00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015
00016
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 int sebchvxx_(real *thresh, char *path)
00028 {
00029
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
00058 integer i__1, i__2, i__3, i__4, i__5, i__6;
00059 real r__1, r__2, r__3;
00060
00061
00062 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
00068 extern 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] , b[36]
00077 , c__[6];
00078 integer i__, j, k;
00079 real m;
00080 integer n;
00081 real r__[6], s[6], x[36] , cwise_bnd__;
00082 char c2[2];
00083 real nwise_bnd__, cwise_err__, nwise_err__, errthresh, ab[66]
00084 , af[36] ;
00085 integer kl, ku;
00086 real condthresh, afb[96] ;
00087 integer lda;
00088 real eps, cwise_rcond__, nwise_rcond__;
00089 integer n_aux_tests__, ldab;
00090 real diff[36] ;
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] ;
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] ;
00110 extern logical lsamen_(integer *, char *, char *);
00111 real params[2], orcond;
00112 extern int slacpy_(char *, integer *, integer *, real *,
00113 integer *, real *, integer *);
00114 real rinorm, tstrat[6], rpvgrw;
00115 extern int slahilb_(integer *, integer *, real *,
00116 integer *, real *, integer *, real *, integer *, real *, integer *
00117 );
00118 real invhilb[36] , normdif;
00119 extern 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
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
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
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
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
00236 r__1 = sqrt((real) n);
00237 m = dmax(r__1,10.f);
00238
00239
00240 slahilb_(&n, &n, a, &lda, invhilb, &lda, b, &lda, work, &info);
00241
00242 slacpy_("ALL", &n, &n, a, &c__6, acopy, &c__6);
00243
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
00254 i__2 = 1, i__3 = j - ku;
00255
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
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
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
00296
00297
00298
00299
00300
00301 } else {
00302
00303
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
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
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
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
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
00387
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
00399 r__2 = (r__1 = invhilb[i__ + k * 6 - 7], dabs(r__1));
00400 normt = dmax(r__2,normt);
00401
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
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
00446
00447 ccond = a[0] / rinorm;
00448
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
00454
00455
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
00481
00482
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
00508 tstrat[2] = berr[k - 1] / eps;
00509
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
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
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
00598 return 0;
00599 }