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
213 EXTERNAL lsame, icamax
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 )
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
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...
subroutine csscal(N, SA, CX, INCX)
CSSCAL