142 SUBROUTINE cpstrf( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
151 INTEGER INFO, LDA, N, RANK
164 parameter( one = 1.0e+0, zero = 0.0e+0 )
166 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
170 REAL AJJ, SSTOP, STEMP
171 INTEGER I, ITEMP, J, JB, K, NB, PVT
177 LOGICAL LSAME, SISNAN
178 EXTERNAL slamch, ilaenv, lsame, sisnan
185 INTRINSIC conjg, max, min,
REAL, 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(
'CPSTRF', -info )
212 nb = ilaenv( 1,
'CPOTRF', uplo, n, -1, -1, -1 )
213 IF( nb.LE.1 .OR. nb.GE.n )
THEN
217 CALL cpstf2( uplo, n, a( 1, 1 ), lda, piv, rank, tol, work,
232 work( i ) =
REAL( A( I, I ) )
234 pvt = maxloc( work( 1:n ), 1 )
235 ajj =
REAL( A( PVT, PVT ) )
236 IF( ajj.EQ.zero.OR.sisnan( ajj ) )
THEN
244 IF( tol.LT.zero )
THEN
245 sstop = n * slamch(
'Epsilon' ) * ajj
259 jb = min( nb, n-k+1 )
268 DO 150 j = k, k + jb - 1
277 work( i ) = work( i ) +
278 $
REAL( CONJG( A( J-1, I ) )*
281 work( n+i ) =
REAL( A( I, I ) ) - WORK( i )
286 itemp = maxloc( work( (n+j):(2*n) ), 1 )
289 IF( ajj.LE.sstop.OR.sisnan( ajj ) )
THEN
299 a( pvt, pvt ) = a( j, j )
300 CALL cswap( j-1, a( 1, j ), 1, a( 1, pvt ), 1 )
302 $
CALL cswap( n-pvt, a( j, pvt+1 ), lda,
303 $ a( pvt, pvt+1 ), lda )
304 DO 140 i = j + 1, pvt - 1
305 ctemp = conjg( a( j, i ) )
306 a( j, i ) = conjg( a( i, pvt ) )
309 a( j, pvt ) = conjg( a( j, pvt ) )
314 work( j ) = work( pvt )
317 piv( pvt ) = piv( j )
327 CALL clacgv( j-1, a( 1, j ), 1 )
328 CALL cgemv(
'Trans', j-k, n-j, -cone, a( k, j+1 ),
329 $ lda, a( k, j ), 1, cone, a( j, j+1 ),
331 CALL clacgv( j-1, a( 1, j ), 1 )
332 CALL csscal( n-j, one / ajj, a( j, j+1 ), lda )
340 CALL cherk(
'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 $
REAL( CONJG( A( I, J-1 ) )*
376 work( n+i ) =
REAL( A( I, I ) ) - WORK( i )
381 itemp = maxloc( work( (n+j):(2*n) ), 1 )
384 IF( ajj.LE.sstop.OR.sisnan( ajj ) )
THEN
394 a( pvt, pvt ) = a( j, j )
395 CALL cswap( j-1, a( j, 1 ), lda, a( pvt, 1 ), lda )
397 $
CALL cswap( n-pvt, a( pvt+1, j ), 1,
398 $ a( pvt+1, pvt ), 1 )
399 DO 190 i = j + 1, pvt - 1
400 ctemp = conjg( a( i, j ) )
401 a( i, j ) = conjg( a( pvt, i ) )
404 a( pvt, j ) = conjg( a( pvt, j ) )
409 work( j ) = work( pvt )
412 piv( pvt ) = piv( j )
422 CALL clacgv( j-1, a( j, 1 ), lda )
423 CALL cgemv(
'No Trans', n-j, j-k, -cone,
424 $ a( j+1, k ), lda, a( j, k ), lda, cone,
426 CALL clacgv( j-1, a( j, 1 ), lda )
427 CALL csscal( n-j, one / ajj, a( j+1, j ), 1 )
435 CALL cherk(
'Lower',
'No Trans', n-j+1, jb, -one,
436 $ a( j, k ), lda, one, a( j, j ), lda )
subroutine cpstf2(UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO)
CPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric or complex Herm...
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
subroutine cpstrf(UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO)
CPSTRF
subroutine csscal(N, SA, CX, INCX)
CSSCAL
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK