Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_tinvit.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008 int tinvit_(integer *nm, integer *n, doublereal *d__,
00009 doublereal *e, doublereal *e2, integer *m, doublereal *w, integer *
00010 ind, doublereal *z__, integer *ierr, doublereal *rv1, doublereal *rv2,
00011 doublereal *rv3, doublereal *rv4, doublereal *rv6)
00012 {
00013
00014 integer z_dim1, z_offset, i__1, i__2, i__3;
00015 doublereal d__1, d__2, d__3, d__4;
00016
00017
00018 double sqrt(doublereal);
00019
00020
00021 static doublereal norm;
00022 static integer i__, j, p, q, r__, s;
00023 static doublereal u, v, order;
00024 static integer group;
00025 static doublereal x0, x1;
00026 static integer ii, jj, ip;
00027 static doublereal uk, xu;
00028 extern doublereal pythag_(doublereal *, doublereal *), epslon_(doublereal
00029 *);
00030 static integer tag, its;
00031 static doublereal eps2, eps3, eps4;
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
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103 --rv6;
00104 --rv4;
00105 --rv3;
00106 --rv2;
00107 --rv1;
00108 --e2;
00109 --e;
00110 --d__;
00111 z_dim1 = *nm;
00112 z_offset = z_dim1 + 1;
00113 z__ -= z_offset;
00114 --ind;
00115 --w;
00116
00117
00118 *ierr = 0;
00119 if (*m == 0) {
00120 goto L1001;
00121 }
00122 tag = 0;
00123 order = 1. - e2[1];
00124 q = 0;
00125
00126 L100:
00127 p = q + 1;
00128
00129 i__1 = *n;
00130 for (q = p; q <= i__1; ++q) {
00131 if (q == *n) {
00132 goto L140;
00133 }
00134 if (e2[q + 1] == 0.) {
00135 goto L140;
00136 }
00137
00138 }
00139
00140 L140:
00141 ++tag;
00142 s = 0;
00143
00144 i__1 = *m;
00145 for (r__ = 1; r__ <= i__1; ++r__) {
00146 if (ind[r__] != tag) {
00147 goto L920;
00148 }
00149 its = 1;
00150 x1 = w[r__];
00151 if (s != 0) {
00152 goto L510;
00153 }
00154
00155 xu = 1.;
00156 if (p != q) {
00157 goto L490;
00158 }
00159 rv6[p] = 1.;
00160 goto L870;
00161 L490:
00162 norm = (d__1 = d__[p], abs(d__1));
00163 ip = p + 1;
00164
00165 i__2 = q;
00166 for (i__ = ip; i__ <= i__2; ++i__) {
00167
00168
00169 d__3 = norm, d__4 = (d__1 = d__[i__], abs(d__1)) + (d__2 = e[i__],
00170 abs(d__2));
00171 norm = max(d__3,d__4);
00172 }
00173
00174
00175
00176
00177
00178 eps2 = norm * .001;
00179 eps3 = epslon_(&norm);
00180 uk = (doublereal) (q - p + 1);
00181 eps4 = uk * eps3;
00182 uk = eps4 / sqrt(uk);
00183 s = p;
00184 L505:
00185 group = 0;
00186 goto L520;
00187
00188 L510:
00189 if ((d__1 = x1 - x0, abs(d__1)) >= eps2) {
00190 goto L505;
00191 }
00192 ++group;
00193 if (order * (x1 - x0) <= 0.) {
00194 x1 = x0 + order * eps3;
00195 }
00196
00197
00198 L520:
00199 v = 0.;
00200
00201 i__2 = q;
00202 for (i__ = p; i__ <= i__2; ++i__) {
00203 rv6[i__] = uk;
00204 if (i__ == p) {
00205 goto L560;
00206 }
00207 if ((d__1 = e[i__], abs(d__1)) < abs(u)) {
00208 goto L540;
00209 }
00210
00211
00212
00213 xu = u / e[i__];
00214 rv4[i__] = xu;
00215 rv1[i__ - 1] = e[i__];
00216 rv2[i__ - 1] = d__[i__] - x1;
00217 rv3[i__ - 1] = 0.;
00218 if (i__ != q) {
00219 rv3[i__ - 1] = e[i__ + 1];
00220 }
00221 u = v - xu * rv2[i__ - 1];
00222 v = -xu * rv3[i__ - 1];
00223 goto L580;
00224 L540:
00225 xu = e[i__] / u;
00226 rv4[i__] = xu;
00227 rv1[i__ - 1] = u;
00228 rv2[i__ - 1] = v;
00229 rv3[i__ - 1] = 0.;
00230 L560:
00231 u = d__[i__] - x1 - xu * v;
00232 if (i__ != q) {
00233 v = e[i__ + 1];
00234 }
00235 L580:
00236 ;
00237 }
00238
00239 if (u == 0.) {
00240 u = eps3;
00241 }
00242 rv1[q] = u;
00243 rv2[q] = 0.;
00244 rv3[q] = 0.;
00245
00246
00247 L600:
00248 i__2 = q;
00249 for (ii = p; ii <= i__2; ++ii) {
00250 i__ = p + q - ii;
00251 rv6[i__] = (rv6[i__] - u * rv2[i__] - v * rv3[i__]) / rv1[i__];
00252 v = u;
00253 u = rv6[i__];
00254
00255 }
00256
00257
00258 if (group == 0) {
00259 goto L700;
00260 }
00261 j = r__;
00262
00263 i__2 = group;
00264 for (jj = 1; jj <= i__2; ++jj) {
00265 L630:
00266 --j;
00267 if (ind[j] != tag) {
00268 goto L630;
00269 }
00270 xu = 0.;
00271
00272 i__3 = q;
00273 for (i__ = p; i__ <= i__3; ++i__) {
00274
00275 xu += rv6[i__] * z__[i__ + j * z_dim1];
00276 }
00277
00278 i__3 = q;
00279 for (i__ = p; i__ <= i__3; ++i__) {
00280
00281 rv6[i__] -= xu * z__[i__ + j * z_dim1];
00282 }
00283
00284
00285 }
00286
00287 L700:
00288 norm = 0.;
00289
00290 i__2 = q;
00291 for (i__ = p; i__ <= i__2; ++i__) {
00292
00293 norm += (d__1 = rv6[i__], abs(d__1));
00294 }
00295
00296 if (norm >= 1.) {
00297 goto L840;
00298 }
00299
00300 if (its == 5) {
00301 goto L830;
00302 }
00303 if (norm != 0.) {
00304 goto L740;
00305 }
00306 rv6[s] = eps4;
00307 ++s;
00308 if (s > q) {
00309 s = p;
00310 }
00311 goto L780;
00312 L740:
00313 xu = eps4 / norm;
00314
00315 i__2 = q;
00316 for (i__ = p; i__ <= i__2; ++i__) {
00317
00318 rv6[i__] *= xu;
00319 }
00320
00321
00322 L780:
00323 i__2 = q;
00324 for (i__ = ip; i__ <= i__2; ++i__) {
00325 u = rv6[i__];
00326
00327
00328
00329 if (rv1[i__ - 1] != e[i__]) {
00330 goto L800;
00331 }
00332 u = rv6[i__ - 1];
00333 rv6[i__ - 1] = rv6[i__];
00334 L800:
00335 rv6[i__] = u - rv4[i__] * rv6[i__ - 1];
00336
00337 }
00338
00339 ++its;
00340 goto L600;
00341
00342 L830:
00343 *ierr = -r__;
00344 xu = 0.;
00345 goto L870;
00346
00347
00348 L840:
00349 u = 0.;
00350
00351 i__2 = q;
00352 for (i__ = p; i__ <= i__2; ++i__) {
00353
00354 u = pythag_(&u, &rv6[i__]);
00355 }
00356
00357 xu = 1. / u;
00358
00359 L870:
00360 i__2 = *n;
00361 for (i__ = 1; i__ <= i__2; ++i__) {
00362
00363 z__[i__ + r__ * z_dim1] = 0.;
00364 }
00365
00366 i__2 = q;
00367 for (i__ = p; i__ <= i__2; ++i__) {
00368
00369 z__[i__ + r__ * z_dim1] = rv6[i__] * xu;
00370 }
00371
00372 x0 = x1;
00373 L920:
00374 ;
00375 }
00376
00377 if (q < *n) {
00378 goto L100;
00379 }
00380 L1001:
00381 return 0;
00382 }
00383