37 s = d_1 * d_1 + d_2 * d_2;
38 *cr = (ars * brs + ais * bis) / s;
39 *ci = (ais * brs - ars * bis) / s;
111 a = 1.3333333333333333;
115 eps = (d_1 = c - 1.,
abs(d_1));
119 ret_val = eps *
abs(*x);
135 d_1 =
abs(*a), d_2 =
abs(*b);
141 d_2 =
abs(*a), d_3 =
abs(*b);
143 d_1 =
min(d_2,d_3) / p;
166 integer t_dim1, t_offset, z_dim1, z_offset, i_1, i_2;
228 t_offset = t_dim1 + 1;
232 z_offset = z_dim1 + 1;
246 for (i = 2; i <= i_1; ++i) {
250 if (t[i + t_dim1] != 0. || t[i - 1 + t_dim1 * 3] != 0.) {
256 e[i] = e[i - 1] * e[i] / t[i - 1 + t_dim1 * 3];
262 for (j = 1; j <= i_1; ++j) {
265 for (i = 2; i <= i_2; ++i) {
266 z[i + j * z_dim1] *= e[i];
275 *ierr = (*n << 1) + i;
284 integer a_dim1, a_offset, i_1, i_2;
354 a_offset = a_dim1 + 1;
373 for (i = 1; i <= i_1; ++i) {
374 f = a[i + j * a_dim1];
375 a[i + j * a_dim1] = a[i + m * a_dim1];
376 a[i + m * a_dim1] = f;
381 for (i = k; i <= i_1; ++i) {
382 f = a[j + i * a_dim1];
383 a[j + i * a_dim1] = a[m + i * a_dim1];
384 a[m + i * a_dim1] = f;
403 for (jj = 1; jj <= i_1; ++jj) {
407 for (i = 1; i <= i_2; ++i) {
411 if (a[j + i * a_dim1] != 0.) {
433 for (j = k; j <= i_1; ++j) {
436 for (i = k; i <= i_2; ++i) {
440 if (a[i + j * a_dim1] != 0.) {
455 for (i = k; i <= i_1; ++i) {
464 for (i = k; i <= i_1; ++i) {
469 for (j = k; j <= i_2; ++j) {
473 c += (d_1 = a[j + i * a_dim1],
abs(d_1));
474 r += (d_1 = a[i + j * a_dim1],
abs(d_1));
480 if (c == 0. || r == 0.) {
504 if ((c + r) / f >= s * .95) {
512 for (j = k; j <= i_2; ++j) {
514 a[i + j * a_dim1] *= g;
518 for (j = 1; j <= i_2; ++j) {
520 a[j + i * a_dim1] *= f;
541 integer z_dim1, z_offset, i_1, i_2;
593 z_offset = z_dim1 + 1;
605 for (i = *low; i <= i_1; ++i) {
611 for (j = 1; j <= i_2; ++j) {
613 z[i + j * z_dim1] *=
s;
622 for (ii = 1; ii <= i_1; ++ii) {
624 if (i >= *low && i <= *igh) {
636 for (j = 1; j <= i_2; ++j) {
637 s = z[i + j * z_dim1];
638 z[i + j * z_dim1] = z[k + j * z_dim1];
639 z[k + j * z_dim1] =
s;
656 integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2, i_3, i_4, i_5,
669 static integer i1, i2, j1, j2, m1, n2, r1;
738 z_offset = z_dim1 + 1;
744 a_offset = a_dim1 + 1;
748 dmin_ = 5.4210108624275222e-20;
749 dminrt = 2.3283064365386963e-10;
752 for (j = 1; j <= i_1; ++j) {
762 for (j = 1; j <= i_1; ++j) {
765 for (k = 1; k <= i_2; ++k) {
767 z[j + k * z_dim1] = 0.;
770 z[j + j * z_dim1] = 1.;
776 if ((i_1 = m1 - 1) < 0) {
778 }
else if (i_1 == 0) {
787 for (k = 1; k <= i_1; ++k) {
789 i_2 = m1, i_3 = *n - k;
793 for (r1 = 2; r1 <= i_2; ++r1) {
797 g = a[kr + mr * a_dim1];
798 a[kr - 1 + a_dim1] = a[kr - 1 + (mr + 1) * a_dim1];
803 for (j = kr; i_4 < 0 ? j >= i_3 : j <= i_3; j += i_4) {
809 b1 = a[j1 + a_dim1] / g;
810 b2 = b1 * d[j1] / d[j];
811 s2 = 1. / (b1 * b2 + 1.);
815 b1 = g / a[j1 + a_dim1];
816 b2 = b1 * d[j] / d[j1];
820 f1 = a[j + m1 * a_dim1] * 2.;
821 f2 = b1 * a[j1 + *mb * a_dim1];
822 a[j + m1 * a_dim1] = -b2 * (b1 * a[j + m1 * a_dim1] - a[j + *
823 mb * a_dim1]) - f2 + a[j + m1 * a_dim1];
824 a[j1 + *mb * a_dim1] = b2 * (b2 * a[j + *mb * a_dim1] + f1) +
825 a[j1 + *mb * a_dim1];
826 a[j + *mb * a_dim1] = b1 * (f2 - f1) + a[j + *mb * a_dim1];
829 for (l = ugl; l <= i_5; ++l) {
831 u = a[j1 + (i2 + 1) * a_dim1] + b2 * a[j + i2 * a_dim1];
832 a[j + i2 * a_dim1] = -b1 * a[j1 + (i2 + 1) * a_dim1] + a[
834 a[j1 + (i2 + 1) * a_dim1] = u;
839 a[j1 + a_dim1] += b2 * g;
844 i_5 = m1, i_6 = *n - j1;
848 for (l = 2; l <= i_5; ++l) {
851 u = a[i1 + i2 * a_dim1] + b2 * a[i1 + (i2 + 1) * a_dim1];
852 a[i1 + (i2 + 1) * a_dim1] = -b1 * a[i1 + i2 * a_dim1] + a[
853 i1 + (i2 + 1) * a_dim1];
854 a[i1 + i2 * a_dim1] = u;
862 g = b2 * a[i1 + a_dim1];
869 for (l = 1; l <= i_5; ++l) {
870 u = z[l + j1 * z_dim1] + b2 * z[l + j * z_dim1];
871 z[l + j * z_dim1] = -b1 * z[l + j1 * z_dim1] + z[l + j *
873 z[l + j1 * z_dim1] = u;
883 f1 = a[j + m1 * a_dim1] * 2.;
884 f2 = b1 * a[j + *mb * a_dim1];
885 u = b1 * (f2 - f1) + a[j1 + *mb * a_dim1];
886 a[j + m1 * a_dim1] = b2 * (b1 * a[j + m1 * a_dim1] - a[j1 + *
887 mb * a_dim1]) + f2 - a[j + m1 * a_dim1];
888 a[j1 + *mb * a_dim1] = b2 * (b2 * a[j1 + *mb * a_dim1] + f1)
889 + a[j + *mb * a_dim1];
890 a[j + *mb * a_dim1] = u;
893 for (l = ugl; l <= i_5; ++l) {
895 u = b2 * a[j1 + (i2 + 1) * a_dim1] + a[j + i2 * a_dim1];
896 a[j + i2 * a_dim1] = -a[j1 + (i2 + 1) * a_dim1] + b1 * a[
898 a[j1 + (i2 + 1) * a_dim1] = u;
903 a[j1 + a_dim1] = b2 * a[j1 + a_dim1] + g;
908 i_5 = m1, i_6 = *n - j1;
912 for (l = 2; l <= i_5; ++l) {
915 u = b2 * a[i1 + i2 * a_dim1] + a[i1 + (i2 + 1) * a_dim1];
916 a[i1 + (i2 + 1) * a_dim1] = -a[i1 + i2 * a_dim1] + b1 * a[
917 i1 + (i2 + 1) * a_dim1];
918 a[i1 + i2 * a_dim1] = u;
927 a[i1 + a_dim1] = b1 * a[i1 + a_dim1];
934 for (l = 1; l <= i_5; ++l) {
935 u = b2 * z[l + j1 * z_dim1] + z[l + j * z_dim1];
936 z[l + j * z_dim1] = -z[l + j1 * z_dim1] + b1 * z[l + j *
938 z[l + j1 * z_dim1] = u;
955 for (j = k; j <= i_2; ++j) {
960 i_4 = 1, i_3 = *mb + 1 - j;
964 for (l = maxl; l <= i_4; ++l) {
966 a[j + l * a_dim1] = dminrt * a[j + l * a_dim1];
973 i_4 = m1, i_3 = *n - j;
977 for (l = 1; l <= i_4; ++l) {
980 a[i1 + i2 * a_dim1] = dminrt * a[i1 + i2 * a_dim1];
990 for (l = 1; l <= i_4; ++l) {
992 z[l + j * z_dim1] = dminrt * z[l + j * z_dim1];
996 a[j + *mb * a_dim1] = dmin_ * a[j + *mb * a_dim1];
1008 for (j = 2; j <= i_1; ++j) {
1018 for (j = 1; j <= i_1; ++j) {
1021 for (k = 2; k <= i_2; ++k) {
1023 z[j + k * z_dim1] = e[k] * z[j + k * z_dim1];
1033 for (j = 2; j <= i_1; ++j) {
1034 a[j + m1 * a_dim1] = u * e[j] * a[j + m1 * a_dim1];
1037 d_1 = a[j + m1 * a_dim1];
1039 a[j + *mb * a_dim1] = d[j] * a[j + *mb * a_dim1];
1040 d[j] = a[j + *mb * a_dim1];
1041 e[j] = a[j + m1 * a_dim1];
1045 d[1] = a[*mb * a_dim1 + 1];
1052 for (j = 1; j <= i_1; ++j) {
1053 d[j] = a[j + *mb * a_dim1];
1068 integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2, i_3, i_4, i_5;
1081 static integer mb, m21, ii, ij, jj, kj;
1192 a_offset = a_dim1 + 1;
1195 z_offset = z_dim1 + 1;
1207 mb = (*mbw + 1) / 2;
1211 order = 1. -
abs(*e21);
1214 for (r = 1; r <= i_1; ++r) {
1224 for (j = 1; j <= i_2; ++j) {
1231 for (i = jj; i <= i_3; ++i) {
1232 v += (d_1 = a[i + j * a_dim1],
abs(d_1));
1236 v += (d_1 = a[ij + kj * a_dim1],
abs(d_1));
1257 eps2 = norm * .001 *
abs(order);
1267 if ((d_1 = x1 - x0,
abs(d_1)) >= eps2) {
1271 if (order * (x1 - x0) <= 0.) {
1272 x1 = x0 + order * eps3;
1278 for (i = 1; i <= i_2; ++i) {
1280 i_3 = 0, i_4 = i - m1;
1281 ij = i +
min(i_3,i_4) * *n;
1289 for (j = 1; j <= i_3; ++j) {
1300 rv[ij] = a[i + j * a_dim1];
1314 rv[kj] = a[ii + jj * a_dim1];
1321 rv[ij] = a[i + mb * a_dim1] - x1;
1324 rv6[i] = z[i + r * z_dim1];
1334 for (i = 1; i <= i_2; ++i) {
1340 i_3 = *n - i, i_4 = m21 - 2;
1341 maxj =
min(i_3,i_4) * *n;
1344 for (k = i; k <= i_3; ++k) {
1351 for (kj = j; i_5 < 0 ? kj >= i_4 : kj <= i_4; kj += i_5) {
1369 i_3 = *n - ii, i_5 = m21 - 2;
1370 maxj =
min(i_3,i_5) * *n;
1373 for (j = i; j <= i_3; ++j) {
1374 if ((d_1 = rv[j],
abs(d_1)) <
abs(u)) {
1392 for (ij = i; i_5 < 0 ? ij >= i_3 : ij <= i_3; ij += i_5) {
1412 for (k = ii; k <= i_5; ++k) {
1418 for (ij = j; i_4 < 0 ? ij >= i_3 : ij <= i_3; ij += i_4) {
1420 rv[kj] -= v * rv[ij];
1425 rv6[k] -= v * rv6[i];
1437 for (ii = 1; ii <= i_2; ++ii) {
1445 jj = j + (maxj - 2) * *n;
1449 for (ij = j; i_4 < 0 ? ij >= i_5 : ij <= i_5; ij += i_4) {
1451 rv6[i] -= rv[ij] * rv6[ij1];
1457 if (
abs(v) >= eps3) {
1482 for (jj = 1; jj <= i_2; ++jj) {
1483 j = r - group - 1 + jj;
1487 for (i = 1; i <= i_4; ++i) {
1489 xu += rv6[i] * z[i + j * z_dim1];
1493 for (i = 1; i <= i_4; ++i) {
1495 rv6[i] -= xu * z[i + j * z_dim1];
1505 for (i = 1; i <= i_2; ++i) {
1507 norm += (d_1 = rv6[i],
abs(d_1));
1519 xu = eps4 / (uk + 1.);
1523 for (i = 2; i <= i_2; ++i) {
1528 rv6[its] -= eps4 * uk;
1541 for (i = 1; i <= i_2; ++i) {
1550 for (i = 1; i <= i_2; ++i) {
1552 z[i + r * z_dim1] = rv6[i] * xu;
1573 static integer i, j, k, l, p, q, r,
s;
1678 for (i = 1; i <= i_1; ++i) {
1682 tst1 = (d_1 = d[i],
abs(d_1)) + (d_2 = d[i - 1],
abs(d_2));
1683 tst2 = tst1 + (d_1 = e[i],
abs(d_1));
1724 for (q = p; q <= i_1; ++q) {
1731 u = (d_1 = e[q + 1],
abs(d_1));
1735 d_1 = d[q] - (x1 + u);
1738 d_1 = d[q] + (x1 + u);
1748 d_2 =
abs(xu), d_3 =
abs(x0);
1758 if (t1 > d[p] || d[p] >= t2) {
1768 d_1 = t1, d_2 = xu - x1;
1771 d_1 = t2, d_2 = x0 + x1;
1791 for (i = m1; i <= i_1; ++i) {
1805 for (ii = m1; ii <= i_1; ++ii) {
1822 x1 = (xu + x0) * .5;
1823 if (x0 - xu <=
abs(*eps1)) {
1826 tst1 = (
abs(xu) +
abs(x0)) * 2.;
1827 tst2 = tst1 + (x0 - xu);
1837 for (i = p; i <= i_1; ++i) {
1841 v = (d_1 = e[i],
abs(d_1)) /
epslon_(&c_b141);
1894 r = r + m2 - m1 + 1;
1899 for (l = 1; l <= i_1; ++l) {
1906 if (rv5[k] >= w[l]) {
1911 for (ii = j; ii <= i_2; ++ii) {
1914 ind[i + 1] = ind[i];
1949 integer a_dim1, a_offset, i_1, i_2, i_3;
1959 static integer imult, m1, m2, m3, m4, m21, m31, ii, ik, jk, kj, jm, kk,
1960 km, ll, mk, mn, ni, mz;
2059 a_offset = a_dim1 + 1;
2077 g = a[*n + *mb * a_dim1];
2084 for (k = 1; k <= i_1; ++k) {
2086 f += (d_1 = a[*n + mk * a_dim1],
abs(d_1));
2090 if (its == 0 && f > *r) {
2103 if (f > *r * .25 && its < 5) {
2106 f = a[*n + (*mb - 1) * a_dim1];
2110 q = (a[*n - 1 + *mb * a_dim1] - g) / (f * 2.);
2112 g -= f / (q +
d_sign(&s, &q));
2117 for (i = 1; i <= i_1; ++i) {
2119 a[i + *mb * a_dim1] -= g;
2124 for (k = m31; k <= i_1; ++k) {
2130 for (ii = 1; ii <= i_1; ++ii) {
2138 i_2 = 1, i_3 = 2 - i;
2142 for (k = 1; k <= i_2; ++k) {
2148 for (k = l; k <= i_2; ++k) {
2151 rv[km] = a[ii + mk * a_dim1];
2161 for (k = 1; k <= i_2; ++k) {
2165 rv[km] = a[ik + mk * a_dim1];
2178 for (j = 1; j <= i_2; ++j) {
2187 for (k = 1; k <= i_3; ++k) {
2190 f += rv[kj] * rv[jk];
2198 for (k = 1; k <= i_3; ++k) {
2201 rv[jk] -= rv[kj] * f;
2220 for (k = m21; k <= i_2; ++k) {
2222 scale += (d_1 = rv[k],
abs(d_1));
2230 for (k = m21; k <= i_2; ++k) {
2233 d_1 = rv[k] / scale;
2237 s = scale * scale *
s;
2242 kj = m4 + m2 * m1 + 1;
2246 for (k = 2; k <= i_2; ++k) {
2255 for (k = l; k <= i_2; ++k) {
2258 a[ii + mk * a_dim1] = rv[km];
2264 i_2 = 1, i_3 = m1 + 1 - i;
2271 for (k = 1; k <= i_2; ++k) {
2277 i_2 = m1, i_3 = ni + m1;
2281 for (kk = 1; kk <= i_2; ++kk) {
2286 rv[km] = a[ik + mk * a_dim1];
2297 for (k = l; k <= i_2; ++k) {
2299 a[i + mk * a_dim1] = rv[k];
2310 for (j = l; j <= i_2; ++j) {
2312 rv[jm] = rv[jm + 1];
2315 for (k = 1; k <= i_3; ++k) {
2332 for (i = 1; i <= i_1; ++i) {
2334 a[i + *mb * a_dim1] -= g;
2338 for (k = 1; k <= i_1; ++k) {
2340 a[*n + mk * a_dim1] = 0.;
2357 integer zr_dim1, zr_offset, zi_dim1, zi_offset, i_1, i_2;
2412 zi_offset = zi_dim1 + 1;
2415 zr_offset = zr_dim1 + 1;
2427 for (i = *low; i <= i_1; ++i) {
2433 for (j = 1; j <= i_2; ++j) {
2434 zr[i + j * zr_dim1] *=
s;
2435 zi[i + j * zi_dim1] *=
s;
2445 for (ii = 1; ii <= i_1; ++ii) {
2447 if (i >= *low && i <= *igh) {
2459 for (j = 1; j <= i_2; ++j) {
2460 s = zr[i + j * zr_dim1];
2461 zr[i + j * zr_dim1] = zr[k + j * zr_dim1];
2462 zr[k + j * zr_dim1] =
s;
2463 s = zi[i + j * zi_dim1];
2464 zi[i + j * zi_dim1] = zi[k + j * zi_dim1];
2465 zi[k + j * zi_dim1] =
s;
2481 integer ar_dim1, ar_offset, ai_dim1, ai_offset, i_1, i_2;
2556 ai_offset = ai_dim1 + 1;
2559 ar_offset = ar_dim1 + 1;
2578 for (i = 1; i <= i_1; ++i) {
2579 f = ar[i + j * ar_dim1];
2580 ar[i + j * ar_dim1] = ar[i + m * ar_dim1];
2581 ar[i + m * ar_dim1] = f;
2582 f = ai[i + j * ai_dim1];
2583 ai[i + j * ai_dim1] = ai[i + m * ai_dim1];
2584 ai[i + m * ai_dim1] = f;
2589 for (i = k; i <= i_1; ++i) {
2590 f = ar[j + i * ar_dim1];
2591 ar[j + i * ar_dim1] = ar[m + i * ar_dim1];
2592 ar[m + i * ar_dim1] = f;
2593 f = ai[j + i * ai_dim1];
2594 ai[j + i * ai_dim1] = ai[m + i * ai_dim1];
2595 ai[m + i * ai_dim1] = f;
2614 for (jj = 1; jj <= i_1; ++jj) {
2618 for (i = 1; i <= i_2; ++i) {
2622 if (ar[j + i * ar_dim1] != 0. || ai[j + i * ai_dim1] != 0.) {
2644 for (j = k; j <= i_1; ++j) {
2647 for (i = k; i <= i_2; ++i) {
2651 if (ar[i + j * ar_dim1] != 0. || ai[i + j * ai_dim1] != 0.) {
2666 for (i = k; i <= i_1; ++i) {
2675 for (i = k; i <= i_1; ++i) {
2680 for (j = k; j <= i_2; ++j) {
2684 c = c + (d_1 = ar[j + i * ar_dim1],
abs(d_1)) + (d_2 = ai[j +
2685 i * ai_dim1],
abs(d_2));
2686 r = r + (d_1 = ar[i + j * ar_dim1],
abs(d_1)) + (d_2 = ai[i +
2687 j * ai_dim1],
abs(d_2));
2693 if (c == 0. || r == 0.) {
2717 if ((c + r) / f >= s * .95) {
2725 for (j = k; j <= i_2; ++j) {
2726 ar[i + j * ar_dim1] *= g;
2727 ai[i + j * ai_dim1] *= g;
2732 for (j = 1; j <= i_2; ++j) {
2733 ar[j + i * ar_dim1] *= f;
2734 ai[j + i * ai_dim1] *= f;
2758 integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset,
2824 zi_offset = zi_dim1 + 1;
2827 zr_offset = zr_dim1 + 1;
2832 ai_offset = ai_dim1 + 1;
2835 ar_offset = ar_dim1 + 1;
2846 cbal_(nm, n, &ar[ar_offset], &ai[ai_offset], &is1, &is2, &fv1[1]);
2847 corth_(nm, n, &is1, &is2, &ar[ar_offset], &ai[ai_offset], &fv2[1], &fv3[1]
2853 comqr_(nm, n, &is1, &is2, &ar[ar_offset], &ai[ai_offset], &wr[1], &wi[1],
2858 comqr2_(nm, n, &is1, &is2, &fv2[1], &fv3[1], &ar[ar_offset], &ai[
2859 ai_offset], &wr[1], &wi[1], &zr[zr_offset], &zi[zi_offset], ierr);
2863 cbabk2_(nm, n, &is1, &is2, &fv1[1], n, &zr[zr_offset], &zi[zi_offset]);
2873 integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset,
2874 zi_dim1, zi_offset, i_1, i_2;
2935 zi_offset = zi_dim1 + 1;
2938 zr_offset = zr_dim1 + 1;
2942 ai_offset = ai_dim1 + 1;
2945 ar_offset = ar_dim1 + 1;
2956 htridi_(nm, n, &ar[ar_offset], &ai[ai_offset], &w[1], &fv1[1], &fv2[1], &
2962 tqlrat_(n, &w[1], &fv2[1], ierr);
2967 for (i = 1; i <= i_1; ++i) {
2970 for (j = 1; j <= i_2; ++j) {
2971 zr[j + i * zr_dim1] = 0.;
2975 zr[i + i * zr_dim1] = 1.;
2979 tql2_(nm, n, &w[1], &fv1[1], &zr[zr_offset], ierr);
2983 htribk_(nm, n, &ar[ar_offset], &ai[ai_offset], &fm1[3], n, &zr[zr_offset],
2996 integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset,
2997 zi_dim1, zi_offset, rm1_dim1, rm1_offset, rm2_dim1, rm2_offset,
3097 rm2_offset = rm2_dim1 + 1;
3100 rm1_offset = rm1_dim1 + 1;
3106 ai_offset = ai_dim1 + 1;
3109 ar_offset = ar_dim1 + 1;
3112 zi_offset = zi_dim1 + 1;
3115 zr_offset = zr_dim1 + 1;
3124 for (k = 1; k <= i_1; ++k) {
3136 for (uk = k; uk <= i_2; ++uk) {
3140 if (ar[uk + 1 + uk * ar_dim1] == 0. && ai[uk + 1 + uk * ai_dim1]
3153 for (i = 1; i <= i_2; ++i) {
3157 for (j = mp; j <= i_3; ++j) {
3159 x +=
pythag_(&ar[i + j * ar_dim1], &ai[i + j * ai_dim1]);
3176 ukroot = sqrt(ukroot);
3177 growto = .1 / ukroot;
3193 for (ii = 1; ii <= i_2; ++ii) {
3195 if (select[i] && (d_1 = wr[i] - rlambd,
abs(d_1)) < eps3 && (
3196 d_2 = wi[i] - ilambd,
abs(d_2)) < eps3) {
3209 for (i = 1; i <= i_2; ++i) {
3212 for (j = mp; j <= i_3; ++j) {
3213 rm1[i + j * rm1_dim1] = ar[i + j * ar_dim1];
3214 rm2[i + j * rm2_dim1] = ai[i + j * ai_dim1];
3218 rm1[i + i * rm1_dim1] -= rlambd;
3219 rm2[i + i * rm2_dim1] -= ilambd;
3231 for (i = 2; i <= i_2; ++i) {
3233 if (
pythag_(&rm1[i + mp * rm1_dim1], &rm2[i + mp * rm2_dim1]) <=
3234 pythag_(&rm1[mp + mp * rm1_dim1], &rm2[mp + mp * rm2_dim1]
3240 for (j = mp; j <= i_3; ++j) {
3241 y = rm1[i + j * rm1_dim1];
3242 rm1[i + j * rm1_dim1] = rm1[mp + j * rm1_dim1];
3243 rm1[mp + j * rm1_dim1] =
y;
3244 y = rm2[i + j * rm2_dim1];
3245 rm2[i + j * rm2_dim1] = rm2[mp + j * rm2_dim1];
3246 rm2[mp + j * rm2_dim1] =
y;
3251 if (rm1[mp + mp * rm1_dim1] == 0. && rm2[mp + mp * rm2_dim1] ==
3253 rm1[mp + mp * rm1_dim1] = eps3;
3255 cdiv_(&rm1[i + mp * rm1_dim1], &rm2[i + mp * rm2_dim1], &rm1[mp +
3256 mp * rm1_dim1], &rm2[mp + mp * rm2_dim1], &x, &y);
3257 if (x == 0. && y == 0.) {
3262 for (j = i; j <= i_3; ++j) {
3263 rm1[i + j * rm1_dim1] = rm1[i + j * rm1_dim1] - x * rm1[mp +
3264 j * rm1_dim1] + y * rm2[mp + j * rm2_dim1];
3265 rm2[i + j * rm2_dim1] = rm2[i + j * rm2_dim1] - x * rm2[mp +
3266 j * rm2_dim1] - y * rm1[mp + j * rm1_dim1];
3275 if (rm1[uk + uk * rm1_dim1] == 0. && rm2[uk + uk * rm2_dim1] == 0.) {
3276 rm1[uk + uk * rm1_dim1] = eps3;
3283 for (ii = 1; ii <= i_2; ++ii) {
3293 for (j = ip1; j <= i_3; ++j) {
3294 x = x - rm1[i + j * rm1_dim1] * rv1[j] + rm2[i + j * rm2_dim1]
3296 y = y - rm1[i + j * rm1_dim1] * rv2[j] - rm2[i + j * rm2_dim1]
3302 cdiv_(&x, &y, &rm1[i + i * rm1_dim1], &rm2[i + i * rm2_dim1], &
3313 for (i = 1; i <= i_2; ++i) {
3314 x =
pythag_(&rv1[i], &rv2[i]);
3325 if (norm < growto) {
3333 for (i = 1; i <= i_2; ++i) {
3334 cdiv_(&rv1[i], &rv2[i], &x, &y, &zr[i + s * zr_dim1], &zi[i + s *
3351 y = eps3 / (x + 1.);
3355 for (i = 2; i <= i_2; ++i) {
3371 for (i = j; i <= i_2; ++i) {
3372 zr[i + s * zr_dim1] = 0.;
3373 zi[i + s * zi_dim1] = 0.;
3391 *ierr = -((*n << 1) + 1);
3403 integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset,
3404 zi_dim1, zi_offset, i_1, i_2, i_3;
3407 static integer i, j, la, mm, mp;
3463 ai_offset = ai_dim1 + 1;
3466 ar_offset = ar_dim1 + 1;
3469 zi_offset = zi_dim1 + 1;
3472 zr_offset = zr_dim1 + 1;
3486 for (mm = kp1; mm <= i_1; ++mm) {
3487 mp = *low + *igh - mm;
3491 for (i = mp1; i <= i_2; ++i) {
3492 xr = ar[i + (mp - 1) * ar_dim1];
3493 xi = ai[i + (mp - 1) * ai_dim1];
3494 if (xr == 0. && xi == 0.) {
3499 for (j = 1; j <= i_3; ++j) {
3500 zr[i + j * zr_dim1] = zr[i + j * zr_dim1] + xr * zr[mp + j *
3501 zr_dim1] - xi * zi[mp + j * zi_dim1];
3502 zi[i + j * zi_dim1] = zi[i + j * zi_dim1] + xr * zi[mp + j *
3503 zi_dim1] + xi * zr[mp + j * zr_dim1];
3517 for (j = 1; j <= i_2; ++j) {
3518 xr = zr[i + j * zr_dim1];
3519 zr[i + j * zr_dim1] = zr[mp + j * zr_dim1];
3520 zr[mp + j * zr_dim1] = xr;
3521 xi = zi[i + j * zi_dim1];
3522 zi[i + j * zi_dim1] = zi[mp + j * zi_dim1];
3523 zi[mp + j * zi_dim1] =
xi;
3539 integer ar_dim1, ar_offset, ai_dim1, ai_offset, i_1, i_2, i_3;
3600 ai_offset = ai_dim1 + 1;
3603 ar_offset = ar_dim1 + 1;
3615 for (m = kp1; m <= i_1; ++m) {
3622 for (j = m; j <= i_2; ++j) {
3623 if ((d_1 = ar[j + mm1 * ar_dim1],
abs(d_1)) + (d_2 = ai[j +
3624 mm1 * ai_dim1],
abs(d_2)) <=
abs(xr) +
abs(xi)) {
3627 xr = ar[j + mm1 * ar_dim1];
3628 xi = ai[j + mm1 * ai_dim1];
3641 for (j = mm1; j <= i_2; ++j) {
3642 yr = ar[i + j * ar_dim1];
3643 ar[i + j * ar_dim1] = ar[m + j * ar_dim1];
3644 ar[m + j * ar_dim1] = yr;
3645 yi = ai[i + j * ai_dim1];
3646 ai[i + j * ai_dim1] = ai[m + j * ai_dim1];
3647 ai[m + j * ai_dim1] = yi;
3652 for (j = 1; j <= i_2; ++j) {
3653 yr = ar[j + i * ar_dim1];
3654 ar[j + i * ar_dim1] = ar[j + m * ar_dim1];
3655 ar[j + m * ar_dim1] = yr;
3656 yi = ai[j + i * ai_dim1];
3657 ai[j + i * ai_dim1] = ai[j + m * ai_dim1];
3658 ai[j + m * ai_dim1] = yi;
3663 if (xr == 0. && xi == 0.) {
3669 for (i = mp1; i <= i_2; ++i) {
3670 yr = ar[i + mm1 * ar_dim1];
3671 yi = ai[i + mm1 * ai_dim1];
3672 if (yr == 0. && yi == 0.) {
3675 cdiv_(&yr, &yi, &xr, &xi, &yr, &yi);
3676 ar[i + mm1 * ar_dim1] = yr;
3677 ai[i + mm1 * ai_dim1] = yi;
3680 for (j = m; j <= i_3; ++j) {
3681 ar[i + j * ar_dim1] = ar[i + j * ar_dim1] - yr * ar[m + j *
3682 ar_dim1] + yi * ai[m + j * ai_dim1];
3683 ai[i + j * ai_dim1] = ai[i + j * ai_dim1] - yr * ai[m + j *
3684 ai_dim1] - yi * ar[m + j * ar_dim1];
3689 for (j = 1; j <= i_3; ++j) {
3690 ar[j + m * ar_dim1] = ar[j + m * ar_dim1] + yr * ar[j + i *
3691 ar_dim1] - yi * ai[j + i * ai_dim1];
3692 ai[j + m * ai_dim1] = ai[j + m * ai_dim1] + yr * ai[j + i *
3693 ai_dim1] + yi * ar[j + i * ar_dim1];
3714 integer hr_dim1, hr_offset, hi_dim1, hi_offset, i_1, i_2;
3720 static integer i, j, l, m, en, ll, mm;
3790 hi_offset = hi_dim1 + 1;
3793 hr_offset = hr_dim1 + 1;
3800 for (i = 1; i <= i_1; ++i) {
3801 if (i >= *low && i <= *igh) {
3804 wr[i] = hr[i + i * hr_dim1];
3805 wi[i] = hi[i + i * hi_dim1];
3825 for (ll = *low; ll <= i_1; ++ll) {
3830 tst1 = (d_1 = hr[l - 1 + (l - 1) * hr_dim1],
abs(d_1)) + (d_2 = hi[
3831 l - 1 + (l - 1) * hi_dim1],
abs(d_2)) + (d_3 = hr[l + l *
3832 hr_dim1],
abs(d_3)) + (d_4 = hi[l + l * hi_dim1],
abs(d_4))
3834 tst2 = tst1 + (d_1 = hr[l + (l - 1) * hr_dim1],
abs(d_1)) + (d_2 =
3835 hi[l + (l - 1) * hi_dim1],
abs(d_2));
3849 if (its == 10 || its == 20) {
3852 sr = hr[en + en * hr_dim1];
3853 si = hi[en + en * hi_dim1];
3854 xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1] - hi[enm1 + en *
3855 hi_dim1] * hi[en + enm1 * hi_dim1];
3856 xi = hr[enm1 + en * hr_dim1] * hi[en + enm1 * hi_dim1] + hi[enm1 + en *
3857 hi_dim1] * hr[en + enm1 * hr_dim1];
3858 if (xr == 0. && xi == 0.) {
3861 yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
3862 yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
3867 d_1 = d_2 * d_2 - d_3 * d_3 + xr;
3868 d_4 = yr * 2. * yi +
xi;
3869 csroot_(&d_1, &d_4, &zzr, &zzi);
3870 if (yr * zzr + yi * zzi >= 0.) {
3878 cdiv_(&xr, &xi, &d_1, &d_2, &xr, &xi);
3884 sr = (d_1 = hr[en + enm1 * hr_dim1],
abs(d_1)) + (d_2 = hr[enm1 + (en
3885 - 2) * hr_dim1],
abs(d_2));
3886 si = (d_1 = hi[en + enm1 * hi_dim1],
abs(d_1)) + (d_2 = hi[enm1 + (en
3887 - 2) * hi_dim1],
abs(d_2));
3891 for (i = *low; i <= i_1; ++i) {
3892 hr[i + i * hr_dim1] -= sr;
3893 hi[i + i * hi_dim1] -= si;
3903 xr = (d_1 = hr[enm1 + enm1 * hr_dim1],
abs(d_1)) + (d_2 = hi[enm1 +
3904 enm1 * hi_dim1],
abs(d_2));
3905 yr = (d_1 = hr[en + enm1 * hr_dim1],
abs(d_1)) + (d_2 = hi[en + enm1 *
3906 hi_dim1],
abs(d_2));
3907 zzr = (d_1 = hr[en + en * hr_dim1],
abs(d_1)) + (d_2 = hi[en + en *
3908 hi_dim1],
abs(d_2));
3911 for (mm = l; mm <= i_1; ++mm) {
3917 yr = (d_1 = hr[m + (m - 1) * hr_dim1],
abs(d_1)) + (d_2 = hi[m + (
3918 m - 1) * hi_dim1],
abs(d_2));
3921 xr = (d_1 = hr[m - 1 + (m - 1) * hr_dim1],
abs(d_1)) + (d_2 = hi[m
3922 - 1 + (m - 1) * hi_dim1],
abs(d_2));
3923 tst1 = zzr / yi * (zzr + xr +
xi);
3935 for (i = mp1; i <= i_1; ++i) {
3937 xr = hr[im1 + im1 * hr_dim1];
3938 xi = hi[im1 + im1 * hi_dim1];
3939 yr = hr[i + im1 * hr_dim1];
3940 yi = hi[i + im1 * hi_dim1];
3946 for (j = im1; j <= i_2; ++j) {
3947 zzr = hr[im1 + j * hr_dim1];
3948 hr[im1 + j * hr_dim1] = hr[i + j * hr_dim1];
3949 hr[i + j * hr_dim1] = zzr;
3950 zzi = hi[im1 + j * hi_dim1];
3951 hi[im1 + j * hi_dim1] = hi[i + j * hi_dim1];
3952 hi[i + j * hi_dim1] = zzi;
3956 cdiv_(&xr, &xi, &yr, &yi, &zzr, &zzi);
3960 cdiv_(&yr, &yi, &xr, &xi, &zzr, &zzi);
3963 hr[i + im1 * hr_dim1] = zzr;
3964 hi[i + im1 * hi_dim1] = zzi;
3967 for (j = i; j <= i_2; ++j) {
3968 hr[i + j * hr_dim1] = hr[i + j * hr_dim1] - zzr * hr[im1 + j *
3969 hr_dim1] + zzi * hi[im1 + j * hi_dim1];
3970 hi[i + j * hi_dim1] = hi[i + j * hi_dim1] - zzr * hi[im1 + j *
3971 hi_dim1] - zzi * hr[im1 + j * hr_dim1];
3979 for (j = mp1; j <= i_1; ++j) {
3980 xr = hr[j + (j - 1) * hr_dim1];
3981 xi = hi[j + (j - 1) * hi_dim1];
3982 hr[j + (j - 1) * hr_dim1] = 0.;
3983 hi[j + (j - 1) * hi_dim1] = 0.;
3991 for (i = l; i <= i_2; ++i) {
3992 zzr = hr[i + (j - 1) * hr_dim1];
3993 hr[i + (j - 1) * hr_dim1] = hr[i + j * hr_dim1];
3994 hr[i + j * hr_dim1] = zzr;
3995 zzi = hi[i + (j - 1) * hi_dim1];
3996 hi[i + (j - 1) * hi_dim1] = hi[i + j * hi_dim1];
3997 hi[i + j * hi_dim1] = zzi;
4003 for (i = l; i <= i_2; ++i) {
4004 hr[i + (j - 1) * hr_dim1] = hr[i + (j - 1) * hr_dim1] + xr * hr[i
4005 + j * hr_dim1] - xi * hi[i + j * hi_dim1];
4006 hi[i + (j - 1) * hi_dim1] = hi[i + (j - 1) * hi_dim1] + xr * hi[i
4007 + j * hi_dim1] + xi * hr[i + j * hr_dim1];
4017 wr[en] = hr[en + en * hr_dim1] + tr;
4018 wi[en] = hi[en + en * hi_dim1] + ti;
4034 integer hr_dim1, hr_offset, hi_dim1, hi_offset, zr_dim1, zr_offset,
4035 zi_dim1, zi_offset, i_1, i_2, i_3;
4043 static integer i, j, k, l, m, ii, en, jj, ll, mm, nn;
4048 static integer ip1, mp1, itn, its;
4126 zi_offset = zi_dim1 + 1;
4129 zr_offset = zr_dim1 + 1;
4134 hi_offset = hi_dim1 + 1;
4137 hr_offset = hr_dim1 + 1;
4145 for (i = 1; i <= i_1; ++i) {
4148 for (j = 1; j <= i_2; ++j) {
4149 zr[i + j * zr_dim1] = 0.;
4150 zi[i + j * zi_dim1] = 0.;
4152 zr[i + j * zr_dim1] = 1.;
4159 iend = *igh - *low - 1;
4165 for (ii = 1; ii <= i_2; ++ii) {
4170 for (k = ip1; k <= i_1; ++k) {
4171 zr[k + i * zr_dim1] = hr[k + (i - 1) * hr_dim1];
4172 zi[k + i * zi_dim1] = hi[k + (i - 1) * hi_dim1];
4182 for (k = i; k <= i_1; ++k) {
4183 zr[i + k * zr_dim1] = zr[j + k * zr_dim1];
4184 zi[i + k * zi_dim1] = zi[j + k * zi_dim1];
4185 zr[j + k * zr_dim1] = 0.;
4186 zi[j + k * zi_dim1] = 0.;
4190 zr[j + i * zr_dim1] = 1.;
4197 for (i = 1; i <= i_2; ++i) {
4198 if (i >= *low && i <= *igh) {
4201 wr[i] = hr[i + i * hr_dim1];
4202 wi[i] = hi[i + i * hi_dim1];
4222 for (ll = *low; ll <= i_2; ++ll) {
4227 tst1 = (d_1 = hr[l - 1 + (l - 1) * hr_dim1],
abs(d_1)) + (d_2 = hi[
4228 l - 1 + (l - 1) * hi_dim1],
abs(d_2)) + (d_3 = hr[l + l *
4229 hr_dim1],
abs(d_3)) + (d_4 = hi[l + l * hi_dim1],
abs(d_4))
4231 tst2 = tst1 + (d_1 = hr[l + (l - 1) * hr_dim1],
abs(d_1)) + (d_2 =
4232 hi[l + (l - 1) * hi_dim1],
abs(d_2));
4246 if (its == 10 || its == 20) {
4249 sr = hr[en + en * hr_dim1];
4250 si = hi[en + en * hi_dim1];
4251 xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1] - hi[enm1 + en *
4252 hi_dim1] * hi[en + enm1 * hi_dim1];
4253 xi = hr[enm1 + en * hr_dim1] * hi[en + enm1 * hi_dim1] + hi[enm1 + en *
4254 hi_dim1] * hr[en + enm1 * hr_dim1];
4255 if (xr == 0. && xi == 0.) {
4258 yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
4259 yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
4264 d_1 = d_2 * d_2 - d_3 * d_3 + xr;
4265 d_4 = yr * 2. * yi +
xi;
4266 csroot_(&d_1, &d_4, &zzr, &zzi);
4267 if (yr * zzr + yi * zzi >= 0.) {
4275 cdiv_(&xr, &xi, &d_1, &d_2, &xr, &xi);
4281 sr = (d_1 = hr[en + enm1 * hr_dim1],
abs(d_1)) + (d_2 = hr[enm1 + (en
4282 - 2) * hr_dim1],
abs(d_2));
4283 si = (d_1 = hi[en + enm1 * hi_dim1],
abs(d_1)) + (d_2 = hi[enm1 + (en
4284 - 2) * hi_dim1],
abs(d_2));
4288 for (i = *low; i <= i_2; ++i) {
4289 hr[i + i * hr_dim1] -= sr;
4290 hi[i + i * hi_dim1] -= si;
4300 xr = (d_1 = hr[enm1 + enm1 * hr_dim1],
abs(d_1)) + (d_2 = hi[enm1 +
4301 enm1 * hi_dim1],
abs(d_2));
4302 yr = (d_1 = hr[en + enm1 * hr_dim1],
abs(d_1)) + (d_2 = hi[en + enm1 *
4303 hi_dim1],
abs(d_2));
4304 zzr = (d_1 = hr[en + en * hr_dim1],
abs(d_1)) + (d_2 = hi[en + en *
4305 hi_dim1],
abs(d_2));
4308 for (mm = l; mm <= i_2; ++mm) {
4314 yr = (d_1 = hr[m + (m - 1) * hr_dim1],
abs(d_1)) + (d_2 = hi[m + (
4315 m - 1) * hi_dim1],
abs(d_2));
4318 xr = (d_1 = hr[m - 1 + (m - 1) * hr_dim1],
abs(d_1)) + (d_2 = hi[m
4319 - 1 + (m - 1) * hi_dim1],
abs(d_2));
4320 tst1 = zzr / yi * (zzr + xr +
xi);
4332 for (i = mp1; i <= i_2; ++i) {
4334 xr = hr[im1 + im1 * hr_dim1];
4335 xi = hi[im1 + im1 * hi_dim1];
4336 yr = hr[i + im1 * hr_dim1];
4337 yi = hi[i + im1 * hi_dim1];
4343 for (j = im1; j <= i_1; ++j) {
4344 zzr = hr[im1 + j * hr_dim1];
4345 hr[im1 + j * hr_dim1] = hr[i + j * hr_dim1];
4346 hr[i + j * hr_dim1] = zzr;
4347 zzi = hi[im1 + j * hi_dim1];
4348 hi[im1 + j * hi_dim1] = hi[i + j * hi_dim1];
4349 hi[i + j * hi_dim1] = zzi;
4353 cdiv_(&xr, &xi, &yr, &yi, &zzr, &zzi);
4357 cdiv_(&yr, &yi, &xr, &xi, &zzr, &zzi);
4360 hr[i + im1 * hr_dim1] = zzr;
4361 hi[i + im1 * hi_dim1] = zzi;
4364 for (j = i; j <= i_1; ++j) {
4365 hr[i + j * hr_dim1] = hr[i + j * hr_dim1] - zzr * hr[im1 + j *
4366 hr_dim1] + zzi * hi[im1 + j * hi_dim1];
4367 hi[i + j * hi_dim1] = hi[i + j * hi_dim1] - zzr * hi[im1 + j *
4368 hi_dim1] - zzi * hr[im1 + j * hr_dim1];
4376 for (j = mp1; j <= i_2; ++j) {
4377 xr = hr[j + (j - 1) * hr_dim1];
4378 xi = hi[j + (j - 1) * hi_dim1];
4379 hr[j + (j - 1) * hr_dim1] = 0.;
4380 hi[j + (j - 1) * hi_dim1] = 0.;
4388 for (i = 1; i <= i_1; ++i) {
4389 zzr = hr[i + (j - 1) * hr_dim1];
4390 hr[i + (j - 1) * hr_dim1] = hr[i + j * hr_dim1];
4391 hr[i + j * hr_dim1] = zzr;
4392 zzi = hi[i + (j - 1) * hi_dim1];
4393 hi[i + (j - 1) * hi_dim1] = hi[i + j * hi_dim1];
4394 hi[i + j * hi_dim1] = zzi;
4399 for (i = *low; i <= i_1; ++i) {
4400 zzr = zr[i + (j - 1) * zr_dim1];
4401 zr[i + (j - 1) * zr_dim1] = zr[i + j * zr_dim1];
4402 zr[i + j * zr_dim1] = zzr;
4403 zzi = zi[i + (j - 1) * zi_dim1];
4404 zi[i + (j - 1) * zi_dim1] = zi[i + j * zi_dim1];
4405 zi[i + j * zi_dim1] = zzi;
4411 for (i = 1; i <= i_1; ++i) {
4412 hr[i + (j - 1) * hr_dim1] = hr[i + (j - 1) * hr_dim1] + xr * hr[i
4413 + j * hr_dim1] - xi * hi[i + j * hi_dim1];
4414 hi[i + (j - 1) * hi_dim1] = hi[i + (j - 1) * hi_dim1] + xr * hi[i
4415 + j * hi_dim1] + xi * hr[i + j * hr_dim1];
4420 for (i = *low; i <= i_1; ++i) {
4421 zr[i + (j - 1) * zr_dim1] = zr[i + (j - 1) * zr_dim1] + xr * zr[i
4422 + j * zr_dim1] - xi * zi[i + j * zi_dim1];
4423 zi[i + (j - 1) * zi_dim1] = zi[i + (j - 1) * zi_dim1] + xr * zi[i
4424 + j * zi_dim1] + xi * zr[i + j * zr_dim1];
4434 hr[en + en * hr_dim1] += tr;
4435 wr[en] = hr[en + en * hr_dim1];
4436 hi[en + en * hi_dim1] += ti;
4437 wi[en] = hi[en + en * hi_dim1];
4446 for (i = 1; i <= i_2; ++i) {
4449 for (j = i; j <= i_1; ++j) {
4450 tr = (d_1 = hr[i + j * hr_dim1],
abs(d_1)) + (d_2 = hi[i + j *
4451 hi_dim1],
abs(d_2));
4459 hr[hr_dim1 + 1] = norm;
4460 if (*n == 1 || norm == 0.) {
4465 for (nn = 2; nn <= i_1; ++nn) {
4469 hr[en + en * hr_dim1] = 1.;
4470 hi[en + en * hi_dim1] = 0.;
4474 for (ii = 1; ii <= i_2; ++ii) {
4481 for (j = ip1; j <= i_3; ++j) {
4482 zzr = zzr + hr[i + j * hr_dim1] * hr[j + en * hr_dim1] - hi[i
4483 + j * hi_dim1] * hi[j + en * hi_dim1];
4484 zzi = zzi + hr[i + j * hr_dim1] * hi[j + en * hi_dim1] + hi[i
4485 + j * hi_dim1] * hr[j + en * hr_dim1];
4491 if (yr != 0. || yi != 0.) {
4503 cdiv_(&zzr, &zzi, &yr, &yi, &hr[i + en * hr_dim1], &hi[i + en *
4506 tr = (d_1 = hr[i + en * hr_dim1],
abs(d_1)) + (d_2 = hi[i + en
4507 * hi_dim1],
abs(d_2));
4512 tst2 = tst1 + 1. / tst1;
4517 for (j = i; j <= i_3; ++j) {
4518 hr[j + en * hr_dim1] /= tr;
4519 hi[j + en * hi_dim1] /= tr;
4533 for (i = 1; i <= i_1; ++i) {
4534 if (i >= *low && i <= *igh) {
4540 for (j = ip1; j <= i_2; ++j) {
4541 zr[i + j * zr_dim1] = hr[i + j * hr_dim1];
4542 zi[i + j * zi_dim1] = hi[i + j * hi_dim1];
4553 for (jj = *low; jj <= i_1; ++jj) {
4558 for (i = *low; i <= i_2; ++i) {
4563 for (k = *low; k <= i_3; ++k) {
4564 zzr = zzr + zr[i + k * zr_dim1] * hr[k + j * hr_dim1] - zi[i
4565 + k * zi_dim1] * hi[k + j * hi_dim1];
4566 zzi = zzi + zr[i + k * zr_dim1] * hi[k + j * hi_dim1] + zi[i
4567 + k * zi_dim1] * hr[k + j * hr_dim1];
4571 zr[i + j * zr_dim1] = zzr;
4572 zi[i + j * zi_dim1] = zzi;
4591 integer hr_dim1, hr_offset, hi_dim1, hi_offset, i_1, i_2;
4598 static integer i, j, l, en, ll;
4672 hi_offset = hi_dim1 + 1;
4675 hr_offset = hr_dim1 + 1;
4687 for (i = l; i <= i_1; ++i) {
4691 if (hi[i + (i - 1) * hi_dim1] == 0.) {
4694 norm =
pythag_(&hr[i + (i - 1) * hr_dim1], &hi[i + (i - 1) * hi_dim1])
4696 yr = hr[i + (i - 1) * hr_dim1] / norm;
4697 yi = hi[i + (i - 1) * hi_dim1] / norm;
4698 hr[i + (i - 1) * hr_dim1] = norm;
4699 hi[i + (i - 1) * hi_dim1] = 0.;
4702 for (j = i; j <= i_2; ++j) {
4703 si = yr * hi[i + j * hi_dim1] - yi * hr[i + j * hr_dim1];
4704 hr[i + j * hr_dim1] = yr * hr[i + j * hr_dim1] + yi * hi[i + j *
4706 hi[i + j * hi_dim1] = si;
4711 for (j = *low; j <= i_2; ++j) {
4712 si = yr * hi[j + i * hi_dim1] + yi * hr[j + i * hr_dim1];
4713 hr[j + i * hr_dim1] = yr * hr[j + i * hr_dim1] - yi * hi[j + i *
4715 hi[j + i * hi_dim1] = si;
4725 for (i = 1; i <= i_1; ++i) {
4726 if (i >= *low && i <= *igh) {
4729 wr[i] = hr[i + i * hr_dim1];
4730 wi[i] = hi[i + i * hi_dim1];
4750 for (ll = *low; ll <= i_1; ++ll) {
4755 tst1 = (d_1 = hr[l - 1 + (l - 1) * hr_dim1],
abs(d_1)) + (d_2 = hi[
4756 l - 1 + (l - 1) * hi_dim1],
abs(d_2)) + (d_3 = hr[l + l *
4757 hr_dim1],
abs(d_3)) + (d_4 = hi[l + l * hi_dim1],
abs(d_4))
4759 tst2 = tst1 + (d_1 = hr[l + (l - 1) * hr_dim1],
abs(d_1));
4773 if (its == 10 || its == 20) {
4776 sr = hr[en + en * hr_dim1];
4777 si = hi[en + en * hi_dim1];
4778 xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1];
4779 xi = hi[enm1 + en * hi_dim1] * hr[en + enm1 * hr_dim1];
4780 if (xr == 0. && xi == 0.) {
4783 yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
4784 yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
4789 d_1 = d_2 * d_2 - d_3 * d_3 + xr;
4790 d_4 = yr * 2. * yi +
xi;
4791 csroot_(&d_1, &d_4, &zzr, &zzi);
4792 if (yr * zzr + yi * zzi >= 0.) {
4800 cdiv_(&xr, &xi, &d_1, &d_2, &xr, &xi);
4806 sr = (d_1 = hr[en + enm1 * hr_dim1],
abs(d_1)) + (d_2 = hr[enm1 + (en
4807 - 2) * hr_dim1],
abs(d_2));
4812 for (i = *low; i <= i_1; ++i) {
4813 hr[i + i * hr_dim1] -= sr;
4814 hi[i + i * hi_dim1] -= si;
4826 for (i = lp1; i <= i_1; ++i) {
4827 sr = hr[i + (i - 1) * hr_dim1];
4828 hr[i + (i - 1) * hr_dim1] = 0.;
4829 d_1 =
pythag_(&hr[i - 1 + (i - 1) * hr_dim1], &hi[i - 1 + (i - 1) *
4832 xr = hr[i - 1 + (i - 1) * hr_dim1] / norm;
4834 xi = hi[i - 1 + (i - 1) * hi_dim1] / norm;
4836 hr[i - 1 + (i - 1) * hr_dim1] = norm;
4837 hi[i - 1 + (i - 1) * hi_dim1] = 0.;
4838 hi[i + (i - 1) * hi_dim1] = sr / norm;
4841 for (j = i; j <= i_2; ++j) {
4842 yr = hr[i - 1 + j * hr_dim1];
4843 yi = hi[i - 1 + j * hi_dim1];
4844 zzr = hr[i + j * hr_dim1];
4845 zzi = hi[i + j * hi_dim1];
4846 hr[i - 1 + j * hr_dim1] = xr * yr + xi * yi + hi[i + (i - 1) *
4848 hi[i - 1 + j * hi_dim1] = xr * yi - xi * yr + hi[i + (i - 1) *
4850 hr[i + j * hr_dim1] = xr * zzr - xi * zzi - hi[i + (i - 1) *
4852 hi[i + j * hi_dim1] = xr * zzi + xi * zzr - hi[i + (i - 1) *
4860 si = hi[en + en * hi_dim1];
4864 norm =
pythag_(&hr[en + en * hr_dim1], &si);
4865 sr = hr[en + en * hr_dim1] / norm;
4867 hr[en + en * hr_dim1] = norm;
4868 hi[en + en * hi_dim1] = 0.;
4872 for (j = lp1; j <= i_1; ++j) {
4877 for (i = l; i <= i_2; ++i) {
4878 yr = hr[i + (j - 1) * hr_dim1];
4880 zzr = hr[i + j * hr_dim1];
4881 zzi = hi[i + j * hi_dim1];
4885 yi = hi[i + (j - 1) * hi_dim1];
4886 hi[i + (j - 1) * hi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) *
4889 hr[i + (j - 1) * hr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) *
4891 hr[i + j * hr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) *
4893 hi[i + j * hi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) *
4906 for (i = l; i <= i_1; ++i) {
4907 yr = hr[i + en * hr_dim1];
4908 yi = hi[i + en * hi_dim1];
4909 hr[i + en * hr_dim1] = sr * yr - si * yi;
4910 hi[i + en * hi_dim1] = sr * yi + si * yr;
4917 wr[en] = hr[en + en * hr_dim1] + tr;
4918 wi[en] = hi[en + en * hi_dim1] + ti;
4935 integer hr_dim1, hr_offset, hi_dim1, hi_offset, zr_dim1, zr_offset,
4936 zi_dim1, zi_offset, i_1, i_2, i_3;
4944 static integer i, j, k, l, m, ii, en, jj, ll, nn;
4949 static integer ip1, lp1, itn, its;
5032 zi_offset = zi_dim1 + 1;
5035 zr_offset = zr_dim1 + 1;
5040 hi_offset = hi_dim1 + 1;
5043 hr_offset = hr_dim1 + 1;
5052 for (j = 1; j <= i_1; ++j) {
5055 for (i = 1; i <= i_2; ++i) {
5056 zr[i + j * zr_dim1] = 0.;
5057 zi[i + j * zi_dim1] = 0.;
5060 zr[j + j * zr_dim1] = 1.;
5065 iend = *igh - *low - 1;
5068 }
else if (iend == 0) {
5076 for (ii = 1; ii <= i_1; ++ii) {
5078 if (ortr[i] == 0. && orti[i] == 0.) {
5081 if (hr[i + (i - 1) * hr_dim1] == 0. && hi[i + (i - 1) * hi_dim1] ==
5087 norm = hr[i + (i - 1) * hr_dim1] * ortr[i] + hi[i + (i - 1) * hi_dim1]
5092 for (k = ip1; k <= i_2; ++k) {
5093 ortr[k] = hr[k + (i - 1) * hr_dim1];
5094 orti[k] = hi[k + (i - 1) * hi_dim1];
5099 for (j = i; j <= i_2; ++j) {
5104 for (k = i; k <= i_3; ++k) {
5105 sr = sr + ortr[k] * zr[k + j * zr_dim1] + orti[k] * zi[k + j *
5107 si = si + ortr[k] * zi[k + j * zi_dim1] - orti[k] * zr[k + j *
5116 for (k = i; k <= i_3; ++k) {
5117 zr[k + j * zr_dim1] = zr[k + j * zr_dim1] + sr * ortr[k] - si
5119 zi[k + j * zi_dim1] = zi[k + j * zi_dim1] + sr * orti[k] + si
5135 for (i = l; i <= i_1; ++i) {
5139 if (hi[i + (i - 1) * hi_dim1] == 0.) {
5142 norm =
pythag_(&hr[i + (i - 1) * hr_dim1], &hi[i + (i - 1) * hi_dim1])
5144 yr = hr[i + (i - 1) * hr_dim1] / norm;
5145 yi = hi[i + (i - 1) * hi_dim1] / norm;
5146 hr[i + (i - 1) * hr_dim1] = norm;
5147 hi[i + (i - 1) * hi_dim1] = 0.;
5150 for (j = i; j <= i_2; ++j) {
5151 si = yr * hi[i + j * hi_dim1] - yi * hr[i + j * hr_dim1];
5152 hr[i + j * hr_dim1] = yr * hr[i + j * hr_dim1] + yi * hi[i + j *
5154 hi[i + j * hi_dim1] = si;
5159 for (j = 1; j <= i_2; ++j) {
5160 si = yr * hi[j + i * hi_dim1] + yi * hr[j + i * hr_dim1];
5161 hr[j + i * hr_dim1] = yr * hr[j + i * hr_dim1] - yi * hi[j + i *
5163 hi[j + i * hi_dim1] = si;
5168 for (j = *low; j <= i_2; ++j) {
5169 si = yr * zi[j + i * zi_dim1] + yi * zr[j + i * zr_dim1];
5170 zr[j + i * zr_dim1] = yr * zr[j + i * zr_dim1] - yi * zi[j + i *
5172 zi[j + i * zi_dim1] = si;
5182 for (i = 1; i <= i_1; ++i) {
5183 if (i >= *low && i <= *igh) {
5186 wr[i] = hr[i + i * hr_dim1];
5187 wi[i] = hi[i + i * hi_dim1];
5207 for (ll = *low; ll <= i_1; ++ll) {
5212 tst1 = (d_1 = hr[l - 1 + (l - 1) * hr_dim1],
abs(d_1)) + (d_2 = hi[
5213 l - 1 + (l - 1) * hi_dim1],
abs(d_2)) + (d_3 = hr[l + l *
5214 hr_dim1],
abs(d_3)) + (d_4 = hi[l + l * hi_dim1],
abs(d_4))
5216 tst2 = tst1 + (d_1 = hr[l + (l - 1) * hr_dim1],
abs(d_1));
5230 if (its == 10 || its == 20) {
5233 sr = hr[en + en * hr_dim1];
5234 si = hi[en + en * hi_dim1];
5235 xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1];
5236 xi = hi[enm1 + en * hi_dim1] * hr[en + enm1 * hr_dim1];
5237 if (xr == 0. && xi == 0.) {
5240 yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
5241 yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
5246 d_1 = d_2 * d_2 - d_3 * d_3 + xr;
5247 d_4 = yr * 2. * yi +
xi;
5248 csroot_(&d_1, &d_4, &zzr, &zzi);
5249 if (yr * zzr + yi * zzi >= 0.) {
5257 cdiv_(&xr, &xi, &d_1, &d_2, &xr, &xi);
5263 sr = (d_1 = hr[en + enm1 * hr_dim1],
abs(d_1)) + (d_2 = hr[enm1 + (en
5264 - 2) * hr_dim1],
abs(d_2));
5269 for (i = *low; i <= i_1; ++i) {
5270 hr[i + i * hr_dim1] -= sr;
5271 hi[i + i * hi_dim1] -= si;
5283 for (i = lp1; i <= i_1; ++i) {
5284 sr = hr[i + (i - 1) * hr_dim1];
5285 hr[i + (i - 1) * hr_dim1] = 0.;
5286 d_1 =
pythag_(&hr[i - 1 + (i - 1) * hr_dim1], &hi[i - 1 + (i - 1) *
5289 xr = hr[i - 1 + (i - 1) * hr_dim1] / norm;
5291 xi = hi[i - 1 + (i - 1) * hi_dim1] / norm;
5293 hr[i - 1 + (i - 1) * hr_dim1] = norm;
5294 hi[i - 1 + (i - 1) * hi_dim1] = 0.;
5295 hi[i + (i - 1) * hi_dim1] = sr / norm;
5298 for (j = i; j <= i_2; ++j) {
5299 yr = hr[i - 1 + j * hr_dim1];
5300 yi = hi[i - 1 + j * hi_dim1];
5301 zzr = hr[i + j * hr_dim1];
5302 zzi = hi[i + j * hi_dim1];
5303 hr[i - 1 + j * hr_dim1] = xr * yr + xi * yi + hi[i + (i - 1) *
5305 hi[i - 1 + j * hi_dim1] = xr * yi - xi * yr + hi[i + (i - 1) *
5307 hr[i + j * hr_dim1] = xr * zzr - xi * zzi - hi[i + (i - 1) *
5309 hi[i + j * hi_dim1] = xr * zzi + xi * zzr - hi[i + (i - 1) *
5317 si = hi[en + en * hi_dim1];
5321 norm =
pythag_(&hr[en + en * hr_dim1], &si);
5322 sr = hr[en + en * hr_dim1] / norm;
5324 hr[en + en * hr_dim1] = norm;
5325 hi[en + en * hi_dim1] = 0.;
5332 for (j = ip1; j <= i_1; ++j) {
5333 yr = hr[en + j * hr_dim1];
5334 yi = hi[en + j * hi_dim1];
5335 hr[en + j * hr_dim1] = sr * yr + si * yi;
5336 hi[en + j * hi_dim1] = sr * yi - si * yr;
5342 for (j = lp1; j <= i_1; ++j) {
5347 for (i = 1; i <= i_2; ++i) {
5348 yr = hr[i + (j - 1) * hr_dim1];
5350 zzr = hr[i + j * hr_dim1];
5351 zzi = hi[i + j * hi_dim1];
5355 yi = hi[i + (j - 1) * hi_dim1];
5356 hi[i + (j - 1) * hi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) *
5359 hr[i + (j - 1) * hr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) *
5361 hr[i + j * hr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) *
5363 hi[i + j * hi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) *
5369 for (i = *low; i <= i_2; ++i) {
5370 yr = zr[i + (j - 1) * zr_dim1];
5371 yi = zi[i + (j - 1) * zi_dim1];
5372 zzr = zr[i + j * zr_dim1];
5373 zzi = zi[i + j * zi_dim1];
5374 zr[i + (j - 1) * zr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) *
5376 zi[i + (j - 1) * zi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) *
5378 zr[i + j * zr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) *
5380 zi[i + j * zi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) *
5393 for (i = 1; i <= i_1; ++i) {
5394 yr = hr[i + en * hr_dim1];
5395 yi = hi[i + en * hi_dim1];
5396 hr[i + en * hr_dim1] = sr * yr - si * yi;
5397 hi[i + en * hi_dim1] = sr * yi + si * yr;
5402 for (i = *low; i <= i_1; ++i) {
5403 yr = zr[i + en * zr_dim1];
5404 yi = zi[i + en * zi_dim1];
5405 zr[i + en * zr_dim1] = sr * yr - si * yi;
5406 zi[i + en * zi_dim1] = sr * yi + si * yr;
5413 hr[en + en * hr_dim1] += tr;
5414 wr[en] = hr[en + en * hr_dim1];
5415 hi[en + en * hi_dim1] += ti;
5416 wi[en] = hi[en + en * hi_dim1];
5425 for (i = 1; i <= i_1; ++i) {
5428 for (j = i; j <= i_2; ++j) {
5429 tr = (d_1 = hr[i + j * hr_dim1],
abs(d_1)) + (d_2 = hi[i + j *
5430 hi_dim1],
abs(d_2));
5438 if (*n == 1 || norm == 0.) {
5443 for (nn = 2; nn <= i_2; ++nn) {
5447 hr[en + en * hr_dim1] = 1.;
5448 hi[en + en * hi_dim1] = 0.;
5452 for (ii = 1; ii <= i_1; ++ii) {
5459 for (j = ip1; j <= i_3; ++j) {
5460 zzr = zzr + hr[i + j * hr_dim1] * hr[j + en * hr_dim1] - hi[i
5461 + j * hi_dim1] * hi[j + en * hi_dim1];
5462 zzi = zzi + hr[i + j * hr_dim1] * hi[j + en * hi_dim1] + hi[i
5463 + j * hi_dim1] * hr[j + en * hr_dim1];
5469 if (yr != 0. || yi != 0.) {
5481 cdiv_(&zzr, &zzi, &yr, &yi, &hr[i + en * hr_dim1], &hi[i + en *
5484 tr = (d_1 = hr[i + en * hr_dim1],
abs(d_1)) + (d_2 = hi[i + en
5485 * hi_dim1],
abs(d_2));
5490 tst2 = tst1 + 1. / tst1;
5495 for (j = i; j <= i_3; ++j) {
5496 hr[j + en * hr_dim1] /= tr;
5497 hi[j + en * hi_dim1] /= tr;
5511 for (i = 1; i <= i_2; ++i) {
5512 if (i >= *low && i <= *igh) {
5518 for (j = ip1; j <= i_1; ++j) {
5519 zr[i + j * zr_dim1] = hr[i + j * hr_dim1];
5520 zi[i + j * zi_dim1] = hi[i + j * hi_dim1];
5531 for (jj = *low; jj <= i_2; ++jj) {
5536 for (i = *low; i <= i_1; ++i) {
5541 for (k = *low; k <= i_3; ++k) {
5542 zzr = zzr + zr[i + k * zr_dim1] * hr[k + j * hr_dim1] - zi[i
5543 + k * zi_dim1] * hi[k + j * hi_dim1];
5544 zzi = zzi + zr[i + k * zr_dim1] * hi[k + j * hi_dim1] + zi[i
5545 + k * zi_dim1] * hr[k + j * hr_dim1];
5549 zr[i + j * zr_dim1] = zzr;
5550 zi[i + j * zi_dim1] = zzi;
5569 integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset,
5570 zi_dim1, zi_offset, i_1, i_2, i_3;
5576 static integer mm, mp, kp1, mp1;
5637 ai_offset = ai_dim1 + 1;
5640 ar_offset = ar_dim1 + 1;
5643 zi_offset = zi_dim1 + 1;
5646 zr_offset = zr_dim1 + 1;
5660 for (mm = kp1; mm <= i_1; ++mm) {
5661 mp = *low + *igh - mm;
5662 if (ar[mp + (mp - 1) * ar_dim1] == 0. && ai[mp + (mp - 1) * ai_dim1]
5668 h = ar[mp + (mp - 1) * ar_dim1] * ortr[mp] + ai[mp + (mp - 1) *
5669 ai_dim1] * orti[mp];
5673 for (i = mp1; i <= i_2; ++i) {
5674 ortr[i] = ar[i + (mp - 1) * ar_dim1];
5675 orti[i] = ai[i + (mp - 1) * ai_dim1];
5680 for (j = 1; j <= i_2; ++j) {
5685 for (i = mp; i <= i_3; ++i) {
5686 gr = gr + ortr[i] * zr[i + j * zr_dim1] + orti[i] * zi[i + j *
5688 gi = gi + ortr[i] * zi[i + j * zi_dim1] - orti[i] * zr[i + j *
5697 for (i = mp; i <= i_3; ++i) {
5698 zr[i + j * zr_dim1] = zr[i + j * zr_dim1] + gr * ortr[i] - gi
5700 zi[i + j * zi_dim1] = zi[i + j * zi_dim1] + gr * orti[i] + gi
5721 integer ar_dim1, ar_offset, ai_dim1, ai_offset, i_1, i_2, i_3;
5790 ai_offset = ai_dim1 + 1;
5793 ar_offset = ar_dim1 + 1;
5806 for (m = kp1; m <= i_1; ++m) {
5814 for (i = m; i <= i_2; ++i) {
5816 scale = scale + (d_1 = ar[i + (m - 1) * ar_dim1],
abs(d_1)) + (
5817 d_2 = ai[i + (m - 1) * ai_dim1],
abs(d_2));
5826 for (ii = m; ii <= i_2; ++ii) {
5828 ortr[i] = ar[i + (m - 1) * ar_dim1] / scale;
5829 orti[i] = ai[i + (m - 1) * ai_dim1] / scale;
5830 h = h + ortr[i] * ortr[i] + orti[i] * orti[i];
5835 f =
pythag_(&ortr[m], &orti[m]);
5841 ortr[m] = (g + 1.) * ortr[m];
5842 orti[m] = (g + 1.) * orti[m];
5847 ar[m + (m - 1) * ar_dim1] = scale;
5851 for (j = m; j <= i_2; ++j) {
5856 for (ii = m; ii <= i_3; ++ii) {
5858 fr = fr + ortr[i] * ar[i + j * ar_dim1] + orti[i] * ai[i + j *
5860 fi = fi + ortr[i] * ai[i + j * ai_dim1] - orti[i] * ar[i + j *
5869 for (i = m; i <= i_3; ++i) {
5870 ar[i + j * ar_dim1] = ar[i + j * ar_dim1] - fr * ortr[i] + fi
5872 ai[i + j * ai_dim1] = ai[i + j * ai_dim1] - fr * orti[i] - fi
5881 for (i = 1; i <= i_2; ++i) {
5886 for (jj = m; jj <= i_3; ++jj) {
5888 fr = fr + ortr[j] * ar[i + j * ar_dim1] - orti[j] * ai[i + j *
5890 fi = fi + ortr[j] * ai[i + j * ai_dim1] + orti[j] * ar[i + j *
5899 for (j = m; j <= i_3; ++j) {
5900 ar[i + j * ar_dim1] = ar[i + j * ar_dim1] - fr * ortr[j] - fi
5902 ai[i + j * ai_dim1] = ai[i + j * ai_dim1] + fr * orti[j] - fi
5910 ortr[m] = scale * ortr[m];
5911 orti[m] = scale * orti[m];
5912 ar[m + (m - 1) * ar_dim1] = -g * ar[m + (m - 1) * ar_dim1];
5913 ai[m + (m - 1) * ai_dim1] = -g * ai[m + (m - 1) * ai_dim1];
5926 integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2, i_3;
5931 static integer la, mm, mp, kp1, mp1;
5983 a_offset = a_dim1 + 1;
5986 z_offset = z_dim1 + 1;
6000 for (mm = kp1; mm <= i_1; ++mm) {
6001 mp = *low + *igh - mm;
6005 for (i = mp1; i <= i_2; ++i) {
6006 x = a[i + (mp - 1) * a_dim1];
6012 for (j = 1; j <= i_3; ++j) {
6014 z[i + j * z_dim1] += x * z[mp + j * z_dim1];
6027 for (j = 1; j <= i_2; ++j) {
6028 x = z[i + j * z_dim1];
6029 z[i + j * z_dim1] = z[mp + j * z_dim1];
6030 z[mp + j * z_dim1] = x;
6046 integer a_dim1, a_offset, i_1, i_2, i_3;
6052 static integer la, mm1, kp1, mp1;
6100 a_offset = a_dim1 + 1;
6112 for (m = kp1; m <= i_1; ++m) {
6118 for (j = m; j <= i_2; ++j) {
6119 if ((d_1 = a[j + mm1 * a_dim1],
abs(d_1)) <=
abs(x)) {
6122 x = a[j + mm1 * a_dim1];
6134 for (j = mm1; j <= i_2; ++j) {
6135 y = a[i + j * a_dim1];
6136 a[i + j * a_dim1] = a[m + j * a_dim1];
6137 a[m + j * a_dim1] =
y;
6142 for (j = 1; j <= i_2; ++j) {
6143 y = a[j + i * a_dim1];
6144 a[j + i * a_dim1] = a[j + m * a_dim1];
6145 a[j + m * a_dim1] =
y;
6156 for (i = mp1; i <= i_2; ++i) {
6157 y = a[i + mm1 * a_dim1];
6162 a[i + mm1 * a_dim1] =
y;
6165 for (j = m; j <= i_3; ++j) {
6167 a[i + j * a_dim1] -= y * a[m + j * a_dim1];
6171 for (j = 1; j <= i_3; ++j) {
6173 a[j + m * a_dim1] += y * a[j + i * a_dim1];
6192 integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2;
6195 static integer i, j, kl, mm, mp, mp1;
6245 z_offset = z_dim1 + 1;
6249 a_offset = a_dim1 + 1;
6254 for (j = 1; j <= i_1; ++j) {
6257 for (i = 1; i <= i_2; ++i) {
6259 z[i + j * z_dim1] = 0.;
6262 z[j + j * z_dim1] = 1.;
6266 kl = *igh - *low - 1;
6272 for (mm = 1; mm <= i_1; ++mm) {
6277 for (i = mp1; i <= i_2; ++i) {
6279 z[i + mp * z_dim1] = a[i + (mp - 1) * a_dim1];
6288 for (j = mp; j <= i_2; ++j) {
6289 z[mp + j * z_dim1] = z[i + j * z_dim1];
6290 z[i + j * z_dim1] = 0.;
6294 z[i + mp * z_dim1] = 1.;
6307 integer t_dim1, t_offset, i_1;
6370 t_offset = t_dim1 + 1;
6380 for (i = 1; i <= i_1; ++i) {
6384 e2[i] = t[i + t_dim1] * t[i - 1 + t_dim1 * 3];
6385 if ((d_1 = e2[i]) < 0.) {
6387 }
else if (d_1 == 0) {
6393 if (t[i + t_dim1] == 0. && t[i - 1 + t_dim1 * 3] == 0.) {
6399 *ierr = -(*n * 3 + i);
6403 d[i] = t[i + (t_dim1 << 1)];
6420 integer t_dim1, t_offset, z_dim1, z_offset, i_1, i_2;
6480 t_offset = t_dim1 + 1;
6483 z_offset = z_dim1 + 1;
6492 for (i = 1; i <= i_1; ++i) {
6495 for (j = 1; j <= i_2; ++j) {
6497 z[i + j * z_dim1] = 0.;
6503 h = t[i + t_dim1] * t[i - 1 + t_dim1 * 3];
6506 }
else if (h == 0) {
6512 if (t[i + t_dim1] != 0. || t[i - 1 + t_dim1 * 3] != 0.) {
6517 z[i + i * z_dim1] = 1.;
6521 z[i + i * z_dim1] = z[i - 1 + (i - 1) * z_dim1] * e[i] / t[i - 1 +
6524 d[i] = t[i + (t_dim1 << 1)];
6537 *ierr = (*n << 1) + i;
6546 integer h_dim1, h_offset, i_1, i_2, i_3;
6556 static integer na, en, ll, mm;
6559 static integer mp2, itn, its, enm2;
6620 h_offset = h_dim1 + 1;
6630 for (i = 1; i <= i_1; ++i) {
6633 for (j = k; j <= i_2; ++j) {
6635 norm += (d_1 = h[i + j * h_dim1],
abs(d_1));
6639 if (i >= *low && i <= *igh) {
6642 wr[i] = h[i + i * h_dim1];
6663 for (ll = *low; ll <= i_1; ++ll) {
6668 s = (d_1 = h[l - 1 + (l - 1) * h_dim1],
abs(d_1)) + (d_2 = h[l + l
6669 * h_dim1],
abs(d_2));
6674 tst2 = tst1 + (d_1 = h[l + (l - 1) * h_dim1],
abs(d_1));
6682 x = h[en + en * h_dim1];
6686 y = h[na + na * h_dim1];
6687 w = h[en + na * h_dim1] * h[na + en * h_dim1];
6694 if (its != 10 && its != 20) {
6701 for (i = *low; i <= i_1; ++i) {
6703 h[i + i * h_dim1] -= x;
6706 s = (d_1 = h[en + na * h_dim1],
abs(d_1)) + (d_2 = h[na + enm2 *
6718 for (mm = l; mm <= i_1; ++mm) {
6720 zz = h[m + m * h_dim1];
6723 p = (r * s - w) / h[m + 1 + m * h_dim1] + h[m + (m + 1) * h_dim1];
6724 q = h[m + 1 + (m + 1) * h_dim1] - zz - r - s;
6725 r = h[m + 2 + (m + 1) * h_dim1];
6733 tst1 =
abs(p) * ((d_1 = h[m - 1 + (m - 1) * h_dim1],
abs(d_1)) +
6734 abs(zz) + (d_2 = h[m + 1 + (m + 1) * h_dim1],
abs(d_2)));
6735 tst2 = tst1 + (d_1 = h[m + (m - 1) * h_dim1],
abs(d_1)) * (
abs(q) +
6747 for (i = mp2; i <= i_1; ++i) {
6748 h[i + (i - 2) * h_dim1] = 0.;
6752 h[i + (i - 3) * h_dim1] = 0.;
6759 for (k = m; k <= i_1; ++k) {
6764 p = h[k + (k - 1) * h_dim1];
6765 q = h[k + 1 + (k - 1) * h_dim1];
6768 r = h[k + 2 + (k - 1) * h_dim1];
6778 d_1 = sqrt(p * p + q * q + r * r);
6783 h[k + (k - 1) * h_dim1] = -s * x;
6787 h[k + (k - 1) * h_dim1] = -h[k + (k - 1) * h_dim1];
6801 for (j = k; j <= i_2; ++j) {
6802 p = h[k + j * h_dim1] + q * h[k + 1 + j * h_dim1];
6803 h[k + j * h_dim1] -= p * x;
6804 h[k + 1 + j * h_dim1] -= p *
y;
6809 i_2 = en, i_3 = k + 3;
6813 for (i = 1; i <= i_2; ++i) {
6814 p = x * h[i + k * h_dim1] + y * h[i + (k + 1) * h_dim1];
6815 h[i + k * h_dim1] -= p;
6816 h[i + (k + 1) * h_dim1] -= p * q;
6823 for (j = k; j <= i_2; ++j) {
6824 p = h[k + j * h_dim1] + q * h[k + 1 + j * h_dim1] + r * h[k + 2 +
6826 h[k + j * h_dim1] -= p * x;
6827 h[k + 1 + j * h_dim1] -= p *
y;
6828 h[k + 2 + j * h_dim1] -= p * zz;
6833 i_2 = en, i_3 = k + 3;
6837 for (i = 1; i <= i_2; ++i) {
6838 p = x * h[i + k * h_dim1] + y * h[i + (k + 1) * h_dim1] + zz * h[
6839 i + (k + 2) * h_dim1];
6840 h[i + k * h_dim1] -= p;
6841 h[i + (k + 1) * h_dim1] -= p * q;
6842 h[i + (k + 2) * h_dim1] -= p * r;
6862 zz = sqrt((
abs(q)));
6868 zz = p +
d_sign(&zz, &p);
6872 wr[en] = x - w / zz;
6899 integer h_dim1, h_offset, z_dim1, z_offset, i_1, i_2, i_3;
6911 static integer na, ii, en, jj;
6916 static integer mp2, itn, its, enm2;
6992 z_offset = z_dim1 + 1;
6997 h_offset = h_dim1 + 1;
7007 for (i = 1; i <= i_1; ++i) {
7010 for (j = k; j <= i_2; ++j) {
7012 norm += (d_1 = h[i + j * h_dim1],
abs(d_1));
7016 if (i >= *low && i <= *igh) {
7019 wr[i] = h[i + i * h_dim1];
7040 for (ll = *low; ll <= i_1; ++ll) {
7045 s = (d_1 = h[l - 1 + (l - 1) * h_dim1],
abs(d_1)) + (d_2 = h[l + l
7046 * h_dim1],
abs(d_2));
7051 tst2 = tst1 + (d_1 = h[l + (l - 1) * h_dim1],
abs(d_1));
7059 x = h[en + en * h_dim1];
7063 y = h[na + na * h_dim1];
7064 w = h[en + na * h_dim1] * h[na + en * h_dim1];
7071 if (its != 10 && its != 20) {
7078 for (i = *low; i <= i_1; ++i) {
7080 h[i + i * h_dim1] -= x;
7083 s = (d_1 = h[en + na * h_dim1],
abs(d_1)) + (d_2 = h[na + enm2 *
7095 for (mm = l; mm <= i_1; ++mm) {
7097 zz = h[m + m * h_dim1];
7100 p = (r * s - w) / h[m + 1 + m * h_dim1] + h[m + (m + 1) * h_dim1];
7101 q = h[m + 1 + (m + 1) * h_dim1] - zz - r - s;
7102 r = h[m + 2 + (m + 1) * h_dim1];
7110 tst1 =
abs(p) * ((d_1 = h[m - 1 + (m - 1) * h_dim1],
abs(d_1)) +
7111 abs(zz) + (d_2 = h[m + 1 + (m + 1) * h_dim1],
abs(d_2)));
7112 tst2 = tst1 + (d_1 = h[m + (m - 1) * h_dim1],
abs(d_1)) * (
abs(q) +
7124 for (i = mp2; i <= i_1; ++i) {
7125 h[i + (i - 2) * h_dim1] = 0.;
7129 h[i + (i - 3) * h_dim1] = 0.;
7136 for (k = m; k <= i_1; ++k) {
7141 p = h[k + (k - 1) * h_dim1];
7142 q = h[k + 1 + (k - 1) * h_dim1];
7145 r = h[k + 2 + (k - 1) * h_dim1];
7155 d_1 = sqrt(p * p + q * q + r * r);
7160 h[k + (k - 1) * h_dim1] = -s * x;
7164 h[k + (k - 1) * h_dim1] = -h[k + (k - 1) * h_dim1];
7178 for (j = k; j <= i_2; ++j) {
7179 p = h[k + j * h_dim1] + q * h[k + 1 + j * h_dim1];
7180 h[k + j * h_dim1] -= p * x;
7181 h[k + 1 + j * h_dim1] -= p *
y;
7186 i_2 = en, i_3 = k + 3;
7190 for (i = 1; i <= i_2; ++i) {
7191 p = x * h[i + k * h_dim1] + y * h[i + (k + 1) * h_dim1];
7192 h[i + k * h_dim1] -= p;
7193 h[i + (k + 1) * h_dim1] -= p * q;
7198 for (i = *low; i <= i_2; ++i) {
7199 p = x * z[i + k * z_dim1] + y * z[i + (k + 1) * z_dim1];
7200 z[i + k * z_dim1] -= p;
7201 z[i + (k + 1) * z_dim1] -= p * q;
7208 for (j = k; j <= i_2; ++j) {
7209 p = h[k + j * h_dim1] + q * h[k + 1 + j * h_dim1] + r * h[k + 2 +
7211 h[k + j * h_dim1] -= p * x;
7212 h[k + 1 + j * h_dim1] -= p *
y;
7213 h[k + 2 + j * h_dim1] -= p * zz;
7218 i_2 = en, i_3 = k + 3;
7222 for (i = 1; i <= i_2; ++i) {
7223 p = x * h[i + k * h_dim1] + y * h[i + (k + 1) * h_dim1] + zz * h[
7224 i + (k + 2) * h_dim1];
7225 h[i + k * h_dim1] -= p;
7226 h[i + (k + 1) * h_dim1] -= p * q;
7227 h[i + (k + 2) * h_dim1] -= p * r;
7232 for (i = *low; i <= i_2; ++i) {
7233 p = x * z[i + k * z_dim1] + y * z[i + (k + 1) * z_dim1] + zz * z[
7234 i + (k + 2) * z_dim1];
7235 z[i + k * z_dim1] -= p;
7236 z[i + (k + 1) * z_dim1] -= p * q;
7237 z[i + (k + 2) * z_dim1] -= p * r;
7249 h[en + en * h_dim1] = x + t;
7250 wr[en] = h[en + en * h_dim1];
7258 zz = sqrt((
abs(q)));
7259 h[en + en * h_dim1] = x + t;
7260 x = h[en + en * h_dim1];
7261 h[na + na * h_dim1] = y + t;
7266 zz = p +
d_sign(&zz, &p);
7270 wr[en] = x - w / zz;
7274 x = h[en + na * h_dim1];
7278 r = sqrt(p * p + q * q);
7283 for (j = na; j <= i_1; ++j) {
7284 zz = h[na + j * h_dim1];
7285 h[na + j * h_dim1] = q * zz + p * h[en + j * h_dim1];
7286 h[en + j * h_dim1] = q * h[en + j * h_dim1] - p * zz;
7291 for (i = 1; i <= i_1; ++i) {
7292 zz = h[i + na * h_dim1];
7293 h[i + na * h_dim1] = q * zz + p * h[i + en * h_dim1];
7294 h[i + en * h_dim1] = q * h[i + en * h_dim1] - p * zz;
7299 for (i = *low; i <= i_1; ++i) {
7300 zz = z[i + na * z_dim1];
7301 z[i + na * z_dim1] = q * zz + p * z[i + en * z_dim1];
7302 z[i + en * z_dim1] = q * z[i + en * z_dim1] - p * zz;
7324 for (nn = 1; nn <= i_1; ++nn) {
7331 }
else if (q == 0) {
7339 h[en + en * h_dim1] = 1.;
7345 for (ii = 1; ii <= i_2; ++ii) {
7347 w = h[i + i * h_dim1] - p;
7351 for (j = m; j <= i_3; ++j) {
7353 r += h[i + j * h_dim1] * h[j + en * h_dim1];
7380 h[i + en * h_dim1] = -r / t;
7384 x = h[i + (i + 1) * h_dim1];
7385 y = h[i + 1 + i * h_dim1];
7386 q = (wr[i] - p) * (wr[i] - p) + wi[i] * wi[i];
7387 t = (x * s - zz * r) / q;
7388 h[i + en * h_dim1] = t;
7392 h[i + 1 + en * h_dim1] = (-r - w * t) / x;
7395 h[i + 1 + en * h_dim1] = (-s - y * t) / zz;
7399 t = (d_1 = h[i + en * h_dim1],
abs(d_1));
7404 tst2 = tst1 + 1. / tst1;
7409 for (j = i; j <= i_3; ++j) {
7410 h[j + en * h_dim1] /= t;
7424 if ((d_1 = h[en + na * h_dim1],
abs(d_1)) <= (d_2 = h[na + en *
7425 h_dim1],
abs(d_2))) {
7428 h[na + na * h_dim1] = q / h[en + na * h_dim1];
7429 h[na + en * h_dim1] = -(h[en + en * h_dim1] - p) / h[en + na * h_dim1]
7433 d_1 = -h[na + en * h_dim1];
7434 d_2 = h[na + na * h_dim1] - p;
7435 cdiv_(&c_b550, &d_1, &d_2, &q, &h[na + na * h_dim1], &h[na + en *
7438 h[en + na * h_dim1] = 0.;
7439 h[en + en * h_dim1] = 1.;
7446 for (ii = 1; ii <= i_2; ++ii) {
7448 w = h[i + i * h_dim1] - p;
7453 for (j = m; j <= i_3; ++j) {
7454 ra += h[i + j * h_dim1] * h[j + na * h_dim1];
7455 sa += h[i + j * h_dim1] * h[j + en * h_dim1];
7473 cdiv_(&d_1, &d_2, &w, &q, &h[i + na * h_dim1], &h[i + en *
7478 x = h[i + (i + 1) * h_dim1];
7479 y = h[i + 1 + i * h_dim1];
7480 vr = (wr[i] - p) * (wr[i] - p) + wi[i] * wi[i] - q * q;
7481 vi = (wr[i] - p) * 2. * q;
7482 if (vr != 0. || vi != 0.) {
7494 d_1 = x * r - zz * ra + q * sa;
7495 d_2 = x * s - zz * sa - q * ra;
7496 cdiv_(&d_1, &d_2, &vr, &vi, &h[i + na * h_dim1], &h[i + en *
7501 h[i + 1 + na * h_dim1] = (-ra - w * h[i + na * h_dim1] + q * h[i
7502 + en * h_dim1]) / x;
7503 h[i + 1 + en * h_dim1] = (-sa - w * h[i + en * h_dim1] - q * h[i
7504 + na * h_dim1]) / x;
7507 d_1 = -r - y * h[i + na * h_dim1];
7508 d_2 = -s - y * h[i + en * h_dim1];
7509 cdiv_(&d_1, &d_2, &zz, &q, &h[i + 1 + na * h_dim1], &h[i + 1 +
7515 d_3 = (d_1 = h[i + na * h_dim1],
abs(d_1)), d_4 = (d_2 = h[i
7516 + en * h_dim1],
abs(d_2));
7522 tst2 = tst1 + 1. / tst1;
7527 for (j = i; j <= i_3; ++j) {
7528 h[j + na * h_dim1] /= t;
7529 h[j + en * h_dim1] /= t;
7543 for (i = 1; i <= i_1; ++i) {
7544 if (i >= *low && i <= *igh) {
7549 for (j = i; j <= i_2; ++j) {
7551 z[i + j * z_dim1] = h[i + j * h_dim1];
7561 for (jj = *low; jj <= i_1; ++jj) {
7566 for (i = *low; i <= i_2; ++i) {
7570 for (k = *low; k <= i_3; ++k) {
7572 zz += z[i + k * z_dim1] * h[k + j * h_dim1];
7575 z[i + j * z_dim1] = zz;
7593 integer a_dim1, a_offset, zr_dim1, zr_offset, zi_dim1, zi_offset, i_1,
7651 a_offset = a_dim1 + 1;
7654 zi_offset = zi_dim1 + 1;
7657 zr_offset = zr_dim1 + 1;
7668 for (k = 1; k <= i_1; ++k) {
7671 for (j = 1; j <= i_2; ++j) {
7672 zi[k + j * zi_dim1] = -zr[k + j * zr_dim1] * tau[(k << 1) + 2];
7673 zr[k + j * zr_dim1] *= tau[(k << 1) + 1];
7683 for (i = 2; i <= i_2; ++i) {
7685 h = a[i + i * a_dim1];
7691 for (j = 1; j <= i_1; ++j) {
7696 for (k = 1; k <= i_3; ++k) {
7697 s = s + a[i + k * a_dim1] * zr[k + j * zr_dim1] - a[k + i *
7698 a_dim1] * zi[k + j * zi_dim1];
7699 si = si + a[i + k * a_dim1] * zi[k + j * zi_dim1] + a[k + i *
7700 a_dim1] * zr[k + j * zr_dim1];
7709 for (k = 1; k <= i_3; ++k) {
7710 zr[k + j * zr_dim1] = zr[k + j * zr_dim1] - s * a[i + k *
7711 a_dim1] - si * a[k + i * a_dim1];
7712 zi[k + j * zi_dim1] = zi[k + j * zi_dim1] - si * a[i + k *
7713 a_dim1] + s * a[k + i * a_dim1];
7733 integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset,
7734 zi_dim1, zi_offset, i_1, i_2, i_3;
7792 ai_offset = ai_dim1 + 1;
7795 ar_offset = ar_dim1 + 1;
7798 zi_offset = zi_dim1 + 1;
7801 zr_offset = zr_dim1 + 1;
7812 for (k = 1; k <= i_1; ++k) {
7815 for (j = 1; j <= i_2; ++j) {
7816 zi[k + j * zi_dim1] = -zr[k + j * zr_dim1] * tau[(k << 1) + 2];
7817 zr[k + j * zr_dim1] *= tau[(k << 1) + 1];
7827 for (i = 2; i <= i_2; ++i) {
7829 h = ai[i + i * ai_dim1];
7835 for (j = 1; j <= i_1; ++j) {
7840 for (k = 1; k <= i_3; ++k) {
7841 s = s + ar[i + k * ar_dim1] * zr[k + j * zr_dim1] - ai[i + k *
7842 ai_dim1] * zi[k + j * zi_dim1];
7843 si = si + ar[i + k * ar_dim1] * zi[k + j * zi_dim1] + ai[i +
7844 k * ai_dim1] * zr[k + j * zr_dim1];
7853 for (k = 1; k <= i_3; ++k) {
7854 zr[k + j * zr_dim1] = zr[k + j * zr_dim1] - s * ar[i + k *
7855 ar_dim1] - si * ai[i + k * ai_dim1];
7856 zi[k + j * zi_dim1] = zi[k + j * zi_dim1] - si * ar[i + k *
7857 ar_dim1] + s * ai[i + k * ai_dim1];
7876 integer a_dim1, a_offset, i_1, i_2, i_3;
7950 a_offset = a_dim1 + 1;
7954 tau[(*n << 1) + 1] = 1.;
7955 tau[(*n << 1) + 2] = 0.;
7958 for (ii = 1; ii <= i_1; ++ii) {
7968 for (k = 1; k <= i_2; ++k) {
7970 scale = scale + (d_1 = a[i + k * a_dim1],
abs(d_1)) + (d_2 = a[
7971 k + i * a_dim1],
abs(d_2));
7977 tau[(l << 1) + 1] = 1.;
7978 tau[(l << 1) + 2] = 0.;
7986 for (k = 1; k <= i_2; ++k) {
7987 a[i + k * a_dim1] /= scale;
7988 a[k + i * a_dim1] /= scale;
7989 h = h + a[i + k * a_dim1] * a[i + k * a_dim1] + a[k + i * a_dim1]
7990 * a[k + i * a_dim1];
7994 e2[i] = scale * scale * h;
7997 f =
pythag_(&a[i + l * a_dim1], &a[l + i * a_dim1]);
8002 tau[(l << 1) + 1] = (a[l + i * a_dim1] * tau[(i << 1) + 2] - a[i + l *
8003 a_dim1] * tau[(i << 1) + 1]) / f;
8004 si = (a[i + l * a_dim1] * tau[(i << 1) + 2] + a[l + i * a_dim1] * tau[
8008 a[i + l * a_dim1] = g * a[i + l * a_dim1];
8009 a[l + i * a_dim1] = g * a[l + i * a_dim1];
8015 tau[(l << 1) + 1] = -tau[(i << 1) + 1];
8016 si = tau[(i << 1) + 2];
8017 a[i + l * a_dim1] = g;
8022 for (j = 1; j <= i_2; ++j) {
8031 for (k = 1; k <= i_3; ++k) {
8032 g = g + a[j + k * a_dim1] * a[i + k * a_dim1] + a[k + j *
8033 a_dim1] * a[k + i * a_dim1];
8034 gi = gi - a[j + k * a_dim1] * a[k + i * a_dim1] + a[k + j *
8035 a_dim1] * a[i + k * a_dim1];
8040 g += a[j + j * a_dim1] * a[i + j * a_dim1];
8041 gi -= a[j + j * a_dim1] * a[j + i * a_dim1];
8048 for (k = jp1; k <= i_3; ++k) {
8049 g = g + a[k + j * a_dim1] * a[i + k * a_dim1] - a[j + k *
8050 a_dim1] * a[k + i * a_dim1];
8051 gi = gi - a[k + j * a_dim1] * a[k + i * a_dim1] - a[j + k *
8052 a_dim1] * a[i + k * a_dim1];
8058 tau[(j << 1) + 2] = gi / h;
8059 f = f + e[j] * a[i + j * a_dim1] - tau[(j << 1) + 2] * a[j + i *
8067 for (j = 1; j <= i_2; ++j) {
8068 f = a[i + j * a_dim1];
8071 fi = -a[j + i * a_dim1];
8072 gi = tau[(j << 1) + 2] - hh * fi;
8073 tau[(j << 1) + 2] = -gi;
8074 a[j + j * a_dim1] -= (f * g + fi * gi) * 2.;
8081 for (k = 1; k <= i_3; ++k) {
8082 a[j + k * a_dim1] = a[j + k * a_dim1] - f * e[k] - g * a[i +
8083 k * a_dim1] + fi * tau[(k << 1) + 2] + gi * a[k + i *
8085 a[k + j * a_dim1] = a[k + j * a_dim1] - f * tau[(k << 1) + 2]
8086 - g * a[k + i * a_dim1] - fi * e[k] - gi * a[i + k *
8097 for (k = 1; k <= i_2; ++k) {
8098 a[i + k * a_dim1] = scale * a[i + k * a_dim1];
8099 a[k + i * a_dim1] = scale * a[k + i * a_dim1];
8103 tau[(l << 1) + 2] = -si;
8105 d[i] = a[i + i * a_dim1];
8106 a[i + i * a_dim1] = scale * sqrt(h);
8118 integer ar_dim1, ar_offset, ai_dim1, ai_offset, i_1, i_2, i_3;
8191 ai_offset = ai_dim1 + 1;
8194 ar_offset = ar_dim1 + 1;
8198 tau[(*n << 1) + 1] = 1.;
8199 tau[(*n << 1) + 2] = 0.;
8202 for (i = 1; i <= i_1; ++i) {
8204 d[i] = ar[i + i * ar_dim1];
8208 for (ii = 1; ii <= i_1; ++ii) {
8218 for (k = 1; k <= i_2; ++k) {
8220 scale = scale + (d_1 = ar[i + k * ar_dim1],
abs(d_1)) + (d_2 =
8221 ai[i + k * ai_dim1],
abs(d_2));
8227 tau[(l << 1) + 1] = 1.;
8228 tau[(l << 1) + 2] = 0.;
8236 for (k = 1; k <= i_2; ++k) {
8237 ar[i + k * ar_dim1] /= scale;
8238 ai[i + k * ai_dim1] /= scale;
8239 h = h + ar[i + k * ar_dim1] * ar[i + k * ar_dim1] + ai[i + k *
8240 ai_dim1] * ai[i + k * ai_dim1];
8244 e2[i] = scale * scale * h;
8247 f =
pythag_(&ar[i + l * ar_dim1], &ai[i + l * ai_dim1]);
8252 tau[(l << 1) + 1] = (ai[i + l * ai_dim1] * tau[(i << 1) + 2] - ar[i +
8253 l * ar_dim1] * tau[(i << 1) + 1]) / f;
8254 si = (ar[i + l * ar_dim1] * tau[(i << 1) + 2] + ai[i + l * ai_dim1] *
8255 tau[(i << 1) + 1]) / f;
8258 ar[i + l * ar_dim1] = g * ar[i + l * ar_dim1];
8259 ai[i + l * ai_dim1] = g * ai[i + l * ai_dim1];
8265 tau[(l << 1) + 1] = -tau[(i << 1) + 1];
8266 si = tau[(i << 1) + 2];
8267 ar[i + l * ar_dim1] = g;
8272 for (j = 1; j <= i_2; ++j) {
8277 for (k = 1; k <= i_3; ++k) {
8278 g = g + ar[j + k * ar_dim1] * ar[i + k * ar_dim1] + ai[j + k *
8279 ai_dim1] * ai[i + k * ai_dim1];
8280 gi = gi - ar[j + k * ar_dim1] * ai[i + k * ai_dim1] + ai[j +
8281 k * ai_dim1] * ar[i + k * ar_dim1];
8291 for (k = jp1; k <= i_3; ++k) {
8292 g = g + ar[k + j * ar_dim1] * ar[i + k * ar_dim1] - ai[k + j *
8293 ai_dim1] * ai[i + k * ai_dim1];
8294 gi = gi - ar[k + j * ar_dim1] * ai[i + k * ai_dim1] - ai[k +
8295 j * ai_dim1] * ar[i + k * ar_dim1];
8301 tau[(j << 1) + 2] = gi / h;
8302 f = f + e[j] * ar[i + j * ar_dim1] - tau[(j << 1) + 2] * ai[i + j
8310 for (j = 1; j <= i_2; ++j) {
8311 f = ar[i + j * ar_dim1];
8314 fi = -ai[i + j * ai_dim1];
8315 gi = tau[(j << 1) + 2] - hh * fi;
8316 tau[(j << 1) + 2] = -gi;
8319 for (k = 1; k <= i_3; ++k) {
8320 ar[j + k * ar_dim1] = ar[j + k * ar_dim1] - f * e[k] - g * ar[
8321 i + k * ar_dim1] + fi * tau[(k << 1) + 2] + gi * ai[i
8323 ai[j + k * ai_dim1] = ai[j + k * ai_dim1] - f * tau[(k << 1)
8324 + 2] - g * ai[i + k * ai_dim1] - fi * e[k] - gi * ar[
8332 for (k = 1; k <= i_3; ++k) {
8333 ar[i + k * ar_dim1] = scale * ar[i + k * ar_dim1];
8334 ai[i + k * ai_dim1] = scale * ai[i + k * ai_dim1];
8338 tau[(l << 1) + 2] = -si;
8341 d[i] = ar[i + i * ar_dim1];
8342 ar[i + i * ar_dim1] = hh;
8343 ai[i + i * ai_dim1] = scale * sqrt(h);
8424 for (i = 2; i <= i_1; ++i) {
8432 for (l = 1; l <= i_1; ++l) {
8437 for (m = l; m <= i_2; ++m) {
8441 tst1 = (d_1 = d[m],
abs(d_1)) + (d_2 = d[m + 1],
abs(d_2));
8442 tst2 = tst1 + (d_1 = e[m],
abs(d_1));
8459 g = (d[l + 1] - p) / (e[l] * 2.);
8461 g = d[m] - p + e[l] / (g +
d_sign(&r, &g));
8468 for (ii = 1; ii <= i_2; ++ii) {
8480 r = (d[i] - g) * s + c * 2. * b;
8503 for (ii = 2; ii <= i_2; ++ii) {
8505 if (p >= d[i - 1]) {
8532 integer z_dim1, z_offset, i_1, i_2, i_3;
8609 z_offset = z_dim1 + 1;
8621 for (i = 2; i <= i_1; ++i) {
8629 for (l = 1; l <= i_1; ++l) {
8634 for (m = l; m <= i_2; ++m) {
8638 tst1 = (d_1 = d[m],
abs(d_1)) + (d_2 = d[m + 1],
abs(d_2));
8639 tst2 = tst1 + (d_1 = e[m],
abs(d_1));
8656 g = (d[l + 1] - p) / (e[l] * 2.);
8658 g = d[m] - p + e[l] / (g +
d_sign(&r, &g));
8665 for (ii = 1; ii <= i_2; ++ii) {
8677 r = (d[i] - g) * s + c * 2. * b;
8683 for (k = 1; k <= i_3; ++k) {
8684 f = z[k + (i + 1) * z_dim1];
8685 z[k + (i + 1) * z_dim1] = s * z[k + i * z_dim1] + c * f;
8686 z[k + i * z_dim1] = c * z[k + i * z_dim1] - s * f;
8707 for (ii = 2; ii <= i_1; ++ii) {
8713 for (j = ii; j <= i_2; ++j) {
8730 for (j = 1; j <= i_2; ++j) {
8731 p = z[j + i * z_dim1];
8732 z[j + i * z_dim1] = z[j + k * z_dim1];
8733 z[j + k * z_dim1] = p;
8846 for (i = 1; i <= i_1; ++i) {
8858 for (l = 1; l <= i_1; ++l) {
8863 for (m = l; m <= i_2; ++m) {
8867 tst1 = (d_1 = w[m],
abs(d_1)) + (d_2 = w[m + 1],
abs(d_2));
8868 tst2 = tst1 + (d_1 = rv1[m],
abs(d_1));
8874 if (e2[m + 1] == 0.) {
8900 g = (w[l + 1] - p) / (rv1[l] * 2.);
8902 g = w[m] - p + rv1[l] / (g +
d_sign(&r, &g));
8909 for (ii = 1; ii <= i_2; ++ii) {
8921 r = (w[i] - g) * s + c * 2. * b;
8944 for (ii = 2; ii <= i_2; ++ii) {
8946 if (p >= w[i - 1]) {
8950 ind[i] = ind[i - 1];
8977 integer a_dim1, a_offset, z_dim1, z_offset, rm1_dim1, rm1_offset, i_1,
8994 static integer ip, mp, ns, uk;
9093 rm1_offset = rm1_dim1 + 1;
9099 a_offset = a_dim1 + 1;
9102 z_offset = z_dim1 + 1;
9116 for (k = 1; k <= i_1; ++k) {
9117 if (wi[k] == 0. || ip < 0) {
9121 if (select[k] && select[k + 1]) {
9139 for (uk = k; uk <= i_2; ++uk) {
9143 if (a[uk + 1 + uk * a_dim1] == 0.) {
9155 for (i = 1; i <= i_2; ++i) {
9159 for (j = mp; j <= i_3; ++j) {
9161 x += (d_1 = a[i + j * a_dim1],
abs(d_1));
9178 ukroot = sqrt(ukroot);
9179 growto = .1 / ukroot;
9195 for (ii = 1; ii <= i_2; ++ii) {
9197 if (select[i] && (d_1 = wr[i] - rlambd,
abs(d_1)) < eps3 && (
9198 d_2 = wi[i] - ilambd,
abs(d_2)) < eps3) {
9214 for (i = 1; i <= i_2; ++i) {
9217 for (j = mp; j <= i_3; ++j) {
9219 rm1[j + i * rm1_dim1] = a[i + j * a_dim1];
9222 rm1[i + i * rm1_dim1] -= rlambd;
9240 for (i = 2; i <= i_2; ++i) {
9242 if ((d_1 = rm1[mp + i * rm1_dim1],
abs(d_1)) <= (d_2 = rm1[mp
9243 + mp * rm1_dim1],
abs(d_2))) {
9248 for (j = mp; j <= i_3; ++j) {
9249 y = rm1[j + i * rm1_dim1];
9250 rm1[j + i * rm1_dim1] = rm1[j + mp * rm1_dim1];
9251 rm1[j + mp * rm1_dim1] =
y;
9256 if (rm1[mp + mp * rm1_dim1] == 0.) {
9257 rm1[mp + mp * rm1_dim1] = eps3;
9259 x = rm1[mp + i * rm1_dim1] / rm1[mp + mp * rm1_dim1];
9265 for (j = i; j <= i_3; ++j) {
9267 rm1[j + i * rm1_dim1] -= x * rm1[j + mp * rm1_dim1];
9275 if (rm1[uk + uk * rm1_dim1] == 0.) {
9276 rm1[uk + uk * rm1_dim1] = eps3;
9282 for (ii = 1; ii <= i_2; ++ii) {
9291 for (j = ip1; j <= i_3; ++j) {
9293 y -= rm1[j + i * rm1_dim1] * rv1[j];
9297 rv1[i] = y / rm1[i + i * rm1_dim1];
9309 z[(s - 1) * z_dim1 + 1] = -ilambd;
9310 z[s * z_dim1 + 1] = 0.;
9314 rm1[rm1_dim1 * 3 + 1] = -ilambd;
9315 z[(s - 1) * z_dim1 + 1] = 0.;
9321 for (i = 4; i <= i_2; ++i) {
9323 rm1[i * rm1_dim1 + 1] = 0.;
9328 for (i = 2; i <= i_2; ++i) {
9330 w = rm1[mp + i * rm1_dim1];
9332 t = rm1[mp + (i + 1) * rm1_dim1];
9335 t = z[mp + (s - 1) * z_dim1];
9337 x = rm1[mp + mp * rm1_dim1] * rm1[mp + mp * rm1_dim1] + t * t;
9341 x = rm1[mp + mp * rm1_dim1] / w;
9343 rm1[mp + mp * rm1_dim1] = w;
9345 rm1[mp + (i + 1) * rm1_dim1] = 0.;
9348 z[mp + (s - 1) * z_dim1] = 0.;
9352 for (j = i; j <= i_3; ++j) {
9353 w = rm1[j + i * rm1_dim1];
9354 rm1[j + i * rm1_dim1] = rm1[j + mp * rm1_dim1] - x * w;
9355 rm1[j + mp * rm1_dim1] = w;
9360 z[i + l * z_dim1] = z[mp + l * z_dim1] - y * w;
9361 z[mp + l * z_dim1] = 0.;
9364 rm1[i + (j + 2) * rm1_dim1] = rm1[mp + (j + 2) * rm1_dim1] -
9366 rm1[mp + (j + 2) * rm1_dim1] = 0.;
9371 rm1[i + i * rm1_dim1] -= y * ilambd;
9376 z[mp + l * z_dim1] = -ilambd;
9377 z[i + l * z_dim1] += x * ilambd;
9380 rm1[mp + (i + 2) * rm1_dim1] = -ilambd;
9381 rm1[i + (i + 2) * rm1_dim1] += x * ilambd;
9387 rm1[mp + mp * rm1_dim1] = eps3;
9389 rm1[mp + (i + 1) * rm1_dim1] = 0.;
9392 z[mp + (s - 1) * z_dim1] = 0.;
9398 x = rm1[mp + mp * rm1_dim1] * w;
9402 for (j = i; j <= i_3; ++j) {
9407 t = z[mp + l * z_dim1];
9408 z[i + l * z_dim1] = -x * t - y * rm1[j + mp * rm1_dim1];
9411 t = rm1[mp + (j + 2) * rm1_dim1];
9412 rm1[i + (j + 2) * rm1_dim1] = -x * t - y * rm1[j + mp *
9415 rm1[j + i * rm1_dim1] = rm1[j + i * rm1_dim1] - x * rm1[j +
9416 mp * rm1_dim1] + y * t;
9424 z[i + l * z_dim1] -= ilambd;
9427 rm1[i + (i + 2) * rm1_dim1] -= ilambd;
9436 t = z[uk + l * z_dim1];
9439 t = rm1[uk + (uk + 2) * rm1_dim1];
9441 if (rm1[uk + uk * rm1_dim1] == 0. && t == 0.) {
9442 rm1[uk + uk * rm1_dim1] = eps3;
9448 for (ii = 1; ii <= i_2; ++ii) {
9458 for (j = ip1; j <= i_3; ++j) {
9463 t = z[i + l * z_dim1];
9466 t = rm1[i + (j + 2) * rm1_dim1];
9468 x = x - rm1[j + i * rm1_dim1] * rv1[j] + t * rv2[j];
9469 y = y - rm1[j + i * rm1_dim1] * rv2[j] - t * rv1[j];
9478 t = z[i + l * z_dim1];
9481 t = rm1[i + (i + 2) * rm1_dim1];
9483 cdiv_(&x, &y, &rm1[i + i * rm1_dim1], &t, &rv1[i], &rv2[i]);
9494 for (i = 1; i <= i_2; ++i) {
9496 x = (d_1 = rv1[i],
abs(d_1));
9499 x =
pythag_(&rv1[i], &rv2[i]);
9511 if (norm < growto) {
9524 for (i = 1; i <= i_2; ++i) {
9528 z[i + s * z_dim1] = rv1[i] * x;
9531 cdiv_(&rv1[i], &rv2[i], &x, &y, &z[i + (s - 1) * z_dim1], &z[i +
9549 y = eps3 / (x + 1.);
9553 for (i = 2; i <= i_2; ++i) {
9572 for (i = j; i <= i_2; ++i) {
9573 z[i + s * z_dim1] = 0.;
9575 z[i + (s - 1) * z_dim1] = 0.;
9600 *ierr = -((*n << 1) + 1);
9603 *m = s - 1 -
abs(ip);
9612 integer a_dim1, a_offset, b_dim1, b_offset, i_1, i_2, i_3;
9622 static integer i1, k1, l1, m1, ii, kk, ll;
9699 a_offset = a_dim1 + 1;
9702 b_offset = b_dim1 + 1;
9713 for (i = 1; i <= i_1; ++i) {
9724 for (k = i; k <= i_2; ++k) {
9726 scale += (d_1 = a[k + i * a_dim1],
abs(d_1));
9734 for (k = i; k <= i_2; ++k) {
9735 a[k + i * a_dim1] /= scale;
9737 d_1 = a[k + i * a_dim1];
9742 f = a[i + i * a_dim1];
9746 a[i + i * a_dim1] = f - g;
9752 for (j = l; j <= i_2; ++j) {
9756 for (k = i; k <= i_3; ++k) {
9758 s += a[k + i * a_dim1] * a[k + j * a_dim1];
9764 for (k = i; k <= i_3; ++k) {
9765 a[k + j * a_dim1] += f * a[k + i * a_dim1];
9776 for (j = 1; j <= i_3; ++j) {
9780 for (k = i; k <= i_2; ++k) {
9782 s += a[k + i * a_dim1] * b[k + j * b_dim1];
9788 for (k = i; k <= i_2; ++k) {
9789 b[k + j * b_dim1] += f * a[k + i * a_dim1];
9796 for (k = i; k <= i_2; ++k) {
9798 a[k + i * a_dim1] = scale * a[k + i * a_dim1];
9806 if (i > *m || i == *n) {
9811 for (k = l; k <= i_2; ++k) {
9813 scale += (d_1 = a[i + k * a_dim1],
abs(d_1));
9821 for (k = l; k <= i_2; ++k) {
9822 a[i + k * a_dim1] /= scale;
9824 d_1 = a[i + k * a_dim1];
9829 f = a[i + l * a_dim1];
9833 a[i + l * a_dim1] = f - g;
9836 for (k = l; k <= i_2; ++k) {
9838 rv1[k] = a[i + k * a_dim1] / h;
9846 for (j = l; j <= i_2; ++j) {
9850 for (k = l; k <= i_3; ++k) {
9852 s += a[j + k * a_dim1] * a[i + k * a_dim1];
9856 for (k = l; k <= i_3; ++k) {
9857 a[j + k * a_dim1] += s * rv1[k];
9864 for (k = l; k <= i_3; ++k) {
9866 a[i + k * a_dim1] = scale * a[i + k * a_dim1];
9871 d_3 = x, d_4 = (d_1 = w[i],
abs(d_1)) + (d_2 = rv1[i],
abs(d_2))
9879 for (ii = 1; ii <= i_1; ++ii) {
9889 for (j = l; j <= i_3; ++j) {
9893 a[j + i * a_dim1] = a[i + j * a_dim1] / a[i + l * a_dim1] / g;
9897 for (j = l; j <= i_3; ++j) {
9901 for (k = l; k <= i_2; ++k) {
9903 s += a[i + k * a_dim1] * a[k + j * a_dim1];
9907 for (k = l; k <= i_2; ++k) {
9908 a[k + j * a_dim1] += s * a[k + i * a_dim1];
9915 for (j = l; j <= i_2; ++j) {
9916 a[i + j * a_dim1] = 0.;
9917 a[j + i * a_dim1] = 0.;
9922 a[i + i * a_dim1] = 1.;
9928 if (*m >= *n || *ip == 0) {
9934 for (i = m1; i <= i_1; ++i) {
9937 for (j = 1; j <= i_2; ++j) {
9938 b[i + j * b_dim1] = 0.;
9947 for (kk = 1; kk <= i_2; ++kk) {
9955 for (ll = 1; ll <= i_1; ++ll) {
9958 tst2 = tst1 + (d_1 = rv1[l],
abs(d_1));
9964 tst2 = tst1 + (d_1 = w[l1],
abs(d_1));
9977 for (i = l; i <= i_1; ++i) {
9979 rv1[i] = c * rv1[i];
9980 tst2 = tst1 +
abs(f);
9994 for (j = 1; j <= i_3; ++j) {
9995 y = b[l1 + j * b_dim1];
9996 z = b[i + j * b_dim1];
9997 b[l1 + j * b_dim1] = y * c + z *
s;
9998 b[i + j * b_dim1] = -y * s + z * c;
10020 f = ((g + z) / h * ((g - z) /
y) + y / h - h / y) * .5;
10022 f = x - z / x * z + h / x * (y / (f +
d_sign(&g, &f)) - h);
10028 for (i1 = l; i1 <= i_1; ++i1) {
10039 g = -x * s + g * c;
10044 for (j = 1; j <= i_3; ++j) {
10045 x = a[j + i1 * a_dim1];
10046 z = a[j + i * a_dim1];
10047 a[j + i1 * a_dim1] = x * c + z *
s;
10048 a[j + i * a_dim1] = -x * s + z * c;
10063 x = -s * g + c *
y;
10069 for (j = 1; j <= i_3; ++j) {
10070 y = b[i1 + j * b_dim1];
10071 z = b[i + j * b_dim1];
10072 b[i1 + j * b_dim1] = y * c + z *
s;
10073 b[i + j * b_dim1] = -y * s + z * c;
10094 for (j = 1; j <= i_1; ++j) {
10096 a[j + k * a_dim1] = -a[j + k * a_dim1];
10116 integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2, i_3;
10120 static integer i, j, la, mm, mp, kp1, mp1;
10176 a_offset = a_dim1 + 1;
10179 z_offset = z_dim1 + 1;
10193 for (mm = kp1; mm <= i_1; ++mm) {
10194 mp = *low + *igh - mm;
10195 if (a[mp + (mp - 1) * a_dim1] == 0.) {
10201 for (i = mp1; i <= i_2; ++i) {
10203 ort[i] = a[i + (mp - 1) * a_dim1];
10207 for (j = 1; j <= i_2; ++j) {
10211 for (i = mp; i <= i_3; ++i) {
10213 g += ort[i] * z[i + j * z_dim1];
10219 g = g / ort[mp] / a[mp + (mp - 1) * a_dim1];
10222 for (i = mp; i <= i_3; ++i) {
10224 z[i + j * z_dim1] += g * ort[i];
10242 integer a_dim1, a_offset, i_1, i_2, i_3;
10252 static integer la, ii, jj, mp, kp1;
10300 a_offset = a_dim1 + 1;
10312 for (m = kp1; m <= i_1; ++m) {
10319 for (i = m; i <= i_2; ++i) {
10321 scale += (d_1 = a[i + (m - 1) * a_dim1],
abs(d_1));
10330 for (ii = m; ii <= i_2; ++ii) {
10332 ort[i] = a[i + (m - 1) * a_dim1] / scale;
10333 h += ort[i] * ort[i];
10338 g = -
d_sign(&d_1, &ort[m]);
10343 for (j = m; j <= i_2; ++j) {
10347 for (ii = m; ii <= i_3; ++ii) {
10349 f += ort[i] * a[i + j * a_dim1];
10356 for (i = m; i <= i_3; ++i) {
10358 a[i + j * a_dim1] -= f * ort[i];
10365 for (i = 1; i <= i_2; ++i) {
10369 for (jj = m; jj <= i_3; ++jj) {
10371 f += ort[j] * a[i + j * a_dim1];
10378 for (j = m; j <= i_3; ++j) {
10380 a[i + j * a_dim1] -= f * ort[j];
10386 ort[m] = scale * ort[m];
10387 a[m + (m - 1) * a_dim1] = scale * g;
10400 integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2, i_3;
10404 static integer i, j, kl, mm, mp, mp1;
10455 z_offset = z_dim1 + 1;
10459 a_offset = a_dim1 + 1;
10464 for (j = 1; j <= i_1; ++j) {
10467 for (i = 1; i <= i_2; ++i) {
10469 z[i + j * z_dim1] = 0.;
10472 z[j + j * z_dim1] = 1.;
10476 kl = *igh - *low - 1;
10482 for (mm = 1; mm <= i_1; ++mm) {
10484 if (a[mp + (mp - 1) * a_dim1] == 0.) {
10490 for (i = mp1; i <= i_2; ++i) {
10492 ort[i] = a[i + (mp - 1) * a_dim1];
10496 for (j = mp; j <= i_2; ++j) {
10500 for (i = mp; i <= i_3; ++i) {
10502 g += ort[i] * z[i + j * z_dim1];
10508 g = g / ort[mp] / a[mp + (mp - 1) * a_dim1];
10511 for (i = mp; i <= i_3; ++i) {
10513 z[i + j * z_dim1] += g * ort[i];
10531 integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i_1, i_2,
10543 static integer lb, nk1, nm1, nm2;
10598 z_offset = z_dim1 + 1;
10601 b_offset = b_dim1 + 1;
10604 a_offset = a_dim1 + 1;
10613 for (j = 1; j <= i_1; ++j) {
10616 for (i = 1; i <= i_2; ++i) {
10617 z[i + j * z_dim1] = 0.;
10621 z[j + j * z_dim1] = 1.;
10632 for (l = 1; l <= i_1; ++l) {
10637 for (i = l1; i <= i_2; ++i) {
10638 s += (d_1 = b[i + l * b_dim1],
abs(d_1));
10645 s += (d_1 = b[l + l * b_dim1],
abs(d_1));
10649 for (i = l; i <= i_2; ++i) {
10650 b[i + l * b_dim1] /=
s;
10652 d_1 = b[i + l * b_dim1];
10658 r =
d_sign(&d_1, &b[l + l * b_dim1]);
10659 b[l + l * b_dim1] += r;
10660 rho = r * b[l + l * b_dim1];
10663 for (j = l1; j <= i_2; ++j) {
10667 for (i = l; i <= i_3; ++i) {
10668 t += b[i + l * b_dim1] * b[i + j * b_dim1];
10675 for (i = l; i <= i_3; ++i) {
10676 b[i + j * b_dim1] += t * b[i + l * b_dim1];
10684 for (j = 1; j <= i_2; ++j) {
10688 for (i = l; i <= i_3; ++i) {
10689 t += b[i + l * b_dim1] * a[i + j * a_dim1];
10696 for (i = l; i <= i_3; ++i) {
10697 a[i + j * a_dim1] += t * b[i + l * b_dim1];
10704 b[l + l * b_dim1] = -s * r;
10707 for (i = l1; i <= i_2; ++i) {
10708 b[i + l * b_dim1] = 0.;
10723 for (k = 1; k <= i_1; ++k) {
10727 for (lb = 1; lb <= i_2; ++lb) {
10731 s = (d_1 = a[l + k * a_dim1],
abs(d_1)) + (d_2 = a[l1 + k *
10732 a_dim1],
abs(d_2));
10736 u1 = a[l + k * a_dim1] /
s;
10737 u2 = a[l1 + k * a_dim1] /
s;
10738 d_1 = sqrt(u1 * u1 + u2 * u2);
10740 v1 = -(u1 + r) / r;
10745 for (j = k; j <= i_3; ++j) {
10746 t = a[l + j * a_dim1] + u2 * a[l1 + j * a_dim1];
10747 a[l + j * a_dim1] += t * v1;
10748 a[l1 + j * a_dim1] += t * v2;
10752 a[l1 + k * a_dim1] = 0.;
10755 for (j = l; j <= i_3; ++j) {
10756 t = b[l + j * b_dim1] + u2 * b[l1 + j * b_dim1];
10757 b[l + j * b_dim1] += t * v1;
10758 b[l1 + j * b_dim1] += t * v2;
10762 s = (d_1 = b[l1 + l1 * b_dim1],
abs(d_1)) + (d_2 = b[l1 + l *
10763 b_dim1],
abs(d_2));
10767 u1 = b[l1 + l1 * b_dim1] /
s;
10768 u2 = b[l1 + l * b_dim1] /
s;
10769 d_1 = sqrt(u1 * u1 + u2 * u2);
10771 v1 = -(u1 + r) / r;
10776 for (i = 1; i <= i_3; ++i) {
10777 t = b[i + l1 * b_dim1] + u2 * b[i + l * b_dim1];
10778 b[i + l1 * b_dim1] += t * v1;
10779 b[i + l * b_dim1] += t * v2;
10783 b[l1 + l * b_dim1] = 0.;
10786 for (i = 1; i <= i_3; ++i) {
10787 t = a[i + l1 * a_dim1] + u2 * a[i + l * a_dim1];
10788 a[i + l1 * a_dim1] += t * v1;
10789 a[i + l * a_dim1] += t * v2;
10798 for (i = 1; i <= i_3; ++i) {
10799 t = z[i + l1 * z_dim1] + u2 * z[i + l * z_dim1];
10800 z[i + l1 * z_dim1] += t * v1;
10801 z[i + l * z_dim1] += t * v2;
10820 integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i_1, i_2,
10834 static doublereal u1, u2, u3, v1, v2, v3, a11, a12, a21, a22, a33, a34,
10835 a43, a44, b11, b12, b22, b33;
10846 static integer ish, itn, its, enm2, lor1;
10923 z_offset = z_dim1 + 1;
10926 b_offset = b_dim1 + 1;
10929 a_offset = a_dim1 + 1;
10939 for (i = 1; i <= i_1; ++i) {
10942 ani = (d_1 = a[i + (i - 1) * a_dim1],
abs(d_1));
10947 for (j = i; j <= i_2; ++j) {
10948 ani += (d_1 = a[i + j * a_dim1],
abs(d_1));
10949 bni += (d_1 = b[i + j * b_dim1],
abs(d_1));
10999 for (ll = 1; ll <= i_1; ++ll) {
11005 if ((d_1 = a[l + lm1 * a_dim1],
abs(d_1)) <= epsa) {
11012 a[l + lm1 * a_dim1] = 0.;
11024 b11 = b[l + l * b_dim1];
11025 if (
abs(b11) > epsb) {
11028 b[l + l * b_dim1] = 0.;
11029 s = (d_1 = a[l + l * a_dim1],
abs(d_1)) + (d_2 = a[l1 + l * a_dim1],
11031 u1 = a[l + l * a_dim1] /
s;
11032 u2 = a[l1 + l * a_dim1] /
s;
11033 d_1 = sqrt(u1 * u1 + u2 * u2);
11035 v1 = -(u1 + r) / r;
11040 for (j = l; j <= i_1; ++j) {
11041 t = a[l + j * a_dim1] + u2 * a[l1 + j * a_dim1];
11042 a[l + j * a_dim1] += t * v1;
11043 a[l1 + j * a_dim1] += t * v2;
11044 t = b[l + j * b_dim1] + u2 * b[l1 + j * b_dim1];
11045 b[l + j * b_dim1] += t * v1;
11046 b[l1 + j * b_dim1] += t * v2;
11051 a[l + lm1 * a_dim1] = -a[l + lm1 * a_dim1];
11057 a11 = a[l + l * a_dim1] / b11;
11058 a21 = a[l1 + l * a_dim1] / b11;
11070 b22 = b[l1 + l1 * b_dim1];
11071 if (
abs(b22) < epsb) {
11074 b33 = b[na + na * b_dim1];
11075 if (
abs(b33) < epsb) {
11078 b44 = b[en + en * b_dim1];
11079 if (
abs(b44) < epsb) {
11082 a33 = a[na + na * a_dim1] / b33;
11083 a34 = a[na + en * a_dim1] / b44;
11084 a43 = a[en + na * a_dim1] / b33;
11085 a44 = a[en + en * a_dim1] / b44;
11086 b34 = b[na + en * b_dim1] / b44;
11087 t = (a43 * b34 - a33 - a44) * .5;
11088 r = t * t + a34 * a43 - a33 * a44;
11097 if ((d_1 = s - a44,
abs(d_1)) < (d_2 = sh - a44,
abs(d_2))) {
11104 for (ll = ld; ll <= i_1; ++ll) {
11105 l = enm2 + ld - ll;
11111 t = a[l + l * a_dim1];
11112 if ((d_1 = b[l + l * b_dim1],
abs(d_1)) > epsb) {
11113 t -= sh * b[l + l * b_dim1];
11115 if ((d_1 = a[l + lm1 * a_dim1],
abs(d_1)) <= (d_2 = t / a[l1 + l *
11116 a_dim1],
abs(d_2)) * epsa) {
11126 a[l + lm1 * a_dim1] = -a[l + lm1 * a_dim1];
11131 a12 = a[l + l1 * a_dim1] / b22;
11132 a22 = a[l1 + l1 * a_dim1] / b22;
11133 b12 = b[l + l1 * b_dim1] / b22;
11134 a1 = ((a33 - a11) * (a44 - a11) - a34 * a43 + a43 * b34 * a11) / a21 +
11136 a2 = a22 - a11 - a21 * b12 - (a33 - a11) - (a44 - a11) + a43 * b34;
11137 a3 = a[l1 + 1 + l1 * a_dim1] / b22;
11152 for (k = l; k <= i_1; ++k) {
11153 notlas = k != na && ish == 2;
11160 i_2 = en, i_3 = k1 + ish;
11169 a1 = a[k + km1 * a_dim1];
11170 a2 = a[k1 + km1 * a_dim1];
11178 d_1 = sqrt(u1 * u1 + u2 * u2);
11180 v1 = -(u1 + r) / r;
11185 for (j = km1; j <= i_2; ++j) {
11186 t = a[k + j * a_dim1] + u2 * a[k1 + j * a_dim1];
11187 a[k + j * a_dim1] += t * v1;
11188 a[k1 + j * a_dim1] += t * v2;
11189 t = b[k + j * b_dim1] + u2 * b[k1 + j * b_dim1];
11190 b[k + j * b_dim1] += t * v1;
11191 b[k1 + j * b_dim1] += t * v2;
11196 a[k1 + km1 * a_dim1] = 0.;
11204 a1 = a[k + km1 * a_dim1];
11205 a2 = a[k1 + km1 * a_dim1];
11206 a3 = a[k2 + km1 * a_dim1];
11215 d_1 = sqrt(u1 * u1 + u2 * u2 + u3 * u3);
11217 v1 = -(u1 + r) / r;
11224 for (j = km1; j <= i_2; ++j) {
11225 t = a[k + j * a_dim1] + u2 * a[k1 + j * a_dim1] + u3 * a[k2 + j *
11227 a[k + j * a_dim1] += t * v1;
11228 a[k1 + j * a_dim1] += t * v2;
11229 a[k2 + j * a_dim1] += t * v3;
11230 t = b[k + j * b_dim1] + u2 * b[k1 + j * b_dim1] + u3 * b[k2 + j *
11232 b[k + j * b_dim1] += t * v1;
11233 b[k1 + j * b_dim1] += t * v2;
11234 b[k2 + j * b_dim1] += t * v3;
11241 a[k1 + km1 * a_dim1] = 0.;
11242 a[k2 + km1 * a_dim1] = 0.;
11245 s = (d_1 = b[k2 + k2 * b_dim1],
abs(d_1)) + (d_2 = b[k2 + k1 *
11246 b_dim1],
abs(d_2)) + (d_3 = b[k2 + k * b_dim1],
abs(d_3));
11250 u1 = b[k2 + k2 * b_dim1] /
s;
11251 u2 = b[k2 + k1 * b_dim1] /
s;
11252 u3 = b[k2 + k * b_dim1] /
s;
11253 d_1 = sqrt(u1 * u1 + u2 * u2 + u3 * u3);
11255 v1 = -(u1 + r) / r;
11262 for (i = lor1; i <= i_2; ++i) {
11263 t = a[i + k2 * a_dim1] + u2 * a[i + k1 * a_dim1] + u3 * a[i + k *
11265 a[i + k2 * a_dim1] += t * v1;
11266 a[i + k1 * a_dim1] += t * v2;
11267 a[i + k * a_dim1] += t * v3;
11268 t = b[i + k2 * b_dim1] + u2 * b[i + k1 * b_dim1] + u3 * b[i + k *
11270 b[i + k2 * b_dim1] += t * v1;
11271 b[i + k1 * b_dim1] += t * v2;
11272 b[i + k * b_dim1] += t * v3;
11276 b[k2 + k * b_dim1] = 0.;
11277 b[k2 + k1 * b_dim1] = 0.;
11283 for (i = 1; i <= i_2; ++i) {
11284 t = z[i + k2 * z_dim1] + u2 * z[i + k1 * z_dim1] + u3 * z[i + k *
11286 z[i + k2 * z_dim1] += t * v1;
11287 z[i + k1 * z_dim1] += t * v2;
11288 z[i + k * z_dim1] += t * v3;
11293 s = (d_1 = b[k1 + k1 * b_dim1],
abs(d_1)) + (d_2 = b[k1 + k *
11294 b_dim1],
abs(d_2));
11298 u1 = b[k1 + k1 * b_dim1] /
s;
11299 u2 = b[k1 + k * b_dim1] /
s;
11300 d_1 = sqrt(u1 * u1 + u2 * u2);
11302 v1 = -(u1 + r) / r;
11307 for (i = lor1; i <= i_2; ++i) {
11308 t = a[i + k1 * a_dim1] + u2 * a[i + k * a_dim1];
11309 a[i + k1 * a_dim1] += t * v1;
11310 a[i + k * a_dim1] += t * v2;
11311 t = b[i + k1 * b_dim1] + u2 * b[i + k * b_dim1];
11312 b[i + k1 * b_dim1] += t * v1;
11313 b[i + k * b_dim1] += t * v2;
11317 b[k1 + k * b_dim1] = 0.;
11323 for (i = 1; i <= i_2; ++i) {
11324 t = z[i + k1 * z_dim1] + u2 * z[i + k * z_dim1];
11325 z[i + k1 * z_dim1] += t * v1;
11326 z[i + k * z_dim1] += t * v2;
11342 b[*n + b_dim1] = epsb;
11352 integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i_1, i_2;
11361 static doublereal r,
s, t, a1, a2, u1, u2, v1, v2, a11, a12, a21, a22,
11362 b11, b12, b22, di, ei;
11368 static doublereal cz, ti, tr, a1i, a2i, a11i, a12i, a22i, a11r, a12r,
11445 z_offset = z_dim1 + 1;
11451 b_offset = b_dim1 + 1;
11454 a_offset = a_dim1 + 1;
11458 epsb = b[*n + b_dim1];
11463 for (nn = 1; nn <= i_1; ++nn) {
11472 if (a[en + na * a_dim1] != 0.) {
11477 alfr[en] = a[en + en * a_dim1];
11478 if (b[en + en * b_dim1] < 0.) {
11479 alfr[en] = -alfr[en];
11481 beta[en] = (d_1 = b[en + en * b_dim1],
abs(d_1));
11486 if ((d_1 = b[na + na * b_dim1],
abs(d_1)) <= epsb) {
11489 if ((d_1 = b[en + en * b_dim1],
abs(d_1)) > epsb) {
11492 a1 = a[en + en * a_dim1];
11493 a2 = a[en + na * a_dim1];
11497 an = (d_1 = a[na + na * a_dim1],
abs(d_1)) + (d_2 = a[na + en *
11498 a_dim1],
abs(d_2)) + (d_3 = a[en + na * a_dim1],
abs(d_3))
11499 + (d_4 = a[en + en * a_dim1],
abs(d_4));
11500 bn = (d_1 = b[na + na * b_dim1],
abs(d_1)) + (d_2 = b[na + en *
11501 b_dim1],
abs(d_2)) + (d_3 = b[en + en * b_dim1],
abs(d_3));
11502 a11 = a[na + na * a_dim1] / an;
11503 a12 = a[na + en * a_dim1] / an;
11504 a21 = a[en + na * a_dim1] / an;
11505 a22 = a[en + en * a_dim1] / an;
11506 b11 = b[na + na * b_dim1] / bn;
11507 b12 = b[na + en * b_dim1] / bn;
11508 b22 = b[en + en * b_dim1] / bn;
11511 s = a21 / (b11 * b22);
11512 t = (a22 - e * b22) / b22;
11513 if (
abs(e) <=
abs(ei)) {
11517 t = (a11 - e * b11) / b11;
11519 c = (t - s * b12) * .5;
11520 d = c * c + s * (a12 - e * b12);
11527 e += c +
d_sign(&d_1, &c);
11545 d_1 = sqrt(u1 * u1 + u2 * u2);
11547 v1 = -(u1 + r) / r;
11552 for (i = 1; i <= i_2; ++i) {
11553 t = a[i + en * a_dim1] + u2 * a[i + na * a_dim1];
11554 a[i + en * a_dim1] += t * v1;
11555 a[i + na * a_dim1] += t * v2;
11556 t = b[i + en * b_dim1] + u2 * b[i + na * b_dim1];
11557 b[i + en * b_dim1] += t * v1;
11558 b[i + na * b_dim1] += t * v2;
11567 for (i = 1; i <= i_2; ++i) {
11568 t = z[i + en * z_dim1] + u2 * z[i + na * z_dim1];
11569 z[i + en * z_dim1] += t * v1;
11570 z[i + na * z_dim1] += t * v2;
11578 if (an <
abs(e) * bn) {
11581 a1 = b[na + na * b_dim1];
11582 a2 = b[en + na * b_dim1];
11585 a1 = a[na + na * a_dim1];
11586 a2 = a[en + na * a_dim1];
11595 d_1 = sqrt(u1 * u1 + u2 * u2);
11597 v1 = -(u1 + r) / r;
11602 for (j = na; j <= i_2; ++j) {
11603 t = a[na + j * a_dim1] + u2 * a[en + j * a_dim1];
11604 a[na + j * a_dim1] += t * v1;
11605 a[en + j * a_dim1] += t * v2;
11606 t = b[na + j * b_dim1] + u2 * b[en + j * b_dim1];
11607 b[na + j * b_dim1] += t * v1;
11608 b[en + j * b_dim1] += t * v2;
11613 a[en + na * a_dim1] = 0.;
11614 b[en + na * b_dim1] = 0.;
11615 alfr[na] = a[na + na * a_dim1];
11616 alfr[en] = a[en + en * a_dim1];
11617 if (b[na + na * b_dim1] < 0.) {
11618 alfr[na] = -alfr[na];
11620 if (b[en + en * b_dim1] < 0.) {
11621 alfr[en] = -alfr[en];
11623 beta[na] = (d_1 = b[na + na * b_dim1],
abs(d_1));
11624 beta[en] = (d_1 = b[en + en * b_dim1],
abs(d_1));
11632 a11r = a11 - e * b11;
11634 a12r = a12 - e * b12;
11636 a22r = a22 - e * b22;
11639 a22r) +
abs(a22i)) {
11654 cz = sqrt(a1 * a1 + a1i * a1i);
11658 szr = (a1 * a2 + a1i * a2i) / cz;
11659 szi = (a1 * a2i - a1i * a2) / cz;
11660 r = sqrt(cz * cz + szr * szr + szi * szi);
11669 if (an < (
abs(e) + ei) * bn) {
11672 a1 = cz * b11 + szr * b12;
11678 a1 = cz * a11 + szr * a12;
11680 a2 = cz * a21 + szr * a22;
11684 cq = sqrt(a1 * a1 + a1i * a1i);
11688 sqr = (a1 * a2 + a1i * a2i) / cq;
11689 sqi = (a1 * a2i - a1i * a2) / cq;
11690 r = sqrt(cq * cq + sqr * sqr + sqi * sqi);
11701 ssr = sqr * szr + sqi * szi;
11702 ssi = sqr * szi - sqi * szr;
11704 tr = cq * cz * a11 + cq * szr * a12 + sqr * cz * a21 + ssr * a22;
11705 ti = cq * szi * a12 - sqi * cz * a21 + ssi * a22;
11706 dr = cq * cz * b11 + cq * szr * b12 + ssr * b22;
11707 di = cq * szi * b12 + ssi * b22;
11711 tr = ssr * a11 - sqr * cz * a12 - cq * szr * a21 + cq * cz * a22;
11712 ti = -ssi * a11 - sqi * cz * a12 + cq * szi * a21;
11713 dr = ssr * b11 - sqr * cz * b12 + cq * cz * b22;
11714 di = -ssi * b11 - sqi * cz * b12;
11716 t = ti * dr - tr * di;
11721 r = sqrt(dr * dr + di * di);
11723 alfr[j] = an * (tr * dr + ti * di) / r;
11724 alfi[j] = an * t / r;
11733 b[*n + b_dim1] = epsb;
11743 integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i_1, i_2,
11751 static doublereal alfm, almi, betm, epsb, almr, d;
11753 static doublereal q, r,
s, t, w, x,
y, t1, t2, w1, x1, z1, di;
11754 static integer na, ii, en, jj;
11830 z_offset = z_dim1 + 1;
11836 b_offset = b_dim1 + 1;
11839 a_offset = a_dim1 + 1;
11843 epsb = b[*n + b_dim1];
11847 for (nn = 1; nn <= i_1; ++nn) {
11853 if (alfi[en] != 0.) {
11858 b[en + en * b_dim1] = 1.;
11866 for (ii = 1; ii <= i_2; ++ii) {
11868 w = betm * a[i + i * a_dim1] - alfm * b[i + i * b_dim1];
11872 for (j = m; j <= i_3; ++j) {
11874 r += (betm * a[i + j * a_dim1] - alfm * b[i + j * b_dim1]) *
11875 b[j + en * b_dim1];
11878 if (i == 1 || isw == 2) {
11881 if (betm * a[i + (i - 1) * a_dim1] == 0.) {
11897 b[i + en * b_dim1] = -r / t;
11901 x = betm * a[i + (i + 1) * a_dim1] - alfm * b[i + (i + 1) *
11903 y = betm * a[i + 1 + i * a_dim1];
11904 q = w * zz - x *
y;
11905 t = (x * s - zz * r) / q;
11906 b[i + en * b_dim1] = t;
11907 if (
abs(x) <=
abs(zz)) {
11910 b[i + 1 + en * b_dim1] = (-r - w * t) / x;
11913 b[i + 1 + en * b_dim1] = (-s - y * t) / zz;
11929 y = betm * a[en + na * a_dim1];
11930 b[na + na * b_dim1] = -almi * b[en + en * b_dim1] /
y;
11931 b[na + en * b_dim1] = (almr * b[en + en * b_dim1] - betm * a[en + en *
11933 b[en + na * b_dim1] = 0.;
11934 b[en + en * b_dim1] = 1.;
11941 for (ii = 1; ii <= i_2; ++ii) {
11943 w = betm * a[i + i * a_dim1] - almr * b[i + i * b_dim1];
11944 w1 = -almi * b[i + i * b_dim1];
11949 for (j = m; j <= i_3; ++j) {
11950 x = betm * a[i + j * a_dim1] - almr * b[i + j * b_dim1];
11951 x1 = -almi * b[i + j * b_dim1];
11952 ra = ra + x * b[j + na * b_dim1] - x1 * b[j + en * b_dim1];
11953 sa = sa + x * b[j + en * b_dim1] + x1 * b[j + na * b_dim1];
11957 if (i == 1 || isw == 2) {
11960 if (betm * a[i + (i - 1) * a_dim1] == 0.) {
11983 if (
abs(di) >
abs(dr)) {
11988 t1 = (tr + ti * rr) / d;
11989 t2 = (ti - tr * rr) / d;
11997 t1 = (tr * rr + ti) / d;
11998 t2 = (ti * rr - tr) / d;
12005 x = betm * a[i + (i + 1) * a_dim1] - almr * b[i + (i + 1) *
12007 x1 = -almi * b[i + (i + 1) * b_dim1];
12008 y = betm * a[i + 1 + i * a_dim1];
12009 tr = y * ra - w * r + w1 *
s;
12010 ti = y * sa - w * s - w1 * r;
12011 dr = w * zz - w1 * z1 - x *
y;
12012 di = w * z1 + w1 * zz - x1 *
y;
12013 if (dr == 0. && di == 0.) {
12018 b[i + 1 + na * b_dim1] = t1;
12019 b[i + 1 + en * b_dim1] = t2;
12024 tr = -ra - x * b[i + 1 + na * b_dim1] + x1 * b[i + 1 + en *
12026 ti = -sa - x * b[i + 1 + en * b_dim1] - x1 * b[i + 1 + na *
12030 t1 = (-r - zz * b[i + 1 + na * b_dim1] + z1 * b[i + 1 + en *
12032 t2 = (-s - zz * b[i + 1 + en * b_dim1] - z1 * b[i + 1 + na *
12035 b[i + na * b_dim1] = t1;
12036 b[i + en * b_dim1] = t2;
12050 for (jj = 1; jj <= i_1; ++jj) {
12054 for (i = 1; i <= i_2; ++i) {
12058 for (k = 1; k <= i_3; ++k) {
12060 zz += z[i + k * z_dim1] * b[k + j * b_dim1];
12063 z[i + j * z_dim1] = zz;
12071 for (j = 1; j <= i_2; ++j) {
12076 if (alfi[j] != 0.) {
12081 for (i = 1; i <= i_1; ++i) {
12082 if ((d_1 = z[i + j * z_dim1],
abs(d_1)) > d) {
12083 d = (d_2 = z[i + j * z_dim1],
abs(d_2));
12089 for (i = 1; i <= i_1; ++i) {
12091 z[i + j * z_dim1] /= d;
12098 for (i = 1; i <= i_1; ++i) {
12099 r = (d_1 = z[i + (j - 1) * z_dim1],
abs(d_1)) + (d_2 = z[i + j
12100 * z_dim1],
abs(d_2));
12103 d_1 = z[i + (j - 1) * z_dim1] / r;
12105 d_2 = z[i + j * z_dim1] / r;
12106 r *= sqrt(d_1 * d_1 + d_2 * d_2);
12115 for (i = 1; i <= i_1; ++i) {
12116 z[i + (j - 1) * z_dim1] /= d;
12117 z[i + j * z_dim1] /= d;
12258 for (i = 1; i <= i_1; ++i) {
12279 for (i = 1; i <= i_1; ++i) {
12284 d_3 = (d_1 = d[i],
abs(d_1)) + (d_2 = d[i - 1],
abs(d_2));
12300 q = (d_1 = e[i + 1],
abs(d_1));
12303 d_1 = w[i] - p - q;
12304 tot =
min(d_1,tot);
12308 if (jdef == 1 && tot < 0.) {
12313 for (i = 1; i <= i_1; ++i) {
12324 for (k = 1; k <= i_1; ++k) {
12334 if (delta > *eps1) {
12337 if (delta < -(*eps1)) {
12349 for (j = k1; j <= i_2; ++j) {
12350 d_2 = w[j] + w[j - 1];
12353 if (bd[j] <= d_1 * d_1) {
12360 f = bd[*n] / delta;
12369 for (ii = 1; ii <= i_2; ++ii) {
12375 w[i + 1] = qp + ep;
12377 if (delta > *eps1) {
12380 if (delta < -(*eps1)) {
12387 bd[i + 1] = qp * ep;
12394 if (tot + s > tot) {
12399 *ierr = *n * 5 + k;
12404 for (j = k; j <= i_2; ++j) {
12405 if (w[j] > delta) {
12416 bd[i + 1] = bd[i] * f / qp;
12425 for (jj = 1; jj <= i_2; ++jj) {
12427 w[j + 1] = w[j] -
s;
12429 ind[j + 1] = ind[j];
12451 for (i = 1; i <= i_1; ++i) {
12459 case 2:
goto L1001;
12463 *ierr = *n * 6 + 1;
12472 integer b_dim1, b_offset, z_dim1, z_offset, i_1, i_2, i_3;
12525 b_offset = b_dim1 + 1;
12528 z_offset = z_dim1 + 1;
12537 for (j = 1; j <= i_1; ++j) {
12540 for (ii = 1; ii <= i_2; ++ii) {
12543 x = z[i + j * z_dim1];
12549 for (k = i1; k <= i_3; ++k) {
12551 x -= b[k + i * b_dim1] * z[k + j * z_dim1];
12555 z[i + j * z_dim1] = x / dl[i];
12568 integer b_dim1, b_offset, z_dim1, z_offset, i_1, i_2, i_3;
12621 b_offset = b_dim1 + 1;
12624 z_offset = z_dim1 + 1;
12633 for (j = 1; j <= i_1; ++j) {
12636 for (ii = 1; ii <= i_2; ++ii) {
12639 x = dl[i] * z[i + j * z_dim1];
12645 for (k = 1; k <= i_3; ++k) {
12647 x += b[i + k * b_dim1] * z[k + j * z_dim1];
12651 z[i + j * z_dim1] = x;
12664 integer a_dim1, a_offset, b_dim1, b_offset, i_1, i_2, i_3;
12731 b_offset = b_dim1 + 1;
12734 a_offset = a_dim1 + 1;
12745 for (i = 1; i <= i_1; ++i) {
12749 for (j = i; j <= i_2; ++j) {
12750 x = b[i + j * b_dim1];
12756 for (k = 1; k <= i_3; ++k) {
12758 x -= b[i + k * b_dim1] * b[j + k * b_dim1];
12772 b[j + i * b_dim1] = x /
y;
12781 for (i = 1; i <= i_2; ++i) {
12786 for (j = i; j <= i_1; ++j) {
12787 x = a[i + j * a_dim1];
12793 for (k = 1; k <= i_3; ++k) {
12795 x -= b[i + k * b_dim1] * a[j + k * a_dim1];
12799 a[j + i * a_dim1] = x /
y;
12805 for (j = 1; j <= i_1; ++j) {
12809 for (i = j; i <= i_2; ++i) {
12810 x = a[i + j * a_dim1];
12817 for (k = j; k <= i_3; ++k) {
12819 x -= a[k + j * a_dim1] * b[i + k * b_dim1];
12828 for (k = 1; k <= i_3; ++k) {
12830 x -= a[j + k * a_dim1] * b[i + k * b_dim1];
12834 a[i + j * a_dim1] = x / dl[i];
12842 *ierr = *n * 7 + 1;
12851 integer a_dim1, a_offset, b_dim1, b_offset, i_1, i_2, i_3;
12919 b_offset = b_dim1 + 1;
12922 a_offset = a_dim1 + 1;
12933 for (i = 1; i <= i_1; ++i) {
12937 for (j = i; j <= i_2; ++j) {
12938 x = b[i + j * b_dim1];
12944 for (k = 1; k <= i_3; ++k) {
12946 x -= b[i + k * b_dim1] * b[j + k * b_dim1];
12960 b[j + i * b_dim1] = x /
y;
12969 for (i = 1; i <= i_2; ++i) {
12973 for (j = 1; j <= i_1; ++j) {
12974 x = a[j + i * a_dim1] * dl[j];
12981 for (k = j1; k <= i_3; ++k) {
12983 x += a[k + i * a_dim1] * b[k + j * b_dim1];
12992 for (k = i1; k <= i_3; ++k) {
12994 x += a[i + k * a_dim1] * b[k + j * b_dim1];
12998 a[i + j * a_dim1] = x;
13004 for (i = 1; i <= i_1; ++i) {
13009 for (j = 1; j <= i_2; ++j) {
13010 x = y * a[i + j * a_dim1];
13016 for (k = i1; k <= i_3; ++k) {
13018 x += a[k + j * a_dim1] * b[k + i * b_dim1];
13022 a[i + j * a_dim1] = x;
13030 *ierr = *n * 7 + 1;
13040 integer a_dim1, a_offset, z_dim1, z_offset;
13110 z_offset = z_dim1 + 1;
13115 a_offset = a_dim1 + 1;
13126 balanc_(nm, n, &a[a_offset], &is1, &is2, &fv1[1]);
13127 elmhes_(nm, n, &is1, &is2, &a[a_offset], &iv1[1]);
13132 hqr_(nm, n, &is1, &is2, &a[a_offset], &wr[1], &wi[1], ierr);
13136 eltran_(nm, n, &is1, &is2, &a[a_offset], &iv1[1], &z[z_offset]);
13137 hqr2_(nm, n, &is1, &is2, &a[a_offset], &wr[1], &wi[1], &z[z_offset], ierr)
13142 balbak_(nm, n, &is1, &is2, &fv1[1], n, &z[z_offset]);
13152 integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset;
13221 z_offset = z_dim1 + 1;
13227 b_offset = b_dim1 + 1;
13230 a_offset = a_dim1 + 1;
13246 qzhes_(nm, n, &a[a_offset], &b[b_offset], &tf, &z[z_offset]);
13247 qzit_(nm, n, &a[a_offset], &b[b_offset], &c_b550, &tf, &z[z_offset], ierr)
13249 qzval_(nm, n, &a[a_offset], &b[b_offset], &alfr[1], &alfi[1], &beta[1], &
13255 qzhes_(nm, n, &a[a_offset], &b[b_offset], &tf, &z[z_offset]);
13256 qzit_(nm, n, &a[a_offset], &b[b_offset], &c_b550, &tf, &z[z_offset], ierr)
13258 qzval_(nm, n, &a[a_offset], &b[b_offset], &alfr[1], &alfi[1], &beta[1], &
13263 qzvec_(nm, n, &a[a_offset], &b[b_offset], &alfr[1], &alfi[1], &beta[1], &
13274 integer a_dim1, a_offset, z_dim1, z_offset;
13330 z_offset = z_dim1 + 1;
13334 a_offset = a_dim1 + 1;
13349 tred1_(nm, n, &a[a_offset], &w[1], &fv1[1], &fv2[1]);
13350 tqlrat_(n, &w[1], &fv2[1], ierr);
13354 tred2_(nm, n, &a[a_offset], &w[1], &fv1[1], &z[z_offset]);
13355 tql2_(nm, n, &w[1], &fv1[1], &z[z_offset], ierr);
13365 integer a_dim1, a_offset, z_dim1, z_offset;
13434 z_offset = z_dim1 + 1;
13438 a_offset = a_dim1 + 1;
13466 bandr_(nm, n, mb, &a[a_offset], &w[1], &fv1[1], &fv2[1], &tf, &z[z_offset]
13468 tqlrat_(n, &w[1], &fv2[1], ierr);
13473 bandr_(nm, n, mb, &a[a_offset], &w[1], &fv1[1], &fv1[1], &tf, &z[z_offset]
13475 tql2_(nm, n, &w[1], &fv1[1], &z[z_offset], ierr);
13485 integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset;
13546 z_offset = z_dim1 + 1;
13550 b_offset = b_dim1 + 1;
13553 a_offset = a_dim1 + 1;
13564 reduc_(nm, n, &a[a_offset], &b[b_offset], &fv2[1], ierr);
13572 tred1_(nm, n, &a[a_offset], &w[1], &fv1[1], &fv2[1]);
13573 tqlrat_(n, &w[1], &fv2[1], ierr);
13577 tred2_(nm, n, &a[a_offset], &w[1], &fv1[1], &z[z_offset]);
13578 tql2_(nm, n, &w[1], &fv1[1], &z[z_offset], ierr);
13582 rebak_(nm, n, &b[b_offset], &fv2[1], n, &z[z_offset]);
13592 integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset;
13653 z_offset = z_dim1 + 1;
13657 b_offset = b_dim1 + 1;
13660 a_offset = a_dim1 + 1;
13671 reduc2_(nm, n, &a[a_offset], &b[b_offset], &fv2[1], ierr);
13679 tred1_(nm, n, &a[a_offset], &w[1], &fv1[1], &fv2[1]);
13680 tqlrat_(n, &w[1], &fv2[1], ierr);
13684 tred2_(nm, n, &a[a_offset], &w[1], &fv1[1], &z[z_offset]);
13685 tql2_(nm, n, &w[1], &fv1[1], &z[z_offset], ierr);
13689 rebak_(nm, n, &b[b_offset], &fv2[1], n, &z[z_offset]);
13699 integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset;
13760 z_offset = z_dim1 + 1;
13764 b_offset = b_dim1 + 1;
13767 a_offset = a_dim1 + 1;
13778 reduc2_(nm, n, &a[a_offset], &b[b_offset], &fv2[1], ierr);
13786 tred1_(nm, n, &a[a_offset], &w[1], &fv1[1], &fv2[1]);
13787 tqlrat_(n, &w[1], &fv2[1], ierr);
13791 tred2_(nm, n, &a[a_offset], &w[1], &fv1[1], &z[z_offset]);
13792 tql2_(nm, n, &w[1], &fv1[1], &z[z_offset], ierr);
13796 rebakb_(nm, n, &b[b_offset], &fv2[1], n, &z[z_offset]);
13806 integer a_dim1, a_offset, z_dim1, z_offset;
13811 static integer k1, k2, k3, k4, k5, k6, k7, k8;
13871 a_offset = a_dim1 + 1;
13874 z_offset = z_dim1 + 1;
13880 if (*n > *nm || *m > *nm) {
13895 tred1_(nm, n, &a[a_offset], &w[1], &fwork[k1], &fwork[k2]);
13896 tqlrat_(n, &w[1], &fwork[k2], ierr);
13900 tred1_(nm, n, &a[a_offset], &fwork[k1], &fwork[k2], &fwork[k3]);
13901 imtqlv_(n, &fwork[k1], &fwork[k2], &fwork[k3], &w[1], &iwork[1], ierr, &
13903 tinvit_(nm, n, &fwork[k1], &fwork[k2], &fwork[k3], m, &w[1], &iwork[1], &
13904 z[z_offset], ierr, &fwork[k4], &fwork[k5], &fwork[k6], &fwork[k7],
13906 trbak1_(nm, n, &a[a_offset], &fwork[k2], m, &z[z_offset]);
13916 integer z_dim1, z_offset, i_1, i_2;
13979 z_offset = z_dim1 + 1;
13991 if (*nv >= *n * (*n + 1) / 2) {
13998 tred3_(n, nv, &a[1], &w[1], &fv1[1], &fv2[1]);
14003 tqlrat_(n, &w[1], &fv2[1], ierr);
14008 for (i = 1; i <= i_1; ++i) {
14011 for (j = 1; j <= i_2; ++j) {
14012 z[j + i * z_dim1] = 0.;
14016 z[i + i * z_dim1] = 1.;
14020 tql2_(nm, n, &w[1], &fv1[1], &z[z_offset], ierr);
14024 trbak3_(nm, n, nv, &a[1], n, &z[z_offset]);
14033 integer z_dim1, z_offset, i_1, i_2;
14087 z_offset = z_dim1 + 1;
14104 imtql1_(n, &w[1], &e[1], ierr);
14109 for (i = 1; i <= i_1; ++i) {
14112 for (j = 1; j <= i_2; ++j) {
14113 z[j + i * z_dim1] = 0.;
14117 z[i + i * z_dim1] = 1.;
14121 imtql2_(nm, n, &w[1], &e[1], &z[z_offset], ierr);
14130 integer a_dim1, a_offset, z_dim1, z_offset;
14189 a_offset = a_dim1 + 1;
14193 z_offset = z_dim1 + 1;
14209 figi_(nm, n, &a[a_offset], &w[1], &fv1[1], &fv1[1], ierr);
14213 imtql1_(n, &w[1], &fv1[1], ierr);
14217 figi2_(nm, n, &a[a_offset], &w[1], &fv1[1], &z[z_offset], ierr);
14221 imtql2_(nm, n, &w[1], &fv1[1], &z[z_offset], ierr);
14231 integer a_dim1, a_offset, u_dim1, u_offset, v_dim1, v_offset, i_1, i_2,
14242 static integer i1, k1, l1, ii, kk, ll, mn;
14319 v_offset = v_dim1 + 1;
14322 u_offset = u_dim1 + 1;
14326 a_offset = a_dim1 + 1;
14333 for (i = 1; i <= i_1; ++i) {
14336 for (j = 1; j <= i_2; ++j) {
14337 u[i + j * u_dim1] = a[i + j * a_dim1];
14347 for (i = 1; i <= i_2; ++i) {
14349 rv1[i] = scale * g;
14358 for (k = i; k <= i_1; ++k) {
14360 scale += (d_1 = u[k + i * u_dim1],
abs(d_1));
14368 for (k = i; k <= i_1; ++k) {
14369 u[k + i * u_dim1] /= scale;
14371 d_1 = u[k + i * u_dim1];
14376 f = u[i + i * u_dim1];
14380 u[i + i * u_dim1] = f - g;
14386 for (j = l; j <= i_1; ++j) {
14390 for (k = i; k <= i_3; ++k) {
14392 s += u[k + i * u_dim1] * u[k + j * u_dim1];
14398 for (k = i; k <= i_3; ++k) {
14399 u[k + j * u_dim1] += f * u[k + i * u_dim1];
14406 for (k = i; k <= i_3; ++k) {
14408 u[k + i * u_dim1] = scale * u[k + i * u_dim1];
14416 if (i > *m || i == *n) {
14421 for (k = l; k <= i_3; ++k) {
14423 scale += (d_1 = u[i + k * u_dim1],
abs(d_1));
14431 for (k = l; k <= i_3; ++k) {
14432 u[i + k * u_dim1] /= scale;
14434 d_1 = u[i + k * u_dim1];
14439 f = u[i + l * u_dim1];
14443 u[i + l * u_dim1] = f - g;
14446 for (k = l; k <= i_3; ++k) {
14448 rv1[k] = u[i + k * u_dim1] / h;
14456 for (j = l; j <= i_3; ++j) {
14460 for (k = l; k <= i_1; ++k) {
14462 s += u[j + k * u_dim1] * u[i + k * u_dim1];
14466 for (k = l; k <= i_1; ++k) {
14467 u[j + k * u_dim1] += s * rv1[k];
14474 for (k = l; k <= i_1; ++k) {
14476 u[i + k * u_dim1] = scale * u[i + k * u_dim1];
14481 d_3 = x, d_4 = (d_1 = w[i],
abs(d_1)) + (d_2 = rv1[i],
abs(d_2))
14492 for (ii = 1; ii <= i_2; ++ii) {
14502 for (j = l; j <= i_1; ++j) {
14506 v[j + i * v_dim1] = u[i + j * u_dim1] / u[i + l * u_dim1] / g;
14510 for (j = l; j <= i_1; ++j) {
14514 for (k = l; k <= i_3; ++k) {
14516 s += u[i + k * u_dim1] * v[k + j * v_dim1];
14520 for (k = l; k <= i_3; ++k) {
14521 v[k + j * v_dim1] += s * v[k + i * v_dim1];
14528 for (j = l; j <= i_3; ++j) {
14529 v[i + j * v_dim1] = 0.;
14530 v[j + i * v_dim1] = 0.;
14535 v[i + i * v_dim1] = 1.;
14552 for (ii = 1; ii <= i_2; ++ii) {
14561 for (j = l; j <= i_3; ++j) {
14563 u[i + j * u_dim1] = 0.;
14575 for (j = l; j <= i_3; ++j) {
14579 for (k = l; k <= i_1; ++k) {
14581 s += u[k + i * u_dim1] * u[k + j * u_dim1];
14585 f = s / u[i + i * u_dim1] / g;
14588 for (k = i; k <= i_1; ++k) {
14589 u[k + j * u_dim1] += f * u[k + i * u_dim1];
14596 for (j = i; j <= i_1; ++j) {
14598 u[j + i * u_dim1] /= g;
14605 for (j = i; j <= i_1; ++j) {
14607 u[j + i * u_dim1] = 0.;
14611 u[i + i * u_dim1] += 1.;
14619 for (kk = 1; kk <= i_2; ++kk) {
14627 for (ll = 1; ll <= i_1; ++ll) {
14630 tst2 = tst1 + (d_1 = rv1[l],
abs(d_1));
14631 if (tst2 == tst1) {
14636 tst2 = tst1 + (d_1 = w[l1],
abs(d_1));
14637 if (tst2 == tst1) {
14649 for (i = l; i <= i_1; ++i) {
14651 rv1[i] = c * rv1[i];
14652 tst2 = tst1 +
abs(f);
14653 if (tst2 == tst1) {
14666 for (j = 1; j <= i_3; ++j) {
14667 y = u[j + l1 * u_dim1];
14668 z = u[j + i * u_dim1];
14669 u[j + l1 * u_dim1] = y * c + z *
s;
14670 u[j + i * u_dim1] = -y * s + z * c;
14692 f = ((g + z) / h * ((g - z) /
y) + y / h - h / y) * .5;
14694 f = x - z / x * z + h / x * (y / (f +
d_sign(&g, &f)) - h);
14700 for (i1 = l; i1 <= i_1; ++i1) {
14711 g = -x * s + g * c;
14719 for (j = 1; j <= i_3; ++j) {
14720 x = v[j + i1 * v_dim1];
14721 z = v[j + i * v_dim1];
14722 v[j + i1 * v_dim1] = x * c + z *
s;
14723 v[j + i * v_dim1] = -x * s + z * c;
14739 x = -s * g + c *
y;
14745 for (j = 1; j <= i_3; ++j) {
14746 y = u[j + i1 * u_dim1];
14747 z = u[j + i * u_dim1];
14748 u[j + i1 * u_dim1] = y * c + z *
s;
14749 u[j + i * u_dim1] = -y * s + z * c;
14773 for (j = 1; j <= i_1; ++j) {
14775 v[j + k * v_dim1] = -v[j + k * v_dim1];
14797 integer z_dim1, z_offset, i_1, i_2, i_3;
14895 z_offset = z_dim1 + 1;
14906 order = 1. - e2[1];
14913 for (q = p; q <= i_1; ++q) {
14917 if (e2[q + 1] == 0.) {
14928 for (r = 1; r <= i_1; ++r) {
14929 if (ind[r] != tag) {
14945 norm = (d_1 = d[p],
abs(d_1));
14949 for (i = ip; i <= i_2; ++i) {
14952 d_3 = norm, d_4 = (d_1 = d[i],
abs(d_1)) + (d_2 = e[i],
abs(
14954 norm =
max(d_3,d_4);
14961 eps2 = norm * .001;
14965 uk = eps4 / sqrt(uk);
14972 if ((d_1 = x1 - x0,
abs(d_1)) >= eps2) {
14976 if (order * (x1 - x0) <= 0.) {
14977 x1 = x0 + order * eps3;
14985 for (i = p; i <= i_2; ++i) {
14990 if ((d_1 = e[i],
abs(d_1)) <
abs(u)) {
14999 rv2[i - 1] = d[i] - x1;
15002 rv3[i - 1] = e[i + 1];
15004 u = v - xu * rv2[i - 1];
15005 v = -xu * rv3[i - 1];
15014 u = d[i] - x1 - xu * v;
15032 for (ii = p; ii <= i_2; ++ii) {
15034 rv6[i] = (rv6[i] - u * rv2[i] - v * rv3[i]) / rv1[i];
15047 for (jj = 1; jj <= i_2; ++jj) {
15050 if (ind[j] != tag) {
15056 for (i = p; i <= i_3; ++i) {
15058 xu += rv6[i] * z[i + j * z_dim1];
15062 for (i = p; i <= i_3; ++i) {
15064 rv6[i] -= xu * z[i + j * z_dim1];
15074 for (i = p; i <= i_2; ++i) {
15076 norm += (d_1 = rv6[i],
abs(d_1));
15099 for (i = p; i <= i_2; ++i) {
15107 for (i = ip; i <= i_2; ++i) {
15112 if (rv1[i - 1] != e[i]) {
15116 rv6[i - 1] = rv6[i];
15118 rv6[i] = u - rv4[i] * rv6[i - 1];
15135 for (i = p; i <= i_2; ++i) {
15144 for (i = 1; i <= i_2; ++i) {
15146 z[i + r * z_dim1] = 0.;
15150 for (i = p; i <= i_2; ++i) {
15152 z[i + r * z_dim1] = rv6[i] * xu;
15244 for (i = 2; i <= i_1; ++i) {
15254 for (l = 1; l <= i_1; ++l) {
15256 h = (d_1 = d[l],
abs(d_1)) + (d_2 = e[l],
abs(d_2));
15262 for (m = l; m <= i_2; ++m) {
15263 tst2 = tst1 + (d_1 = e[m],
abs(d_1));
15264 if (tst2 == tst1) {
15285 p = (d[l1] - g) / (e[l] * 2.);
15287 d[l] = e[l] / (p +
d_sign(&r, &p));
15288 d[l1] = e[l] * (p +
d_sign(&r, &p));
15296 for (i = l2; i <= i_2; ++i) {
15312 for (ii = 1; ii <= i_2; ++ii) {
15323 p = c * d[i] - s * g;
15324 d[i + 1] = h + s * (c * g + s * d[i]);
15328 p = -s * s2 * c3 * el1 * e[l] / dl1;
15331 tst2 = tst1 + (d_1 = e[l],
abs(d_1));
15343 for (ii = 2; ii <= i_2; ++ii) {
15345 if (p >= d[i - 1]) {
15372 integer z_dim1, z_offset, i_1, i_2, i_3;
15380 static integer i, j, k, l, m;
15452 z_offset = z_dim1 + 1;
15464 for (i = 2; i <= i_1; ++i) {
15474 for (l = 1; l <= i_1; ++l) {
15476 h = (d_1 = d[l],
abs(d_1)) + (d_2 = e[l],
abs(d_2));
15482 for (m = l; m <= i_2; ++m) {
15483 tst2 = tst1 + (d_1 = e[m],
abs(d_1));
15484 if (tst2 == tst1) {
15505 p = (d[l1] - g) / (e[l] * 2.);
15507 d[l] = e[l] / (p +
d_sign(&r, &p));
15508 d[l1] = e[l] * (p +
d_sign(&r, &p));
15516 for (i = l2; i <= i_2; ++i) {
15532 for (ii = 1; ii <= i_2; ++ii) {
15543 p = c * d[i] - s * g;
15544 d[i + 1] = h + s * (c * g + s * d[i]);
15547 for (k = 1; k <= i_3; ++k) {
15548 h = z[k + (i + 1) * z_dim1];
15549 z[k + (i + 1) * z_dim1] = s * z[k + i * z_dim1] + c * h;
15550 z[k + i * z_dim1] = c * z[k + i * z_dim1] - s * h;
15557 p = -s * s2 * c3 * el1 * e[l] / dl1;
15560 tst2 = tst1 + (d_1 = e[l],
abs(d_1));
15570 for (ii = 2; ii <= i_1; ++ii) {
15576 for (j = ii; j <= i_2; ++j) {
15593 for (j = 1; j <= i_2; ++j) {
15594 p = z[j + i * z_dim1];
15595 z[j + i * z_dim1] = z[j + k * z_dim1];
15596 z[j + k * z_dim1] = p;
15686 for (i = 2; i <= i_1; ++i) {
15696 for (l = 1; l <= i_1; ++l) {
15698 h = (d_1 = d[l],
abs(d_1)) + sqrt(e2[l]);
15709 for (m = l; m <= i_2; ++m) {
15731 p = (d[l1] - g) / (s * 2.);
15733 d[l] = s / (p +
d_sign(&r, &p));
15737 for (i = l1; i <= i_2; ++i) {
15753 for (ii = 1; ii <= i_2; ++ii) {
15759 d[i + 1] = h + s * (h + d[i]);
15760 g = d[i] - e2[i] / g;
15775 if ((d_1 = e2[l],
abs(d_1)) <= (d_2 = c / h,
abs(d_2))) {
15790 for (ii = 2; ii <= i_2; ++ii) {
15792 if (p >= d[i - 1]) {
15819 integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2, i_3;
15874 a_offset = a_dim1 + 1;
15877 z_offset = z_dim1 + 1;
15889 for (i = 2; i <= i_1; ++i) {
15896 for (j = 1; j <= i_2; ++j) {
15900 for (k = 1; k <= i_3; ++k) {
15902 s += a[i + k * a_dim1] * z[k + j * z_dim1];
15908 s = s / a[i + l * a_dim1] / e[i];
15911 for (k = 1; k <= i_3; ++k) {
15913 z[k + j * z_dim1] += s * a[i + k * a_dim1];
15931 integer z_dim1, z_offset, i_1, i_2, i_3;
15988 z_offset = z_dim1 + 1;
16000 for (i = 2; i <= i_1; ++i) {
16010 for (j = 1; j <= i_2; ++j) {
16015 for (k = 1; k <= i_3; ++k) {
16017 s += a[ik] * z[k + j * z_dim1];
16026 for (k = 1; k <= i_3; ++k) {
16028 z[k + j * z_dim1] -= s * a[ik];
16047 integer a_dim1, a_offset, i_1, i_2, i_3;
16108 a_offset = a_dim1 + 1;
16113 for (i = 1; i <= i_1; ++i) {
16114 d[i] = a[*n + i * a_dim1];
16115 a[*n + i * a_dim1] = a[i + i * a_dim1];
16120 for (ii = 1; ii <= i_1; ++ii) {
16130 for (k = 1; k <= i_2; ++k) {
16132 scale += (d_1 = d[k],
abs(d_1));
16140 for (j = 1; j <= i_2; ++j) {
16141 d[j] = a[l + j * a_dim1];
16142 a[l + j * a_dim1] = a[i + j * a_dim1];
16143 a[i + j * a_dim1] = 0.;
16154 for (k = 1; k <= i_2; ++k) {
16160 e2[i] = scale * scale * h;
16172 for (j = 1; j <= i_2; ++j) {
16178 for (j = 1; j <= i_2; ++j) {
16180 g = e[j] + a[j + j * a_dim1] * f;
16187 for (k = jp1; k <= i_3; ++k) {
16188 g += a[k + j * a_dim1] * d[k];
16189 e[k] += a[k + j * a_dim1] * f;
16201 for (j = 1; j <= i_2; ++j) {
16210 for (j = 1; j <= i_2; ++j) {
16216 for (j = 1; j <= i_2; ++j) {
16221 for (k = j; k <= i_3; ++k) {
16223 a[k + j * a_dim1] = a[k + j * a_dim1] - f * e[k] - g * d[k];
16231 for (j = 1; j <= i_2; ++j) {
16233 d[j] = a[l + j * a_dim1];
16234 a[l + j * a_dim1] = a[i + j * a_dim1];
16235 a[i + j * a_dim1] = f * scale;
16250 integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2, i_3;
16306 z_offset = z_dim1 + 1;
16311 a_offset = a_dim1 + 1;
16316 for (i = 1; i <= i_1; ++i) {
16319 for (j = i; j <= i_2; ++j) {
16321 z[j + i * z_dim1] = a[j + i * a_dim1];
16324 d[i] = a[*n + i * a_dim1];
16333 for (ii = 2; ii <= i_1; ++ii) {
16343 for (k = 1; k <= i_2; ++k) {
16345 scale += (d_1 = d[k],
abs(d_1));
16355 for (j = 1; j <= i_2; ++j) {
16356 d[j] = z[l + j * z_dim1];
16357 z[i + j * z_dim1] = 0.;
16358 z[j + i * z_dim1] = 0.;
16366 for (k = 1; k <= i_2; ++k) {
16380 for (j = 1; j <= i_2; ++j) {
16386 for (j = 1; j <= i_2; ++j) {
16388 z[j + i * z_dim1] = f;
16389 g = e[j] + z[j + j * z_dim1] * f;
16396 for (k = jp1; k <= i_3; ++k) {
16397 g += z[k + j * z_dim1] * d[k];
16398 e[k] += z[k + j * z_dim1] * f;
16410 for (j = 1; j <= i_2; ++j) {
16419 for (j = 1; j <= i_2; ++j) {
16425 for (j = 1; j <= i_2; ++j) {
16430 for (k = j; k <= i_3; ++k) {
16432 z[k + j * z_dim1] = z[k + j * z_dim1] - f * e[k] - g * d[k];
16435 d[j] = z[l + j * z_dim1];
16436 z[i + j * z_dim1] = 0.;
16446 for (i = 2; i <= i_1; ++i) {
16448 z[*n + l * z_dim1] = z[l + l * z_dim1];
16449 z[l + l * z_dim1] = 1.;
16456 for (k = 1; k <= i_2; ++k) {
16458 d[k] = z[k + i * z_dim1] / h;
16462 for (j = 1; j <= i_2; ++j) {
16466 for (k = 1; k <= i_3; ++k) {
16468 g += z[k + i * z_dim1] * z[k + j * z_dim1];
16472 for (k = 1; k <= i_3; ++k) {
16473 z[k + j * z_dim1] -= g * d[k];
16480 for (k = 1; k <= i_3; ++k) {
16482 z[k + i * z_dim1] = 0.;
16490 for (i = 1; i <= i_1; ++i) {
16491 d[i] = z[*n + i * z_dim1];
16492 z[*n + i * z_dim1] = 0.;
16496 z[*n + *n * z_dim1] = 1.;
16515 static integer ii, jk, iz, jm1;
16569 for (ii = 1; ii <= i_1; ++ii) {
16580 for (k = 1; k <= i_2; ++k) {
16583 scale += (d_1 = d[k],
abs(d_1));
16597 for (k = 1; k <= i_2; ++k) {
16603 e2[i] = scale * scale * h;
16610 a[iz] = scale * d[l];
16617 for (j = 1; j <= i_2; ++j) {
16626 for (k = 1; k <= i_3; ++k) {
16634 e[j] = g + a[jk] * f;
16642 for (j = 1; j <= i_2; ++j) {
16651 for (j = 1; j <= i_2; ++j) {
16659 for (j = 1; j <= i_2; ++j) {
16664 for (k = 1; k <= i_3; ++k) {
16665 a[jk] = a[jk] - f * e[k] - g * d[k];
16675 a[iz + 1] = scale * sqrt(h);
16692 static integer i, j, k, l, p, q, r,
s;
16799 for (i = 1; i <= i_1; ++i) {
16803 u = (d_1 = e[i + 1],
abs(d_1));
16806 d_1 = d[i] - (x1 + u);
16809 d_1 = d[i] + (x1 + u);
16814 tst1 = (d_1 = d[i],
abs(d_1)) + (d_2 = d[i - 1],
abs(d_2));
16815 tst2 = tst1 + (d_1 = e[i],
abs(d_1));
16827 d_2 =
abs(xu), d_3 =
abs(x0);
16828 d_1 =
max(d_2,d_3);
16845 x1 = xu + (x0 - xu) * .5;
16851 if ((i_1 = s - m1) < 0) {
16853 }
else if (i_1 == 0) {
16876 if ((i_1 = s - m22) < 0) {
16878 }
else if (i_1 == 0) {
16901 for (q = p; q <= i_1; ++q) {
16908 u = (d_1 = e[q + 1],
abs(d_1));
16912 d_1 = d[q] - (x1 + u);
16915 d_1 = d[q] + (x1 + u);
16925 d_2 =
abs(xu), d_3 =
abs(x0);
16926 d_1 =
max(d_2,d_3);
16935 if (t1 > d[p] || d[p] >= t2) {
16945 d_1 = t1, d_2 = xu - x1;
16946 *lb =
max(d_1,d_2);
16948 d_1 = t2, d_2 = x0 + x1;
16949 *ub =
min(d_1,d_2);
16968 for (i = m1; i <= i_1; ++i) {
16982 for (ii = m1; ii <= i_1; ++ii) {
16984 if (xu >= rv4[i]) {
16999 x1 = (xu + x0) * .5;
17000 if (x0 - xu <=
abs(*eps1)) {
17003 tst1 = (
abs(xu) +
abs(x0)) * 2.;
17004 tst2 = tst1 + (x0 - xu);
17005 if (tst2 == tst1) {
17014 for (i = p; i <= i_1; ++i) {
17018 v = (d_1 = e[i],
abs(d_1)) /
epslon_(&c_b141);
17071 r = r + m2 - m1 + 1;
17076 for (l = 1; l <= i_1; ++l) {
17083 if (rv5[k] >= w[l]) {
17088 for (ii = j; ii <= i_2; ++ii) {
17091 ind[i + 1] = ind[i];
17114 *ierr = *n * 3 + isturm;
17128 integer z_dim1, z_offset, i_1, i_2, i_3;
17136 static integer i, j, k, p, q, r,
s;
17145 static doublereal eps2, eps3, eps4, tst1, tst2;
17247 z_offset = z_dim1 + 1;
17257 for (i = 1; i <= i_1; ++i) {
17261 tst1 = (d_1 = d[i],
abs(d_1)) + (d_2 = d[i - 1],
abs(d_2));
17262 tst2 = tst1 + (d_1 = e[i],
abs(d_1));
17302 for (q = p; q <= i_1; ++q) {
17309 u = (d_1 = e[q + 1],
abs(d_1));
17313 d_1 = d[q] - (x1 + u);
17316 d_1 = d[q] + (x1 + u);
17326 d_2 =
abs(xu), d_3 =
abs(x0);
17327 d_1 =
max(d_2,d_3);
17336 if (t1 > d[p] || d[p] >= t2) {
17342 for (i = 1; i <= i_1; ++i) {
17344 z[i + r * z_dim1] = 0.;
17348 z[p + r * z_dim1] = 1.;
17354 d_1 = t1, d_2 = xu - x1;
17355 *lb =
max(d_1,d_2);
17357 d_1 = t2, d_2 = x0 + x1;
17358 *ub =
min(d_1,d_2);
17377 for (i = m1; i <= i_1; ++i) {
17391 for (ii = m1; ii <= i_1; ++ii) {
17393 if (xu >= rv4[i]) {
17408 x1 = (xu + x0) * .5;
17409 if (x0 - xu <=
abs(*eps1)) {
17412 tst1 = (
abs(xu) +
abs(x0)) * 2.;
17413 tst2 = tst1 + (x0 - xu);
17414 if (tst2 == tst1) {
17423 for (i = p; i <= i_1; ++i) {
17427 v = (d_1 = e[i],
abs(d_1)) /
epslon_(&c_b141);
17477 norm = (d_1 = d[p],
abs(d_1));
17481 for (i = ip; i <= i_1; ++i) {
17484 d_3 = norm, d_4 = (d_1 = d[i],
abs(d_1)) + (d_2 = e[i],
abs(d_2)
17486 norm =
max(d_3,d_4);
17492 eps2 = norm * .001;
17496 uk = eps4 / sqrt(uk);
17501 for (k = m1; k <= i_1; ++k) {
17510 if (x1 - x0 >= eps2) {
17523 for (i = p; i <= i_2; ++i) {
17528 if ((d_1 = e[i],
abs(d_1)) <
abs(u)) {
17534 rv2[i - 1] = d[i] - x1;
17537 rv3[i - 1] = e[i + 1];
17539 u = v - xu * rv2[i - 1];
17540 v = -xu * rv3[i - 1];
17549 u = d[i] - x1 - xu * v;
17567 for (ii = p; ii <= i_2; ++ii) {
17569 rv6[i] = (rv6[i] - u * rv2[i] - v * rv3[i]) / rv1[i];
17581 for (jj = 1; jj <= i_2; ++jj) {
17582 j = r - group - 1 + jj;
17586 for (i = p; i <= i_3; ++i) {
17588 xu += rv6[i] * z[i + j * z_dim1];
17592 for (i = p; i <= i_3; ++i) {
17594 rv6[i] -= xu * z[i + j * z_dim1];
17604 for (i = p; i <= i_2; ++i) {
17606 norm += (d_1 = rv6[i],
abs(d_1));
17629 for (i = p; i <= i_2; ++i) {
17637 for (i = ip; i <= i_2; ++i) {
17642 if (rv1[i - 1] != e[i]) {
17646 rv6[i - 1] = rv6[i];
17648 rv6[i] = u - rv4[i] * rv6[i - 1];
17660 for (i = p; i <= i_2; ++i) {
17668 for (i = 1; i <= i_2; ++i) {
17670 z[i + r * z_dim1] = 0.;
17674 for (i = p; i <= i_2; ++i) {
17676 z[i + r * z_dim1] = rv6[i] * xu;
17690 *ierr = (*n << 2) + r;
17695 *ierr = *n * 3 + 1;
int tql2_(integer *nm, integer *n, doublereal *d, doublereal *e, doublereal *z, integer *ierr)
int rg_(integer *nm, integer *n, doublereal *a, doublereal *wr, doublereal *wi, integer *matz, doublereal *z, integer *iv1, doublereal *fv1, integer *ierr)
int rt_(integer *nm, integer *n, doublereal *a, doublereal *w, integer *matz, doublereal *z, doublereal *fv1, integer *ierr)
int qzhes_(integer *nm, integer *n, doublereal *a, doublereal *b, logical *matz, doublereal *z)
int svd_(integer *nm, integer *m, integer *n, doublereal *a, doublereal *w, logical *matu, doublereal *u, logical *matv, doublereal *v, integer *ierr, doublereal *rv1)
int eltran_(integer *nm, integer *n, integer *low, integer *igh, doublereal *a, integer *int_, doublereal *z)
int trbak1_(integer *nm, integer *n, doublereal *a, doublereal *e, integer *m, doublereal *z)
int comqr2_(integer *nm, integer *n, integer *low, integer *igh, doublereal *ortr, doublereal *orti, doublereal *hr, doublereal *hi, doublereal *wr, doublereal *wi, doublereal *zr, doublereal *zi, integer *ierr)
int rst_(integer *nm, integer *n, doublereal *w, doublereal *e, integer *matz, doublereal *z, integer *ierr)
int rsp_(integer *nm, integer *n, integer *nv, doublereal *a, doublereal *w, integer *matz, doublereal *z, doublereal *fv1, doublereal *fv2, integer *ierr)
int cbal_(integer *nm, integer *n, doublereal *ar, doublereal *ai, integer *low, integer *igh, doublereal *scale)
int imtqlv_(integer *n, doublereal *d, doublereal *e, doublereal *e2, doublereal *w, integer *ind, integer *ierr, doublereal *rv1)
int rebak_(integer *nm, integer *n, doublereal *b, doublereal *dl, integer *m, doublereal *z)
int balbak_(integer *nm, integer *n, integer *low, integer *igh, doublereal *scale, integer *m, doublereal *z)
int bandr_(integer *nm, integer *n, integer *mb, doublereal *a, doublereal *d, doublereal *e, doublereal *e2, logical *matz, doublereal *z)
int figi2_(integer *nm, integer *n, doublereal *t, doublereal *d, doublereal *e, doublereal *z, integer *ierr)
int comqr_(integer *nm, integer *n, integer *low, integer *igh, doublereal *hr, doublereal *hi, doublereal *wr, doublereal *wi, integer *ierr)
int cortb_(integer *nm, integer *low, integer *igh, doublereal *ar, doublereal *ai, doublereal *ortr, doublereal *orti, integer *m, doublereal *zr, doublereal *zi)
int rs_(integer *nm, integer *n, doublereal *a, doublereal *w, integer *matz, doublereal *z, doublereal *fv1, doublereal *fv2, integer *ierr)
int tql1_(integer *n, doublereal *d, doublereal *e, integer *ierr)
int bakvec_(integer *nm, integer *n, doublereal *t, doublereal *e, integer *m, doublereal *z, integer *ierr)
int cbabk2_(integer *nm, integer *n, integer *low, integer *igh, doublereal *scale, integer *m, doublereal *zr, doublereal *zi)
int hqr_(integer *nm, integer *n, integer *low, integer *igh, doublereal *h, doublereal *wr, doublereal *wi, integer *ierr)
int comlr_(integer *nm, integer *n, integer *low, integer *igh, doublereal *hr, doublereal *hi, doublereal *wr, doublereal *wi, integer *ierr)
int qzval_(integer *nm, integer *n, doublereal *a, doublereal *b, doublereal *alfr, doublereal *alfi, doublereal *beta, logical *matz, doublereal *z)
int bqr_(integer *nm, integer *n, integer *mb, doublereal *a, doublereal *t, doublereal *r, integer *ierr, integer *, doublereal *rv)
int tred3_(integer *n, integer *, doublereal *a, doublereal *d, doublereal *e, doublereal *e2)
int csroot_(doublereal *xr, doublereal *xi, doublereal *yr, doublereal *yi)
int bisect_(integer *n, doublereal *eps1, doublereal *d, doublereal *e, doublereal *e2, doublereal *lb, doublereal *ub, integer *mm, integer *m, doublereal *w, integer *ind, integer *ierr, doublereal *rv4, doublereal *rv5)
int tridib_(integer *n, doublereal *eps1, doublereal *d, doublereal *e, doublereal *e2, doublereal *lb, doublereal *ub, integer *m11, integer *m, doublereal *w, integer *ind, integer *ierr, doublereal *rv4, doublereal *rv5)
int balanc_(integer *nm, integer *n, doublereal *a, integer *low, integer *igh, doublereal *scale)
int tinvit_(integer *nm, integer *n, doublereal *d, doublereal *e, doublereal *e2, integer *m, doublereal *w, integer *ind, doublereal *z, integer *ierr, doublereal *rv1, doublereal *rv2, doublereal *rv3, doublereal *rv4, doublereal *rv6)
int elmhes_(integer *nm, integer *n, integer *low, integer *igh, doublereal *a, integer *int_)
int orthes_(integer *nm, integer *n, integer *low, integer *igh, doublereal *a, doublereal *ort)
int ratqr_(integer *n, doublereal *eps1, doublereal *d, doublereal *e, doublereal *e2, integer *m, doublereal *w, integer *ind, doublereal *bd, logical *type, integer *idef, integer *ierr)
doublereal epslon_(doublereal *x)
int trbak3_(integer *nm, integer *n, integer *, doublereal *a, integer *m, doublereal *z)
int ch_(integer *nm, integer *n, doublereal *ar, doublereal *ai, doublereal *w, integer *matz, doublereal *zr, doublereal *zi, doublereal *fv1, doublereal *fv2, doublereal *fm1, integer *ierr)
int ortran_(integer *nm, integer *n, integer *low, integer *igh, doublereal *a, doublereal *ort, doublereal *z)
int htribk_(integer *nm, integer *n, doublereal *ar, doublereal *ai, doublereal *tau, integer *m, doublereal *zr, doublereal *zi)
int cinvit_(integer *nm, integer *n, doublereal *ar, doublereal *ai, doublereal *wr, doublereal *wi, logical *select, integer *mm, integer *m, doublereal *zr, doublereal *zi, integer *ierr, doublereal *rm1, doublereal *rm2, doublereal *rv1, doublereal *rv2)
int bandv_(integer *nm, integer *n, integer *mbw, doublereal *a, doublereal *e21, integer *m, doublereal *w, doublereal *z, integer *ierr, integer *, doublereal *rv, doublereal *rv6)
int rsb_(integer *nm, integer *n, integer *mb, doublereal *a, doublereal *w, integer *matz, doublereal *z, doublereal *fv1, doublereal *fv2, integer *ierr)
int rgg_(integer *nm, integer *n, doublereal *a, doublereal *b, doublereal *alfr, doublereal *alfi, doublereal *beta, integer *matz, doublereal *z, integer *ierr)
int minfit_(integer *nm, integer *m, integer *n, doublereal *a, doublereal *w, integer *ip, doublereal *b, integer *ierr, doublereal *rv1)
int htrid3_(integer *nm, integer *n, doublereal *a, doublereal *d, doublereal *e, doublereal *e2, doublereal *tau)
int hqr2_(integer *nm, integer *n, integer *low, integer *igh, doublereal *h, doublereal *wr, doublereal *wi, doublereal *z, integer *ierr)
int cdiv_(doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, doublereal *cr, doublereal *ci)
int figi_(integer *nm, integer *n, doublereal *t, doublereal *d, doublereal *e, doublereal *e2, integer *ierr)
int tsturm_(integer *nm, integer *n, doublereal *eps1, doublereal *d, doublereal *e, doublereal *e2, doublereal *lb, doublereal *ub, integer *mm, integer *m, doublereal *w, doublereal *z, integer *ierr, doublereal *rv1, doublereal *rv2, doublereal *rv3, doublereal *rv4, doublereal *rv5, doublereal *rv6)
double d_sign(doublereal *a, doublereal *b)
int qzvec_(integer *nm, integer *n, doublereal *a, doublereal *b, doublereal *alfr, doublereal *alfi, doublereal *beta, doublereal *z)
int tqlrat_(integer *n, doublereal *d, doublereal *e2, integer *ierr)
int reduc_(integer *nm, integer *n, doublereal *a, doublereal *b, doublereal *dl, integer *ierr)
int htridi_(integer *nm, integer *n, doublereal *ar, doublereal *ai, doublereal *d, doublereal *e, doublereal *e2, doublereal *tau)
int ortbak_(integer *nm, integer *low, integer *igh, doublereal *a, doublereal *ort, integer *m, doublereal *z)
int rsgba_(integer *nm, integer *n, doublereal *a, doublereal *b, doublereal *w, integer *matz, doublereal *z, doublereal *fv1, doublereal *fv2, integer *ierr)
int corth_(integer *nm, integer *n, integer *low, integer *igh, doublereal *ar, doublereal *ai, doublereal *ortr, doublereal *orti)
int qzit_(integer *nm, integer *n, doublereal *a, doublereal *b, doublereal *eps1, logical *matz, doublereal *z, integer *ierr)
int combak_(integer *nm, integer *low, integer *igh, doublereal *ar, doublereal *ai, integer *int_, integer *m, doublereal *zr, doublereal *zi)
int rebakb_(integer *nm, integer *n, doublereal *b, doublereal *dl, integer *m, doublereal *z)
int comlr2_(integer *nm, integer *n, integer *low, integer *igh, integer *int_, doublereal *hr, doublereal *hi, doublereal *wr, doublereal *wi, doublereal *zr, doublereal *zi, integer *ierr)
int htrib3_(integer *nm, integer *n, doublereal *a, doublereal *tau, integer *m, doublereal *zr, doublereal *zi)
int rsm_(integer *nm, integer *n, doublereal *a, doublereal *w, integer *m, doublereal *z, doublereal *fwork, integer *iwork, integer *ierr)
int rsgab_(integer *nm, integer *n, doublereal *a, doublereal *b, doublereal *w, integer *matz, doublereal *z, doublereal *fv1, doublereal *fv2, integer *ierr)
int elmbak_(integer *nm, integer *low, integer *igh, doublereal *a, integer *int_, integer *m, doublereal *z)
int cg_(integer *nm, integer *n, doublereal *ar, doublereal *ai, doublereal *wr, doublereal *wi, integer *matz, doublereal *zr, doublereal *zi, doublereal *fv1, doublereal *fv2, doublereal *fv3, integer *ierr)
int tred1_(integer *nm, integer *n, doublereal *a, doublereal *d, doublereal *e, doublereal *e2)
int imtql1_(integer *n, doublereal *d, doublereal *e, integer *ierr)
doublereal pythag_(doublereal *a, doublereal *b)
int tred2_(integer *nm, integer *n, doublereal *a, doublereal *d, doublereal *e, doublereal *z)
int invit_(integer *nm, integer *n, doublereal *a, doublereal *wr, doublereal *wi, logical *select, integer *mm, integer *m, doublereal *z, integer *ierr, doublereal *rm1, doublereal *rv1, doublereal *rv2)
int comhes_(integer *nm, integer *n, integer *low, integer *igh, doublereal *ar, doublereal *ai, integer *int_)
int imtql2_(integer *nm, integer *n, doublereal *d, doublereal *e, doublereal *z, integer *ierr)
int reduc2_(integer *nm, integer *n, doublereal *a, doublereal *b, doublereal *dl, integer *ierr)
int rsg_(integer *nm, integer *n, doublereal *a, doublereal *b, doublereal *w, integer *matz, doublereal *z, doublereal *fv1, doublereal *fv2, integer *ierr)
GB_write_int const char s