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__3 = 3;
00019 static integer c__1 = 1;
00020 static integer c__6 = 6;
00021 static integer c__4 = 4;
00022 static integer c__20 = 20;
00023 static integer c__1200 = 1200;
00024 static integer c__0 = 0;
00025
00026 int cget37_(real *rmax, integer *lmax, integer *ninfo,
00027 integer *knt, integer *nin)
00028 {
00029
00030 integer i__1, i__2, i__3;
00031 real r__1, r__2;
00032
00033
00034 double sqrt(doublereal);
00035 integer s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
00036 e_rsle(void);
00037 double r_imag(complex *);
00038
00039
00040 integer i__, j, m, n;
00041 real s[20];
00042 complex t[400] ;
00043 real v;
00044 complex w[20], le[400] , re[400]
00045 ;
00046 real val[3], dum[1], eps, sep[20], sin__[20], tol;
00047 complex tmp[400] ;
00048 integer icmp;
00049 complex cdum[1];
00050 integer iscl, info, lcmp[3], kmin;
00051 real wiin[20], vmin, vmax, tnrm;
00052 integer isrt;
00053 real wrin[20], vmul, stmp[20];
00054 complex work[1200], wtmp[20];
00055 real wsrt[20];
00056 extern int sscal_(integer *, real *, real *, integer *);
00057 real vcmin;
00058 extern int ccopy_(integer *, complex *, integer *,
00059 complex *, integer *);
00060 real sepin[20], tolin;
00061 extern int scopy_(integer *, real *, integer *, real *,
00062 integer *);
00063 real rwork[40];
00064 extern int slabad_(real *, real *);
00065 extern doublereal clange_(char *, integer *, integer *, complex *,
00066 integer *, real *);
00067 extern int cgehrd_(integer *, integer *, integer *,
00068 complex *, integer *, complex *, complex *, integer *, integer *);
00069 extern doublereal slamch_(char *);
00070 extern int csscal_(integer *, real *, complex *, integer
00071 *), clacpy_(char *, integer *, integer *, complex *, integer *,
00072 complex *, integer *);
00073 logical select[20];
00074 real bignum;
00075 extern int chseqr_(char *, char *, integer *, integer *,
00076 integer *, complex *, integer *, complex *, complex *, integer *,
00077 complex *, integer *, integer *), ctrevc_(char *,
00078 char *, logical *, integer *, complex *, integer *, complex *,
00079 integer *, complex *, integer *, integer *, integer *, complex *,
00080 real *, integer *), ctrsna_(char *, char *,
00081 logical *, integer *, complex *, integer *, complex *, integer *,
00082 complex *, integer *, real *, real *, integer *, integer *,
00083 complex *, integer *, real *, integer *);
00084 real septmp[20], smlnum;
00085
00086
00087 static cilist io___5 = { 0, 0, 0, 0, 0 };
00088 static cilist io___9 = { 0, 0, 0, 0, 0 };
00089 static cilist io___12 = { 0, 0, 0, 0, 0 };
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
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 --ninfo;
00158 --lmax;
00159 --rmax;
00160
00161
00162 eps = slamch_("P");
00163 smlnum = slamch_("S") / eps;
00164 bignum = 1.f / smlnum;
00165 slabad_(&smlnum, &bignum);
00166
00167
00168
00169 eps = dmax(eps,5.9605e-8f);
00170 rmax[1] = 0.f;
00171 rmax[2] = 0.f;
00172 rmax[3] = 0.f;
00173 lmax[1] = 0;
00174 lmax[2] = 0;
00175 lmax[3] = 0;
00176 *knt = 0;
00177 ninfo[1] = 0;
00178 ninfo[2] = 0;
00179 ninfo[3] = 0;
00180 val[0] = sqrt(smlnum);
00181 val[1] = 1.f;
00182 val[2] = sqrt(bignum);
00183
00184
00185
00186
00187
00188 L10:
00189 io___5.ciunit = *nin;
00190 s_rsle(&io___5);
00191 do_lio(&c__3, &c__1, (char *)&n, (ftnlen)sizeof(integer));
00192 do_lio(&c__3, &c__1, (char *)&isrt, (ftnlen)sizeof(integer));
00193 e_rsle();
00194 if (n == 0) {
00195 return 0;
00196 }
00197 i__1 = n;
00198 for (i__ = 1; i__ <= i__1; ++i__) {
00199 io___9.ciunit = *nin;
00200 s_rsle(&io___9);
00201 i__2 = n;
00202 for (j = 1; j <= i__2; ++j) {
00203 do_lio(&c__6, &c__1, (char *)&tmp[i__ + j * 20 - 21], (ftnlen)
00204 sizeof(complex));
00205 }
00206 e_rsle();
00207
00208 }
00209 i__1 = n;
00210 for (i__ = 1; i__ <= i__1; ++i__) {
00211 io___12.ciunit = *nin;
00212 s_rsle(&io___12);
00213 do_lio(&c__4, &c__1, (char *)&wrin[i__ - 1], (ftnlen)sizeof(real));
00214 do_lio(&c__4, &c__1, (char *)&wiin[i__ - 1], (ftnlen)sizeof(real));
00215 do_lio(&c__4, &c__1, (char *)&sin__[i__ - 1], (ftnlen)sizeof(real));
00216 do_lio(&c__4, &c__1, (char *)&sepin[i__ - 1], (ftnlen)sizeof(real));
00217 e_rsle();
00218
00219 }
00220 tnrm = clange_("M", &n, &n, tmp, &c__20, rwork);
00221 for (iscl = 1; iscl <= 3; ++iscl) {
00222
00223
00224
00225 ++(*knt);
00226 clacpy_("F", &n, &n, tmp, &c__20, t, &c__20);
00227 vmul = val[iscl - 1];
00228 i__1 = n;
00229 for (i__ = 1; i__ <= i__1; ++i__) {
00230 csscal_(&n, &vmul, &t[i__ * 20 - 20], &c__1);
00231
00232 }
00233 if (tnrm == 0.f) {
00234 vmul = 1.f;
00235 }
00236
00237
00238
00239 i__1 = 1200 - n;
00240 cgehrd_(&n, &c__1, &n, t, &c__20, work, &work[n], &i__1, &info);
00241 if (info != 0) {
00242 lmax[1] = *knt;
00243 ++ninfo[1];
00244 goto L260;
00245 }
00246 i__1 = n - 2;
00247 for (j = 1; j <= i__1; ++j) {
00248 i__2 = n;
00249 for (i__ = j + 2; i__ <= i__2; ++i__) {
00250 i__3 = i__ + j * 20 - 21;
00251 t[i__3].r = 0.f, t[i__3].i = 0.f;
00252
00253 }
00254
00255 }
00256
00257
00258
00259 chseqr_("S", "N", &n, &c__1, &n, t, &c__20, w, cdum, &c__1, work, &
00260 c__1200, &info);
00261 if (info != 0) {
00262 lmax[2] = *knt;
00263 ++ninfo[2];
00264 goto L260;
00265 }
00266
00267
00268
00269 i__1 = n;
00270 for (i__ = 1; i__ <= i__1; ++i__) {
00271 select[i__ - 1] = TRUE_;
00272
00273 }
00274 ctrevc_("B", "A", select, &n, t, &c__20, le, &c__20, re, &c__20, &n, &
00275 m, work, rwork, &info);
00276
00277
00278
00279 ctrsna_("B", "A", select, &n, t, &c__20, le, &c__20, re, &c__20, s,
00280 sep, &n, &m, work, &n, rwork, &info);
00281 if (info != 0) {
00282 lmax[3] = *knt;
00283 ++ninfo[3];
00284 goto L260;
00285 }
00286
00287
00288
00289
00290 ccopy_(&n, w, &c__1, wtmp, &c__1);
00291 if (isrt == 0) {
00292
00293
00294
00295 i__1 = n;
00296 for (i__ = 1; i__ <= i__1; ++i__) {
00297 i__2 = i__ - 1;
00298 wsrt[i__ - 1] = w[i__2].r;
00299
00300 }
00301 } else {
00302
00303
00304
00305 i__1 = n;
00306 for (i__ = 1; i__ <= i__1; ++i__) {
00307 wsrt[i__ - 1] = r_imag(&w[i__ - 1]);
00308
00309 }
00310 }
00311 scopy_(&n, s, &c__1, stmp, &c__1);
00312 scopy_(&n, sep, &c__1, septmp, &c__1);
00313 r__1 = 1.f / vmul;
00314 sscal_(&n, &r__1, septmp, &c__1);
00315 i__1 = n - 1;
00316 for (i__ = 1; i__ <= i__1; ++i__) {
00317 kmin = i__;
00318 vmin = wsrt[i__ - 1];
00319 i__2 = n;
00320 for (j = i__ + 1; j <= i__2; ++j) {
00321 if (wsrt[j - 1] < vmin) {
00322 kmin = j;
00323 vmin = wsrt[j - 1];
00324 }
00325
00326 }
00327 wsrt[kmin - 1] = wsrt[i__ - 1];
00328 wsrt[i__ - 1] = vmin;
00329 i__2 = i__ - 1;
00330 vcmin = wtmp[i__2].r;
00331 i__2 = i__ - 1;
00332 i__3 = kmin - 1;
00333 wtmp[i__2].r = w[i__3].r, wtmp[i__2].i = w[i__3].i;
00334 i__2 = kmin - 1;
00335 wtmp[i__2].r = vcmin, wtmp[i__2].i = 0.f;
00336 vmin = stmp[kmin - 1];
00337 stmp[kmin - 1] = stmp[i__ - 1];
00338 stmp[i__ - 1] = vmin;
00339 vmin = septmp[kmin - 1];
00340 septmp[kmin - 1] = septmp[i__ - 1];
00341 septmp[i__ - 1] = vmin;
00342
00343 }
00344
00345
00346
00347
00348
00349 r__1 = (real) n * 2.f * eps * tnrm;
00350 v = dmax(r__1,smlnum);
00351 if (tnrm == 0.f) {
00352 v = 1.f;
00353 }
00354 i__1 = n;
00355 for (i__ = 1; i__ <= i__1; ++i__) {
00356 if (v > septmp[i__ - 1]) {
00357 tol = 1.f;
00358 } else {
00359 tol = v / septmp[i__ - 1];
00360 }
00361 if (v > sepin[i__ - 1]) {
00362 tolin = 1.f;
00363 } else {
00364 tolin = v / sepin[i__ - 1];
00365 }
00366
00367 r__1 = tol, r__2 = smlnum / eps;
00368 tol = dmax(r__1,r__2);
00369
00370 r__1 = tolin, r__2 = smlnum / eps;
00371 tolin = dmax(r__1,r__2);
00372 if (eps * (sin__[i__ - 1] - tolin) > stmp[i__ - 1] + tol) {
00373 vmax = 1.f / eps;
00374 } else if (sin__[i__ - 1] - tolin > stmp[i__ - 1] + tol) {
00375 vmax = (sin__[i__ - 1] - tolin) / (stmp[i__ - 1] + tol);
00376 } else if (sin__[i__ - 1] + tolin < eps * (stmp[i__ - 1] - tol)) {
00377 vmax = 1.f / eps;
00378 } else if (sin__[i__ - 1] + tolin < stmp[i__ - 1] - tol) {
00379 vmax = (stmp[i__ - 1] - tol) / (sin__[i__ - 1] + tolin);
00380 } else {
00381 vmax = 1.f;
00382 }
00383 if (vmax > rmax[2]) {
00384 rmax[2] = vmax;
00385 if (ninfo[2] == 0) {
00386 lmax[2] = *knt;
00387 }
00388 }
00389
00390 }
00391
00392
00393
00394
00395 i__1 = n;
00396 for (i__ = 1; i__ <= i__1; ++i__) {
00397 if (v > septmp[i__ - 1] * stmp[i__ - 1]) {
00398 tol = septmp[i__ - 1];
00399 } else {
00400 tol = v / stmp[i__ - 1];
00401 }
00402 if (v > sepin[i__ - 1] * sin__[i__ - 1]) {
00403 tolin = sepin[i__ - 1];
00404 } else {
00405 tolin = v / sin__[i__ - 1];
00406 }
00407
00408 r__1 = tol, r__2 = smlnum / eps;
00409 tol = dmax(r__1,r__2);
00410
00411 r__1 = tolin, r__2 = smlnum / eps;
00412 tolin = dmax(r__1,r__2);
00413 if (eps * (sepin[i__ - 1] - tolin) > septmp[i__ - 1] + tol) {
00414 vmax = 1.f / eps;
00415 } else if (sepin[i__ - 1] - tolin > septmp[i__ - 1] + tol) {
00416 vmax = (sepin[i__ - 1] - tolin) / (septmp[i__ - 1] + tol);
00417 } else if (sepin[i__ - 1] + tolin < eps * (septmp[i__ - 1] - tol))
00418 {
00419 vmax = 1.f / eps;
00420 } else if (sepin[i__ - 1] + tolin < septmp[i__ - 1] - tol) {
00421 vmax = (septmp[i__ - 1] - tol) / (sepin[i__ - 1] + tolin);
00422 } else {
00423 vmax = 1.f;
00424 }
00425 if (vmax > rmax[2]) {
00426 rmax[2] = vmax;
00427 if (ninfo[2] == 0) {
00428 lmax[2] = *knt;
00429 }
00430 }
00431
00432 }
00433
00434
00435
00436
00437 i__1 = n;
00438 for (i__ = 1; i__ <= i__1; ++i__) {
00439 if (sin__[i__ - 1] <= (real) (n << 1) * eps && stmp[i__ - 1] <= (
00440 real) (n << 1) * eps) {
00441 vmax = 1.f;
00442 } else if (eps * sin__[i__ - 1] > stmp[i__ - 1]) {
00443 vmax = 1.f / eps;
00444 } else if (sin__[i__ - 1] > stmp[i__ - 1]) {
00445 vmax = sin__[i__ - 1] / stmp[i__ - 1];
00446 } else if (sin__[i__ - 1] < eps * stmp[i__ - 1]) {
00447 vmax = 1.f / eps;
00448 } else if (sin__[i__ - 1] < stmp[i__ - 1]) {
00449 vmax = stmp[i__ - 1] / sin__[i__ - 1];
00450 } else {
00451 vmax = 1.f;
00452 }
00453 if (vmax > rmax[3]) {
00454 rmax[3] = vmax;
00455 if (ninfo[3] == 0) {
00456 lmax[3] = *knt;
00457 }
00458 }
00459
00460 }
00461
00462
00463
00464
00465 i__1 = n;
00466 for (i__ = 1; i__ <= i__1; ++i__) {
00467 if (sepin[i__ - 1] <= v && septmp[i__ - 1] <= v) {
00468 vmax = 1.f;
00469 } else if (eps * sepin[i__ - 1] > septmp[i__ - 1]) {
00470 vmax = 1.f / eps;
00471 } else if (sepin[i__ - 1] > septmp[i__ - 1]) {
00472 vmax = sepin[i__ - 1] / septmp[i__ - 1];
00473 } else if (sepin[i__ - 1] < eps * septmp[i__ - 1]) {
00474 vmax = 1.f / eps;
00475 } else if (sepin[i__ - 1] < septmp[i__ - 1]) {
00476 vmax = septmp[i__ - 1] / sepin[i__ - 1];
00477 } else {
00478 vmax = 1.f;
00479 }
00480 if (vmax > rmax[3]) {
00481 rmax[3] = vmax;
00482 if (ninfo[3] == 0) {
00483 lmax[3] = *knt;
00484 }
00485 }
00486
00487 }
00488
00489
00490
00491 vmax = 0.f;
00492 dum[0] = -1.f;
00493 scopy_(&n, dum, &c__0, stmp, &c__1);
00494 scopy_(&n, dum, &c__0, septmp, &c__1);
00495 ctrsna_("E", "A", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp,
00496 septmp, &n, &m, work, &n, rwork, &info)
00497 ;
00498 if (info != 0) {
00499 lmax[3] = *knt;
00500 ++ninfo[3];
00501 goto L260;
00502 }
00503 i__1 = n;
00504 for (i__ = 1; i__ <= i__1; ++i__) {
00505 if (stmp[i__ - 1] != s[i__ - 1]) {
00506 vmax = 1.f / eps;
00507 }
00508 if (septmp[i__ - 1] != dum[0]) {
00509 vmax = 1.f / eps;
00510 }
00511
00512 }
00513
00514
00515
00516 scopy_(&n, dum, &c__0, stmp, &c__1);
00517 scopy_(&n, dum, &c__0, septmp, &c__1);
00518 ctrsna_("V", "A", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp,
00519 septmp, &n, &m, work, &n, rwork, &info)
00520 ;
00521 if (info != 0) {
00522 lmax[3] = *knt;
00523 ++ninfo[3];
00524 goto L260;
00525 }
00526 i__1 = n;
00527 for (i__ = 1; i__ <= i__1; ++i__) {
00528 if (stmp[i__ - 1] != dum[0]) {
00529 vmax = 1.f / eps;
00530 }
00531 if (septmp[i__ - 1] != sep[i__ - 1]) {
00532 vmax = 1.f / eps;
00533 }
00534
00535 }
00536
00537
00538
00539 i__1 = n;
00540 for (i__ = 1; i__ <= i__1; ++i__) {
00541 select[i__ - 1] = TRUE_;
00542
00543 }
00544 scopy_(&n, dum, &c__0, stmp, &c__1);
00545 scopy_(&n, dum, &c__0, septmp, &c__1);
00546 ctrsna_("B", "S", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp,
00547 septmp, &n, &m, work, &n, rwork, &info)
00548 ;
00549 if (info != 0) {
00550 lmax[3] = *knt;
00551 ++ninfo[3];
00552 goto L260;
00553 }
00554 i__1 = n;
00555 for (i__ = 1; i__ <= i__1; ++i__) {
00556 if (septmp[i__ - 1] != sep[i__ - 1]) {
00557 vmax = 1.f / eps;
00558 }
00559 if (stmp[i__ - 1] != s[i__ - 1]) {
00560 vmax = 1.f / eps;
00561 }
00562
00563 }
00564
00565
00566
00567 scopy_(&n, dum, &c__0, stmp, &c__1);
00568 scopy_(&n, dum, &c__0, septmp, &c__1);
00569 ctrsna_("E", "S", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp,
00570 septmp, &n, &m, work, &n, rwork, &info)
00571 ;
00572 if (info != 0) {
00573 lmax[3] = *knt;
00574 ++ninfo[3];
00575 goto L260;
00576 }
00577 i__1 = n;
00578 for (i__ = 1; i__ <= i__1; ++i__) {
00579 if (stmp[i__ - 1] != s[i__ - 1]) {
00580 vmax = 1.f / eps;
00581 }
00582 if (septmp[i__ - 1] != dum[0]) {
00583 vmax = 1.f / eps;
00584 }
00585
00586 }
00587
00588
00589
00590 scopy_(&n, dum, &c__0, stmp, &c__1);
00591 scopy_(&n, dum, &c__0, septmp, &c__1);
00592 ctrsna_("V", "S", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp,
00593 septmp, &n, &m, work, &n, rwork, &info)
00594 ;
00595 if (info != 0) {
00596 lmax[3] = *knt;
00597 ++ninfo[3];
00598 goto L260;
00599 }
00600 i__1 = n;
00601 for (i__ = 1; i__ <= i__1; ++i__) {
00602 if (stmp[i__ - 1] != dum[0]) {
00603 vmax = 1.f / eps;
00604 }
00605 if (septmp[i__ - 1] != sep[i__ - 1]) {
00606 vmax = 1.f / eps;
00607 }
00608
00609 }
00610 if (vmax > rmax[1]) {
00611 rmax[1] = vmax;
00612 if (ninfo[1] == 0) {
00613 lmax[1] = *knt;
00614 }
00615 }
00616
00617
00618
00619 i__1 = n;
00620 for (i__ = 1; i__ <= i__1; ++i__) {
00621 select[i__ - 1] = FALSE_;
00622
00623 }
00624 icmp = 0;
00625 if (n > 1) {
00626 icmp = 1;
00627 lcmp[0] = 2;
00628 select[1] = TRUE_;
00629 ccopy_(&n, &re[20], &c__1, re, &c__1);
00630 ccopy_(&n, &le[20], &c__1, le, &c__1);
00631 }
00632 if (n > 3) {
00633 icmp = 2;
00634 lcmp[1] = n - 1;
00635 select[n - 2] = TRUE_;
00636 ccopy_(&n, &re[(n - 1) * 20 - 20], &c__1, &re[20], &c__1);
00637 ccopy_(&n, &le[(n - 1) * 20 - 20], &c__1, &le[20], &c__1);
00638 }
00639
00640
00641
00642 scopy_(&icmp, dum, &c__0, stmp, &c__1);
00643 scopy_(&icmp, dum, &c__0, septmp, &c__1);
00644 ctrsna_("B", "S", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp,
00645 septmp, &n, &m, work, &n, rwork, &info)
00646 ;
00647 if (info != 0) {
00648 lmax[3] = *knt;
00649 ++ninfo[3];
00650 goto L260;
00651 }
00652 i__1 = icmp;
00653 for (i__ = 1; i__ <= i__1; ++i__) {
00654 j = lcmp[i__ - 1];
00655 if (septmp[i__ - 1] != sep[j - 1]) {
00656 vmax = 1.f / eps;
00657 }
00658 if (stmp[i__ - 1] != s[j - 1]) {
00659 vmax = 1.f / eps;
00660 }
00661
00662 }
00663
00664
00665
00666 scopy_(&icmp, dum, &c__0, stmp, &c__1);
00667 scopy_(&icmp, dum, &c__0, septmp, &c__1);
00668 ctrsna_("E", "S", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp,
00669 septmp, &n, &m, work, &n, rwork, &info)
00670 ;
00671 if (info != 0) {
00672 lmax[3] = *knt;
00673 ++ninfo[3];
00674 goto L260;
00675 }
00676 i__1 = icmp;
00677 for (i__ = 1; i__ <= i__1; ++i__) {
00678 j = lcmp[i__ - 1];
00679 if (stmp[i__ - 1] != s[j - 1]) {
00680 vmax = 1.f / eps;
00681 }
00682 if (septmp[i__ - 1] != dum[0]) {
00683 vmax = 1.f / eps;
00684 }
00685
00686 }
00687
00688
00689
00690 scopy_(&icmp, dum, &c__0, stmp, &c__1);
00691 scopy_(&icmp, dum, &c__0, septmp, &c__1);
00692 ctrsna_("V", "S", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp,
00693 septmp, &n, &m, work, &n, rwork, &info)
00694 ;
00695 if (info != 0) {
00696 lmax[3] = *knt;
00697 ++ninfo[3];
00698 goto L260;
00699 }
00700 i__1 = icmp;
00701 for (i__ = 1; i__ <= i__1; ++i__) {
00702 j = lcmp[i__ - 1];
00703 if (stmp[i__ - 1] != dum[0]) {
00704 vmax = 1.f / eps;
00705 }
00706 if (septmp[i__ - 1] != sep[j - 1]) {
00707 vmax = 1.f / eps;
00708 }
00709
00710 }
00711 if (vmax > rmax[1]) {
00712 rmax[1] = vmax;
00713 if (ninfo[1] == 0) {
00714 lmax[1] = *knt;
00715 }
00716 }
00717 L260:
00718 ;
00719 }
00720 goto L10;
00721
00722
00723
00724 }