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