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 )
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
LOGICAL function lsame(CA, CB)
LSAME
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
subroutine zlasyf_rook(UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
ZLASYF_ROOK computes a partial factorization of a complex symmetric matrix using the bounded Bunch-Ka...
INTEGER function izamax(N, ZX, INCX)
IZAMAX
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL