194 INTEGER info, kb, lda, ldw, n, nb
198 COMPLEX a( lda, * ), w( ldw, * )
205 parameter ( zero = 0.0e+0, one = 1.0e+0 )
207 parameter ( cone = ( 1.0e+0, 0.0e+0 ) )
209 parameter ( eight = 8.0e+0, sevten = 17.0e+0 )
213 INTEGER imax, itemp, ii, j, jb, jj, jmax, jp1, jp2, k,
214 $ kk, kkw, kp, kstep, kw, p
215 REAL absakk, alpha, colmax, stemp, r1, rowmax, t,
217 COMPLEX d11, d21, d22, z
229 INTRINSIC abs, conjg, aimag, max, min,
REAL, sqrt
235 cabs1( z ) = abs(
REAL( Z ) ) + abs( aimag( 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 )
275 $
CALL ccopy( k-1, a( 1, k ), 1, w( 1, kw ), 1 )
276 w( k, kw ) =
REAL( A( K, K ) )
278 CALL cgemv(
'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 ) =
REAL( W( K, KW ) )
286 absakk = abs(
REAL( W( K, KW ) ) )
293 imax =
icamax( k-1, w( 1, kw ), 1 )
294 colmax = cabs1( w( imax, kw ) )
299 IF( max( absakk, colmax ).EQ.zero )
THEN
306 a( k, k ) =
REAL( W( K, KW ) )
308 $
CALL ccopy( k-1, w( 1, kw ), 1, a( 1, k ), 1 )
318 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
338 $
CALL ccopy( imax-1, a( 1, imax ), 1, w( 1, kw-1 ),
340 w( imax, kw-1 ) =
REAL( A( IMAX, IMAX ) )
342 CALL ccopy( k-imax, a( imax, imax+1 ), lda,
343 $ w( imax+1, kw-1 ), 1 )
344 CALL clacgv( k-imax, w( imax+1, kw-1 ), 1 )
347 CALL cgemv(
'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 ) =
REAL( W( IMAX, KW-1 ) )
358 jmax = imax +
icamax( k-imax, w( imax+1, kw-1 ),
360 rowmax = cabs1( w( jmax, kw-1 ) )
366 itemp =
icamax( imax-1, w( 1, kw-1 ), 1 )
367 stemp = cabs1( w( itemp, kw-1 ) )
368 IF( stemp.GT.rowmax )
THEN
379 IF( .NOT.( abs(
REAL( W( IMAX,KW-1 ) ) )
380 $ .LT.alpha*rowmax ) )
THEN
389 CALL ccopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
397 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
418 CALL ccopy( 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 ) =
REAL( A( K, K ) )
452 CALL ccopy( k-1-p, a( p+1, k ), 1, a( p, p+1 ),
454 CALL clacgv( k-1-p, a( p, p+1 ), lda )
456 $
CALL ccopy( p-1, a( 1, k ), 1, a( 1, p ), 1 )
464 $
CALL cswap( n-k, a( k, k+1 ), lda, a( p, k+1 ),
466 CALL cswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ),
480 a( kp, kp ) =
REAL( A( KK, KK ) )
481 CALL ccopy( kk-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
483 CALL clacgv( kk-1-kp, a( kp, kp+1 ), lda )
485 $
CALL ccopy( kp-1, a( 1, kk ), 1, a( 1, kp ), 1 )
493 $
CALL cswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),
495 CALL cswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
499 IF( kstep.EQ.1 )
THEN
517 CALL ccopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
526 t =
REAL( A( K, K ) )
527 IF( abs( t ).GE.sfmin )
THEN
529 CALL csscal( k-1, r1, a( 1, k ), 1 )
532 a( ii, k ) = a( ii, k ) / t
538 CALL clacgv( k-1, w( 1, kw ), 1 )
606 d11 = w( k, kw ) / conjg( d21 )
607 d22 = w( k-1, kw-1 ) / d21
608 t = one / (
REAL( 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 clacgv( k-1, w( 1, kw ), 1 )
631 CALL clacgv( 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 ) =
REAL( A( JJ, JJ ) )
667 CALL cgemv(
'No transpose', jj-j+1, n-k, -cone,
668 $ a( j, k+1 ), lda, w( jj, kw+1 ), ldw, cone,
670 a( jj, jj ) =
REAL( A( JJ, JJ ) )
676 $
CALL cgemm(
'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 cswap( 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 cswap( 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 ) =
REAL( A( K, K ) )
740 $
CALL ccopy( n-k, a( k+1, k ), 1, w( k+1, k ), 1 )
742 CALL cgemv(
'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 ) =
REAL( W( K, K ) )
750 absakk = abs(
REAL( W( K, K ) ) )
757 imax = k +
icamax( 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 ) =
REAL( W( K, K ) )
772 $
CALL ccopy( n-k, w( k+1, k ), 1, a( k+1, k ), 1 )
783 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
802 CALL ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
803 CALL clacgv( imax-k, w( k, k+1 ), 1 )
804 w( imax, k+1 ) =
REAL( A( IMAX, IMAX ) )
807 $
CALL ccopy( n-imax, a( imax+1, imax ), 1,
808 $ w( imax+1, k+1 ), 1 )
811 CALL cgemv(
'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 ) =
REAL( W( IMAX, K+1 ) )
822 jmax = k - 1 +
icamax( imax-k, w( k, k+1 ), 1 )
823 rowmax = cabs1( w( jmax, k+1 ) )
829 itemp = imax +
icamax( n-imax, w( imax+1, k+1 ), 1)
830 stemp = cabs1( w( itemp, k+1 ) )
831 IF( stemp.GT.rowmax )
THEN
842 IF( .NOT.( abs(
REAL( W( IMAX,K+1 ) ) )
843 $ .LT.alpha*rowmax ) )
THEN
852 CALL ccopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
860 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
881 CALL ccopy( 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 ) =
REAL( A( K, K ) )
911 CALL ccopy( p-k-1, a( k+1, k ), 1, a( p, k+1 ), lda )
912 CALL clacgv( p-k-1, a( p, k+1 ), lda )
914 $
CALL ccopy( n-p, a( p+1, k ), 1, a( p+1, p ), 1 )
922 $
CALL cswap( k-1, a( k, 1 ), lda, a( p, 1 ), lda )
923 CALL cswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
936 a( kp, kp ) =
REAL( A( KK, KK ) )
937 CALL ccopy( kp-kk-1, a( kk+1, kk ), 1, a( kp, kk+1 ),
939 CALL clacgv( kp-kk-1, a( kp, kk+1 ), lda )
941 $
CALL ccopy( n-kp, a( kp+1, kk ), 1, a( kp+1, kp ), 1 )
949 $
CALL cswap( k-1, a( kk, 1 ), lda, a( kp, 1 ), lda )
950 CALL cswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
953 IF( kstep.EQ.1 )
THEN
971 CALL ccopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
980 t =
REAL( A( K, K ) )
981 IF( abs( t ).GE.sfmin )
THEN
983 CALL csscal( n-k, r1, a( k+1, k ), 1 )
986 a( ii, k ) = a( ii, k ) / t
992 CALL clacgv( n-k, w( k+1, k ), 1 )
1060 d11 = w( k+1, k+1 ) / d21
1061 d22 = w( k, k ) / conjg( d21 )
1062 t = one / (
REAL( 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 clacgv( n-k, w( k+1, k ), 1 )
1085 CALL clacgv( 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 ) =
REAL( A( JJ, JJ ) )
1121 CALL cgemv(
'No transpose', j+jb-jj, k-1, -cone,
1122 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, cone,
1124 a( jj, jj ) =
REAL( A( JJ, JJ ) )
1130 $
CALL cgemm(
'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 cswap( j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
1162 IF( kstep.EQ.2 .AND. jp1.NE.jj .AND. j.GE.1 )
1163 $
CALL cswap( j, a( jp1, 1 ), lda, a( jj, 1 ), lda )
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
integer function icamax(N, CX, INCX)
ICAMAX
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
real function slamch(CMACH)
SLAMCH
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
logical function lsame(CA, CB)
LSAME
subroutine csscal(N, SA, CX, INCX)
CSSCAL