00001 /* minfit.f -- translated by f2c (version 19961017).
00002    You must link the resulting object file with the libraries:
00003         -lf2c -lm   (in that order)
00004 */
00005 
00006 #include "f2c.h"
00007 
00008 /* Table of constant values */
00009 
00010 static doublereal c_b39 = 1.;
00011 
00012 /* Subroutine */ int minfit_(integer *nm, integer *m, integer *n, doublereal *
00013         a, doublereal *w, integer *ip, doublereal *b, integer *ierr, 
00014         doublereal *rv1)
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_ */
00517