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
175 DOUBLE PRECISION DLAMCH
177 LOGICAL LSAME, DISNAN
178 EXTERNAL dlamch, ilaenv, lsame, disnan
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 )
327 CALL zlacgv( j-1, a( 1, j ), 1 )
328 CALL zgemv(
'Trans', j-k, n-j, -cone, a( k, j+1 ),
329 $ lda, a( k, j ), 1, cone, a( j, j+1 ),
331 CALL zlacgv( j-1, a( 1, 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 )
423 CALL zlacgv( j-1, a( j, 1 ), lda )
424 CALL zgemv(
'No Trans', n-j, j-k, -cone,
425 $ a( j+1, k ), lda, a( j, k ), lda, cone,
427 CALL zlacgv( j-1, a( j, 1 ), lda )
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
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.
subroutine zpstrf(UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO)
ZPSTRF