223 SUBROUTINE zhpt21( ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP,
224 $ tau, work, rwork, result )
233 INTEGER ITYPE, KBAND, LDU, N
236 DOUBLE PRECISION D( * ), E( * ), RESULT( 2 ), RWORK( * )
237 COMPLEX*16 AP( * ), TAU( * ), U( ldu, * ), VP( * ),
244 DOUBLE PRECISION ZERO, ONE, TEN
245 parameter ( zero = 0.0d+0, one = 1.0d+0, ten = 10.0d+0 )
246 DOUBLE PRECISION HALF
247 parameter ( half = 1.0d+0 / 2.0d+0 )
248 COMPLEX*16 CZERO, CONE
249 parameter ( czero = ( 0.0d+0, 0.0d+0 ),
250 $ cone = ( 1.0d+0, 0.0d+0 ) )
255 INTEGER IINFO, J, JP, JP1, JR, LAP
256 DOUBLE PRECISION ANORM, ULP, UNFL, WNORM
257 COMPLEX*16 TEMP, VSAVE
261 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHP
263 EXTERNAL lsame, dlamch, zlange, zlanhp, zdotc
270 INTRINSIC dble, dcmplx, max, min
282 lap = ( n*( n+1 ) ) / 2
284 IF( lsame( uplo,
'U' ) )
THEN
292 unfl = dlamch(
'Safe minimum' )
293 ulp = dlamch(
'Epsilon' )*dlamch(
'Base' )
297 IF( itype.LT.1 .OR. itype.GT.3 )
THEN
298 result( 1 ) = ten / ulp
306 IF( itype.EQ.3 )
THEN
309 anorm = max( zlanhp(
'1', cuplo, n, ap, rwork ), unfl )
314 IF( itype.EQ.1 )
THEN
318 CALL zlaset(
'Full', n, n, czero, czero, work, n )
319 CALL zcopy( lap, ap, 1, work, 1 )
322 CALL zhpr( cuplo, n, -d( j ), u( 1, j ), 1, work )
325 IF( n.GT.1 .AND. kband.EQ.1 )
THEN
327 CALL zhpr2( cuplo, n, -dcmplx( e( j ) ), u( 1, j ), 1,
328 $ u( 1, j-1 ), 1, work )
331 wnorm = zlanhp(
'1', cuplo, n, work, rwork )
333 ELSE IF( itype.EQ.2 )
THEN
337 CALL zlaset(
'Full', n, n, czero, czero, work, n )
341 DO 40 j = n - 1, 1, -1
342 jp = ( ( 2*n-j )*( j-1 ) ) / 2
344 IF( kband.EQ.1 )
THEN
345 work( jp+j+1 ) = ( cone-tau( j ) )*e( j )
347 work( jp+jr ) = -tau( j )*e( j )*vp( jp+jr )
351 IF( tau( j ).NE.czero )
THEN
354 CALL zhpmv(
'L', n-j, cone, work( jp1+j+1 ),
355 $ vp( jp+j+1 ), 1, czero, work( lap+1 ), 1 )
356 temp = -half*tau( j )*zdotc( n-j, work( lap+1 ), 1,
358 CALL zaxpy( n-j, temp, vp( jp+j+1 ), 1, work( lap+1 ),
360 CALL zhpr2(
'L', n-j, -tau( j ), vp( jp+j+1 ), 1,
361 $ work( lap+1 ), 1, work( jp1+j+1 ) )
365 work( jp+j ) = d( j )
370 jp = ( j*( j-1 ) ) / 2
372 IF( kband.EQ.1 )
THEN
373 work( jp1+j ) = ( cone-tau( j ) )*e( j )
375 work( jp1+jr ) = -tau( j )*e( j )*vp( jp1+jr )
379 IF( tau( j ).NE.czero )
THEN
382 CALL zhpmv(
'U', j, cone, work, vp( jp1+1 ), 1, czero,
384 temp = -half*tau( j )*zdotc( j, work( lap+1 ), 1,
386 CALL zaxpy( j, temp, vp( jp1+1 ), 1, work( lap+1 ),
388 CALL zhpr2(
'U', j, -tau( j ), vp( jp1+1 ), 1,
389 $ work( lap+1 ), 1, work )
392 work( jp1+j+1 ) = d( j+1 )
397 work( j ) = work( j ) - ap( j )
399 wnorm = zlanhp(
'1', cuplo, n, work, rwork )
401 ELSE IF( itype.EQ.3 )
THEN
407 CALL zlacpy(
' ', n, n, u, ldu, work, n )
408 CALL zupmtr(
'R', cuplo,
'C', n, n, vp, tau, work, n,
409 $ work( n**2+1 ), iinfo )
410 IF( iinfo.NE.0 )
THEN
411 result( 1 ) = ten / ulp
416 work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - cone
419 wnorm = zlange(
'1', n, n, work, n, rwork )
422 IF( anorm.GT.wnorm )
THEN
423 result( 1 ) = ( wnorm / anorm ) / ( n*ulp )
425 IF( anorm.LT.one )
THEN
426 result( 1 ) = ( min( wnorm, n*anorm ) / anorm ) / ( n*ulp )
428 result( 1 ) = min( wnorm / anorm, dble( n ) ) / ( n*ulp )
436 IF( itype.EQ.1 )
THEN
437 CALL zgemm(
'N',
'C', n, n, n, cone, u, ldu, u, ldu, czero,
441 work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - cone
444 result( 2 ) = min( zlange(
'1', n, n, work, n, rwork ),
445 $ dble( n ) ) / ( n*ulp )
subroutine zhpr2(UPLO, N, ALPHA, X, INCX, Y, INCY, AP)
ZHPR2
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zhpt21(ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP, TAU, WORK, RWORK, RESULT)
ZHPT21
subroutine zupmtr(SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, INFO)
ZUPMTR
subroutine zhpr(UPLO, N, ALPHA, X, INCX, AP)
ZHPR
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine zhpmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
ZHPMV
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY