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 )