135 SUBROUTINE dbdt03( UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK,
145 INTEGER kd, ldu, ldvt, n
146 DOUBLE PRECISION resid
149 DOUBLE PRECISION d( * ), e( * ), s( * ), u( ldu, * ),
150 $ vt( ldvt, * ), work( * )
156 DOUBLE PRECISION zero, one
157 parameter( zero = 0.0d+0, one = 1.0d+0 )
161 DOUBLE PRECISION bnorm, eps
173 INTRINSIC abs, dble, max, min
190 IF(
lsame( uplo,
'U' ) )
THEN
196 work( n+i ) = s( i )*vt( i, j )
198 CALL
dgemv(
'No transpose', n, n, -one, u, ldu,
199 $ work( n+1 ), 1, zero, work, 1 )
200 work( j ) = work( j ) + d( j )
202 work( j-1 ) = work( j-1 ) + e( j-1 )
203 bnorm = max( bnorm, abs( d( j ) )+abs( e( j-1 ) ) )
205 bnorm = max( bnorm, abs( d( j ) ) )
207 resid = max( resid,
dasum( n, work, 1 ) )
215 work( n+i ) = s( i )*vt( i, j )
217 CALL
dgemv(
'No transpose', n, n, -one, u, ldu,
218 $ work( n+1 ), 1, zero, work, 1 )
219 work( j ) = work( j ) + d( j )
221 work( j+1 ) = work( j+1 ) + e( j )
222 bnorm = max( bnorm, abs( d( j ) )+abs( e( j ) ) )
224 bnorm = max( bnorm, abs( d( j ) ) )
226 resid = max( resid,
dasum( n, work, 1 ) )
235 work( n+i ) = s( i )*vt( i, j )
237 CALL
dgemv(
'No transpose', n, n, -one, u, ldu, work( n+1 ),
239 work( j ) = work( j ) + d( j )
240 resid = max( resid,
dasum( n, work, 1 ) )
243 bnorm = abs( d( j ) )
248 eps =
dlamch(
'Precision' )
250 IF( bnorm.LE.zero )
THEN
254 IF( bnorm.GE.resid )
THEN
255 resid = ( resid / bnorm ) / ( dble( n )*eps )
257 IF( bnorm.LT.one )
THEN
258 resid = ( min( resid, dble( n )*bnorm ) / bnorm ) /
261 resid = min( resid / bnorm, dble( n ) ) /