Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_tql1.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008
00009
00010 static doublereal c_b10 = 1.;
00011
00012 int tql1_(integer *n, doublereal *d__, doublereal *e,
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 c__, f, g, h__;
00024 static integer i__, j, l, m;
00025 static doublereal p, r__, s, c2, c3;
00026 static integer l1, l2;
00027 static doublereal s2;
00028 static integer ii;
00029 extern doublereal pythag_(doublereal *, doublereal *);
00030 static doublereal dl1, el1;
00031 static integer mml;
00032 static doublereal tst1, tst2;
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
00076
00077
00078
00079 --e;
00080 --d__;
00081
00082
00083 *ierr = 0;
00084 if (*n == 1) {
00085 goto L1001;
00086 }
00087
00088 i__1 = *n;
00089 for (i__ = 2; i__ <= i__1; ++i__) {
00090
00091 e[i__ - 1] = e[i__];
00092 }
00093
00094 f = 0.;
00095 tst1 = 0.;
00096 e[*n] = 0.;
00097
00098 i__1 = *n;
00099 for (l = 1; l <= i__1; ++l) {
00100 j = 0;
00101 h__ = (d__1 = d__[l], abs(d__1)) + (d__2 = e[l], abs(d__2));
00102 if (tst1 < h__) {
00103 tst1 = h__;
00104 }
00105
00106 i__2 = *n;
00107 for (m = l; m <= i__2; ++m) {
00108 tst2 = tst1 + (d__1 = e[m], abs(d__1));
00109 if (tst2 == tst1) {
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 l2 = l1 + 1;
00129 g = d__[l];
00130 p = (d__[l1] - g) / (e[l] * 2.);
00131 r__ = pythag_(&p, &c_b10);
00132 d__[l] = e[l] / (p + d_sign(&r__, &p));
00133 d__[l1] = e[l] * (p + d_sign(&r__, &p));
00134 dl1 = d__[l1];
00135 h__ = g - d__[l];
00136 if (l2 > *n) {
00137 goto L145;
00138 }
00139
00140 i__2 = *n;
00141 for (i__ = l2; i__ <= i__2; ++i__) {
00142
00143 d__[i__] -= h__;
00144 }
00145
00146 L145:
00147 f += h__;
00148
00149 p = d__[m];
00150 c__ = 1.;
00151 c2 = c__;
00152 el1 = e[l1];
00153 s = 0.;
00154 mml = m - l;
00155
00156 i__2 = mml;
00157 for (ii = 1; ii <= i__2; ++ii) {
00158 c3 = c2;
00159 c2 = c__;
00160 s2 = s;
00161 i__ = m - ii;
00162 g = c__ * e[i__];
00163 h__ = c__ * p;
00164 r__ = pythag_(&p, &e[i__]);
00165 e[i__ + 1] = s * r__;
00166 s = e[i__] / r__;
00167 c__ = p / r__;
00168 p = c__ * d__[i__] - s * g;
00169 d__[i__ + 1] = h__ + s * (c__ * g + s * d__[i__]);
00170
00171 }
00172
00173 p = -s * s2 * c3 * el1 * e[l] / dl1;
00174 e[l] = s * p;
00175 d__[l] = c__ * p;
00176 tst2 = tst1 + (d__1 = e[l], abs(d__1));
00177 if (tst2 > tst1) {
00178 goto L130;
00179 }
00180 L210:
00181 p = d__[l] + f;
00182
00183 if (l == 1) {
00184 goto L250;
00185 }
00186
00187 i__2 = l;
00188 for (ii = 2; ii <= i__2; ++ii) {
00189 i__ = l + 2 - ii;
00190 if (p >= d__[i__ - 1]) {
00191 goto L270;
00192 }
00193 d__[i__] = d__[i__ - 1];
00194
00195 }
00196
00197 L250:
00198 i__ = 1;
00199 L270:
00200 d__[i__] = p;
00201
00202 }
00203
00204 goto L1001;
00205
00206
00207 L1000:
00208 *ierr = l;
00209 L1001:
00210 return 0;
00211 }
00212