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 )
LOGICAL function lsame(CA, CB)
LSAME
REAL function slamch(CMACH)
SLAMCH
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
subroutine clahef_rook(UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
Download CLAHEF_ROOK + dependencies <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clahef_rook.f"> [TGZ]</a> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clahef_rook.f"> [ZIP]</a> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clahef_rook.f"> [TXT]</a>
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine csscal(N, SA, CX, INCX)
CSSCAL
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
INTEGER function icamax(N, CX, INCX)
ICAMAX