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 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