cgerfsx.c
Go to the documentation of this file.
00001 /* cgerfsx.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 logical c_true = TRUE_;
00019 static logical c_false = FALSE_;
00020 
00021 /* Subroutine */ int cgerfsx_(char *trans, char *equed, integer *n, integer *
00022         nrhs, complex *a, integer *lda, complex *af, integer *ldaf, integer *
00023         ipiv, real *r__, real *c__, complex *b, integer *ldb, complex *x, 
00024         integer *ldx, real *rcond, real *berr, integer *n_err_bnds__, real *
00025         err_bnds_norm__, real *err_bnds_comp__, integer *nparams, real *
00026         params, complex *work, real *rwork, integer *info)
00027 {
00028     /* System generated locals */
00029     integer a_dim1, a_offset, af_dim1, af_offset, b_dim1, b_offset, x_dim1, 
00030             x_offset, err_bnds_norm_dim1, err_bnds_norm_offset, 
00031             err_bnds_comp_dim1, err_bnds_comp_offset, i__1;
00032     real r__1, r__2;
00033 
00034     /* Builtin functions */
00035     double sqrt(doublereal);
00036 
00037     /* Local variables */
00038     real illrcond_thresh__, unstable_thresh__, err_lbnd__;
00039     integer ref_type__;
00040     extern integer ilatrans_(char *);
00041     integer j;
00042     real rcond_tmp__;
00043     integer prec_type__, trans_type__;
00044     real cwise_wrong__;
00045     extern /* Subroutine */ int cla_gerfsx_extended__(integer *, integer *, 
00046             integer *, integer *, complex *, integer *, complex *, integer *, 
00047             integer *, logical *, real *, complex *, integer *, complex *, 
00048             integer *, real *, integer *, real *, real *, complex *, real *, 
00049             complex *, complex *, real *, integer *, real *, real *, logical *
00050             , integer *);
00051     char norm[1];
00052     logical ignore_cwise__;
00053     extern doublereal cla_gercond_c__(char *, integer *, complex *, integer *,
00054              complex *, integer *, integer *, real *, logical *, integer *, 
00055             complex *, real *, ftnlen);
00056     extern logical lsame_(char *, char *);
00057     real anorm;
00058     extern doublereal cla_gercond_x__(char *, integer *, complex *, integer *,
00059              complex *, integer *, integer *, complex *, integer *, complex *,
00060              real *, ftnlen), clange_(char *, integer *, integer *, complex *, 
00061              integer *, real *);
00062     extern /* Subroutine */ int cgecon_(char *, integer *, complex *, integer 
00063             *, real *, real *, complex *, real *, integer *);
00064     extern doublereal slamch_(char *);
00065     extern /* Subroutine */ int xerbla_(char *, integer *);
00066     logical colequ, notran, rowequ;
00067     extern integer ilaprec_(char *);
00068     integer ithresh, n_norms__;
00069     real rthresh;
00070 
00071 
00072 /*     -- LAPACK routine (version 3.2.1)                                 -- */
00073 /*     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- */
00074 /*     -- Jason Riedy of Univ. of California Berkeley.                 -- */
00075 /*     -- April 2009                                                   -- */
00076 
00077 /*     -- LAPACK is a software package provided by Univ. of Tennessee, -- */
00078 /*     -- Univ. of California Berkeley and NAG Ltd.                    -- */
00079 
00080 /*     .. */
00081 /*     .. Scalar Arguments .. */
00082 /*     .. */
00083 /*     .. Array Arguments .. */
00084 /*     .. */
00085 
00086 /*     Purpose */
00087 /*     ======= */
00088 
00089 /*     CGERFSX improves the computed solution to a system of linear */
00090 /*     equations and provides error bounds and backward error estimates */
00091 /*     for the solution.  In addition to normwise error bound, the code */
00092 /*     provides maximum componentwise error bound if possible.  See */
00093 /*     comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the */
00094 /*     error bounds. */
00095 
00096 /*     The original system of linear equations may have been equilibrated */
00097 /*     before calling this routine, as described by arguments EQUED, R */
00098 /*     and C below. In this case, the solution and error bounds returned */
00099 /*     are for the original unequilibrated system. */
00100 
00101 /*     Arguments */
00102 /*     ========= */
00103 
00104 /*     Some optional parameters are bundled in the PARAMS array.  These */
00105 /*     settings determine how refinement is performed, but often the */
00106 /*     defaults are acceptable.  If the defaults are acceptable, users */
00107 /*     can pass NPARAMS = 0 which prevents the source code from accessing */
00108 /*     the PARAMS argument. */
00109 
00110 /*     TRANS   (input) CHARACTER*1 */
00111 /*     Specifies the form of the system of equations: */
00112 /*       = 'N':  A * X = B     (No transpose) */
00113 /*       = 'T':  A**T * X = B  (Transpose) */
00114 /*       = 'C':  A**H * X = B  (Conjugate transpose = Transpose) */
00115 
00116 /*     EQUED   (input) CHARACTER*1 */
00117 /*     Specifies the form of equilibration that was done to A */
00118 /*     before calling this routine. This is needed to compute */
00119 /*     the solution and error bounds correctly. */
00120 /*       = 'N':  No equilibration */
00121 /*       = 'R':  Row equilibration, i.e., A has been premultiplied by */
00122 /*               diag(R). */
00123 /*       = 'C':  Column equilibration, i.e., A has been postmultiplied */
00124 /*               by diag(C). */
00125 /*       = 'B':  Both row and column equilibration, i.e., A has been */
00126 /*               replaced by diag(R) * A * diag(C). */
00127 /*               The right hand side B has been changed accordingly. */
00128 
00129 /*     N       (input) INTEGER */
00130 /*     The order of the matrix A.  N >= 0. */
00131 
00132 /*     NRHS    (input) INTEGER */
00133 /*     The number of right hand sides, i.e., the number of columns */
00134 /*     of the matrices B and X.  NRHS >= 0. */
00135 
00136 /*     A       (input) COMPLEX array, dimension (LDA,N) */
00137 /*     The original N-by-N matrix A. */
00138 
00139 /*     LDA     (input) INTEGER */
00140 /*     The leading dimension of the array A.  LDA >= max(1,N). */
00141 
00142 /*     AF      (input) COMPLEX array, dimension (LDAF,N) */
00143 /*     The factors L and U from the factorization A = P*L*U */
00144 /*     as computed by CGETRF. */
00145 
00146 /*     LDAF    (input) INTEGER */
00147 /*     The leading dimension of the array AF.  LDAF >= max(1,N). */
00148 
00149 /*     IPIV    (input) INTEGER array, dimension (N) */
00150 /*     The pivot indices from CGETRF; for 1<=i<=N, row i of the */
00151 /*     matrix was interchanged with row IPIV(i). */
00152 
00153 /*     R       (input or output) REAL array, dimension (N) */
00154 /*     The row scale factors for A.  If EQUED = 'R' or 'B', A is */
00155 /*     multiplied on the left by diag(R); if EQUED = 'N' or 'C', R */
00156 /*     is not accessed.  R is an input argument if FACT = 'F'; */
00157 /*     otherwise, R is an output argument.  If FACT = 'F' and */
00158 /*     EQUED = 'R' or 'B', each element of R must be positive. */
00159 /*     If R is output, each element of R is a power of the radix. */
00160 /*     If R is input, each element of R should be a power of the radix */
00161 /*     to ensure a reliable solution and error estimates. Scaling by */
00162 /*     powers of the radix does not cause rounding errors unless the */
00163 /*     result underflows or overflows. Rounding errors during scaling */
00164 /*     lead to refining with a matrix that is not equivalent to the */
00165 /*     input matrix, producing error estimates that may not be */
00166 /*     reliable. */
00167 
00168 /*     C       (input or output) REAL array, dimension (N) */
00169 /*     The column scale factors for A.  If EQUED = 'C' or 'B', A is */
00170 /*     multiplied on the right by diag(C); if EQUED = 'N' or 'R', C */
00171 /*     is not accessed.  C is an input argument if FACT = 'F'; */
00172 /*     otherwise, C is an output argument.  If FACT = 'F' and */
00173 /*     EQUED = 'C' or 'B', each element of C must be positive. */
00174 /*     If C is output, each element of C is a power of the radix. */
00175 /*     If C is input, each element of C should be a power of the radix */
00176 /*     to ensure a reliable solution and error estimates. Scaling by */
00177 /*     powers of the radix does not cause rounding errors unless the */
00178 /*     result underflows or overflows. Rounding errors during scaling */
00179 /*     lead to refining with a matrix that is not equivalent to the */
00180 /*     input matrix, producing error estimates that may not be */
00181 /*     reliable. */
00182 
00183 /*     B       (input) COMPLEX array, dimension (LDB,NRHS) */
00184 /*     The right hand side matrix B. */
00185 
00186 /*     LDB     (input) INTEGER */
00187 /*     The leading dimension of the array B.  LDB >= max(1,N). */
00188 
00189 /*     X       (input/output) COMPLEX array, dimension (LDX,NRHS) */
00190 /*     On entry, the solution matrix X, as computed by CGETRS. */
00191 /*     On exit, the improved solution matrix X. */
00192 
00193 /*     LDX     (input) INTEGER */
00194 /*     The leading dimension of the array X.  LDX >= max(1,N). */
00195 
00196 /*     RCOND   (output) REAL */
00197 /*     Reciprocal scaled condition number.  This is an estimate of the */
00198 /*     reciprocal Skeel condition number of the matrix A after */
00199 /*     equilibration (if done).  If this is less than the machine */
00200 /*     precision (in particular, if it is zero), the matrix is singular */
00201 /*     to working precision.  Note that the error may still be small even */
00202 /*     if this number is very small and the matrix appears ill- */
00203 /*     conditioned. */
00204 
00205 /*     BERR    (output) REAL array, dimension (NRHS) */
00206 /*     Componentwise relative backward error.  This is the */
00207 /*     componentwise relative backward error of each solution vector X(j) */
00208 /*     (i.e., the smallest relative change in any element of A or B that */
00209 /*     makes X(j) an exact solution). */
00210 
00211 /*     N_ERR_BNDS (input) INTEGER */
00212 /*     Number of error bounds to return for each right hand side */
00213 /*     and each type (normwise or componentwise).  See ERR_BNDS_NORM and */
00214 /*     ERR_BNDS_COMP below. */
00215 
00216 /*     ERR_BNDS_NORM  (output) REAL array, dimension (NRHS, N_ERR_BNDS) */
00217 /*     For each right-hand side, this array contains information about */
00218 /*     various error bounds and condition numbers corresponding to the */
00219 /*     normwise relative error, which is defined as follows: */
00220 
00221 /*     Normwise relative error in the ith solution vector: */
00222 /*             max_j (abs(XTRUE(j,i) - X(j,i))) */
00223 /*            ------------------------------ */
00224 /*                  max_j abs(X(j,i)) */
00225 
00226 /*     The array is indexed by the type of error information as described */
00227 /*     below. There currently are up to three pieces of information */
00228 /*     returned. */
00229 
00230 /*     The first index in ERR_BNDS_NORM(i,:) corresponds to the ith */
00231 /*     right-hand side. */
00232 
00233 /*     The second index in ERR_BNDS_NORM(:,err) contains the following */
00234 /*     three fields: */
00235 /*     err = 1 "Trust/don't trust" boolean. Trust the answer if the */
00236 /*              reciprocal condition number is less than the threshold */
00237 /*              sqrt(n) * slamch('Epsilon'). */
00238 
00239 /*     err = 2 "Guaranteed" error bound: The estimated forward error, */
00240 /*              almost certainly within a factor of 10 of the true error */
00241 /*              so long as the next entry is greater than the threshold */
00242 /*              sqrt(n) * slamch('Epsilon'). This error bound should only */
00243 /*              be trusted if the previous boolean is true. */
00244 
00245 /*     err = 3  Reciprocal condition number: Estimated normwise */
00246 /*              reciprocal condition number.  Compared with the threshold */
00247 /*              sqrt(n) * slamch('Epsilon') to determine if the error */
00248 /*              estimate is "guaranteed". These reciprocal condition */
00249 /*              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some */
00250 /*              appropriately scaled matrix Z. */
00251 /*              Let Z = S*A, where S scales each row by a power of the */
00252 /*              radix so all absolute row sums of Z are approximately 1. */
00253 
00254 /*     See Lapack Working Note 165 for further details and extra */
00255 /*     cautions. */
00256 
00257 /*     ERR_BNDS_COMP  (output) REAL array, dimension (NRHS, N_ERR_BNDS) */
00258 /*     For each right-hand side, this array contains information about */
00259 /*     various error bounds and condition numbers corresponding to the */
00260 /*     componentwise relative error, which is defined as follows: */
00261 
00262 /*     Componentwise relative error in the ith solution vector: */
00263 /*                    abs(XTRUE(j,i) - X(j,i)) */
00264 /*             max_j ---------------------- */
00265 /*                         abs(X(j,i)) */
00266 
00267 /*     The array is indexed by the right-hand side i (on which the */
00268 /*     componentwise relative error depends), and the type of error */
00269 /*     information as described below. There currently are up to three */
00270 /*     pieces of information returned for each right-hand side. If */
00271 /*     componentwise accuracy is not requested (PARAMS(3) = 0.0), then */
00272 /*     ERR_BNDS_COMP is not accessed.  If N_ERR_BNDS .LT. 3, then at most */
00273 /*     the first (:,N_ERR_BNDS) entries are returned. */
00274 
00275 /*     The first index in ERR_BNDS_COMP(i,:) corresponds to the ith */
00276 /*     right-hand side. */
00277 
00278 /*     The second index in ERR_BNDS_COMP(:,err) contains the following */
00279 /*     three fields: */
00280 /*     err = 1 "Trust/don't trust" boolean. Trust the answer if the */
00281 /*              reciprocal condition number is less than the threshold */
00282 /*              sqrt(n) * slamch('Epsilon'). */
00283 
00284 /*     err = 2 "Guaranteed" error bound: The estimated forward error, */
00285 /*              almost certainly within a factor of 10 of the true error */
00286 /*              so long as the next entry is greater than the threshold */
00287 /*              sqrt(n) * slamch('Epsilon'). This error bound should only */
00288 /*              be trusted if the previous boolean is true. */
00289 
00290 /*     err = 3  Reciprocal condition number: Estimated componentwise */
00291 /*              reciprocal condition number.  Compared with the threshold */
00292 /*              sqrt(n) * slamch('Epsilon') to determine if the error */
00293 /*              estimate is "guaranteed". These reciprocal condition */
00294 /*              numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some */
00295 /*              appropriately scaled matrix Z. */
00296 /*              Let Z = S*(A*diag(x)), where x is the solution for the */
00297 /*              current right-hand side and S scales each row of */
00298 /*              A*diag(x) by a power of the radix so all absolute row */
00299 /*              sums of Z are approximately 1. */
00300 
00301 /*     See Lapack Working Note 165 for further details and extra */
00302 /*     cautions. */
00303 
00304 /*     NPARAMS (input) INTEGER */
00305 /*     Specifies the number of parameters set in PARAMS.  If .LE. 0, the */
00306 /*     PARAMS array is never referenced and default values are used. */
00307 
00308 /*     PARAMS  (input / output) REAL array, dimension NPARAMS */
00309 /*     Specifies algorithm parameters.  If an entry is .LT. 0.0, then */
00310 /*     that entry will be filled with default value used for that */
00311 /*     parameter.  Only positions up to NPARAMS are accessed; defaults */
00312 /*     are used for higher-numbered parameters. */
00313 
00314 /*       PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative */
00315 /*            refinement or not. */
00316 /*         Default: 1.0 */
00317 /*            = 0.0 : No refinement is performed, and no error bounds are */
00318 /*                    computed. */
00319 /*            = 1.0 : Use the double-precision refinement algorithm, */
00320 /*                    possibly with doubled-single computations if the */
00321 /*                    compilation environment does not support DOUBLE */
00322 /*                    PRECISION. */
00323 /*              (other values are reserved for future use) */
00324 
00325 /*       PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual */
00326 /*            computations allowed for refinement. */
00327 /*         Default: 10 */
00328 /*         Aggressive: Set to 100 to permit convergence using approximate */
00329 /*                     factorizations or factorizations other than LU. If */
00330 /*                     the factorization uses a technique other than */
00331 /*                     Gaussian elimination, the guarantees in */
00332 /*                     err_bnds_norm and err_bnds_comp may no longer be */
00333 /*                     trustworthy. */
00334 
00335 /*       PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code */
00336 /*            will attempt to find a solution with small componentwise */
00337 /*            relative error in the double-precision algorithm.  Positive */
00338 /*            is true, 0.0 is false. */
00339 /*         Default: 1.0 (attempt componentwise convergence) */
00340 
00341 /*     WORK    (workspace) COMPLEX array, dimension (2*N) */
00342 
00343 /*     RWORK   (workspace) REAL array, dimension (2*N) */
00344 
00345 /*     INFO    (output) INTEGER */
00346 /*       = 0:  Successful exit. The solution to every right-hand side is */
00347 /*         guaranteed. */
00348 /*       < 0:  If INFO = -i, the i-th argument had an illegal value */
00349 /*       > 0 and <= N:  U(INFO,INFO) is exactly zero.  The factorization */
00350 /*         has been completed, but the factor U is exactly singular, so */
00351 /*         the solution and error bounds could not be computed. RCOND = 0 */
00352 /*         is returned. */
00353 /*       = N+J: The solution corresponding to the Jth right-hand side is */
00354 /*         not guaranteed. The solutions corresponding to other right- */
00355 /*         hand sides K with K > J may not be guaranteed as well, but */
00356 /*         only the first such right-hand side is reported. If a small */
00357 /*         componentwise error is not requested (PARAMS(3) = 0.0) then */
00358 /*         the Jth right-hand side is the first with a normwise error */
00359 /*         bound that is not guaranteed (the smallest J such */
00360 /*         that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0) */
00361 /*         the Jth right-hand side is the first with either a normwise or */
00362 /*         componentwise error bound that is not guaranteed (the smallest */
00363 /*         J such that either ERR_BNDS_NORM(J,1) = 0.0 or */
00364 /*         ERR_BNDS_COMP(J,1) = 0.0). See the definition of */
00365 /*         ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information */
00366 /*         about all of the right-hand sides check ERR_BNDS_NORM or */
00367 /*         ERR_BNDS_COMP. */
00368 
00369 /*     ================================================================== */
00370 
00371 /*     .. Parameters .. */
00372 /*     .. */
00373 /*     .. Local Scalars .. */
00374 /*     .. */
00375 /*     .. External Subroutines .. */
00376 /*     .. */
00377 /*     .. Intrinsic Functions .. */
00378 /*     .. */
00379 /*     .. External Functions .. */
00380 /*     .. */
00381 /*     .. Executable Statements .. */
00382 
00383 /*     Check the input parameters. */
00384 
00385     /* Parameter adjustments */
00386     err_bnds_comp_dim1 = *nrhs;
00387     err_bnds_comp_offset = 1 + err_bnds_comp_dim1;
00388     err_bnds_comp__ -= err_bnds_comp_offset;
00389     err_bnds_norm_dim1 = *nrhs;
00390     err_bnds_norm_offset = 1 + err_bnds_norm_dim1;
00391     err_bnds_norm__ -= err_bnds_norm_offset;
00392     a_dim1 = *lda;
00393     a_offset = 1 + a_dim1;
00394     a -= a_offset;
00395     af_dim1 = *ldaf;
00396     af_offset = 1 + af_dim1;
00397     af -= af_offset;
00398     --ipiv;
00399     --r__;
00400     --c__;
00401     b_dim1 = *ldb;
00402     b_offset = 1 + b_dim1;
00403     b -= b_offset;
00404     x_dim1 = *ldx;
00405     x_offset = 1 + x_dim1;
00406     x -= x_offset;
00407     --berr;
00408     --params;
00409     --work;
00410     --rwork;
00411 
00412     /* Function Body */
00413     *info = 0;
00414     trans_type__ = ilatrans_(trans);
00415     ref_type__ = 1;
00416     if (*nparams >= 1) {
00417         if (params[1] < 0.f) {
00418             params[1] = 1.f;
00419         } else {
00420             ref_type__ = params[1];
00421         }
00422     }
00423 
00424 /*     Set default parameters. */
00425 
00426     illrcond_thresh__ = (real) (*n) * slamch_("Epsilon");
00427     ithresh = 10;
00428     rthresh = .5f;
00429     unstable_thresh__ = .25f;
00430     ignore_cwise__ = FALSE_;
00431 
00432     if (*nparams >= 2) {
00433         if (params[2] < 0.f) {
00434             params[2] = (real) ithresh;
00435         } else {
00436             ithresh = (integer) params[2];
00437         }
00438     }
00439     if (*nparams >= 3) {
00440         if (params[3] < 0.f) {
00441             if (ignore_cwise__) {
00442                 params[3] = 0.f;
00443             } else {
00444                 params[3] = 1.f;
00445             }
00446         } else {
00447             ignore_cwise__ = params[3] == 0.f;
00448         }
00449     }
00450     if (ref_type__ == 0 || *n_err_bnds__ == 0) {
00451         n_norms__ = 0;
00452     } else if (ignore_cwise__) {
00453         n_norms__ = 1;
00454     } else {
00455         n_norms__ = 2;
00456     }
00457 
00458     notran = lsame_(trans, "N");
00459     rowequ = lsame_(equed, "R") || lsame_(equed, "B");
00460     colequ = lsame_(equed, "C") || lsame_(equed, "B");
00461 
00462 /*     Test input parameters. */
00463 
00464     if (trans_type__ == -1) {
00465         *info = -1;
00466     } else if (! rowequ && ! colequ && ! lsame_(equed, "N")) {
00467         *info = -2;
00468     } else if (*n < 0) {
00469         *info = -3;
00470     } else if (*nrhs < 0) {
00471         *info = -4;
00472     } else if (*lda < max(1,*n)) {
00473         *info = -6;
00474     } else if (*ldaf < max(1,*n)) {
00475         *info = -8;
00476     } else if (*ldb < max(1,*n)) {
00477         *info = -13;
00478     } else if (*ldx < max(1,*n)) {
00479         *info = -15;
00480     }
00481     if (*info != 0) {
00482         i__1 = -(*info);
00483         xerbla_("CGERFSX", &i__1);
00484         return 0;
00485     }
00486 
00487 /*     Quick return if possible. */
00488 
00489     if (*n == 0 || *nrhs == 0) {
00490         *rcond = 1.f;
00491         i__1 = *nrhs;
00492         for (j = 1; j <= i__1; ++j) {
00493             berr[j] = 0.f;
00494             if (*n_err_bnds__ >= 1) {
00495                 err_bnds_norm__[j + err_bnds_norm_dim1] = 1.f;
00496                 err_bnds_comp__[j + err_bnds_comp_dim1] = 1.f;
00497             } else if (*n_err_bnds__ >= 2) {
00498                 err_bnds_norm__[j + (err_bnds_norm_dim1 << 1)] = 0.f;
00499                 err_bnds_comp__[j + (err_bnds_comp_dim1 << 1)] = 0.f;
00500             } else if (*n_err_bnds__ >= 3) {
00501                 err_bnds_norm__[j + err_bnds_norm_dim1 * 3] = 1.f;
00502                 err_bnds_comp__[j + err_bnds_comp_dim1 * 3] = 1.f;
00503             }
00504         }
00505         return 0;
00506     }
00507 
00508 /*     Default to failure. */
00509 
00510     *rcond = 0.f;
00511     i__1 = *nrhs;
00512     for (j = 1; j <= i__1; ++j) {
00513         berr[j] = 1.f;
00514         if (*n_err_bnds__ >= 1) {
00515             err_bnds_norm__[j + err_bnds_norm_dim1] = 1.f;
00516             err_bnds_comp__[j + err_bnds_comp_dim1] = 1.f;
00517         } else if (*n_err_bnds__ >= 2) {
00518             err_bnds_norm__[j + (err_bnds_norm_dim1 << 1)] = 1.f;
00519             err_bnds_comp__[j + (err_bnds_comp_dim1 << 1)] = 1.f;
00520         } else if (*n_err_bnds__ >= 3) {
00521             err_bnds_norm__[j + err_bnds_norm_dim1 * 3] = 0.f;
00522             err_bnds_comp__[j + err_bnds_comp_dim1 * 3] = 0.f;
00523         }
00524     }
00525 
00526 /*     Compute the norm of A and the reciprocal of the condition */
00527 /*     number of A. */
00528 
00529     if (notran) {
00530         *(unsigned char *)norm = 'I';
00531     } else {
00532         *(unsigned char *)norm = '1';
00533     }
00534     anorm = clange_(norm, n, n, &a[a_offset], lda, &rwork[1]);
00535     cgecon_(norm, n, &af[af_offset], ldaf, &anorm, rcond, &work[1], &rwork[1], 
00536              info);
00537 
00538 /*     Perform refinement on each right-hand side */
00539 
00540     if (ref_type__ != 0) {
00541         prec_type__ = ilaprec_("D");
00542         if (notran) {
00543             cla_gerfsx_extended__(&prec_type__, &trans_type__, n, nrhs, &a[
00544                     a_offset], lda, &af[af_offset], ldaf, &ipiv[1], &colequ, &
00545                     c__[1], &b[b_offset], ldb, &x[x_offset], ldx, &berr[1], &
00546                     n_norms__, &err_bnds_norm__[err_bnds_norm_offset], &
00547                     err_bnds_comp__[err_bnds_comp_offset], &work[1], &rwork[1]
00548                     , &work[*n + 1], (complex *)(&rwork[1]), rcond, &ithresh, &rthresh, &
00549                     unstable_thresh__, &ignore_cwise__, info);
00550         } else {
00551             cla_gerfsx_extended__(&prec_type__, &trans_type__, n, nrhs, &a[
00552                     a_offset], lda, &af[af_offset], ldaf, &ipiv[1], &rowequ, &
00553                     r__[1], &b[b_offset], ldb, &x[x_offset], ldx, &berr[1], &
00554                     n_norms__, &err_bnds_norm__[err_bnds_norm_offset], &
00555                     err_bnds_comp__[err_bnds_comp_offset], &work[1], &rwork[1]
00556                     , &work[*n + 1], (complex *)(&rwork[1]), rcond, &ithresh, &rthresh, &
00557                     unstable_thresh__, &ignore_cwise__, info);
00558         }
00559     }
00560 /* Computing MAX */
00561     r__1 = 10.f, r__2 = sqrt((real) (*n));
00562     err_lbnd__ = dmax(r__1,r__2) * slamch_("Epsilon");
00563     if (*n_err_bnds__ >= 1 && n_norms__ >= 1) {
00564 
00565 /*     Compute scaled normwise condition number cond(A*C). */
00566 
00567         if (colequ && notran) {
00568             rcond_tmp__ = cla_gercond_c__(trans, n, &a[a_offset], lda, &af[
00569                     af_offset], ldaf, &ipiv[1], &c__[1], &c_true, info, &work[
00570                     1], &rwork[1], (ftnlen)1);
00571         } else if (rowequ && ! notran) {
00572             rcond_tmp__ = cla_gercond_c__(trans, n, &a[a_offset], lda, &af[
00573                     af_offset], ldaf, &ipiv[1], &r__[1], &c_true, info, &work[
00574                     1], &rwork[1], (ftnlen)1);
00575         } else {
00576             rcond_tmp__ = cla_gercond_c__(trans, n, &a[a_offset], lda, &af[
00577                     af_offset], ldaf, &ipiv[1], &c__[1], &c_false, info, &
00578                     work[1], &rwork[1], (ftnlen)1);
00579         }
00580         i__1 = *nrhs;
00581         for (j = 1; j <= i__1; ++j) {
00582 
00583 /*     Cap the error at 1.0. */
00584 
00585             if (*n_err_bnds__ >= 2 && err_bnds_norm__[j + (err_bnds_norm_dim1 
00586                     << 1)] > 1.f) {
00587                 err_bnds_norm__[j + (err_bnds_norm_dim1 << 1)] = 1.f;
00588             }
00589 
00590 /*     Threshold the error (see LAWN). */
00591 
00592             if (rcond_tmp__ < illrcond_thresh__) {
00593                 err_bnds_norm__[j + (err_bnds_norm_dim1 << 1)] = 1.f;
00594                 err_bnds_norm__[j + err_bnds_norm_dim1] = 0.f;
00595                 if (*info <= *n) {
00596                     *info = *n + j;
00597                 }
00598             } else if (err_bnds_norm__[j + (err_bnds_norm_dim1 << 1)] < 
00599                     err_lbnd__) {
00600                 err_bnds_norm__[j + (err_bnds_norm_dim1 << 1)] = err_lbnd__;
00601                 err_bnds_norm__[j + err_bnds_norm_dim1] = 1.f;
00602             }
00603 
00604 /*     Save the condition number. */
00605 
00606             if (*n_err_bnds__ >= 3) {
00607                 err_bnds_norm__[j + err_bnds_norm_dim1 * 3] = rcond_tmp__;
00608             }
00609         }
00610     }
00611     if (*n_err_bnds__ >= 1 && n_norms__ >= 2) {
00612 
00613 /*     Compute componentwise condition number cond(A*diag(Y(:,J))) for */
00614 /*     each right-hand side using the current solution as an estimate of */
00615 /*     the true solution.  If the componentwise error estimate is too */
00616 /*     large, then the solution is a lousy estimate of truth and the */
00617 /*     estimated RCOND may be too optimistic.  To avoid misleading users, */
00618 /*     the inverse condition number is set to 0.0 when the estimated */
00619 /*     cwise error is at least CWISE_WRONG. */
00620 
00621         cwise_wrong__ = sqrt(slamch_("Epsilon"));
00622         i__1 = *nrhs;
00623         for (j = 1; j <= i__1; ++j) {
00624             if (err_bnds_comp__[j + (err_bnds_comp_dim1 << 1)] < 
00625                     cwise_wrong__) {
00626                 rcond_tmp__ = cla_gercond_x__(trans, n, &a[a_offset], lda, &
00627                         af[af_offset], ldaf, &ipiv[1], &x[j * x_dim1 + 1], 
00628                         info, &work[1], &rwork[1], (ftnlen)1);
00629             } else {
00630                 rcond_tmp__ = 0.f;
00631             }
00632 
00633 /*     Cap the error at 1.0. */
00634 
00635             if (*n_err_bnds__ >= 2 && err_bnds_comp__[j + (err_bnds_comp_dim1 
00636                     << 1)] > 1.f) {
00637                 err_bnds_comp__[j + (err_bnds_comp_dim1 << 1)] = 1.f;
00638             }
00639 
00640 /*     Threshold the error (see LAWN). */
00641 
00642             if (rcond_tmp__ < illrcond_thresh__) {
00643                 err_bnds_comp__[j + (err_bnds_comp_dim1 << 1)] = 1.f;
00644                 err_bnds_comp__[j + err_bnds_comp_dim1] = 0.f;
00645                 if (params[3] == 1.f && *info < *n + j) {
00646                     *info = *n + j;
00647                 }
00648             } else if (err_bnds_comp__[j + (err_bnds_comp_dim1 << 1)] < 
00649                     err_lbnd__) {
00650                 err_bnds_comp__[j + (err_bnds_comp_dim1 << 1)] = err_lbnd__;
00651                 err_bnds_comp__[j + err_bnds_comp_dim1] = 1.f;
00652             }
00653 
00654 /*     Save the condition number. */
00655 
00656             if (*n_err_bnds__ >= 3) {
00657                 err_bnds_comp__[j + err_bnds_comp_dim1 * 3] = rcond_tmp__;
00658             }
00659         }
00660     }
00661 
00662     return 0;
00663 
00664 /*     End of CGERFSX */
00665 
00666 } /* cgerfsx_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:26