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 )