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