131 SUBROUTINE dbdt04( UPLO, N, D, E, S, NS, U, LDU, VT, LDVT, WORK,
141 INTEGER LDU, LDVT, N, NS
142 DOUBLE PRECISION RESID
145 DOUBLE PRECISION D( * ), E( * ), S( * ), U( ldu, * ),
146 $ vt( ldvt, * ), work( * )
152 DOUBLE PRECISION ZERO, ONE
153 parameter ( zero = 0.0d+0, one = 1.0d+0 )
157 DOUBLE PRECISION BNORM, EPS
162 DOUBLE PRECISION DASUM, DLAMCH
163 EXTERNAL lsame, idamax, dasum, dlamch
169 INTRINSIC abs, dble, max, min
176 IF( n.LE.0 .OR. ns.LE.0 )
179 eps = dlamch(
'Precision' )
185 IF( lsame( uplo,
'U' ) )
THEN
193 work( k ) = d( j )*vt( i, j ) + e( j )*vt( i, j+1 )
196 work( k ) = d( n )*vt( i, n )
198 bnorm = abs( d( 1 ) )
200 bnorm = max( bnorm, abs( d( i ) )+abs( e( i-1 ) ) )
209 work( k ) = d( 1 )*vt( i, 1 )
212 work( k ) = e( j )*vt( i, j ) + d( j+1 )*vt( i, j+1 )
215 bnorm = abs( d( n ) )
217 bnorm = max( bnorm, abs( d( i ) )+abs( e( i ) ) )
221 CALL dgemm(
'T',
'N', ns, ns, n, -one, u, ldu, work( 1 ),
222 $ n, zero, work( 1+n*ns ), ns )
228 work( k+i ) = work( k+i ) + s( i )
229 resid = max( resid, dasum( ns, work( k+1 ), 1 ) )
233 IF( bnorm.LE.zero )
THEN
237 IF( bnorm.GE.resid )
THEN
238 resid = ( resid / bnorm ) / ( dble( n )*eps )
240 IF( bnorm.LT.one )
THEN
241 resid = ( min( resid, dble( n )*bnorm ) / bnorm ) /
244 resid = min( resid / bnorm, dble( n ) ) /
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dbdt04(UPLO, N, D, E, S, NS, U, LDU, VT, LDVT, WORK, RESID)
DBDT04