142 SUBROUTINE spstrf( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
151 INTEGER INFO, LDA, N, RANK
155 REAL A( lda, * ), WORK( 2*n )
163 parameter ( one = 1.0e+0, zero = 0.0e+0 )
166 REAL AJJ, SSTOP, STEMP
167 INTEGER I, ITEMP, J, JB, K, NB, PVT
173 LOGICAL LSAME, SISNAN
174 EXTERNAL slamch, ilaenv, lsame, sisnan
180 INTRINSIC max, min, sqrt, maxloc
187 upper = lsame( uplo,
'U' )
188 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
190 ELSE IF( n.LT.0 )
THEN
192 ELSE IF( lda.LT.max( 1, n ) )
THEN
196 CALL xerbla(
'SPSTRF', -info )
207 nb = ilaenv( 1,
'SPOTRF', uplo, n, -1, -1, -1 )
208 IF( nb.LE.1 .OR. nb.GE.n )
THEN
212 CALL spstf2( uplo, n, a( 1, 1 ), lda, piv, rank, tol, work,
229 IF( a( i, i ).GT.ajj )
THEN
234 IF( ajj.LE.zero.OR.sisnan( ajj ) )
THEN
242 IF( tol.LT.zero )
THEN
243 sstop = n * slamch(
'Epsilon' ) * ajj
257 jb = min( nb, n-k+1 )
266 DO 130 j = k, k + jb - 1
275 work( i ) = work( i ) + a( j-1, i )**2
277 work( n+i ) = a( i, i ) - work( i )
282 itemp = maxloc( work( (n+j):(2*n) ), 1 )
285 IF( ajj.LE.sstop.OR.sisnan( ajj ) )
THEN
295 a( pvt, pvt ) = a( j, j )
296 CALL sswap( j-1, a( 1, j ), 1, a( 1, pvt ), 1 )
298 $
CALL sswap( n-pvt, a( j, pvt+1 ), lda,
299 $ a( pvt, pvt+1 ), lda )
300 CALL sswap( pvt-j-1, a( j, j+1 ), lda,
306 work( j ) = work( pvt )
309 piv( pvt ) = piv( j )
319 CALL sgemv(
'Trans', j-k, n-j, -one, a( k, j+1 ),
320 $ lda, a( k, j ), 1, one, a( j, j+1 ),
322 CALL sscal( n-j, one / ajj, a( j, j+1 ), lda )
330 CALL ssyrk(
'Upper',
'Trans', n-j+1, jb, -one,
331 $ a( k, j ), lda, one, a( j, j ), lda )
344 jb = min( nb, n-k+1 )
353 DO 170 j = k, k + jb - 1
362 work( i ) = work( i ) + a( i, j-1 )**2
364 work( n+i ) = a( i, i ) - work( i )
369 itemp = maxloc( work( (n+j):(2*n) ), 1 )
372 IF( ajj.LE.sstop.OR.sisnan( ajj ) )
THEN
382 a( pvt, pvt ) = a( j, j )
383 CALL sswap( j-1, a( j, 1 ), lda, a( pvt, 1 ), lda )
385 $
CALL sswap( n-pvt, a( pvt+1, j ), 1,
386 $ a( pvt+1, pvt ), 1 )
387 CALL sswap( pvt-j-1, a( j+1, j ), 1, a( pvt, j+1 ),
393 work( j ) = work( pvt )
396 piv( pvt ) = piv( j )
406 CALL sgemv(
'No Trans', n-j, j-k, -one,
407 $ a( j+1, k ), lda, a( j, k ), lda, one,
409 CALL sscal( n-j, one / ajj, a( j+1, j ), 1 )
417 CALL ssyrk(
'Lower',
'No Trans', n-j+1, jb, -one,
418 $ a( j, k ), lda, one, a( j, j ), lda )
subroutine spstrf(UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO)
SPSTRF computes the Cholesky factorization with complete pivoting of a real symmetric positive semide...
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
subroutine spstf2(UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO)
SPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric positive semide...
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP