Chris@25: #include "nnls.h"
Chris@25: 
Chris@35: /*
Chris@35:   NNLS-Chroma / Chordino
Chris@35: 
Chris@35:   This file is converted from the Netlib FORTRAN code NNLS.FOR,
Chris@35:   developed by Charles L. Lawson and Richard J. Hanson at Jet
Chris@35:   Propulsion Laboratory 1973 JUN 15, and published in the book
Chris@35:   "SOLVING LEAST SQUARES PROBLEMS", Prentice-Hall, 1974.
Chris@35: 
Chris@35:   Refer to nnls.f for the original code and comments.
Chris@35: */
Chris@35: 
Chris@29: #include <math.h>
Chris@25: 
Chris@29: #define nnls_max(a,b) ((a) >= (b) ? (a) : (b))
Chris@29: #define nnls_abs(x) ((x) >= 0 ? (x) : -(x))
Chris@25: 
Chris@29: float d_sign(float a, float b)
Chris@25: {
Chris@29:   float x;
Chris@29:   x = (a >= 0 ? a : - a);
Chris@29:   return (b >= 0 ? x : -x);
Chris@25: }
Chris@25: 
Chris@25: /* Table of constant values */
Chris@25: 
Chris@29: int c__1 = 1;
Chris@29: int c__0 = 0;
Chris@29: int c__2 = 2;
Chris@25: 
Chris@29: 
Chris@29: int g1(float* a, float* b, float* cterm, float* sterm, float* sig)
Chris@25: {
Chris@29:   /* System generated locals */
Chris@29:   float d;
Chris@25: 
Chris@29:   float xr, yr;
Chris@25: 
Chris@25: 
Chris@29:   if (nnls_abs(*a) > nnls_abs(*b)) {
Chris@29:     xr = *b / *a;
Chris@29:     /* Computing 2nd power */
Chris@29:     d = xr;
Chris@29:     yr = sqrt(d * d + 1.);
Chris@29:     d = 1. / yr;
Chris@29:     *cterm = d_sign(d, *a);
Chris@29:     *sterm = *cterm * xr;
Chris@29:     *sig = nnls_abs(*a) * yr;
Chris@25:     return 0;
Chris@29:   }
Chris@29:   if (*b != 0.) {
Chris@29:     xr = *a / *b;
Chris@29:     /* Computing 2nd power */
Chris@29:     d = xr;
Chris@29:     yr = sqrt(d * d + 1.);
Chris@29:     d = 1. / yr;
Chris@29:     *sterm = d_sign(d, *b);
Chris@29:     *cterm = *sterm * xr;
Chris@29:     *sig = nnls_abs(*b) * yr;
Chris@29:     return 0;
Chris@29:   }
Chris@29:   *sig = 0.;
Chris@29:   *cterm = 0.;
Chris@29:   *sterm = 1.;
Chris@29:   return 0;
Chris@25: } /* g1_ */
Chris@25: 
Chris@29: int h12(int mode, int* lpivot, int* l1,
Chris@29:         int m, float* u, int* iue, float* up, float* c__,
Chris@29:         int* ice, int* icv, int* ncv)
Chris@29: {
Chris@29:   /* System generated locals */
Chris@29:   int u_dim1, u_offset, idx1, idx2;
Chris@29:   float d, d2;
Chris@25: 
Chris@29:   /* Local variables */
Chris@29:   int incr;
Chris@29:   float b;
Chris@29:   int i__, j;
Chris@29:   float clinv;
Chris@29:   int i2, i3, i4;
Chris@29:   float cl, sm;
Chris@25: 
Chris@29:   /*     ------------------------------------------------------------------
Chris@29:    */
Chris@29:   /*     float precision U(IUE,M) */
Chris@29:   /*     ------------------------------------------------------------------
Chris@29:    */
Chris@29:   /* Parameter adjustments */
Chris@29:   u_dim1 = *iue;
Chris@29:   u_offset = u_dim1 + 1;
Chris@29:   u -= u_offset;
Chris@29:   --c__;
Chris@25: 
Chris@29:   /* Function Body */
Chris@29:   if (0 >= *lpivot || *lpivot >= *l1 || *l1 > m) {
Chris@29:     return 0;
Chris@29:   }
Chris@29:   cl = (d = u[*lpivot * u_dim1 + 1], nnls_abs(d));
Chris@29:   if (mode == 2) {
Chris@29:     goto L60;
Chris@29:   }
Chris@29:   /*                            ****** CONSTRUCT THE TRANSFORMATION. ******
Chris@29:    */
Chris@29:   idx1 = m;
Chris@29:   for (j = *l1; j <= idx1; ++j) {
Chris@29:     /* L10: */
Chris@29:     /* Computing MAX */
Chris@29:     d2 = (d = u[j * u_dim1 + 1], nnls_abs(d));
Chris@29:     cl = nnls_max(d2,cl);
Chris@29:   }
Chris@29:   if (cl <= 0.) {
Chris@29:     goto L130;
Chris@29:   } else {
Chris@29:     goto L20;
Chris@29:   }
Chris@29:  L20:
Chris@29:   clinv = 1. / cl;
Chris@29:   /* Computing 2nd power */
Chris@29:   d = u[*lpivot * u_dim1 + 1] * clinv;
Chris@29:   sm = d * d;
Chris@29:   idx1 = m;
Chris@29:   for (j = *l1; j <= idx1; ++j) {
Chris@29:     /* L30: */
Chris@29:     /* Computing 2nd power */
Chris@29:     d = u[j * u_dim1 + 1] * clinv;
Chris@29:     sm += d * d;
Chris@29:   }
Chris@29:   cl *= sqrt(sm);
Chris@29:   if (u[*lpivot * u_dim1 + 1] <= 0.) {
Chris@29:     goto L50;
Chris@29:   } else {
Chris@29:     goto L40;
Chris@29:   }
Chris@29:  L40:
Chris@29:   cl = -cl;
Chris@29:  L50:
Chris@29:   *up = u[*lpivot * u_dim1 + 1] - cl;
Chris@29:   u[*lpivot * u_dim1 + 1] = cl;
Chris@29:   goto L70;
Chris@29:   /*            ****** APPLY THE TRANSFORMATION  I+U*(U**T)/B  TO C. ******
Chris@29:    */
Chris@29: 
Chris@29:  L60:
Chris@29:   if (cl <= 0.) {
Chris@29:     goto L130;
Chris@29:   } else {
Chris@29:     goto L70;
Chris@29:   }
Chris@29:  L70:
Chris@29:   if (*ncv <= 0) {
Chris@29:     return 0;
Chris@29:   }
Chris@29:   b = *up * u[*lpivot * u_dim1 + 1];
Chris@29:   /*                       B  MUST BE NONPOSITIVE HERE.  IF B = 0., RETURN.
Chris@29:    */
Chris@29: 
Chris@29:   if (b >= 0.) {
Chris@29:     goto L130;
Chris@29:   } else {
Chris@29:     goto L80;
Chris@29:   }
Chris@29:  L80:
Chris@29:   b = 1. / b;
Chris@29:   i2 = 1 - *icv + *ice * (*lpivot - 1);
Chris@29:   incr = *ice * (*l1 - *lpivot);
Chris@29:   idx1 = *ncv;
Chris@29:   for (j = 1; j <= idx1; ++j) {
Chris@29:     i2 += *icv;
Chris@29:     i3 = i2 + incr;
Chris@29:     i4 = i3;
Chris@29:     sm = c__[i2] * *up;
Chris@29:     idx2 = m;
Chris@29:     for (i__ = *l1; i__ <= idx2; ++i__) {
Chris@29:       sm += c__[i3] * u[i__ * u_dim1 + 1];
Chris@29:       /* L90: */
Chris@29:       i3 += *ice;
Chris@29:     }
Chris@29:     if (sm != 0.) {
Chris@29:       goto L100;
Chris@29:     } else {
Chris@29:       goto L120;
Chris@29:     }
Chris@29:   L100:
Chris@29:     sm *= b;
Chris@29:     c__[i2] += sm * *up;
Chris@29:     idx2 = m;
Chris@29:     for (i__ = *l1; i__ <= idx2; ++i__) {
Chris@29:       c__[i4] += sm * u[i__ * u_dim1 + 1];
Chris@29:       /* L110: */
Chris@29:       i4 += *ice;
Chris@29:     }
Chris@29:   L120:
Chris@29:     ;
Chris@29:   }
Chris@29:  L130:
Chris@29:   return 0;
Chris@29: } /* h12 */
Chris@29: 
Chris@29: int nnls(float* a,  int mda,  int m,  int n, float* b,
Chris@29:          float* x, float* rnorm, float* w, float* zz, int* index,
Chris@29:          int* mode)
Chris@25: {
Chris@29:   /* System generated locals */
Chris@29:   int a_dim1, a_offset, idx1, idx2;
Chris@29:   float d1, d2;
Chris@25: 
Chris@25: 
Chris@29:   /* Local variables */
Chris@29:   int iter;
Chris@29:   float temp, wmax;
Chris@29:   int i__, j, l;
Chris@29:   float t, alpha, asave;
Chris@176:   int itmax, izmax = 0, nsetp;
Chris@29:   float unorm, ztest, cc;
Chris@29:   float dummy[2];
Chris@29:   int ii, jj, ip;
Chris@29:   float sm;
Chris@29:   int iz, jz;
Chris@29:   float up, ss;
Chris@29:   int rtnkey, iz1, iz2, npp1;
Chris@25: 
Chris@29:   /*     ------------------------------------------------------------------
Chris@29:    */
Chris@29:   /*     int INDEX(N) */
Chris@29:   /*     float precision A(MDA,N), B(M), W(N), X(N), ZZ(M) */
Chris@29:   /*     ------------------------------------------------------------------
Chris@29:    */
Chris@29:   /* Parameter adjustments */
Chris@29:   a_dim1 = mda;
Chris@29:   a_offset = a_dim1 + 1;
Chris@29:   a -= a_offset;
Chris@29:   --b;
Chris@29:   --x;
Chris@29:   --w;
Chris@29:   --zz;
Chris@29:   --index;
Chris@25: 
Chris@29:   /* Function Body */
Chris@29:   *mode = 1;
Chris@29:   if (m <= 0 || n <= 0) {
Chris@29:     *mode = 2;
Chris@29:     return 0;
Chris@29:   }
Chris@29:   iter = 0;
Chris@29:   itmax = n * 3;
Chris@29: 
Chris@29:   /*                    INITIALIZE THE ARRAYS INDEX() AND X(). */
Chris@29: 
Chris@29:   idx1 = n;
Chris@29:   for (i__ = 1; i__ <= idx1; ++i__) {
Chris@29:     x[i__] = 0.;
Chris@29:     /* L20: */
Chris@29:     index[i__] = i__;
Chris@29:   }
Chris@29: 
Chris@29:   iz2 = n;
Chris@29:   iz1 = 1;
Chris@29:   nsetp = 0;
Chris@29:   npp1 = 1;
Chris@29:   /*                             ******  MAIN LOOP BEGINS HERE  ****** */
Chris@29:  L30:
Chris@29:   /*                  QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION.
Chris@29:    */
Chris@29:   /*                        OR IF M COLS OF A HAVE BEEN TRIANGULARIZED. */
Chris@29: 
Chris@29:   if (iz1 > iz2 || nsetp >= m) {
Chris@29:     goto L350;
Chris@29:   }
Chris@29: 
Chris@29:   /*         COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W().
Chris@29:    */
Chris@29: 
Chris@29:   idx1 = iz2;
Chris@29:   for (iz = iz1; iz <= idx1; ++iz) {
Chris@29:     j = index[iz];
Chris@29:     sm = 0.;
Chris@29:     idx2 = m;
Chris@29:     for (l = npp1; l <= idx2; ++l) {
Chris@29:       /* L40: */
Chris@29:       sm += a[l + j * a_dim1] * b[l];
Chris@25:     }
Chris@29:     w[j] = sm;
Chris@29:     /* L50: */
Chris@29:   }
Chris@29:   /*                                   FIND LARGEST POSITIVE W(J). */
Chris@29:  L60:
Chris@29:   wmax = 0.;
Chris@29:   idx1 = iz2;
Chris@29:   for (iz = iz1; iz <= idx1; ++iz) {
Chris@29:     j = index[iz];
Chris@29:     if (w[j] > wmax) {
Chris@29:       wmax = w[j];
Chris@29:       izmax = iz;
Chris@25:     }
Chris@29:     /* L70: */
Chris@29:   }
Chris@29: 
Chris@29:   /*             IF WMAX .LE. 0. GO TO TERMINATION. */
Chris@29:   /*             THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS.
Chris@29:    */
Chris@29: 
Chris@29:   if (wmax <= 0.) {
Chris@29:     goto L350;
Chris@29:   }
Chris@29:   iz = izmax;
Chris@29:   j = index[iz];
Chris@29: 
Chris@29:   /*     THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P. */
Chris@29:   /*     BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID */
Chris@29:   /*     NEAR LINEAR DEPENDENCE. */
Chris@29: 
Chris@29:   asave = a[npp1 + j * a_dim1];
Chris@29:   idx1 = npp1 + 1;
Chris@29:   h12(c__1, &npp1, &idx1, m, &a[j * a_dim1 + 1], &c__1, &up, dummy, &
Chris@29:       c__1, &c__1, &c__0);
Chris@29:   unorm = 0.;
Chris@29:   if (nsetp != 0) {
Chris@29:     idx1 = nsetp;
Chris@29:     for (l = 1; l <= idx1; ++l) {
Chris@29:       /* L90: */
Chris@29:       /* Computing 2nd power */
Chris@29:       d1 = a[l + j * a_dim1];
Chris@29:       unorm += d1 * d1;
Chris@25:     }
Chris@29:   }
Chris@29:   unorm = sqrt(unorm);
Chris@29:   d2 = unorm + (d1 = a[npp1 + j * a_dim1], nnls_abs(d1)) * .01;
Chris@29:   if ((d2- unorm) > 0.) {
Chris@29: 
Chris@29:     /*        COL J IS SUFFICIENTLY INDEPENDENT.  COPY B INTO ZZ, UPDATE Z
Chris@29:               Z */
Chris@29:     /*        AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ). */
Chris@29: 
Chris@29:     idx1 = m;
Chris@29:     for (l = 1; l <= idx1; ++l) {
Chris@29:       /* L120: */
Chris@29:       zz[l] = b[l];
Chris@25:     }
Chris@29:     idx1 = npp1 + 1;
Chris@29:     h12(c__2, &npp1, &idx1, m, &a[j * a_dim1 + 1], &c__1, &up, (zz+1), &
Chris@29:         c__1, &c__1, &c__1);
Chris@29:     ztest = zz[npp1] / a[npp1 + j * a_dim1];
Chris@29: 
Chris@29:     /*                                     SEE IF ZTEST IS POSITIVE */
Chris@29: 
Chris@29:     if (ztest > 0.) {
Chris@29:       goto L140;
Chris@25:     }
Chris@29:   }
Chris@29: 
Chris@29:   /*     REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P. */
Chris@29:   /*     RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL */
Chris@29:   /*     COEFFS AGAIN. */
Chris@29: 
Chris@29:   a[npp1 + j * a_dim1] = asave;
Chris@29:   w[j] = 0.;
Chris@29:   goto L60;
Chris@29: 
Chris@29:   /*     THE INDEX  J=INDEX(IZ)  HAS BEEN SELECTED TO BE MOVED FROM */
Chris@29:   /*     SET Z TO SET P.    UPDATE B,  UPDATE INDICES,  APPLY HOUSEHOLDER */
Chris@29:   /*     TRANSFORMATIONS TO COLS IN NEW SET Z,  ZERO SUBDIAGONAL ELTS IN */
Chris@29:   /*     COL J,  SET W(J)=0. */
Chris@29: 
Chris@29:  L140:
Chris@29:   idx1 = m;
Chris@29:   for (l = 1; l <= idx1; ++l) {
Chris@29:     /* L150: */
Chris@29:     b[l] = zz[l];
Chris@29:   }
Chris@29: 
Chris@29:   index[iz] = index[iz1];
Chris@29:   index[iz1] = j;
Chris@29:   ++iz1;
Chris@29:   nsetp = npp1;
Chris@29:   ++npp1;
Chris@29: 
Chris@29:   if (iz1 <= iz2) {
Chris@29:     idx1 = iz2;
Chris@29:     for (jz = iz1; jz <= idx1; ++jz) {
Chris@29:       jj = index[jz];
Chris@29:       h12(c__2, &nsetp, &npp1, m,
Chris@29:           &a[j * a_dim1 + 1], &c__1, &up,
Chris@29:           &a[jj * a_dim1 + 1], &c__1, &mda, &c__1);
Chris@29:       /* L160: */
Chris@25:     }
Chris@29:   }
Chris@25: 
Chris@29:   if (nsetp != m) {
Chris@29:     idx1 = m;
Chris@29:     for (l = npp1; l <= idx1; ++l) {
Chris@29:       /* L180: */
Chris@29:       // SS: CHECK THIS DAMAGE....
Chris@29:       a[l + j * a_dim1] = 0.;
Chris@25:     }
Chris@29:   }
Chris@29: 
Chris@29:   w[j] = 0.;
Chris@29:   /*                                SOLVE THE TRIANGULAR SYSTEM. */
Chris@29:   /*                                STORE THE SOLUTION TEMPORARILY IN ZZ().
Chris@29:    */
Chris@29:   rtnkey = 1;
Chris@29:   goto L400;
Chris@29:  L200:
Chris@29: 
Chris@29:   /*                       ******  SECONDARY LOOP BEGINS HERE ****** */
Chris@29: 
Chris@29:   /*                          ITERATION COUNTER. */
Chris@29: 
Chris@29:  L210:
Chris@29:   ++iter;
Chris@29:   if (iter > itmax) {
Chris@29:     *mode = 3;
Chris@29:     goto L350;
Chris@29:   }
Chris@29: 
Chris@29:   /*                    SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE. */
Chris@29:   /*                                  IF NOT COMPUTE ALPHA. */
Chris@29: 
Chris@29:   alpha = 2.;
Chris@29:   idx1 = nsetp;
Chris@29:   for (ip = 1; ip <= idx1; ++ip) {
Chris@29:     l = index[ip];
Chris@29:     if (zz[ip] <= 0.) {
Chris@29:       t = -x[l] / (zz[ip] - x[l]);
Chris@29:       if (alpha > t) {
Chris@29:         alpha = t;
Chris@29:         jj = ip;
Chris@29:       }
Chris@25:     }
Chris@29:     /* L240: */
Chris@29:   }
Chris@25: 
Chris@29:   /*          IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL */
Chris@29:   /*          STILL = 2.    IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP. */
Chris@29: 
Chris@29:   if (alpha == 2.) {
Chris@29:     goto L330;
Chris@29:   }
Chris@29: 
Chris@29:   /*          OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO */
Chris@29:   /*          INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ. */
Chris@29: 
Chris@29:   idx1 = nsetp;
Chris@29:   for (ip = 1; ip <= idx1; ++ip) {
Chris@29:     l = index[ip];
Chris@29:     x[l] += alpha * (zz[ip] - x[l]);
Chris@29:     /* L250: */
Chris@29:   }
Chris@29: 
Chris@29:   /*        MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I */
Chris@29:   /*        FROM SET P TO SET Z. */
Chris@29: 
Chris@29:   i__ = index[jj];
Chris@29:  L260:
Chris@29:   x[i__] = 0.;
Chris@29: 
Chris@29:   if (jj != nsetp) {
Chris@29:     ++jj;
Chris@29:     idx1 = nsetp;
Chris@29:     for (j = jj; j <= idx1; ++j) {
Chris@29:       ii = index[j];
Chris@29:       index[j - 1] = ii;
Chris@29:       g1(&a[j - 1 + ii * a_dim1], &a[j + ii * a_dim1],
Chris@29:          &cc, &ss, &a[j - 1 + ii * a_dim1]);
Chris@29:       // SS: CHECK THIS DAMAGE...
Chris@29:       a[j + ii * a_dim1] = 0.;
Chris@29:       idx2 = n;
Chris@29:       for (l = 1; l <= idx2; ++l) {
Chris@29:         if (l != ii) {
Chris@29: 
Chris@29:           /*                 Apply procedure G2 (CC,SS,A(J-1,L),A(J,
Chris@29:                              L)) */
Chris@29: 
Chris@29:           temp = a[j - 1 + l * a_dim1];
Chris@29:           // SS: CHECK THIS DAMAGE
Chris@29:           a[j - 1 + l * a_dim1] = cc * temp + ss * a[j + l * a_dim1];
Chris@29:           a[j + l * a_dim1] = -ss * temp + cc * a[j + l * a_dim1];
Chris@29:         }
Chris@29:         /* L270: */
Chris@29:       }
Chris@29: 
Chris@29:       /*                 Apply procedure G2 (CC,SS,B(J-1),B(J)) */
Chris@29: 
Chris@29:       temp = b[j - 1];
Chris@29:       b[j - 1] = cc * temp + ss * b[j];
Chris@29:       b[j] = -ss * temp + cc * b[j];
Chris@29:       /* L280: */
Chris@25:     }
Chris@29:   }
Chris@29: 
Chris@29:   npp1 = nsetp;
Chris@29:   --nsetp;
Chris@29:   --iz1;
Chris@29:   index[iz1] = i__;
Chris@29: 
Chris@29:   /*        SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE.  THEY SHOULD
Chris@29:    */
Chris@29:   /*        BE BECAUSE OF THE WAY ALPHA WAS DETERMINED. */
Chris@29:   /*        IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR.  ANY */
Chris@29:   /*        THAT ARE NONPOSITIVE WILL BE SET TO ZERO */
Chris@29:   /*        AND MOVED FROM SET P TO SET Z. */
Chris@29: 
Chris@29:   idx1 = nsetp;
Chris@29:   for (jj = 1; jj <= idx1; ++jj) {
Chris@29:     i__ = index[jj];
Chris@29:     if (x[i__] <= 0.) {
Chris@29:       goto L260;
Chris@25:     }
Chris@29:     /* L300: */
Chris@29:   }
Chris@25: 
Chris@29:   /*         COPY B( ) INTO ZZ( ).  THEN SOLVE AGAIN AND LOOP BACK. */
Chris@25: 
Chris@29:   idx1 = m;
Chris@29:   for (i__ = 1; i__ <= idx1; ++i__) {
Chris@29:     /* L310: */
Chris@29:     zz[i__] = b[i__];
Chris@29:   }
Chris@29:   rtnkey = 2;
Chris@29:   goto L400;
Chris@29:  L320:
Chris@29:   goto L210;
Chris@29:   /*                      ******  END OF SECONDARY LOOP  ****** */
Chris@25: 
Chris@29:  L330:
Chris@29:   idx1 = nsetp;
Chris@29:   for (ip = 1; ip <= idx1; ++ip) {
Chris@29:     i__ = index[ip];
Chris@29:     /* L340: */
Chris@29:     x[i__] = zz[ip];
Chris@29:   }
Chris@29:   /*        ALL NEW COEFFS ARE POSITIVE.  LOOP BACK TO BEGINNING. */
Chris@29:   goto L30;
Chris@25: 
Chris@29:   /*                        ******  END OF MAIN LOOP  ****** */
Chris@25: 
Chris@29:   /*                        COME TO HERE FOR TERMINATION. */
Chris@29:   /*                     COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR. */
Chris@25: 
Chris@29:  L350:
Chris@29:   sm = 0.;
Chris@29:   if (npp1 <= m) {
Chris@29:     idx1 = m;
Chris@29:     for (i__ = npp1; i__ <= idx1; ++i__) {
Chris@29:       /* L360: */
Chris@29:       /* Computing 2nd power */
Chris@29:       d1 = b[i__];
Chris@29:       sm += d1 * d1;
Chris@29:     }
Chris@29:   } else {
Chris@29:     idx1 = n;
Chris@29:     for (j = 1; j <= idx1; ++j) {
Chris@29:       /* L380: */
Chris@29:       w[j] = 0.;
Chris@29:     }
Chris@29:   }
Chris@29:   *rnorm = sqrt(sm);
Chris@29:   return 0;
Chris@25: 
Chris@29:   /*     THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE */
Chris@29:   /*     TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ(). */
Chris@25: 
Chris@29:  L400:
Chris@29:   idx1 = nsetp;
Chris@29:   for (l = 1; l <= idx1; ++l) {
Chris@29:     ip = nsetp + 1 - l;
Chris@29:     if (l != 1) {
Chris@29:       idx2 = ip;
Chris@29:       for (ii = 1; ii <= idx2; ++ii) {
Chris@29:         zz[ii] -= a[ii + jj * a_dim1] * zz[ip + 1];
Chris@29:         /* L410: */
Chris@29:       }
Chris@29:     }
Chris@29:     jj = index[ip];
Chris@29:     zz[ip] /= a[ip + jj * a_dim1];
Chris@29:     /* L430: */
Chris@29:   }
Chris@29:   switch ((int)rtnkey) {
Chris@29:   case 1:  goto L200;
Chris@29:   case 2:  goto L320;
Chris@29:   }
Chris@25: 
Chris@29:   return 0;
Chris@25: 
Chris@25: } /* nnls_ */
Chris@25: 
Chris@29: