|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO ) 00002 * 00003 * -- LAPACK PROTOTYPE 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 * CPSTF2 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' * A * P = U' * U , if UPLO = 'U', 00026 * P' * A * P = L * L', 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 2 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 * PIV (output) INTEGER array, dimension (N) 00058 * PIV is such that the nonzero entries are P( PIV(K), K ) = 1. 00059 * 00060 * RANK (output) INTEGER 00061 * The rank of A given by the number of steps the algorithm 00062 * completed. 00063 * 00064 * TOL (input) REAL 00065 * User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) ) 00066 * will be used. The algorithm terminates at the (K-1)st step 00067 * if the pivot <= TOL. 00068 * 00069 * LDA (input) INTEGER 00070 * The leading dimension of the array A. LDA >= max(1,N). 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, PVT 00094 LOGICAL UPPER 00095 * .. 00096 * .. External Functions .. 00097 REAL SLAMCH 00098 LOGICAL LSAME, SISNAN 00099 EXTERNAL SLAMCH, LSAME, SISNAN 00100 * .. 00101 * .. External Subroutines .. 00102 EXTERNAL CGEMV, CLACGV, CSSCAL, CSWAP, XERBLA 00103 * .. 00104 * .. Intrinsic Functions .. 00105 INTRINSIC CONJG, MAX, REAL, SQRT 00106 * .. 00107 * .. Executable Statements .. 00108 * 00109 * Test the input parameters 00110 * 00111 INFO = 0 00112 UPPER = LSAME( UPLO, 'U' ) 00113 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00114 INFO = -1 00115 ELSE IF( N.LT.0 ) THEN 00116 INFO = -2 00117 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00118 INFO = -4 00119 END IF 00120 IF( INFO.NE.0 ) THEN 00121 CALL XERBLA( 'CPSTF2', -INFO ) 00122 RETURN 00123 END IF 00124 * 00125 * Quick return if possible 00126 * 00127 IF( N.EQ.0 ) 00128 $ RETURN 00129 * 00130 * Initialize PIV 00131 * 00132 DO 100 I = 1, N 00133 PIV( I ) = I 00134 100 CONTINUE 00135 * 00136 * Compute stopping value 00137 * 00138 DO 110 I = 1, N 00139 WORK( I ) = REAL( A( I, I ) ) 00140 110 CONTINUE 00141 PVT = MAXLOC( WORK( 1:N ), 1 ) 00142 AJJ = REAL ( A( PVT, PVT ) ) 00143 IF( AJJ.EQ.ZERO.OR.SISNAN( AJJ ) ) THEN 00144 RANK = 0 00145 INFO = 1 00146 GO TO 200 00147 END IF 00148 * 00149 * Compute stopping value if not supplied 00150 * 00151 IF( TOL.LT.ZERO ) THEN 00152 SSTOP = N * SLAMCH( 'Epsilon' ) * AJJ 00153 ELSE 00154 SSTOP = TOL 00155 END IF 00156 * 00157 * Set first half of WORK to zero, holds dot products 00158 * 00159 DO 120 I = 1, N 00160 WORK( I ) = 0 00161 120 CONTINUE 00162 * 00163 IF( UPPER ) THEN 00164 * 00165 * Compute the Cholesky factorization P' * A * P = U' * U 00166 * 00167 DO 150 J = 1, N 00168 * 00169 * Find pivot, test for exit, else swap rows and columns 00170 * Update dot products, compute possible pivots which are 00171 * stored in the second half of WORK 00172 * 00173 DO 130 I = J, N 00174 * 00175 IF( J.GT.1 ) THEN 00176 WORK( I ) = WORK( I ) + 00177 $ REAL( CONJG( A( J-1, I ) )* 00178 $ A( J-1, I ) ) 00179 END IF 00180 WORK( N+I ) = REAL( A( I, I ) ) - WORK( I ) 00181 * 00182 130 CONTINUE 00183 * 00184 IF( J.GT.1 ) THEN 00185 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 ) 00186 PVT = ITEMP + J - 1 00187 AJJ = WORK( N+PVT ) 00188 IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN 00189 A( J, J ) = AJJ 00190 GO TO 190 00191 END IF 00192 END IF 00193 * 00194 IF( J.NE.PVT ) THEN 00195 * 00196 * Pivot OK, so can now swap pivot rows and columns 00197 * 00198 A( PVT, PVT ) = A( J, J ) 00199 CALL CSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 ) 00200 IF( PVT.LT.N ) 00201 $ CALL CSWAP( N-PVT, A( J, PVT+1 ), LDA, 00202 $ A( PVT, PVT+1 ), LDA ) 00203 DO 140 I = J + 1, PVT - 1 00204 CTEMP = CONJG( A( J, I ) ) 00205 A( J, I ) = CONJG( A( I, PVT ) ) 00206 A( I, PVT ) = CTEMP 00207 140 CONTINUE 00208 A( J, PVT ) = CONJG( A( J, PVT ) ) 00209 * 00210 * Swap dot products and PIV 00211 * 00212 STEMP = WORK( J ) 00213 WORK( J ) = WORK( PVT ) 00214 WORK( PVT ) = STEMP 00215 ITEMP = PIV( PVT ) 00216 PIV( PVT ) = PIV( J ) 00217 PIV( J ) = ITEMP 00218 END IF 00219 * 00220 AJJ = SQRT( AJJ ) 00221 A( J, J ) = AJJ 00222 * 00223 * Compute elements J+1:N of row J 00224 * 00225 IF( J.LT.N ) THEN 00226 CALL CLACGV( J-1, A( 1, J ), 1 ) 00227 CALL CGEMV( 'Trans', J-1, N-J, -CONE, A( 1, J+1 ), LDA, 00228 $ A( 1, J ), 1, CONE, A( J, J+1 ), LDA ) 00229 CALL CLACGV( J-1, A( 1, J ), 1 ) 00230 CALL CSSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA ) 00231 END IF 00232 * 00233 150 CONTINUE 00234 * 00235 ELSE 00236 * 00237 * Compute the Cholesky factorization P' * A * P = L * L' 00238 * 00239 DO 180 J = 1, N 00240 * 00241 * Find pivot, test for exit, else swap rows and columns 00242 * Update dot products, compute possible pivots which are 00243 * stored in the second half of WORK 00244 * 00245 DO 160 I = J, N 00246 * 00247 IF( J.GT.1 ) THEN 00248 WORK( I ) = WORK( I ) + 00249 $ REAL( CONJG( A( I, J-1 ) )* 00250 $ A( I, J-1 ) ) 00251 END IF 00252 WORK( N+I ) = REAL( A( I, I ) ) - WORK( I ) 00253 * 00254 160 CONTINUE 00255 * 00256 IF( J.GT.1 ) THEN 00257 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 ) 00258 PVT = ITEMP + J - 1 00259 AJJ = WORK( N+PVT ) 00260 IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN 00261 A( J, J ) = AJJ 00262 GO TO 190 00263 END IF 00264 END IF 00265 * 00266 IF( J.NE.PVT ) THEN 00267 * 00268 * Pivot OK, so can now swap pivot rows and columns 00269 * 00270 A( PVT, PVT ) = A( J, J ) 00271 CALL CSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA ) 00272 IF( PVT.LT.N ) 00273 $ CALL CSWAP( N-PVT, A( PVT+1, J ), 1, A( PVT+1, PVT ), 00274 $ 1 ) 00275 DO 170 I = J + 1, PVT - 1 00276 CTEMP = CONJG( A( I, J ) ) 00277 A( I, J ) = CONJG( A( PVT, I ) ) 00278 A( PVT, I ) = CTEMP 00279 170 CONTINUE 00280 A( PVT, J ) = CONJG( A( PVT, J ) ) 00281 * 00282 * Swap dot products and PIV 00283 * 00284 STEMP = WORK( J ) 00285 WORK( J ) = WORK( PVT ) 00286 WORK( PVT ) = STEMP 00287 ITEMP = PIV( PVT ) 00288 PIV( PVT ) = PIV( J ) 00289 PIV( J ) = ITEMP 00290 END IF 00291 * 00292 AJJ = SQRT( AJJ ) 00293 A( J, J ) = AJJ 00294 * 00295 * Compute elements J+1:N of column J 00296 * 00297 IF( J.LT.N ) THEN 00298 CALL CLACGV( J-1, A( J, 1 ), LDA ) 00299 CALL CGEMV( 'No Trans', N-J, J-1, -CONE, A( J+1, 1 ), 00300 $ LDA, A( J, 1 ), LDA, CONE, A( J+1, J ), 1 ) 00301 CALL CLACGV( J-1, A( J, 1 ), LDA ) 00302 CALL CSSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 ) 00303 END IF 00304 * 00305 180 CONTINUE 00306 * 00307 END IF 00308 * 00309 * Ran to completion, A has full rank 00310 * 00311 RANK = N 00312 * 00313 GO TO 200 00314 190 CONTINUE 00315 * 00316 * Rank is number of steps completed. Set INFO = 1 to signal 00317 * that the factorization cannot be used to solve a system. 00318 * 00319 RANK = J - 1 00320 INFO = 1 00321 * 00322 200 CONTINUE 00323 RETURN 00324 * 00325 * End of CPSTF2 00326 * 00327 END
1.7.3 
