Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
eis_minfit.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include "f2c.h"
00007
00008
00009
00010 static doublereal c_b39 = 1.;
00011
00012 int minfit_(integer *nm, integer *m, integer *n, doublereal *
00013 a, doublereal *w, integer *ip, doublereal *b, integer *ierr,
00014 doublereal *rv1)
00015 {
00016
00017 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
00018 doublereal d__1, d__2, d__3, d__4;
00019
00020
00021 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
00022
00023
00024 static doublereal c__, f, g, h__;
00025 static integer i__, j, k, l;
00026 static doublereal s, x, y, z__, scale;
00027 static integer i1, k1, l1, m1, ii, kk, ll;
00028 extern doublereal pythag_(doublereal *, doublereal *);
00029 static integer its;
00030 static doublereal tst1, tst2;
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
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 --rv1;
00102 --w;
00103 a_dim1 = *nm;
00104 a_offset = a_dim1 + 1;
00105 a -= a_offset;
00106 b_dim1 = *nm;
00107 b_offset = b_dim1 + 1;
00108 b -= b_offset;
00109
00110
00111 *ierr = 0;
00112
00113 g = 0.;
00114 scale = 0.;
00115 x = 0.;
00116
00117 i__1 = *n;
00118 for (i__ = 1; i__ <= i__1; ++i__) {
00119 l = i__ + 1;
00120 rv1[i__] = scale * g;
00121 g = 0.;
00122 s = 0.;
00123 scale = 0.;
00124 if (i__ > *m) {
00125 goto L210;
00126 }
00127
00128 i__2 = *m;
00129 for (k = i__; k <= i__2; ++k) {
00130
00131 scale += (d__1 = a[k + i__ * a_dim1], abs(d__1));
00132 }
00133
00134 if (scale == 0.) {
00135 goto L210;
00136 }
00137
00138 i__2 = *m;
00139 for (k = i__; k <= i__2; ++k) {
00140 a[k + i__ * a_dim1] /= scale;
00141
00142 d__1 = a[k + i__ * a_dim1];
00143 s += d__1 * d__1;
00144
00145 }
00146
00147 f = a[i__ + i__ * a_dim1];
00148 d__1 = sqrt(s);
00149 g = -d_sign(&d__1, &f);
00150 h__ = f * g - s;
00151 a[i__ + i__ * a_dim1] = f - g;
00152 if (i__ == *n) {
00153 goto L160;
00154 }
00155
00156 i__2 = *n;
00157 for (j = l; j <= i__2; ++j) {
00158 s = 0.;
00159
00160 i__3 = *m;
00161 for (k = i__; k <= i__3; ++k) {
00162
00163 s += a[k + i__ * a_dim1] * a[k + j * a_dim1];
00164 }
00165
00166 f = s / h__;
00167
00168 i__3 = *m;
00169 for (k = i__; k <= i__3; ++k) {
00170 a[k + j * a_dim1] += f * a[k + i__ * a_dim1];
00171
00172 }
00173 }
00174
00175 L160:
00176 if (*ip == 0) {
00177 goto L190;
00178 }
00179
00180 i__3 = *ip;
00181 for (j = 1; j <= i__3; ++j) {
00182 s = 0.;
00183
00184 i__2 = *m;
00185 for (k = i__; k <= i__2; ++k) {
00186
00187 s += a[k + i__ * a_dim1] * b[k + j * b_dim1];
00188 }
00189
00190 f = s / h__;
00191
00192 i__2 = *m;
00193 for (k = i__; k <= i__2; ++k) {
00194 b[k + j * b_dim1] += f * a[k + i__ * a_dim1];
00195
00196 }
00197 }
00198
00199 L190:
00200 i__2 = *m;
00201 for (k = i__; k <= i__2; ++k) {
00202
00203 a[k + i__ * a_dim1] = scale * a[k + i__ * a_dim1];
00204 }
00205
00206 L210:
00207 w[i__] = scale * g;
00208 g = 0.;
00209 s = 0.;
00210 scale = 0.;
00211 if (i__ > *m || i__ == *n) {
00212 goto L290;
00213 }
00214
00215 i__2 = *n;
00216 for (k = l; k <= i__2; ++k) {
00217
00218 scale += (d__1 = a[i__ + k * a_dim1], abs(d__1));
00219 }
00220
00221 if (scale == 0.) {
00222 goto L290;
00223 }
00224
00225 i__2 = *n;
00226 for (k = l; k <= i__2; ++k) {
00227 a[i__ + k * a_dim1] /= scale;
00228
00229 d__1 = a[i__ + k * a_dim1];
00230 s += d__1 * d__1;
00231
00232 }
00233
00234 f = a[i__ + l * a_dim1];
00235 d__1 = sqrt(s);
00236 g = -d_sign(&d__1, &f);
00237 h__ = f * g - s;
00238 a[i__ + l * a_dim1] = f - g;
00239
00240 i__2 = *n;
00241 for (k = l; k <= i__2; ++k) {
00242
00243 rv1[k] = a[i__ + k * a_dim1] / h__;
00244 }
00245
00246 if (i__ == *m) {
00247 goto L270;
00248 }
00249
00250 i__2 = *m;
00251 for (j = l; j <= i__2; ++j) {
00252 s = 0.;
00253
00254 i__3 = *n;
00255 for (k = l; k <= i__3; ++k) {
00256
00257 s += a[j + k * a_dim1] * a[i__ + k * a_dim1];
00258 }
00259
00260 i__3 = *n;
00261 for (k = l; k <= i__3; ++k) {
00262 a[j + k * a_dim1] += s * rv1[k];
00263
00264 }
00265 }
00266
00267 L270:
00268 i__3 = *n;
00269 for (k = l; k <= i__3; ++k) {
00270
00271 a[i__ + k * a_dim1] = scale * a[i__ + k * a_dim1];
00272 }
00273
00274 L290:
00275
00276 d__3 = x, d__4 = (d__1 = w[i__], abs(d__1)) + (d__2 = rv1[i__], abs(
00277 d__2));
00278 x = max(d__3,d__4);
00279
00280 }
00281
00282
00283 i__1 = *n;
00284 for (ii = 1; ii <= i__1; ++ii) {
00285 i__ = *n + 1 - ii;
00286 if (i__ == *n) {
00287 goto L390;
00288 }
00289 if (g == 0.) {
00290 goto L360;
00291 }
00292
00293 i__3 = *n;
00294 for (j = l; j <= i__3; ++j) {
00295
00296
00297
00298 a[j + i__ * a_dim1] = a[i__ + j * a_dim1] / a[i__ + l * a_dim1] /
00299 g;
00300 }
00301
00302 i__3 = *n;
00303 for (j = l; j <= i__3; ++j) {
00304 s = 0.;
00305
00306 i__2 = *n;
00307 for (k = l; k <= i__2; ++k) {
00308
00309 s += a[i__ + k * a_dim1] * a[k + j * a_dim1];
00310 }
00311
00312 i__2 = *n;
00313 for (k = l; k <= i__2; ++k) {
00314 a[k + j * a_dim1] += s * a[k + i__ * a_dim1];
00315
00316 }
00317 }
00318
00319 L360:
00320 i__2 = *n;
00321 for (j = l; j <= i__2; ++j) {
00322 a[i__ + j * a_dim1] = 0.;
00323 a[j + i__ * a_dim1] = 0.;
00324
00325 }
00326
00327 L390:
00328 a[i__ + i__ * a_dim1] = 1.;
00329 g = rv1[i__];
00330 l = i__;
00331
00332 }
00333
00334 if (*m >= *n || *ip == 0) {
00335 goto L510;
00336 }
00337 m1 = *m + 1;
00338
00339 i__1 = *n;
00340 for (i__ = m1; i__ <= i__1; ++i__) {
00341
00342 i__2 = *ip;
00343 for (j = 1; j <= i__2; ++j) {
00344 b[i__ + j * b_dim1] = 0.;
00345
00346 }
00347 }
00348
00349 L510:
00350 tst1 = x;
00351
00352 i__2 = *n;
00353 for (kk = 1; kk <= i__2; ++kk) {
00354 k1 = *n - kk;
00355 k = k1 + 1;
00356 its = 0;
00357
00358
00359 L520:
00360 i__1 = k;
00361 for (ll = 1; ll <= i__1; ++ll) {
00362 l1 = k - ll;
00363 l = l1 + 1;
00364 tst2 = tst1 + (d__1 = rv1[l], abs(d__1));
00365 if (tst2 == tst1) {
00366 goto L565;
00367 }
00368
00369
00370 tst2 = tst1 + (d__1 = w[l1], abs(d__1));
00371 if (tst2 == tst1) {
00372 goto L540;
00373 }
00374
00375 }
00376
00377
00378 L540:
00379 c__ = 0.;
00380 s = 1.;
00381
00382 i__1 = k;
00383 for (i__ = l; i__ <= i__1; ++i__) {
00384 f = s * rv1[i__];
00385 rv1[i__] = c__ * rv1[i__];
00386 tst2 = tst1 + abs(f);
00387 if (tst2 == tst1) {
00388 goto L565;
00389 }
00390 g = w[i__];
00391 h__ = pythag_(&f, &g);
00392 w[i__] = h__;
00393 c__ = g / h__;
00394 s = -f / h__;
00395 if (*ip == 0) {
00396 goto L560;
00397 }
00398
00399 i__3 = *ip;
00400 for (j = 1; j <= i__3; ++j) {
00401 y = b[l1 + j * b_dim1];
00402 z__ = b[i__ + j * b_dim1];
00403 b[l1 + j * b_dim1] = y * c__ + z__ * s;
00404 b[i__ + j * b_dim1] = -y * s + z__ * c__;
00405
00406 }
00407
00408 L560:
00409 ;
00410 }
00411
00412 L565:
00413 z__ = w[k];
00414 if (l == k) {
00415 goto L650;
00416 }
00417
00418 if (its == 30) {
00419 goto L1000;
00420 }
00421 ++its;
00422 x = w[l];
00423 y = w[k1];
00424 g = rv1[k1];
00425 h__ = rv1[k];
00426 f = ((g + z__) / h__ * ((g - z__) / y) + y / h__ - h__ / y) * .5;
00427 g = pythag_(&f, &c_b39);
00428 f = x - z__ / x * z__ + h__ / x * (y / (f + d_sign(&g, &f)) - h__);
00429
00430 c__ = 1.;
00431 s = 1.;
00432
00433 i__1 = k1;
00434 for (i1 = l; i1 <= i__1; ++i1) {
00435 i__ = i1 + 1;
00436 g = rv1[i__];
00437 y = w[i__];
00438 h__ = s * g;
00439 g = c__ * g;
00440 z__ = pythag_(&f, &h__);
00441 rv1[i1] = z__;
00442 c__ = f / z__;
00443 s = h__ / z__;
00444 f = x * c__ + g * s;
00445 g = -x * s + g * c__;
00446 h__ = y * s;
00447 y *= c__;
00448
00449 i__3 = *n;
00450 for (j = 1; j <= i__3; ++j) {
00451 x = a[j + i1 * a_dim1];
00452 z__ = a[j + i__ * a_dim1];
00453 a[j + i1 * a_dim1] = x * c__ + z__ * s;
00454 a[j + i__ * a_dim1] = -x * s + z__ * c__;
00455
00456 }
00457
00458 z__ = pythag_(&f, &h__);
00459 w[i1] = z__;
00460
00461
00462 if (z__ == 0.) {
00463 goto L580;
00464 }
00465 c__ = f / z__;
00466 s = h__ / z__;
00467 L580:
00468 f = c__ * g + s * y;
00469 x = -s * g + c__ * y;
00470 if (*ip == 0) {
00471 goto L600;
00472 }
00473
00474 i__3 = *ip;
00475 for (j = 1; j <= i__3; ++j) {
00476 y = b[i1 + j * b_dim1];
00477 z__ = b[i__ + j * b_dim1];
00478 b[i1 + j * b_dim1] = y * c__ + z__ * s;
00479 b[i__ + j * b_dim1] = -y * s + z__ * c__;
00480
00481 }
00482
00483 L600:
00484 ;
00485 }
00486
00487 rv1[l] = 0.;
00488 rv1[k] = f;
00489 w[k] = x;
00490 goto L520;
00491
00492 L650:
00493 if (z__ >= 0.) {
00494 goto L700;
00495 }
00496
00497 w[k] = -z__;
00498
00499 i__1 = *n;
00500 for (j = 1; j <= i__1; ++j) {
00501
00502 a[j + k * a_dim1] = -a[j + k * a_dim1];
00503 }
00504
00505 L700:
00506 ;
00507 }
00508
00509 goto L1001;
00510
00511
00512 L1000:
00513 *ierr = k;
00514 L1001:
00515 return 0;
00516 }
00517