141 SUBROUTINE spstrf( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
150 INTEGER info, lda, n, rank
154 REAL a( lda, * ), work( 2*n )
162 parameter( one = 1.0e+0, zero = 0.0e+0 )
165 REAL ajj, sstop, stemp
166 INTEGER i, itemp,
j, jb, k, nb, pvt
179 INTRINSIC max, min, sqrt, maxloc
186 upper =
lsame( uplo,
'U' )
187 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
189 ELSE IF( n.LT.0 )
THEN
191 ELSE IF( lda.LT.max( 1, n ) )
THEN
195 CALL
xerbla(
'SPSTRF', -info )
206 nb =
ilaenv( 1,
'SPOTRF', uplo, n, -1, -1, -1 )
207 IF( nb.LE.1 .OR. nb.GE.n )
THEN
211 CALL
spstf2( uplo, n, a( 1, 1 ), lda, piv, rank, tol, work,
228 IF( a( i, i ).GT.ajj )
THEN
233 IF( ajj.EQ.zero.OR.
sisnan( ajj ) )
THEN
241 IF( tol.LT.zero )
THEN
242 sstop = n *
slamch(
'Epsilon' ) * ajj
256 jb = min( nb, n-k+1 )
265 DO 130
j = k, k + jb - 1
274 work( i ) = work( i ) + a(
j-1, i )**2
276 work( n+i ) = a( i, i ) - work( i )
281 itemp = maxloc( work( (n+
j):(2*n) ), 1 )
284 IF( ajj.LE.sstop.OR.
sisnan( ajj ) )
THEN
294 a( pvt, pvt ) = a(
j,
j )
295 CALL
sswap(
j-1, a( 1,
j ), 1, a( 1, pvt ), 1 )
297 $ CALL
sswap( n-pvt, a(
j, pvt+1 ), lda,
298 $ a( pvt, pvt+1 ), lda )
299 CALL
sswap( pvt-
j-1, a(
j,
j+1 ), lda,
305 work(
j ) = work( pvt )
308 piv( pvt ) = piv(
j )
318 CALL
sgemv(
'Trans',
j-k, n-
j, -one, a( k,
j+1 ),
319 $ lda, a( k,
j ), 1, one, a(
j,
j+1 ),
321 CALL
sscal( n-
j, one / ajj, a(
j,
j+1 ), lda )
329 CALL
ssyrk(
'Upper',
'Trans', n-
j+1, jb, -one,
330 $ a( k,
j ), lda, one, a(
j,
j ), lda )
343 jb = min( nb, n-k+1 )
352 DO 170
j = k, k + jb - 1
361 work( i ) = work( i ) + a( i,
j-1 )**2
363 work( n+i ) = a( i, i ) - work( i )
368 itemp = maxloc( work( (n+
j):(2*n) ), 1 )
371 IF( ajj.LE.sstop.OR.
sisnan( ajj ) )
THEN
381 a( pvt, pvt ) = a(
j,
j )
382 CALL
sswap(
j-1, a(
j, 1 ), lda, a( pvt, 1 ), lda )
384 $ CALL
sswap( n-pvt, a( pvt+1,
j ), 1,
385 $ a( pvt+1, pvt ), 1 )
386 CALL
sswap( pvt-
j-1, a(
j+1,
j ), 1, a( pvt,
j+1 ),
392 work(
j ) = work( pvt )
395 piv( pvt ) = piv(
j )
405 CALL
sgemv(
'No Trans', n-
j,
j-k, -one,
406 $ a(
j+1, k ), lda, a(
j, k ), lda, one,
408 CALL
sscal( n-
j, one / ajj, a(
j+1,
j ), 1 )
416 CALL
ssyrk(
'Lower',
'No Trans', n-
j+1, jb, -one,
417 $ a(
j, k ), lda, one, a(
j,
j ), lda )
subroutine spstf2(UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO)
SPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric or complex Herm...
LOGICAL function lsame(CA, CB)
LSAME
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
REAL function slamch(CMACH)
SLAMCH
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
subroutine spstrf(UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO)
SPSTRF
subroutine xerbla(SRNAME, INFO)
XERBLA
LOGICAL function sisnan(SIN)
SISNAN tests input for NaN.
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine sscal(N, SA, SX, INCX)
SSCAL