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