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