nnls.c

Go to the documentation of this file.
00001 /* Distributed with ASYNCHRONOUS PARALLEL PATTERN SEARCH (APPS) */
00002 
00003 /* The routines in this file have been translated from Fortran to C by
00004    f2c. Additional modifications have been made to remove the
00005    dependencies on the f2c header file and library. The original
00006    Fortran 77 code accompanies the SIAM Publications printing of
00007    "Solving Least Squares Problems," by C. Lawson and R. Hanson and is
00008    freely available at www.netlib.org/lawson-hanson/all. */
00009 
00010 /* nnls.F -- translated by f2c (version 19970805).
00011    You must link the resulting object file with the libraries:
00012     -lf2c -lm   (in that order)
00013 */
00014 
00015 /* The next line was removed after the f2c translation */
00016 /* #include "f2c.h" */
00017 
00018 /* The next lines were added after the f2c translation. Also swapped
00019    abs for nnls_abs and max for nnls_max to avoid confusion with some
00020    compilers. */
00021 #include <stdio.h>
00022 #include <math.h>
00023 #define nnls_max(a,b) ((a) >= (b) ? (a) : (b))
00024 #define nnls_abs(x) ((x) >= 0 ? (x) : -(x))
00025 typedef int integer;
00026 typedef double doublereal;
00027 
00028 /* The following subroutine was added after the f2c translation */
00029 double d_sign(double *a, double *b)
00030 {
00031   double x;
00032   x = (*a >= 0 ? *a : - *a);
00033   return( *b >= 0 ? x : -x);
00034 }
00035 
00036 /* Table of constant values */
00037 
00038 static integer c__1 = 1;
00039 static integer c__0 = 0;
00040 static integer c__2 = 2;
00041 
00042 /*     SUBROUTINE NNLS  (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) */
00043 
00044 /*  Algorithm NNLS: NONNEGATIVE LEAST SQUARES */
00045 
00046 /*  The original version of this code was developed by */
00047 /*  Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
00048 /*  1973 JUN 15, and published in the book */
00049 /*  "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
00050 /*  Revised FEB 1995 to accompany reprinting of the book by SIAM. */
00051 
00052 /*     GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B,  COMPUTE AN */
00053 /*     N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM */
00054 
00055 /*                      A * X = B  SUBJECT TO X .GE. 0 */
00056 /*     ------------------------------------------------------------------ */
00057 /*                     Subroutine Arguments */
00058 
00059 /*     A(),MDA,M,N     MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE */
00060 /*                     ARRAY, A().   ON ENTRY A() CONTAINS THE M BY N */
00061 /*                     MATRIX, A.           ON EXIT A() CONTAINS */
00062 /*                     THE PRODUCT MATRIX, Q*A , WHERE Q IS AN */
00063 /*                     M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY */
00064 /*                     THIS SUBROUTINE. */
00065 /*     B()     ON ENTRY B() CONTAINS THE M-VECTOR, B.   ON EXIT B() CON- */
00066 /*             TAINS Q*B. */
00067 /*     X()     ON ENTRY X() NEED NOT BE INITIALIZED.  ON EXIT X() WILL */
00068 /*             CONTAIN THE SOLUTION VECTOR. */
00069 /*     RNORM   ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE */
00070 /*             RESIDUAL VECTOR. */
00071 /*     W()     AN N-ARRAY OF WORKING SPACE.  ON EXIT W() WILL CONTAIN */
00072 /*             THE DUAL SOLUTION VECTOR.   W WILL SATISFY W(I) = 0. */
00073 /*             FOR ALL I IN SET P  AND W(I) .LE. 0. FOR ALL I IN SET Z */
00074 /*     ZZ()     AN M-ARRAY OF WORKING SPACE. */
00075 /*     INDEX()     AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N. */
00076 /*                 ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS */
00077 /*                 P AND Z AS FOLLOWS.. */
00078 
00079 /*                 INDEX(1)   THRU INDEX(NSETP) = SET P. */
00080 /*                 INDEX(IZ1) THRU INDEX(IZ2)   = SET Z. */
00081 /*                 IZ1 = NSETP + 1 = NPP1 */
00082 /*                 IZ2 = N */
00083 /*     MODE    THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING */
00084 /*             MEANINGS. */
00085 /*             1     THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY. */
00086 /*             2     THE DIMENSIONS OF THE PROBLEM ARE BAD. */
00087 /*                   EITHER M .LE. 0 OR N .LE. 0. */
00088 /*             3    ITERATION COUNT EXCEEDED.  MORE THAN 3*N ITERATIONS. */
00089 
00090 /*     ------------------------------------------------------------------ */
00091 /* Subroutine */ int nnls_(a, mda, m, n, b, x, rnorm, w, zz, index, mode)
00092 doublereal *a;
00093 integer *mda, *m, *n;
00094 doublereal *b, *x, *rnorm, *w, *zz;
00095 integer *index, *mode;
00096 {
00097     /* System generated locals */
00098     integer a_dim1, a_offset, i__1, i__2;
00099     doublereal d__1, d__2;
00100 
00101     /* Builtin functions */
00102     /* The following lines were commented out after the f2c translation */
00103     /* double sqrt(); */
00104     /* integer s_wsfe(), do_fio(), e_wsfe(); */
00105 
00106     /* Local variables */
00107     extern doublereal diff_();
00108     static integer iter;
00109     static doublereal temp, wmax;
00110     static integer i__, j, l;
00111     static doublereal t, alpha, asave;
00112     static integer itmax, izmax, nsetp;
00113     extern /* Subroutine */ int g1_();
00114     static doublereal dummy, unorm, ztest, cc;
00115     extern /* Subroutine */ int h12_();
00116     static integer ii, jj, ip;
00117     static doublereal sm;
00118     static integer iz, jz;
00119     static doublereal up, ss;
00120     static integer rtnkey, iz1, iz2, npp1;
00121 
00122     /* Fortran I/O blocks */
00123     /* The following line was commented out after the f2c translation */
00124     /* static cilist io___22 = { 0, 6, 0, "(/a)", 0 }; */
00125 
00126 
00127 /*     ------------------------------------------------------------------ 
00128 */
00129 /*     integer INDEX(N) */
00130 /*     double precision A(MDA,N), B(M), W(N), X(N), ZZ(M) */
00131 /*     ------------------------------------------------------------------ 
00132 */
00133     /* Parameter adjustments */
00134     a_dim1 = *mda;
00135     a_offset = a_dim1 + 1;
00136     a -= a_offset;
00137     --b;
00138     --x;
00139     --w;
00140     --zz;
00141     --index;
00142 
00143     /* Function Body */
00144     *mode = 1;
00145     if (*m <= 0 || *n <= 0) {
00146     *mode = 2;
00147     return 0;
00148     }
00149     iter = 0;
00150     itmax = *n * 3;
00151 
00152 /*                    INITIALIZE THE ARRAYS INDEX() AND X(). */
00153 
00154     i__1 = *n;
00155     for (i__ = 1; i__ <= i__1; ++i__) {
00156     x[i__] = 0.;
00157 /* L20: */
00158     index[i__] = i__;
00159     }
00160 
00161     iz2 = *n;
00162     iz1 = 1;
00163     nsetp = 0;
00164     npp1 = 1;
00165 /*                             ******  MAIN LOOP BEGINS HERE  ****** */
00166 L30:
00167 /*                  QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION. 
00168 */
00169 /*                        OR IF M COLS OF A HAVE BEEN TRIANGULARIZED. */
00170 
00171     if (iz1 > iz2 || nsetp >= *m) {
00172     goto L350;
00173     }
00174 
00175 /*         COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W(). 
00176 */
00177 
00178     i__1 = iz2;
00179     for (iz = iz1; iz <= i__1; ++iz) {
00180     j = index[iz];
00181     sm = 0.;
00182     i__2 = *m;
00183     for (l = npp1; l <= i__2; ++l) {
00184 /* L40: */
00185         sm += a[l + j * a_dim1] * b[l];
00186     }
00187     w[j] = sm;
00188 /* L50: */
00189     }
00190 /*                                   FIND LARGEST POSITIVE W(J). */
00191 L60:
00192     wmax = 0.;
00193     i__1 = iz2;
00194     for (iz = iz1; iz <= i__1; ++iz) {
00195     j = index[iz];
00196     if (w[j] > wmax) {
00197         wmax = w[j];
00198         izmax = iz;
00199     }
00200 /* L70: */
00201     }
00202 
00203 /*             IF WMAX .LE. 0. GO TO TERMINATION. */
00204 /*             THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS. 
00205 */
00206 
00207     if (wmax <= 0.) {
00208     goto L350;
00209     }
00210     iz = izmax;
00211     j = index[iz];
00212 
00213 /*     THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P. */
00214 /*     BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID */
00215 /*     NEAR LINEAR DEPENDENCE. */
00216 
00217     asave = a[npp1 + j * a_dim1];
00218     i__1 = npp1 + 1;
00219     h12_(&c__1, &npp1, &i__1, m, &a[j * a_dim1 + 1], &c__1, &up, &dummy, &
00220         c__1, &c__1, &c__0);
00221     unorm = 0.;
00222     if (nsetp != 0) {
00223     i__1 = nsetp;
00224     for (l = 1; l <= i__1; ++l) {
00225 /* L90: */
00226 /* Computing 2nd power */
00227         d__1 = a[l + j * a_dim1];
00228         unorm += d__1 * d__1;
00229     }
00230     }
00231     unorm = sqrt(unorm);
00232     d__2 = unorm + (d__1 = a[npp1 + j * a_dim1], nnls_abs(d__1)) * .01;
00233     if (diff_(&d__2, &unorm) > 0.) {
00234 
00235 /*        COL J IS SUFFICIENTLY INDEPENDENT.  COPY B INTO ZZ, UPDATE Z
00236 Z */
00237 /*        AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ). */
00238 
00239     i__1 = *m;
00240     for (l = 1; l <= i__1; ++l) {
00241 /* L120: */
00242         zz[l] = b[l];
00243     }
00244     i__1 = npp1 + 1;
00245     h12_(&c__2, &npp1, &i__1, m, &a[j * a_dim1 + 1], &c__1, &up, &zz[1], &
00246         c__1, &c__1, &c__1);
00247     ztest = zz[npp1] / a[npp1 + j * a_dim1];
00248 
00249 /*                                     SEE IF ZTEST IS POSITIVE */
00250 
00251     if (ztest > 0.) {
00252         goto L140;
00253     }
00254     }
00255 
00256 /*     REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P. */
00257 /*     RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL */
00258 /*     COEFFS AGAIN. */
00259 
00260     a[npp1 + j * a_dim1] = asave;
00261     w[j] = 0.;
00262     goto L60;
00263 
00264 /*     THE INDEX  J=INDEX(IZ)  HAS BEEN SELECTED TO BE MOVED FROM */
00265 /*     SET Z TO SET P.    UPDATE B,  UPDATE INDICES,  APPLY HOUSEHOLDER */
00266 /*     TRANSFORMATIONS TO COLS IN NEW SET Z,  ZERO SUBDIAGONAL ELTS IN */
00267 /*     COL J,  SET W(J)=0. */
00268 
00269 L140:
00270     i__1 = *m;
00271     for (l = 1; l <= i__1; ++l) {
00272 /* L150: */
00273     b[l] = zz[l];
00274     }
00275 
00276     index[iz] = index[iz1];
00277     index[iz1] = j;
00278     ++iz1;
00279     nsetp = npp1;
00280     ++npp1;
00281 
00282     if (iz1 <= iz2) {
00283     i__1 = iz2;
00284     for (jz = iz1; jz <= i__1; ++jz) {
00285         jj = index[jz];
00286         h12_(&c__2, &nsetp, &npp1, m, &a[j * a_dim1 + 1], &c__1, &up, &a[
00287             jj * a_dim1 + 1], &c__1, mda, &c__1);
00288 /* L160: */
00289     }
00290     }
00291 
00292     if (nsetp != *m) {
00293     i__1 = *m;
00294     for (l = npp1; l <= i__1; ++l) {
00295 /* L180: */
00296         a[l + j * a_dim1] = 0.;
00297     }
00298     }
00299 
00300     w[j] = 0.;
00301 /*                                SOLVE THE TRIANGULAR SYSTEM. */
00302 /*                                STORE THE SOLUTION TEMPORARILY IN ZZ(). 
00303 */
00304     rtnkey = 1;
00305     goto L400;
00306 L200:
00307 
00308 /*                       ******  SECONDARY LOOP BEGINS HERE ****** */
00309 
00310 /*                          ITERATION COUNTER. */
00311 
00312 L210:
00313     ++iter;
00314     if (iter > itmax) {
00315     *mode = 3;
00316     /* The following lines were replaced after the f2c translation */
00317     /* s_wsfe(&io___22); */
00318     /* do_fio(&c__1, " NNLS quitting on iteration count.", 34L); */
00319     /* e_wsfe(); */
00320     fprintf(stdout, "\n NNLS quitting on iteration count.\n");
00321     fflush(stdout);
00322     goto L350;
00323     }
00324 
00325 /*                    SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE. */
00326 /*                                  IF NOT COMPUTE ALPHA. */
00327 
00328     alpha = 2.;
00329     i__1 = nsetp;
00330     for (ip = 1; ip <= i__1; ++ip) {
00331     l = index[ip];
00332     if (zz[ip] <= 0.) {
00333         t = -x[l] / (zz[ip] - x[l]);
00334         if (alpha > t) {
00335         alpha = t;
00336         jj = ip;
00337         }
00338     }
00339 /* L240: */
00340     }
00341 
00342 /*          IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL */
00343 /*          STILL = 2.    IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP. */
00344 
00345     if (alpha == 2.) {
00346     goto L330;
00347     }
00348 
00349 /*          OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO */
00350 /*          INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ. */
00351 
00352     i__1 = nsetp;
00353     for (ip = 1; ip <= i__1; ++ip) {
00354     l = index[ip];
00355     x[l] += alpha * (zz[ip] - x[l]);
00356 /* L250: */
00357     }
00358 
00359 /*        MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I */
00360 /*        FROM SET P TO SET Z. */
00361 
00362     i__ = index[jj];
00363 L260:
00364     x[i__] = 0.;
00365 
00366     if (jj != nsetp) {
00367     ++jj;
00368     i__1 = nsetp;
00369     for (j = jj; j <= i__1; ++j) {
00370         ii = index[j];
00371         index[j - 1] = ii;
00372         g1_(&a[j - 1 + ii * a_dim1], &a[j + ii * a_dim1], &cc, &ss, &a[j 
00373             - 1 + ii * a_dim1]);
00374         a[j + ii * a_dim1] = 0.;
00375         i__2 = *n;
00376         for (l = 1; l <= i__2; ++l) {
00377         if (l != ii) {
00378 
00379 /*                 Apply procedure G2 (CC,SS,A(J-1,L),A(J,
00380 L)) */
00381 
00382             temp = a[j - 1 + l * a_dim1];
00383             a[j - 1 + l * a_dim1] = cc * temp + ss * a[j + l * a_dim1]
00384                 ;
00385             a[j + l * a_dim1] = -ss * temp + cc * a[j + l * a_dim1];
00386         }
00387 /* L270: */
00388         }
00389 
00390 /*                 Apply procedure G2 (CC,SS,B(J-1),B(J)) */
00391 
00392         temp = b[j - 1];
00393         b[j - 1] = cc * temp + ss * b[j];
00394         b[j] = -ss * temp + cc * b[j];
00395 /* L280: */
00396     }
00397     }
00398 
00399     npp1 = nsetp;
00400     --nsetp;
00401     --iz1;
00402     index[iz1] = i__;
00403 
00404 /*        SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE.  THEY SHOULD 
00405 */
00406 /*        BE BECAUSE OF THE WAY ALPHA WAS DETERMINED. */
00407 /*        IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR.  ANY */
00408 /*        THAT ARE NONPOSITIVE WILL BE SET TO ZERO */
00409 /*        AND MOVED FROM SET P TO SET Z. */
00410 
00411     i__1 = nsetp;
00412     for (jj = 1; jj <= i__1; ++jj) {
00413     i__ = index[jj];
00414     if (x[i__] <= 0.) {
00415         goto L260;
00416     }
00417 /* L300: */
00418     }
00419 
00420 /*         COPY B( ) INTO ZZ( ).  THEN SOLVE AGAIN AND LOOP BACK. */
00421 
00422     i__1 = *m;
00423     for (i__ = 1; i__ <= i__1; ++i__) {
00424 /* L310: */
00425     zz[i__] = b[i__];
00426     }
00427     rtnkey = 2;
00428     goto L400;
00429 L320:
00430     goto L210;
00431 /*                      ******  END OF SECONDARY LOOP  ****** */
00432 
00433 L330:
00434     i__1 = nsetp;
00435     for (ip = 1; ip <= i__1; ++ip) {
00436     i__ = index[ip];
00437 /* L340: */
00438     x[i__] = zz[ip];
00439     }
00440 /*        ALL NEW COEFFS ARE POSITIVE.  LOOP BACK TO BEGINNING. */
00441     goto L30;
00442 
00443 /*                        ******  END OF MAIN LOOP  ****** */
00444 
00445 /*                        COME TO HERE FOR TERMINATION. */
00446 /*                     COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR. */
00447 
00448 L350:
00449     sm = 0.;
00450     if (npp1 <= *m) {
00451     i__1 = *m;
00452     for (i__ = npp1; i__ <= i__1; ++i__) {
00453 /* L360: */
00454 /* Computing 2nd power */
00455         d__1 = b[i__];
00456         sm += d__1 * d__1;
00457     }
00458     } else {
00459     i__1 = *n;
00460     for (j = 1; j <= i__1; ++j) {
00461 /* L380: */
00462         w[j] = 0.;
00463     }
00464     }
00465     *rnorm = sqrt(sm);
00466     return 0;
00467 
00468 /*     THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE */
00469 /*     TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ(). */
00470 
00471 L400:
00472     i__1 = nsetp;
00473     for (l = 1; l <= i__1; ++l) {
00474     ip = nsetp + 1 - l;
00475     if (l != 1) {
00476         i__2 = ip;
00477         for (ii = 1; ii <= i__2; ++ii) {
00478         zz[ii] -= a[ii + jj * a_dim1] * zz[ip + 1];
00479 /* L410: */
00480         }
00481     }
00482     jj = index[ip];
00483     zz[ip] /= a[ip + jj * a_dim1];
00484 /* L430: */
00485     }
00486     switch ((int)rtnkey) {
00487     case 1:  goto L200;
00488     case 2:  goto L320;
00489     }
00490 
00491     /* The next line was added after the f2c translation to keep
00492        compilers from complaining about a void return from a non-void
00493        function. */
00494     return 0;
00495 
00496 } /* nnls_ */
00497 
00498 /* Subroutine */ int g1_(a, b, cterm, sterm, sig)
00499 doublereal *a, *b, *cterm, *sterm, *sig;
00500 {
00501     /* System generated locals */
00502     doublereal d__1;
00503 
00504     /* Builtin functions */
00505     /* The following line was commented out after the f2c translation */
00506     /* double sqrt(), d_sign(); */
00507 
00508     /* Local variables */
00509     static doublereal xr, yr;
00510 
00511 
00512 /*     COMPUTE ORTHOGONAL ROTATION MATRIX.. */
00513 
00514 /*  The original version of this code was developed by */
00515 /*  Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory 
00516 */
00517 /*  1973 JUN 12, and published in the book */
00518 /*  "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
00519 /*  Revised FEB 1995 to accompany reprinting of the book by SIAM. */
00520 
00521 /*     COMPUTE.. MATRIX   (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2)) */
00522 /*                        (-S,C)         (-S,C)(B)   (   0          ) */
00523 /*     COMPUTE SIG = SQRT(A**2+B**2) */
00524 /*        SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT */
00525 /*        SIG MAY BE IN THE SAME LOCATION AS A OR B . */
00526 /*     ------------------------------------------------------------------ 
00527 */
00528 /*     ------------------------------------------------------------------ 
00529 */
00530     if (nnls_abs(*a) > nnls_abs(*b)) {
00531     xr = *b / *a;
00532 /* Computing 2nd power */
00533     d__1 = xr;
00534     yr = sqrt(d__1 * d__1 + 1.);
00535     d__1 = 1. / yr;
00536     *cterm = d_sign(&d__1, a);
00537     *sterm = *cterm * xr;
00538     *sig = nnls_abs(*a) * yr;
00539     return 0;
00540     }
00541     if (*b != 0.) {
00542     xr = *a / *b;
00543 /* Computing 2nd power */
00544     d__1 = xr;
00545     yr = sqrt(d__1 * d__1 + 1.);
00546     d__1 = 1. / yr;
00547     *sterm = d_sign(&d__1, b);
00548     *cterm = *sterm * xr;
00549     *sig = nnls_abs(*b) * yr;
00550     return 0;
00551     }
00552     *sig = 0.;
00553     *cterm = 0.;
00554     *sterm = 1.;
00555     return 0;
00556 } /* g1_ */
00557 
00558 /*     SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) */
00559 
00560 /*  CONSTRUCTION AND/OR APPLICATION OF A SINGLE */
00561 /*  HOUSEHOLDER TRANSFORMATION..     Q = I + U*(U**T)/B */
00562 
00563 /*  The original version of this code was developed by */
00564 /*  Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
00565 /*  1973 JUN 12, and published in the book */
00566 /*  "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
00567 /*  Revised FEB 1995 to accompany reprinting of the book by SIAM. */
00568 /*     ------------------------------------------------------------------ */
00569 /*                     Subroutine Arguments */
00570 
00571 /*     MODE   = 1 OR 2   Selects Algorithm H1 to construct and apply a */
00572 /*            Householder transformation, or Algorithm H2 to apply a */
00573 /*            previously constructed transformation. */
00574 /*     LPIVOT IS THE INDEX OF THE PIVOT ELEMENT. */
00575 /*     L1,M   IF L1 .LE. M   THE TRANSFORMATION WILL BE CONSTRUCTED TO */
00576 /*            ZERO ELEMENTS INDEXED FROM L1 THROUGH M.   IF L1 GT. M */
00577 /*            THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION. */
00578 /*     U(),IUE,UP    On entry with MODE = 1, U() contains the pivot */
00579 /*            vector.  IUE is the storage increment between elements. */
00580 /*            On exit when MODE = 1, U() and UP contain quantities */
00581 /*            defining the vector U of the Householder transformation. */
00582 /*            on entry with MODE = 2, U() and UP should contain */
00583 /*            quantities previously computed with MODE = 1.  These will */
00584 /*            not be modified during the entry with MODE = 2. */
00585 /*     C()    ON ENTRY with MODE = 1 or 2, C() CONTAINS A MATRIX WHICH */
00586 /*            WILL BE REGARDED AS A SET OF VECTORS TO WHICH THE */
00587 /*            HOUSEHOLDER TRANSFORMATION IS TO BE APPLIED. */
00588 /*            ON EXIT C() CONTAINS THE SET OF TRANSFORMED VECTORS. */
00589 /*     ICE    STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C(). */
00590 /*     ICV    STORAGE INCREMENT BETWEEN VECTORS IN C(). */
00591 /*     NCV    NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0 */
00592 /*            NO OPERATIONS WILL BE DONE ON C(). */
00593 /*     ------------------------------------------------------------------ */
00594 /* Subroutine */ int h12_(mode, lpivot, l1, m, u, iue, up, c__, ice, icv, ncv)
00595 integer *mode, *lpivot, *l1, *m;
00596 doublereal *u;
00597 integer *iue;
00598 doublereal *up, *c__;
00599 integer *ice, *icv, *ncv;
00600 {
00601     /* System generated locals */
00602     integer u_dim1, u_offset, i__1, i__2;
00603     doublereal d__1, d__2;
00604 
00605     /* Builtin functions */
00606     /* The following line was commented out after the f2c translation */
00607     /* double sqrt(); */
00608 
00609     /* Local variables */
00610     static integer incr;
00611     static doublereal b;
00612     static integer i__, j;
00613     static doublereal clinv;
00614     static integer i2, i3, i4;
00615     static doublereal cl, sm;
00616 
00617 /*     ------------------------------------------------------------------ 
00618 */
00619 /*     double precision U(IUE,M) */
00620 /*     ------------------------------------------------------------------ 
00621 */
00622     /* Parameter adjustments */
00623     u_dim1 = *iue;
00624     u_offset = u_dim1 + 1;
00625     u -= u_offset;
00626     --c__;
00627 
00628     /* Function Body */
00629     if (0 >= *lpivot || *lpivot >= *l1 || *l1 > *m) {
00630     return 0;
00631     }
00632     cl = (d__1 = u[*lpivot * u_dim1 + 1], nnls_abs(d__1));
00633     if (*mode == 2) {
00634     goto L60;
00635     }
00636 /*                            ****** CONSTRUCT THE TRANSFORMATION. ****** 
00637 */
00638     i__1 = *m;
00639     for (j = *l1; j <= i__1; ++j) {
00640 /* L10: */
00641 /* Computing MAX */
00642     d__2 = (d__1 = u[j * u_dim1 + 1], nnls_abs(d__1));
00643     cl = nnls_max(d__2,cl);
00644     }
00645     if (cl <= 0.) {
00646     goto L130;
00647     } else {
00648     goto L20;
00649     }
00650 L20:
00651     clinv = 1. / cl;
00652 /* Computing 2nd power */
00653     d__1 = u[*lpivot * u_dim1 + 1] * clinv;
00654     sm = d__1 * d__1;
00655     i__1 = *m;
00656     for (j = *l1; j <= i__1; ++j) {
00657 /* L30: */
00658 /* Computing 2nd power */
00659     d__1 = u[j * u_dim1 + 1] * clinv;
00660     sm += d__1 * d__1;
00661     }
00662     cl *= sqrt(sm);
00663     if (u[*lpivot * u_dim1 + 1] <= 0.) {
00664     goto L50;
00665     } else {
00666     goto L40;
00667     }
00668 L40:
00669     cl = -cl;
00670 L50:
00671     *up = u[*lpivot * u_dim1 + 1] - cl;
00672     u[*lpivot * u_dim1 + 1] = cl;
00673     goto L70;
00674 /*            ****** APPLY THE TRANSFORMATION  I+U*(U**T)/B  TO C. ****** 
00675 */
00676 
00677 L60:
00678     if (cl <= 0.) {
00679     goto L130;
00680     } else {
00681     goto L70;
00682     }
00683 L70:
00684     if (*ncv <= 0) {
00685     return 0;
00686     }
00687     b = *up * u[*lpivot * u_dim1 + 1];
00688 /*                       B  MUST BE NONPOSITIVE HERE.  IF B = 0., RETURN. 
00689 */
00690 
00691     if (b >= 0.) {
00692     goto L130;
00693     } else {
00694     goto L80;
00695     }
00696 L80:
00697     b = 1. / b;
00698     i2 = 1 - *icv + *ice * (*lpivot - 1);
00699     incr = *ice * (*l1 - *lpivot);
00700     i__1 = *ncv;
00701     for (j = 1; j <= i__1; ++j) {
00702     i2 += *icv;
00703     i3 = i2 + incr;
00704     i4 = i3;
00705     sm = c__[i2] * *up;
00706     i__2 = *m;
00707     for (i__ = *l1; i__ <= i__2; ++i__) {
00708         sm += c__[i3] * u[i__ * u_dim1 + 1];
00709 /* L90: */
00710         i3 += *ice;
00711     }
00712     if (sm != 0.) {
00713         goto L100;
00714     } else {
00715         goto L120;
00716     }
00717 L100:
00718     sm *= b;
00719     c__[i2] += sm * *up;
00720     i__2 = *m;
00721     for (i__ = *l1; i__ <= i__2; ++i__) {
00722         c__[i4] += sm * u[i__ * u_dim1 + 1];
00723 /* L110: */
00724         i4 += *ice;
00725     }
00726 L120:
00727     ;
00728     }
00729 L130:
00730     return 0;
00731 } /* h12_ */
00732 
00733 doublereal diff_(x, y)
00734 doublereal *x, *y;
00735 {
00736     /* System generated locals */
00737     doublereal ret_val;
00738 
00739 
00740 /*  Function used in tests that depend on machine precision. */
00741 
00742 /*  The original version of this code was developed by */
00743 /*  Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory 
00744 */
00745 /*  1973 JUN 7, and published in the book */
00746 /*  "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
00747 /*  Revised FEB 1995 to accompany reprinting of the book by SIAM. */
00748 
00749     ret_val = *x - *y;
00750     return ret_val;
00751 } /* diff_ */
00752 
00753 
00754 /* The following subroutine was added after the f2c translation */
00755 int nnls_c(double* a, const int* mda, const int* m, const int* n, double* b, 
00756      double* x, double* rnorm, double* w, double* zz, int* index, 
00757      int* mode)
00758 {
00759   return (nnls_(a, mda, m, n, b, x, rnorm, w, zz, index, mode));
00760 }