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__1 = 1;
00019 static integer c__0 = 0;
00020 static integer c__5 = 5;
00021 static logical c_true = TRUE_;
00022 static logical c_false = FALSE_;
00023 static integer c__3 = 3;
00024
00025 int ddrgvx_(integer *nsize, doublereal *thresh, integer *nin,
00026 integer *nout, doublereal *a, integer *lda, doublereal *b,
00027 doublereal *ai, doublereal *bi, doublereal *alphar, doublereal *
00028 alphai, doublereal *beta, doublereal *vl, doublereal *vr, integer *
00029 ilo, integer *ihi, doublereal *lscale, doublereal *rscale, doublereal
00030 *s, doublereal *dtru, doublereal *dif, doublereal *diftru, doublereal
00031 *work, integer *lwork, integer *iwork, integer *liwork, doublereal *
00032 result, logical *bwork, integer *info)
00033 {
00034
00035 static char fmt_9999[] = "(\002 DDRGVX: \002,a,\002 returned INFO=\002,i"
00036 "6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002)\002)";
00037 static char fmt_9998[] = "(\002 DDRGVX: \002,a,\002 Eigenvectors from"
00038 " \002,a,\002 incorrectly \002,\002normalized.\002,/\002 Bits of "
00039 "error=\002,0p,g10.3,\002,\002,9x,\002N=\002,i6,\002, JTYPE=\002,"
00040 "i6,\002, IWA=\002,i5,\002, IWB=\002,i5,\002, IWX=\002,i5,\002, I"
00041 "WY=\002,i5)";
00042 static char fmt_9997[] = "(/1x,a3,\002 -- Real Expert Eigenvalue/vecto"
00043 "r\002,\002 problem driver\002)";
00044 static char fmt_9995[] = "(\002 Matrix types: \002,/)";
00045 static char fmt_9994[] = "(\002 TYPE 1: Da is diagonal, Db is identity,"
00046 " \002,/\002 A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) \002,/"
00047 "\002 YH and X are left and right eigenvectors. \002,/)";
00048 static char fmt_9993[] = "(\002 TYPE 2: Da is quasi-diagonal, Db is iden"
00049 "tity, \002,/\002 A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1)"
00050 " \002,/\002 YH and X are left and right eigenvectors. \002,/)"
00051 ;
00052 static char fmt_9992[] = "(/\002 Tests performed: \002,/4x,\002 a is al"
00053 "pha, b is beta, l is a left eigenvector, \002,/4x,\002 r is a ri"
00054 "ght eigenvector and \002,a,\002 means \002,a,\002.\002,/\002 1 ="
00055 " max | ( b A - a B )\002,a,\002 l | / const.\002,/\002 2 = max |"
00056 " ( b A - a B ) r | / const.\002,/\002 3 = max ( Sest/Stru, Stru/"
00057 "Sest ) \002,\002 over all eigenvalues\002,/\002 4 = max( DIFest/"
00058 "DIFtru, DIFtru/DIFest ) \002,\002 over the 1st and 5th eigenvect"
00059 "ors\002,/)";
00060 static char fmt_9991[] = "(\002 Type=\002,i2,\002,\002,\002 IWA=\002,i2"
00061 ",\002, IWB=\002,i2,\002, IWX=\002,i2,\002, IWY=\002,i2,\002, res"
00062 "ult \002,i2,\002 is\002,0p,f8.2)";
00063 static char fmt_9990[] = "(\002 Type=\002,i2,\002,\002,\002 IWA=\002,i2"
00064 ",\002, IWB=\002,i2,\002, IWX=\002,i2,\002, IWY=\002,i2,\002, res"
00065 "ult \002,i2,\002 is\002,1p,d10.3)";
00066 static char fmt_9987[] = "(\002 DDRGVX: \002,a,\002 returned INFO=\002,i"
00067 "6,\002.\002,/9x,\002N=\002,i6,\002, Input example #\002,i2,\002"
00068 ")\002)";
00069 static char fmt_9986[] = "(\002 DDRGVX: \002,a,\002 Eigenvectors from"
00070 " \002,a,\002 incorrectly \002,\002normalized.\002,/\002 Bits of "
00071 "error=\002,0p,g10.3,\002,\002,9x,\002N=\002,i6,\002, Input Examp"
00072 "le #\002,i2,\002)\002)";
00073 static char fmt_9996[] = "(\002 Input Example\002)";
00074 static char fmt_9989[] = "(\002 Input example #\002,i2,\002, matrix orde"
00075 "r=\002,i4,\002,\002,\002 result \002,i2,\002 is\002,0p,f8.2)";
00076 static char fmt_9988[] = "(\002 Input example #\002,i2,\002, matrix orde"
00077 "r=\002,i4,\002,\002,\002 result \002,i2,\002 is\002,1p,d10.3)";
00078
00079
00080 integer a_dim1, a_offset, ai_dim1, ai_offset, b_dim1, b_offset, bi_dim1,
00081 bi_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1, i__2;
00082 doublereal d__1, d__2, d__3, d__4;
00083
00084
00085 double sqrt(doublereal);
00086 integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void),
00087 s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
00088 e_rsle(void);
00089
00090
00091 integer i__, j, n, iwa, iwb;
00092 doublereal ulp;
00093 integer iwx, iwy, nmax;
00094 extern int dget52_(logical *, integer *, doublereal *,
00095 integer *, doublereal *, integer *, doublereal *, integer *,
00096 doublereal *, doublereal *, doublereal *, doublereal *,
00097 doublereal *);
00098 integer linfo;
00099 doublereal anorm, bnorm;
00100 integer nerrs;
00101 extern int dlatm6_(integer *, integer *, doublereal *,
00102 integer *, doublereal *, doublereal *, integer *, doublereal *,
00103 integer *, doublereal *, doublereal *, doublereal *, doublereal *,
00104 doublereal *, doublereal *);
00105 doublereal ratio1, ratio2, thrsh2;
00106 extern doublereal dlamch_(char *), dlange_(char *, integer *,
00107 integer *, doublereal *, integer *, doublereal *);
00108 extern int dlacpy_(char *, integer *, integer *,
00109 doublereal *, integer *, doublereal *, integer *),
00110 xerbla_(char *, integer *);
00111 doublereal abnorm;
00112 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
00113 integer *, integer *);
00114 extern int alasvm_(char *, integer *, integer *, integer
00115 *, integer *), dggevx_(char *, char *, char *, char *,
00116 integer *, doublereal *, integer *, doublereal *, integer *,
00117 doublereal *, doublereal *, doublereal *, doublereal *, integer *,
00118 doublereal *, integer *, integer *, integer *, doublereal *,
00119 doublereal *, doublereal *, doublereal *, doublereal *,
00120 doublereal *, doublereal *, integer *, integer *, logical *,
00121 integer *);
00122 doublereal weight[5];
00123 integer minwrk, maxwrk, iptype;
00124 doublereal ulpinv;
00125 integer nptknt, ntestt;
00126
00127
00128 static cilist io___20 = { 0, 0, 0, fmt_9999, 0 };
00129 static cilist io___22 = { 0, 0, 0, fmt_9998, 0 };
00130 static cilist io___23 = { 0, 0, 0, fmt_9998, 0 };
00131 static cilist io___28 = { 0, 0, 0, fmt_9997, 0 };
00132 static cilist io___29 = { 0, 0, 0, fmt_9995, 0 };
00133 static cilist io___30 = { 0, 0, 0, fmt_9994, 0 };
00134 static cilist io___31 = { 0, 0, 0, fmt_9993, 0 };
00135 static cilist io___32 = { 0, 0, 0, fmt_9992, 0 };
00136 static cilist io___33 = { 0, 0, 0, fmt_9991, 0 };
00137 static cilist io___34 = { 0, 0, 0, fmt_9990, 0 };
00138 static cilist io___35 = { 0, 0, 1, 0, 0 };
00139 static cilist io___36 = { 0, 0, 0, 0, 0 };
00140 static cilist io___37 = { 0, 0, 0, 0, 0 };
00141 static cilist io___38 = { 0, 0, 0, 0, 0 };
00142 static cilist io___39 = { 0, 0, 0, 0, 0 };
00143 static cilist io___40 = { 0, 0, 0, fmt_9987, 0 };
00144 static cilist io___41 = { 0, 0, 0, fmt_9986, 0 };
00145 static cilist io___42 = { 0, 0, 0, fmt_9986, 0 };
00146 static cilist io___43 = { 0, 0, 0, fmt_9997, 0 };
00147 static cilist io___44 = { 0, 0, 0, fmt_9996, 0 };
00148 static cilist io___45 = { 0, 0, 0, fmt_9992, 0 };
00149 static cilist io___46 = { 0, 0, 0, fmt_9989, 0 };
00150 static cilist io___47 = { 0, 0, 0, fmt_9988, 0 };
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
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342 vr_dim1 = *lda;
00343 vr_offset = 1 + vr_dim1;
00344 vr -= vr_offset;
00345 vl_dim1 = *lda;
00346 vl_offset = 1 + vl_dim1;
00347 vl -= vl_offset;
00348 bi_dim1 = *lda;
00349 bi_offset = 1 + bi_dim1;
00350 bi -= bi_offset;
00351 ai_dim1 = *lda;
00352 ai_offset = 1 + ai_dim1;
00353 ai -= ai_offset;
00354 b_dim1 = *lda;
00355 b_offset = 1 + b_dim1;
00356 b -= b_offset;
00357 a_dim1 = *lda;
00358 a_offset = 1 + a_dim1;
00359 a -= a_offset;
00360 --alphar;
00361 --alphai;
00362 --beta;
00363 --lscale;
00364 --rscale;
00365 --s;
00366 --dtru;
00367 --dif;
00368 --diftru;
00369 --work;
00370 --iwork;
00371 --result;
00372 --bwork;
00373
00374
00375 *info = 0;
00376
00377 nmax = 5;
00378
00379 if (*nsize < 0) {
00380 *info = -1;
00381 } else if (*thresh < 0.) {
00382 *info = -2;
00383 } else if (*nin <= 0) {
00384 *info = -3;
00385 } else if (*nout <= 0) {
00386 *info = -4;
00387 } else if (*lda < 1 || *lda < nmax) {
00388 *info = -6;
00389 } else if (*liwork < nmax + 6) {
00390 *info = -26;
00391 }
00392
00393
00394
00395
00396
00397
00398
00399
00400 minwrk = 1;
00401 if (*info == 0 && *lwork >= 1) {
00402 minwrk = (nmax << 1) * nmax + nmax * 12 + 16;
00403 maxwrk = nmax * 6 + nmax * ilaenv_(&c__1, "DGEQRF", " ", &nmax, &c__1,
00404 &nmax, &c__0);
00405
00406 i__1 = maxwrk, i__2 = (nmax << 1) * nmax + nmax * 12 + 16;
00407 maxwrk = max(i__1,i__2);
00408 work[1] = (doublereal) maxwrk;
00409 }
00410
00411 if (*lwork < minwrk) {
00412 *info = -24;
00413 }
00414
00415 if (*info != 0) {
00416 i__1 = -(*info);
00417 xerbla_("DDRGVX", &i__1);
00418 return 0;
00419 }
00420
00421 n = 5;
00422 ulp = dlamch_("P");
00423 ulpinv = 1. / ulp;
00424 thrsh2 = *thresh * 10.;
00425 nerrs = 0;
00426 nptknt = 0;
00427 ntestt = 0;
00428
00429 if (*nsize == 0) {
00430 goto L90;
00431 }
00432
00433
00434
00435 weight[0] = sqrt(sqrt(ulp));
00436 weight[1] = .1;
00437 weight[2] = 1.;
00438 weight[3] = 1. / weight[1];
00439 weight[4] = 1. / weight[0];
00440
00441 for (iptype = 1; iptype <= 2; ++iptype) {
00442 for (iwa = 1; iwa <= 5; ++iwa) {
00443 for (iwb = 1; iwb <= 5; ++iwb) {
00444 for (iwx = 1; iwx <= 5; ++iwx) {
00445 for (iwy = 1; iwy <= 5; ++iwy) {
00446
00447
00448
00449 dlatm6_(&iptype, &c__5, &a[a_offset], lda, &b[
00450 b_offset], &vr[vr_offset], lda, &vl[vl_offset]
00451 , lda, &weight[iwa - 1], &weight[iwb - 1], &
00452 weight[iwx - 1], &weight[iwy - 1], &dtru[1], &
00453 diftru[1]);
00454
00455
00456
00457
00458
00459 dlacpy_("F", &n, &n, &a[a_offset], lda, &ai[ai_offset]
00460 , lda);
00461 dlacpy_("F", &n, &n, &b[b_offset], lda, &bi[bi_offset]
00462 , lda);
00463
00464 dggevx_("N", "V", "V", "B", &n, &ai[ai_offset], lda, &
00465 bi[bi_offset], lda, &alphar[1], &alphai[1], &
00466 beta[1], &vl[vl_offset], lda, &vr[vr_offset],
00467 lda, ilo, ihi, &lscale[1], &rscale[1], &anorm,
00468 &bnorm, &s[1], &dif[1], &work[1], lwork, &
00469 iwork[1], &bwork[1], &linfo);
00470 if (linfo != 0) {
00471 result[1] = ulpinv;
00472 io___20.ciunit = *nout;
00473 s_wsfe(&io___20);
00474 do_fio(&c__1, "DGGEVX", (ftnlen)6);
00475 do_fio(&c__1, (char *)&linfo, (ftnlen)sizeof(
00476 integer));
00477 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer))
00478 ;
00479 do_fio(&c__1, (char *)&iptype, (ftnlen)sizeof(
00480 integer));
00481 e_wsfe();
00482 goto L30;
00483 }
00484
00485
00486
00487 dlacpy_("Full", &n, &n, &ai[ai_offset], lda, &work[1],
00488 &n);
00489 dlacpy_("Full", &n, &n, &bi[bi_offset], lda, &work[n *
00490 n + 1], &n);
00491 i__1 = n << 1;
00492 abnorm = dlange_("Fro", &n, &i__1, &work[1], &n, &
00493 work[1]);
00494
00495
00496
00497 result[1] = 0.;
00498 dget52_(&c_true, &n, &a[a_offset], lda, &b[b_offset],
00499 lda, &vl[vl_offset], lda, &alphar[1], &alphai[
00500 1], &beta[1], &work[1], &result[1]);
00501 if (result[2] > *thresh) {
00502 io___22.ciunit = *nout;
00503 s_wsfe(&io___22);
00504 do_fio(&c__1, "Left", (ftnlen)4);
00505 do_fio(&c__1, "DGGEVX", (ftnlen)6);
00506 do_fio(&c__1, (char *)&result[2], (ftnlen)sizeof(
00507 doublereal));
00508 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer))
00509 ;
00510 do_fio(&c__1, (char *)&iptype, (ftnlen)sizeof(
00511 integer));
00512 do_fio(&c__1, (char *)&iwa, (ftnlen)sizeof(
00513 integer));
00514 do_fio(&c__1, (char *)&iwb, (ftnlen)sizeof(
00515 integer));
00516 do_fio(&c__1, (char *)&iwx, (ftnlen)sizeof(
00517 integer));
00518 do_fio(&c__1, (char *)&iwy, (ftnlen)sizeof(
00519 integer));
00520 e_wsfe();
00521 }
00522
00523 result[2] = 0.;
00524 dget52_(&c_false, &n, &a[a_offset], lda, &b[b_offset],
00525 lda, &vr[vr_offset], lda, &alphar[1], &
00526 alphai[1], &beta[1], &work[1], &result[2]);
00527 if (result[3] > *thresh) {
00528 io___23.ciunit = *nout;
00529 s_wsfe(&io___23);
00530 do_fio(&c__1, "Right", (ftnlen)5);
00531 do_fio(&c__1, "DGGEVX", (ftnlen)6);
00532 do_fio(&c__1, (char *)&result[3], (ftnlen)sizeof(
00533 doublereal));
00534 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer))
00535 ;
00536 do_fio(&c__1, (char *)&iptype, (ftnlen)sizeof(
00537 integer));
00538 do_fio(&c__1, (char *)&iwa, (ftnlen)sizeof(
00539 integer));
00540 do_fio(&c__1, (char *)&iwb, (ftnlen)sizeof(
00541 integer));
00542 do_fio(&c__1, (char *)&iwx, (ftnlen)sizeof(
00543 integer));
00544 do_fio(&c__1, (char *)&iwy, (ftnlen)sizeof(
00545 integer));
00546 e_wsfe();
00547 }
00548
00549
00550
00551 result[3] = 0.;
00552 i__1 = n;
00553 for (i__ = 1; i__ <= i__1; ++i__) {
00554 if (s[i__] == 0.) {
00555 if (dtru[i__] > abnorm * ulp) {
00556 result[3] = ulpinv;
00557 }
00558 } else if (dtru[i__] == 0.) {
00559 if (s[i__] > abnorm * ulp) {
00560 result[3] = ulpinv;
00561 }
00562 } else {
00563
00564 d__3 = (d__1 = dtru[i__] / s[i__], abs(d__1)),
00565 d__4 = (d__2 = s[i__] / dtru[i__],
00566 abs(d__2));
00567 work[i__] = max(d__3,d__4);
00568
00569 d__1 = result[3], d__2 = work[i__];
00570 result[3] = max(d__1,d__2);
00571 }
00572
00573 }
00574
00575
00576
00577 result[4] = 0.;
00578 if (dif[1] == 0.) {
00579 if (diftru[1] > abnorm * ulp) {
00580 result[4] = ulpinv;
00581 }
00582 } else if (diftru[1] == 0.) {
00583 if (dif[1] > abnorm * ulp) {
00584 result[4] = ulpinv;
00585 }
00586 } else if (dif[5] == 0.) {
00587 if (diftru[5] > abnorm * ulp) {
00588 result[4] = ulpinv;
00589 }
00590 } else if (diftru[5] == 0.) {
00591 if (dif[5] > abnorm * ulp) {
00592 result[4] = ulpinv;
00593 }
00594 } else {
00595
00596 d__3 = (d__1 = diftru[1] / dif[1], abs(d__1)),
00597 d__4 = (d__2 = dif[1] / diftru[1], abs(
00598 d__2));
00599 ratio1 = max(d__3,d__4);
00600
00601 d__3 = (d__1 = diftru[5] / dif[5], abs(d__1)),
00602 d__4 = (d__2 = dif[5] / diftru[5], abs(
00603 d__2));
00604 ratio2 = max(d__3,d__4);
00605 result[4] = max(ratio1,ratio2);
00606 }
00607
00608 ntestt += 4;
00609
00610
00611
00612 for (j = 1; j <= 4; ++j) {
00613 if (result[j] >= thrsh2 && j >= 4 || result[j] >=
00614 *thresh && j <= 3) {
00615
00616
00617
00618
00619 if (nerrs == 0) {
00620 io___28.ciunit = *nout;
00621 s_wsfe(&io___28);
00622 do_fio(&c__1, "DXV", (ftnlen)3);
00623 e_wsfe();
00624
00625
00626
00627
00628
00629 io___29.ciunit = *nout;
00630 s_wsfe(&io___29);
00631 e_wsfe();
00632 io___30.ciunit = *nout;
00633 s_wsfe(&io___30);
00634 e_wsfe();
00635 io___31.ciunit = *nout;
00636 s_wsfe(&io___31);
00637 e_wsfe();
00638
00639
00640
00641 io___32.ciunit = *nout;
00642 s_wsfe(&io___32);
00643 do_fio(&c__1, "'", (ftnlen)1);
00644 do_fio(&c__1, "transpose", (ftnlen)9);
00645 do_fio(&c__1, "'", (ftnlen)1);
00646 e_wsfe();
00647
00648 }
00649 ++nerrs;
00650 if (result[j] < 1e4) {
00651 io___33.ciunit = *nout;
00652 s_wsfe(&io___33);
00653 do_fio(&c__1, (char *)&iptype, (ftnlen)
00654 sizeof(integer));
00655 do_fio(&c__1, (char *)&iwa, (ftnlen)
00656 sizeof(integer));
00657 do_fio(&c__1, (char *)&iwb, (ftnlen)
00658 sizeof(integer));
00659 do_fio(&c__1, (char *)&iwx, (ftnlen)
00660 sizeof(integer));
00661 do_fio(&c__1, (char *)&iwy, (ftnlen)
00662 sizeof(integer));
00663 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(
00664 integer));
00665 do_fio(&c__1, (char *)&result[j], (ftnlen)
00666 sizeof(doublereal));
00667 e_wsfe();
00668 } else {
00669 io___34.ciunit = *nout;
00670 s_wsfe(&io___34);
00671 do_fio(&c__1, (char *)&iptype, (ftnlen)
00672 sizeof(integer));
00673 do_fio(&c__1, (char *)&iwa, (ftnlen)
00674 sizeof(integer));
00675 do_fio(&c__1, (char *)&iwb, (ftnlen)
00676 sizeof(integer));
00677 do_fio(&c__1, (char *)&iwx, (ftnlen)
00678 sizeof(integer));
00679 do_fio(&c__1, (char *)&iwy, (ftnlen)
00680 sizeof(integer));
00681 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(
00682 integer));
00683 do_fio(&c__1, (char *)&result[j], (ftnlen)
00684 sizeof(doublereal));
00685 e_wsfe();
00686 }
00687 }
00688
00689 }
00690
00691 L30:
00692
00693
00694 ;
00695 }
00696
00697 }
00698
00699 }
00700
00701 }
00702
00703 }
00704
00705 goto L150;
00706
00707 L90:
00708
00709
00710
00711
00712 io___35.ciunit = *nin;
00713 i__1 = s_rsle(&io___35);
00714 if (i__1 != 0) {
00715 goto L150;
00716 }
00717 i__1 = do_lio(&c__3, &c__1, (char *)&n, (ftnlen)sizeof(integer));
00718 if (i__1 != 0) {
00719 goto L150;
00720 }
00721 i__1 = e_rsle();
00722 if (i__1 != 0) {
00723 goto L150;
00724 }
00725 if (n == 0) {
00726 goto L150;
00727 }
00728 i__1 = n;
00729 for (i__ = 1; i__ <= i__1; ++i__) {
00730 io___36.ciunit = *nin;
00731 s_rsle(&io___36);
00732 i__2 = n;
00733 for (j = 1; j <= i__2; ++j) {
00734 do_lio(&c__5, &c__1, (char *)&a[i__ + j * a_dim1], (ftnlen)sizeof(
00735 doublereal));
00736 }
00737 e_rsle();
00738
00739 }
00740 i__1 = n;
00741 for (i__ = 1; i__ <= i__1; ++i__) {
00742 io___37.ciunit = *nin;
00743 s_rsle(&io___37);
00744 i__2 = n;
00745 for (j = 1; j <= i__2; ++j) {
00746 do_lio(&c__5, &c__1, (char *)&b[i__ + j * b_dim1], (ftnlen)sizeof(
00747 doublereal));
00748 }
00749 e_rsle();
00750
00751 }
00752 io___38.ciunit = *nin;
00753 s_rsle(&io___38);
00754 i__1 = n;
00755 for (i__ = 1; i__ <= i__1; ++i__) {
00756 do_lio(&c__5, &c__1, (char *)&dtru[i__], (ftnlen)sizeof(doublereal));
00757 }
00758 e_rsle();
00759 io___39.ciunit = *nin;
00760 s_rsle(&io___39);
00761 i__1 = n;
00762 for (i__ = 1; i__ <= i__1; ++i__) {
00763 do_lio(&c__5, &c__1, (char *)&diftru[i__], (ftnlen)sizeof(doublereal))
00764 ;
00765 }
00766 e_rsle();
00767
00768 ++nptknt;
00769
00770
00771
00772
00773
00774 dlacpy_("F", &n, &n, &a[a_offset], lda, &ai[ai_offset], lda);
00775 dlacpy_("F", &n, &n, &b[b_offset], lda, &bi[bi_offset], lda);
00776
00777 dggevx_("N", "V", "V", "B", &n, &ai[ai_offset], lda, &bi[bi_offset], lda,
00778 &alphar[1], &alphai[1], &beta[1], &vl[vl_offset], lda, &vr[
00779 vr_offset], lda, ilo, ihi, &lscale[1], &rscale[1], &anorm, &bnorm,
00780 &s[1], &dif[1], &work[1], lwork, &iwork[1], &bwork[1], &linfo);
00781
00782 if (linfo != 0) {
00783 result[1] = ulpinv;
00784 io___40.ciunit = *nout;
00785 s_wsfe(&io___40);
00786 do_fio(&c__1, "DGGEVX", (ftnlen)6);
00787 do_fio(&c__1, (char *)&linfo, (ftnlen)sizeof(integer));
00788 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00789 do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
00790 e_wsfe();
00791 goto L140;
00792 }
00793
00794
00795
00796 dlacpy_("Full", &n, &n, &ai[ai_offset], lda, &work[1], &n);
00797 dlacpy_("Full", &n, &n, &bi[bi_offset], lda, &work[n * n + 1], &n);
00798 i__1 = n << 1;
00799 abnorm = dlange_("Fro", &n, &i__1, &work[1], &n, &work[1]);
00800
00801
00802
00803 result[1] = 0.;
00804 dget52_(&c_true, &n, &a[a_offset], lda, &b[b_offset], lda, &vl[vl_offset],
00805 lda, &alphar[1], &alphai[1], &beta[1], &work[1], &result[1]);
00806 if (result[2] > *thresh) {
00807 io___41.ciunit = *nout;
00808 s_wsfe(&io___41);
00809 do_fio(&c__1, "Left", (ftnlen)4);
00810 do_fio(&c__1, "DGGEVX", (ftnlen)6);
00811 do_fio(&c__1, (char *)&result[2], (ftnlen)sizeof(doublereal));
00812 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00813 do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
00814 e_wsfe();
00815 }
00816
00817 result[2] = 0.;
00818 dget52_(&c_false, &n, &a[a_offset], lda, &b[b_offset], lda, &vr[vr_offset]
00819 , lda, &alphar[1], &alphai[1], &beta[1], &work[1], &result[2]);
00820 if (result[3] > *thresh) {
00821 io___42.ciunit = *nout;
00822 s_wsfe(&io___42);
00823 do_fio(&c__1, "Right", (ftnlen)5);
00824 do_fio(&c__1, "DGGEVX", (ftnlen)6);
00825 do_fio(&c__1, (char *)&result[3], (ftnlen)sizeof(doublereal));
00826 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00827 do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
00828 e_wsfe();
00829 }
00830
00831
00832
00833 result[3] = 0.;
00834 i__1 = n;
00835 for (i__ = 1; i__ <= i__1; ++i__) {
00836 if (s[i__] == 0.) {
00837 if (dtru[i__] > abnorm * ulp) {
00838 result[3] = ulpinv;
00839 }
00840 } else if (dtru[i__] == 0.) {
00841 if (s[i__] > abnorm * ulp) {
00842 result[3] = ulpinv;
00843 }
00844 } else {
00845
00846 d__3 = (d__1 = dtru[i__] / s[i__], abs(d__1)), d__4 = (d__2 = s[
00847 i__] / dtru[i__], abs(d__2));
00848 work[i__] = max(d__3,d__4);
00849
00850 d__1 = result[3], d__2 = work[i__];
00851 result[3] = max(d__1,d__2);
00852 }
00853
00854 }
00855
00856
00857
00858 result[4] = 0.;
00859 if (dif[1] == 0.) {
00860 if (diftru[1] > abnorm * ulp) {
00861 result[4] = ulpinv;
00862 }
00863 } else if (diftru[1] == 0.) {
00864 if (dif[1] > abnorm * ulp) {
00865 result[4] = ulpinv;
00866 }
00867 } else if (dif[5] == 0.) {
00868 if (diftru[5] > abnorm * ulp) {
00869 result[4] = ulpinv;
00870 }
00871 } else if (diftru[5] == 0.) {
00872 if (dif[5] > abnorm * ulp) {
00873 result[4] = ulpinv;
00874 }
00875 } else {
00876
00877 d__3 = (d__1 = diftru[1] / dif[1], abs(d__1)), d__4 = (d__2 = dif[1] /
00878 diftru[1], abs(d__2));
00879 ratio1 = max(d__3,d__4);
00880
00881 d__3 = (d__1 = diftru[5] / dif[5], abs(d__1)), d__4 = (d__2 = dif[5] /
00882 diftru[5], abs(d__2));
00883 ratio2 = max(d__3,d__4);
00884 result[4] = max(ratio1,ratio2);
00885 }
00886
00887 ntestt += 4;
00888
00889
00890
00891 for (j = 1; j <= 4; ++j) {
00892 if (result[j] >= thrsh2) {
00893
00894
00895
00896
00897 if (nerrs == 0) {
00898 io___43.ciunit = *nout;
00899 s_wsfe(&io___43);
00900 do_fio(&c__1, "DXV", (ftnlen)3);
00901 e_wsfe();
00902
00903
00904
00905
00906
00907 io___44.ciunit = *nout;
00908 s_wsfe(&io___44);
00909 e_wsfe();
00910
00911
00912
00913 io___45.ciunit = *nout;
00914 s_wsfe(&io___45);
00915 do_fio(&c__1, "'", (ftnlen)1);
00916 do_fio(&c__1, "transpose", (ftnlen)9);
00917 do_fio(&c__1, "'", (ftnlen)1);
00918 e_wsfe();
00919
00920 }
00921 ++nerrs;
00922 if (result[j] < 1e4) {
00923 io___46.ciunit = *nout;
00924 s_wsfe(&io___46);
00925 do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
00926 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00927 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
00928 do_fio(&c__1, (char *)&result[j], (ftnlen)sizeof(doublereal));
00929 e_wsfe();
00930 } else {
00931 io___47.ciunit = *nout;
00932 s_wsfe(&io___47);
00933 do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
00934 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00935 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
00936 do_fio(&c__1, (char *)&result[j], (ftnlen)sizeof(doublereal));
00937 e_wsfe();
00938 }
00939 }
00940
00941 }
00942
00943 L140:
00944
00945 goto L90;
00946 L150:
00947
00948
00949
00950 alasvm_("DXV", nout, &nerrs, &ntestt, &c__0);
00951
00952 work[1] = (doublereal) maxwrk;
00953
00954 return 0;
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969 }