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 )
integer function icamax(N, CX, INCX)
ICAMAX
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
logical function lsame(CA, CB)
LSAME
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