139 SUBROUTINE dstt22( N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK,
148 INTEGER kband, ldu, ldwork, m, n
151 DOUBLE PRECISION ad( * ), ae( * ), result( 2 ), sd( * ),
152 $ se( * ), u( ldu, * ), work( ldwork, * )
158 DOUBLE PRECISION zero, one
159 parameter( zero = 0.0d0, one = 1.0d0 )
163 DOUBLE PRECISION anorm, aukj, ulp, unfl, wnorm
173 INTRINSIC abs, dble, max, min
179 IF( n.LE.0 .OR. m.LE.0 )
182 unfl =
dlamch(
'Safe minimum' )
190 anorm = abs( ad( 1 ) ) + abs( ae( 1 ) )
192 anorm = max( anorm, abs( ad(
j ) )+abs( ae(
j ) )+
195 anorm = max( anorm, abs( ad( n ) )+abs( ae( n-1 ) ) )
197 anorm = abs( ad( 1 ) )
199 anorm = max( anorm, unfl )
207 aukj = ad( k )*u( k,
j )
209 $ aukj = aukj + ae( k )*u( k+1,
j )
211 $ aukj = aukj + ae( k-1 )*u( k-1,
j )
212 work( i,
j ) = work( i,
j ) + u( k, i )*aukj
215 work( i, i ) = work( i, i ) - sd( i )
216 IF( kband.EQ.1 )
THEN
218 $ work( i, i-1 ) = work( i, i-1 ) - se( i-1 )
220 $ work( i, i+1 ) = work( i, i+1 ) - se( i )
224 wnorm =
dlansy(
'1',
'L', m, work, m, work( 1, m+1 ) )
226 IF( anorm.GT.wnorm )
THEN
227 result( 1 ) = ( wnorm / anorm ) / ( m*ulp )
229 IF( anorm.LT.one )
THEN
230 result( 1 ) = ( min( wnorm, m*anorm ) / anorm ) / ( m*ulp )
232 result( 1 ) = min( wnorm / anorm, dble( m ) ) / ( m*ulp )
240 CALL
dgemm(
'T',
'N', m, m, n, one, u, ldu, u, ldu, zero, work,
244 work(
j,
j ) = work(
j,
j ) - one
247 result( 2 ) = min( dble( m ),
dlange(
'1', m, m, work, m, work( 1,
248 $ m+1 ) ) ) / ( m*ulp )
DOUBLE PRECISION function dlansy(NORM, UPLO, N, A, LDA, WORK)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix.
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dstt22(N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK, LDWORK, RESULT)
DSTT22
DOUBLE PRECISION function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH