215 parameter( zero = 0.0e+0, one = 1.0e+0 )
217 parameter( eight = 8.0e+0, sevten = 17.0e+0 )
221 INTEGER i, ii, imax, itemp,
j, jmax, k, kk, kp, kstep,
223 REAL absakk, alpha, colmax, d, d11, d22, r1, stemp,
225 COMPLEX d12, d21, t, wk, wkm1, wkp1, z
238 INTRINSIC abs, aimag, cmplx, conjg, max,
REAL, sqrt
244 cabs1( z ) = abs(
REAL( Z ) ) + abs( aimag( z ) )
251 upper =
lsame( uplo,
'U' )
252 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
254 ELSE IF( n.LT.0 )
THEN
256 ELSE IF( lda.LT.max( 1, n ) )
THEN
260 CALL
xerbla(
'CHETF2_ROOK', -info )
266 alpha = ( one+sqrt( sevten ) ) / eight
292 absakk = abs(
REAL( A( K, K ) ) )
299 imax =
icamax( k-1, a( 1, k ), 1 )
300 colmax = cabs1( a( imax, k ) )
305 IF( ( max( absakk, colmax ).EQ.zero ) )
THEN
312 a( k, k ) =
REAL( A( K, K ) )
323 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
345 jmax = imax +
icamax( k-imax, a( imax, imax+1 ),
347 rowmax = cabs1( a( imax, jmax ) )
353 itemp =
icamax( imax-1, a( 1, imax ), 1 )
354 stemp = cabs1( a( itemp, imax ) )
355 IF( stemp.GT.rowmax )
THEN
366 IF( .NOT.( abs(
REAL( A( IMAX, IMAX ) ) )
367 $ .LT.alpha*rowmax ) )
THEN
379 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
401 IF( .NOT.done ) goto 12
416 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
419 $ CALL
cswap( p-1, a( 1, k ), 1, a( 1, p ), 1 )
421 DO 14
j = p + 1, k - 1
422 t = conjg( a(
j, k ) )
423 a(
j, k ) = conjg( a( p,
j ) )
427 a( p, k ) = conjg( a( p, k ) )
429 r1 =
REAL( A( K, K ) )
430 a( k, k ) =
REAL( A( P, P ) )
440 $ CALL
cswap( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
442 DO 15
j = kp + 1, kk - 1
443 t = conjg( a(
j, kk ) )
444 a(
j, kk ) = conjg( a( kp,
j ) )
448 a( kp, kk ) = conjg( a( kp, kk ) )
450 r1 =
REAL( A( KK, KK ) )
451 a( kk, kk ) =
REAL( A( KP, KP ) )
454 IF( kstep.EQ.2 )
THEN
456 a( k, k ) =
REAL( A( K, K ) )
459 a( k-1, k ) = a( kp, k )
464 a( k, k ) =
REAL( A( K, K ) )
466 $ a( k-1, k-1 ) =
REAL( A( K-1, K-1 ) )
471 IF( kstep.EQ.1 )
THEN
484 IF( abs(
REAL( A( K, K ) ) ).GE.sfmin ) then
490 d11 = one /
REAL( A( K, K ) )
491 CALL
cher( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
495 CALL
csscal( k-1, d11, a( 1, k ), 1 )
500 d11 =
REAL( A( K, K ) )
502 a( ii, k ) = a( ii, k ) / d11
510 CALL
cher( uplo, k-1, -d11, a( 1, k ), 1, a, lda )
532 d =
slapy2(
REAL( A( K-1, K ) ),
533 $ aimag( a( k-1, k ) ) )
535 d22 = a( k-1, k-1 ) / d
536 d12 = a( k-1, k ) / d
537 tt = one / ( d11*d22-one )
539 DO 30
j = k - 2, 1, -1
543 wkm1 = tt*( d11*a(
j, k-1 )-conjg( d12 )*
545 wk = tt*( d22*a(
j, k )-d12*a(
j, k-1 ) )
550 a( i,
j ) = a( i,
j ) -
551 $ ( a( i, k ) / d )*conjg( wk ) -
552 $ ( a( i, k-1 ) / d )*conjg( wkm1 )
558 a(
j, k-1 ) = wkm1 / d
560 a(
j,
j ) = cmplx(
REAL( A( J, J ) ), zero )
572 IF( kstep.EQ.1 )
THEN
604 absakk = abs(
REAL( A( K, K ) ) )
611 imax = k +
icamax( n-k, a( k+1, k ), 1 )
612 colmax = cabs1( a( imax, k ) )
617 IF( max( absakk, colmax ).EQ.zero )
THEN
624 a( k, k ) =
REAL( A( K, K ) )
635 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
657 jmax = k - 1 +
icamax( imax-k, a( imax, k ), lda )
658 rowmax = cabs1( a( imax, jmax ) )
664 itemp = imax +
icamax( n-imax, a( imax+1, imax ),
666 stemp = cabs1( a( itemp, imax ) )
667 IF( stemp.GT.rowmax )
THEN
678 IF( .NOT.( abs(
REAL( A( IMAX, IMAX ) ) )
679 $ .LT.alpha*rowmax ) )
THEN
691 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
714 IF( .NOT.done ) goto 42
729 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
732 $ CALL
cswap( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
734 DO 44
j = k + 1, p - 1
735 t = conjg( a(
j, k ) )
736 a(
j, k ) = conjg( a( p,
j ) )
740 a( p, k ) = conjg( a( p, k ) )
742 r1 =
REAL( A( K, K ) )
743 a( k, k ) =
REAL( A( P, P ) )
753 $ CALL
cswap( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
755 DO 45
j = kk + 1, kp - 1
756 t = conjg( a(
j, kk ) )
757 a(
j, kk ) = conjg( a( kp,
j ) )
761 a( kp, kk ) = conjg( a( kp, kk ) )
763 r1 =
REAL( A( KK, KK ) )
764 a( kk, kk ) =
REAL( A( KP, KP ) )
767 IF( kstep.EQ.2 )
THEN
769 a( k, k ) =
REAL( A( K, K ) )
772 a( k+1, k ) = a( kp, k )
777 a( k, k ) =
REAL( A( K, K ) )
779 $ a( k+1, k+1 ) =
REAL( A( K+1, K+1 ) )
784 IF( kstep.EQ.1 )
THEN
799 IF( abs(
REAL( A( K, K ) ) ).GE.sfmin ) then
805 d11 = one /
REAL( A( K, K ) )
806 CALL
cher( uplo, n-k, -d11, a( k+1, k ), 1,
807 $ a( k+1, k+1 ), lda )
811 CALL
csscal( n-k, d11, a( k+1, k ), 1 )
816 d11 =
REAL( A( K, K ) )
818 a( ii, k ) = a( ii, k ) / d11
826 CALL
cher( uplo, n-k, -d11, a( k+1, k ), 1,
827 $ a( k+1, k+1 ), lda )
850 d =
slapy2(
REAL( A( K+1, K ) ),
851 $ aimag( a( k+1, k ) ) )
852 d11 =
REAL( A( K+1, K+1 ) ) / d
853 d22 =
REAL( A( K, K ) ) / d
854 d21 = a( k+1, k ) / d
855 tt = one / ( d11*d22-one )
861 wk = tt*( d11*a(
j, k )-d21*a(
j, k+1 ) )
862 wkp1 = tt*( d22*a(
j, k+1 )-conjg( d21 )*
868 a( i,
j ) = a( i,
j ) -
869 $ ( a( i, k ) / d )*conjg( wk ) -
870 $ ( a( i, k+1 ) / d )*conjg( wkp1 )
876 a(
j, k+1 ) = wkp1 / d
878 a(
j,
j ) = cmplx(
REAL( A( J, J ) ), zero )
890 IF( kstep.EQ.1 )
THEN
LOGICAL function lsame(CA, CB)
LSAME
subroutine chetf2_rook(UPLO, N, A, LDA, IPIV, INFO)
CHETF2_ROOK computes the factorization of a complex Hermitian indefinite matrix using the bounded Bun...
REAL function slamch(CMACH)
SLAMCH
REAL function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cher(UPLO, N, ALPHA, X, INCX, A, LDA)
CHER
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine csscal(N, SA, CX, INCX)
CSSCAL
INTEGER function icamax(N, CX, INCX)
ICAMAX