156 SUBROUTINE dget54( N, A, LDA, B, LDB, S, LDS, T, LDT, U, LDU, V,
157 $ ldv, work, result )
165 INTEGER lda, ldb, lds, ldt, ldu, ldv, n
166 DOUBLE PRECISION result
169 DOUBLE PRECISION a( lda, * ), b( ldb, * ), s( lds, * ),
170 $ t( ldt, * ), u( ldu, * ), v( ldv, * ),
177 DOUBLE PRECISION zero, one
178 parameter( zero = 0.0d+0, one = 1.0d+0 )
181 DOUBLE PRECISION abnorm, ulp, unfl, wnorm
184 DOUBLE PRECISION dum( 1 )
194 INTRINSIC dble, max, min
204 unfl =
dlamch(
'Safe minimum' )
209 CALL
dlacpy(
'Full', n, n, a, lda, work, n )
210 CALL
dlacpy(
'Full', n, n, b, ldb, work( n*n+1 ), n )
211 abnorm = max(
dlange(
'1', n, 2*n, work, n, dum ), unfl )
215 CALL
dlacpy(
' ', n, n, a, lda, work, n )
216 CALL
dgemm(
'N',
'N', n, n, n, one, u, ldu, s, lds, zero,
219 CALL
dgemm(
'N',
'C', n, n, n, -one, work( n*n+1 ), n, v, ldv,
224 CALL
dlacpy(
' ', n, n, b, ldb, work( n*n+1 ), n )
225 CALL
dgemm(
'N',
'N', n, n, n, one, u, ldu, t, ldt, zero,
226 $ work( 2*n*n+1 ), n )
228 CALL
dgemm(
'N',
'C', n, n, n, -one, work( 2*n*n+1 ), n, v, ldv,
229 $ one, work( n*n+1 ), n )
233 wnorm =
dlange(
'1', n, 2*n, work, n, dum )
235 IF( abnorm.GT.wnorm )
THEN
236 result = ( wnorm / abnorm ) / ( 2*n*ulp )
238 IF( abnorm.LT.one )
THEN
239 result = ( min( wnorm, 2*n*abnorm ) / abnorm ) / ( 2*n*ulp )
241 result = min( wnorm / abnorm, dble( 2*n ) ) / ( 2*n*ulp )