143 SUBROUTINE dpstrf( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
152 INTEGER INFO, LDA, N, RANK
156 DOUBLE PRECISION A( lda, * ), WORK( 2*n )
163 DOUBLE PRECISION ONE, ZERO
164 parameter ( one = 1.0d+0, zero = 0.0d+0 )
167 DOUBLE PRECISION AJJ, DSTOP, DTEMP
168 INTEGER I, ITEMP, J, JB, K, NB, PVT
172 DOUBLE PRECISION DLAMCH
174 LOGICAL LSAME, DISNAN
175 EXTERNAL dlamch, ilaenv, lsame, disnan
181 INTRINSIC max, min, sqrt, maxloc
188 upper = lsame( uplo,
'U' )
189 IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
191 ELSE IF( n.LT.0 )
THEN
193 ELSE IF( lda.LT.max( 1, n ) )
THEN
197 CALL xerbla(
'DPSTRF', -info )
208 nb = ilaenv( 1,
'DPOTRF', uplo, n, -1, -1, -1 )
209 IF( nb.LE.1 .OR. nb.GE.n )
THEN
213 CALL dpstf2( uplo, n, a( 1, 1 ), lda, piv, rank, tol, work,
230 IF( a( i, i ).GT.ajj )
THEN
235 IF( ajj.LE.zero.OR.disnan( ajj ) )
THEN
243 IF( tol.LT.zero )
THEN
244 dstop = n * dlamch(
'Epsilon' ) * ajj
258 jb = min( nb, n-k+1 )
267 DO 130 j = k, k + jb - 1
276 work( i ) = work( i ) + a( j-1, i )**2
278 work( n+i ) = a( i, i ) - work( i )
283 itemp = maxloc( work( (n+j):(2*n) ), 1 )
286 IF( ajj.LE.dstop.OR.disnan( ajj ) )
THEN
296 a( pvt, pvt ) = a( j, j )
297 CALL dswap( j-1, a( 1, j ), 1, a( 1, pvt ), 1 )
299 $
CALL dswap( n-pvt, a( j, pvt+1 ), lda,
300 $ a( pvt, pvt+1 ), lda )
301 CALL dswap( pvt-j-1, a( j, j+1 ), lda,
307 work( j ) = work( pvt )
310 piv( pvt ) = piv( j )
320 CALL dgemv(
'Trans', j-k, n-j, -one, a( k, j+1 ),
321 $ lda, a( k, j ), 1, one, a( j, j+1 ),
323 CALL dscal( n-j, one / ajj, a( j, j+1 ), lda )
331 CALL dsyrk(
'Upper',
'Trans', n-j+1, jb, -one,
332 $ a( k, j ), lda, one, a( j, j ), lda )
345 jb = min( nb, n-k+1 )
354 DO 170 j = k, k + jb - 1
363 work( i ) = work( i ) + a( i, j-1 )**2
365 work( n+i ) = a( i, i ) - work( i )
370 itemp = maxloc( work( (n+j):(2*n) ), 1 )
373 IF( ajj.LE.dstop.OR.disnan( ajj ) )
THEN
383 a( pvt, pvt ) = a( j, j )
384 CALL dswap( j-1, a( j, 1 ), lda, a( pvt, 1 ), lda )
386 $
CALL dswap( n-pvt, a( pvt+1, j ), 1,
387 $ a( pvt+1, pvt ), 1 )
388 CALL dswap( pvt-j-1, a( j+1, j ), 1, a( pvt, j+1 ),
394 work( j ) = work( pvt )
397 piv( pvt ) = piv( j )
407 CALL dgemv(
'No Trans', n-j, j-k, -one,
408 $ a( j+1, k ), lda, a( j, k ), lda, one,
410 CALL dscal( n-j, one / ajj, a( j+1, j ), 1 )
418 CALL dsyrk(
'Lower',
'No Trans', n-j+1, jb, -one,
419 $ a( j, k ), lda, one, a( j, j ), lda )
subroutine dpstf2(UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO)
DPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric positive semide...
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
subroutine dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine dpstrf(UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO)
DPSTRF computes the Cholesky factorization with complete pivoting of a real symmetric positive semide...
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dscal(N, DA, DX, INCX)
DSCAL