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