135 SUBROUTINE cbdt03( UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK,
145 INTEGER kd, ldu, ldvt, n
149 REAL d( * ), e( * ), s( * )
150 COMPLEX u( ldu, * ), vt( ldvt, * ), work( * )
157 parameter( zero = 0.0e+0, one = 1.0e+0 )
173 INTRINSIC abs, cmplx, max, min, real
190 IF(
lsame( uplo,
'U' ) )
THEN
196 work( n+i ) = s( i )*vt( i,
j )
198 CALL
cgemv(
'No transpose', n, n, -cmplx( one ), u, ldu,
199 $ work( n+1 ), 1, cmplx( 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,
scasum( n, work, 1 ) )
215 work( n+i ) = s( i )*vt( i,
j )
217 CALL
cgemv(
'No transpose', n, n, -cmplx( one ), u, ldu,
218 $ work( n+1 ), 1, cmplx( 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,
scasum( n, work, 1 ) )
235 work( n+i ) = s( i )*vt( i,
j )
237 CALL
cgemv(
'No transpose', n, n, -cmplx( one ), u, ldu,
238 $ work( n+1 ), 1, cmplx( zero ), work, 1 )
239 work(
j ) = work(
j ) + d(
j )
240 resid = max( resid,
scasum( n, work, 1 ) )
243 bnorm = abs( d(
j ) )
248 eps =
slamch(
'Precision' )
250 IF( bnorm.LE.zero )
THEN
254 IF( bnorm.GE.resid )
THEN
255 resid = ( resid / bnorm ) / (
REAL( n )*eps )
257 IF( bnorm.LT.one )
THEN
258 resid = ( min( resid,
REAL( n )*bnorm ) / bnorm ) /
261 resid = min( resid / bnorm,
REAL( N ) ) /
LOGICAL function lsame(CA, CB)
LSAME
INTEGER function isamax(N, SX, INCX)
ISAMAX
REAL function slamch(CMACH)
SLAMCH
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
REAL function scasum(N, CX, INCX)
SCASUM
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine cbdt03(UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK, RESID)
CBDT03