Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_tqlrat.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008
00009
00010 static doublereal c_b11 = 1.;
00011
00012 int tqlrat_(integer *n, doublereal *d__, doublereal *e2,
00013 integer *ierr)
00014 {
00015
00016 integer i__1, i__2;
00017 doublereal d__1, d__2;
00018
00019
00020 double d_sign(doublereal *, doublereal *);
00021
00022
00023 static doublereal b, c__, f, g, h__;
00024 static integer i__, j, l, m;
00025 static doublereal p, r__, s, t;
00026 static integer l1, ii;
00027 extern doublereal pythag_(doublereal *, doublereal *), epslon_(doublereal
00028 *);
00029 static integer mml;
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075 --e2;
00076 --d__;
00077
00078
00079 *ierr = 0;
00080 if (*n == 1) {
00081 goto L1001;
00082 }
00083
00084 i__1 = *n;
00085 for (i__ = 2; i__ <= i__1; ++i__) {
00086
00087 e2[i__ - 1] = e2[i__];
00088 }
00089
00090 f = 0.;
00091 t = 0.;
00092 e2[*n] = 0.;
00093
00094 i__1 = *n;
00095 for (l = 1; l <= i__1; ++l) {
00096 j = 0;
00097 h__ = (d__1 = d__[l], abs(d__1)) + sqrt(e2[l]);
00098 if (t > h__) {
00099 goto L105;
00100 }
00101 t = h__;
00102 b = epslon_(&t);
00103 c__ = b * b;
00104
00105
00106 L105:
00107 i__2 = *n;
00108 for (m = l; m <= i__2; ++m) {
00109 if (e2[m] <= c__) {
00110 goto L120;
00111 }
00112
00113
00114
00115 }
00116
00117 L120:
00118 if (m == l) {
00119 goto L210;
00120 }
00121 L130:
00122 if (j == 30) {
00123 goto L1000;
00124 }
00125 ++j;
00126
00127 l1 = l + 1;
00128 s = sqrt(e2[l]);
00129 g = d__[l];
00130 p = (d__[l1] - g) / (s * 2.);
00131 r__ = pythag_(&p, &c_b11);
00132 d__[l] = s / (p + d_sign(&r__, &p));
00133 h__ = g - d__[l];
00134
00135 i__2 = *n;
00136 for (i__ = l1; i__ <= i__2; ++i__) {
00137
00138 d__[i__] -= h__;
00139 }
00140
00141 f += h__;
00142
00143 g = d__[m];
00144 if (g == 0.) {
00145 g = b;
00146 }
00147 h__ = g;
00148 s = 0.;
00149 mml = m - l;
00150
00151 i__2 = mml;
00152 for (ii = 1; ii <= i__2; ++ii) {
00153 i__ = m - ii;
00154 p = g * h__;
00155 r__ = p + e2[i__];
00156 e2[i__ + 1] = s * r__;
00157 s = e2[i__] / r__;
00158 d__[i__ + 1] = h__ + s * (h__ + d__[i__]);
00159 g = d__[i__] - e2[i__] / g;
00160 if (g == 0.) {
00161 g = b;
00162 }
00163 h__ = g * p / r__;
00164
00165 }
00166
00167 e2[l] = s * g;
00168 d__[l] = h__;
00169
00170
00171 if (h__ == 0.) {
00172 goto L210;
00173 }
00174 if ((d__1 = e2[l], abs(d__1)) <= (d__2 = c__ / h__, abs(d__2))) {
00175 goto L210;
00176 }
00177 e2[l] = h__ * e2[l];
00178 if (e2[l] != 0.) {
00179 goto L130;
00180 }
00181 L210:
00182 p = d__[l] + f;
00183
00184 if (l == 1) {
00185 goto L250;
00186 }
00187
00188 i__2 = l;
00189 for (ii = 2; ii <= i__2; ++ii) {
00190 i__ = l + 2 - ii;
00191 if (p >= d__[i__ - 1]) {
00192 goto L270;
00193 }
00194 d__[i__] = d__[i__ - 1];
00195
00196 }
00197
00198 L250:
00199 i__ = 1;
00200 L270:
00201 d__[i__] = p;
00202
00203 }
00204
00205 goto L1001;
00206
00207
00208 L1000:
00209 *ierr = l;
00210 L1001:
00211 return 0;
00212 }
00213