![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CPSTF2 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CPSTF2 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpstf2.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpstf2.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpstf2.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CPSTF2( 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 *> CPSTF2 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 2 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[out] PIV 00087 *> \verbatim 00088 *> PIV is INTEGER array, dimension (N) 00089 *> PIV is such that the nonzero entries are P( PIV(K), K ) = 1. 00090 *> \endverbatim 00091 *> 00092 *> \param[out] RANK 00093 *> \verbatim 00094 *> RANK is INTEGER 00095 *> The rank of A given by the number of steps the algorithm 00096 *> completed. 00097 *> \endverbatim 00098 *> 00099 *> \param[in] TOL 00100 *> \verbatim 00101 *> TOL is REAL 00102 *> User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) ) 00103 *> will be used. The algorithm terminates at the (K-1)st step 00104 *> if the pivot <= TOL. 00105 *> \endverbatim 00106 *> 00107 *> \param[in] LDA 00108 *> \verbatim 00109 *> LDA is INTEGER 00110 *> The leading dimension of the array A. LDA >= max(1,N). 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 CPSTF2( 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, PVT 00172 LOGICAL UPPER 00173 * .. 00174 * .. External Functions .. 00175 REAL SLAMCH 00176 LOGICAL LSAME, SISNAN 00177 EXTERNAL SLAMCH, LSAME, SISNAN 00178 * .. 00179 * .. External Subroutines .. 00180 EXTERNAL CGEMV, CLACGV, CSSCAL, CSWAP, XERBLA 00181 * .. 00182 * .. Intrinsic Functions .. 00183 INTRINSIC CONJG, MAX, REAL, SQRT 00184 * .. 00185 * .. Executable Statements .. 00186 * 00187 * Test the input parameters 00188 * 00189 INFO = 0 00190 UPPER = LSAME( UPLO, 'U' ) 00191 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00192 INFO = -1 00193 ELSE IF( N.LT.0 ) THEN 00194 INFO = -2 00195 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00196 INFO = -4 00197 END IF 00198 IF( INFO.NE.0 ) THEN 00199 CALL XERBLA( 'CPSTF2', -INFO ) 00200 RETURN 00201 END IF 00202 * 00203 * Quick return if possible 00204 * 00205 IF( N.EQ.0 ) 00206 $ RETURN 00207 * 00208 * Initialize PIV 00209 * 00210 DO 100 I = 1, N 00211 PIV( I ) = I 00212 100 CONTINUE 00213 * 00214 * Compute stopping value 00215 * 00216 DO 110 I = 1, N 00217 WORK( I ) = REAL( A( I, I ) ) 00218 110 CONTINUE 00219 PVT = MAXLOC( WORK( 1:N ), 1 ) 00220 AJJ = REAL ( A( PVT, PVT ) ) 00221 IF( AJJ.EQ.ZERO.OR.SISNAN( AJJ ) ) THEN 00222 RANK = 0 00223 INFO = 1 00224 GO TO 200 00225 END IF 00226 * 00227 * Compute stopping value if not supplied 00228 * 00229 IF( TOL.LT.ZERO ) THEN 00230 SSTOP = N * SLAMCH( 'Epsilon' ) * AJJ 00231 ELSE 00232 SSTOP = TOL 00233 END IF 00234 * 00235 * Set first half of WORK to zero, holds dot products 00236 * 00237 DO 120 I = 1, N 00238 WORK( I ) = 0 00239 120 CONTINUE 00240 * 00241 IF( UPPER ) THEN 00242 * 00243 * Compute the Cholesky factorization P**T * A * P = U**H * U 00244 * 00245 DO 150 J = 1, N 00246 * 00247 * Find pivot, test for exit, else swap rows and columns 00248 * Update dot products, compute possible pivots which are 00249 * stored in the second half of WORK 00250 * 00251 DO 130 I = J, N 00252 * 00253 IF( J.GT.1 ) THEN 00254 WORK( I ) = WORK( I ) + 00255 $ REAL( CONJG( A( J-1, I ) )* 00256 $ A( J-1, I ) ) 00257 END IF 00258 WORK( N+I ) = REAL( A( I, I ) ) - WORK( I ) 00259 * 00260 130 CONTINUE 00261 * 00262 IF( J.GT.1 ) THEN 00263 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 ) 00264 PVT = ITEMP + J - 1 00265 AJJ = WORK( N+PVT ) 00266 IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN 00267 A( J, J ) = AJJ 00268 GO TO 190 00269 END IF 00270 END IF 00271 * 00272 IF( J.NE.PVT ) THEN 00273 * 00274 * Pivot OK, so can now swap pivot rows and columns 00275 * 00276 A( PVT, PVT ) = A( J, J ) 00277 CALL CSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 ) 00278 IF( PVT.LT.N ) 00279 $ CALL CSWAP( N-PVT, A( J, PVT+1 ), LDA, 00280 $ A( PVT, PVT+1 ), LDA ) 00281 DO 140 I = J + 1, PVT - 1 00282 CTEMP = CONJG( A( J, I ) ) 00283 A( J, I ) = CONJG( A( I, PVT ) ) 00284 A( I, PVT ) = CTEMP 00285 140 CONTINUE 00286 A( J, PVT ) = CONJG( A( J, PVT ) ) 00287 * 00288 * Swap dot products and PIV 00289 * 00290 STEMP = WORK( J ) 00291 WORK( J ) = WORK( PVT ) 00292 WORK( PVT ) = STEMP 00293 ITEMP = PIV( PVT ) 00294 PIV( PVT ) = PIV( J ) 00295 PIV( J ) = ITEMP 00296 END IF 00297 * 00298 AJJ = SQRT( AJJ ) 00299 A( J, J ) = AJJ 00300 * 00301 * Compute elements J+1:N of row J 00302 * 00303 IF( J.LT.N ) THEN 00304 CALL CLACGV( J-1, A( 1, J ), 1 ) 00305 CALL CGEMV( 'Trans', J-1, N-J, -CONE, A( 1, J+1 ), LDA, 00306 $ A( 1, J ), 1, CONE, A( J, J+1 ), LDA ) 00307 CALL CLACGV( J-1, A( 1, J ), 1 ) 00308 CALL CSSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA ) 00309 END IF 00310 * 00311 150 CONTINUE 00312 * 00313 ELSE 00314 * 00315 * Compute the Cholesky factorization P**T * A * P = L * L**H 00316 * 00317 DO 180 J = 1, N 00318 * 00319 * Find pivot, test for exit, else swap rows and columns 00320 * Update dot products, compute possible pivots which are 00321 * stored in the second half of WORK 00322 * 00323 DO 160 I = J, N 00324 * 00325 IF( J.GT.1 ) THEN 00326 WORK( I ) = WORK( I ) + 00327 $ REAL( CONJG( A( I, J-1 ) )* 00328 $ A( I, J-1 ) ) 00329 END IF 00330 WORK( N+I ) = REAL( A( I, I ) ) - WORK( I ) 00331 * 00332 160 CONTINUE 00333 * 00334 IF( J.GT.1 ) THEN 00335 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 ) 00336 PVT = ITEMP + J - 1 00337 AJJ = WORK( N+PVT ) 00338 IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN 00339 A( J, J ) = AJJ 00340 GO TO 190 00341 END IF 00342 END IF 00343 * 00344 IF( J.NE.PVT ) THEN 00345 * 00346 * Pivot OK, so can now swap pivot rows and columns 00347 * 00348 A( PVT, PVT ) = A( J, J ) 00349 CALL CSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA ) 00350 IF( PVT.LT.N ) 00351 $ CALL CSWAP( N-PVT, A( PVT+1, J ), 1, A( PVT+1, PVT ), 00352 $ 1 ) 00353 DO 170 I = J + 1, PVT - 1 00354 CTEMP = CONJG( A( I, J ) ) 00355 A( I, J ) = CONJG( A( PVT, I ) ) 00356 A( PVT, I ) = CTEMP 00357 170 CONTINUE 00358 A( PVT, J ) = CONJG( A( PVT, J ) ) 00359 * 00360 * Swap dot products and PIV 00361 * 00362 STEMP = WORK( J ) 00363 WORK( J ) = WORK( PVT ) 00364 WORK( PVT ) = STEMP 00365 ITEMP = PIV( PVT ) 00366 PIV( PVT ) = PIV( J ) 00367 PIV( J ) = ITEMP 00368 END IF 00369 * 00370 AJJ = SQRT( AJJ ) 00371 A( J, J ) = AJJ 00372 * 00373 * Compute elements J+1:N of column J 00374 * 00375 IF( J.LT.N ) THEN 00376 CALL CLACGV( J-1, A( J, 1 ), LDA ) 00377 CALL CGEMV( 'No Trans', N-J, J-1, -CONE, A( J+1, 1 ), 00378 $ LDA, A( J, 1 ), LDA, CONE, A( J+1, J ), 1 ) 00379 CALL CLACGV( J-1, A( J, 1 ), LDA ) 00380 CALL CSSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 ) 00381 END IF 00382 * 00383 180 CONTINUE 00384 * 00385 END IF 00386 * 00387 * Ran to completion, A has full rank 00388 * 00389 RANK = N 00390 * 00391 GO TO 200 00392 190 CONTINUE 00393 * 00394 * Rank is number of steps completed. Set INFO = 1 to signal 00395 * that the factorization cannot be used to solve a system. 00396 * 00397 RANK = J - 1 00398 INFO = 1 00399 * 00400 200 CONTINUE 00401 RETURN 00402 * 00403 * End of CPSTF2 00404 * 00405 END