178 SUBROUTINE zlahef( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
187 INTEGER info, kb, lda, ldw, n, nb
191 COMPLEX*16 a( lda, * ), w( ldw, * )
197 DOUBLE PRECISION zero, one
198 parameter( zero = 0.0d+0, one = 1.0d+0 )
200 parameter( cone = ( 1.0d+0, 0.0d+0 ) )
201 DOUBLE PRECISION eight, sevten
202 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
205 INTEGER imax,
j, jb, jj, jmax, jp, k, kk, kkw, kp,
207 DOUBLE PRECISION absakk, alpha, colmax, r1, rowmax, t
208 COMPLEX*16 d11, d21, d22, z
219 INTRINSIC abs, dble, dconjg, dimag, max, min, sqrt
222 DOUBLE PRECISION cabs1
225 cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
233 alpha = ( one+sqrt( sevten ) ) / eight
235 IF(
lsame( uplo,
'U' ) )
THEN
251 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
258 CALL
zcopy( k-1, a( 1, k ), 1, w( 1, kw ), 1 )
259 w( k, kw ) = dble( a( k, k ) )
261 CALL
zgemv(
'No transpose', k, n-k, -cone, a( 1, k+1 ), lda,
262 $ w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
263 w( k, kw ) = dble( w( k, kw ) )
269 absakk = abs( dble( w( k, kw ) ) )
276 imax =
izamax( k-1, w( 1, kw ), 1 )
277 colmax = cabs1( w( imax, kw ) )
282 IF( max( absakk, colmax ).EQ.zero )
THEN
289 a( k, k ) = dble( a( k, k ) )
297 IF( absakk.GE.alpha*colmax )
THEN
309 CALL
zcopy( imax-1, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
310 w( imax, kw-1 ) = dble( a( imax, imax ) )
311 CALL
zcopy( k-imax, a( imax, imax+1 ), lda,
312 $ w( imax+1, kw-1 ), 1 )
313 CALL
zlacgv( k-imax, w( imax+1, kw-1 ), 1 )
315 CALL
zgemv(
'No transpose', k, n-k, -cone,
316 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
317 $ cone, w( 1, kw-1 ), 1 )
318 w( imax, kw-1 ) = dble( w( imax, kw-1 ) )
325 jmax = imax +
izamax( k-imax, w( imax+1, kw-1 ), 1 )
326 rowmax = cabs1( w( jmax, kw-1 ) )
328 jmax =
izamax( imax-1, w( 1, kw-1 ), 1 )
329 rowmax = max( rowmax, cabs1( w( jmax, kw-1 ) ) )
333 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
340 ELSE IF( abs( dble( w( imax, kw-1 ) ) ).GE.alpha*rowmax )
350 CALL
zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
389 a( kp, kp ) = dble( a( kk, kk ) )
390 CALL
zcopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
392 CALL
zlacgv( kk-1-kp, a( kp, kp+1 ), lda )
394 $ CALL
zcopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
402 $ CALL
zswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
404 CALL
zswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
408 IF( kstep.EQ.1 )
THEN
426 CALL
zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
433 r1 = one / dble( a( k, k ) )
434 CALL
zdscal( k-1, r1, a( 1, k ), 1 )
438 CALL
zlacgv( k-1, w( 1, kw ), 1 )
503 d11 = w( k, kw ) / dconjg( d21 )
504 d22 = w( k-1, kw-1 ) / d21
505 t = one / ( dble( d11*d22 )-one )
513 a(
j, k-1 ) = d21*( d11*w(
j, kw-1 )-w(
j, kw ) )
514 a(
j, k ) = dconjg( d21 )*
515 $ ( d22*w(
j, kw )-w(
j, kw-1 ) )
521 a( k-1, k-1 ) = w( k-1, kw-1 )
522 a( k-1, k ) = w( k-1, kw )
523 a( k, k ) = w( k, kw )
527 CALL
zlacgv( k-1, w( 1, kw ), 1 )
528 CALL
zlacgv( k-2, w( 1, kw-1 ), 1 )
536 IF( kstep.EQ.1 )
THEN
557 DO 50
j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
558 jb = min( nb, k-
j+1 )
562 DO 40 jj =
j,
j + jb - 1
563 a( jj, jj ) = dble( a( jj, jj ) )
564 CALL
zgemv(
'No transpose', jj-
j+1, n-k, -cone,
565 $ a(
j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
567 a( jj, jj ) = dble( a( jj, jj ) )
572 CALL
zgemm(
'No transpose',
'Transpose',
j-1, jb, n-k,
573 $ -cone, a( 1, k+1 ), lda, w(
j, kw+1 ), ldw,
574 $ cone, a( 1,
j ), lda )
597 IF( jp.NE.jj .AND.
j.LE.n )
598 $ CALL
zswap( n-
j+1, a( jp,
j ), lda, a( jj,
j ), lda )
619 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
626 w( k, k ) = dble( a( k, k ) )
628 $ CALL
zcopy( n-k, a( k+1, k ), 1, w( k+1, k ), 1 )
629 CALL
zgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ), lda,
630 $ w( k, 1 ), ldw, cone, w( k, k ), 1 )
631 w( k, k ) = dble( w( k, k ) )
636 absakk = abs( dble( w( k, k ) ) )
643 imax = k +
izamax( n-k, w( k+1, k ), 1 )
644 colmax = cabs1( w( imax, k ) )
649 IF( max( absakk, colmax ).EQ.zero )
THEN
656 a( k, k ) = dble( a( k, k ) )
664 IF( absakk.GE.alpha*colmax )
THEN
676 CALL
zcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1 )
677 CALL
zlacgv( imax-k, w( k, k+1 ), 1 )
678 w( imax, k+1 ) = dble( a( imax, imax ) )
680 $ CALL
zcopy( n-imax, a( imax+1, imax ), 1,
681 $ w( imax+1, k+1 ), 1 )
682 CALL
zgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
683 $ lda, w( imax, 1 ), ldw, cone, w( k, k+1 ),
685 w( imax, k+1 ) = dble( w( imax, k+1 ) )
691 jmax = k - 1 +
izamax( imax-k, w( k, k+1 ), 1 )
692 rowmax = cabs1( w( jmax, k+1 ) )
694 jmax = imax +
izamax( n-imax, w( imax+1, k+1 ), 1 )
695 rowmax = max( rowmax, cabs1( w( jmax, k+1 ) ) )
699 IF( absakk.GE.alpha*colmax*( colmax / rowmax ) )
THEN
706 ELSE IF( abs( dble( w( imax, k+1 ) ) ).GE.alpha*rowmax )
716 CALL
zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
751 a( kp, kp ) = dble( a( kk, kk ) )
752 CALL
zcopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
754 CALL
zlacgv( kp-kk-1, a( kp, kk+1 ), lda )
756 $ CALL
zcopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
764 $ CALL
zswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
765 CALL
zswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
768 IF( kstep.EQ.1 )
THEN
786 CALL
zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
793 r1 = one / dble( a( k, k ) )
794 CALL
zdscal( n-k, r1, a( k+1, k ), 1 )
798 CALL
zlacgv( n-k, w( k+1, k ), 1 )
863 d11 = w( k+1, k+1 ) / d21
864 d22 = w( k, k ) / dconjg( d21 )
865 t = one / ( dble( d11*d22 )-one )
873 a(
j, k ) = dconjg( d21 )*
874 $ ( d11*w(
j, k )-w(
j, k+1 ) )
875 a(
j, k+1 ) = d21*( d22*w(
j, k+1 )-w(
j, k ) )
881 a( k, k ) = w( k, k )
882 a( k+1, k ) = w( k+1, k )
883 a( k+1, k+1 ) = w( k+1, k+1 )
887 CALL
zlacgv( n-k, w( k+1, k ), 1 )
888 CALL
zlacgv( n-k-1, w( k+2, k+1 ), 1 )
896 IF( kstep.EQ.1 )
THEN
918 jb = min( nb, n-
j+1 )
922 DO 100 jj =
j,
j + jb - 1
923 a( jj, jj ) = dble( a( jj, jj ) )
924 CALL
zgemv(
'No transpose',
j+jb-jj, k-1, -cone,
925 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
927 a( jj, jj ) = dble( a( jj, jj ) )
933 $ CALL
zgemm(
'No transpose',
'Transpose', n-
j-jb+1, jb,
934 $ k-1, -cone, a(
j+jb, 1 ), lda, w(
j, 1 ),
935 $ ldw, cone, a(
j+jb,
j ), lda )
958 IF( jp.NE.jj .AND.
j.GE.1 )
959 $ CALL
zswap(
j, a( jp, 1 ), lda, a( jj, 1 ), lda )