150 SUBROUTINE zhbt21( UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK,
160 INTEGER ka, ks, lda, ldu, n
163 DOUBLE PRECISION d( * ), e( * ), result( 2 ), rwork( * )
164 COMPLEX*16 a( lda, * ), u( ldu, * ), work( * )
170 COMPLEX*16 czero, cone
171 parameter( czero = ( 0.0d+0, 0.0d+0 ),
172 $ cone = ( 1.0d+0, 0.0d+0 ) )
173 DOUBLE PRECISION zero, one
174 parameter( zero = 0.0d+0, one = 1.0d+0 )
179 INTEGER ika, j, jc, jr
180 DOUBLE PRECISION anorm, ulp, unfl, wnorm
191 INTRINSIC dble, dcmplx, max, min
202 ika = max( 0, min( n-1, ka ) )
204 IF(
lsame( uplo,
'U' ) )
THEN
212 unfl =
dlamch(
'Safe minimum' )
221 anorm = max(
zlanhb(
'1', cuplo, n, ika, a, lda, rwork ), unfl )
230 DO 10 jr = 1, min( ika+1, n+1-jc )
232 work( j ) = a( jr, jc )
234 DO 20 jr = ika + 2, n + 1 - jc
239 DO 30 jr = ika + 2, jc
243 DO 40 jr = min( ika, jc-1 ), 0, -1
245 work( j ) = a( ika+1-jr, jc )
251 CALL
zhpr( cuplo, n, -d( j ), u( 1, j ), 1, work )
254 IF( n.GT.1 .AND. ks.EQ.1 )
THEN
256 CALL
zhpr2( cuplo, n, -dcmplx( e( j ) ), u( 1, j ), 1,
257 $ u( 1, j+1 ), 1, work )
260 wnorm =
zlanhp(
'1', cuplo, n, work, rwork )
262 IF( anorm.GT.wnorm )
THEN
263 result( 1 ) = ( wnorm / anorm ) / ( n*ulp )
265 IF( anorm.LT.one )
THEN
266 result( 1 ) = ( min( wnorm, n*anorm ) / anorm ) / ( n*ulp )
268 result( 1 ) = min( wnorm / anorm, dble( n ) ) / ( n*ulp )
276 CALL
zgemm(
'N',
'C', n, n, n, cone, u, ldu, u, ldu, czero, work,
280 work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - cone
283 result( 2 ) = min(
zlange(
'1', n, n, work, n, rwork ),
284 $ dble( n ) ) / ( n*ulp )