178 SUBROUTINE clahef( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
187 INTEGER info, kb, lda, ldw, n, nb
191 COMPLEX a( lda, * ), w( ldw, * )
198 parameter( zero = 0.0e+0, one = 1.0e+0 )
200 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
202 parameter( eight = 8.0e+0, sevten = 17.0e+0 )
205 INTEGER imax,
j, jb, jj, jmax, jp, k, kk, kkw, kp,
207 REAL absakk, alpha, colmax, r1, rowmax, t
208 COMPLEX d11, d21, d22, z
219 INTRINSIC abs, aimag, conjg, max, min,
REAL, sqrt
225 cabs1( z ) = abs(
REAL( Z ) ) + abs( aimag( z ) )
233 alpha = ( one+sqrt( sevten ) ) / eight
235 IF(
lsame( uplo,
'U' ) )
THEN
252 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
259 CALL
ccopy( k-1, a( 1, k ), 1, w( 1, kw ), 1 )
260 w( k, kw ) =
REAL( A( K, K ) )
262 CALL
cgemv(
'No transpose', k, n-k, -cone, a( 1, k+1 ), lda,
263 $ w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
264 w( k, kw ) =
REAL( W( K, KW ) )
270 absakk = abs(
REAL( W( K, KW ) ) )
277 imax =
icamax( k-1, w( 1, kw ), 1 )
278 colmax = cabs1( w( imax, kw ) )
283 IF( max( absakk, colmax ).EQ.zero )
THEN
290 a( k, k ) =
REAL( A( K, K ) )
298 IF( absakk.GE.alpha*colmax )
THEN
310 CALL
ccopy( imax-1, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
311 w( imax, kw-1 ) =
REAL( A( IMAX, IMAX ) )
312 CALL
ccopy( k-imax, a( imax, imax+1 ), lda,
313 $ w( imax+1, kw-1 ), 1 )
314 CALL
clacgv( k-imax, w( imax+1, kw-1 ), 1 )
316 CALL
cgemv(
'No transpose', k, n-k, -cone,
317 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
318 $ cone, w( 1, kw-1 ), 1 )
319 w( imax, kw-1 ) =
REAL( W( IMAX, KW-1 ) )
326 jmax = imax +
icamax( k-imax, w( imax+1, kw-1 ), 1 )
327 rowmax = cabs1( w( jmax, kw-1 ) )
329 jmax =
icamax( imax-1, w( 1, kw-1 ), 1 )
330 rowmax = max( rowmax, cabs1( w( jmax, kw-1 ) ) )
334 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
341 ELSE IF( abs(
REAL( W( IMAX, KW-1 ) ) ).GE.alpha*rowmax )
351 CALL
ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
390 a( kp, kp ) =
REAL( A( KK, KK ) )
391 CALL
ccopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
393 CALL
clacgv( kk-1-kp, a( kp, kp+1 ), lda )
395 $ CALL
ccopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
403 $ CALL
cswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
405 CALL
cswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
409 IF( kstep.EQ.1 )
THEN
427 CALL
ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
434 r1 = one /
REAL( A( K, K ) )
435 CALL
csscal( k-1, r1, a( 1, k ), 1 )
439 CALL
clacgv( k-1, w( 1, kw ), 1 )
504 d11 = w( k, kw ) / conjg( d21 )
505 d22 = w( k-1, kw-1 ) / d21
506 t = one / (
REAL( d11*d22 )-one )
514 a(
j, k-1 ) = d21*( d11*w(
j, kw-1 )-w(
j, kw ) )
515 a(
j, k ) = conjg( d21 )*
516 $ ( d22*w(
j, kw )-w(
j, kw-1 ) )
522 a( k-1, k-1 ) = w( k-1, kw-1 )
523 a( k-1, k ) = w( k-1, kw )
524 a( k, k ) = w( k, kw )
528 CALL
clacgv( k-1, w( 1, kw ), 1 )
529 CALL
clacgv( k-2, w( 1, kw-1 ), 1 )
537 IF( kstep.EQ.1 )
THEN
558 DO 50
j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
559 jb = min( nb, k-
j+1 )
563 DO 40 jj =
j,
j + jb - 1
564 a( jj, jj ) =
REAL( A( JJ, JJ ) )
565 CALL
cgemv(
'No transpose', jj-
j+1, n-k, -cone,
566 $ a(
j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
568 a( jj, jj ) =
REAL( A( JJ, JJ ) )
573 CALL
cgemm(
'No transpose',
'Transpose',
j-1, jb, n-k,
574 $ -cone, a( 1, k+1 ), lda, w(
j, kw+1 ), ldw,
575 $ cone, a( 1,
j ), lda )
598 IF( jp.NE.jj .AND.
j.LE.n )
599 $ CALL
cswap( n-
j+1, a( jp,
j ), lda, a( jj,
j ), lda )
620 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
627 w( k, k ) =
REAL( A( K, K ) )
629 $ CALL
ccopy( n-k, a( k+1, k ), 1, w( k+1, k ), 1 )
630 CALL
cgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ), lda,
631 $ w( k, 1 ), ldw, cone, w( k, k ), 1 )
632 w( k, k ) =
REAL( W( K, K ) )
637 absakk = abs(
REAL( W( K, K ) ) )
644 imax = k +
icamax( n-k, w( k+1, k ), 1 )
645 colmax = cabs1( w( imax, k ) )
650 IF( max( absakk, colmax ).EQ.zero )
THEN
657 a( k, k ) =
REAL( A( K, K ) )
665 IF( absakk.GE.alpha*colmax )
THEN
677 CALL
ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1 )
678 CALL
clacgv( imax-k, w( k, k+1 ), 1 )
679 w( imax, k+1 ) =
REAL( A( IMAX, IMAX ) )
681 $ CALL
ccopy( n-imax, a( imax+1, imax ), 1,
682 $ w( imax+1, k+1 ), 1 )
683 CALL
cgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
684 $ lda, w( imax, 1 ), ldw, cone, w( k, k+1 ),
686 w( imax, k+1 ) =
REAL( W( IMAX, K+1 ) )
692 jmax = k - 1 +
icamax( imax-k, w( k, k+1 ), 1 )
693 rowmax = cabs1( w( jmax, k+1 ) )
695 jmax = imax +
icamax( n-imax, w( imax+1, k+1 ), 1 )
696 rowmax = max( rowmax, cabs1( w( jmax, k+1 ) ) )
700 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
707 ELSE IF( abs(
REAL( W( IMAX, K+1 ) ) ).GE.alpha*rowmax )
717 CALL
ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
752 a( kp, kp ) =
REAL( A( KK, KK ) )
753 CALL
ccopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
755 CALL
clacgv( kp-kk-1, a( kp, kk+1 ), lda )
757 $ CALL
ccopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
765 $ CALL
cswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
766 CALL
cswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
769 IF( kstep.EQ.1 )
THEN
787 CALL
ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
794 r1 = one /
REAL( A( K, K ) )
795 CALL
csscal( n-k, r1, a( k+1, k ), 1 )
799 CALL
clacgv( n-k, w( k+1, k ), 1 )
864 d11 = w( k+1, k+1 ) / d21
865 d22 = w( k, k ) / conjg( d21 )
866 t = one / (
REAL( d11*d22 )-one )
874 a(
j, k ) = conjg( d21 )*
875 $ ( d11*w(
j, k )-w(
j, k+1 ) )
876 a(
j, k+1 ) = d21*( d22*w(
j, k+1 )-w(
j, k ) )
882 a( k, k ) = w( k, k )
883 a( k+1, k ) = w( k+1, k )
884 a( k+1, k+1 ) = w( k+1, k+1 )
888 CALL
clacgv( n-k, w( k+1, k ), 1 )
889 CALL
clacgv( n-k-1, w( k+2, k+1 ), 1 )
897 IF( kstep.EQ.1 )
THEN
919 jb = min( nb, n-
j+1 )
923 DO 100 jj =
j,
j + jb - 1
924 a( jj, jj ) =
REAL( A( JJ, JJ ) )
925 CALL
cgemv(
'No transpose',
j+jb-jj, k-1, -cone,
926 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
928 a( jj, jj ) =
REAL( A( JJ, JJ ) )
934 $ CALL
cgemm(
'No transpose',
'Transpose', n-
j-jb+1, jb,
935 $ k-1, -cone, a(
j+jb, 1 ), lda, w(
j, 1 ),
936 $ ldw, cone, a(
j+jb,
j ), lda )
959 IF( jp.NE.jj .AND.
j.GE.1 )
960 $ CALL
cswap(
j, a( jp, 1 ), lda, a( jj, 1 ), lda )
LOGICAL function lsame(CA, CB)
LSAME
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
subroutine clahef(UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
CLAHEF computes a partial factorization of a complex Hermitian indefinite matrix using the Bunch-Kauf...
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine csscal(N, SA, CX, INCX)
CSSCAL
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
INTEGER function icamax(N, CX, INCX)
ICAMAX