192 SUBROUTINE zhetf2( UPLO, N, A, LDA, IPIV, INFO )
205 COMPLEX*16 a( lda, * )
211 DOUBLE PRECISION zero, one
212 parameter( zero = 0.0d+0, one = 1.0d+0 )
213 DOUBLE PRECISION eight, sevten
214 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
218 INTEGER i, imax,
j, jmax, k, kk, kp, kstep
219 DOUBLE PRECISION absakk, alpha, colmax, d, d11, d22, r1, rowmax,
221 COMPLEX*16 d12, d21, t, wk, wkm1, wkp1, zdum
233 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, sqrt
236 DOUBLE PRECISION cabs1
239 cabs1( zdum ) = abs( dble( zdum ) ) + abs( dimag( zdum ) )
246 upper =
lsame( uplo,
'U' )
247 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
249 ELSE IF( n.LT.0 )
THEN
251 ELSE IF( lda.LT.max( 1, n ) )
THEN
255 CALL
xerbla(
'ZHETF2', -info )
261 alpha = ( one+sqrt( sevten ) ) / eight
282 absakk = abs( dble( a( k, k ) ) )
289 imax =
izamax( k-1, a( 1, k ), 1 )
290 colmax = cabs1( a( imax, k ) )
295 IF( (max( absakk, colmax ).EQ.zero) .OR.
disnan(absakk) )
THEN
303 a( k, k ) = dble( a( k, k ) )
310 IF( absakk.GE.alpha*colmax )
THEN
321 jmax = imax +
izamax( k-imax, a( imax, imax+1 ), lda )
322 rowmax = cabs1( a( imax, jmax ) )
324 jmax =
izamax( imax-1, a( 1, imax ), 1 )
325 rowmax = max( rowmax, cabs1( a( jmax, imax ) ) )
328 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
334 ELSE IF( abs( dble( a( imax, imax ) ) ).GE.alpha*rowmax )
360 CALL
zswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
361 DO 20
j = kp + 1, kk - 1
362 t = dconjg( a(
j, kk ) )
363 a(
j, kk ) = dconjg( a( kp,
j ) )
366 a( kp, kk ) = dconjg( a( kp, kk ) )
367 r1 = dble( a( kk, kk ) )
368 a( kk, kk ) = dble( a( kp, kp ) )
370 IF( kstep.EQ.2 )
THEN
371 a( k, k ) = dble( a( k, k ) )
373 a( k-1, k ) = a( kp, k )
377 a( k, k ) = dble( a( k, k ) )
379 $ a( k-1, k-1 ) = dble( a( k-1, k-1 ) )
384 IF( kstep.EQ.1 )
THEN
396 r1 = one / dble( a( k, k ) )
397 CALL
zher( uplo, k-1, -r1, a( 1, k ), 1, a, lda )
401 CALL
zdscal( k-1, r1, a( 1, k ), 1 )
418 d =
dlapy2( dble( a( k-1, k ) ),
419 $ dimag( a( k-1, k ) ) )
420 d22 = dble( a( k-1, k-1 ) ) / d
421 d11 = dble( a( k, k ) ) / d
422 tt = one / ( d11*d22-one )
423 d12 = a( k-1, k ) / d
426 DO 40
j = k - 2, 1, -1
427 wkm1 = d*( d11*a(
j, k-1 )-dconjg( d12 )*
429 wk = d*( d22*a(
j, k )-d12*a(
j, k-1 ) )
431 a( i,
j ) = a( i,
j ) - a( i, k )*dconjg( wk ) -
432 $ a( i, k-1 )*dconjg( wkm1 )
436 a(
j,
j ) = dcmplx( dble( a(
j,
j ) ), 0.0d+0 )
446 IF( kstep.EQ.1 )
THEN
477 absakk = abs( dble( a( k, k ) ) )
484 imax = k +
izamax( n-k, a( k+1, k ), 1 )
485 colmax = cabs1( a( imax, k ) )
490 IF( (max( absakk, colmax ).EQ.zero) .OR.
disnan(absakk) )
THEN
498 a( k, k ) = dble( a( k, k ) )
505 IF( absakk.GE.alpha*colmax )
THEN
516 jmax = k - 1 +
izamax( imax-k, a( imax, k ), lda )
517 rowmax = cabs1( a( imax, jmax ) )
519 jmax = imax +
izamax( n-imax, a( imax+1, imax ), 1 )
520 rowmax = max( rowmax, cabs1( a( jmax, imax ) ) )
523 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
529 ELSE IF( abs( dble( a( imax, imax ) ) ).GE.alpha*rowmax )
556 $ CALL
zswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
557 DO 60
j = kk + 1, kp - 1
558 t = dconjg( a(
j, kk ) )
559 a(
j, kk ) = dconjg( a( kp,
j ) )
562 a( kp, kk ) = dconjg( a( kp, kk ) )
563 r1 = dble( a( kk, kk ) )
564 a( kk, kk ) = dble( a( kp, kp ) )
566 IF( kstep.EQ.2 )
THEN
567 a( k, k ) = dble( a( k, k ) )
569 a( k+1, k ) = a( kp, k )
573 a( k, k ) = dble( a( k, k ) )
575 $ a( k+1, k+1 ) = dble( a( k+1, k+1 ) )
580 IF( kstep.EQ.1 )
THEN
594 r1 = one / dble( a( k, k ) )
595 CALL
zher( uplo, n-k, -r1, a( k+1, k ), 1,
596 $ a( k+1, k+1 ), lda )
600 CALL
zdscal( n-k, r1, a( k+1, k ), 1 )
616 d =
dlapy2( dble( a( k+1, k ) ),
617 $ dimag( a( k+1, k ) ) )
618 d11 = dble( a( k+1, k+1 ) ) / d
619 d22 = dble( a( k, k ) ) / d
620 tt = one / ( d11*d22-one )
621 d21 = a( k+1, k ) / d
625 wk = d*( d11*a(
j, k )-d21*a(
j, k+1 ) )
626 wkp1 = d*( d22*a(
j, k+1 )-dconjg( d21 )*
629 a( i,
j ) = a( i,
j ) - a( i, k )*dconjg( wk ) -
630 $ a( i, k+1 )*dconjg( wkp1 )
634 a(
j,
j ) = dcmplx( dble( a(
j,
j ) ), 0.0d+0 )
642 IF( kstep.EQ.1 )
THEN
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
integer function izamax(N, ZX, INCX)
IZAMAX
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine zhetf2(UPLO, N, A, LDA, IPIV, INFO)
ZHETF2 computes the factorization of a complex Hermitian matrix, using the diagonal pivoting method (...
subroutine zher(UPLO, N, ALPHA, X, INCX, A, LDA)
ZHER
logical function disnan(DIN)
DISNAN tests input for NaN.