00001 /* svd.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_b47 = 1.;
00011 
00012 /* Subroutine */ int svd_(integer *m, integer *n, integer *lda, doublereal *a,
00013          doublereal *w, logical *matu, integer *ldu, doublereal *u, logical *
00014         matv, integer *ldv, doublereal *v, integer *ierr, doublereal *rv1)
00015 {
00016     /* System generated locals */
00017     integer a_dim1, a_offset, u_dim1, u_offset, v_dim1, v_offset, i__1, i__2, 
00018             i__3;
00019     doublereal d__1, d__2, d__3, d__4;
00020 
00021     /* Builtin functions */
00022     double sqrt(doublereal), d_sign(doublereal *, doublereal *);
00023 
00024     /* Local variables */
00025     static doublereal c__, f, g, h__;
00026     static integer i__, j, k, l;
00027     static doublereal s, x, y, z__, scale;
00028     static integer i1, k1, l1, ii, kk, ll, mn;
00029     extern doublereal pythag_(doublereal *, doublereal *);
00030     static integer its;
00031     static doublereal tst1, tst2;
00032 
00033 
00034 
00035 /*     this subroutine is a translation of the algol procedure svd, */
00036 /*     num. math. 14, 403-420(1970) by golub and reinsch. */
00037 /*     handbook for auto. comp., vol ii-linear algebra, 134-151(1971). */
00038 
00039 /*     this subroutine determines the singular value decomposition */
00040 /*          t */
00041 /*     a=usv  of a real m by n rectangular matrix.  householder */
00042 /*     bidiagonalization and a variant of the qr algorithm are used. */
00043 
00044 /*     on input */
00045 
00046 /*        nm must be set to the row dimension of two-dimensional */
00047 /*          array parameters as declared in the calling program */
00048 /*          dimension statement.  note that nm must be at least */
00049 /*          as large as the maximum of m and n. */
00050 
00051 /*        m is the number of rows of a (and u). */
00052 
00053 /*        n is the number of columns of a, u, and v */
00054 
00055 /*        a contains the rectangular input matrix to be decomposed. */
00056 
00057 /*        matu should be set to .true. if the u matrix in the */
00058 /*          decomposition is desired, and to .false. otherwise. */
00059 
00060 /*        matv should be set to .true. if the v matrix in the */
00061 /*          decomposition is desired, and to .false. otherwise. */
00062 
00063 /*        lda, ldu, ldv: are the leading dimensions of matrices */
00064 /*          a, u, and v (respectively);  must have */
00065 /*           lda >= m ; ldu >= m ; ldv >= n */
00066 
00067 /*     on output */
00068 
00069 /*        a is unaltered (unless overwritten by u or v). */
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 /*        u contains the matrix u (orthogonal column vectors) of the */
00077 /*          decomposition if matu has been set to .true.  otherwise */
00078 /*          u is used as a temporary array.  u may coincide with a. */
00079 /*          if an error exit is made, the columns of u corresponding */
00080 /*          to indices of correct singular values should be correct. */
00081 
00082 /*        v contains the matrix v (orthogonal) of the decomposition if */
00083 /*          matv has been set to .true.  otherwise v is not referenced. */
00084 /*          v may also coincide with a if u is not needed.  if an error */
00085 /*          exit is made, the columns of v corresponding to indices of */
00086 /*          correct singular values should be correct. */
00087 
00088 /*        ierr is set to */
00089 /*          zero       for normal return, */
00090 /*          k          if the k-th singular value has not been */
00091 /*                     determined after 30 iterations. */
00092 
00093 /*        rv1 is a temporary storage array. */
00094 
00095 /*     calls pythag for  dsqrt(a*a + b*b) . */
00096 
00097 /*     questions and comments should be directed to burton s. garbow, */
00098 /*     mathematics and computer science div, argonne national laboratory 
00099 */
00100 
00101 /*     this version dated august 1983. */
00102 
00103 /*     ------------------------------------------------------------------ 
00104 */
00105 
00106     /* Parameter adjustments */
00107     --rv1;
00108     --w;
00109     a_dim1 = *lda;
00110     a_offset = a_dim1 + 1;
00111     a -= a_offset;
00112     u_dim1 = *ldu;
00113     u_offset = u_dim1 + 1;
00114     u -= u_offset;
00115     v_dim1 = *ldv;
00116     v_offset = v_dim1 + 1;
00117     if( v != (doublereal *)0 ) v -= v_offset;
00118 
00119     /* Function Body */
00120     *ierr = 0;
00121 
00122     i__1 = *m;
00123     for (i__ = 1; i__ <= i__1; ++i__) {
00124 
00125         i__2 = *n;
00126         for (j = 1; j <= i__2; ++j) {
00127             u[i__ + j * u_dim1] = a[i__ + j * a_dim1];
00128 /* L100: */
00129         }
00130     }
00131 /*     .......... householder reduction to bidiagonal form .......... */
00132     g = 0.;
00133     scale = 0.;
00134     x = 0.;
00135 
00136     i__2 = *n;
00137     for (i__ = 1; i__ <= i__2; ++i__) {
00138         l = i__ + 1;
00139         rv1[i__] = scale * g;
00140         g = 0.;
00141         s = 0.;
00142         scale = 0.;
00143         if (i__ > *m) {
00144             goto L210;
00145         }
00146 
00147         i__1 = *m;
00148         for (k = i__; k <= i__1; ++k) {
00149 /* L120: */
00150             scale += (d__1 = u[k + i__ * u_dim1], abs(d__1));
00151         }
00152 
00153         if (scale == 0.) {
00154             goto L210;
00155         }
00156 
00157         i__1 = *m;
00158         for (k = i__; k <= i__1; ++k) {
00159             u[k + i__ * u_dim1] /= scale;
00160 /* Computing 2nd power */
00161             d__1 = u[k + i__ * u_dim1];
00162             s += d__1 * d__1;
00163 /* L130: */
00164         }
00165 
00166         f = u[i__ + i__ * u_dim1];
00167         d__1 = sqrt(s);
00168         g = -d_sign(&d__1, &f);
00169         h__ = f * g - s;
00170         u[i__ + i__ * u_dim1] = f - g;
00171         if (i__ == *n) {
00172             goto L190;
00173         }
00174 
00175         i__1 = *n;
00176         for (j = l; j <= i__1; ++j) {
00177             s = 0.;
00178 
00179             i__3 = *m;
00180             for (k = i__; k <= i__3; ++k) {
00181 /* L140: */
00182                 s += u[k + i__ * u_dim1] * u[k + j * u_dim1];
00183             }
00184 
00185             f = s / h__;
00186 
00187             i__3 = *m;
00188             for (k = i__; k <= i__3; ++k) {
00189                 u[k + j * u_dim1] += f * u[k + i__ * u_dim1];
00190 /* L150: */
00191             }
00192         }
00193 
00194 L190:
00195         i__3 = *m;
00196         for (k = i__; k <= i__3; ++k) {
00197 /* L200: */
00198             u[k + i__ * u_dim1] = scale * u[k + i__ * u_dim1];
00199         }
00200 
00201 L210:
00202         w[i__] = scale * g;
00203         g = 0.;
00204         s = 0.;
00205         scale = 0.;
00206         if (i__ > *m || i__ == *n) {
00207             goto L290;
00208         }
00209 
00210         i__3 = *n;
00211         for (k = l; k <= i__3; ++k) {
00212 /* L220: */
00213             scale += (d__1 = u[i__ + k * u_dim1], abs(d__1));
00214         }
00215 
00216         if (scale == 0.) {
00217             goto L290;
00218         }
00219 
00220         i__3 = *n;
00221         for (k = l; k <= i__3; ++k) {
00222             u[i__ + k * u_dim1] /= scale;
00223 /* Computing 2nd power */
00224             d__1 = u[i__ + k * u_dim1];
00225             s += d__1 * d__1;
00226 /* L230: */
00227         }
00228 
00229         f = u[i__ + l * u_dim1];
00230         d__1 = sqrt(s);
00231         g = -d_sign(&d__1, &f);
00232         h__ = f * g - s;
00233         u[i__ + l * u_dim1] = f - g;
00234 
00235         i__3 = *n;
00236         for (k = l; k <= i__3; ++k) {
00237 /* L240: */
00238             rv1[k] = u[i__ + k * u_dim1] / h__;
00239         }
00240 
00241         if (i__ == *m) {
00242             goto L270;
00243         }
00244 
00245         i__3 = *m;
00246         for (j = l; j <= i__3; ++j) {
00247             s = 0.;
00248 
00249             i__1 = *n;
00250             for (k = l; k <= i__1; ++k) {
00251 /* L250: */
00252                 s += u[j + k * u_dim1] * u[i__ + k * u_dim1];
00253             }
00254 
00255             i__1 = *n;
00256             for (k = l; k <= i__1; ++k) {
00257                 u[j + k * u_dim1] += s * rv1[k];
00258 /* L260: */
00259             }
00260         }
00261 
00262 L270:
00263         i__1 = *n;
00264         for (k = l; k <= i__1; ++k) {
00265 /* L280: */
00266             u[i__ + k * u_dim1] = scale * u[i__ + k * u_dim1];
00267         }
00268 
00269 L290:
00270 /* Computing MAX */
00271         d__3 = x, d__4 = (d__1 = w[i__], abs(d__1)) + (d__2 = rv1[i__], abs(
00272                 d__2));
00273         x = max(d__3,d__4);
00274 /* L300: */
00275     }
00276 /*     .......... accumulation of right-hand transformations .......... */
00277     if (! (*matv)) {
00278         goto L410;
00279     }
00280 /*     .......... for i=n step -1 until 1 do -- .......... */
00281     i__2 = *n;
00282     for (ii = 1; ii <= i__2; ++ii) {
00283         i__ = *n + 1 - ii;
00284         if (i__ == *n) {
00285             goto L390;
00286         }
00287         if (g == 0.) {
00288             goto L360;
00289         }
00290 
00291         i__1 = *n;
00292         for (j = l; j <= i__1; ++j) {
00293 /*     .......... double division avoids possible underflow ......
00294 .... */
00295 /* L320: */
00296             v[j + i__ * v_dim1] = u[i__ + j * u_dim1] / u[i__ + l * u_dim1] / 
00297                     g;
00298         }
00299 
00300         i__1 = *n;
00301         for (j = l; j <= i__1; ++j) {
00302             s = 0.;
00303 
00304             i__3 = *n;
00305             for (k = l; k <= i__3; ++k) {
00306 /* L340: */
00307                 s += u[i__ + k * u_dim1] * v[k + j * v_dim1];
00308             }
00309 
00310             i__3 = *n;
00311             for (k = l; k <= i__3; ++k) {
00312                 v[k + j * v_dim1] += s * v[k + i__ * v_dim1];
00313 /* L350: */
00314             }
00315         }
00316 
00317 L360:
00318         i__3 = *n;
00319         for (j = l; j <= i__3; ++j) {
00320             v[i__ + j * v_dim1] = 0.;
00321             v[j + i__ * v_dim1] = 0.;
00322 /* L380: */
00323         }
00324 
00325 L390:
00326         v[i__ + i__ * v_dim1] = 1.;
00327         g = rv1[i__];
00328         l = i__;
00329 /* L400: */
00330     }
00331 /*     .......... accumulation of left-hand transformations .......... */
00332 L410:
00333     if (! (*matu)) {
00334         goto L510;
00335     }
00336 /*     ..........for i=min(m,n) step -1 until 1 do -- .......... */
00337     mn = *n;
00338     if (*m < *n) {
00339         mn = *m;
00340     }
00341 
00342     i__2 = mn;
00343     for (ii = 1; ii <= i__2; ++ii) {
00344         i__ = mn + 1 - ii;
00345         l = i__ + 1;
00346         g = w[i__];
00347         if (i__ == *n) {
00348             goto L430;
00349         }
00350 
00351         i__3 = *n;
00352         for (j = l; j <= i__3; ++j) {
00353 /* L420: */
00354             u[i__ + j * u_dim1] = 0.;
00355         }
00356 
00357 L430:
00358         if (g == 0.) {
00359             goto L475;
00360         }
00361         if (i__ == mn) {
00362             goto L460;
00363         }
00364 
00365         i__3 = *n;
00366         for (j = l; j <= i__3; ++j) {
00367             s = 0.;
00368 
00369             i__1 = *m;
00370             for (k = l; k <= i__1; ++k) {
00371 /* L440: */
00372                 s += u[k + i__ * u_dim1] * u[k + j * u_dim1];
00373             }
00374 /*     .......... double division avoids possible underflow ......
00375 .... */
00376             f = s / u[i__ + i__ * u_dim1] / g;
00377 
00378             i__1 = *m;
00379             for (k = i__; k <= i__1; ++k) {
00380                 u[k + j * u_dim1] += f * u[k + i__ * u_dim1];
00381 /* L450: */
00382             }
00383         }
00384 
00385 L460:
00386         i__1 = *m;
00387         for (j = i__; j <= i__1; ++j) {
00388 /* L470: */
00389             u[j + i__ * u_dim1] /= g;
00390         }
00391 
00392         goto L490;
00393 
00394 L475:
00395         i__1 = *m;
00396         for (j = i__; j <= i__1; ++j) {
00397 /* L480: */
00398             u[j + i__ * u_dim1] = 0.;
00399         }
00400 
00401 L490:
00402         u[i__ + i__ * u_dim1] += 1.;
00403 /* L500: */
00404     }
00405 /*     .......... diagonalization of the bidiagonal form .......... */
00406 L510:
00407     tst1 = x;
00408 /*     .......... for k=n step -1 until 1 do -- .......... */
00409     i__2 = *n;
00410     for (kk = 1; kk <= i__2; ++kk) {
00411         k1 = *n - kk;
00412         k = k1 + 1;
00413         its = 0;
00414 /*     .......... test for splitting. */
00415 /*                for l=k step -1 until 1 do -- .......... */
00416 L520:
00417         i__1 = k;
00418         for (ll = 1; ll <= i__1; ++ll) {
00419             l1 = k - ll;
00420             l = l1 + 1;
00421             tst2 = tst1 + (d__1 = rv1[l], abs(d__1));
00422             if (tst2 == tst1) {
00423                 goto L565;
00424             }
00425 /*     .......... rv1(1) is always zero, so there is no exit */
00426 /*                through the bottom of the loop .......... */
00427             tst2 = tst1 + (d__1 = w[l1], abs(d__1));
00428             if (tst2 == tst1) {
00429                 goto L540;
00430             }
00431 /* L530: */
00432         }
00433 /*     .......... cancellation of rv1(l) if l greater than 1 .........
00434 . */
00435 L540:
00436         c__ = 0.;
00437         s = 1.;
00438 
00439         i__1 = k;
00440         for (i__ = l; i__ <= i__1; ++i__) {
00441             f = s * rv1[i__];
00442             rv1[i__] = c__ * rv1[i__];
00443             tst2 = tst1 + abs(f);
00444             if (tst2 == tst1) {
00445                 goto L565;
00446             }
00447             g = w[i__];
00448             h__ = pythag_(&f, &g);
00449             w[i__] = h__;
00450             c__ = g / h__;
00451             s = -f / h__;
00452             if (! (*matu)) {
00453                 goto L560;
00454             }
00455 
00456             i__3 = *m;
00457             for (j = 1; j <= i__3; ++j) {
00458                 y = u[j + l1 * u_dim1];
00459                 z__ = u[j + i__ * u_dim1];
00460                 u[j + l1 * u_dim1] = y * c__ + z__ * s;
00461                 u[j + i__ * u_dim1] = -y * s + z__ * c__;
00462 /* L550: */
00463             }
00464 
00465 L560:
00466             ;
00467         }
00468 /*     .......... test for convergence .......... */
00469 L565:
00470         z__ = w[k];
00471         if (l == k) {
00472             goto L650;
00473         }
00474 /*     .......... shift from bottom 2 by 2 minor .......... */
00475         if (its == 30) {
00476             goto L1000;
00477         }
00478         ++its;
00479         x = w[l];
00480         y = w[k1];
00481         g = rv1[k1];
00482         h__ = rv1[k];
00483         f = ((g + z__) / h__ * ((g - z__) / y) + y / h__ - h__ / y) * .5;
00484         g = pythag_(&f, &c_b47);
00485         f = x - z__ / x * z__ + h__ / x * (y / (f + d_sign(&g, &f)) - h__);
00486 /*     .......... next qr transformation .......... */
00487         c__ = 1.;
00488         s = 1.;
00489 
00490         i__1 = k1;
00491         for (i1 = l; i1 <= i__1; ++i1) {
00492             i__ = i1 + 1;
00493             g = rv1[i__];
00494             y = w[i__];
00495             h__ = s * g;
00496             g = c__ * g;
00497             z__ = pythag_(&f, &h__);
00498             rv1[i1] = z__;
00499             c__ = f / z__;
00500             s = h__ / z__;
00501             f = x * c__ + g * s;
00502             g = -x * s + g * c__;
00503             h__ = y * s;
00504             y *= c__;
00505             if (! (*matv)) {
00506                 goto L575;
00507             }
00508 
00509             i__3 = *n;
00510             for (j = 1; j <= i__3; ++j) {
00511                 x = v[j + i1 * v_dim1];
00512                 z__ = v[j + i__ * v_dim1];
00513                 v[j + i1 * v_dim1] = x * c__ + z__ * s;
00514                 v[j + i__ * v_dim1] = -x * s + z__ * c__;
00515 /* L570: */
00516             }
00517 
00518 L575:
00519             z__ = pythag_(&f, &h__);
00520             w[i1] = z__;
00521 /*     .......... rotation can be arbitrary if z is zero .........
00522 . */
00523             if (z__ == 0.) {
00524                 goto L580;
00525             }
00526             c__ = f / z__;
00527             s = h__ / z__;
00528 L580:
00529             f = c__ * g + s * y;
00530             x = -s * g + c__ * y;
00531             if (! (*matu)) {
00532                 goto L600;
00533             }
00534 
00535             i__3 = *m;
00536             for (j = 1; j <= i__3; ++j) {
00537                 y = u[j + i1 * u_dim1];
00538                 z__ = u[j + i__ * u_dim1];
00539                 u[j + i1 * u_dim1] = y * c__ + z__ * s;
00540                 u[j + i__ * u_dim1] = -y * s + z__ * c__;
00541 /* L590: */
00542             }
00543 
00544 L600:
00545             ;
00546         }
00547 
00548         rv1[l] = 0.;
00549         rv1[k] = f;
00550         w[k] = x;
00551         goto L520;
00552 /*     .......... convergence .......... */
00553 L650:
00554         if (z__ >= 0.) {
00555             goto L700;
00556         }
00557 /*     .......... w(k) is made non-negative .......... */
00558         w[k] = -z__;
00559         if (! (*matv)) {
00560             goto L700;
00561         }
00562 
00563         i__1 = *n;
00564         for (j = 1; j <= i__1; ++j) {
00565 /* L690: */
00566             v[j + k * v_dim1] = -v[j + k * v_dim1];
00567         }
00568 
00569 L700:
00570         ;
00571     }
00572 
00573     goto L1001;
00574 /*     .......... set error -- no convergence to a */
00575 /*                singular value after 30 iterations .......... */
00576 L1000:
00577     *ierr = k;
00578 L1001:
00579     return 0;
00580 } /* svd_ */
00581