142 SUBROUTINE zpstrf( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
151 INTEGER info, lda, n, rank
155 COMPLEX*16 a( lda, * )
156 DOUBLE PRECISION work( 2*n )
163 DOUBLE PRECISION one, zero
164 parameter( one = 1.0d+0, zero = 0.0d+0 )
166 parameter( cone = ( 1.0d+0, 0.0d+0 ) )
170 DOUBLE PRECISION ajj, dstop, dtemp
171 INTEGER i, itemp,
j, jb, k, nb, pvt
185 INTRINSIC dble, dconjg, max, min, sqrt, maxloc
192 upper =
lsame( uplo,
'U' )
193 IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
195 ELSE IF( n.LT.0 )
THEN
197 ELSE IF( lda.LT.max( 1, n ) )
THEN
201 CALL
xerbla(
'ZPSTRF', -info )
212 nb =
ilaenv( 1,
'ZPOTRF', uplo, n, -1, -1, -1 )
213 IF( nb.LE.1 .OR. nb.GE.n )
THEN
217 CALL
zpstf2( uplo, n, a( 1, 1 ), lda, piv, rank, tol, work,
232 work( i ) = dble( a( i, i ) )
234 pvt = maxloc( work( 1:n ), 1 )
235 ajj = dble( a( pvt, pvt ) )
236 IF( ajj.EQ.zero.OR.
disnan( ajj ) )
THEN
244 IF( tol.LT.zero )
THEN
245 dstop = n *
dlamch(
'Epsilon' ) * ajj
259 jb = min( nb, n-k+1 )
268 DO 150
j = k, k + jb - 1
277 work( i ) = work( i ) +
278 $ dble( dconjg( a(
j-1, i ) )*
281 work( n+i ) = dble( a( i, i ) ) - work( i )
286 itemp = maxloc( work( (n+
j):(2*n) ), 1 )
289 IF( ajj.LE.dstop.OR.
disnan( ajj ) )
THEN
299 a( pvt, pvt ) = a(
j,
j )
300 CALL
zswap(
j-1, a( 1,
j ), 1, a( 1, pvt ), 1 )
302 $ CALL
zswap( n-pvt, a(
j, pvt+1 ), lda,
303 $ a( pvt, pvt+1 ), lda )
304 DO 140 i =
j + 1, pvt - 1
305 ztemp = dconjg( a(
j, i ) )
306 a(
j, i ) = dconjg( a( i, pvt ) )
309 a(
j, pvt ) = dconjg( a(
j, pvt ) )
314 work(
j ) = work( pvt )
317 piv( pvt ) = piv(
j )
328 CALL
zgemv(
'Trans',
j-k, n-
j, -cone, a( k,
j+1 ),
329 $ lda, a( k,
j ), 1, cone, a(
j,
j+1 ),
332 CALL
zdscal( n-
j, one / ajj, a(
j,
j+1 ), lda )
340 CALL
zherk(
'Upper',
'Conj Trans', n-
j+1, jb, -one,
341 $ a( k,
j ), lda, one, a(
j,
j ), lda )
354 jb = min( nb, n-k+1 )
363 DO 200
j = k, k + jb - 1
372 work( i ) = work( i ) +
373 $ dble( dconjg( a( i,
j-1 ) )*
376 work( n+i ) = dble( a( i, i ) ) - work( i )
381 itemp = maxloc( work( (n+
j):(2*n) ), 1 )
384 IF( ajj.LE.dstop.OR.
disnan( ajj ) )
THEN
394 a( pvt, pvt ) = a(
j,
j )
395 CALL
zswap(
j-1, a(
j, 1 ), lda, a( pvt, 1 ), lda )
397 $ CALL
zswap( n-pvt, a( pvt+1,
j ), 1,
398 $ a( pvt+1, pvt ), 1 )
399 DO 190 i =
j + 1, pvt - 1
400 ztemp = dconjg( a( i,
j ) )
401 a( i,
j ) = dconjg( a( pvt, i ) )
404 a( pvt,
j ) = dconjg( a( pvt,
j ) )
410 work(
j ) = work( pvt )
413 piv( pvt ) = piv(
j )
424 CALL
zgemv(
'No Trans', n-
j,
j-k, -cone,
425 $ a(
j+1, k ), lda, a(
j, k ), lda, cone,
428 CALL
zdscal( n-
j, one / ajj, a(
j+1,
j ), 1 )
436 CALL
zherk(
'Lower',
'No Trans', n-
j+1, jb, -one,
437 $ a(
j, k ), lda, one, a(
j,
j ), lda )