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 )
206 DOUBLE PRECISION eight, sevten
207 parameter( eight = 8.0d+0, sevten = 17.0d+0 )
208 COMPLEX*16 cone, czero
209 parameter( cone = ( 1.0d+0, 0.0d+0 ),
210 $ czero = ( 0.0d+0, 0.0d+0 ) )
214 INTEGER imax, itemp,
j, jb, jj, jmax, jp1, jp2, k, kk,
215 $ kw, kkw, kp, kstep, p, ii
216 DOUBLE PRECISION absakk, alpha, colmax, rowmax, dtemp, sfmin
217 COMPLEX*16 d11, d12, d21, d22, r1, t, z
229 INTRINSIC abs, max, min, sqrt, dimag, dble
232 DOUBLE PRECISION cabs1
235 cabs1( z ) = abs( dble( z ) ) + abs( dimag( z ) )
243 alpha = ( one+sqrt( sevten ) ) / eight
249 IF(
lsame( uplo,
'U' ) )
THEN
266 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
274 CALL
zcopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
276 $ CALL
zgemv(
'No transpose', k, n-k, -cone, a( 1, k+1 ),
277 $ lda, w( k, kw+1 ), ldw, cone, w( 1, kw ), 1 )
282 absakk = cabs1( w( k, kw ) )
289 imax =
izamax( k-1, w( 1, kw ), 1 )
290 colmax = cabs1( w( imax, kw ) )
295 IF( max( absakk, colmax ).EQ.zero )
THEN
302 CALL
zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
312 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
331 CALL
zcopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
332 CALL
zcopy( k-imax, a( imax, imax+1 ), lda,
333 $ w( imax+1, kw-1 ), 1 )
336 $ CALL
zgemv(
'No transpose', k, n-k, -cone,
337 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
338 $ cone, w( 1, kw-1 ), 1 )
345 jmax = imax +
izamax( k-imax, w( imax+1, kw-1 ),
347 rowmax = cabs1( w( jmax, kw-1 ) )
353 itemp =
izamax( imax-1, w( 1, kw-1 ), 1 )
354 dtemp = cabs1( w( itemp, kw-1 ) )
355 IF( dtemp.GT.rowmax )
THEN
365 IF( .NOT.(cabs1( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
375 CALL
zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
382 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
401 CALL
zcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
407 IF( .NOT. done ) goto 12
419 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
423 CALL
zcopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
424 CALL
zcopy( p, a( 1, k ), 1, a( 1, p ), 1 )
429 CALL
zswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
430 CALL
zswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
439 a( kp, k ) = a( kk, k )
440 CALL
zcopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
442 CALL
zcopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
447 CALL
zswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
448 CALL
zswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
452 IF( kstep.EQ.1 )
THEN
462 CALL
zcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
464 IF( cabs1( a( k, k ) ).GE.sfmin )
THEN
465 r1 = cone / a( k, k )
466 CALL
zscal( k-1, r1, a( 1, k ), 1 )
467 ELSE IF( a( k, k ).NE.czero )
THEN
469 a( ii, k ) = a( ii, k ) / a( k, k )
489 d11 = w( k, kw ) / d12
490 d22 = w( k-1, kw-1 ) / d12
491 t = cone / ( d11*d22-cone )
493 a(
j, k-1 ) = t*( (d11*w(
j, kw-1 )-w(
j, kw ) ) /
495 a(
j, k ) = t*( ( d22*w(
j, kw )-w(
j, kw-1 ) ) /
502 a( k-1, k-1 ) = w( k-1, kw-1 )
503 a( k-1, k ) = w( k-1, kw )
504 a( k, k ) = w( k, kw )
510 IF( kstep.EQ.1 )
THEN
530 DO 50
j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
531 jb = min( nb, k-
j+1 )
535 DO 40 jj =
j,
j + jb - 1
536 CALL
zgemv(
'No transpose', jj-
j+1, n-k, -cone,
537 $ a(
j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
544 $ CALL
zgemm(
'No transpose',
'Transpose',
j-1, jb,
545 $ n-k, -cone, a( 1, k+1 ), lda, w(
j, kw+1 ), ldw,
546 $ cone, a( 1,
j ), lda )
567 IF( jp2.NE.jj .AND.
j.LE.n )
568 $ CALL
zswap( n-
j+1, a( jp2,
j ), lda, a( jj,
j ), lda )
570 IF( jp1.NE.jj .AND. kstep.EQ.2 )
571 $ CALL
zswap( n-
j+1, a( jp1,
j ), lda, a( jj,
j ), lda )
592 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
600 CALL
zcopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
602 $ CALL
zgemv(
'No transpose', n-k+1, k-1, -cone, a( k, 1 ),
603 $ lda, w( k, 1 ), ldw, cone, w( k, k ), 1 )
608 absakk = cabs1( w( k, k ) )
615 imax = k +
izamax( n-k, w( k+1, k ), 1 )
616 colmax = cabs1( w( imax, k ) )
621 IF( max( absakk, colmax ).EQ.zero )
THEN
628 CALL
zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
638 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
657 CALL
zcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
658 CALL
zcopy( n-imax+1, a( imax, imax ), 1,
659 $ w( imax, k+1 ), 1 )
661 $ CALL
zgemv(
'No transpose', n-k+1, k-1, -cone,
662 $ a( k, 1 ), lda, w( imax, 1 ), ldw,
663 $ cone, w( k, k+1 ), 1 )
670 jmax = k - 1 +
izamax( imax-k, w( k, k+1 ), 1 )
671 rowmax = cabs1( w( jmax, k+1 ) )
677 itemp = imax +
izamax( n-imax, w( imax+1, k+1 ), 1)
678 dtemp = cabs1( w( itemp, k+1 ) )
679 IF( dtemp.GT.rowmax )
THEN
689 IF( .NOT.( cabs1( w( imax, k+1 ) ).LT.alpha*rowmax ) )
699 CALL
zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
706 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
725 CALL
zcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
731 IF( .NOT. done ) goto 72
739 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
743 CALL
zcopy( p-k, a( k, k ), 1, a( p, k ), lda )
744 CALL
zcopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
749 CALL
zswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
750 CALL
zswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
759 a( kp, k ) = a( kk, k )
760 CALL
zcopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
761 CALL
zcopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
765 CALL
zswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
766 CALL
zswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
769 IF( kstep.EQ.1 )
THEN
779 CALL
zcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
781 IF( cabs1( a( k, k ) ).GE.sfmin )
THEN
782 r1 = cone / a( k, k )
783 CALL
zscal( n-k, r1, a( k+1, k ), 1 )
784 ELSE IF( a( k, k ).NE.czero )
THEN
786 a( ii, k ) = a( ii, k ) / a( k, k )
805 d11 = w( k+1, k+1 ) / d21
806 d22 = w( k, k ) / d21
807 t = cone / ( d11*d22-cone )
809 a(
j, k ) = t*( ( d11*w(
j, k )-w(
j, k+1 ) ) /
811 a(
j, k+1 ) = t*( ( d22*w(
j, k+1 )-w(
j, k ) ) /
818 a( k, k ) = w( k, k )
819 a( k+1, k ) = w( k+1, k )
820 a( k+1, k+1 ) = w( k+1, k+1 )
826 IF( kstep.EQ.1 )
THEN
847 jb = min( nb, n-
j+1 )
851 DO 100 jj =
j,
j + jb - 1
852 CALL
zgemv(
'No transpose',
j+jb-jj, k-1, -cone,
853 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
860 $ CALL
zgemm(
'No transpose',
'Transpose', n-
j-jb+1, jb,
861 $ k-1, -cone, a(
j+jb, 1 ), lda, w(
j, 1 ), ldw,
862 $ cone, a(
j+jb,
j ), lda )
883 IF( jp2.NE.jj .AND.
j.GE.1 )
884 $ CALL
zswap(
j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
886 IF( jp1.NE.jj .AND. kstep.EQ.2 )
887 $ CALL
zswap(
j, a( jp1, 1 ), lda, a( jj, 1 ), lda )