146 SUBROUTINE dglmts( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U,
147 $ work, lwork, rwork, result )
155 INTEGER lda, ldb, lwork, m, n, p
156 DOUBLE PRECISION result
162 DOUBLE PRECISION a( lda, * ), af( lda, * ), b( ldb, * ),
163 $ bf( ldb, * ), d( * ), df( * ), rwork( * ),
164 $ u( * ), work( lwork ), x( * )
167 DOUBLE PRECISION zero, one
168 parameter( zero = 0.0d+0, one = 1.0d+0 )
172 DOUBLE PRECISION anorm, bnorm, dnorm, eps, unfl, xnorm, ynorm
188 unfl =
dlamch(
'Safe minimum' )
189 anorm = max(
dlange(
'1', n, m, a, lda, rwork ), unfl )
190 bnorm = max(
dlange(
'1', n, p, b, ldb, rwork ), unfl )
195 CALL
dlacpy(
'Full', n, m, a, lda, af, lda )
196 CALL
dlacpy(
'Full', n, p, b, ldb, bf, ldb )
197 CALL
dcopy( n, d, 1, df, 1 )
201 CALL
dggglm( n, m, p, af, lda, bf, ldb, df, x, u, work, lwork,
210 CALL
dcopy( n, d, 1, df, 1 )
211 CALL
dgemv(
'No transpose', n, m, -one, a, lda, x, 1, one, df, 1 )
213 CALL
dgemv(
'No transpose', n, p, -one, b, ldb, u, 1, one, df, 1 )
215 dnorm =
dasum( n, df, 1 )
217 ynorm = anorm + bnorm
219 IF( xnorm.LE.zero )
THEN
222 result = ( ( dnorm / ynorm ) / xnorm ) / eps