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 m, n, mplusn, k;
00020 logical fs;
00021 } mn_;
00022
00023 #define mn_1 mn_
00024
00025
00026
00027 static integer c__1 = 1;
00028 static integer c__0 = 0;
00029 static integer c_n1 = -1;
00030 static real c_b26 = 0.f;
00031 static integer c__3 = 3;
00032 static integer c__4 = 4;
00033
00034 int sdrgsx_(integer *nsize, integer *ncmax, real *thresh,
00035 integer *nin, integer *nout, real *a, integer *lda, real *b, real *ai,
00036 real *bi, real *z__, real *q, real *alphar, real *alphai, real *beta,
00037 real *c__, integer *ldc, real *s, real *work, integer *lwork,
00038 integer *iwork, integer *liwork, logical *bwork, integer *info)
00039 {
00040
00041 static char fmt_9999[] = "(\002 SDRGSX: \002,a,\002 returned INFO=\002,i"
00042 "6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002)\002)";
00043 static char fmt_9997[] = "(\002 SDRGSX: SGET53 returned INFO=\002,i1,"
00044 "\002 for eigenvalue \002,i6,\002.\002,/9x,\002N=\002,i6,\002, JT"
00045 "YPE=\002,i6,\002)\002)";
00046 static char fmt_9996[] = "(\002 SDRGSX: S not in Schur form at eigenvalu"
00047 "e \002,i6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002"
00048 ")\002)";
00049 static char fmt_9995[] = "(/1x,a3,\002 -- Real Expert Generalized Schur "
00050 "form\002,\002 problem driver\002)";
00051 static char fmt_9993[] = "(\002 Matrix types: \002,/\002 1: A is a blo"
00052 "ck diagonal matrix of Jordan blocks \002,\002and B is the identi"
00053 "ty \002,/\002 matrix, \002,/\002 2: A and B are upper tri"
00054 "angular matrices, \002,/\002 3: A and B are as type 2, but eac"
00055 "h second diagonal \002,\002block in A_11 and \002,/\002 eac"
00056 "h third diaongal block in A_22 are 2x2 blocks,\002,/\002 4: A "
00057 "and B are block diagonal matrices, \002,/\002 5: (A,B) has pot"
00058 "entially close or common \002,\002eigenvalues.\002,/)";
00059 static char fmt_9992[] = "(/\002 Tests performed: (S is Schur, T is tri"
00060 "angular, \002,\002Q and Z are \002,a,\002,\002,/19x,\002 a is al"
00061 "pha, b is beta, and \002,a,\002 means \002,a,\002.)\002,/\002 1"
00062 " = | A - Q S Z\002,a,\002 | / ( |A| n ulp ) 2 = | B - Q T "
00063 "Z\002,a,\002 | / ( |B| n ulp )\002,/\002 3 = | I - QQ\002,a,"
00064 "\002 | / ( n ulp ) 4 = | I - ZZ\002,a,\002 | / ( n u"
00065 "lp )\002,/\002 5 = 1/ULP if A is not in \002,\002Schur form "
00066 "S\002,/\002 6 = difference between (alpha,beta)\002,\002 and di"
00067 "agonals of (S,T)\002,/\002 7 = 1/ULP if SDIM is not the correc"
00068 "t number of \002,\002selected eigenvalues\002,/\002 8 = 1/ULP "
00069 "if DIFEST/DIFTRU > 10*THRESH or \002,\002DIFTRU/DIFEST > 10*THRE"
00070 "SH\002,/\002 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B"
00071 ") \002,\002when reordering fails\002,/\002 10 = 1/ULP if PLEST/"
00072 "PLTRU > THRESH or \002,\002PLTRU/PLEST > THRESH\002,/\002 ( T"
00073 "est 10 is only for input examples )\002,/)";
00074 static char fmt_9991[] = "(\002 Matrix order=\002,i2,\002, type=\002,i2"
00075 ",\002, a=\002,e10.4,\002, order(A_11)=\002,i2,\002, result \002,"
00076 "i2,\002 is \002,0p,f8.2)";
00077 static char fmt_9990[] = "(\002 Matrix order=\002,i2,\002, type=\002,i2"
00078 ",\002, a=\002,e10.4,\002, order(A_11)=\002,i2,\002, result \002,"
00079 "i2,\002 is \002,0p,e10.4)";
00080 static char fmt_9998[] = "(\002 SDRGSX: \002,a,\002 returned INFO=\002,i"
00081 "6,\002.\002,/9x,\002N=\002,i6,\002, Input Example #\002,i2,\002"
00082 ")\002)";
00083 static char fmt_9994[] = "(\002Input Example\002)";
00084 static char fmt_9989[] = "(\002 Input example #\002,i2,\002, matrix orde"
00085 "r=\002,i4,\002,\002,\002 result \002,i2,\002 is\002,0p,f8.2)";
00086 static char fmt_9988[] = "(\002 Input example #\002,i2,\002, matrix orde"
00087 "r=\002,i4,\002,\002,\002 result \002,i2,\002 is\002,1p,e10.3)";
00088
00089
00090 integer a_dim1, a_offset, ai_dim1, ai_offset, b_dim1, b_offset, bi_dim1,
00091 bi_offset, c_dim1, c_offset, q_dim1, q_offset, z_dim1, z_offset,
00092 i__1, i__2, i__3, i__4;
00093 real r__1, r__2, r__3, r__4, r__5, r__6, r__7, r__8, r__9, r__10;
00094
00095
00096 double sqrt(doublereal);
00097 integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void),
00098 s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
00099 e_rsle(void);
00100
00101
00102 integer i__, j, i1, mm;
00103 real pl[2];
00104 integer mn2, qba, qbb;
00105 real ulp, temp1, temp2, abnrm;
00106 integer ifunc, iinfo, linfo;
00107 extern int sget51_(integer *, integer *, real *, integer
00108 *, real *, integer *, real *, integer *, real *, integer *, real *
00109 , real *), sget53_(real *, integer *, real *, integer *, real *,
00110 real *, real *, real *, integer *);
00111 char sense[1];
00112 integer nerrs, ntest;
00113 real pltru;
00114 extern int slakf2_(integer *, integer *, real *, integer
00115 *, real *, real *, real *, real *, integer *), slatm5_(integer *,
00116 integer *, integer *, real *, integer *, real *, integer *, real *
00117 , integer *, real *, integer *, real *, integer *, real *,
00118 integer *, real *, integer *, real *, integer *, real *, integer *
00119 , integer *);
00120 real thrsh2;
00121 logical ilabad;
00122 extern int slabad_(real *, real *);
00123 integer bdspac;
00124 extern doublereal slamch_(char *), slange_(char *, integer *,
00125 integer *, real *, integer *, real *);
00126 extern int xerbla_(char *, integer *);
00127 real difest[2];
00128 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
00129 integer *, integer *);
00130 real bignum;
00131 extern int alasvm_(char *, integer *, integer *, integer
00132 *, integer *);
00133 real weight;
00134 extern int sgesvd_(char *, char *, integer *, integer *,
00135 real *, integer *, real *, real *, integer *, real *, integer *,
00136 real *, integer *, integer *), slacpy_(char *,
00137 integer *, integer *, real *, integer *, real *, integer *);
00138 real diftru;
00139 extern int slaset_(char *, integer *, integer *, real *,
00140 real *, real *, integer *), sggesx_(char *, char *, char *
00141 , L_fp, char *, integer *, real *, integer *, real *, integer *,
00142 integer *, real *, real *, real *, real *, integer *, real *,
00143 integer *, real *, real *, real *, integer *, integer *, integer *
00144 , logical *, integer *);
00145 integer minwrk, maxwrk;
00146 real smlnum, ulpinv;
00147 integer nptknt;
00148 real result[10];
00149 extern logical slctsx_();
00150 integer ntestt, prtype;
00151
00152
00153 static cilist io___22 = { 0, 0, 0, fmt_9999, 0 };
00154 static cilist io___31 = { 0, 0, 0, fmt_9997, 0 };
00155 static cilist io___32 = { 0, 0, 0, fmt_9996, 0 };
00156 static cilist io___35 = { 0, 0, 0, fmt_9995, 0 };
00157 static cilist io___36 = { 0, 0, 0, fmt_9993, 0 };
00158 static cilist io___37 = { 0, 0, 0, fmt_9992, 0 };
00159 static cilist io___39 = { 0, 0, 0, fmt_9991, 0 };
00160 static cilist io___40 = { 0, 0, 0, fmt_9990, 0 };
00161 static cilist io___42 = { 0, 0, 1, 0, 0 };
00162 static cilist io___43 = { 0, 0, 1, 0, 0 };
00163 static cilist io___44 = { 0, 0, 0, 0, 0 };
00164 static cilist io___45 = { 0, 0, 0, 0, 0 };
00165 static cilist io___46 = { 0, 0, 0, 0, 0 };
00166 static cilist io___48 = { 0, 0, 0, fmt_9998, 0 };
00167 static cilist io___49 = { 0, 0, 0, fmt_9997, 0 };
00168 static cilist io___50 = { 0, 0, 0, fmt_9996, 0 };
00169 static cilist io___51 = { 0, 0, 0, fmt_9995, 0 };
00170 static cilist io___52 = { 0, 0, 0, fmt_9994, 0 };
00171 static cilist io___53 = { 0, 0, 0, fmt_9992, 0 };
00172 static cilist io___54 = { 0, 0, 0, fmt_9989, 0 };
00173 static cilist io___55 = { 0, 0, 0, fmt_9988, 0 };
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
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446 q_dim1 = *lda;
00447 q_offset = 1 + q_dim1;
00448 q -= q_offset;
00449 z_dim1 = *lda;
00450 z_offset = 1 + z_dim1;
00451 z__ -= z_offset;
00452 bi_dim1 = *lda;
00453 bi_offset = 1 + bi_dim1;
00454 bi -= bi_offset;
00455 ai_dim1 = *lda;
00456 ai_offset = 1 + ai_dim1;
00457 ai -= ai_offset;
00458 b_dim1 = *lda;
00459 b_offset = 1 + b_dim1;
00460 b -= b_offset;
00461 a_dim1 = *lda;
00462 a_offset = 1 + a_dim1;
00463 a -= a_offset;
00464 --alphar;
00465 --alphai;
00466 --beta;
00467 c_dim1 = *ldc;
00468 c_offset = 1 + c_dim1;
00469 c__ -= c_offset;
00470 --s;
00471 --work;
00472 --iwork;
00473 --bwork;
00474
00475
00476 if (*nsize < 0) {
00477 *info = -1;
00478 } else if (*thresh < 0.f) {
00479 *info = -2;
00480 } else if (*nin <= 0) {
00481 *info = -3;
00482 } else if (*nout <= 0) {
00483 *info = -4;
00484 } else if (*lda < 1 || *lda < *nsize) {
00485 *info = -6;
00486 } else if (*ldc < 1 || *ldc < *nsize * *nsize / 2) {
00487 *info = -17;
00488 } else if (*liwork < *nsize + 6) {
00489 *info = -21;
00490 }
00491
00492
00493
00494
00495
00496
00497
00498
00499 minwrk = 1;
00500 if (*info == 0 && *lwork >= 1) {
00501
00502
00503 i__1 = (*nsize + 1) * 10, i__2 = *nsize * 5 * *nsize / 2;
00504 minwrk = max(i__1,i__2);
00505
00506
00507
00508 maxwrk = (*nsize + 1) * 9 + *nsize * ilaenv_(&c__1, "SGEQRF", " ",
00509 nsize, &c__1, nsize, &c__0);
00510
00511 i__1 = maxwrk, i__2 = (*nsize + 1) * 9 + *nsize * ilaenv_(&c__1,
00512 "SORGQR", " ", nsize, &c__1, nsize, &c_n1);
00513 maxwrk = max(i__1,i__2);
00514
00515
00516
00517 bdspac = *nsize * 5 * *nsize / 2;
00518
00519 i__3 = *nsize * *nsize / 2;
00520 i__4 = *nsize * *nsize / 2;
00521 i__1 = maxwrk, i__2 = *nsize * 3 * *nsize / 2 + *nsize * *nsize *
00522 ilaenv_(&c__1, "SGEBRD", " ", &i__3, &i__4, &c_n1, &c_n1);
00523 maxwrk = max(i__1,i__2);
00524 maxwrk = max(maxwrk,bdspac);
00525
00526 maxwrk = max(maxwrk,minwrk);
00527
00528 work[1] = (real) maxwrk;
00529 }
00530
00531 if (*lwork < minwrk) {
00532 *info = -19;
00533 }
00534
00535 if (*info != 0) {
00536 i__1 = -(*info);
00537 xerbla_("SDRGSX", &i__1);
00538 return 0;
00539 }
00540
00541
00542
00543 ulp = slamch_("P");
00544 ulpinv = 1.f / ulp;
00545 smlnum = slamch_("S") / ulp;
00546 bignum = 1.f / smlnum;
00547 slabad_(&smlnum, &bignum);
00548 thrsh2 = *thresh * 10.f;
00549 ntestt = 0;
00550 nerrs = 0;
00551
00552
00553
00554 ifunc = 0;
00555 if (*nsize == 0) {
00556 goto L70;
00557 }
00558
00559
00560
00561
00562
00563 prtype = 0;
00564 qba = 3;
00565 qbb = 4;
00566 weight = sqrt(ulp);
00567
00568 for (ifunc = 0; ifunc <= 3; ++ifunc) {
00569 for (prtype = 1; prtype <= 5; ++prtype) {
00570 i__1 = *nsize - 1;
00571 for (mn_1.m = 1; mn_1.m <= i__1; ++mn_1.m) {
00572 i__2 = *nsize - mn_1.m;
00573 for (mn_1.n = 1; mn_1.n <= i__2; ++mn_1.n) {
00574
00575 weight = 1.f / weight;
00576 mn_1.mplusn = mn_1.m + mn_1.n;
00577
00578
00579
00580 mn_1.fs = TRUE_;
00581 mn_1.k = 0;
00582
00583 slaset_("Full", &mn_1.mplusn, &mn_1.mplusn, &c_b26, &
00584 c_b26, &ai[ai_offset], lda);
00585 slaset_("Full", &mn_1.mplusn, &mn_1.mplusn, &c_b26, &
00586 c_b26, &bi[bi_offset], lda);
00587
00588 slatm5_(&prtype, &mn_1.m, &mn_1.n, &ai[ai_offset], lda, &
00589 ai[mn_1.m + 1 + (mn_1.m + 1) * ai_dim1], lda, &ai[
00590 (mn_1.m + 1) * ai_dim1 + 1], lda, &bi[bi_offset],
00591 lda, &bi[mn_1.m + 1 + (mn_1.m + 1) * bi_dim1],
00592 lda, &bi[(mn_1.m + 1) * bi_dim1 + 1], lda, &q[
00593 q_offset], lda, &z__[z_offset], lda, &weight, &
00594 qba, &qbb);
00595
00596
00597
00598
00599
00600
00601 if (ifunc == 0) {
00602 *(unsigned char *)sense = 'N';
00603 } else if (ifunc == 1) {
00604 *(unsigned char *)sense = 'E';
00605 } else if (ifunc == 2) {
00606 *(unsigned char *)sense = 'V';
00607 } else if (ifunc == 3) {
00608 *(unsigned char *)sense = 'B';
00609 }
00610
00611 slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &ai[ai_offset]
00612 , lda, &a[a_offset], lda);
00613 slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &bi[bi_offset]
00614 , lda, &b[b_offset], lda);
00615
00616 sggesx_("V", "V", "S", (L_fp)slctsx_, sense, &mn_1.mplusn,
00617 &ai[ai_offset], lda, &bi[bi_offset], lda, &mm, &
00618 alphar[1], &alphai[1], &beta[1], &q[q_offset],
00619 lda, &z__[z_offset], lda, pl, difest, &work[1],
00620 lwork, &iwork[1], liwork, &bwork[1], &linfo);
00621
00622 if (linfo != 0 && linfo != mn_1.mplusn + 2) {
00623 result[0] = ulpinv;
00624 io___22.ciunit = *nout;
00625 s_wsfe(&io___22);
00626 do_fio(&c__1, "SGGESX", (ftnlen)6);
00627 do_fio(&c__1, (char *)&linfo, (ftnlen)sizeof(integer))
00628 ;
00629 do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)sizeof(
00630 integer));
00631 do_fio(&c__1, (char *)&prtype, (ftnlen)sizeof(integer)
00632 );
00633 e_wsfe();
00634 *info = linfo;
00635 goto L30;
00636 }
00637
00638
00639
00640 slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &ai[ai_offset]
00641 , lda, &work[1], &mn_1.mplusn);
00642 slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &bi[bi_offset]
00643 , lda, &work[mn_1.mplusn * mn_1.mplusn + 1], &
00644 mn_1.mplusn);
00645 i__3 = mn_1.mplusn << 1;
00646 abnrm = slange_("Fro", &mn_1.mplusn, &i__3, &work[1], &
00647 mn_1.mplusn, &work[1]);
00648
00649
00650
00651 sget51_(&c__1, &mn_1.mplusn, &a[a_offset], lda, &ai[
00652 ai_offset], lda, &q[q_offset], lda, &z__[z_offset]
00653 , lda, &work[1], result);
00654 sget51_(&c__1, &mn_1.mplusn, &b[b_offset], lda, &bi[
00655 bi_offset], lda, &q[q_offset], lda, &z__[z_offset]
00656 , lda, &work[1], &result[1]);
00657 sget51_(&c__3, &mn_1.mplusn, &b[b_offset], lda, &bi[
00658 bi_offset], lda, &q[q_offset], lda, &q[q_offset],
00659 lda, &work[1], &result[2]);
00660 sget51_(&c__3, &mn_1.mplusn, &b[b_offset], lda, &bi[
00661 bi_offset], lda, &z__[z_offset], lda, &z__[
00662 z_offset], lda, &work[1], &result[3]);
00663 ntest = 4;
00664
00665
00666
00667
00668 temp1 = 0.f;
00669 result[4] = 0.f;
00670 result[5] = 0.f;
00671
00672 i__3 = mn_1.mplusn;
00673 for (j = 1; j <= i__3; ++j) {
00674 ilabad = FALSE_;
00675 if (alphai[j] == 0.f) {
00676
00677 r__7 = smlnum, r__8 = (r__2 = alphar[j], dabs(
00678 r__2)), r__7 = max(r__7,r__8), r__8 = (
00679 r__3 = ai[j + j * ai_dim1], dabs(r__3));
00680
00681 r__9 = smlnum, r__10 = (r__5 = beta[j], dabs(r__5)
00682 ), r__9 = max(r__9,r__10), r__10 = (r__6 =
00683 bi[j + j * bi_dim1], dabs(r__6));
00684 temp2 = ((r__1 = alphar[j] - ai[j + j * ai_dim1],
00685 dabs(r__1)) / dmax(r__7,r__8) + (r__4 =
00686 beta[j] - bi[j + j * bi_dim1], dabs(r__4))
00687 / dmax(r__9,r__10)) / ulp;
00688 if (j < mn_1.mplusn) {
00689 if (ai[j + 1 + j * ai_dim1] != 0.f) {
00690 ilabad = TRUE_;
00691 result[4] = ulpinv;
00692 }
00693 }
00694 if (j > 1) {
00695 if (ai[j + (j - 1) * ai_dim1] != 0.f) {
00696 ilabad = TRUE_;
00697 result[4] = ulpinv;
00698 }
00699 }
00700 } else {
00701 if (alphai[j] > 0.f) {
00702 i1 = j;
00703 } else {
00704 i1 = j - 1;
00705 }
00706 if (i1 <= 0 || i1 >= mn_1.mplusn) {
00707 ilabad = TRUE_;
00708 } else if (i1 < mn_1.mplusn - 1) {
00709 if (ai[i1 + 2 + (i1 + 1) * ai_dim1] != 0.f) {
00710 ilabad = TRUE_;
00711 result[4] = ulpinv;
00712 }
00713 } else if (i1 > 1) {
00714 if (ai[i1 + (i1 - 1) * ai_dim1] != 0.f) {
00715 ilabad = TRUE_;
00716 result[4] = ulpinv;
00717 }
00718 }
00719 if (! ilabad) {
00720 sget53_(&ai[i1 + i1 * ai_dim1], lda, &bi[i1 +
00721 i1 * bi_dim1], lda, &beta[j], &alphar[
00722 j], &alphai[j], &temp2, &iinfo);
00723 if (iinfo >= 3) {
00724 io___31.ciunit = *nout;
00725 s_wsfe(&io___31);
00726 do_fio(&c__1, (char *)&iinfo, (ftnlen)
00727 sizeof(integer));
00728 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(
00729 integer));
00730 do_fio(&c__1, (char *)&mn_1.mplusn, (
00731 ftnlen)sizeof(integer));
00732 do_fio(&c__1, (char *)&prtype, (ftnlen)
00733 sizeof(integer));
00734 e_wsfe();
00735 *info = abs(iinfo);
00736 }
00737 } else {
00738 temp2 = ulpinv;
00739 }
00740 }
00741 temp1 = dmax(temp1,temp2);
00742 if (ilabad) {
00743 io___32.ciunit = *nout;
00744 s_wsfe(&io___32);
00745 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer))
00746 ;
00747 do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)
00748 sizeof(integer));
00749 do_fio(&c__1, (char *)&prtype, (ftnlen)sizeof(
00750 integer));
00751 e_wsfe();
00752 }
00753
00754 }
00755 result[5] = temp1;
00756 ntest += 2;
00757
00758
00759
00760 result[6] = 0.f;
00761 if (linfo == mn_1.mplusn + 3) {
00762 result[6] = ulpinv;
00763 } else if (mm != mn_1.n) {
00764 result[6] = ulpinv;
00765 }
00766 ++ntest;
00767
00768
00769
00770
00771 result[7] = 0.f;
00772 mn2 = mm * (mn_1.mplusn - mm) << 1;
00773 if (ifunc >= 2 && mn2 <= *ncmax * *ncmax) {
00774
00775
00776
00777
00778 i__3 = mn_1.mplusn - mm;
00779 slakf2_(&mm, &i__3, &ai[ai_offset], lda, &ai[mm + 1 +
00780 (mm + 1) * ai_dim1], &bi[bi_offset], &bi[mm +
00781 1 + (mm + 1) * bi_dim1], &c__[c_offset], ldc);
00782
00783 i__3 = *lwork - 2;
00784 sgesvd_("N", "N", &mn2, &mn2, &c__[c_offset], ldc, &s[
00785 1], &work[1], &c__1, &work[2], &c__1, &work[3]
00786 , &i__3, info);
00787 diftru = s[mn2];
00788
00789 if (difest[1] == 0.f) {
00790 if (diftru > abnrm * ulp) {
00791 result[7] = ulpinv;
00792 }
00793 } else if (diftru == 0.f) {
00794 if (difest[1] > abnrm * ulp) {
00795 result[7] = ulpinv;
00796 }
00797 } else if (diftru > thrsh2 * difest[1] || diftru *
00798 thrsh2 < difest[1]) {
00799
00800 r__1 = diftru / difest[1], r__2 = difest[1] /
00801 diftru;
00802 result[7] = dmax(r__1,r__2);
00803 }
00804 ++ntest;
00805 }
00806
00807
00808
00809 result[8] = 0.f;
00810 if (linfo == mn_1.mplusn + 2) {
00811 if (diftru > abnrm * ulp) {
00812 result[8] = ulpinv;
00813 }
00814 if (ifunc > 1 && difest[1] != 0.f) {
00815 result[8] = ulpinv;
00816 }
00817 if (ifunc == 1 && pl[0] != 0.f) {
00818 result[8] = ulpinv;
00819 }
00820 ++ntest;
00821 }
00822
00823 ntestt += ntest;
00824
00825
00826
00827 for (j = 1; j <= 9; ++j) {
00828 if (result[j - 1] >= *thresh) {
00829
00830
00831
00832
00833 if (nerrs == 0) {
00834 io___35.ciunit = *nout;
00835 s_wsfe(&io___35);
00836 do_fio(&c__1, "SGX", (ftnlen)3);
00837 e_wsfe();
00838
00839
00840
00841 io___36.ciunit = *nout;
00842 s_wsfe(&io___36);
00843 e_wsfe();
00844
00845
00846
00847 io___37.ciunit = *nout;
00848 s_wsfe(&io___37);
00849 do_fio(&c__1, "orthogonal", (ftnlen)10);
00850 do_fio(&c__1, "'", (ftnlen)1);
00851 do_fio(&c__1, "transpose", (ftnlen)9);
00852 for (i__ = 1; i__ <= 4; ++i__) {
00853 do_fio(&c__1, "'", (ftnlen)1);
00854 }
00855 e_wsfe();
00856
00857 }
00858 ++nerrs;
00859 if (result[j - 1] < 1e4f) {
00860 io___39.ciunit = *nout;
00861 s_wsfe(&io___39);
00862 do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)
00863 sizeof(integer));
00864 do_fio(&c__1, (char *)&prtype, (ftnlen)sizeof(
00865 integer));
00866 do_fio(&c__1, (char *)&weight, (ftnlen)sizeof(
00867 real));
00868 do_fio(&c__1, (char *)&mn_1.m, (ftnlen)sizeof(
00869 integer));
00870 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(
00871 integer));
00872 do_fio(&c__1, (char *)&result[j - 1], (ftnlen)
00873 sizeof(real));
00874 e_wsfe();
00875 } else {
00876 io___40.ciunit = *nout;
00877 s_wsfe(&io___40);
00878 do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)
00879 sizeof(integer));
00880 do_fio(&c__1, (char *)&prtype, (ftnlen)sizeof(
00881 integer));
00882 do_fio(&c__1, (char *)&weight, (ftnlen)sizeof(
00883 real));
00884 do_fio(&c__1, (char *)&mn_1.m, (ftnlen)sizeof(
00885 integer));
00886 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(
00887 integer));
00888 do_fio(&c__1, (char *)&result[j - 1], (ftnlen)
00889 sizeof(real));
00890 e_wsfe();
00891 }
00892 }
00893
00894 }
00895
00896 L30:
00897 ;
00898 }
00899
00900 }
00901
00902 }
00903
00904 }
00905
00906 goto L150;
00907
00908 L70:
00909
00910
00911
00912
00913 nptknt = 0;
00914
00915 L80:
00916 io___42.ciunit = *nin;
00917 i__1 = s_rsle(&io___42);
00918 if (i__1 != 0) {
00919 goto L140;
00920 }
00921 i__1 = do_lio(&c__3, &c__1, (char *)&mn_1.mplusn, (ftnlen)sizeof(integer))
00922 ;
00923 if (i__1 != 0) {
00924 goto L140;
00925 }
00926 i__1 = e_rsle();
00927 if (i__1 != 0) {
00928 goto L140;
00929 }
00930 if (mn_1.mplusn == 0) {
00931 goto L140;
00932 }
00933 io___43.ciunit = *nin;
00934 i__1 = s_rsle(&io___43);
00935 if (i__1 != 0) {
00936 goto L140;
00937 }
00938 i__1 = do_lio(&c__3, &c__1, (char *)&mn_1.n, (ftnlen)sizeof(integer));
00939 if (i__1 != 0) {
00940 goto L140;
00941 }
00942 i__1 = e_rsle();
00943 if (i__1 != 0) {
00944 goto L140;
00945 }
00946 i__1 = mn_1.mplusn;
00947 for (i__ = 1; i__ <= i__1; ++i__) {
00948 io___44.ciunit = *nin;
00949 s_rsle(&io___44);
00950 i__2 = mn_1.mplusn;
00951 for (j = 1; j <= i__2; ++j) {
00952 do_lio(&c__4, &c__1, (char *)&ai[i__ + j * ai_dim1], (ftnlen)
00953 sizeof(real));
00954 }
00955 e_rsle();
00956
00957 }
00958 i__1 = mn_1.mplusn;
00959 for (i__ = 1; i__ <= i__1; ++i__) {
00960 io___45.ciunit = *nin;
00961 s_rsle(&io___45);
00962 i__2 = mn_1.mplusn;
00963 for (j = 1; j <= i__2; ++j) {
00964 do_lio(&c__4, &c__1, (char *)&bi[i__ + j * bi_dim1], (ftnlen)
00965 sizeof(real));
00966 }
00967 e_rsle();
00968
00969 }
00970 io___46.ciunit = *nin;
00971 s_rsle(&io___46);
00972 do_lio(&c__4, &c__1, (char *)&pltru, (ftnlen)sizeof(real));
00973 do_lio(&c__4, &c__1, (char *)&diftru, (ftnlen)sizeof(real));
00974 e_rsle();
00975
00976 ++nptknt;
00977 mn_1.fs = TRUE_;
00978 mn_1.k = 0;
00979 mn_1.m = mn_1.mplusn - mn_1.n;
00980
00981 slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &ai[ai_offset], lda, &a[
00982 a_offset], lda);
00983 slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &bi[bi_offset], lda, &b[
00984 b_offset], lda);
00985
00986
00987
00988
00989 sggesx_("V", "V", "S", (L_fp)slctsx_, "B", &mn_1.mplusn, &ai[ai_offset],
00990 lda, &bi[bi_offset], lda, &mm, &alphar[1], &alphai[1], &beta[1], &
00991 q[q_offset], lda, &z__[z_offset], lda, pl, difest, &work[1],
00992 lwork, &iwork[1], liwork, &bwork[1], &linfo);
00993
00994 if (linfo != 0 && linfo != mn_1.mplusn + 2) {
00995 result[0] = ulpinv;
00996 io___48.ciunit = *nout;
00997 s_wsfe(&io___48);
00998 do_fio(&c__1, "SGGESX", (ftnlen)6);
00999 do_fio(&c__1, (char *)&linfo, (ftnlen)sizeof(integer));
01000 do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)sizeof(integer));
01001 do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
01002 e_wsfe();
01003 goto L130;
01004 }
01005
01006
01007
01008
01009 slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &ai[ai_offset], lda, &work[1],
01010 &mn_1.mplusn);
01011 slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &bi[bi_offset], lda, &work[
01012 mn_1.mplusn * mn_1.mplusn + 1], &mn_1.mplusn);
01013 i__1 = mn_1.mplusn << 1;
01014 abnrm = slange_("Fro", &mn_1.mplusn, &i__1, &work[1], &mn_1.mplusn, &work[
01015 1]);
01016
01017
01018
01019 sget51_(&c__1, &mn_1.mplusn, &a[a_offset], lda, &ai[ai_offset], lda, &q[
01020 q_offset], lda, &z__[z_offset], lda, &work[1], result);
01021 sget51_(&c__1, &mn_1.mplusn, &b[b_offset], lda, &bi[bi_offset], lda, &q[
01022 q_offset], lda, &z__[z_offset], lda, &work[1], &result[1]);
01023 sget51_(&c__3, &mn_1.mplusn, &b[b_offset], lda, &bi[bi_offset], lda, &q[
01024 q_offset], lda, &q[q_offset], lda, &work[1], &result[2]);
01025 sget51_(&c__3, &mn_1.mplusn, &b[b_offset], lda, &bi[bi_offset], lda, &z__[
01026 z_offset], lda, &z__[z_offset], lda, &work[1], &result[3]);
01027
01028
01029
01030
01031 ntest = 6;
01032 temp1 = 0.f;
01033 result[4] = 0.f;
01034 result[5] = 0.f;
01035
01036 i__1 = mn_1.mplusn;
01037 for (j = 1; j <= i__1; ++j) {
01038 ilabad = FALSE_;
01039 if (alphai[j] == 0.f) {
01040
01041 r__7 = smlnum, r__8 = (r__2 = alphar[j], dabs(r__2)), r__7 = max(
01042 r__7,r__8), r__8 = (r__3 = ai[j + j * ai_dim1], dabs(r__3)
01043 );
01044
01045 r__9 = smlnum, r__10 = (r__5 = beta[j], dabs(r__5)), r__9 = max(
01046 r__9,r__10), r__10 = (r__6 = bi[j + j * bi_dim1], dabs(
01047 r__6));
01048 temp2 = ((r__1 = alphar[j] - ai[j + j * ai_dim1], dabs(r__1)) /
01049 dmax(r__7,r__8) + (r__4 = beta[j] - bi[j + j * bi_dim1],
01050 dabs(r__4)) / dmax(r__9,r__10)) / ulp;
01051 if (j < mn_1.mplusn) {
01052 if (ai[j + 1 + j * ai_dim1] != 0.f) {
01053 ilabad = TRUE_;
01054 result[4] = ulpinv;
01055 }
01056 }
01057 if (j > 1) {
01058 if (ai[j + (j - 1) * ai_dim1] != 0.f) {
01059 ilabad = TRUE_;
01060 result[4] = ulpinv;
01061 }
01062 }
01063 } else {
01064 if (alphai[j] > 0.f) {
01065 i1 = j;
01066 } else {
01067 i1 = j - 1;
01068 }
01069 if (i1 <= 0 || i1 >= mn_1.mplusn) {
01070 ilabad = TRUE_;
01071 } else if (i1 < mn_1.mplusn - 1) {
01072 if (ai[i1 + 2 + (i1 + 1) * ai_dim1] != 0.f) {
01073 ilabad = TRUE_;
01074 result[4] = ulpinv;
01075 }
01076 } else if (i1 > 1) {
01077 if (ai[i1 + (i1 - 1) * ai_dim1] != 0.f) {
01078 ilabad = TRUE_;
01079 result[4] = ulpinv;
01080 }
01081 }
01082 if (! ilabad) {
01083 sget53_(&ai[i1 + i1 * ai_dim1], lda, &bi[i1 + i1 * bi_dim1],
01084 lda, &beta[j], &alphar[j], &alphai[j], &temp2, &iinfo)
01085 ;
01086 if (iinfo >= 3) {
01087 io___49.ciunit = *nout;
01088 s_wsfe(&io___49);
01089 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01090 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
01091 do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)sizeof(
01092 integer));
01093 do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
01094 e_wsfe();
01095 *info = abs(iinfo);
01096 }
01097 } else {
01098 temp2 = ulpinv;
01099 }
01100 }
01101 temp1 = dmax(temp1,temp2);
01102 if (ilabad) {
01103 io___50.ciunit = *nout;
01104 s_wsfe(&io___50);
01105 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
01106 do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)sizeof(integer));
01107 do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
01108 e_wsfe();
01109 }
01110
01111 }
01112 result[5] = temp1;
01113
01114
01115
01116 ntest = 7;
01117 result[6] = 0.f;
01118 if (linfo == mn_1.mplusn + 3) {
01119 result[6] = ulpinv;
01120 }
01121
01122
01123
01124 ntest = 8;
01125 result[7] = 0.f;
01126 if (difest[1] == 0.f) {
01127 if (diftru > abnrm * ulp) {
01128 result[7] = ulpinv;
01129 }
01130 } else if (diftru == 0.f) {
01131 if (difest[1] > abnrm * ulp) {
01132 result[7] = ulpinv;
01133 }
01134 } else if (diftru > thrsh2 * difest[1] || diftru * thrsh2 < difest[1]) {
01135
01136 r__1 = diftru / difest[1], r__2 = difest[1] / diftru;
01137 result[7] = dmax(r__1,r__2);
01138 }
01139
01140
01141
01142 ntest = 9;
01143 result[8] = 0.f;
01144 if (linfo == mn_1.mplusn + 2) {
01145 if (diftru > abnrm * ulp) {
01146 result[8] = ulpinv;
01147 }
01148 if (ifunc > 1 && difest[1] != 0.f) {
01149 result[8] = ulpinv;
01150 }
01151 if (ifunc == 1 && pl[0] != 0.f) {
01152 result[8] = ulpinv;
01153 }
01154 }
01155
01156
01157
01158 ntest = 10;
01159 result[9] = 0.f;
01160 if (pl[0] == 0.f) {
01161 if (pltru > abnrm * ulp) {
01162 result[9] = ulpinv;
01163 }
01164 } else if (pltru == 0.f) {
01165 if (pl[0] > abnrm * ulp) {
01166 result[9] = ulpinv;
01167 }
01168 } else if (pltru > *thresh * pl[0] || pltru * *thresh < pl[0]) {
01169 result[9] = ulpinv;
01170 }
01171
01172 ntestt += ntest;
01173
01174
01175
01176 i__1 = ntest;
01177 for (j = 1; j <= i__1; ++j) {
01178 if (result[j - 1] >= *thresh) {
01179
01180
01181
01182
01183 if (nerrs == 0) {
01184 io___51.ciunit = *nout;
01185 s_wsfe(&io___51);
01186 do_fio(&c__1, "SGX", (ftnlen)3);
01187 e_wsfe();
01188
01189
01190
01191 io___52.ciunit = *nout;
01192 s_wsfe(&io___52);
01193 e_wsfe();
01194
01195
01196
01197 io___53.ciunit = *nout;
01198 s_wsfe(&io___53);
01199 do_fio(&c__1, "orthogonal", (ftnlen)10);
01200 do_fio(&c__1, "'", (ftnlen)1);
01201 do_fio(&c__1, "transpose", (ftnlen)9);
01202 for (i__ = 1; i__ <= 4; ++i__) {
01203 do_fio(&c__1, "'", (ftnlen)1);
01204 }
01205 e_wsfe();
01206
01207 }
01208 ++nerrs;
01209 if (result[j - 1] < 1e4f) {
01210 io___54.ciunit = *nout;
01211 s_wsfe(&io___54);
01212 do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
01213 do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)sizeof(integer));
01214 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
01215 do_fio(&c__1, (char *)&result[j - 1], (ftnlen)sizeof(real));
01216 e_wsfe();
01217 } else {
01218 io___55.ciunit = *nout;
01219 s_wsfe(&io___55);
01220 do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
01221 do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)sizeof(integer));
01222 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
01223 do_fio(&c__1, (char *)&result[j - 1], (ftnlen)sizeof(real));
01224 e_wsfe();
01225 }
01226 }
01227
01228
01229 }
01230
01231 L130:
01232 goto L80;
01233 L140:
01234
01235 L150:
01236
01237
01238
01239 alasvm_("SGX", nout, &nerrs, &ntestt, &c__0);
01240
01241 work[1] = (real) maxwrk;
01242
01243 return 0;
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255 }