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 zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zlahef_rook(UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY