00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015
00016
00017
00018 struct {
00019 integer infot, nunit;
00020 logical ok, lerr;
00021 } infoc_;
00022
00023 #define infoc_1 infoc_
00024
00025 struct {
00026 char srnamt[32];
00027 } srnamc_;
00028
00029 #define srnamc_1 srnamc_
00030
00031
00032
00033 static integer c__2 = 2;
00034 static integer c__0 = 0;
00035 static integer c_n1 = -1;
00036 static integer c__1 = 1;
00037 static real c_b50 = 0.f;
00038 static real c_b51 = 1.f;
00039 static integer c__7 = 7;
00040
00041 int schkpb_(logical *dotype, integer *nn, integer *nval,
00042 integer *nnb, integer *nbval, integer *nns, integer *nsval, real *
00043 thresh, logical *tsterr, integer *nmax, real *a, real *afac, real *
00044 ainv, real *b, real *x, real *xact, real *work, real *rwork, integer *
00045 iwork, integer *nout)
00046 {
00047
00048
00049 static integer iseedy[4] = { 1988,1989,1990,1991 };
00050
00051
00052 static char fmt_9999[] = "(\002 UPLO='\002,a1,\002', N=\002,i5,\002, KD"
00053 "=\002,i5,\002, NB=\002,i4,\002, type \002,i2,\002, test \002,i2"
00054 ",\002, ratio= \002,g12.5)";
00055 static char fmt_9998[] = "(\002 UPLO='\002,a1,\002', N=\002,i5,\002, KD"
00056 "=\002,i5,\002, NRHS=\002,i3,\002, type \002,i2,\002, test(\002,i"
00057 "2,\002) = \002,g12.5)";
00058 static char fmt_9997[] = "(\002 UPLO='\002,a1,\002', N=\002,i5,\002, KD"
00059 "=\002,i5,\002,\002,10x,\002 type \002,i2,\002, test(\002,i2,\002"
00060 ") = \002,g12.5)";
00061
00062
00063 integer i__1, i__2, i__3, i__4, i__5, i__6;
00064
00065
00066 int s_copy(char *, char *, ftnlen, ftnlen);
00067 integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
00068
00069
00070 integer i__, k, n, i1, i2, kd, nb, in, kl, iw, ku, lda, ikd, inb, nkd,
00071 ldab, ioff, mode, koff, imat, info;
00072 char path[3], dist[1];
00073 integer irhs, nrhs;
00074 char uplo[1], type__[1];
00075 integer nrun;
00076 extern int alahd_(integer *, char *);
00077 integer nfail, iseed[4], kdval[4];
00078 real rcond;
00079 extern int sget04_(integer *, integer *, real *, integer
00080 *, real *, integer *, real *, real *);
00081 integer nimat;
00082 extern doublereal sget06_(real *, real *);
00083 extern int spbt01_(char *, integer *, integer *, real *,
00084 integer *, real *, integer *, real *, real *), spbt02_(
00085 char *, integer *, integer *, integer *, real *, integer *, real *
00086 , integer *, real *, integer *, real *, real *);
00087 real anorm;
00088 extern int spbt05_(char *, integer *, integer *, integer
00089 *, real *, integer *, real *, integer *, real *, integer *, real *
00090 , integer *, real *, real *, real *);
00091 integer iuplo, izero, nerrs;
00092 extern int scopy_(integer *, real *, integer *, real *,
00093 integer *), sswap_(integer *, real *, integer *, real *, integer *
00094 );
00095 logical zerot;
00096 char xtype[1];
00097 extern int slatb4_(char *, integer *, integer *, integer
00098 *, char *, integer *, integer *, real *, integer *, real *, char *
00099 ), alaerh_(char *, char *, integer *,
00100 integer *, char *, integer *, integer *, integer *, integer *,
00101 integer *, integer *, integer *, integer *, integer *);
00102 real rcondc;
00103 extern doublereal slange_(char *, integer *, integer *, real *, integer *,
00104 real *);
00105 char packit[1];
00106 extern doublereal slansb_(char *, char *, integer *, integer *, real *,
00107 integer *, real *);
00108 real cndnum;
00109 extern int alasum_(char *, integer *, integer *, integer
00110 *, integer *), spbcon_(char *, integer *, integer *, real
00111 *, integer *, real *, real *, real *, integer *, integer *);
00112 real ainvnm;
00113 extern int slacpy_(char *, integer *, integer *, real *,
00114 integer *, real *, integer *), slarhs_(char *, char *,
00115 char *, char *, integer *, integer *, integer *, integer *,
00116 integer *, real *, integer *, real *, integer *, real *, integer *
00117 , integer *, integer *), slaset_(
00118 char *, integer *, integer *, real *, real *, real *, integer *), spbrfs_(char *, integer *, integer *, integer *, real *,
00119 integer *, real *, integer *, real *, integer *, real *, integer *
00120 , real *, real *, real *, integer *, integer *), spbtrf_(
00121 char *, integer *, integer *, real *, integer *, integer *), xlaenv_(integer *, integer *), slatms_(integer *,
00122 integer *, char *, integer *, char *, real *, integer *, real *,
00123 real *, integer *, integer *, char *, real *, integer *, real *,
00124 integer *), serrpo_(char *, integer *), spbtrs_(char *, integer *, integer *, integer *, real *,
00125 integer *, real *, integer *, integer *);
00126 real result[7];
00127
00128
00129 static cilist io___40 = { 0, 0, 0, fmt_9999, 0 };
00130 static cilist io___46 = { 0, 0, 0, fmt_9998, 0 };
00131 static cilist io___48 = { 0, 0, 0, fmt_9997, 0 };
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
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
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231 --iwork;
00232 --rwork;
00233 --work;
00234 --xact;
00235 --x;
00236 --b;
00237 --ainv;
00238 --afac;
00239 --a;
00240 --nsval;
00241 --nbval;
00242 --nval;
00243 --dotype;
00244
00245
00246
00247
00248
00249
00250
00251 s_copy(path, "Single precision", (ftnlen)1, (ftnlen)16);
00252 s_copy(path + 1, "PB", (ftnlen)2, (ftnlen)2);
00253 nrun = 0;
00254 nfail = 0;
00255 nerrs = 0;
00256 for (i__ = 1; i__ <= 4; ++i__) {
00257 iseed[i__ - 1] = iseedy[i__ - 1];
00258
00259 }
00260
00261
00262
00263 if (*tsterr) {
00264 serrpo_(path, nout);
00265 }
00266 infoc_1.infot = 0;
00267 xlaenv_(&c__2, &c__2);
00268 kdval[0] = 0;
00269
00270
00271
00272 i__1 = *nn;
00273 for (in = 1; in <= i__1; ++in) {
00274 n = nval[in];
00275 lda = max(n,1);
00276 *(unsigned char *)xtype = 'N';
00277
00278
00279
00280
00281 i__2 = 1, i__3 = min(n,4);
00282 nkd = max(i__2,i__3);
00283 nimat = 8;
00284 if (n == 0) {
00285 nimat = 1;
00286 }
00287
00288 kdval[1] = n + (n + 1) / 4;
00289 kdval[2] = (n * 3 - 1) / 4;
00290 kdval[3] = (n + 1) / 4;
00291
00292 i__2 = nkd;
00293 for (ikd = 1; ikd <= i__2; ++ikd) {
00294
00295
00296
00297
00298
00299 kd = kdval[ikd - 1];
00300 ldab = kd + 1;
00301
00302
00303
00304 for (iuplo = 1; iuplo <= 2; ++iuplo) {
00305 koff = 1;
00306 if (iuplo == 1) {
00307 *(unsigned char *)uplo = 'U';
00308
00309 i__3 = 1, i__4 = kd + 2 - n;
00310 koff = max(i__3,i__4);
00311 *(unsigned char *)packit = 'Q';
00312 } else {
00313 *(unsigned char *)uplo = 'L';
00314 *(unsigned char *)packit = 'B';
00315 }
00316
00317 i__3 = nimat;
00318 for (imat = 1; imat <= i__3; ++imat) {
00319
00320
00321
00322 if (! dotype[imat]) {
00323 goto L60;
00324 }
00325
00326
00327
00328 zerot = imat >= 2 && imat <= 4;
00329 if (zerot && n < imat - 1) {
00330 goto L60;
00331 }
00332
00333 if (! zerot || ! dotype[1]) {
00334
00335
00336
00337
00338 slatb4_(path, &imat, &n, &n, type__, &kl, &ku, &anorm,
00339 &mode, &cndnum, dist);
00340
00341 s_copy(srnamc_1.srnamt, "SLATMS", (ftnlen)32, (ftnlen)
00342 6);
00343 slatms_(&n, &n, dist, iseed, type__, &rwork[1], &mode,
00344 &cndnum, &anorm, &kd, &kd, packit, &a[koff],
00345 &ldab, &work[1], &info);
00346
00347
00348
00349 if (info != 0) {
00350 alaerh_(path, "SLATMS", &info, &c__0, uplo, &n, &
00351 n, &kd, &kd, &c_n1, &imat, &nfail, &nerrs,
00352 nout);
00353 goto L60;
00354 }
00355 } else if (izero > 0) {
00356
00357
00358
00359
00360 iw = (lda << 1) + 1;
00361 if (iuplo == 1) {
00362 ioff = (izero - 1) * ldab + kd + 1;
00363 i__4 = izero - i1;
00364 scopy_(&i__4, &work[iw], &c__1, &a[ioff - izero +
00365 i1], &c__1);
00366 iw = iw + izero - i1;
00367 i__4 = i2 - izero + 1;
00368
00369 i__6 = ldab - 1;
00370 i__5 = max(i__6,1);
00371 scopy_(&i__4, &work[iw], &c__1, &a[ioff], &i__5);
00372 } else {
00373 ioff = (i1 - 1) * ldab + 1;
00374 i__4 = izero - i1;
00375
00376 i__6 = ldab - 1;
00377 i__5 = max(i__6,1);
00378 scopy_(&i__4, &work[iw], &c__1, &a[ioff + izero -
00379 i1], &i__5);
00380 ioff = (izero - 1) * ldab + 1;
00381 iw = iw + izero - i1;
00382 i__4 = i2 - izero + 1;
00383 scopy_(&i__4, &work[iw], &c__1, &a[ioff], &c__1);
00384 }
00385 }
00386
00387
00388
00389
00390 izero = 0;
00391 if (zerot) {
00392 if (imat == 2) {
00393 izero = 1;
00394 } else if (imat == 3) {
00395 izero = n;
00396 } else {
00397 izero = n / 2 + 1;
00398 }
00399
00400
00401
00402 iw = lda << 1;
00403
00404 i__5 = (kd << 1) + 1;
00405 i__4 = min(i__5,n);
00406 for (i__ = 1; i__ <= i__4; ++i__) {
00407 work[iw + i__] = 0.f;
00408
00409 }
00410 ++iw;
00411
00412 i__4 = izero - kd;
00413 i1 = max(i__4,1);
00414
00415 i__4 = izero + kd;
00416 i2 = min(i__4,n);
00417
00418 if (iuplo == 1) {
00419 ioff = (izero - 1) * ldab + kd + 1;
00420 i__4 = izero - i1;
00421 sswap_(&i__4, &a[ioff - izero + i1], &c__1, &work[
00422 iw], &c__1);
00423 iw = iw + izero - i1;
00424 i__4 = i2 - izero + 1;
00425
00426 i__6 = ldab - 1;
00427 i__5 = max(i__6,1);
00428 sswap_(&i__4, &a[ioff], &i__5, &work[iw], &c__1);
00429 } else {
00430 ioff = (i1 - 1) * ldab + 1;
00431 i__4 = izero - i1;
00432
00433 i__6 = ldab - 1;
00434 i__5 = max(i__6,1);
00435 sswap_(&i__4, &a[ioff + izero - i1], &i__5, &work[
00436 iw], &c__1);
00437 ioff = (izero - 1) * ldab + 1;
00438 iw = iw + izero - i1;
00439 i__4 = i2 - izero + 1;
00440 sswap_(&i__4, &a[ioff], &c__1, &work[iw], &c__1);
00441 }
00442 }
00443
00444
00445
00446 i__4 = *nnb;
00447 for (inb = 1; inb <= i__4; ++inb) {
00448 nb = nbval[inb];
00449 xlaenv_(&c__1, &nb);
00450
00451
00452
00453
00454 i__5 = kd + 1;
00455 slacpy_("Full", &i__5, &n, &a[1], &ldab, &afac[1], &
00456 ldab);
00457 s_copy(srnamc_1.srnamt, "SPBTRF", (ftnlen)32, (ftnlen)
00458 6);
00459 spbtrf_(uplo, &n, &kd, &afac[1], &ldab, &info);
00460
00461
00462
00463 if (info != izero) {
00464 alaerh_(path, "SPBTRF", &info, &izero, uplo, &n, &
00465 n, &kd, &kd, &nb, &imat, &nfail, &nerrs,
00466 nout);
00467 goto L50;
00468 }
00469
00470
00471
00472 if (info != 0) {
00473 goto L50;
00474 }
00475
00476
00477
00478
00479
00480 i__5 = kd + 1;
00481 slacpy_("Full", &i__5, &n, &afac[1], &ldab, &ainv[1],
00482 &ldab);
00483 spbt01_(uplo, &n, &kd, &a[1], &ldab, &ainv[1], &ldab,
00484 &rwork[1], result);
00485
00486
00487
00488 if (result[0] >= *thresh) {
00489 if (nfail == 0 && nerrs == 0) {
00490 alahd_(nout, path);
00491 }
00492 io___40.ciunit = *nout;
00493 s_wsfe(&io___40);
00494 do_fio(&c__1, uplo, (ftnlen)1);
00495 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer))
00496 ;
00497 do_fio(&c__1, (char *)&kd, (ftnlen)sizeof(integer)
00498 );
00499 do_fio(&c__1, (char *)&nb, (ftnlen)sizeof(integer)
00500 );
00501 do_fio(&c__1, (char *)&imat, (ftnlen)sizeof(
00502 integer));
00503 do_fio(&c__1, (char *)&c__1, (ftnlen)sizeof(
00504 integer));
00505 do_fio(&c__1, (char *)&result[0], (ftnlen)sizeof(
00506 real));
00507 e_wsfe();
00508 ++nfail;
00509 }
00510 ++nrun;
00511
00512
00513
00514 if (inb > 1) {
00515 goto L50;
00516 }
00517
00518
00519
00520
00521 slaset_("Full", &n, &n, &c_b50, &c_b51, &ainv[1], &
00522 lda);
00523 s_copy(srnamc_1.srnamt, "SPBTRS", (ftnlen)32, (ftnlen)
00524 6);
00525 spbtrs_(uplo, &n, &kd, &n, &afac[1], &ldab, &ainv[1],
00526 &lda, &info);
00527
00528
00529
00530 anorm = slansb_("1", uplo, &n, &kd, &a[1], &ldab, &
00531 rwork[1]);
00532 ainvnm = slange_("1", &n, &n, &ainv[1], &lda, &rwork[
00533 1]);
00534 if (anorm <= 0.f || ainvnm <= 0.f) {
00535 rcondc = 1.f;
00536 } else {
00537 rcondc = 1.f / anorm / ainvnm;
00538 }
00539
00540 i__5 = *nns;
00541 for (irhs = 1; irhs <= i__5; ++irhs) {
00542 nrhs = nsval[irhs];
00543
00544
00545
00546
00547 s_copy(srnamc_1.srnamt, "SLARHS", (ftnlen)32, (
00548 ftnlen)6);
00549 slarhs_(path, xtype, uplo, " ", &n, &n, &kd, &kd,
00550 &nrhs, &a[1], &ldab, &xact[1], &lda, &b[1]
00551 , &lda, iseed, &info);
00552 slacpy_("Full", &n, &nrhs, &b[1], &lda, &x[1], &
00553 lda);
00554
00555 s_copy(srnamc_1.srnamt, "SPBTRS", (ftnlen)32, (
00556 ftnlen)6);
00557 spbtrs_(uplo, &n, &kd, &nrhs, &afac[1], &ldab, &x[
00558 1], &lda, &info);
00559
00560
00561
00562 if (info != 0) {
00563 alaerh_(path, "SPBTRS", &info, &c__0, uplo, &
00564 n, &n, &kd, &kd, &nrhs, &imat, &nfail,
00565 &nerrs, nout);
00566 }
00567
00568 slacpy_("Full", &n, &nrhs, &b[1], &lda, &work[1],
00569 &lda);
00570 spbt02_(uplo, &n, &kd, &nrhs, &a[1], &ldab, &x[1],
00571 &lda, &work[1], &lda, &rwork[1], &result[
00572 1]);
00573
00574
00575
00576
00577 sget04_(&n, &nrhs, &x[1], &lda, &xact[1], &lda, &
00578 rcondc, &result[2]);
00579
00580
00581
00582
00583 s_copy(srnamc_1.srnamt, "SPBRFS", (ftnlen)32, (
00584 ftnlen)6);
00585 spbrfs_(uplo, &n, &kd, &nrhs, &a[1], &ldab, &afac[
00586 1], &ldab, &b[1], &lda, &x[1], &lda, &
00587 rwork[1], &rwork[nrhs + 1], &work[1], &
00588 iwork[1], &info);
00589
00590
00591
00592 if (info != 0) {
00593 alaerh_(path, "SPBRFS", &info, &c__0, uplo, &
00594 n, &n, &kd, &kd, &nrhs, &imat, &nfail,
00595 &nerrs, nout);
00596 }
00597
00598 sget04_(&n, &nrhs, &x[1], &lda, &xact[1], &lda, &
00599 rcondc, &result[3]);
00600 spbt05_(uplo, &n, &kd, &nrhs, &a[1], &ldab, &b[1],
00601 &lda, &x[1], &lda, &xact[1], &lda, &
00602 rwork[1], &rwork[nrhs + 1], &result[4]);
00603
00604
00605
00606
00607 for (k = 2; k <= 6; ++k) {
00608 if (result[k - 1] >= *thresh) {
00609 if (nfail == 0 && nerrs == 0) {
00610 alahd_(nout, path);
00611 }
00612 io___46.ciunit = *nout;
00613 s_wsfe(&io___46);
00614 do_fio(&c__1, uplo, (ftnlen)1);
00615 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(
00616 integer));
00617 do_fio(&c__1, (char *)&kd, (ftnlen)sizeof(
00618 integer));
00619 do_fio(&c__1, (char *)&nrhs, (ftnlen)
00620 sizeof(integer));
00621 do_fio(&c__1, (char *)&imat, (ftnlen)
00622 sizeof(integer));
00623 do_fio(&c__1, (char *)&k, (ftnlen)sizeof(
00624 integer));
00625 do_fio(&c__1, (char *)&result[k - 1], (
00626 ftnlen)sizeof(real));
00627 e_wsfe();
00628 ++nfail;
00629 }
00630
00631 }
00632 nrun += 5;
00633
00634 }
00635
00636
00637
00638
00639 s_copy(srnamc_1.srnamt, "SPBCON", (ftnlen)32, (ftnlen)
00640 6);
00641 spbcon_(uplo, &n, &kd, &afac[1], &ldab, &anorm, &
00642 rcond, &work[1], &iwork[1], &info);
00643
00644
00645
00646 if (info != 0) {
00647 alaerh_(path, "SPBCON", &info, &c__0, uplo, &n, &
00648 n, &kd, &kd, &c_n1, &imat, &nfail, &nerrs,
00649 nout);
00650 }
00651
00652 result[6] = sget06_(&rcond, &rcondc);
00653
00654
00655
00656 if (result[6] >= *thresh) {
00657 if (nfail == 0 && nerrs == 0) {
00658 alahd_(nout, path);
00659 }
00660 io___48.ciunit = *nout;
00661 s_wsfe(&io___48);
00662 do_fio(&c__1, uplo, (ftnlen)1);
00663 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer))
00664 ;
00665 do_fio(&c__1, (char *)&kd, (ftnlen)sizeof(integer)
00666 );
00667 do_fio(&c__1, (char *)&imat, (ftnlen)sizeof(
00668 integer));
00669 do_fio(&c__1, (char *)&c__7, (ftnlen)sizeof(
00670 integer));
00671 do_fio(&c__1, (char *)&result[6], (ftnlen)sizeof(
00672 real));
00673 e_wsfe();
00674 ++nfail;
00675 }
00676 ++nrun;
00677 L50:
00678 ;
00679 }
00680 L60:
00681 ;
00682 }
00683
00684 }
00685
00686 }
00687
00688 }
00689
00690
00691
00692 alasum_(path, nout, &nfail, &nrun, &nerrs);
00693
00694 return 0;
00695
00696
00697
00698 }