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