194 INTEGER info, kb, lda, ldw, n, nb
198 DOUBLE PRECISION 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 )
211 INTEGER imax, itemp,
j, jb, jj, jmax, jp1, jp2, k, kk,
212 $ kw, kkw, kp, kstep, p, ii
214 DOUBLE PRECISION absakk, alpha, colmax, d11, d12, d21, d22,
215 $ dtemp, 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
dcopy( k, a( 1, k ), 1, w( 1, kw ), 1 )
268 $ CALL
dgemv(
'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 =
idamax( k-1, w( 1, kw ), 1 )
282 colmax = abs( w( imax, kw ) )
287 IF( max( absakk, colmax ).EQ.zero )
THEN
294 CALL
dcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
304 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
323 CALL
dcopy( imax, a( 1, imax ), 1, w( 1, kw-1 ), 1 )
324 CALL
dcopy( k-imax, a( imax, imax+1 ), lda,
325 $ w( imax+1, kw-1 ), 1 )
328 $ CALL
dgemv(
'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 +
idamax( k-imax, w( imax+1, kw-1 ),
339 rowmax = abs( w( jmax, kw-1 ) )
345 itemp =
idamax( imax-1, w( 1, kw-1 ), 1 )
346 dtemp = abs( w( itemp, kw-1 ) )
347 IF( dtemp.GT.rowmax )
THEN
357 IF( .NOT.(abs( w( imax, kw-1 ) ).LT.alpha*rowmax ) )
367 CALL
dcopy( k, w( 1, kw-1 ), 1, w( 1, kw ), 1 )
374 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
393 CALL
dcopy( 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
dcopy( k-p, a( p+1, k ), 1, a( p, p+1 ), lda )
416 CALL
dcopy( p, a( 1, k ), 1, a( 1, p ), 1 )
421 CALL
dswap( n-k+1, a( k, k ), lda, a( p, k ), lda )
422 CALL
dswap( n-kk+1, w( k, kkw ), ldw, w( p, kkw ), ldw )
431 a( kp, k ) = a( kk, k )
432 CALL
dcopy( k-1-kp, a( kp+1, kk ), 1, a( kp, kp+1 ),
434 CALL
dcopy( kp, a( 1, kk ), 1, a( 1, kp ), 1 )
439 CALL
dswap( n-kk+1, a( kk, kk ), lda, a( kp, kk ), lda )
440 CALL
dswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),
444 IF( kstep.EQ.1 )
THEN
454 CALL
dcopy( k, w( 1, kw ), 1, a( 1, k ), 1 )
456 IF( abs( a( k, k ) ).GE.sfmin )
THEN
458 CALL
dscal( 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
dgemv(
'No transpose', jj-
j+1, n-k, -one,
529 $ a(
j, k+1 ), lda, w( jj, kw+1 ), ldw, one,
536 $ CALL
dgemm(
'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
dswap( n-
j+1, a( jp2,
j ), lda, a( jj,
j ), lda )
562 IF( jp1.NE.jj .AND. kstep.EQ.2 )
563 $ CALL
dswap( 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
dcopy( n-k+1, a( k, k ), 1, w( k, k ), 1 )
594 $ CALL
dgemv(
'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 +
idamax( n-k, w( k+1, k ), 1 )
608 colmax = abs( w( imax, k ) )
613 IF( max( absakk, colmax ).EQ.zero )
THEN
620 CALL
dcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
630 IF( .NOT.( absakk.LT.alpha*colmax ) )
THEN
649 CALL
dcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1)
650 CALL
dcopy( n-imax+1, a( imax, imax ), 1,
651 $ w( imax, k+1 ), 1 )
653 $ CALL
dgemv(
'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 +
idamax( imax-k, w( k, k+1 ), 1 )
663 rowmax = abs( w( jmax, k+1 ) )
669 itemp = imax +
idamax( n-imax, w( imax+1, k+1 ), 1)
670 dtemp = abs( w( itemp, k+1 ) )
671 IF( dtemp.GT.rowmax )
THEN
681 IF( .NOT.( abs( w( imax, k+1 ) ).LT.alpha*rowmax ) )
691 CALL
dcopy( n-k+1, w( k, k+1 ), 1, w( k, k ), 1 )
698 ELSE IF( ( p.EQ.jmax ) .OR. ( rowmax.LE.colmax ) )
717 CALL
dcopy( 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
dcopy( p-k, a( k, k ), 1, a( p, k ), lda )
736 CALL
dcopy( n-p+1, a( p, k ), 1, a( p, p ), 1 )
741 CALL
dswap( k, a( k, 1 ), lda, a( p, 1 ), lda )
742 CALL
dswap( kk, w( k, 1 ), ldw, w( p, 1 ), ldw )
751 a( kp, k ) = a( kk, k )
752 CALL
dcopy( kp-k-1, a( k+1, kk ), 1, a( kp, k+1 ), lda )
753 CALL
dcopy( n-kp+1, a( kp, kk ), 1, a( kp, kp ), 1 )
757 CALL
dswap( kk, a( kk, 1 ), lda, a( kp, 1 ), lda )
758 CALL
dswap( kk, w( kk, 1 ), ldw, w( kp, 1 ), ldw )
761 IF( kstep.EQ.1 )
THEN
771 CALL
dcopy( n-k+1, w( k, k ), 1, a( k, k ), 1 )
773 IF( abs( a( k, k ) ).GE.sfmin )
THEN
775 CALL
dscal( 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
dgemv(
'No transpose',
j+jb-jj, k-1, -one,
845 $ a( jj, 1 ), lda, w( jj, 1 ), ldw, one,
852 $ CALL
dgemm(
'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
dswap(
j, a( jp2, 1 ), lda, a( jj, 1 ), lda )
878 IF( jp1.NE.jj .AND. kstep.EQ.2 )
879 $ CALL
dswap(
j, a( jp1, 1 ), lda, a( jj, 1 ), lda )
subroutine dlasyf_rook(UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO)
DLASYF_ROOK *> DLASYF_ROOK computes a partial factorization of a real symmetric matrix using the boun...
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
logical function lsame(CA, CB)
LSAME
subroutine dscal(N, DA, DX, INCX)
DSCAL
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
integer function idamax(N, DX, INCX)
IDAMAX
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV