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__1 = 1;
00034 static integer c__2 = 2;
00035 static integer c__0 = 0;
00036 static integer c_n1 = -1;
00037 static doublereal c_b50 = 0.;
00038
00039 int ddrvpo_(logical *dotype, integer *nn, integer *nval,
00040 integer *nrhs, doublereal *thresh, logical *tsterr, integer *nmax,
00041 doublereal *a, doublereal *afac, doublereal *asav, doublereal *b,
00042 doublereal *bsav, doublereal *x, doublereal *xact, doublereal *s,
00043 doublereal *work, doublereal *rwork, integer *iwork, integer *nout)
00044 {
00045
00046
00047 static integer iseedy[4] = { 1988,1989,1990,1991 };
00048 static char uplos[1*2] = "U" "L";
00049 static char facts[1*3] = "F" "N" "E";
00050 static char equeds[1*2] = "N" "Y";
00051
00052
00053 static char fmt_9999[] = "(1x,a,\002, UPLO='\002,a1,\002', N =\002,i5"
00054 ",\002, type \002,i1,\002, test(\002,i1,\002)=\002,g12.5)";
00055 static char fmt_9997[] = "(1x,a,\002, FACT='\002,a1,\002', UPLO='\002,"
00056 "a1,\002', N=\002,i5,\002, EQUED='\002,a1,\002', type \002,i1,"
00057 "\002, test(\002,i1,\002) =\002,g12.5)";
00058 static char fmt_9998[] = "(1x,a,\002, FACT='\002,a1,\002', UPLO='\002,"
00059 "a1,\002', N=\002,i5,\002, type \002,i1,\002, test(\002,i1,\002)"
00060 "=\002,g12.5)";
00061
00062
00063 address a__1[2];
00064 integer i__1, i__2, i__3, i__4, i__5[2];
00065 char ch__1[2];
00066
00067
00068 int s_copy(char *, char *, ftnlen, ftnlen);
00069 integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
00070 int s_cat(char *, char **, integer *, integer *, ftnlen);
00071
00072
00073 integer i__, k, n, k1, nb, in, kl, ku, nt, lda;
00074 char fact[1];
00075 integer ioff, mode;
00076 doublereal amax;
00077 char path[3];
00078 integer imat, info;
00079 char dist[1], uplo[1], type__[1];
00080 integer nrun, ifact;
00081 extern int dget04_(integer *, integer *, doublereal *,
00082 integer *, doublereal *, integer *, doublereal *, doublereal *);
00083 integer nfail, iseed[4], nfact;
00084 extern doublereal dget06_(doublereal *, doublereal *);
00085 extern logical lsame_(char *, char *);
00086 char equed[1];
00087 integer nbmin;
00088 doublereal rcond, roldc, scond;
00089 integer nimat;
00090 extern int dpot01_(char *, integer *, doublereal *,
00091 integer *, doublereal *, integer *, doublereal *, doublereal *), dpot02_(char *, integer *, integer *, doublereal *,
00092 integer *, doublereal *, integer *, doublereal *, integer *,
00093 doublereal *, doublereal *), dpot05_(char *, integer *,
00094 integer *, doublereal *, integer *, doublereal *, integer *,
00095 doublereal *, integer *, doublereal *, integer *, doublereal *,
00096 doublereal *, doublereal *);
00097 doublereal anorm;
00098 logical equil;
00099 integer iuplo, izero, nerrs;
00100 extern int dposv_(char *, integer *, integer *,
00101 doublereal *, integer *, doublereal *, integer *, integer *);
00102 logical zerot;
00103 char xtype[1];
00104 extern int dlatb4_(char *, integer *, integer *, integer
00105 *, char *, integer *, integer *, doublereal *, integer *,
00106 doublereal *, char *), aladhd_(integer *,
00107 char *), alaerh_(char *, char *, integer *, integer *,
00108 char *, integer *, integer *, integer *, integer *, integer *,
00109 integer *, integer *, integer *, integer *);
00110 logical prefac;
00111 doublereal rcondc;
00112 logical nofact;
00113 integer iequed;
00114 extern int dlacpy_(char *, integer *, integer *,
00115 doublereal *, integer *, doublereal *, integer *),
00116 dlarhs_(char *, char *, char *, char *, integer *, integer *,
00117 integer *, integer *, integer *, doublereal *, integer *,
00118 doublereal *, integer *, doublereal *, integer *, integer *,
00119 integer *), dlaset_(char *,
00120 integer *, integer *, doublereal *, doublereal *, doublereal *,
00121 integer *), alasvm_(char *, integer *, integer *, integer
00122 *, integer *);
00123 doublereal cndnum;
00124 extern int dlatms_(integer *, integer *, char *, integer
00125 *, char *, doublereal *, integer *, doublereal *, doublereal *,
00126 integer *, integer *, char *, doublereal *, integer *, doublereal
00127 *, integer *);
00128 doublereal ainvnm;
00129 extern doublereal dlansy_(char *, char *, integer *, doublereal *,
00130 integer *, doublereal *);
00131 extern int dlaqsy_(char *, integer *, doublereal *,
00132 integer *, doublereal *, doublereal *, doublereal *, char *), dpoequ_(integer *, doublereal *, integer *,
00133 doublereal *, doublereal *, doublereal *, integer *), dpotrf_(
00134 char *, integer *, doublereal *, integer *, integer *),
00135 dpotri_(char *, integer *, doublereal *, integer *, integer *), xlaenv_(integer *, integer *), derrvx_(char *, integer *);
00136 doublereal result[6];
00137 extern int dposvx_(char *, char *, integer *, integer *,
00138 doublereal *, integer *, doublereal *, integer *, char *,
00139 doublereal *, doublereal *, integer *, doublereal *, integer *,
00140 doublereal *, doublereal *, doublereal *, doublereal *, integer *,
00141 integer *);
00142
00143
00144 static cilist io___48 = { 0, 0, 0, fmt_9999, 0 };
00145 static cilist io___51 = { 0, 0, 0, fmt_9997, 0 };
00146 static cilist io___52 = { 0, 0, 0, fmt_9998, 0 };
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
00232
00233
00234
00235
00236
00237
00238
00239
00240 --iwork;
00241 --rwork;
00242 --work;
00243 --s;
00244 --xact;
00245 --x;
00246 --bsav;
00247 --b;
00248 --asav;
00249 --afac;
00250 --a;
00251 --nval;
00252 --dotype;
00253
00254
00255
00256
00257
00258
00259
00260 s_copy(path, "Double precision", (ftnlen)1, (ftnlen)16);
00261 s_copy(path + 1, "PO", (ftnlen)2, (ftnlen)2);
00262 nrun = 0;
00263 nfail = 0;
00264 nerrs = 0;
00265 for (i__ = 1; i__ <= 4; ++i__) {
00266 iseed[i__ - 1] = iseedy[i__ - 1];
00267
00268 }
00269
00270
00271
00272 if (*tsterr) {
00273 derrvx_(path, nout);
00274 }
00275 infoc_1.infot = 0;
00276
00277
00278
00279 nb = 1;
00280 nbmin = 2;
00281 xlaenv_(&c__1, &nb);
00282 xlaenv_(&c__2, &nbmin);
00283
00284
00285
00286 i__1 = *nn;
00287 for (in = 1; in <= i__1; ++in) {
00288 n = nval[in];
00289 lda = max(n,1);
00290 *(unsigned char *)xtype = 'N';
00291 nimat = 9;
00292 if (n <= 0) {
00293 nimat = 1;
00294 }
00295
00296 i__2 = nimat;
00297 for (imat = 1; imat <= i__2; ++imat) {
00298
00299
00300
00301 if (! dotype[imat]) {
00302 goto L120;
00303 }
00304
00305
00306
00307 zerot = imat >= 3 && imat <= 5;
00308 if (zerot && n < imat - 2) {
00309 goto L120;
00310 }
00311
00312
00313
00314 for (iuplo = 1; iuplo <= 2; ++iuplo) {
00315 *(unsigned char *)uplo = *(unsigned char *)&uplos[iuplo - 1];
00316
00317
00318
00319
00320 dlatb4_(path, &imat, &n, &n, type__, &kl, &ku, &anorm, &mode,
00321 &cndnum, dist);
00322
00323 s_copy(srnamc_1.srnamt, "DLATMS", (ftnlen)32, (ftnlen)6);
00324 dlatms_(&n, &n, dist, iseed, type__, &rwork[1], &mode, &
00325 cndnum, &anorm, &kl, &ku, uplo, &a[1], &lda, &work[1],
00326 &info);
00327
00328
00329
00330 if (info != 0) {
00331 alaerh_(path, "DLATMS", &info, &c__0, uplo, &n, &n, &c_n1,
00332 &c_n1, &c_n1, &imat, &nfail, &nerrs, nout);
00333 goto L110;
00334 }
00335
00336
00337
00338
00339 if (zerot) {
00340 if (imat == 3) {
00341 izero = 1;
00342 } else if (imat == 4) {
00343 izero = n;
00344 } else {
00345 izero = n / 2 + 1;
00346 }
00347 ioff = (izero - 1) * lda;
00348
00349
00350
00351 if (iuplo == 1) {
00352 i__3 = izero - 1;
00353 for (i__ = 1; i__ <= i__3; ++i__) {
00354 a[ioff + i__] = 0.;
00355
00356 }
00357 ioff += izero;
00358 i__3 = n;
00359 for (i__ = izero; i__ <= i__3; ++i__) {
00360 a[ioff] = 0.;
00361 ioff += lda;
00362
00363 }
00364 } else {
00365 ioff = izero;
00366 i__3 = izero - 1;
00367 for (i__ = 1; i__ <= i__3; ++i__) {
00368 a[ioff] = 0.;
00369 ioff += lda;
00370
00371 }
00372 ioff -= izero;
00373 i__3 = n;
00374 for (i__ = izero; i__ <= i__3; ++i__) {
00375 a[ioff + i__] = 0.;
00376
00377 }
00378 }
00379 } else {
00380 izero = 0;
00381 }
00382
00383
00384
00385 dlacpy_(uplo, &n, &n, &a[1], &lda, &asav[1], &lda);
00386
00387 for (iequed = 1; iequed <= 2; ++iequed) {
00388 *(unsigned char *)equed = *(unsigned char *)&equeds[
00389 iequed - 1];
00390 if (iequed == 1) {
00391 nfact = 3;
00392 } else {
00393 nfact = 1;
00394 }
00395
00396 i__3 = nfact;
00397 for (ifact = 1; ifact <= i__3; ++ifact) {
00398 *(unsigned char *)fact = *(unsigned char *)&facts[
00399 ifact - 1];
00400 prefac = lsame_(fact, "F");
00401 nofact = lsame_(fact, "N");
00402 equil = lsame_(fact, "E");
00403
00404 if (zerot) {
00405 if (prefac) {
00406 goto L90;
00407 }
00408 rcondc = 0.;
00409
00410 } else if (! lsame_(fact, "N"))
00411 {
00412
00413
00414
00415
00416
00417
00418 dlacpy_(uplo, &n, &n, &asav[1], &lda, &afac[1], &
00419 lda);
00420 if (equil || iequed > 1) {
00421
00422
00423
00424
00425 dpoequ_(&n, &afac[1], &lda, &s[1], &scond, &
00426 amax, &info);
00427 if (info == 0 && n > 0) {
00428 if (iequed > 1) {
00429 scond = 0.;
00430 }
00431
00432
00433
00434 dlaqsy_(uplo, &n, &afac[1], &lda, &s[1], &
00435 scond, &amax, equed);
00436 }
00437 }
00438
00439
00440
00441
00442 if (equil) {
00443 roldc = rcondc;
00444 }
00445
00446
00447
00448 anorm = dlansy_("1", uplo, &n, &afac[1], &lda, &
00449 rwork[1]);
00450
00451
00452
00453 dpotrf_(uplo, &n, &afac[1], &lda, &info);
00454
00455
00456
00457 dlacpy_(uplo, &n, &n, &afac[1], &lda, &a[1], &lda);
00458 dpotri_(uplo, &n, &a[1], &lda, &info);
00459
00460
00461
00462 ainvnm = dlansy_("1", uplo, &n, &a[1], &lda, &
00463 rwork[1]);
00464 if (anorm <= 0. || ainvnm <= 0.) {
00465 rcondc = 1.;
00466 } else {
00467 rcondc = 1. / anorm / ainvnm;
00468 }
00469 }
00470
00471
00472
00473 dlacpy_(uplo, &n, &n, &asav[1], &lda, &a[1], &lda);
00474
00475
00476
00477 s_copy(srnamc_1.srnamt, "DLARHS", (ftnlen)32, (ftnlen)
00478 6);
00479 dlarhs_(path, xtype, uplo, " ", &n, &n, &kl, &ku,
00480 nrhs, &a[1], &lda, &xact[1], &lda, &b[1], &
00481 lda, iseed, &info);
00482 *(unsigned char *)xtype = 'C';
00483 dlacpy_("Full", &n, nrhs, &b[1], &lda, &bsav[1], &lda);
00484
00485 if (nofact) {
00486
00487
00488
00489
00490
00491
00492 dlacpy_(uplo, &n, &n, &a[1], &lda, &afac[1], &lda);
00493 dlacpy_("Full", &n, nrhs, &b[1], &lda, &x[1], &
00494 lda);
00495
00496 s_copy(srnamc_1.srnamt, "DPOSV ", (ftnlen)32, (
00497 ftnlen)6);
00498 dposv_(uplo, &n, nrhs, &afac[1], &lda, &x[1], &
00499 lda, &info);
00500
00501
00502
00503 if (info != izero) {
00504 alaerh_(path, "DPOSV ", &info, &izero, uplo, &
00505 n, &n, &c_n1, &c_n1, nrhs, &imat, &
00506 nfail, &nerrs, nout);
00507 goto L70;
00508 } else if (info != 0) {
00509 goto L70;
00510 }
00511
00512
00513
00514
00515 dpot01_(uplo, &n, &a[1], &lda, &afac[1], &lda, &
00516 rwork[1], result);
00517
00518
00519
00520 dlacpy_("Full", &n, nrhs, &b[1], &lda, &work[1], &
00521 lda);
00522 dpot02_(uplo, &n, nrhs, &a[1], &lda, &x[1], &lda,
00523 &work[1], &lda, &rwork[1], &result[1]);
00524
00525
00526
00527 dget04_(&n, nrhs, &x[1], &lda, &xact[1], &lda, &
00528 rcondc, &result[2]);
00529 nt = 3;
00530
00531
00532
00533
00534 i__4 = nt;
00535 for (k = 1; k <= i__4; ++k) {
00536 if (result[k - 1] >= *thresh) {
00537 if (nfail == 0 && nerrs == 0) {
00538 aladhd_(nout, path);
00539 }
00540 io___48.ciunit = *nout;
00541 s_wsfe(&io___48);
00542 do_fio(&c__1, "DPOSV ", (ftnlen)6);
00543 do_fio(&c__1, uplo, (ftnlen)1);
00544 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(
00545 integer));
00546 do_fio(&c__1, (char *)&imat, (ftnlen)
00547 sizeof(integer));
00548 do_fio(&c__1, (char *)&k, (ftnlen)sizeof(
00549 integer));
00550 do_fio(&c__1, (char *)&result[k - 1], (
00551 ftnlen)sizeof(doublereal));
00552 e_wsfe();
00553 ++nfail;
00554 }
00555
00556 }
00557 nrun += nt;
00558 L70:
00559 ;
00560 }
00561
00562
00563
00564 if (! prefac) {
00565 dlaset_(uplo, &n, &n, &c_b50, &c_b50, &afac[1], &
00566 lda);
00567 }
00568 dlaset_("Full", &n, nrhs, &c_b50, &c_b50, &x[1], &lda);
00569 if (iequed > 1 && n > 0) {
00570
00571
00572
00573
00574 dlaqsy_(uplo, &n, &a[1], &lda, &s[1], &scond, &
00575 amax, equed);
00576 }
00577
00578
00579
00580
00581 s_copy(srnamc_1.srnamt, "DPOSVX", (ftnlen)32, (ftnlen)
00582 6);
00583 dposvx_(fact, uplo, &n, nrhs, &a[1], &lda, &afac[1], &
00584 lda, equed, &s[1], &b[1], &lda, &x[1], &lda, &
00585 rcond, &rwork[1], &rwork[*nrhs + 1], &work[1],
00586 &iwork[1], &info);
00587
00588
00589
00590 if (info != izero) {
00591
00592 i__5[0] = 1, a__1[0] = fact;
00593 i__5[1] = 1, a__1[1] = uplo;
00594 s_cat(ch__1, a__1, i__5, &c__2, (ftnlen)2);
00595 alaerh_(path, "DPOSVX", &info, &izero, ch__1, &n,
00596 &n, &c_n1, &c_n1, nrhs, &imat, &nfail, &
00597 nerrs, nout);
00598 goto L90;
00599 }
00600
00601 if (info == 0) {
00602 if (! prefac) {
00603
00604
00605
00606
00607 dpot01_(uplo, &n, &a[1], &lda, &afac[1], &lda,
00608 &rwork[(*nrhs << 1) + 1], result);
00609 k1 = 1;
00610 } else {
00611 k1 = 2;
00612 }
00613
00614
00615
00616 dlacpy_("Full", &n, nrhs, &bsav[1], &lda, &work[1]
00617 , &lda);
00618 dpot02_(uplo, &n, nrhs, &asav[1], &lda, &x[1], &
00619 lda, &work[1], &lda, &rwork[(*nrhs << 1)
00620 + 1], &result[1]);
00621
00622
00623
00624 if (nofact || prefac && lsame_(equed, "N")) {
00625 dget04_(&n, nrhs, &x[1], &lda, &xact[1], &lda,
00626 &rcondc, &result[2]);
00627 } else {
00628 dget04_(&n, nrhs, &x[1], &lda, &xact[1], &lda,
00629 &roldc, &result[2]);
00630 }
00631
00632
00633
00634
00635 dpot05_(uplo, &n, nrhs, &asav[1], &lda, &b[1], &
00636 lda, &x[1], &lda, &xact[1], &lda, &rwork[
00637 1], &rwork[*nrhs + 1], &result[3]);
00638 } else {
00639 k1 = 6;
00640 }
00641
00642
00643
00644
00645 result[5] = dget06_(&rcond, &rcondc);
00646
00647
00648
00649
00650 for (k = k1; k <= 6; ++k) {
00651 if (result[k - 1] >= *thresh) {
00652 if (nfail == 0 && nerrs == 0) {
00653 aladhd_(nout, path);
00654 }
00655 if (prefac) {
00656 io___51.ciunit = *nout;
00657 s_wsfe(&io___51);
00658 do_fio(&c__1, "DPOSVX", (ftnlen)6);
00659 do_fio(&c__1, fact, (ftnlen)1);
00660 do_fio(&c__1, uplo, (ftnlen)1);
00661 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(
00662 integer));
00663 do_fio(&c__1, equed, (ftnlen)1);
00664 do_fio(&c__1, (char *)&imat, (ftnlen)
00665 sizeof(integer));
00666 do_fio(&c__1, (char *)&k, (ftnlen)sizeof(
00667 integer));
00668 do_fio(&c__1, (char *)&result[k - 1], (
00669 ftnlen)sizeof(doublereal));
00670 e_wsfe();
00671 } else {
00672 io___52.ciunit = *nout;
00673 s_wsfe(&io___52);
00674 do_fio(&c__1, "DPOSVX", (ftnlen)6);
00675 do_fio(&c__1, fact, (ftnlen)1);
00676 do_fio(&c__1, uplo, (ftnlen)1);
00677 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(
00678 integer));
00679 do_fio(&c__1, (char *)&imat, (ftnlen)
00680 sizeof(integer));
00681 do_fio(&c__1, (char *)&k, (ftnlen)sizeof(
00682 integer));
00683 do_fio(&c__1, (char *)&result[k - 1], (
00684 ftnlen)sizeof(doublereal));
00685 e_wsfe();
00686 }
00687 ++nfail;
00688 }
00689
00690 }
00691 nrun = nrun + 7 - k1;
00692 L90:
00693 ;
00694 }
00695
00696 }
00697 L110:
00698 ;
00699 }
00700 L120:
00701 ;
00702 }
00703
00704 }
00705
00706
00707
00708 alasvm_(path, nout, &nfail, &nrun, &nerrs);
00709
00710 return 0;
00711
00712
00713
00714 }