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