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