262 SUBROUTINE clahef_rk( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
272 INTEGER INFO, KB, LDA, LDW, N, NB
276 COMPLEX A( lda, * ), W( ldw, * ), E( * )
283 parameter ( zero = 0.0e+0, one = 1.0e+0 )
285 parameter ( eight = 8.0e+0, sevten = 17.0e+0 )
287 parameter ( cone = ( 1.0e+0, 0.0e+0 ),
288 $ czero = ( 0.0e+0, 0.0e+0 ) )
292 INTEGER IMAX, ITEMP, II, J, JB, JJ, JMAX, K, KK, KKW,
294 REAL ABSAKK, ALPHA, COLMAX, STEMP, R1, ROWMAX, T,
296 COMPLEX D11, D21, D22, Z
302 EXTERNAL lsame, icamax, slamch
308 INTRINSIC abs, conjg, aimag, max, min,
REAL, SQRT
314 cabs1( z ) = abs(
REAL( Z ) ) + abs( AIMAG( z ) )
322 alpha = ( one+sqrt( sevten ) ) / eight
326 sfmin = slamch(
'S' )
328 IF( lsame( uplo,
'U' ) )
THEN
350 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
359 $
CALL ccopy( k-1, a( 1, k ), 1, w( 1, kw ), 1 )
360 w( k, kw ) =
REAL( A( K, K ) )
362 CALL cgemv(
'No transpose', k, n-k, -cone, a( 1, k+1 ), lda,
363 $ w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
364 w( k, kw ) =
REAL( W( K, KW ) )
370 absakk = abs(
REAL( W( K, KW ) ) )
377 imax = icamax( k-1, w( 1, kw ), 1 )
378 colmax = cabs1( w( imax, kw ) )
383 IF( max( absakk, colmax ).EQ.zero )
THEN
390 a( k, k ) =
REAL( W( K, KW ) )
392 $
CALL ccopy( k-1, w( 1, kw ), 1, a( 1, k ), 1 )
408 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
428 $
CALL ccopy( imax-1, a( 1, imax ), 1, w( 1, kw-1 ),
430 w( imax, kw-1 ) =
REAL( A( IMAX, IMAX ) )
432 CALL ccopy( k-imax, a( imax, imax+1 ), lda,
433 $ w( imax+1, kw-1 ), 1 )
434 CALL clacgv( k-imax, w( imax+1, kw-1 ), 1 )
437 CALL cgemv(
'No transpose', k, n-k, -cone,
438 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
439 $ cone, w( 1, kw-1 ), 1 )
440 w( imax, kw-1 ) =
REAL( W( IMAX, KW-1 ) )
448 jmax = imax + icamax( k-imax, w( imax+1, kw-1 ),
450 rowmax = cabs1( w( jmax, kw-1 ) )
456 itemp = icamax( imax-1, w( 1, kw-1 ), 1 )
457 stemp = cabs1( w( itemp, kw-1 ) )
458 IF( stemp.GT.rowmax )
THEN
469 IF( .NOT.( abs(
REAL( W( IMAX,KW-1 ) ) )
470 $ .LT.alpha*rowmax ) )
THEN
479 CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
487 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
508 CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
515 IF( .NOT.done )
GOTO 12
534 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
541 a( p, p ) =
REAL( A( K, K ) )
542 CALL ccopy( k-1-p, a( p+1, k ), 1, a( p, p+1 ),
544 CALL clacgv( k-1-p, a( p, p+1 ), lda )
546 $
CALL ccopy( p-1, a( 1, k ), 1, a( 1, p ), 1 )
554 $
CALL cswap( n-k, a( k, k+1 ), lda, a( p, k+1 ),
556 CALL cswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ),
570 a( kp, kp ) =
REAL( A( KK, KK ) )
571 CALL ccopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
573 CALL clacgv( kk-1-kp, a( kp, kp+1 ), lda )
575 $
CALL ccopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
583 $
CALL cswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
585 CALL cswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
589 IF( kstep.EQ.1 )
THEN
607 CALL ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
616 t =
REAL( A( K, K ) )
617 IF( abs( t ).GE.sfmin )
THEN
619 CALL csscal( k-1, r1, a( 1, k ), 1 )
622 a( ii, k ) = a( ii, k ) / t
628 CALL clacgv( k-1, w( 1, kw ), 1 )
701 d11 = w( k, kw ) / conjg( d21 )
702 d22 = w( k-1, kw-1 ) / d21
703 t = one / (
REAL( d11*d22 )-ONE )
710 a( j, k-1 ) = t*( ( d11*w( j, kw-1 )-w( j, kw ) ) /
712 a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
721 a( k-1, k-1 ) = w( k-1, kw-1 )
723 a( k, k ) = w( k, kw )
724 e( k ) = w( k-1, kw )
729 CALL clacgv( k-1, w( 1, kw ), 1 )
730 CALL clacgv( k-2, w( 1, kw-1 ), 1 )
740 IF( kstep.EQ.1 )
THEN
761 DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
762 jb = min( nb, k-j+1 )
766 DO 40 jj = j, j + jb - 1
767 a( jj, jj ) =
REAL( A( JJ, JJ ) )
768 CALL cgemv(
'No transpose', jj-j+1, n-k, -cone,
769 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
771 a( jj, jj ) =
REAL( A( JJ, JJ ) )
777 $
CALL cgemm(
'No transpose',
'Transpose', j-1, jb, n-k,
778 $ -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
779 $ cone, a( 1, j ), lda )
803 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
811 w( k, k ) =
REAL( A( K, K ) )
813 $
CALL ccopy( n-k, a( k+1, k ), 1, w( k+1, k ), 1 )
815 CALL cgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
816 $ lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
817 w( k, k ) =
REAL( W( K, K ) )
823 absakk = abs(
REAL( W( K, K ) ) )
830 imax = k + icamax( n-k, w( k+1, k ), 1 )
831 colmax = cabs1( w( imax, k ) )
836 IF( max( absakk, colmax ).EQ.zero )
THEN
843 a( k, k ) =
REAL( W( K, K ) )
845 $
CALL ccopy( n-k, w( k+1, k ), 1, a( k+1, k ), 1 )
862 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
881 CALL ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
882 CALL clacgv( imax-k, w( k, k+1 ), 1 )
883 w( imax, k+1 ) =
REAL( A( IMAX, IMAX ) )
886 $
CALL ccopy( n-imax, a( imax+1, imax ), 1,
887 $ w( imax+1, k+1 ), 1 )
890 CALL cgemv(
'No transpose', n-k+1, k-1, -cone,
891 $ a( k, 1 ), lda, w( imax, 1 ), ldw,
892 $ cone, w( k, k+1 ), 1 )
893 w( imax, k+1 ) =
REAL( W( IMAX, K+1 ) )
901 jmax = k - 1 + icamax( imax-k, w( k, k+1 ), 1 )
902 rowmax = cabs1( w( jmax, k+1 ) )
908 itemp = imax + icamax( n-imax, w( imax+1, k+1 ), 1)
909 stemp = cabs1( w( itemp, k+1 ) )
910 IF( stemp.GT.rowmax )
THEN
921 IF( .NOT.( abs(
REAL( W( IMAX,K+1 ) ) )
922 $ .LT.alpha*rowmax ) )
THEN
931 CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
939 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
960 CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
967 IF( .NOT.done )
GOTO 72
982 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
989 a( p, p ) =
REAL( A( K, K ) )
990 CALL ccopy( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
991 CALL clacgv( p-k-1, a( p, k+1 ), lda )
993 $
CALL ccopy( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
1001 $
CALL cswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
1002 CALL cswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
1015 a( kp, kp ) =
REAL( A( KK, KK ) )
1016 CALL ccopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
1018 CALL clacgv( kp-kk-1, a( kp, kk+1 ), lda )
1020 $
CALL ccopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
1028 $
CALL cswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
1029 CALL cswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
1032 IF( kstep.EQ.1 )
THEN
1050 CALL ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
1059 t =
REAL( A( K, K ) )
1060 IF( abs( t ).GE.sfmin )
THEN
1062 CALL csscal( n-k, r1, a( k+1, k ), 1 )
1065 a( ii, k ) = a( ii, k ) / t
1071 CALL clacgv( n-k, w( k+1, k ), 1 )
1144 d11 = w( k+1, k+1 ) / d21
1145 d22 = w( k, k ) / conjg( d21 )
1146 t = one / (
REAL( d11*d22 )-ONE )
1153 a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
1155 a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
1164 a( k, k ) = w( k, k )
1166 a( k+1, k+1 ) = w( k+1, k+1 )
1167 e( k ) = w( k+1, k )
1172 CALL clacgv( n-k, w( k+1, k ), 1 )
1173 CALL clacgv( n-k-1, w( k+2, k+1 ), 1 )
1183 IF( kstep.EQ.1 )
THEN
1205 jb = min( nb, n-j+1 )
1209 DO 100 jj = j, j + jb - 1
1210 a( jj, jj ) =
REAL( A( JJ, JJ ) )
1211 CALL cgemv(
'No transpose', j+jb-jj, k-1, -cone,
1212 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
1214 a( jj, jj ) =
REAL( A( JJ, JJ ) )
1220 $
CALL cgemm(
'No transpose',
'Transpose', n-j-jb+1, jb,
1221 $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
1222 $ ldw, cone, a( j+jb, j ), lda )
subroutine clahef_rk(UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW, INFO)
CLAHEF_RK computes a partial factorization of a complex Hermitian indefinite matrix using bounded Bun...
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 csscal(N, SA, CX, INCX)
CSSCAL