184 SUBROUTINE zlahef_rook( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW,
194 INTEGER INFO, KB, LDA, LDW, N, NB
198 COMPLEX*16 A( lda, * ), W( ldw, * )
204 DOUBLE PRECISION ZERO, ONE
205 parameter ( zero = 0.0d+0, one = 1.0d+0 )
207 parameter ( cone = ( 1.0d+0, 0.0d+0 ) )
208 DOUBLE PRECISION EIGHT, SEVTEN
209 parameter ( eight = 8.0d+0, sevten = 17.0d+0 )
213 INTEGER IMAX, ITEMP, II, J, JB, JJ, JMAX, JP1, JP2, K,
214 $ kk, kkw, kp, kstep, kw, p
215 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, DTEMP, R1, ROWMAX, T,
217 COMPLEX*16 D11, D21, D22, Z
222 DOUBLE PRECISION DLAMCH
223 EXTERNAL lsame, izamax, dlamch
229 INTRINSIC abs, dble, dconjg, dimag, max, min, sqrt
232 DOUBLE PRECISION CABS1
235 cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
243 alpha = ( one+sqrt( sevten ) ) / eight
247 sfmin = dlamch(
'S' )
249 IF( lsame( uplo,
'U' ) )
THEN
266 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
275 $
CALL zcopy( k-1, a( 1, k ), 1, w( 1, kw ), 1 )
276 w( k, kw ) = dble( a( k, k ) )
278 CALL zgemv(
'No transpose', k, n-k, -cone, a( 1, k+1 ), lda,
279 $ w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
280 w( k, kw ) = dble( w( k, kw ) )
286 absakk = abs( dble( w( k, kw ) ) )
293 imax = izamax( k-1, w( 1, kw ), 1 )
294 colmax = cabs1( w( imax, kw ) )
299 IF( max( absakk, colmax ).EQ.zero )
THEN
306 a( k, k ) = dble( w( k, kw ) )
308 $
CALL zcopy( k-1, w( 1, kw ), 1, a( 1, k ), 1 )
318 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
338 $
CALL zcopy( imax-1, a( 1, imax ), 1, w( 1, kw-1 ),
340 w( imax, kw-1 ) = dble( a( imax, imax ) )
342 CALL zcopy( k-imax, a( imax, imax+1 ), lda,
343 $ w( imax+1, kw-1 ), 1 )
344 CALL zlacgv( k-imax, w( imax+1, kw-1 ), 1 )
347 CALL zgemv(
'No transpose', k, n-k, -cone,
348 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
349 $ cone, w( 1, kw-1 ), 1 )
350 w( imax, kw-1 ) = dble( w( imax, kw-1 ) )
358 jmax = imax + izamax( k-imax, w( imax+1, kw-1 ),
360 rowmax = cabs1( w( jmax, kw-1 ) )
366 itemp = izamax( imax-1, w( 1, kw-1 ), 1 )
367 dtemp = cabs1( w( itemp, kw-1 ) )
368 IF( dtemp.GT.rowmax )
THEN
379 IF( .NOT.( abs( dble( w( imax,kw-1 ) ) )
380 $ .LT.alpha*rowmax ) )
THEN
389 CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
397 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
418 CALL zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
425 IF( .NOT.done )
GOTO 12
444 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
451 a( p, p ) = dble( a( k, k ) )
452 CALL zcopy( k-1-p, a( p+1, k ), 1, a( p, p+1 ),
454 CALL zlacgv( k-1-p, a( p, p+1 ), lda )
456 $
CALL zcopy( p-1, a( 1, k ), 1, a( 1, p ), 1 )
464 $
CALL zswap( n-k, a( k, k+1 ), lda, a( p, k+1 ),
466 CALL zswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ),
480 a( kp, kp ) = dble( a( kk, kk ) )
481 CALL zcopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
483 CALL zlacgv( kk-1-kp, a( kp, kp+1 ), lda )
485 $
CALL zcopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
493 $
CALL zswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
495 CALL zswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
499 IF( kstep.EQ.1 )
THEN
517 CALL zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
526 t = dble( a( k, k ) )
527 IF( abs( t ).GE.sfmin )
THEN
529 CALL zdscal( k-1, r1, a( 1, k ), 1 )
532 a( ii, k ) = a( ii, k ) / t
538 CALL zlacgv( k-1, w( 1, kw ), 1 )
606 d11 = w( k, kw ) / dconjg( d21 )
607 d22 = w( k-1, kw-1 ) / d21
608 t = one / ( dble( d11*d22 )-one )
615 a( j, k-1 ) = t*( ( d11*w( j, kw-1 )-w( j, kw ) ) /
617 a( j, k ) = t*( ( d22*w( j, kw )-w( j, kw-1 ) ) /
624 a( k-1, k-1 ) = w( k-1, kw-1 )
625 a( k-1, k ) = w( k-1, kw )
626 a( k, k ) = w( k, kw )
630 CALL zlacgv( k-1, w( 1, kw ), 1 )
631 CALL zlacgv( k-2, w( 1, kw-1 ), 1 )
639 IF( kstep.EQ.1 )
THEN
660 DO 50 j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
661 jb = min( nb, k-j+1 )
665 DO 40 jj = j, j + jb - 1
666 a( jj, jj ) = dble( a( jj, jj ) )
667 CALL zgemv(
'No transpose', jj-j+1, n-k, -cone,
668 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
670 a( jj, jj ) = dble( a( jj, jj ) )
676 $
CALL zgemm(
'No transpose',
'Transpose', j-1, jb, n-k,
677 $ -cone, a( 1, k+1 ), lda, w( j, kw+1 ), ldw,
678 $ cone, a( 1, j ), lda )
705 IF( jp2.NE.jj .AND. j.LE.n )
706 $
CALL zswap( n-j+1, a( jp2, j ), lda, a( jj, j ), lda )
708 IF( kstep.EQ.2 .AND. jp1.NE.jj .AND. j.LE.n )
709 $
CALL zswap( n-j+1, a( jp1, j ), lda, a( jj, j ), lda )
730 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
738 w( k, k ) = dble( a( k, k ) )
740 $
CALL zcopy( n-k, a( k+1, k ), 1, w( k+1, k ), 1 )
742 CALL zgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
743 $ lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
744 w( k, k ) = dble( w( k, k ) )
750 absakk = abs( dble( w( k, k ) ) )
757 imax = k + izamax( n-k, w( k+1, k ), 1 )
758 colmax = cabs1( w( imax, k ) )
763 IF( max( absakk, colmax ).EQ.zero )
THEN
770 a( k, k ) = dble( w( k, k ) )
772 $
CALL zcopy( n-k, w( k+1, k ), 1, a( k+1, k ), 1 )
783 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
802 CALL zcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
803 CALL zlacgv( imax-k, w( k, k+1 ), 1 )
804 w( imax, k+1 ) = dble( a( imax, imax ) )
807 $
CALL zcopy( n-imax, a( imax+1, imax ), 1,
808 $ w( imax+1, k+1 ), 1 )
811 CALL zgemv(
'No transpose', n-k+1, k-1, -cone,
812 $ a( k, 1 ), lda, w( imax, 1 ), ldw,
813 $ cone, w( k, k+1 ), 1 )
814 w( imax, k+1 ) = dble( w( imax, k+1 ) )
822 jmax = k - 1 + izamax( imax-k, w( k, k+1 ), 1 )
823 rowmax = cabs1( w( jmax, k+1 ) )
829 itemp = imax + izamax( n-imax, w( imax+1, k+1 ), 1)
830 dtemp = cabs1( w( itemp, k+1 ) )
831 IF( dtemp.GT.rowmax )
THEN
842 IF( .NOT.( abs( dble( w( imax,k+1 ) ) )
843 $ .LT.alpha*rowmax ) )
THEN
852 CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
860 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
881 CALL zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
888 IF( .NOT.done )
GOTO 72
903 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
910 a( p, p ) = dble( a( k, k ) )
911 CALL zcopy( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
912 CALL zlacgv( p-k-1, a( p, k+1 ), lda )
914 $
CALL zcopy( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
922 $
CALL zswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
923 CALL zswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
936 a( kp, kp ) = dble( a( kk, kk ) )
937 CALL zcopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
939 CALL zlacgv( kp-kk-1, a( kp, kk+1 ), lda )
941 $
CALL zcopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
949 $
CALL zswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
950 CALL zswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
953 IF( kstep.EQ.1 )
THEN
971 CALL zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
980 t = dble( a( k, k ) )
981 IF( abs( t ).GE.sfmin )
THEN
983 CALL zdscal( n-k, r1, a( k+1, k ), 1 )
986 a( ii, k ) = a( ii, k ) / t
992 CALL zlacgv( n-k, w( k+1, k ), 1 )
1060 d11 = w( k+1, k+1 ) / d21
1061 d22 = w( k, k ) / dconjg( d21 )
1062 t = one / ( dble( d11*d22 )-one )
1069 a( j, k ) = t*( ( d11*w( j, k )-w( j, k+1 ) ) /
1071 a( j, k+1 ) = t*( ( d22*w( j, k+1 )-w( j, k ) ) /
1078 a( k, k ) = w( k, k )
1079 a( k+1, k ) = w( k+1, k )
1080 a( k+1, k+1 ) = w( k+1, k+1 )
1084 CALL zlacgv( n-k, w( k+1, k ), 1 )
1085 CALL zlacgv( n-k-1, w( k+2, k+1 ), 1 )
1093 IF( kstep.EQ.1 )
THEN
1115 jb = min( nb, n-j+1 )
1119 DO 100 jj = j, j + jb - 1
1120 a( jj, jj ) = dble( a( jj, jj ) )
1121 CALL zgemv(
'No transpose', j+jb-jj, k-1, -cone,
1122 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
1124 a( jj, jj ) = dble( a( jj, jj ) )
1130 $
CALL zgemm(
'No transpose',
'Transpose', n-j-jb+1, jb,
1131 $ k-1, -cone, a( j+jb, 1 ), lda, w( j, 1 ),
1132 $ ldw, cone, a( j+jb, j ), lda )
1159 IF( jp2.NE.jj .AND. j.GE.1 )
1160 $
CALL zswap( j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
1162 IF( kstep.EQ.2 .AND. jp1.NE.jj .AND. j.GE.1 )
1163 $
CALL zswap( j, a( jp1, 1 ), lda, a( jj, 1 ), lda )
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
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
subroutine zlahef_rook(UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)