sebchvxx.c
Go to the documentation of this file.
00001 /* sebchvxx.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 integer c__6 = 6;
00019 static integer c__2 = 2;
00020 static integer c__3 = 3;
00021 static integer c__1 = 1;
00022 static integer c__4 = 4;
00023 static integer c__5 = 5;
00024 static integer c__7 = 7;
00025 static integer c__8 = 8;
00026 
00027 /* Subroutine */ int sebchvxx_(real *thresh, char *path)
00028 {
00029     /* Format strings */
00030     static char fmt_8000[] = "(\002 S\002,a2,\002SVXX: N =\002,i2,\002, INFO"
00031             " = \002,i3,\002, ORCOND = \002,g12.5,\002, real RCOND = \002,g12"
00032             ".5)";
00033     static char fmt_9996[] = "(3x,i2,\002: Normwise guaranteed forward erro"
00034             "r\002,/5x,\002Guaranteed case: if norm ( abs( Xc - Xt )\002,\002"
00035             " / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ), then\002,/5x"
00036             ",\002ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS\002)"
00037             ;
00038     static char fmt_9995[] = "(3x,i2,\002: Componentwise guaranteed forward "
00039             "error\002)";
00040     static char fmt_9994[] = "(3x,i2,\002: Backwards error\002)";
00041     static char fmt_9993[] = "(3x,i2,\002: Reciprocal condition number\002)";
00042     static char fmt_9992[] = "(3x,i2,\002: Reciprocal normwise condition num"
00043             "ber\002)";
00044     static char fmt_9991[] = "(3x,i2,\002: Raw normwise error estimate\002)";
00045     static char fmt_9990[] = "(3x,i2,\002: Reciprocal componentwise conditio"
00046             "n number\002)";
00047     static char fmt_9989[] = "(3x,i2,\002: Raw componentwise error estimat"
00048             "e\002)";
00049     static char fmt_9999[] = "(\002 S\002,a2,\002SVXX: N =\002,i2,\002, RHS "
00050             "= \002,i2,\002, NWISE GUAR. = \002,a,\002, CWISE GUAR. = \002,a"
00051             ",\002 test(\002,i1,\002) =\002,g12.5)";
00052     static char fmt_9998[] = "(\002 S\002,a2,\002SVXX: \002,i6,\002 out of"
00053             " \002,i6,\002 tests failed to pass the threshold\002)";
00054     static char fmt_9997[] = "(\002 S\002,a2,\002SVXX passed the tests of er"
00055             "ror bounds\002)";
00056 
00057     /* System generated locals */
00058     integer i__1, i__2, i__3, i__4, i__5, i__6;
00059     real r__1, r__2, r__3;
00060 
00061     /* Builtin functions */
00062     /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
00063     double sqrt(doublereal);
00064     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void),
00065              s_wsle(cilist *), e_wsle(void);
00066 
00067     /* Local variables */
00068     extern /* Subroutine */ int sposvxx_(char *, char *, integer *, integer *, 
00069              real *, integer *, real *, integer *, char *, real *, real *, 
00070             integer *, real *, integer *, real *, real *, real *, integer *, 
00071             real *, real *, integer *, real *, real *, integer *, integer *), ssysvxx_(char *, char *, integer *, 
00072             integer *, real *, integer *, real *, integer *, integer *, char *
00073 , real *, real *, integer *, real *, integer *, real *, real *, 
00074             real *, integer *, real *, real *, integer *, real *, real *, 
00075             integer *, integer *);
00076     real errbnd_c__[18], errbnd_n__[18], a[36]  /* was [6][6] */, b[36] /* 
00077             was [6][6] */, c__[6];
00078     integer i__, j, k;
00079     real m;
00080     integer n;
00081     real r__[6], s[6], x[36]    /* was [6][6] */, cwise_bnd__;
00082     char c2[2];
00083     real nwise_bnd__, cwise_err__, nwise_err__, errthresh, ab[66]       /* 
00084             was [11][6] */, af[36]      /* was [6][6] */;
00085     integer kl, ku;
00086     real condthresh, afb[96]    /* was [16][6] */;
00087     integer lda;
00088     real eps, cwise_rcond__, nwise_rcond__;
00089     integer n_aux_tests__, ldab;
00090     real diff[36]       /* was [6][6] */;
00091     char fact[1];
00092     real berr[6];
00093     integer info, ipiv[6], nrhs;
00094     real rinv[6];
00095     char uplo[1];
00096     real work[30], sumr;
00097     integer ldafb;
00098     real ccond;
00099     integer nfail;
00100     char cguar[3];
00101     real ncond;
00102     char equed[1];
00103     real rcond, acopy[36]       /* was [6][6] */;
00104     char nguar[3], trans[1];
00105     integer iwork[6];
00106     real rnorm, normt, sumri;
00107     logical printed_guide__;
00108     extern doublereal slamch_(char *);
00109     real abcopy[66]     /* was [11][6] */;
00110     extern logical lsamen_(integer *, char *, char *);
00111     real params[2], orcond;
00112     extern /* Subroutine */ int slacpy_(char *, integer *, integer *, real *, 
00113             integer *, real *, integer *);
00114     real rinorm, tstrat[6], rpvgrw;
00115     extern /* Subroutine */ int slahilb_(integer *, integer *, real *, 
00116             integer *, real *, integer *, real *, integer *, real *, integer *
00117 );
00118     real invhilb[36]    /* was [6][6] */, normdif;
00119     extern /* Subroutine */ int sgbsvxx_(char *, char *, integer *, integer *, 
00120              integer *, integer *, real *, integer *, real *, integer *, 
00121             integer *, char *, real *, real *, real *, integer *, real *, 
00122             integer *, real *, real *, real *, integer *, real *, real *, 
00123             integer *, real *, real *, integer *, integer *), sgesvxx_(char *, char *, integer *, integer *, real *, 
00124             integer *, real *, integer *, integer *, char *, real *, real *, 
00125             real *, integer *, real *, integer *, real *, real *, real *, 
00126             integer *, real *, real *, integer *, real *, real *, integer *, 
00127             integer *);
00128 
00129     /* Fortran I/O blocks */
00130     static cilist io___42 = { 0, 6, 0, fmt_8000, 0 };
00131     static cilist io___66 = { 0, 6, 0, 0, 0 };
00132     static cilist io___67 = { 0, 6, 0, fmt_9996, 0 };
00133     static cilist io___68 = { 0, 6, 0, fmt_9995, 0 };
00134     static cilist io___69 = { 0, 6, 0, fmt_9994, 0 };
00135     static cilist io___70 = { 0, 6, 0, fmt_9993, 0 };
00136     static cilist io___71 = { 0, 6, 0, fmt_9992, 0 };
00137     static cilist io___72 = { 0, 6, 0, fmt_9991, 0 };
00138     static cilist io___73 = { 0, 6, 0, fmt_9990, 0 };
00139     static cilist io___74 = { 0, 6, 0, fmt_9989, 0 };
00140     static cilist io___75 = { 0, 6, 0, 0, 0 };
00141     static cilist io___76 = { 0, 6, 0, fmt_9999, 0 };
00142     static cilist io___77 = { 0, 6, 0, 0, 0 };
00143     static cilist io___78 = { 0, 6, 0, fmt_9998, 0 };
00144     static cilist io___79 = { 0, 6, 0, fmt_9997, 0 };
00145 
00146 
00147 /*     .. Scalar Arguments .. */
00148 
00149 /*  Purpose */
00150 /*  ====== */
00151 
00152 /*  SEBCHVXX will run S**SVXX on a series of Hilbert matrices and then */
00153 /*  compare the error bounds returned by SGESVXX to see if the returned */
00154 /*  answer indeed falls within those bounds. */
00155 
00156 /*  Eight test ratios will be computed.  The tests will pass if they are .LT. */
00157 /*  THRESH.  There are two cases that are determined by 1 / (SQRT( N ) * EPS). */
00158 /*  If that value is .LE. to the component wise reciprocal condition number, */
00159 /*  it uses the guaranteed case, other wise it uses the unguaranteed case. */
00160 
00161 /*  Test ratios: */
00162 /*     Let Xc be X_computed and Xt be X_truth. */
00163 /*     The norm used is the infinity norm. */
00164 /*     Let A be the guaranteed case and B be the unguaranteed case. */
00165 
00166 /*       1. Normwise guaranteed forward error bound. */
00167 /*       A: norm ( abs( Xc - Xt ) / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ) and */
00168 /*          ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N),10) * EPS. */
00169 /*          If these conditions are met, the test ratio is set to be */
00170 /*          ERRBND( *, nwise_i, bnd_i ) / MAX(SQRT(N), 10).  Otherwise it is 1/EPS. */
00171 /*       B: For this case, SGESVXX should just return 1.  If it is less than */
00172 /*          one, treat it the same as in 1A.  Otherwise it fails. (Set test */
00173 /*          ratio to ERRBND( *, nwise_i, bnd_i ) * THRESH?) */
00174 
00175 /*       2. Componentwise guaranteed forward error bound. */
00176 /*       A: norm ( abs( Xc(j) - Xt(j) ) ) / norm (Xt(j)) .LE. ERRBND( *, cwise_i, bnd_i ) */
00177 /*          for all j .AND. ERRBND( *, cwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS. */
00178 /*          If these conditions are met, the test ratio is set to be */
00179 /*          ERRBND( *, cwise_i, bnd_i ) / MAX(SQRT(N), 10).  Otherwise it is 1/EPS. */
00180 /*       B: Same as normwise test ratio. */
00181 
00182 /*       3. Backwards error. */
00183 /*       A: The test ratio is set to BERR/EPS. */
00184 /*       B: Same test ratio. */
00185 
00186 /*       4. Reciprocal condition number. */
00187 /*       A: A condition number is computed with Xt and compared with the one */
00188 /*          returned from SGESVXX.  Let RCONDc be the RCOND returned by SGESVXX */
00189 /*          and RCONDt be the RCOND from the truth value.  Test ratio is set to */
00190 /*          MAX(RCONDc/RCONDt, RCONDt/RCONDc). */
00191 /*       B: Test ratio is set to 1 / (EPS * RCONDc). */
00192 
00193 /*       5. Reciprocal normwise condition number. */
00194 /*       A: The test ratio is set to */
00195 /*          MAX(ERRBND( *, nwise_i, cond_i ) / NCOND, NCOND / ERRBND( *, nwise_i, cond_i )). */
00196 /*       B: Test ratio is set to 1 / (EPS * ERRBND( *, nwise_i, cond_i )). */
00197 
00198 /*       7. Reciprocal componentwise condition number. */
00199 /*       A: Test ratio is set to */
00200 /*          MAX(ERRBND( *, cwise_i, cond_i ) / CCOND, CCOND / ERRBND( *, cwise_i, cond_i )). */
00201 /*       B: Test ratio is set to 1 / (EPS * ERRBND( *, cwise_i, cond_i )). */
00202 
00203 /*     .. Parameters .. */
00204 /*     NMAX is determined by the largest number in the inverse of the Hilbert */
00205 /*     matrix.  Precision is exhausted when the largest entry in it is greater */
00206 /*     than 2 to the power of the number of bits in the fraction of the data */
00207 /*     type used plus one, which is 24 for single precision. */
00208 /*     NMAX should be 6 for single and 11 for double. */
00209 /*     .. Local Scalars .. */
00210 /*     .. Local Arrays .. */
00211 /*     .. External Functions .. */
00212 /*     .. External Subroutines .. */
00213 /*     .. Intrinsic Functions .. */
00214 /*     .. Parameters .. */
00215 /*     Create the loop to test out the Hilbert matrices */
00216     *(unsigned char *)fact = 'E';
00217     *(unsigned char *)uplo = 'U';
00218     *(unsigned char *)trans = 'N';
00219     *(unsigned char *)equed = 'N';
00220     eps = slamch_("Epsilon");
00221     nfail = 0;
00222     n_aux_tests__ = 0;
00223     lda = 6;
00224     ldab = 11;
00225     ldafb = 16;
00226     s_copy(c2, path + 1, (ftnlen)2, (ftnlen)2);
00227 /*     Main loop to test the different Hilbert Matrices. */
00228     printed_guide__ = FALSE_;
00229     for (n = 1; n <= 6; ++n) {
00230         params[0] = -1.f;
00231         params[1] = -1.f;
00232         kl = n - 1;
00233         ku = n - 1;
00234         nrhs = n;
00235 /* Computing MAX */
00236         r__1 = sqrt((real) n);
00237         m = dmax(r__1,10.f);
00238 /*        Generate the Hilbert matrix, its inverse, and the */
00239 /*        right hand side, all scaled by the LCM(1,..,2N-1). */
00240         slahilb_(&n, &n, a, &lda, invhilb, &lda, b, &lda, work, &info);
00241 /*        Copy A into ACOPY. */
00242         slacpy_("ALL", &n, &n, a, &c__6, acopy, &c__6);
00243 /*        Store A in band format for GB tests */
00244         i__1 = n;
00245         for (j = 1; j <= i__1; ++j) {
00246             i__2 = kl + ku + 1;
00247             for (i__ = 1; i__ <= i__2; ++i__) {
00248                 ab[i__ + j * 11 - 12] = 0.f;
00249             }
00250         }
00251         i__1 = n;
00252         for (j = 1; j <= i__1; ++j) {
00253 /* Computing MAX */
00254             i__2 = 1, i__3 = j - ku;
00255 /* Computing MIN */
00256             i__5 = n, i__6 = j + kl;
00257             i__4 = min(i__5,i__6);
00258             for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
00259                 ab[ku + 1 + i__ - j + j * 11 - 12] = a[i__ + j * 6 - 7];
00260             }
00261         }
00262 /*        Copy AB into ABCOPY. */
00263         i__1 = n;
00264         for (j = 1; j <= i__1; ++j) {
00265             i__4 = kl + ku + 1;
00266             for (i__ = 1; i__ <= i__4; ++i__) {
00267                 abcopy[i__ + j * 11 - 12] = 0.f;
00268             }
00269         }
00270         i__1 = kl + ku + 1;
00271         slacpy_("ALL", &i__1, &n, ab, &ldab, abcopy, &ldab);
00272 /*        Call S**SVXX with default PARAMS and N_ERR_BND = 3. */
00273         if (lsamen_(&c__2, c2, "SY")) {
00274             ssysvxx_(fact, uplo, &n, &nrhs, acopy, &lda, af, &lda, ipiv, 
00275                     equed, s, b, &lda, x, &lda, &orcond, &rpvgrw, berr, &c__3, 
00276                      errbnd_n__, errbnd_c__, &c__2, params, work, iwork, &
00277                     info);
00278         } else if (lsamen_(&c__2, c2, "PO")) {
00279             sposvxx_(fact, uplo, &n, &nrhs, acopy, &lda, af, &lda, equed, s, 
00280                     b, &lda, x, &lda, &orcond, &rpvgrw, berr, &c__3, 
00281                     errbnd_n__, errbnd_c__, &c__2, params, work, iwork, &info);
00282         } else if (lsamen_(&c__2, c2, "GB")) {
00283             sgbsvxx_(fact, trans, &n, &kl, &ku, &nrhs, abcopy, &ldab, afb, &
00284                     ldafb, ipiv, equed, r__, c__, b, &lda, x, &lda, &orcond, &
00285                     rpvgrw, berr, &c__3, errbnd_n__, errbnd_c__, &c__2, 
00286                     params, work, iwork, &info);
00287         } else {
00288             sgesvxx_(fact, trans, &n, &nrhs, acopy, &lda, af, &lda, ipiv, 
00289                     equed, r__, c__, b, &lda, x, &lda, &orcond, &rpvgrw, berr, 
00290                      &c__3, errbnd_n__, errbnd_c__, &c__2, params, work, 
00291                     iwork, &info);
00292         }
00293         ++n_aux_tests__;
00294         if (orcond < eps) {
00295 /*        Either factorization failed or the matrix is flagged, and 1 <= */
00296 /*        INFO <= N+1. We don't decide based on rcond anymore. */
00297 /*            IF (INFO .EQ. 0 .OR. INFO .GT. N+1) THEN */
00298 /*               NFAIL = NFAIL + 1 */
00299 /*               WRITE (*, FMT=8000) N, INFO, ORCOND, RCOND */
00300 /*            END IF */
00301         } else {
00302 /*        Either everything succeeded (INFO == 0) or some solution failed */
00303 /*        to converge (INFO > N+1). */
00304             if (info > 0 && info <= n + 1) {
00305                 ++nfail;
00306                 s_wsfe(&io___42);
00307                 do_fio(&c__1, c2, (ftnlen)2);
00308                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00309                 do_fio(&c__1, (char *)&info, (ftnlen)sizeof(integer));
00310                 do_fio(&c__1, (char *)&orcond, (ftnlen)sizeof(real));
00311                 do_fio(&c__1, (char *)&rcond, (ftnlen)sizeof(real));
00312                 e_wsfe();
00313             }
00314         }
00315 /*        Calculating the difference between S**SVXX's X and the true X. */
00316         i__1 = n;
00317         for (i__ = 1; i__ <= i__1; ++i__) {
00318             i__4 = nrhs;
00319             for (j = 1; j <= i__4; ++j) {
00320                 diff[i__ + j * 6 - 7] = x[i__ + j * 6 - 7] - invhilb[i__ + j *
00321                          6 - 7];
00322             }
00323         }
00324 /*        Calculating the RCOND */
00325         rnorm = 0.f;
00326         rinorm = 0.f;
00327         if (lsamen_(&c__2, c2, "PO") || lsamen_(&c__2, 
00328                 c2, "SY")) {
00329             i__1 = n;
00330             for (i__ = 1; i__ <= i__1; ++i__) {
00331                 sumr = 0.f;
00332                 sumri = 0.f;
00333                 i__4 = n;
00334                 for (j = 1; j <= i__4; ++j) {
00335                     sumr += (r__1 = s[i__ - 1] * a[i__ + j * 6 - 7] * s[j - 1]
00336                             , dabs(r__1));
00337                     sumri += (r__1 = invhilb[i__ + j * 6 - 7] / s[j - 1] / s[
00338                             i__ - 1], dabs(r__1));
00339                 }
00340                 rnorm = dmax(rnorm,sumr);
00341                 rinorm = dmax(rinorm,sumri);
00342             }
00343         } else if (lsamen_(&c__2, c2, "GE") || lsamen_(&
00344                 c__2, c2, "GB")) {
00345             i__1 = n;
00346             for (i__ = 1; i__ <= i__1; ++i__) {
00347                 sumr = 0.f;
00348                 sumri = 0.f;
00349                 i__4 = n;
00350                 for (j = 1; j <= i__4; ++j) {
00351                     sumr += (r__1 = r__[i__ - 1] * a[i__ + j * 6 - 7] * c__[j 
00352                             - 1], dabs(r__1));
00353                     sumri += (r__1 = invhilb[i__ + j * 6 - 7] / r__[j - 1] / 
00354                             c__[i__ - 1], dabs(r__1));
00355                 }
00356                 rnorm = dmax(rnorm,sumr);
00357                 rinorm = dmax(rinorm,sumri);
00358             }
00359         }
00360         rnorm /= a[0];
00361         rcond = 1.f / (rnorm * rinorm);
00362 /*        Calculating the R for normwise rcond. */
00363         i__1 = n;
00364         for (i__ = 1; i__ <= i__1; ++i__) {
00365             rinv[i__ - 1] = 0.f;
00366         }
00367         i__1 = n;
00368         for (j = 1; j <= i__1; ++j) {
00369             i__4 = n;
00370             for (i__ = 1; i__ <= i__4; ++i__) {
00371                 rinv[i__ - 1] += (r__1 = a[i__ + j * 6 - 7], dabs(r__1));
00372             }
00373         }
00374 /*        Calculating the Normwise rcond. */
00375         rinorm = 0.f;
00376         i__1 = n;
00377         for (i__ = 1; i__ <= i__1; ++i__) {
00378             sumri = 0.f;
00379             i__4 = n;
00380             for (j = 1; j <= i__4; ++j) {
00381                 sumri += (r__1 = invhilb[i__ + j * 6 - 7] * rinv[j - 1], dabs(
00382                         r__1));
00383             }
00384             rinorm = dmax(rinorm,sumri);
00385         }
00386 /*        invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm */
00387 /*        by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix) */
00388         ncond = a[0] / rinorm;
00389         condthresh = m * eps;
00390         errthresh = m * eps;
00391         i__1 = nrhs;
00392         for (k = 1; k <= i__1; ++k) {
00393             normt = 0.f;
00394             normdif = 0.f;
00395             cwise_err__ = 0.f;
00396             i__4 = n;
00397             for (i__ = 1; i__ <= i__4; ++i__) {
00398 /* Computing MAX */
00399                 r__2 = (r__1 = invhilb[i__ + k * 6 - 7], dabs(r__1));
00400                 normt = dmax(r__2,normt);
00401 /* Computing MAX */
00402                 r__2 = (r__1 = x[i__ + k * 6 - 7] - invhilb[i__ + k * 6 - 7], 
00403                         dabs(r__1));
00404                 normdif = dmax(r__2,normdif);
00405                 if (invhilb[i__ + k * 6 - 7] != 0.f) {
00406 /* Computing MAX */
00407                     r__3 = (r__1 = x[i__ + k * 6 - 7] - invhilb[i__ + k * 6 - 
00408                             7], dabs(r__1)) / (r__2 = invhilb[i__ + k * 6 - 7]
00409                             , dabs(r__2));
00410                     cwise_err__ = dmax(r__3,cwise_err__);
00411                 } else if (x[i__ + k * 6 - 7] != 0.f) {
00412                     cwise_err__ = slamch_("OVERFLOW");
00413                 }
00414             }
00415             if (normt != 0.f) {
00416                 nwise_err__ = normdif / normt;
00417             } else if (normdif != 0.f) {
00418                 nwise_err__ = slamch_("OVERFLOW");
00419             } else {
00420                 nwise_err__ = 0.f;
00421             }
00422             i__4 = n;
00423             for (i__ = 1; i__ <= i__4; ++i__) {
00424                 rinv[i__ - 1] = 0.f;
00425             }
00426             i__4 = n;
00427             for (j = 1; j <= i__4; ++j) {
00428                 i__2 = n;
00429                 for (i__ = 1; i__ <= i__2; ++i__) {
00430                     rinv[i__ - 1] += (r__1 = a[i__ + j * 6 - 7] * invhilb[j + 
00431                             k * 6 - 7], dabs(r__1));
00432                 }
00433             }
00434             rinorm = 0.f;
00435             i__4 = n;
00436             for (i__ = 1; i__ <= i__4; ++i__) {
00437                 sumri = 0.f;
00438                 i__2 = n;
00439                 for (j = 1; j <= i__2; ++j) {
00440                     sumri += (r__1 = invhilb[i__ + j * 6 - 7] * rinv[j - 1] / 
00441                             invhilb[i__ + k * 6 - 7], dabs(r__1));
00442                 }
00443                 rinorm = dmax(rinorm,sumri);
00444             }
00445 /*        invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm */
00446 /*        by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix) */
00447             ccond = a[0] / rinorm;
00448 /*        Forward error bound tests */
00449             nwise_bnd__ = errbnd_n__[k + nrhs - 1];
00450             cwise_bnd__ = errbnd_c__[k + nrhs - 1];
00451             nwise_rcond__ = errbnd_n__[k + (nrhs << 1) - 1];
00452             cwise_rcond__ = errbnd_c__[k + (nrhs << 1) - 1];
00453 /*            write (*,*) 'nwise : ', n, k, ncond, nwise_rcond, */
00454 /*     $           condthresh, ncond.ge.condthresh */
00455 /*            write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh */
00456             if (ncond >= condthresh) {
00457                 s_copy(nguar, "YES", (ftnlen)3, (ftnlen)3);
00458                 if (nwise_bnd__ > errthresh) {
00459                     tstrat[0] = 1 / (eps * 2.f);
00460                 } else {
00461                     if (nwise_bnd__ != 0.f) {
00462                         tstrat[0] = nwise_err__ / nwise_bnd__;
00463                     } else if (nwise_err__ != 0.f) {
00464                         tstrat[0] = 1 / (eps * 16.f);
00465                     } else {
00466                         tstrat[0] = 0.f;
00467                     }
00468                     if (tstrat[0] > 1.f) {
00469                         tstrat[0] = 1 / (eps * 4.f);
00470                     }
00471                 }
00472             } else {
00473                 s_copy(nguar, "NO", (ftnlen)3, (ftnlen)2);
00474                 if (nwise_bnd__ < 1.f) {
00475                     tstrat[0] = 1 / (eps * 8.f);
00476                 } else {
00477                     tstrat[0] = 1.f;
00478                 }
00479             }
00480 /*            write (*,*) 'cwise : ', n, k, ccond, cwise_rcond, */
00481 /*     $           condthresh, ccond.ge.condthresh */
00482 /*            write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh */
00483             if (ccond >= condthresh) {
00484                 s_copy(cguar, "YES", (ftnlen)3, (ftnlen)3);
00485                 if (cwise_bnd__ > errthresh) {
00486                     tstrat[1] = 1 / (eps * 2.f);
00487                 } else {
00488                     if (cwise_bnd__ != 0.f) {
00489                         tstrat[1] = cwise_err__ / cwise_bnd__;
00490                     } else if (cwise_err__ != 0.f) {
00491                         tstrat[1] = 1 / (eps * 16.f);
00492                     } else {
00493                         tstrat[1] = 0.f;
00494                     }
00495                     if (tstrat[1] > 1.f) {
00496                         tstrat[1] = 1 / (eps * 4.f);
00497                     }
00498                 }
00499             } else {
00500                 s_copy(cguar, "NO", (ftnlen)3, (ftnlen)2);
00501                 if (cwise_bnd__ < 1.f) {
00502                     tstrat[1] = 1 / (eps * 8.f);
00503                 } else {
00504                     tstrat[1] = 1.f;
00505                 }
00506             }
00507 /*     Backwards error test */
00508             tstrat[2] = berr[k - 1] / eps;
00509 /*     Condition number tests */
00510             tstrat[3] = rcond / orcond;
00511             if (rcond >= condthresh && tstrat[3] < 1.f) {
00512                 tstrat[3] = 1.f / tstrat[3];
00513             }
00514             tstrat[4] = ncond / nwise_rcond__;
00515             if (ncond >= condthresh && tstrat[4] < 1.f) {
00516                 tstrat[4] = 1.f / tstrat[4];
00517             }
00518             tstrat[5] = ccond / nwise_rcond__;
00519             if (ccond >= condthresh && tstrat[5] < 1.f) {
00520                 tstrat[5] = 1.f / tstrat[5];
00521             }
00522             for (i__ = 1; i__ <= 6; ++i__) {
00523                 if (tstrat[i__ - 1] > *thresh) {
00524                     if (! printed_guide__) {
00525                         s_wsle(&io___66);
00526                         e_wsle();
00527                         s_wsfe(&io___67);
00528                         do_fio(&c__1, (char *)&c__1, (ftnlen)sizeof(integer));
00529                         e_wsfe();
00530                         s_wsfe(&io___68);
00531                         do_fio(&c__1, (char *)&c__2, (ftnlen)sizeof(integer));
00532                         e_wsfe();
00533                         s_wsfe(&io___69);
00534                         do_fio(&c__1, (char *)&c__3, (ftnlen)sizeof(integer));
00535                         e_wsfe();
00536                         s_wsfe(&io___70);
00537                         do_fio(&c__1, (char *)&c__4, (ftnlen)sizeof(integer));
00538                         e_wsfe();
00539                         s_wsfe(&io___71);
00540                         do_fio(&c__1, (char *)&c__5, (ftnlen)sizeof(integer));
00541                         e_wsfe();
00542                         s_wsfe(&io___72);
00543                         do_fio(&c__1, (char *)&c__6, (ftnlen)sizeof(integer));
00544                         e_wsfe();
00545                         s_wsfe(&io___73);
00546                         do_fio(&c__1, (char *)&c__7, (ftnlen)sizeof(integer));
00547                         e_wsfe();
00548                         s_wsfe(&io___74);
00549                         do_fio(&c__1, (char *)&c__8, (ftnlen)sizeof(integer));
00550                         e_wsfe();
00551                         s_wsle(&io___75);
00552                         e_wsle();
00553                         printed_guide__ = TRUE_;
00554                     }
00555                     s_wsfe(&io___76);
00556                     do_fio(&c__1, c2, (ftnlen)2);
00557                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00558                     do_fio(&c__1, (char *)&k, (ftnlen)sizeof(integer));
00559                     do_fio(&c__1, nguar, (ftnlen)3);
00560                     do_fio(&c__1, cguar, (ftnlen)3);
00561                     do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer));
00562                     do_fio(&c__1, (char *)&tstrat[i__ - 1], (ftnlen)sizeof(
00563                             real));
00564                     e_wsfe();
00565                     ++nfail;
00566                 }
00567             }
00568         }
00569 /* $$$         WRITE(*,*) */
00570 /* $$$         WRITE(*,*) 'Normwise Error Bounds' */
00571 /* $$$         WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,nwise_i,bnd_i) */
00572 /* $$$         WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,nwise_i,cond_i) */
00573 /* $$$         WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,nwise_i,rawbnd_i) */
00574 /* $$$         WRITE(*,*) */
00575 /* $$$         WRITE(*,*) 'Componentwise Error Bounds' */
00576 /* $$$         WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,cwise_i,bnd_i) */
00577 /* $$$         WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,cwise_i,cond_i) */
00578 /* $$$         WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,cwise_i,rawbnd_i) */
00579 /* $$$         print *, 'Info: ', info */
00580 /* $$$         WRITE(*,*) */
00581 /*         WRITE(*,*) 'TSTRAT: ',TSTRAT */
00582     }
00583     s_wsle(&io___77);
00584     e_wsle();
00585     if (nfail > 0) {
00586         s_wsfe(&io___78);
00587         do_fio(&c__1, c2, (ftnlen)2);
00588         do_fio(&c__1, (char *)&nfail, (ftnlen)sizeof(integer));
00589         i__1 = n * 6 + n_aux_tests__;
00590         do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer));
00591         e_wsfe();
00592     } else {
00593         s_wsfe(&io___79);
00594         do_fio(&c__1, c2, (ftnlen)2);
00595         e_wsfe();
00596     }
00597 /*     Test ratios. */
00598     return 0;
00599 } /* sebchvxx_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:01