141 SUBROUTINE dpstf2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
150 INTEGER info, lda, n, rank
154 DOUBLE PRECISION a( lda, * ), work( 2*n )
161 DOUBLE PRECISION one, zero
162 parameter( one = 1.0d+0, zero = 0.0d+0 )
165 DOUBLE PRECISION ajj, dstop, dtemp
166 INTEGER i, itemp,
j, pvt
178 INTRINSIC max, sqrt, maxloc
185 upper =
lsame( uplo,
'U' )
186 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
188 ELSE IF( n.LT.0 )
THEN
190 ELSE IF( lda.LT.max( 1, n ) )
THEN
194 CALL
xerbla(
'DPSTF2', -info )
214 IF( a( i, i ).GT.ajj )
THEN
219 IF( ajj.EQ.zero.OR.
disnan( ajj ) )
THEN
227 IF( tol.LT.zero )
THEN
228 dstop = n *
dlamch(
'Epsilon' ) * ajj
252 work( i ) = work( i ) + a(
j-1, i )**2
254 work( n+i ) = a( i, i ) - work( i )
259 itemp = maxloc( work( (n+
j):(2*n) ), 1 )
262 IF( ajj.LE.dstop.OR.
disnan( ajj ) )
THEN
272 a( pvt, pvt ) = a(
j,
j )
273 CALL
dswap(
j-1, a( 1,
j ), 1, a( 1, pvt ), 1 )
275 $ CALL
dswap( n-pvt, a(
j, pvt+1 ), lda,
276 $ a( pvt, pvt+1 ), lda )
277 CALL
dswap( pvt-
j-1, a(
j,
j+1 ), lda, a(
j+1, pvt ), 1 )
282 work(
j ) = work( pvt )
285 piv( pvt ) = piv(
j )
295 CALL
dgemv(
'Trans',
j-1, n-
j, -one, a( 1,
j+1 ), lda,
296 $ a( 1,
j ), 1, one, a(
j,
j+1 ), lda )
297 CALL
dscal( n-
j, one / ajj, a(
j,
j+1 ), lda )
315 work( i ) = work( i ) + a( i,
j-1 )**2
317 work( n+i ) = a( i, i ) - work( i )
322 itemp = maxloc( work( (n+
j):(2*n) ), 1 )
325 IF( ajj.LE.dstop.OR.
disnan( ajj ) )
THEN
335 a( pvt, pvt ) = a(
j,
j )
336 CALL
dswap(
j-1, a(
j, 1 ), lda, a( pvt, 1 ), lda )
338 $ CALL
dswap( n-pvt, a( pvt+1,
j ), 1, a( pvt+1, pvt ),
340 CALL
dswap( pvt-
j-1, a(
j+1,
j ), 1, a( pvt,
j+1 ), lda )
345 work(
j ) = work( pvt )
348 piv( pvt ) = piv(
j )
358 CALL
dgemv(
'No Trans', n-
j,
j-1, -one, a(
j+1, 1 ), lda,
359 $ a(
j, 1 ), lda, one, a(
j+1,
j ), 1 )
360 CALL
dscal( n-
j, one / ajj, a(
j+1,
j ), 1 )
LOGICAL function lsame(CA, CB)
LSAME
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dscal(N, DA, DX, INCX)
DSCAL
LOGICAL function disnan(DIN)
DISNAN tests input for NaN.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
subroutine dpstf2(UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO)
DPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric or complex Herm...