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