135 SUBROUTINE zbdt03( 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( * )
150 COMPLEX*16 u( ldu, * ), 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, dcmplx, max, min
190 IF(
lsame( uplo,
'U' ) )
THEN
196 work( n+i ) = s( i )*vt( i,
j )
198 CALL
zgemv(
'No transpose', n, n, -dcmplx( one ), u, ldu,
199 $ work( n+1 ), 1, dcmplx( 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,
dzasum( n, work, 1 ) )
215 work( n+i ) = s( i )*vt( i,
j )
217 CALL
zgemv(
'No transpose', n, n, -dcmplx( one ), u, ldu,
218 $ work( n+1 ), 1, dcmplx( 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,
dzasum( n, work, 1 ) )
235 work( n+i ) = s( i )*vt( i,
j )
237 CALL
zgemv(
'No transpose', n, n, -dcmplx( one ), u, ldu,
238 $ work( n+1 ), 1, dcmplx( zero ), work, 1 )
239 work(
j ) = work(
j ) + d(
j )
240 resid = max( resid,
dzasum( 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 ) ) /
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
logical function lsame(CA, CB)
LSAME
double precision function dzasum(N, ZX, INCX)
DZASUM
subroutine zbdt03(UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK, RESID)
ZBDT03
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
integer function idamax(N, DX, INCX)
IDAMAX