![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> DSPOSV computes the solution to system of linear equations A * X = B for PO matrices</b> 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DSPOSV + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsposv.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsposv.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsposv.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DSPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK, 00022 * SWORK, ITER, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER UPLO 00026 * INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS 00027 * .. 00028 * .. Array Arguments .. 00029 * REAL SWORK( * ) 00030 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( N, * ), 00031 * $ X( LDX, * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> DSPOSV computes the solution to a real system of linear equations 00041 *> A * X = B, 00042 *> where A is an N-by-N symmetric positive definite matrix and X and B 00043 *> are N-by-NRHS matrices. 00044 *> 00045 *> DSPOSV first attempts to factorize the matrix in SINGLE PRECISION 00046 *> and use this factorization within an iterative refinement procedure 00047 *> to produce a solution with DOUBLE PRECISION normwise backward error 00048 *> quality (see below). If the approach fails the method switches to a 00049 *> DOUBLE PRECISION factorization and solve. 00050 *> 00051 *> The iterative refinement is not going to be a winning strategy if 00052 *> the ratio SINGLE PRECISION performance over DOUBLE PRECISION 00053 *> performance is too small. A reasonable strategy should take the 00054 *> number of right-hand sides and the size of the matrix into account. 00055 *> This might be done with a call to ILAENV in the future. Up to now, we 00056 *> always try iterative refinement. 00057 *> 00058 *> The iterative refinement process is stopped if 00059 *> ITER > ITERMAX 00060 *> or for all the RHS we have: 00061 *> RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX 00062 *> where 00063 *> o ITER is the number of the current iteration in the iterative 00064 *> refinement process 00065 *> o RNRM is the infinity-norm of the residual 00066 *> o XNRM is the infinity-norm of the solution 00067 *> o ANRM is the infinity-operator-norm of the matrix A 00068 *> o EPS is the machine epsilon returned by DLAMCH('Epsilon') 00069 *> The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00 00070 *> respectively. 00071 *> \endverbatim 00072 * 00073 * Arguments: 00074 * ========== 00075 * 00076 *> \param[in] UPLO 00077 *> \verbatim 00078 *> UPLO is CHARACTER*1 00079 *> = 'U': Upper triangle of A is stored; 00080 *> = 'L': Lower triangle of A is stored. 00081 *> \endverbatim 00082 *> 00083 *> \param[in] N 00084 *> \verbatim 00085 *> N is INTEGER 00086 *> The number of linear equations, i.e., the order of the 00087 *> matrix A. N >= 0. 00088 *> \endverbatim 00089 *> 00090 *> \param[in] NRHS 00091 *> \verbatim 00092 *> NRHS is INTEGER 00093 *> The number of right hand sides, i.e., the number of columns 00094 *> of the matrix B. NRHS >= 0. 00095 *> \endverbatim 00096 *> 00097 *> \param[in,out] A 00098 *> \verbatim 00099 *> A is DOUBLE PRECISION array, 00100 *> dimension (LDA,N) 00101 *> On entry, the symmetric matrix A. If UPLO = 'U', the leading 00102 *> N-by-N upper triangular part of A contains the upper 00103 *> triangular part of the matrix A, and the strictly lower 00104 *> triangular part of A is not referenced. If UPLO = 'L', the 00105 *> leading N-by-N lower triangular part of A contains the lower 00106 *> triangular part of the matrix A, and the strictly upper 00107 *> triangular part of A is not referenced. 00108 *> On exit, if iterative refinement has been successfully used 00109 *> (INFO.EQ.0 and ITER.GE.0, see description below), then A is 00110 *> unchanged, if double precision factorization has been used 00111 *> (INFO.EQ.0 and ITER.LT.0, see description below), then the 00112 *> array A contains the factor U or L from the Cholesky 00113 *> factorization A = U**T*U or A = L*L**T. 00114 *> \endverbatim 00115 *> 00116 *> \param[in] LDA 00117 *> \verbatim 00118 *> LDA is INTEGER 00119 *> The leading dimension of the array A. LDA >= max(1,N). 00120 *> \endverbatim 00121 *> 00122 *> \param[in] B 00123 *> \verbatim 00124 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS) 00125 *> The N-by-NRHS right hand side matrix B. 00126 *> \endverbatim 00127 *> 00128 *> \param[in] LDB 00129 *> \verbatim 00130 *> LDB is INTEGER 00131 *> The leading dimension of the array B. LDB >= max(1,N). 00132 *> \endverbatim 00133 *> 00134 *> \param[out] X 00135 *> \verbatim 00136 *> X is DOUBLE PRECISION array, dimension (LDX,NRHS) 00137 *> If INFO = 0, the N-by-NRHS solution matrix X. 00138 *> \endverbatim 00139 *> 00140 *> \param[in] LDX 00141 *> \verbatim 00142 *> LDX is INTEGER 00143 *> The leading dimension of the array X. LDX >= max(1,N). 00144 *> \endverbatim 00145 *> 00146 *> \param[out] WORK 00147 *> \verbatim 00148 *> WORK is DOUBLE PRECISION array, dimension (N,NRHS) 00149 *> This array is used to hold the residual vectors. 00150 *> \endverbatim 00151 *> 00152 *> \param[out] SWORK 00153 *> \verbatim 00154 *> SWORK is REAL array, dimension (N*(N+NRHS)) 00155 *> This array is used to use the single precision matrix and the 00156 *> right-hand sides or solutions in single precision. 00157 *> \endverbatim 00158 *> 00159 *> \param[out] ITER 00160 *> \verbatim 00161 *> ITER is INTEGER 00162 *> < 0: iterative refinement has failed, double precision 00163 *> factorization has been performed 00164 *> -1 : the routine fell back to full precision for 00165 *> implementation- or machine-specific reasons 00166 *> -2 : narrowing the precision induced an overflow, 00167 *> the routine fell back to full precision 00168 *> -3 : failure of SPOTRF 00169 *> -31: stop the iterative refinement after the 30th 00170 *> iterations 00171 *> > 0: iterative refinement has been sucessfully used. 00172 *> Returns the number of iterations 00173 *> \endverbatim 00174 *> 00175 *> \param[out] INFO 00176 *> \verbatim 00177 *> INFO is INTEGER 00178 *> = 0: successful exit 00179 *> < 0: if INFO = -i, the i-th argument had an illegal value 00180 *> > 0: if INFO = i, the leading minor of order i of (DOUBLE 00181 *> PRECISION) A is not positive definite, so the 00182 *> factorization could not be completed, and the solution 00183 *> has not been computed. 00184 *> \endverbatim 00185 * 00186 * Authors: 00187 * ======== 00188 * 00189 *> \author Univ. of Tennessee 00190 *> \author Univ. of California Berkeley 00191 *> \author Univ. of Colorado Denver 00192 *> \author NAG Ltd. 00193 * 00194 *> \date November 2011 00195 * 00196 *> \ingroup doublePOsolve 00197 * 00198 * ===================================================================== 00199 SUBROUTINE DSPOSV( UPLO, N, NRHS, A, LDA, B, LDB, X, LDX, WORK, 00200 $ SWORK, ITER, INFO ) 00201 * 00202 * -- LAPACK driver routine (version 3.4.0) -- 00203 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00204 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00205 * November 2011 00206 * 00207 * .. Scalar Arguments .. 00208 CHARACTER UPLO 00209 INTEGER INFO, ITER, LDA, LDB, LDX, N, NRHS 00210 * .. 00211 * .. Array Arguments .. 00212 REAL SWORK( * ) 00213 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( N, * ), 00214 $ X( LDX, * ) 00215 * .. 00216 * 00217 * ===================================================================== 00218 * 00219 * .. Parameters .. 00220 LOGICAL DOITREF 00221 PARAMETER ( DOITREF = .TRUE. ) 00222 * 00223 INTEGER ITERMAX 00224 PARAMETER ( ITERMAX = 30 ) 00225 * 00226 DOUBLE PRECISION BWDMAX 00227 PARAMETER ( BWDMAX = 1.0E+00 ) 00228 * 00229 DOUBLE PRECISION NEGONE, ONE 00230 PARAMETER ( NEGONE = -1.0D+0, ONE = 1.0D+0 ) 00231 * 00232 * .. Local Scalars .. 00233 INTEGER I, IITER, PTSA, PTSX 00234 DOUBLE PRECISION ANRM, CTE, EPS, RNRM, XNRM 00235 * 00236 * .. External Subroutines .. 00237 EXTERNAL DAXPY, DSYMM, DLACPY, DLAT2S, DLAG2S, SLAG2D, 00238 $ SPOTRF, SPOTRS, XERBLA 00239 * .. 00240 * .. External Functions .. 00241 INTEGER IDAMAX 00242 DOUBLE PRECISION DLAMCH, DLANSY 00243 LOGICAL LSAME 00244 EXTERNAL IDAMAX, DLAMCH, DLANSY, LSAME 00245 * .. 00246 * .. Intrinsic Functions .. 00247 INTRINSIC ABS, DBLE, MAX, SQRT 00248 * .. 00249 * .. Executable Statements .. 00250 * 00251 INFO = 0 00252 ITER = 0 00253 * 00254 * Test the input parameters. 00255 * 00256 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00257 INFO = -1 00258 ELSE IF( N.LT.0 ) THEN 00259 INFO = -2 00260 ELSE IF( NRHS.LT.0 ) THEN 00261 INFO = -3 00262 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00263 INFO = -5 00264 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00265 INFO = -7 00266 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN 00267 INFO = -9 00268 END IF 00269 IF( INFO.NE.0 ) THEN 00270 CALL XERBLA( 'DSPOSV', -INFO ) 00271 RETURN 00272 END IF 00273 * 00274 * Quick return if (N.EQ.0). 00275 * 00276 IF( N.EQ.0 ) 00277 $ RETURN 00278 * 00279 * Skip single precision iterative refinement if a priori slower 00280 * than double precision factorization. 00281 * 00282 IF( .NOT.DOITREF ) THEN 00283 ITER = -1 00284 GO TO 40 00285 END IF 00286 * 00287 * Compute some constants. 00288 * 00289 ANRM = DLANSY( 'I', UPLO, N, A, LDA, WORK ) 00290 EPS = DLAMCH( 'Epsilon' ) 00291 CTE = ANRM*EPS*SQRT( DBLE( N ) )*BWDMAX 00292 * 00293 * Set the indices PTSA, PTSX for referencing SA and SX in SWORK. 00294 * 00295 PTSA = 1 00296 PTSX = PTSA + N*N 00297 * 00298 * Convert B from double precision to single precision and store the 00299 * result in SX. 00300 * 00301 CALL DLAG2S( N, NRHS, B, LDB, SWORK( PTSX ), N, INFO ) 00302 * 00303 IF( INFO.NE.0 ) THEN 00304 ITER = -2 00305 GO TO 40 00306 END IF 00307 * 00308 * Convert A from double precision to single precision and store the 00309 * result in SA. 00310 * 00311 CALL DLAT2S( UPLO, N, A, LDA, SWORK( PTSA ), N, INFO ) 00312 * 00313 IF( INFO.NE.0 ) THEN 00314 ITER = -2 00315 GO TO 40 00316 END IF 00317 * 00318 * Compute the Cholesky factorization of SA. 00319 * 00320 CALL SPOTRF( UPLO, N, SWORK( PTSA ), N, INFO ) 00321 * 00322 IF( INFO.NE.0 ) THEN 00323 ITER = -3 00324 GO TO 40 00325 END IF 00326 * 00327 * Solve the system SA*SX = SB. 00328 * 00329 CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N, 00330 $ INFO ) 00331 * 00332 * Convert SX back to double precision 00333 * 00334 CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, X, LDX, INFO ) 00335 * 00336 * Compute R = B - AX (R is WORK). 00337 * 00338 CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N ) 00339 * 00340 CALL DSYMM( 'Left', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE, 00341 $ WORK, N ) 00342 * 00343 * Check whether the NRHS normwise backward errors satisfy the 00344 * stopping criterion. If yes, set ITER=0 and return. 00345 * 00346 DO I = 1, NRHS 00347 XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) ) 00348 RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) ) 00349 IF( RNRM.GT.XNRM*CTE ) 00350 $ GO TO 10 00351 END DO 00352 * 00353 * If we are here, the NRHS normwise backward errors satisfy the 00354 * stopping criterion. We are good to exit. 00355 * 00356 ITER = 0 00357 RETURN 00358 * 00359 10 CONTINUE 00360 * 00361 DO 30 IITER = 1, ITERMAX 00362 * 00363 * Convert R (in WORK) from double precision to single precision 00364 * and store the result in SX. 00365 * 00366 CALL DLAG2S( N, NRHS, WORK, N, SWORK( PTSX ), N, INFO ) 00367 * 00368 IF( INFO.NE.0 ) THEN 00369 ITER = -2 00370 GO TO 40 00371 END IF 00372 * 00373 * Solve the system SA*SX = SR. 00374 * 00375 CALL SPOTRS( UPLO, N, NRHS, SWORK( PTSA ), N, SWORK( PTSX ), N, 00376 $ INFO ) 00377 * 00378 * Convert SX back to double precision and update the current 00379 * iterate. 00380 * 00381 CALL SLAG2D( N, NRHS, SWORK( PTSX ), N, WORK, N, INFO ) 00382 * 00383 DO I = 1, NRHS 00384 CALL DAXPY( N, ONE, WORK( 1, I ), 1, X( 1, I ), 1 ) 00385 END DO 00386 * 00387 * Compute R = B - AX (R is WORK). 00388 * 00389 CALL DLACPY( 'All', N, NRHS, B, LDB, WORK, N ) 00390 * 00391 CALL DSYMM( 'L', UPLO, N, NRHS, NEGONE, A, LDA, X, LDX, ONE, 00392 $ WORK, N ) 00393 * 00394 * Check whether the NRHS normwise backward errors satisfy the 00395 * stopping criterion. If yes, set ITER=IITER>0 and return. 00396 * 00397 DO I = 1, NRHS 00398 XNRM = ABS( X( IDAMAX( N, X( 1, I ), 1 ), I ) ) 00399 RNRM = ABS( WORK( IDAMAX( N, WORK( 1, I ), 1 ), I ) ) 00400 IF( RNRM.GT.XNRM*CTE ) 00401 $ GO TO 20 00402 END DO 00403 * 00404 * If we are here, the NRHS normwise backward errors satisfy the 00405 * stopping criterion, we are good to exit. 00406 * 00407 ITER = IITER 00408 * 00409 RETURN 00410 * 00411 20 CONTINUE 00412 * 00413 30 CONTINUE 00414 * 00415 * If we are at this place of the code, this is because we have 00416 * performed ITER=ITERMAX iterations and never satisified the 00417 * stopping criterion, set up the ITER flag accordingly and follow 00418 * up on double precision routine. 00419 * 00420 ITER = -ITERMAX - 1 00421 * 00422 40 CONTINUE 00423 * 00424 * Single-precision iterative refinement failed to converge to a 00425 * satisfactory solution, so we resort to double precision. 00426 * 00427 CALL DPOTRF( UPLO, N, A, LDA, INFO ) 00428 * 00429 IF( INFO.NE.0 ) 00430 $ RETURN 00431 * 00432 CALL DLACPY( 'All', N, NRHS, B, LDB, X, LDX ) 00433 CALL DPOTRS( UPLO, N, NRHS, A, LDA, X, LDX, INFO ) 00434 * 00435 RETURN 00436 * 00437 * End of DSPOSV. 00438 * 00439 END