142 SUBROUTINE dpstf2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
151 INTEGER INFO, LDA, N, RANK
155 DOUBLE PRECISION A( lda, * ), WORK( 2*n )
162 DOUBLE PRECISION ONE, ZERO
163 parameter ( one = 1.0d+0, zero = 0.0d+0 )
166 DOUBLE PRECISION AJJ, DSTOP, DTEMP
167 INTEGER I, ITEMP, J, PVT
171 DOUBLE PRECISION DLAMCH
172 LOGICAL LSAME, DISNAN
173 EXTERNAL dlamch, lsame, disnan
179 INTRINSIC max, 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(
'DPSTF2', -info )
215 IF( a( i, i ).GT.ajj )
THEN
220 IF( ajj.LE.zero.OR.disnan( ajj ) )
THEN
228 IF( tol.LT.zero )
THEN
229 dstop = n * dlamch(
'Epsilon' ) * ajj
253 work( i ) = work( i ) + a( j-1, i )**2
255 work( n+i ) = a( i, i ) - work( i )
260 itemp = maxloc( work( (n+j):(2*n) ), 1 )
263 IF( ajj.LE.dstop.OR.disnan( ajj ) )
THEN
273 a( pvt, pvt ) = a( j, j )
274 CALL dswap( j-1, a( 1, j ), 1, a( 1, pvt ), 1 )
276 $
CALL dswap( n-pvt, a( j, pvt+1 ), lda,
277 $ a( pvt, pvt+1 ), lda )
278 CALL dswap( pvt-j-1, a( j, j+1 ), lda, a( j+1, pvt ), 1 )
283 work( j ) = work( pvt )
286 piv( pvt ) = piv( j )
296 CALL dgemv(
'Trans', j-1, n-j, -one, a( 1, j+1 ), lda,
297 $ a( 1, j ), 1, one, a( j, j+1 ), lda )
298 CALL dscal( n-j, one / ajj, a( j, j+1 ), lda )
316 work( i ) = work( i ) + a( i, j-1 )**2
318 work( n+i ) = a( i, i ) - work( i )
323 itemp = maxloc( work( (n+j):(2*n) ), 1 )
326 IF( ajj.LE.dstop.OR.disnan( ajj ) )
THEN
336 a( pvt, pvt ) = a( j, j )
337 CALL dswap( j-1, a( j, 1 ), lda, a( pvt, 1 ), lda )
339 $
CALL dswap( n-pvt, a( pvt+1, j ), 1, a( pvt+1, pvt ),
341 CALL dswap( pvt-j-1, a( j+1, j ), 1, a( pvt, j+1 ), lda )
346 work( j ) = work( pvt )
349 piv( pvt ) = piv( j )
359 CALL dgemv(
'No Trans', n-j, j-1, -one, a( j+1, 1 ), lda,
360 $ a( j, 1 ), lda, one, a( j+1, j ), 1 )
361 CALL dscal( n-j, one / ajj, a( j+1, j ), 1 )
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 dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dscal(N, DA, DX, INCX)
DSCAL