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
270 INTRINSIC dble, dcmplx, max, min
282 lap = ( n*( n+1 ) ) / 2
284 IF(
lsame( uplo,
'U' ) )
THEN
292 unfl =
dlamch(
'Safe minimum' )
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 )
LOGICAL function lsame(CA, CB)
LSAME
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
DOUBLE PRECISION function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zupmtr(SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, INFO)
ZUPMTR
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
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 zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
DOUBLE PRECISION function zlanhp(NORM, UPLO, N, AP, WORK)
ZLANHP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix supplied in packed form.
subroutine zhpr2(UPLO, N, ALPHA, X, INCX, Y, INCY, AP)
ZHPR2
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
COMPLEX *16 function zdotc(N, ZX, INCX, ZY, INCY)
ZDOTC
subroutine zhpt21(ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP, TAU, WORK, RWORK, RESULT)
ZHPT21
subroutine zhpr(UPLO, N, ALPHA, X, INCX, AP)
ZHPR
subroutine zhpmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
ZHPMV