132 SUBROUTINE zstt21( N, KBAND, AD, AE, SD, SE, U, LDU, WORK, RWORK,
141 INTEGER kband, ldu, n
144 DOUBLE PRECISION ad( * ), ae( * ), result( 2 ), rwork( * ),
146 COMPLEX*16 u( ldu, * ), work( * )
152 DOUBLE PRECISION zero, one
153 parameter( zero = 0.0d+0, one = 1.0d+0 )
154 COMPLEX*16 czero, cone
155 parameter( czero = ( 0.0d+0, 0.0d+0 ),
156 $ cone = ( 1.0d+0, 0.0d+0 ) )
160 DOUBLE PRECISION anorm, temp1, temp2, ulp, unfl, wnorm
170 INTRINSIC abs, dble, dcmplx, max, min
181 unfl =
dlamch(
'Safe minimum' )
182 ulp =
dlamch(
'Precision' )
188 CALL
zlaset(
'Full', n, n, czero, czero, work, n )
194 work( ( n+1 )*( j-1 )+1 ) = ad( j )
195 work( ( n+1 )*( j-1 )+2 ) = ae( j )
196 temp2 = abs( ae( j ) )
197 anorm = max( anorm, abs( ad( j ) )+temp1+temp2 )
201 work( n**2 ) = ad( n )
202 anorm = max( anorm, abs( ad( n ) )+temp1, unfl )
207 CALL
zher(
'L', n, -sd( j ), u( 1, j ), 1, work, n )
210 IF( n.GT.1 .AND. kband.EQ.1 )
THEN
212 CALL
zher2(
'L', n, -dcmplx( se( j ) ), u( 1, j ), 1,
213 $ u( 1, j+1 ), 1, work, n )
217 wnorm =
zlanhe(
'1',
'L', n, work, n, rwork )
219 IF( anorm.GT.wnorm )
THEN
220 result( 1 ) = ( wnorm / anorm ) / ( n*ulp )
222 IF( anorm.LT.one )
THEN
223 result( 1 ) = ( min( wnorm, n*anorm ) / anorm ) / ( n*ulp )
225 result( 1 ) = min( wnorm / anorm, dble( n ) ) / ( n*ulp )
233 CALL
zgemm(
'N',
'C', n, n, n, cone, u, ldu, u, ldu, czero, work,
237 work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - cone
240 result( 2 ) = min( dble( n ),
zlange(
'1', n, n, work, n,
241 $ rwork ) ) / ( n*ulp )