00001
00002
00003
00004
00005 #include "f2c.h"
00006
00007 #ifdef HAVE_CONFIG
00008 #include "config.h"
00009 #else
00010 extern doublereal slamch_(char *);
00011 #define EPSILON slamch_("Epsilon")
00012 #define SAFEMINIMUM slamch_("Safe minimum")
00013 #define PRECISION slamch_("Precision")
00014 #define BASE slamch_("Base")
00015 #endif
00016
00017
00018 extern doublereal slapy2_(real *, real *);
00019
00020
00021
00022
00023
00024 static integer c__0 = 0;
00025 static real c_b163 = 0.f;
00026 static real c_b164 = 1.f;
00027 static integer c__1 = 1;
00028 static real c_b181 = -1.f;
00029 static integer c_n1 = -1;
00030
00031 integer ieeeck_(integer *ispec, real *zero, real *one)
00032 {
00033
00034 integer ret_val;
00035
00036
00037 static real nan1, nan2, nan3, nan4, nan5, nan6, neginf, posinf, negzro,
00038 newzro;
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078 ret_val = 1;
00079
00080 posinf = *one / *zero;
00081 if (posinf <= *one) {
00082 ret_val = 0;
00083 return ret_val;
00084 }
00085
00086 neginf = -(*one) / *zero;
00087 if (neginf >= *zero) {
00088 ret_val = 0;
00089 return ret_val;
00090 }
00091
00092 negzro = *one / (neginf + *one);
00093 if (negzro != *zero) {
00094 ret_val = 0;
00095 return ret_val;
00096 }
00097
00098 neginf = *one / negzro;
00099 if (neginf >= *zero) {
00100 ret_val = 0;
00101 return ret_val;
00102 }
00103
00104 newzro = negzro + *zero;
00105 if (newzro != *zero) {
00106 ret_val = 0;
00107 return ret_val;
00108 }
00109
00110 posinf = *one / newzro;
00111 if (posinf <= *one) {
00112 ret_val = 0;
00113 return ret_val;
00114 }
00115
00116 neginf *= posinf;
00117 if (neginf >= *zero) {
00118 ret_val = 0;
00119 return ret_val;
00120 }
00121
00122 posinf *= posinf;
00123 if (posinf <= *one) {
00124 ret_val = 0;
00125 return ret_val;
00126 }
00127
00128
00129
00130
00131 if (*ispec == 0) {
00132 return ret_val;
00133 }
00134
00135 nan1 = posinf + neginf;
00136
00137 nan2 = posinf / neginf;
00138
00139 nan3 = posinf / posinf;
00140
00141 nan4 = posinf * *zero;
00142
00143 nan5 = neginf * negzro;
00144
00145 nan6 = nan5 * 0.f;
00146
00147 if (nan1 == nan1) {
00148 ret_val = 0;
00149 return ret_val;
00150 }
00151
00152 if (nan2 == nan2) {
00153 ret_val = 0;
00154 return ret_val;
00155 }
00156
00157 if (nan3 == nan3) {
00158 ret_val = 0;
00159 return ret_val;
00160 }
00161
00162 if (nan4 == nan4) {
00163 ret_val = 0;
00164 return ret_val;
00165 }
00166
00167 if (nan5 == nan5) {
00168 ret_val = 0;
00169 return ret_val;
00170 }
00171
00172 if (nan6 == nan6) {
00173 ret_val = 0;
00174 return ret_val;
00175 }
00176
00177 return ret_val;
00178 }
00179
00180 integer ilaenv_(integer *ispec, char *name__, char *opts, integer *n1,
00181 integer *n2, integer *n3, integer *n4, ftnlen name_len, ftnlen
00182 opts_len)
00183 {
00184
00185 integer ret_val;
00186
00187
00188 int s_copy(char *, char *, ftnlen, ftnlen);
00189 integer s_cmp(char *, char *, ftnlen, ftnlen);
00190
00191
00192 static integer i__;
00193 static char c1[1], c2[2], c3[3], c4[2];
00194 static integer ic, nb, iz, nx;
00195 static logical cname, sname;
00196 static integer nbmin;
00197 extern integer ieeeck_(integer *, real *, real *);
00198 static char subnam[6];
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 switch (*ispec) {
00301 case 1: goto L100;
00302 case 2: goto L100;
00303 case 3: goto L100;
00304 case 4: goto L400;
00305 case 5: goto L500;
00306 case 6: goto L600;
00307 case 7: goto L700;
00308 case 8: goto L800;
00309 case 9: goto L900;
00310 case 10: goto L1000;
00311 case 11: goto L1100;
00312 }
00313
00314
00315
00316 ret_val = -1;
00317 return ret_val;
00318
00319 L100:
00320
00321
00322
00323 ret_val = 1;
00324 s_copy(subnam, name__, (ftnlen)6, name_len);
00325 ic = *(unsigned char *)subnam;
00326 iz = 'Z';
00327 if (iz == 90 || iz == 122) {
00328
00329
00330
00331 if (ic >= 97 && ic <= 122) {
00332 *(unsigned char *)subnam = (char) (ic - 32);
00333 for (i__ = 2; i__ <= 6; ++i__) {
00334 ic = *(unsigned char *)&subnam[i__ - 1];
00335 if (ic >= 97 && ic <= 122) {
00336 *(unsigned char *)&subnam[i__ - 1] = (char) (ic - 32);
00337 }
00338
00339 }
00340 }
00341
00342 } else if (iz == 233 || iz == 169) {
00343
00344
00345
00346 if (ic >= 129 && ic <= 137 || ic >= 145 && ic <= 153 || ic >= 162 &&
00347 ic <= 169) {
00348 *(unsigned char *)subnam = (char) (ic + 64);
00349 for (i__ = 2; i__ <= 6; ++i__) {
00350 ic = *(unsigned char *)&subnam[i__ - 1];
00351 if (ic >= 129 && ic <= 137 || ic >= 145 && ic <= 153 || ic >=
00352 162 && ic <= 169) {
00353 *(unsigned char *)&subnam[i__ - 1] = (char) (ic + 64);
00354 }
00355
00356 }
00357 }
00358
00359 } else if (iz == 218 || iz == 250) {
00360
00361
00362
00363 if (ic >= 225 && ic <= 250) {
00364 *(unsigned char *)subnam = (char) (ic - 32);
00365 for (i__ = 2; i__ <= 6; ++i__) {
00366 ic = *(unsigned char *)&subnam[i__ - 1];
00367 if (ic >= 225 && ic <= 250) {
00368 *(unsigned char *)&subnam[i__ - 1] = (char) (ic - 32);
00369 }
00370
00371 }
00372 }
00373 }
00374
00375 *(unsigned char *)c1 = *(unsigned char *)subnam;
00376 sname = *(unsigned char *)c1 == 'S' || *(unsigned char *)c1 == 'D';
00377 cname = *(unsigned char *)c1 == 'C' || *(unsigned char *)c1 == 'Z';
00378 if (! (cname || sname)) {
00379 return ret_val;
00380 }
00381 s_copy(c2, subnam + 1, (ftnlen)2, (ftnlen)2);
00382 s_copy(c3, subnam + 3, (ftnlen)3, (ftnlen)3);
00383 s_copy(c4, c3 + 1, (ftnlen)2, (ftnlen)2);
00384
00385 switch (*ispec) {
00386 case 1: goto L110;
00387 case 2: goto L200;
00388 case 3: goto L300;
00389 }
00390
00391 L110:
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401 nb = 1;
00402
00403 if (s_cmp(c2, "GE", (ftnlen)2, (ftnlen)2) == 0) {
00404 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
00405 if (sname) {
00406 nb = 64;
00407 } else {
00408 nb = 64;
00409 }
00410 } else if (s_cmp(c3, "QRF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3,
00411 "RQF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)
00412 3, (ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)3, (ftnlen)3)
00413 == 0) {
00414 if (sname) {
00415 nb = 32;
00416 } else {
00417 nb = 32;
00418 }
00419 } else if (s_cmp(c3, "HRD", (ftnlen)3, (ftnlen)3) == 0) {
00420 if (sname) {
00421 nb = 32;
00422 } else {
00423 nb = 32;
00424 }
00425 } else if (s_cmp(c3, "BRD", (ftnlen)3, (ftnlen)3) == 0) {
00426 if (sname) {
00427 nb = 32;
00428 } else {
00429 nb = 32;
00430 }
00431 } else if (s_cmp(c3, "TRI", (ftnlen)3, (ftnlen)3) == 0) {
00432 if (sname) {
00433 nb = 64;
00434 } else {
00435 nb = 64;
00436 }
00437 }
00438 } else if (s_cmp(c2, "PO", (ftnlen)2, (ftnlen)2) == 0) {
00439 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
00440 if (sname) {
00441 nb = 64;
00442 } else {
00443 nb = 64;
00444 }
00445 }
00446 } else if (s_cmp(c2, "SY", (ftnlen)2, (ftnlen)2) == 0) {
00447 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
00448 if (sname) {
00449 nb = 64;
00450 } else {
00451 nb = 64;
00452 }
00453 } else if (sname && s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
00454 nb = 32;
00455 } else if (sname && s_cmp(c3, "GST", (ftnlen)3, (ftnlen)3) == 0) {
00456 nb = 64;
00457 }
00458 } else if (cname && s_cmp(c2, "HE", (ftnlen)2, (ftnlen)2) == 0) {
00459 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
00460 nb = 64;
00461 } else if (s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
00462 nb = 32;
00463 } else if (s_cmp(c3, "GST", (ftnlen)3, (ftnlen)3) == 0) {
00464 nb = 64;
00465 }
00466 } else if (sname && s_cmp(c2, "OR", (ftnlen)2, (ftnlen)2) == 0) {
00467 if (*(unsigned char *)c3 == 'G') {
00468 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
00469 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
00470 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
00471 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
00472 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
00473 ftnlen)2, (ftnlen)2) == 0) {
00474 nb = 32;
00475 }
00476 } else if (*(unsigned char *)c3 == 'M') {
00477 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
00478 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
00479 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
00480 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
00481 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
00482 ftnlen)2, (ftnlen)2) == 0) {
00483 nb = 32;
00484 }
00485 }
00486 } else if (cname && s_cmp(c2, "UN", (ftnlen)2, (ftnlen)2) == 0) {
00487 if (*(unsigned char *)c3 == 'G') {
00488 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
00489 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
00490 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
00491 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
00492 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
00493 ftnlen)2, (ftnlen)2) == 0) {
00494 nb = 32;
00495 }
00496 } else if (*(unsigned char *)c3 == 'M') {
00497 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
00498 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
00499 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
00500 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
00501 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
00502 ftnlen)2, (ftnlen)2) == 0) {
00503 nb = 32;
00504 }
00505 }
00506 } else if (s_cmp(c2, "GB", (ftnlen)2, (ftnlen)2) == 0) {
00507 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
00508 if (sname) {
00509 if (*n4 <= 64) {
00510 nb = 1;
00511 } else {
00512 nb = 32;
00513 }
00514 } else {
00515 if (*n4 <= 64) {
00516 nb = 1;
00517 } else {
00518 nb = 32;
00519 }
00520 }
00521 }
00522 } else if (s_cmp(c2, "PB", (ftnlen)2, (ftnlen)2) == 0) {
00523 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
00524 if (sname) {
00525 if (*n2 <= 64) {
00526 nb = 1;
00527 } else {
00528 nb = 32;
00529 }
00530 } else {
00531 if (*n2 <= 64) {
00532 nb = 1;
00533 } else {
00534 nb = 32;
00535 }
00536 }
00537 }
00538 } else if (s_cmp(c2, "TR", (ftnlen)2, (ftnlen)2) == 0) {
00539 if (s_cmp(c3, "TRI", (ftnlen)3, (ftnlen)3) == 0) {
00540 if (sname) {
00541 nb = 64;
00542 } else {
00543 nb = 64;
00544 }
00545 }
00546 } else if (s_cmp(c2, "LA", (ftnlen)2, (ftnlen)2) == 0) {
00547 if (s_cmp(c3, "UUM", (ftnlen)3, (ftnlen)3) == 0) {
00548 if (sname) {
00549 nb = 64;
00550 } else {
00551 nb = 64;
00552 }
00553 }
00554 } else if (sname && s_cmp(c2, "ST", (ftnlen)2, (ftnlen)2) == 0) {
00555 if (s_cmp(c3, "EBZ", (ftnlen)3, (ftnlen)3) == 0) {
00556 nb = 1;
00557 }
00558 }
00559 ret_val = nb;
00560 return ret_val;
00561
00562 L200:
00563
00564
00565
00566 nbmin = 2;
00567 if (s_cmp(c2, "GE", (ftnlen)2, (ftnlen)2) == 0) {
00568 if (s_cmp(c3, "QRF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "RQF", (
00569 ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)3, (
00570 ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)3, (ftnlen)3) == 0)
00571 {
00572 if (sname) {
00573 nbmin = 2;
00574 } else {
00575 nbmin = 2;
00576 }
00577 } else if (s_cmp(c3, "HRD", (ftnlen)3, (ftnlen)3) == 0) {
00578 if (sname) {
00579 nbmin = 2;
00580 } else {
00581 nbmin = 2;
00582 }
00583 } else if (s_cmp(c3, "BRD", (ftnlen)3, (ftnlen)3) == 0) {
00584 if (sname) {
00585 nbmin = 2;
00586 } else {
00587 nbmin = 2;
00588 }
00589 } else if (s_cmp(c3, "TRI", (ftnlen)3, (ftnlen)3) == 0) {
00590 if (sname) {
00591 nbmin = 2;
00592 } else {
00593 nbmin = 2;
00594 }
00595 }
00596 } else if (s_cmp(c2, "SY", (ftnlen)2, (ftnlen)2) == 0) {
00597 if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
00598 if (sname) {
00599 nbmin = 8;
00600 } else {
00601 nbmin = 8;
00602 }
00603 } else if (sname && s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
00604 nbmin = 2;
00605 }
00606 } else if (cname && s_cmp(c2, "HE", (ftnlen)2, (ftnlen)2) == 0) {
00607 if (s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
00608 nbmin = 2;
00609 }
00610 } else if (sname && s_cmp(c2, "OR", (ftnlen)2, (ftnlen)2) == 0) {
00611 if (*(unsigned char *)c3 == 'G') {
00612 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
00613 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
00614 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
00615 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
00616 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
00617 ftnlen)2, (ftnlen)2) == 0) {
00618 nbmin = 2;
00619 }
00620 } else if (*(unsigned char *)c3 == 'M') {
00621 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
00622 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
00623 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
00624 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
00625 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
00626 ftnlen)2, (ftnlen)2) == 0) {
00627 nbmin = 2;
00628 }
00629 }
00630 } else if (cname && s_cmp(c2, "UN", (ftnlen)2, (ftnlen)2) == 0) {
00631 if (*(unsigned char *)c3 == 'G') {
00632 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
00633 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
00634 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
00635 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
00636 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
00637 ftnlen)2, (ftnlen)2) == 0) {
00638 nbmin = 2;
00639 }
00640 } else if (*(unsigned char *)c3 == 'M') {
00641 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
00642 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
00643 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
00644 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
00645 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
00646 ftnlen)2, (ftnlen)2) == 0) {
00647 nbmin = 2;
00648 }
00649 }
00650 }
00651 ret_val = nbmin;
00652 return ret_val;
00653
00654 L300:
00655
00656
00657
00658 nx = 0;
00659 if (s_cmp(c2, "GE", (ftnlen)2, (ftnlen)2) == 0) {
00660 if (s_cmp(c3, "QRF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "RQF", (
00661 ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)3, (
00662 ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)3, (ftnlen)3) == 0)
00663 {
00664 if (sname) {
00665 nx = 128;
00666 } else {
00667 nx = 128;
00668 }
00669 } else if (s_cmp(c3, "HRD", (ftnlen)3, (ftnlen)3) == 0) {
00670 if (sname) {
00671 nx = 128;
00672 } else {
00673 nx = 128;
00674 }
00675 } else if (s_cmp(c3, "BRD", (ftnlen)3, (ftnlen)3) == 0) {
00676 if (sname) {
00677 nx = 128;
00678 } else {
00679 nx = 128;
00680 }
00681 }
00682 } else if (s_cmp(c2, "SY", (ftnlen)2, (ftnlen)2) == 0) {
00683 if (sname && s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
00684 nx = 32;
00685 }
00686 } else if (cname && s_cmp(c2, "HE", (ftnlen)2, (ftnlen)2) == 0) {
00687 if (s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
00688 nx = 32;
00689 }
00690 } else if (sname && s_cmp(c2, "OR", (ftnlen)2, (ftnlen)2) == 0) {
00691 if (*(unsigned char *)c3 == 'G') {
00692 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
00693 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
00694 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
00695 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
00696 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
00697 ftnlen)2, (ftnlen)2) == 0) {
00698 nx = 128;
00699 }
00700 }
00701 } else if (cname && s_cmp(c2, "UN", (ftnlen)2, (ftnlen)2) == 0) {
00702 if (*(unsigned char *)c3 == 'G') {
00703 if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
00704 (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
00705 ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
00706 0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
00707 c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
00708 ftnlen)2, (ftnlen)2) == 0) {
00709 nx = 128;
00710 }
00711 }
00712 }
00713 ret_val = nx;
00714 return ret_val;
00715
00716 L400:
00717
00718
00719
00720 ret_val = 6;
00721 return ret_val;
00722
00723 L500:
00724
00725
00726
00727 ret_val = 2;
00728 return ret_val;
00729
00730 L600:
00731
00732
00733
00734 ret_val = (integer) ((real) min(*n1,*n2) * 1.6f);
00735 return ret_val;
00736
00737 L700:
00738
00739
00740
00741 ret_val = 1;
00742 return ret_val;
00743
00744 L800:
00745
00746
00747
00748 ret_val = 50;
00749 return ret_val;
00750
00751 L900:
00752
00753
00754
00755
00756
00757
00758
00759 ret_val = 25;
00760 return ret_val;
00761
00762 L1000:
00763
00764
00765
00766
00767
00768
00769 ret_val = 1;
00770 if (ret_val == 1) {
00771 ret_val = ieeeck_(&c__0, &c_b163, &c_b164);
00772 }
00773 return ret_val;
00774
00775 L1100:
00776
00777
00778
00779
00780
00781
00782 ret_val = 1;
00783 if (ret_val == 1) {
00784 ret_val = ieeeck_(&c__1, &c_b163, &c_b164);
00785 }
00786 return ret_val;
00787
00788
00789
00790 }
00791
00792 int sposv_(char *uplo, integer *n, integer *nrhs, real *a,
00793 integer *lda, real *b, integer *ldb, integer *info)
00794 {
00795
00796 integer a_dim1, a_offset, b_dim1, b_offset, i__1;
00797
00798
00799 extern logical lsame_(char *, char *);
00800 extern int xerbla_(char *, integer *), spotrf_(
00801 char *, integer *, real *, integer *, integer *), spotrs_(
00802 char *, integer *, integer *, real *, integer *, real *, integer *
00803 , integer *);
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879 a_dim1 = *lda;
00880 a_offset = 1 + a_dim1;
00881 a -= a_offset;
00882 b_dim1 = *ldb;
00883 b_offset = 1 + b_dim1;
00884 b -= b_offset;
00885
00886
00887 *info = 0;
00888 if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
00889 *info = -1;
00890 } else if (*n < 0) {
00891 *info = -2;
00892 } else if (*nrhs < 0) {
00893 *info = -3;
00894 } else if (*lda < max(1,*n)) {
00895 *info = -5;
00896 } else if (*ldb < max(1,*n)) {
00897 *info = -7;
00898 }
00899 if (*info != 0) {
00900 i__1 = -(*info);
00901 xerbla_("SPOSV ", &i__1);
00902 return 0;
00903 }
00904
00905
00906
00907 spotrf_(uplo, n, &a[a_offset], lda, info);
00908 if (*info == 0) {
00909
00910
00911
00912 spotrs_(uplo, n, nrhs, &a[a_offset], lda, &b[b_offset], ldb, info);
00913
00914 }
00915 return 0;
00916
00917
00918
00919 }
00920
00921 int spotf2_(char *uplo, integer *n, real *a, integer *lda,
00922 integer *info)
00923 {
00924
00925 integer a_dim1, a_offset, i__1, i__2, i__3;
00926 real r__1;
00927
00928
00929 double sqrt(doublereal);
00930
00931
00932 static integer j;
00933 static real ajj;
00934 extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
00935 extern logical lsame_(char *, char *);
00936 extern int sscal_(integer *, real *, real *, integer *),
00937 sgemv_(char *, integer *, integer *, real *, real *, integer *,
00938 real *, integer *, real *, real *, integer *);
00939 static logical upper;
00940 extern int xerbla_(char *, integer *);
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004 a_dim1 = *lda;
01005 a_offset = 1 + a_dim1;
01006 a -= a_offset;
01007
01008
01009 *info = 0;
01010 upper = lsame_(uplo, "U");
01011 if (! upper && ! lsame_(uplo, "L")) {
01012 *info = -1;
01013 } else if (*n < 0) {
01014 *info = -2;
01015 } else if (*lda < max(1,*n)) {
01016 *info = -4;
01017 }
01018 if (*info != 0) {
01019 i__1 = -(*info);
01020 xerbla_("SPOTF2", &i__1);
01021 return 0;
01022 }
01023
01024
01025
01026 if (*n == 0) {
01027 return 0;
01028 }
01029
01030 if (upper) {
01031
01032
01033
01034 i__1 = *n;
01035 for (j = 1; j <= i__1; ++j) {
01036
01037
01038
01039 i__2 = j - 1;
01040 ajj = a[j + j * a_dim1] - sdot_(&i__2, &a[j * a_dim1 + 1], &c__1,
01041 &a[j * a_dim1 + 1], &c__1);
01042 if (ajj <= 0.f) {
01043 a[j + j * a_dim1] = ajj;
01044 goto L30;
01045 }
01046 ajj = sqrt(ajj);
01047 a[j + j * a_dim1] = ajj;
01048
01049
01050
01051 if (j < *n) {
01052 i__2 = j - 1;
01053 i__3 = *n - j;
01054 sgemv_("Transpose", &i__2, &i__3, &c_b181, &a[(j + 1) *
01055 a_dim1 + 1], lda, &a[j * a_dim1 + 1], &c__1, &c_b164,
01056 &a[j + (j + 1) * a_dim1], lda);
01057 i__2 = *n - j;
01058 r__1 = 1.f / ajj;
01059 sscal_(&i__2, &r__1, &a[j + (j + 1) * a_dim1], lda);
01060 }
01061
01062 }
01063 } else {
01064
01065
01066
01067 i__1 = *n;
01068 for (j = 1; j <= i__1; ++j) {
01069
01070
01071
01072 i__2 = j - 1;
01073 ajj = a[j + j * a_dim1] - sdot_(&i__2, &a[j + a_dim1], lda, &a[j
01074 + a_dim1], lda);
01075 if (ajj <= 0.f) {
01076 a[j + j * a_dim1] = ajj;
01077 goto L30;
01078 }
01079 ajj = sqrt(ajj);
01080 a[j + j * a_dim1] = ajj;
01081
01082
01083
01084 if (j < *n) {
01085 i__2 = *n - j;
01086 i__3 = j - 1;
01087 sgemv_("No transpose", &i__2, &i__3, &c_b181, &a[j + 1 +
01088 a_dim1], lda, &a[j + a_dim1], lda, &c_b164, &a[j + 1
01089 + j * a_dim1], &c__1);
01090 i__2 = *n - j;
01091 r__1 = 1.f / ajj;
01092 sscal_(&i__2, &r__1, &a[j + 1 + j * a_dim1], &c__1);
01093 }
01094
01095 }
01096 }
01097 goto L40;
01098
01099 L30:
01100 *info = j;
01101
01102 L40:
01103 return 0;
01104
01105
01106
01107 }
01108
01109 int spotrf_(char *uplo, integer *n, real *a, integer *lda,
01110 integer *info)
01111 {
01112
01113 integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
01114
01115
01116 static integer j, jb, nb;
01117 extern logical lsame_(char *, char *);
01118 extern int sgemm_(char *, char *, integer *, integer *,
01119 integer *, real *, real *, integer *, real *, integer *, real *,
01120 real *, integer *);
01121 static logical upper;
01122 extern int strsm_(char *, char *, char *, char *,
01123 integer *, integer *, real *, real *, integer *, real *, integer *
01124 ), ssyrk_(char *, char *, integer
01125 *, integer *, real *, real *, integer *, real *, real *, integer *
01126 ), spotf2_(char *, integer *, real *, integer *,
01127 integer *), xerbla_(char *, integer *);
01128 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
01129 integer *, integer *, ftnlen, ftnlen);
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191 a_dim1 = *lda;
01192 a_offset = 1 + a_dim1;
01193 a -= a_offset;
01194
01195
01196 *info = 0;
01197 upper = lsame_(uplo, "U");
01198 if (! upper && ! lsame_(uplo, "L")) {
01199 *info = -1;
01200 } else if (*n < 0) {
01201 *info = -2;
01202 } else if (*lda < max(1,*n)) {
01203 *info = -4;
01204 }
01205 if (*info != 0) {
01206 i__1 = -(*info);
01207 xerbla_("SPOTRF", &i__1);
01208 return 0;
01209 }
01210
01211
01212
01213 if (*n == 0) {
01214 return 0;
01215 }
01216
01217
01218
01219 nb = ilaenv_(&c__1, "SPOTRF", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
01220 ftnlen)1);
01221 if (nb <= 1 || nb >= *n) {
01222
01223
01224
01225 spotf2_(uplo, n, &a[a_offset], lda, info);
01226 } else {
01227
01228
01229
01230 if (upper) {
01231
01232
01233
01234 i__1 = *n;
01235 i__2 = nb;
01236 for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
01237
01238
01239
01240
01241
01242
01243
01244 i__3 = nb, i__4 = *n - j + 1;
01245 jb = min(i__3,i__4);
01246 i__3 = j - 1;
01247 ssyrk_("Upper", "Transpose", &jb, &i__3, &c_b181, &a[j *
01248 a_dim1 + 1], lda, &c_b164, &a[j + j * a_dim1], lda);
01249 spotf2_("Upper", &jb, &a[j + j * a_dim1], lda, info);
01250 if (*info != 0) {
01251 goto L30;
01252 }
01253 if (j + jb <= *n) {
01254
01255
01256
01257 i__3 = *n - j - jb + 1;
01258 i__4 = j - 1;
01259 sgemm_("Transpose", "No transpose", &jb, &i__3, &i__4, &
01260 c_b181, &a[j * a_dim1 + 1], lda, &a[(j + jb) *
01261 a_dim1 + 1], lda, &c_b164, &a[j + (j + jb) *
01262 a_dim1], lda);
01263 i__3 = *n - j - jb + 1;
01264 strsm_("Left", "Upper", "Transpose", "Non-unit", &jb, &
01265 i__3, &c_b164, &a[j + j * a_dim1], lda, &a[j + (j
01266 + jb) * a_dim1], lda);
01267 }
01268
01269 }
01270
01271 } else {
01272
01273
01274
01275 i__2 = *n;
01276 i__1 = nb;
01277 for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
01278
01279
01280
01281
01282
01283
01284
01285 i__3 = nb, i__4 = *n - j + 1;
01286 jb = min(i__3,i__4);
01287 i__3 = j - 1;
01288 ssyrk_("Lower", "No transpose", &jb, &i__3, &c_b181, &a[j +
01289 a_dim1], lda, &c_b164, &a[j + j * a_dim1], lda);
01290 spotf2_("Lower", &jb, &a[j + j * a_dim1], lda, info);
01291 if (*info != 0) {
01292 goto L30;
01293 }
01294 if (j + jb <= *n) {
01295
01296
01297
01298 i__3 = *n - j - jb + 1;
01299 i__4 = j - 1;
01300 sgemm_("No transpose", "Transpose", &i__3, &jb, &i__4, &
01301 c_b181, &a[j + jb + a_dim1], lda, &a[j + a_dim1],
01302 lda, &c_b164, &a[j + jb + j * a_dim1], lda);
01303 i__3 = *n - j - jb + 1;
01304 strsm_("Right", "Lower", "Transpose", "Non-unit", &i__3, &
01305 jb, &c_b164, &a[j + j * a_dim1], lda, &a[j + jb +
01306 j * a_dim1], lda);
01307 }
01308
01309 }
01310 }
01311 }
01312 goto L40;
01313
01314 L30:
01315 *info = *info + j - 1;
01316
01317 L40:
01318 return 0;
01319
01320
01321
01322 }
01323
01324 int spotrs_(char *uplo, integer *n, integer *nrhs, real *a,
01325 integer *lda, real *b, integer *ldb, integer *info)
01326 {
01327
01328 integer a_dim1, a_offset, b_dim1, b_offset, i__1;
01329
01330
01331 extern logical lsame_(char *, char *);
01332 static logical upper;
01333 extern int strsm_(char *, char *, char *, char *,
01334 integer *, integer *, real *, real *, integer *, real *, integer *
01335 ), xerbla_(char *, integer *);
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391 a_dim1 = *lda;
01392 a_offset = 1 + a_dim1;
01393 a -= a_offset;
01394 b_dim1 = *ldb;
01395 b_offset = 1 + b_dim1;
01396 b -= b_offset;
01397
01398
01399 *info = 0;
01400 upper = lsame_(uplo, "U");
01401 if (! upper && ! lsame_(uplo, "L")) {
01402 *info = -1;
01403 } else if (*n < 0) {
01404 *info = -2;
01405 } else if (*nrhs < 0) {
01406 *info = -3;
01407 } else if (*lda < max(1,*n)) {
01408 *info = -5;
01409 } else if (*ldb < max(1,*n)) {
01410 *info = -7;
01411 }
01412 if (*info != 0) {
01413 i__1 = -(*info);
01414 xerbla_("SPOTRS", &i__1);
01415 return 0;
01416 }
01417
01418
01419
01420 if (*n == 0 || *nrhs == 0) {
01421 return 0;
01422 }
01423
01424 if (upper) {
01425
01426
01427
01428
01429
01430
01431
01432 strsm_("Left", "Upper", "Transpose", "Non-unit", n, nrhs, &c_b164, &a[
01433 a_offset], lda, &b[b_offset], ldb);
01434
01435
01436
01437 strsm_("Left", "Upper", "No transpose", "Non-unit", n, nrhs, &c_b164,
01438 &a[a_offset], lda, &b[b_offset], ldb);
01439 } else {
01440
01441
01442
01443
01444
01445
01446
01447 strsm_("Left", "Lower", "No transpose", "Non-unit", n, nrhs, &c_b164,
01448 &a[a_offset], lda, &b[b_offset], ldb);
01449
01450
01451
01452 strsm_("Left", "Lower", "Transpose", "Non-unit", n, nrhs, &c_b164, &a[
01453 a_offset], lda, &b[b_offset], ldb);
01454 }
01455
01456 return 0;
01457
01458
01459
01460 }
01461