1 *> \brief \b CPSTRF computes the Cholesky factorization with complete pivoting of complex Hermitian positive semidefinite matrix.
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download CPSTRF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpstrf.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpstrf.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpstrf.f">
21 * SUBROUTINE CPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
23 * .. Scalar Arguments ..
25 * INTEGER INFO, LDA, N, RANK
28 * .. Array Arguments ..
40 *> CPSTRF computes the Cholesky factorization with complete
41 *> pivoting of a complex Hermitian positive semidefinite matrix A.
43 *> The factorization has the form
44 *> P**T * A * P = U**H * U , if UPLO = 'U',
45 *> P**T * A * P = L * L**H, if UPLO = 'L',
46 *> where U is an upper triangular matrix and L is lower triangular, and
47 *> P is stored as vector PIV.
49 *> This algorithm does not attempt to check that A is positive
50 *> semidefinite. This version of the algorithm calls level 3 BLAS.
58 *> UPLO is CHARACTER*1
59 *> Specifies whether the upper or lower triangular part of the
60 *> symmetric matrix A is stored.
61 *> = 'U': Upper triangular
62 *> = 'L': Lower triangular
68 *> The order of the matrix A. N >= 0.
73 *> A is COMPLEX array, dimension (LDA,N)
74 *> On entry, the symmetric matrix A. If UPLO = 'U', the leading
75 *> n by n upper triangular part of A contains the upper
76 *> triangular part of the matrix A, and the strictly lower
77 *> triangular part of A is not referenced. If UPLO = 'L', the
78 *> leading n by n lower triangular part of A contains the lower
79 *> triangular part of the matrix A, and the strictly upper
80 *> triangular part of A is not referenced.
82 *> On exit, if INFO = 0, the factor U or L from the Cholesky
83 *> factorization as above.
89 *> The leading dimension of the array A. LDA >= max(1,N).
94 *> PIV is INTEGER array, dimension (N)
95 *> PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
101 *> The rank of A given by the number of steps the algorithm
108 *> User defined tolerance. If TOL < 0, then N*U*MAX( A(K,K) )
109 *> will be used. The algorithm terminates at the (K-1)st step
110 *> if the pivot <= TOL.
115 *> WORK is REAL array, dimension (2*N)
122 *> < 0: If INFO = -K, the K-th argument had an illegal value,
123 *> = 0: algorithm completed successfully, and
124 *> > 0: the matrix A is either rank deficient with computed rank
125 *> as returned in RANK, or is not positive semidefinite. See
126 *> Section 7 of LAPACK Working Note #161 for further
133 *> \author Univ. of Tennessee
134 *> \author Univ. of California Berkeley
135 *> \author Univ. of Colorado Denver
138 *> \date November 2015
140 *> \ingroup complexOTHERcomputational
142 * =====================================================================
143 SUBROUTINE CPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
145 * -- LAPACK computational routine (version 3.6.0) --
146 * -- LAPACK is a software package provided by Univ. of Tennessee, --
147 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
150 * .. Scalar Arguments ..
152 INTEGER INFO, LDA, N, RANK
155 * .. Array Arguments ..
161 * =====================================================================
165 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
167 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) )
169 * .. Local Scalars ..
171 REAL AJJ, SSTOP, STEMP
172 INTEGER I, ITEMP, J, JB, K, NB, PVT
175 * .. External Functions ..
178 LOGICAL LSAME, SISNAN
179 EXTERNAL SLAMCH, ILAENV, LSAME, SISNAN
181 * .. External Subroutines ..
182 EXTERNAL CGEMV, CHERK, CLACGV, CPSTF2, CSSCAL, CSWAP,
185 * .. Intrinsic Functions ..
186 INTRINSIC CONJG, MAX, MIN, REAL, SQRT, MAXLOC
188 * .. Executable Statements ..
190 * Test the input parameters.
193 UPPER = LSAME( UPLO, 'U' )
194 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
196 ELSE IF( N.LT.0 ) THEN
198 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
202 CALL XERBLA( 'CPSTRF', -INFO )
206 * Quick return if possible
213 NB = ILAENV( 1, 'CPOTRF', UPLO, N, -1, -1, -1 )
214 IF( NB.LE.1 .OR. NB.GE.N ) THEN
218 CALL CPSTF2( UPLO, N, A( 1, 1 ), LDA, PIV, RANK, TOL, WORK,
230 * Compute stopping value
233 WORK( I ) = REAL( A( I, I ) )
235 PVT = MAXLOC( WORK( 1:N ), 1 )
236 AJJ = REAL( A( PVT, PVT ) )
237 IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN
243 * Compute stopping value if not supplied
245 IF( TOL.LT.ZERO ) THEN
246 SSTOP = N * SLAMCH( 'Epsilon' ) * AJJ
254 * Compute the Cholesky factorization P**T * A * P = U**H * U
258 * Account for last block not being NB wide
260 JB = MIN( NB, N-K+1 )
262 * Set relevant part of first half of WORK to zero,
269 DO 150 J = K, K + JB - 1
271 * Find pivot, test for exit, else swap rows and columns
272 * Update dot products, compute possible pivots which are
273 * stored in the second half of WORK
278 WORK( I ) = WORK( I ) +
279 $ REAL( CONJG( A( J-1, I ) )*
282 WORK( N+I ) = REAL( A( I, I ) ) - WORK( I )
287 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
290 IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN
298 * Pivot OK, so can now swap pivot rows and columns
300 A( PVT, PVT ) = A( J, J )
301 CALL CSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
303 $ CALL CSWAP( N-PVT, A( J, PVT+1 ), LDA,
304 $ A( PVT, PVT+1 ), LDA )
305 DO 140 I = J + 1, PVT - 1
306 CTEMP = CONJG( A( J, I ) )
307 A( J, I ) = CONJG( A( I, PVT ) )
310 A( J, PVT ) = CONJG( A( J, PVT ) )
312 * Swap dot products and PIV
315 WORK( J ) = WORK( PVT )
318 PIV( PVT ) = PIV( J )
325 * Compute elements J+1:N of row J.
328 CALL CLACGV( J-1, A( 1, J ), 1 )
329 CALL CGEMV( 'Trans', J-K, N-J, -CONE, A( K, J+1 ),
330 $ LDA, A( K, J ), 1, CONE, A( J, J+1 ),
332 CALL CLACGV( J-1, A( 1, J ), 1 )
333 CALL CSSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
338 * Update trailing matrix, J already incremented
341 CALL CHERK( 'Upper', 'Conj Trans', N-J+1, JB, -ONE,
342 $ A( K, J ), LDA, ONE, A( J, J ), LDA )
349 * Compute the Cholesky factorization P**T * A * P = L * L**H
353 * Account for last block not being NB wide
355 JB = MIN( NB, N-K+1 )
357 * Set relevant part of first half of WORK to zero,
364 DO 200 J = K, K + JB - 1
366 * Find pivot, test for exit, else swap rows and columns
367 * Update dot products, compute possible pivots which are
368 * stored in the second half of WORK
373 WORK( I ) = WORK( I ) +
374 $ REAL( CONJG( A( I, J-1 ) )*
377 WORK( N+I ) = REAL( A( I, I ) ) - WORK( I )
382 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
385 IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN
393 * Pivot OK, so can now swap pivot rows and columns
395 A( PVT, PVT ) = A( J, J )
396 CALL CSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
398 $ CALL CSWAP( N-PVT, A( PVT+1, J ), 1,
399 $ A( PVT+1, PVT ), 1 )
400 DO 190 I = J + 1, PVT - 1
401 CTEMP = CONJG( A( I, J ) )
402 A( I, J ) = CONJG( A( PVT, I ) )
405 A( PVT, J ) = CONJG( A( PVT, J ) )
407 * Swap dot products and PIV
410 WORK( J ) = WORK( PVT )
413 PIV( PVT ) = PIV( J )
420 * Compute elements J+1:N of column J.
423 CALL CLACGV( J-1, A( J, 1 ), LDA )
424 CALL CGEMV( 'No Trans', N-J, J-K, -CONE,
425 $ A( J+1, K ), LDA, A( J, K ), LDA, CONE,
427 CALL CLACGV( J-1, A( J, 1 ), LDA )
428 CALL CSSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
433 * Update trailing matrix, J already incremented
436 CALL CHERK( 'Lower', 'No Trans', N-J+1, JB, -ONE,
437 $ A( J, K ), LDA, ONE, A( J, J ), LDA )
445 * Ran to completion, A has full rank
452 * Rank is the number of steps completed. Set INFO = 1 to signal
453 * that the factorization cannot be used to solve a system.