![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CPSTRF 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CPSTRF + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpstrf.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpstrf.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpstrf.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * REAL TOL 00025 * INTEGER INFO, LDA, N, RANK 00026 * CHARACTER UPLO 00027 * .. 00028 * .. Array Arguments .. 00029 * COMPLEX A( LDA, * ) 00030 * REAL WORK( 2*N ) 00031 * INTEGER PIV( N ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> CPSTRF computes the Cholesky factorization with complete 00041 *> pivoting of a complex Hermitian positive semidefinite matrix A. 00042 *> 00043 *> The factorization has the form 00044 *> P**T * A * P = U**H * U , if UPLO = 'U', 00045 *> P**T * A * P = L * L**H, if UPLO = 'L', 00046 *> where U is an upper triangular matrix and L is lower triangular, and 00047 *> P is stored as vector PIV. 00048 *> 00049 *> This algorithm does not attempt to check that A is positive 00050 *> semidefinite. This version of the algorithm calls level 3 BLAS. 00051 *> \endverbatim 00052 * 00053 * Arguments: 00054 * ========== 00055 * 00056 *> \param[in] UPLO 00057 *> \verbatim 00058 *> UPLO is CHARACTER*1 00059 *> Specifies whether the upper or lower triangular part of the 00060 *> symmetric matrix A is stored. 00061 *> = 'U': Upper triangular 00062 *> = 'L': Lower triangular 00063 *> \endverbatim 00064 *> 00065 *> \param[in] N 00066 *> \verbatim 00067 *> N is INTEGER 00068 *> The order of the matrix A. N >= 0. 00069 *> \endverbatim 00070 *> 00071 *> \param[in,out] A 00072 *> \verbatim 00073 *> A is COMPLEX array, dimension (LDA,N) 00074 *> On entry, the symmetric matrix A. If UPLO = 'U', the leading 00075 *> n by n upper triangular part of A contains the upper 00076 *> triangular part of the matrix A, and the strictly lower 00077 *> triangular part of A is not referenced. If UPLO = 'L', the 00078 *> leading n by n lower triangular part of A contains the lower 00079 *> triangular part of the matrix A, and the strictly upper 00080 *> triangular part of A is not referenced. 00081 *> 00082 *> On exit, if INFO = 0, the factor U or L from the Cholesky 00083 *> factorization as above. 00084 *> \endverbatim 00085 *> 00086 *> \param[in] LDA 00087 *> \verbatim 00088 *> LDA is INTEGER 00089 *> The leading dimension of the array A. LDA >= max(1,N). 00090 *> \endverbatim 00091 *> 00092 *> \param[out] PIV 00093 *> \verbatim 00094 *> PIV is INTEGER array, dimension (N) 00095 *> PIV is such that the nonzero entries are P( PIV(K), K ) = 1. 00096 *> \endverbatim 00097 *> 00098 *> \param[out] RANK 00099 *> \verbatim 00100 *> RANK is INTEGER 00101 *> The rank of A given by the number of steps the algorithm 00102 *> completed. 00103 *> \endverbatim 00104 *> 00105 *> \param[in] TOL 00106 *> \verbatim 00107 *> TOL is REAL 00108 *> User defined tolerance. If TOL < 0, then N*U*MAX( A(K,K) ) 00109 *> will be used. The algorithm terminates at the (K-1)st step 00110 *> if the pivot <= TOL. 00111 *> \endverbatim 00112 *> 00113 *> \param[out] WORK 00114 *> \verbatim 00115 *> WORK is REAL array, dimension (2*N) 00116 *> Work space. 00117 *> \endverbatim 00118 *> 00119 *> \param[out] INFO 00120 *> \verbatim 00121 *> INFO is INTEGER 00122 *> < 0: If INFO = -K, the K-th argument had an illegal value, 00123 *> = 0: algorithm completed successfully, and 00124 *> > 0: the matrix A is either rank deficient with computed rank 00125 *> as returned in RANK, or is indefinite. See Section 7 of 00126 *> LAPACK Working Note #161 for further information. 00127 *> \endverbatim 00128 * 00129 * Authors: 00130 * ======== 00131 * 00132 *> \author Univ. of Tennessee 00133 *> \author Univ. of California Berkeley 00134 *> \author Univ. of Colorado Denver 00135 *> \author NAG Ltd. 00136 * 00137 *> \date November 2011 00138 * 00139 *> \ingroup complexOTHERcomputational 00140 * 00141 * ===================================================================== 00142 SUBROUTINE CPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO ) 00143 * 00144 * -- LAPACK computational routine (version 3.4.0) -- 00145 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00146 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00147 * November 2011 00148 * 00149 * .. Scalar Arguments .. 00150 REAL TOL 00151 INTEGER INFO, LDA, N, RANK 00152 CHARACTER UPLO 00153 * .. 00154 * .. Array Arguments .. 00155 COMPLEX A( LDA, * ) 00156 REAL WORK( 2*N ) 00157 INTEGER PIV( N ) 00158 * .. 00159 * 00160 * ===================================================================== 00161 * 00162 * .. Parameters .. 00163 REAL ONE, ZERO 00164 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00165 COMPLEX CONE 00166 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00167 * .. 00168 * .. Local Scalars .. 00169 COMPLEX CTEMP 00170 REAL AJJ, SSTOP, STEMP 00171 INTEGER I, ITEMP, J, JB, K, NB, PVT 00172 LOGICAL UPPER 00173 * .. 00174 * .. External Functions .. 00175 REAL SLAMCH 00176 INTEGER ILAENV 00177 LOGICAL LSAME, SISNAN 00178 EXTERNAL SLAMCH, ILAENV, LSAME, SISNAN 00179 * .. 00180 * .. External Subroutines .. 00181 EXTERNAL CGEMV, CHERK, CLACGV, CPSTF2, CSSCAL, CSWAP, 00182 $ XERBLA 00183 * .. 00184 * .. Intrinsic Functions .. 00185 INTRINSIC CONJG, MAX, MIN, REAL, SQRT, MAXLOC 00186 * .. 00187 * .. Executable Statements .. 00188 * 00189 * Test the input parameters. 00190 * 00191 INFO = 0 00192 UPPER = LSAME( UPLO, 'U' ) 00193 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00194 INFO = -1 00195 ELSE IF( N.LT.0 ) THEN 00196 INFO = -2 00197 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00198 INFO = -4 00199 END IF 00200 IF( INFO.NE.0 ) THEN 00201 CALL XERBLA( 'CPSTRF', -INFO ) 00202 RETURN 00203 END IF 00204 * 00205 * Quick return if possible 00206 * 00207 IF( N.EQ.0 ) 00208 $ RETURN 00209 * 00210 * Get block size 00211 * 00212 NB = ILAENV( 1, 'CPOTRF', UPLO, N, -1, -1, -1 ) 00213 IF( NB.LE.1 .OR. NB.GE.N ) THEN 00214 * 00215 * Use unblocked code 00216 * 00217 CALL CPSTF2( UPLO, N, A( 1, 1 ), LDA, PIV, RANK, TOL, WORK, 00218 $ INFO ) 00219 GO TO 230 00220 * 00221 ELSE 00222 * 00223 * Initialize PIV 00224 * 00225 DO 100 I = 1, N 00226 PIV( I ) = I 00227 100 CONTINUE 00228 * 00229 * Compute stopping value 00230 * 00231 DO 110 I = 1, N 00232 WORK( I ) = REAL( A( I, I ) ) 00233 110 CONTINUE 00234 PVT = MAXLOC( WORK( 1:N ), 1 ) 00235 AJJ = REAL( A( PVT, PVT ) ) 00236 IF( AJJ.EQ.ZERO.OR.SISNAN( AJJ ) ) THEN 00237 RANK = 0 00238 INFO = 1 00239 GO TO 230 00240 END IF 00241 * 00242 * Compute stopping value if not supplied 00243 * 00244 IF( TOL.LT.ZERO ) THEN 00245 SSTOP = N * SLAMCH( 'Epsilon' ) * AJJ 00246 ELSE 00247 SSTOP = TOL 00248 END IF 00249 * 00250 * 00251 IF( UPPER ) THEN 00252 * 00253 * Compute the Cholesky factorization P**T * A * P = U**H * U 00254 * 00255 DO 160 K = 1, N, NB 00256 * 00257 * Account for last block not being NB wide 00258 * 00259 JB = MIN( NB, N-K+1 ) 00260 * 00261 * Set relevant part of first half of WORK to zero, 00262 * holds dot products 00263 * 00264 DO 120 I = K, N 00265 WORK( I ) = 0 00266 120 CONTINUE 00267 * 00268 DO 150 J = K, K + JB - 1 00269 * 00270 * Find pivot, test for exit, else swap rows and columns 00271 * Update dot products, compute possible pivots which are 00272 * stored in the second half of WORK 00273 * 00274 DO 130 I = J, N 00275 * 00276 IF( J.GT.K ) THEN 00277 WORK( I ) = WORK( I ) + 00278 $ REAL( CONJG( A( J-1, I ) )* 00279 $ A( J-1, I ) ) 00280 END IF 00281 WORK( N+I ) = REAL( A( I, I ) ) - WORK( I ) 00282 * 00283 130 CONTINUE 00284 * 00285 IF( J.GT.1 ) THEN 00286 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 ) 00287 PVT = ITEMP + J - 1 00288 AJJ = WORK( N+PVT ) 00289 IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN 00290 A( J, J ) = AJJ 00291 GO TO 220 00292 END IF 00293 END IF 00294 * 00295 IF( J.NE.PVT ) THEN 00296 * 00297 * Pivot OK, so can now swap pivot rows and columns 00298 * 00299 A( PVT, PVT ) = A( J, J ) 00300 CALL CSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 ) 00301 IF( PVT.LT.N ) 00302 $ CALL CSWAP( N-PVT, A( J, PVT+1 ), LDA, 00303 $ A( PVT, PVT+1 ), LDA ) 00304 DO 140 I = J + 1, PVT - 1 00305 CTEMP = CONJG( A( J, I ) ) 00306 A( J, I ) = CONJG( A( I, PVT ) ) 00307 A( I, PVT ) = CTEMP 00308 140 CONTINUE 00309 A( J, PVT ) = CONJG( A( J, PVT ) ) 00310 * 00311 * Swap dot products and PIV 00312 * 00313 STEMP = WORK( J ) 00314 WORK( J ) = WORK( PVT ) 00315 WORK( PVT ) = STEMP 00316 ITEMP = PIV( PVT ) 00317 PIV( PVT ) = PIV( J ) 00318 PIV( J ) = ITEMP 00319 END IF 00320 * 00321 AJJ = SQRT( AJJ ) 00322 A( J, J ) = AJJ 00323 * 00324 * Compute elements J+1:N of row J. 00325 * 00326 IF( J.LT.N ) THEN 00327 CALL CLACGV( J-1, A( 1, J ), 1 ) 00328 CALL CGEMV( 'Trans', J-K, N-J, -CONE, A( K, J+1 ), 00329 $ LDA, A( K, J ), 1, CONE, A( J, J+1 ), 00330 $ LDA ) 00331 CALL CLACGV( J-1, A( 1, J ), 1 ) 00332 CALL CSSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA ) 00333 END IF 00334 * 00335 150 CONTINUE 00336 * 00337 * Update trailing matrix, J already incremented 00338 * 00339 IF( K+JB.LE.N ) THEN 00340 CALL CHERK( 'Upper', 'Conj Trans', N-J+1, JB, -ONE, 00341 $ A( K, J ), LDA, ONE, A( J, J ), LDA ) 00342 END IF 00343 * 00344 160 CONTINUE 00345 * 00346 ELSE 00347 * 00348 * Compute the Cholesky factorization P**T * A * P = L * L**H 00349 * 00350 DO 210 K = 1, N, NB 00351 * 00352 * Account for last block not being NB wide 00353 * 00354 JB = MIN( NB, N-K+1 ) 00355 * 00356 * Set relevant part of first half of WORK to zero, 00357 * holds dot products 00358 * 00359 DO 170 I = K, N 00360 WORK( I ) = 0 00361 170 CONTINUE 00362 * 00363 DO 200 J = K, K + JB - 1 00364 * 00365 * Find pivot, test for exit, else swap rows and columns 00366 * Update dot products, compute possible pivots which are 00367 * stored in the second half of WORK 00368 * 00369 DO 180 I = J, N 00370 * 00371 IF( J.GT.K ) THEN 00372 WORK( I ) = WORK( I ) + 00373 $ REAL( CONJG( A( I, J-1 ) )* 00374 $ A( I, J-1 ) ) 00375 END IF 00376 WORK( N+I ) = REAL( A( I, I ) ) - WORK( I ) 00377 * 00378 180 CONTINUE 00379 * 00380 IF( J.GT.1 ) THEN 00381 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 ) 00382 PVT = ITEMP + J - 1 00383 AJJ = WORK( N+PVT ) 00384 IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN 00385 A( J, J ) = AJJ 00386 GO TO 220 00387 END IF 00388 END IF 00389 * 00390 IF( J.NE.PVT ) THEN 00391 * 00392 * Pivot OK, so can now swap pivot rows and columns 00393 * 00394 A( PVT, PVT ) = A( J, J ) 00395 CALL CSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA ) 00396 IF( PVT.LT.N ) 00397 $ CALL CSWAP( N-PVT, A( PVT+1, J ), 1, 00398 $ A( PVT+1, PVT ), 1 ) 00399 DO 190 I = J + 1, PVT - 1 00400 CTEMP = CONJG( A( I, J ) ) 00401 A( I, J ) = CONJG( A( PVT, I ) ) 00402 A( PVT, I ) = CTEMP 00403 190 CONTINUE 00404 A( PVT, J ) = CONJG( A( PVT, J ) ) 00405 * 00406 * Swap dot products and PIV 00407 * 00408 STEMP = WORK( J ) 00409 WORK( J ) = WORK( PVT ) 00410 WORK( PVT ) = STEMP 00411 ITEMP = PIV( PVT ) 00412 PIV( PVT ) = PIV( J ) 00413 PIV( J ) = ITEMP 00414 END IF 00415 * 00416 AJJ = SQRT( AJJ ) 00417 A( J, J ) = AJJ 00418 * 00419 * Compute elements J+1:N of column J. 00420 * 00421 IF( J.LT.N ) THEN 00422 CALL CLACGV( J-1, A( J, 1 ), LDA ) 00423 CALL CGEMV( 'No Trans', N-J, J-K, -CONE, 00424 $ A( J+1, K ), LDA, A( J, K ), LDA, CONE, 00425 $ A( J+1, J ), 1 ) 00426 CALL CLACGV( J-1, A( J, 1 ), LDA ) 00427 CALL CSSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 ) 00428 END IF 00429 * 00430 200 CONTINUE 00431 * 00432 * Update trailing matrix, J already incremented 00433 * 00434 IF( K+JB.LE.N ) THEN 00435 CALL CHERK( 'Lower', 'No Trans', N-J+1, JB, -ONE, 00436 $ A( J, K ), LDA, ONE, A( J, J ), LDA ) 00437 END IF 00438 * 00439 210 CONTINUE 00440 * 00441 END IF 00442 END IF 00443 * 00444 * Ran to completion, A has full rank 00445 * 00446 RANK = N 00447 * 00448 GO TO 230 00449 220 CONTINUE 00450 * 00451 * Rank is the number of steps completed. Set INFO = 1 to signal 00452 * that the factorization cannot be used to solve a system. 00453 * 00454 RANK = J - 1 00455 INFO = 1 00456 * 00457 230 CONTINUE 00458 RETURN 00459 * 00460 * End of CPSTRF 00461 * 00462 END