262 SUBROUTINE zlahef_rk( UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW,
272 INTEGER INFO, KB, LDA, LDW, N, NB
276 COMPLEX*16 A( lda, * ), W( ldw, * ), E( * )
282 DOUBLE PRECISION ZERO, ONE
283 parameter ( zero = 0.0d+0, one = 1.0d+0 )
285 parameter ( cone = ( 1.0d+0, 0.0d+0 ) )
286 DOUBLE PRECISION EIGHT, SEVTEN
287 parameter ( eight = 8.0d+0, sevten = 17.0d+0 )
289 parameter ( czero = ( 0.0d+0, 0.0d+0 ) )
293 INTEGER IMAX, ITEMP, II, J, JB, JJ, JMAX, K, KK, KKW,
295 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, DTEMP, R1, ROWMAX, T,
297 COMPLEX*16 D11, D21, D22, Z
302 DOUBLE PRECISION DLAMCH
303 EXTERNAL lsame, izamax, dlamch
309 INTRINSIC abs, dble, dconjg, dimag, max, min, sqrt
312 DOUBLE PRECISION CABS1
315 cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
323 alpha = ( one+sqrt( sevten ) ) / eight
327 sfmin = dlamch(
'S' )
329 IF( lsame( uplo,
'U' ) )
THEN
350 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
359 $
CALL zcopy( k-1, a( 1, k ), 1, w( 1, kw ), 1 )
360 w( k, kw ) = dble( a( k, k ) )
362 CALL zgemv(
'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 ) = dble( w( k, kw ) )
370 absakk = abs( dble( w( k, kw ) ) )
377 imax = izamax( k-1, w( 1, kw ), 1 )
378 colmax = cabs1( w( imax, kw ) )
383 IF( max( absakk, colmax ).EQ.zero )
THEN
390 a( k, k ) = dble( w( k, kw ) )
392 $
CALL zcopy( k-1, w( 1, kw ), 1, a( 1, k ), 1 )
408 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
428 $
CALL zcopy( imax-1, a( 1, imax ), 1, w( 1, kw-1 ),
430 w( imax, kw-1 ) = dble( a( imax, imax ) )
432 CALL zcopy( k-imax, a( imax, imax+1 ), lda,
433 $ w( imax+1, kw-1 ), 1 )
434 CALL zlacgv( k-imax, w( imax+1, kw-1 ), 1 )
437 CALL zgemv(
'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 ) = dble( w( imax, kw-1 ) )
448 jmax = imax + izamax( k-imax, w( imax+1, kw-1 ),
450 rowmax = cabs1( w( jmax, kw-1 ) )
456 itemp = izamax( imax-1, w( 1, kw-1 ), 1 )
457 dtemp = cabs1( w( itemp, kw-1 ) )
458 IF( dtemp.GT.rowmax )
THEN
469 IF( .NOT.( abs( dble( w( imax,kw-1 ) ) )
470 $ .LT.alpha*rowmax ) )
THEN
479 CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
487 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
508 CALL zcopy( 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 ) = dble( a( k, k ) )
542 CALL zcopy( k-1-p, a( p+1, k ), 1, a( p, p+1 ),
544 CALL zlacgv( k-1-p, a( p, p+1 ), lda )
546 $
CALL zcopy( p-1, a( 1, k ), 1, a( 1, p ), 1 )
554 $
CALL zswap( n-k, a( k, k+1 ), lda, a( p, k+1 ),
556 CALL zswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ),
570 a( kp, kp ) = dble( a( kk, kk ) )
571 CALL zcopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
573 CALL zlacgv( kk-1-kp, a( kp, kp+1 ), lda )
575 $
CALL zcopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
583 $
CALL zswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
585 CALL zswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
589 IF( kstep.EQ.1 )
THEN
607 CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
616 t = dble( a( k, k ) )
617 IF( abs( t ).GE.sfmin )
THEN
619 CALL zdscal( k-1, r1, a( 1, k ), 1 )
622 a( ii, k ) = a( ii, k ) / t
628 CALL zlacgv( k-1, w( 1, kw ), 1 )
701 d11 = w( k, kw ) / dconjg( d21 )
702 d22 = w( k-1, kw-1 ) / d21
703 t = one / ( dble( 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 zlacgv( k-1, w( 1, kw ), 1 )
730 CALL zlacgv( 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 ) = dble( a( jj, jj ) )
768 CALL zgemv(
'No transpose', jj-j+1, n-k, -cone,
769 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
771 a( jj, jj ) = dble( a( jj, jj ) )
777 $
CALL zgemm(
'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 ) = dble( a( k, k ) )
813 $
CALL zcopy( n-k, a( k+1, k ), 1, w( k+1, k ), 1 )
815 CALL zgemv(
'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 ) = dble( w( k, k ) )
823 absakk = abs( dble( w( k, k ) ) )
830 imax = k + izamax( 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 ) = dble( w( k, k ) )
845 $
CALL zcopy( n-k, w( k+1, k ), 1, a( k+1, k ), 1 )
862 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
881 CALL zcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
882 CALL zlacgv( imax-k, w( k, k+1 ), 1 )
883 w( imax, k+1 ) = dble( a( imax, imax ) )
886 $
CALL zcopy( n-imax, a( imax+1, imax ), 1,
887 $ w( imax+1, k+1 ), 1 )
890 CALL zgemv(
'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 ) = dble( w( imax, k+1 ) )
901 jmax = k - 1 + izamax( imax-k, w( k, k+1 ), 1 )
902 rowmax = cabs1( w( jmax, k+1 ) )
908 itemp = imax + izamax( n-imax, w( imax+1, k+1 ), 1)
909 dtemp = cabs1( w( itemp, k+1 ) )
910 IF( dtemp.GT.rowmax )
THEN
921 IF( .NOT.( abs( dble( w( imax,k+1 ) ) )
922 $ .LT.alpha*rowmax ) )
THEN
931 CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
939 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
960 CALL zcopy( 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 ) = dble( a( k, k ) )
990 CALL zcopy( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
991 CALL zlacgv( p-k-1, a( p, k+1 ), lda )
993 $
CALL zcopy( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
1001 $
CALL zswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
1002 CALL zswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
1015 a( kp, kp ) = dble( a( kk, kk ) )
1016 CALL zcopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
1018 CALL zlacgv( kp-kk-1, a( kp, kk+1 ), lda )
1020 $
CALL zcopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
1028 $
CALL zswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
1029 CALL zswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
1032 IF( kstep.EQ.1 )
THEN
1050 CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
1059 t = dble( a( k, k ) )
1060 IF( abs( t ).GE.sfmin )
THEN
1062 CALL zdscal( n-k, r1, a( k+1, k ), 1 )
1065 a( ii, k ) = a( ii, k ) / t
1071 CALL zlacgv( n-k, w( k+1, k ), 1 )
1144 d11 = w( k+1, k+1 ) / d21
1145 d22 = w( k, k ) / dconjg( d21 )
1146 t = one / ( dble( 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 zlacgv( n-k, w( k+1, k ), 1 )
1173 CALL zlacgv( 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 ) = dble( a( jj, jj ) )
1211 CALL zgemv(
'No transpose', j+jb-jj, k-1, -cone,
1212 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
1214 a( jj, jj ) = dble( a( jj, jj ) )
1220 $
CALL zgemm(
'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 zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
subroutine zlahef_rk(UPLO, N, NB, KB, A, LDA, E, IPIV, W, LDW, INFO)
ZLAHEF_RK computes a partial factorization of a complex Hermitian indefinite matrix using bounded Bun...
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL