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 )
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
LOGICAL function lsame(CA, CB)
LSAME
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
subroutine zpstf2(UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO)
ZPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric or complex Herm...
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
LOGICAL function disnan(DIN)
DISNAN tests input for NaN.
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
subroutine zpstrf(UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO)
ZPSTRF