00001 /* dsgesv.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static doublereal c_b10 = -1.; 00019 static doublereal c_b11 = 1.; 00020 static integer c__1 = 1; 00021 00022 /* Subroutine */ int dsgesv_(integer *n, integer *nrhs, doublereal *a, 00023 integer *lda, integer *ipiv, doublereal *b, integer *ldb, doublereal * 00024 x, integer *ldx, doublereal *work, real *swork, integer *iter, 00025 integer *info) 00026 { 00027 /* System generated locals */ 00028 integer a_dim1, a_offset, b_dim1, b_offset, work_dim1, work_offset, 00029 x_dim1, x_offset, i__1; 00030 doublereal d__1; 00031 00032 /* Builtin functions */ 00033 double sqrt(doublereal); 00034 00035 /* Local variables */ 00036 integer i__; 00037 doublereal cte, eps, anrm; 00038 integer ptsa; 00039 doublereal rnrm, xnrm; 00040 integer ptsx; 00041 extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *, 00042 integer *, doublereal *, doublereal *, integer *, doublereal *, 00043 integer *, doublereal *, doublereal *, integer *); 00044 integer iiter; 00045 extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 00046 integer *, doublereal *, integer *), dlag2s_(integer *, integer *, 00047 doublereal *, integer *, real *, integer *, integer *), slag2d_( 00048 integer *, integer *, real *, integer *, doublereal *, integer *, 00049 integer *); 00050 extern doublereal dlamch_(char *), dlange_(char *, integer *, 00051 integer *, doublereal *, integer *, doublereal *); 00052 extern integer idamax_(integer *, doublereal *, integer *); 00053 extern /* Subroutine */ int dlacpy_(char *, integer *, integer *, 00054 doublereal *, integer *, doublereal *, integer *), 00055 xerbla_(char *, integer *), dgetrf_(integer *, integer *, 00056 doublereal *, integer *, integer *, integer *), dgetrs_(char *, 00057 integer *, integer *, doublereal *, integer *, integer *, 00058 doublereal *, integer *, integer *), sgetrf_(integer *, 00059 integer *, real *, integer *, integer *, integer *), sgetrs_(char 00060 *, integer *, integer *, real *, integer *, integer *, real *, 00061 integer *, integer *); 00062 00063 00064 /* -- LAPACK PROTOTYPE driver routine (version 3.2) -- */ 00065 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00066 /* February 2007 */ 00067 00068 /* .. */ 00069 /* .. Scalar Arguments .. */ 00070 /* .. */ 00071 /* .. Array Arguments .. */ 00072 /* .. */ 00073 00074 /* Purpose */ 00075 /* ======= */ 00076 00077 /* DSGESV computes the solution to a real system of linear equations */ 00078 /* A * X = B, */ 00079 /* where A is an N-by-N matrix and X and B are N-by-NRHS matrices. */ 00080 00081 /* DSGESV first attempts to factorize the matrix in SINGLE PRECISION */ 00082 /* and use this factorization within an iterative refinement procedure */ 00083 /* to produce a solution with DOUBLE PRECISION normwise backward error */ 00084 /* quality (see below). If the approach fails the method switches to a */ 00085 /* DOUBLE PRECISION factorization and solve. */ 00086 00087 /* The iterative refinement is not going to be a winning strategy if */ 00088 /* the ratio SINGLE PRECISION performance over DOUBLE PRECISION */ 00089 /* performance is too small. A reasonable strategy should take the */ 00090 /* number of right-hand sides and the size of the matrix into account. */ 00091 /* This might be done with a call to ILAENV in the future. Up to now, we */ 00092 /* always try iterative refinement. */ 00093 00094 /* The iterative refinement process is stopped if */ 00095 /* ITER > ITERMAX */ 00096 /* or for all the RHS we have: */ 00097 /* RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX */ 00098 /* where */ 00099 /* o ITER is the number of the current iteration in the iterative */ 00100 /* refinement process */ 00101 /* o RNRM is the infinity-norm of the residual */ 00102 /* o XNRM is the infinity-norm of the solution */ 00103 /* o ANRM is the infinity-operator-norm of the matrix A */ 00104 /* o EPS is the machine epsilon returned by DLAMCH('Epsilon') */ 00105 /* The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00 */ 00106 /* respectively. */ 00107 00108 /* Arguments */ 00109 /* ========= */ 00110 00111 /* N (input) INTEGER */ 00112 /* The number of linear equations, i.e., the order of the */ 00113 /* matrix A. N >= 0. */ 00114 00115 /* NRHS (input) INTEGER */ 00116 /* The number of right hand sides, i.e., the number of columns */ 00117 /* of the matrix B. NRHS >= 0. */ 00118 00119 /* A (input or input/ouptut) DOUBLE PRECISION array, */ 00120 /* dimension (LDA,N) */ 00121 /* On entry, the N-by-N coefficient matrix A. */ 00122 /* On exit, if iterative refinement has been successfully used */ 00123 /* (INFO.EQ.0 and ITER.GE.0, see description below), then A is */ 00124 /* unchanged, if double precision factorization has been used */ 00125 /* (INFO.EQ.0 and ITER.LT.0, see description below), then the */ 00126 /* array A contains the factors L and U from the factorization */ 00127 /* A = P*L*U; the unit diagonal elements of L are not stored. */ 00128 00129 /* LDA (input) INTEGER */ 00130 /* The leading dimension of the array A. LDA >= max(1,N). */ 00131 00132 /* IPIV (output) INTEGER array, dimension (N) */ 00133 /* The pivot indices that define the permutation matrix P; */ 00134 /* row i of the matrix was interchanged with row IPIV(i). */ 00135 /* Corresponds either to the single precision factorization */ 00136 /* (if INFO.EQ.0 and ITER.GE.0) or the double precision */ 00137 /* factorization (if INFO.EQ.0 and ITER.LT.0). */ 00138 00139 /* B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) */ 00140 /* The N-by-NRHS right hand side matrix B. */ 00141 00142 /* LDB (input) INTEGER */ 00143 /* The leading dimension of the array B. LDB >= max(1,N). */ 00144 00145 /* X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) */ 00146 /* If INFO = 0, the N-by-NRHS solution matrix X. */ 00147 00148 /* LDX (input) INTEGER */ 00149 /* The leading dimension of the array X. LDX >= max(1,N). */ 00150 00151 /* WORK (workspace) DOUBLE PRECISION array, dimension (N*NRHS) */ 00152 /* This array is used to hold the residual vectors. */ 00153 00154 /* SWORK (workspace) 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 00158 /* ITER (output) INTEGER */ 00159 /* < 0: iterative refinement has failed, double precision */ 00160 /* factorization has been performed */ 00161 /* -1 : the routine fell back to full precision for */ 00162 /* implementation- or machine-specific reasons */ 00163 /* -2 : narrowing the precision induced an overflow, */ 00164 /* the routine fell back to full precision */ 00165 /* -3 : failure of SGETRF */ 00166 /* -31: stop the iterative refinement after the 30th */ 00167 /* iterations */ 00168 /* > 0: iterative refinement has been sucessfully used. */ 00169 /* Returns the number of iterations */ 00170 00171 /* INFO (output) INTEGER */ 00172 /* = 0: successful exit */ 00173 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00174 /* > 0: if INFO = i, U(i,i) computed in DOUBLE PRECISION is */ 00175 /* exactly zero. The factorization has been completed, */ 00176 /* but the factor U is exactly singular, so the solution */ 00177 /* could not be computed. */ 00178 00179 /* ========= */ 00180 00181 /* .. Parameters .. */ 00182 00183 00184 00185 00186 /* .. Local Scalars .. */ 00187 00188 /* .. External Subroutines .. */ 00189 /* .. */ 00190 /* .. External Functions .. */ 00191 /* .. */ 00192 /* .. Intrinsic Functions .. */ 00193 /* .. */ 00194 /* .. Executable Statements .. */ 00195 00196 /* Parameter adjustments */ 00197 work_dim1 = *n; 00198 work_offset = 1 + work_dim1; 00199 work -= work_offset; 00200 a_dim1 = *lda; 00201 a_offset = 1 + a_dim1; 00202 a -= a_offset; 00203 --ipiv; 00204 b_dim1 = *ldb; 00205 b_offset = 1 + b_dim1; 00206 b -= b_offset; 00207 x_dim1 = *ldx; 00208 x_offset = 1 + x_dim1; 00209 x -= x_offset; 00210 --swork; 00211 00212 /* Function Body */ 00213 *info = 0; 00214 *iter = 0; 00215 00216 /* Test the input parameters. */ 00217 00218 if (*n < 0) { 00219 *info = -1; 00220 } else if (*nrhs < 0) { 00221 *info = -2; 00222 } else if (*lda < max(1,*n)) { 00223 *info = -4; 00224 } else if (*ldb < max(1,*n)) { 00225 *info = -7; 00226 } else if (*ldx < max(1,*n)) { 00227 *info = -9; 00228 } 00229 if (*info != 0) { 00230 i__1 = -(*info); 00231 xerbla_("DSGESV", &i__1); 00232 return 0; 00233 } 00234 00235 /* Quick return if (N.EQ.0). */ 00236 00237 if (*n == 0) { 00238 return 0; 00239 } 00240 00241 /* Skip single precision iterative refinement if a priori slower */ 00242 /* than double precision factorization. */ 00243 00244 if (FALSE_) { 00245 *iter = -1; 00246 goto L40; 00247 } 00248 00249 /* Compute some constants. */ 00250 00251 anrm = dlange_("I", n, n, &a[a_offset], lda, &work[work_offset]); 00252 eps = dlamch_("Epsilon"); 00253 cte = anrm * eps * sqrt((doublereal) (*n)) * 1.; 00254 00255 /* Set the indices PTSA, PTSX for referencing SA and SX in SWORK. */ 00256 00257 ptsa = 1; 00258 ptsx = ptsa + *n * *n; 00259 00260 /* Convert B from double precision to single precision and store the */ 00261 /* result in SX. */ 00262 00263 dlag2s_(n, nrhs, &b[b_offset], ldb, &swork[ptsx], n, info); 00264 00265 if (*info != 0) { 00266 *iter = -2; 00267 goto L40; 00268 } 00269 00270 /* Convert A from double precision to single precision and store the */ 00271 /* result in SA. */ 00272 00273 dlag2s_(n, n, &a[a_offset], lda, &swork[ptsa], n, info); 00274 00275 if (*info != 0) { 00276 *iter = -2; 00277 goto L40; 00278 } 00279 00280 /* Compute the LU factorization of SA. */ 00281 00282 sgetrf_(n, n, &swork[ptsa], n, &ipiv[1], info); 00283 00284 if (*info != 0) { 00285 *iter = -3; 00286 goto L40; 00287 } 00288 00289 /* Solve the system SA*SX = SB. */ 00290 00291 sgetrs_("No transpose", n, nrhs, &swork[ptsa], n, &ipiv[1], &swork[ptsx], 00292 n, info); 00293 00294 /* Convert SX back to double precision */ 00295 00296 slag2d_(n, nrhs, &swork[ptsx], n, &x[x_offset], ldx, info); 00297 00298 /* Compute R = B - AX (R is WORK). */ 00299 00300 dlacpy_("All", n, nrhs, &b[b_offset], ldb, &work[work_offset], n); 00301 00302 dgemm_("No Transpose", "No Transpose", n, nrhs, n, &c_b10, &a[a_offset], 00303 lda, &x[x_offset], ldx, &c_b11, &work[work_offset], n); 00304 00305 /* Check whether the NRHS normwise backward errors satisfy the */ 00306 /* stopping criterion. If yes, set ITER=0 and return. */ 00307 00308 i__1 = *nrhs; 00309 for (i__ = 1; i__ <= i__1; ++i__) { 00310 xnrm = (d__1 = x[idamax_(n, &x[i__ * x_dim1 + 1], &c__1) + i__ * 00311 x_dim1], abs(d__1)); 00312 rnrm = (d__1 = work[idamax_(n, &work[i__ * work_dim1 + 1], &c__1) + 00313 i__ * work_dim1], abs(d__1)); 00314 if (rnrm > xnrm * cte) { 00315 goto L10; 00316 } 00317 } 00318 00319 /* If we are here, the NRHS normwise backward errors satisfy the */ 00320 /* stopping criterion. We are good to exit. */ 00321 00322 *iter = 0; 00323 return 0; 00324 00325 L10: 00326 00327 for (iiter = 1; iiter <= 30; ++iiter) { 00328 00329 /* Convert R (in WORK) from double precision to single precision */ 00330 /* and store the result in SX. */ 00331 00332 dlag2s_(n, nrhs, &work[work_offset], n, &swork[ptsx], n, info); 00333 00334 if (*info != 0) { 00335 *iter = -2; 00336 goto L40; 00337 } 00338 00339 /* Solve the system SA*SX = SR. */ 00340 00341 sgetrs_("No transpose", n, nrhs, &swork[ptsa], n, &ipiv[1], &swork[ 00342 ptsx], n, info); 00343 00344 /* Convert SX back to double precision and update the current */ 00345 /* iterate. */ 00346 00347 slag2d_(n, nrhs, &swork[ptsx], n, &work[work_offset], n, info); 00348 00349 i__1 = *nrhs; 00350 for (i__ = 1; i__ <= i__1; ++i__) { 00351 daxpy_(n, &c_b11, &work[i__ * work_dim1 + 1], &c__1, &x[i__ * 00352 x_dim1 + 1], &c__1); 00353 } 00354 00355 /* Compute R = B - AX (R is WORK). */ 00356 00357 dlacpy_("All", n, nrhs, &b[b_offset], ldb, &work[work_offset], n); 00358 00359 dgemm_("No Transpose", "No Transpose", n, nrhs, n, &c_b10, &a[ 00360 a_offset], lda, &x[x_offset], ldx, &c_b11, &work[work_offset], 00361 n); 00362 00363 /* Check whether the NRHS normwise backward errors satisfy the */ 00364 /* stopping criterion. If yes, set ITER=IITER>0 and return. */ 00365 00366 i__1 = *nrhs; 00367 for (i__ = 1; i__ <= i__1; ++i__) { 00368 xnrm = (d__1 = x[idamax_(n, &x[i__ * x_dim1 + 1], &c__1) + i__ * 00369 x_dim1], abs(d__1)); 00370 rnrm = (d__1 = work[idamax_(n, &work[i__ * work_dim1 + 1], &c__1) 00371 + i__ * work_dim1], abs(d__1)); 00372 if (rnrm > xnrm * cte) { 00373 goto L20; 00374 } 00375 } 00376 00377 /* If we are here, the NRHS normwise backward errors satisfy the */ 00378 /* stopping criterion, we are good to exit. */ 00379 00380 *iter = iiter; 00381 00382 return 0; 00383 00384 L20: 00385 00386 /* L30: */ 00387 ; 00388 } 00389 00390 /* If we are at this place of the code, this is because we have */ 00391 /* performed ITER=ITERMAX iterations and never satisified the */ 00392 /* stopping criterion, set up the ITER flag accordingly and follow up */ 00393 /* on double precision routine. */ 00394 00395 *iter = -31; 00396 00397 L40: 00398 00399 /* Single-precision iterative refinement failed to converge to a */ 00400 /* satisfactory solution, so we resort to double precision. */ 00401 00402 dgetrf_(n, n, &a[a_offset], lda, &ipiv[1], info); 00403 00404 if (*info != 0) { 00405 return 0; 00406 } 00407 00408 dlacpy_("All", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx); 00409 dgetrs_("No transpose", n, nrhs, &a[a_offset], lda, &ipiv[1], &x[x_offset] 00410 , ldx, info); 00411 00412 return 0; 00413 00414 /* End of DSGESV. */ 00415 00416 } /* dsgesv_ */