194 INTEGER info, kb, lda, ldw, n, nb
198 REAL a( lda, * ), w( ldw, * )
205 parameter( zero = 0.0e+0, one = 1.0e+0 )
207 parameter( eight = 8.0e+0, sevten = 17.0e+0 )
211 INTEGER imax, itemp,
j, jb, jj, jmax, jp1, jp2, k, kk,
212 $ kw, kkw, kp, kstep, p, ii
214 REAL absakk, alpha, colmax, d11, d12, d21, d22,
215 $ stemp, r1, rowmax, t, sfmin
227 INTRINSIC abs, max, min, sqrt
235 alpha = ( one+sqrt( sevten ) ) / eight
241 IF(
lsame( uplo,
'U' ) )
THEN
258 IF( ( k.LE.n-nb+1 .AND. nb.LT.n ) .OR. k.LT.1 )
266 CALL
scopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
268 $ CALL
sgemv(
'No transpose', k, n-k, -one, a( 1, k+1 ),
269 $ lda, w( k, kw+1 ), ldw, one, w( 1, kw ), 1 )
274 absakk = abs( w( k, kw ) )
281 imax =
isamax( k-1, w( 1, kw ), 1 )
282 colmax = abs( w( imax, kw ) )
287 IF( max( absakk, colmax ).EQ.zero )
THEN
294 CALL
scopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
304 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
323 CALL
scopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
324 CALL
scopy( k-imax, a( imax, imax+1 ), lda,
325 $ w( imax+1, kw-1 ), 1 )
328 $ CALL
sgemv(
'No transpose', k, n-k, -one,
329 $ a( 1, k+1 ), lda, w( imax, kw+1 ), ldw,
330 $ one, w( 1, kw-1 ), 1 )
337 jmax = imax +
isamax( k-imax, w( imax+1, kw-1 ),
339 rowmax = abs( w( jmax, kw-1 ) )
345 itemp =
isamax( imax-1, w( 1, kw-1 ), 1 )
346 stemp = abs( w( itemp, kw-1 ) )
347 IF( stemp.GT.rowmax )
THEN
357 IF( .NOT.(abs( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
367 CALL
scopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
374 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
393 CALL
scopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
399 IF( .NOT. done ) goto 12
411 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
415 CALL
scopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
416 CALL
scopy( p, a( 1, k ), 1, a( 1, p ), 1 )
421 CALL
sswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
422 CALL
sswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
431 a( kp, k ) = a( kk, k )
432 CALL
scopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
434 CALL
scopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
439 CALL
sswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
440 CALL
sswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
444 IF( kstep.EQ.1 )
THEN
454 CALL
scopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
456 IF( abs( a( k, k ) ).GE.sfmin )
THEN
458 CALL
sscal( k-1, r1, a( 1, k ), 1 )
459 ELSE IF( a( k, k ).NE.zero )
THEN
461 a( ii, k ) = a( ii, k ) / a( k, k )
481 d11 = w( k, kw ) / d12
482 d22 = w( k-1, kw-1 ) / d12
483 t = one / ( d11*d22-one )
485 a(
j, k-1 ) = t*( (d11*w(
j, kw-1 )-w(
j, kw ) ) /
487 a(
j, k ) = t*( ( d22*w(
j, kw )-w(
j, kw-1 ) ) /
494 a( k-1, k-1 ) = w( k-1, kw-1 )
495 a( k-1, k ) = w( k-1, kw )
496 a( k, k ) = w( k, kw )
502 IF( kstep.EQ.1 )
THEN
522 DO 50
j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
523 jb = min( nb, k-
j+1 )
527 DO 40 jj =
j,
j + jb - 1
528 CALL
sgemv(
'No transpose', jj-
j+1, n-k, -one,
529 $ a(
j, k+1 ), lda, w( jj, kw+1 ), ldw, one,
536 $ CALL
sgemm(
'No transpose',
'Transpose',
j-1, jb,
537 $ n-k, -one, a( 1, k+1 ), lda, w(
j, kw+1 ), ldw,
538 $ one, a( 1,
j ), lda )
559 IF( jp2.NE.jj .AND.
j.LE.n )
560 $ CALL
sswap( n-
j+1, a( jp2,
j ), lda, a( jj,
j ), lda )
562 IF( jp1.NE.jj .AND. kstep.EQ.2 )
563 $ CALL
sswap( n-
j+1, a( jp1,
j ), lda, a( jj,
j ), lda )
584 IF( ( k.GE.nb .AND. nb.LT.n ) .OR. k.GT.n )
592 CALL
scopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
594 $ CALL
sgemv(
'No transpose', n-k+1, k-1, -one, a( k, 1 ),
595 $ lda, w( k, 1 ), ldw, one, w( k, k ), 1 )
600 absakk = abs( w( k, k ) )
607 imax = k +
isamax( n-k, w( k+1, k ), 1 )
608 colmax = abs( w( imax, k ) )
613 IF( max( absakk, colmax ).EQ.zero )
THEN
620 CALL
scopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
630 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
649 CALL
scopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
650 CALL
scopy( n-imax+1, a( imax, imax ), 1,
651 $ w( imax, k+1 ), 1 )
653 $ CALL
sgemv(
'No transpose', n-k+1, k-1, -one,
654 $ a( k, 1 ), lda, w( imax, 1 ), ldw,
655 $ one, w( k, k+1 ), 1 )
662 jmax = k - 1 +
isamax( imax-k, w( k, k+1 ), 1 )
663 rowmax = abs( w( jmax, k+1 ) )
669 itemp = imax +
isamax( n-imax, w( imax+1, k+1 ), 1)
670 stemp = abs( w( itemp, k+1 ) )
671 IF( stemp.GT.rowmax )
THEN
681 IF( .NOT.( abs( w( imax, k+1 ) ).LT.alpha*rowmax ) )
691 CALL
scopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
698 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
717 CALL
scopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
723 IF( .NOT. done ) goto 72
731 IF( ( kstep.EQ.2 ) .AND. ( p.NE.k ) )
THEN
735 CALL
scopy( p-k, a( k, k ), 1, a( p, k ), lda )
736 CALL
scopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
741 CALL
sswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
742 CALL
sswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
751 a( kp, k ) = a( kk, k )
752 CALL
scopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
753 CALL
scopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
757 CALL
sswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
758 CALL
sswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
761 IF( kstep.EQ.1 )
THEN
771 CALL
scopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
773 IF( abs( a( k, k ) ).GE.sfmin )
THEN
775 CALL
sscal( n-k, r1, a( k+1, k ), 1 )
776 ELSE IF( a( k, k ).NE.zero )
THEN
778 a( ii, k ) = a( ii, k ) / a( k, k )
797 d11 = w( k+1, k+1 ) / d21
798 d22 = w( k, k ) / d21
799 t = one / ( d11*d22-one )
801 a(
j, k ) = t*( ( d11*w(
j, k )-w(
j, k+1 ) ) /
803 a(
j, k+1 ) = t*( ( d22*w(
j, k+1 )-w(
j, k ) ) /
810 a( k, k ) = w( k, k )
811 a( k+1, k ) = w( k+1, k )
812 a( k+1, k+1 ) = w( k+1, k+1 )
818 IF( kstep.EQ.1 )
THEN
839 jb = min( nb, n-
j+1 )
843 DO 100 jj =
j,
j + jb - 1
844 CALL
sgemv(
'No transpose',
j+jb-jj, k-1, -one,
845 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, one,
852 $ CALL
sgemm(
'No transpose',
'Transpose', n-
j-jb+1, jb,
853 $ k-1, -one, a(
j+jb, 1 ), lda, w(
j, 1 ), ldw,
854 $ one, a(
j+jb,
j ), lda )
875 IF( jp2.NE.jj .AND.
j.GE.1 )
876 $ CALL
sswap(
j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
878 IF( jp1.NE.jj .AND. kstep.EQ.2 )
879 $ CALL
sswap(
j, a( jp1, 1 ), lda, a( jj, 1 ), lda )
LOGICAL function lsame(CA, CB)
LSAME
INTEGER function isamax(N, SX, INCX)
ISAMAX
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
REAL function slamch(CMACH)
SLAMCH
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
subroutine slasyf_rook(UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
SLASYF_ROOK computes a partial factorization of a real symmetric matrix using the bounded Bunch-Kaufm...
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine sscal(N, SA, SX, INCX)
SSCAL