Doxygen Source Code Documentation
eis_minfit.c File Reference
#include "f2c.h"Go to the source code of this file.
Functions | |
| int | minfit_ (integer *nm, integer *m, integer *n, doublereal *a, doublereal *w, integer *ip, doublereal *b, integer *ierr, doublereal *rv1) |
Variables | |
| doublereal | c_b39 = 1. |
Function Documentation
|
||||||||||||||||||||||||||||||||||||||||
|
Definition at line 12 of file eis_minfit.c. References a, abs, c_b39, d_sign(), i1, l, m1, max, pythag_(), and scale.
00015 {
00016 /* System generated locals */
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 /* Builtin functions */
00021 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
00022
00023 /* Local variables */
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 /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE MINFIT, */
00035 /* NUM. MATH. 14, 403-420(1970) BY GOLUB AND REINSCH. */
00036 /* HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971). */
00037
00038 /* THIS SUBROUTINE DETERMINES, TOWARDS THE SOLUTION OF THE LINEAR */
00039 /* T */
00040 /* SYSTEM AX=B, THE SINGULAR VALUE DECOMPOSITION A=USV OF A REAL */
00041 /* T */
00042 /* M BY N RECTANGULAR MATRIX, FORMING U B RATHER THAN U. HOUSEHOLDER
00043 */
00044 /* BIDIAGONALIZATION AND A VARIANT OF THE QR ALGORITHM ARE USED. */
00045
00046 /* ON INPUT */
00047
00048 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00049 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00050 /* DIMENSION STATEMENT. NOTE THAT NM MUST BE AT LEAST */
00051 /* AS LARGE AS THE MAXIMUM OF M AND N. */
00052
00053 /* M IS THE NUMBER OF ROWS OF A AND B. */
00054
00055 /* N IS THE NUMBER OF COLUMNS OF A AND THE ORDER OF V. */
00056
00057 /* A CONTAINS THE RECTANGULAR COEFFICIENT MATRIX OF THE SYSTEM. */
00058
00059 /* IP IS THE NUMBER OF COLUMNS OF B. IP CAN BE ZERO. */
00060
00061 /* B CONTAINS THE CONSTANT COLUMN MATRIX OF THE SYSTEM */
00062 /* IF IP IS NOT ZERO. OTHERWISE B IS NOT REFERENCED. */
00063
00064 /* ON OUTPUT */
00065
00066 /* A HAS BEEN OVERWRITTEN BY THE MATRIX V (ORTHOGONAL) OF THE */
00067 /* DECOMPOSITION IN ITS FIRST N ROWS AND COLUMNS. IF AN */
00068 /* ERROR EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO */
00069 /* INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT. */
00070
00071 /* W CONTAINS THE N (NON-NEGATIVE) SINGULAR VALUES OF A (THE */
00072 /* DIAGONAL ELEMENTS OF S). THEY ARE UNORDERED. IF AN */
00073 /* ERROR EXIT IS MADE, THE SINGULAR VALUES SHOULD BE CORRECT */
00074 /* FOR INDICES IERR+1,IERR+2,...,N. */
00075
00076 /* T */
00077 /* B HAS BEEN OVERWRITTEN BY U B. IF AN ERROR EXIT IS MADE, */
00078 /* T */
00079 /* THE ROWS OF U B CORRESPONDING TO INDICES OF CORRECT */
00080 /* SINGULAR VALUES SHOULD BE CORRECT. */
00081
00082 /* IERR IS SET TO */
00083 /* ZERO FOR NORMAL RETURN, */
00084 /* K IF THE K-TH SINGULAR VALUE HAS NOT BEEN */
00085 /* DETERMINED AFTER 30 ITERATIONS. */
00086
00087 /* RV1 IS A TEMPORARY STORAGE ARRAY. */
00088
00089 /* CALLS PYTHAG FOR DSQRT(A*A + B*B) . */
00090
00091 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00092 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
00093 */
00094
00095 /* THIS VERSION DATED AUGUST 1983. */
00096
00097 /* ------------------------------------------------------------------
00098 */
00099
00100 /* Parameter adjustments */
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 /* Function Body */
00111 *ierr = 0;
00112 /* .......... HOUSEHOLDER REDUCTION TO BIDIAGONAL FORM .......... */
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 /* L120: */
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 /* Computing 2nd power */
00142 d__1 = a[k + i__ * a_dim1];
00143 s += d__1 * d__1;
00144 /* L130: */
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 /* L140: */
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 /* L150: */
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 /* L170: */
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 /* L180: */
00196 }
00197 }
00198
00199 L190:
00200 i__2 = *m;
00201 for (k = i__; k <= i__2; ++k) {
00202 /* L200: */
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 /* L220: */
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 /* Computing 2nd power */
00229 d__1 = a[i__ + k * a_dim1];
00230 s += d__1 * d__1;
00231 /* L230: */
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 /* L240: */
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 /* L250: */
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 /* L260: */
00264 }
00265 }
00266
00267 L270:
00268 i__3 = *n;
00269 for (k = l; k <= i__3; ++k) {
00270 /* L280: */
00271 a[i__ + k * a_dim1] = scale * a[i__ + k * a_dim1];
00272 }
00273
00274 L290:
00275 /* Computing MAX */
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 /* L300: */
00280 }
00281 /* .......... ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS. */
00282 /* FOR I=N STEP -1 UNTIL 1 DO -- .......... */
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 /* .......... DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW ......
00296 .... */
00297 /* L320: */
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 /* L340: */
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 /* L350: */
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 /* L380: */
00325 }
00326
00327 L390:
00328 a[i__ + i__ * a_dim1] = 1.;
00329 g = rv1[i__];
00330 l = i__;
00331 /* L400: */
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 /* L500: */
00346 }
00347 }
00348 /* .......... DIAGONALIZATION OF THE BIDIAGONAL FORM .......... */
00349 L510:
00350 tst1 = x;
00351 /* .......... FOR K=N STEP -1 UNTIL 1 DO -- .......... */
00352 i__2 = *n;
00353 for (kk = 1; kk <= i__2; ++kk) {
00354 k1 = *n - kk;
00355 k = k1 + 1;
00356 its = 0;
00357 /* .......... TEST FOR SPLITTING. */
00358 /* FOR L=K STEP -1 UNTIL 1 DO -- .......... */
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 /* .......... RV1(1) IS ALWAYS ZERO, SO THERE IS NO EXIT */
00369 /* THROUGH THE BOTTOM OF THE LOOP .......... */
00370 tst2 = tst1 + (d__1 = w[l1], abs(d__1));
00371 if (tst2 == tst1) {
00372 goto L540;
00373 }
00374 /* L530: */
00375 }
00376 /* .......... CANCELLATION OF RV1(L) IF L GREATER THAN 1 .........
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 /* L550: */
00406 }
00407
00408 L560:
00409 ;
00410 }
00411 /* .......... TEST FOR CONVERGENCE .......... */
00412 L565:
00413 z__ = w[k];
00414 if (l == k) {
00415 goto L650;
00416 }
00417 /* .......... SHIFT FROM BOTTOM 2 BY 2 MINOR .......... */
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 /* .......... NEXT QR TRANSFORMATION .......... */
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 /* L570: */
00456 }
00457
00458 z__ = pythag_(&f, &h__);
00459 w[i1] = z__;
00460 /* .......... ROTATION CAN BE ARBITRARY IF Z IS ZERO .........
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 /* L590: */
00481 }
00482
00483 L600:
00484 ;
00485 }
00486
00487 rv1[l] = 0.;
00488 rv1[k] = f;
00489 w[k] = x;
00490 goto L520;
00491 /* .......... CONVERGENCE .......... */
00492 L650:
00493 if (z__ >= 0.) {
00494 goto L700;
00495 }
00496 /* .......... W(K) IS MADE NON-NEGATIVE .......... */
00497 w[k] = -z__;
00498
00499 i__1 = *n;
00500 for (j = 1; j <= i__1; ++j) {
00501 /* L690: */
00502 a[j + k * a_dim1] = -a[j + k * a_dim1];
00503 }
00504
00505 L700:
00506 ;
00507 }
00508
00509 goto L1001;
00510 /* .......... SET ERROR -- NO CONVERGENCE TO A */
00511 /* SINGULAR VALUE AFTER 30 ITERATIONS .......... */
00512 L1000:
00513 *ierr = k;
00514 L1001:
00515 return 0;
00516 } /* minfit_ */
|
Variable Documentation
|
|
Definition at line 10 of file eis_minfit.c. Referenced by minfit_(). |
