141 SUBROUTINE dpstrf( 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, 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(
'DPSTRF', -info )
206 nb =
ilaenv( 1,
'DPOTRF', uplo, n, -1, -1, -1 )
207 IF( nb.LE.1 .OR. nb.GE.n )
THEN
211 CALL
dpstf2( uplo, n, a( 1, 1 ), lda, piv, rank, tol, work,
228 IF( a( i, i ).GT.ajj )
THEN
233 IF( ajj.EQ.zero.OR.
disnan( ajj ) )
THEN
241 IF( tol.LT.zero )
THEN
242 dstop = n *
dlamch(
'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.dstop.OR.
disnan( ajj ) )
THEN
294 a( pvt, pvt ) = a(
j,
j )
295 CALL
dswap(
j-1, a( 1,
j ), 1, a( 1, pvt ), 1 )
297 $ CALL
dswap( n-pvt, a(
j, pvt+1 ), lda,
298 $ a( pvt, pvt+1 ), lda )
299 CALL
dswap( pvt-
j-1, a(
j,
j+1 ), lda,
305 work(
j ) = work( pvt )
308 piv( pvt ) = piv(
j )
318 CALL
dgemv(
'Trans',
j-k, n-
j, -one, a( k,
j+1 ),
319 $ lda, a( k,
j ), 1, one, a(
j,
j+1 ),
321 CALL
dscal( n-
j, one / ajj, a(
j,
j+1 ), lda )
329 CALL
dsyrk(
'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.dstop.OR.
disnan( ajj ) )
THEN
381 a( pvt, pvt ) = a(
j,
j )
382 CALL
dswap(
j-1, a(
j, 1 ), lda, a( pvt, 1 ), lda )
384 $ CALL
dswap( n-pvt, a( pvt+1,
j ), 1,
385 $ a( pvt+1, pvt ), 1 )
386 CALL
dswap( pvt-
j-1, a(
j+1,
j ), 1, a( pvt,
j+1 ),
392 work(
j ) = work( pvt )
395 piv( pvt ) = piv(
j )
405 CALL
dgemv(
'No Trans', n-
j,
j-k, -one,
406 $ a(
j+1, k ), lda, a(
j, k ), lda, one,
408 CALL
dscal( n-
j, one / ajj, a(
j+1,
j ), 1 )
416 CALL
dsyrk(
'Lower',
'No Trans', n-
j+1, jb, -one,
417 $ a(
j, k ), lda, one, a(
j,
j ), lda )
subroutine dpstrf(UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO)
DPSTRF
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
logical function disnan(DIN)
DISNAN tests input for NaN.
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...