Chris@25: 
Chris@25: #include "nnls.h"
Chris@25: 
Chris@25: typedef int integer;
Chris@25: typedef unsigned int uinteger;
Chris@25: typedef char *address;
Chris@25: typedef short int shortint;
Chris@25: typedef float real;
Chris@25: typedef double doublereal;
Chris@25: 
Chris@25: #define abs(x) ((x) >= 0 ? (x) : -(x))
Chris@25: #define dabs(x) (doublereal)abs(x)
Chris@25: #define min(a,b) ((a) <= (b) ? (a) : (b))
Chris@25: #define max(a,b) ((a) >= (b) ? (a) : (b))
Chris@25: #define dmin(a,b) (doublereal)min(a,b)
Chris@25: #define dmax(a,b) (doublereal)max(a,b)
Chris@25: #define bit_test(a,b)	((a) >> (b) & 1)
Chris@25: #define bit_clear(a,b)	((a) & ~((uinteger)1 << (b)))
Chris@25: #define bit_set(a,b)	((a) |  ((uinteger)1 << (b)))
Chris@25: 
Chris@25: double d_sign(doublereal *a, doublereal *b)
Chris@25: {
Chris@25:     double x;
Chris@25:     x = (*a >= 0 ? *a : - *a);
Chris@25:     return (*b >= 0 ? x : -x);
Chris@25: }
Chris@25: 
Chris@25: /* Table of constant values */
Chris@25: 
Chris@25: static integer c__1 = 1;
Chris@25: static integer c__0 = 0;
Chris@25: static integer c__2 = 2;
Chris@25: 
Chris@25: doublereal diff_(doublereal *x, doublereal *y)
Chris@25: {
Chris@25:     /* System generated locals */
Chris@25:     doublereal ret_val;
Chris@25: 
Chris@25: 
Chris@25: /*  Function used in tests that depend on machine precision. */
Chris@25: 
Chris@25: /*  The original version of this code was developed by */
Chris@25: /*  Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
Chris@25: /*  1973 JUN 7, and published in the book */
Chris@25: /*  "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
Chris@25: /*  Revised FEB 1995 to accompany reprinting of the book by SIAM. */
Chris@25: 
Chris@25:     ret_val = *x - *y;
Chris@25:     return ret_val;
Chris@25: } /* diff_ */
Chris@25: 
Chris@25: /* Subroutine */ int g1_(doublereal *a, doublereal *b, doublereal *cterm, 
Chris@25: 	doublereal *sterm, doublereal *sig)
Chris@25: {
Chris@25:     /* System generated locals */
Chris@25:     doublereal d__1;
Chris@25: 
Chris@25:     /* Builtin functions */
Chris@25:     double sqrt(doublereal), d_sign(doublereal *, doublereal *);
Chris@25: 
Chris@25:     /* Local variables */
Chris@25:     static doublereal xr, yr;
Chris@25: 
Chris@25: 
Chris@25: /*     COMPUTE ORTHOGONAL ROTATION MATRIX.. */
Chris@25: 
Chris@25: /*  The original version of this code was developed by */
Chris@25: /*  Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
Chris@25: /*  1973 JUN 12, and published in the book */
Chris@25: /*  "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
Chris@25: /*  Revised FEB 1995 to accompany reprinting of the book by SIAM. */
Chris@25: 
Chris@25: /*     COMPUTE.. MATRIX   (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2)) */
Chris@25: /*                        (-S,C)         (-S,C)(B)   (   0          ) */
Chris@25: /*     COMPUTE SIG = SQRT(A**2+B**2) */
Chris@25: /*        SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT */
Chris@25: /*        SIG MAY BE IN THE SAME LOCATION AS A OR B . */
Chris@25: /*     ------------------------------------------------------------------ */
Chris@25: /*     ------------------------------------------------------------------ */
Chris@25:     if (abs(*a) > abs(*b)) {
Chris@25: 	xr = *b / *a;
Chris@25: /* Computing 2nd power */
Chris@25: 	d__1 = xr;
Chris@25: 	yr = sqrt(d__1 * d__1 + 1.);
Chris@25: 	d__1 = 1. / yr;
Chris@25: 	*cterm = d_sign(&d__1, a);
Chris@25: 	*sterm = *cterm * xr;
Chris@25: 	*sig = abs(*a) * yr;
Chris@25: 	return 0;
Chris@25:     }
Chris@25:     if (*b != 0.) {
Chris@25: 	xr = *a / *b;
Chris@25: /* Computing 2nd power */
Chris@25: 	d__1 = xr;
Chris@25: 	yr = sqrt(d__1 * d__1 + 1.);
Chris@25: 	d__1 = 1. / yr;
Chris@25: 	*sterm = d_sign(&d__1, b);
Chris@25: 	*cterm = *sterm * xr;
Chris@25: 	*sig = abs(*b) * yr;
Chris@25: 	return 0;
Chris@25:     }
Chris@25:     *sig = 0.;
Chris@25:     *cterm = 0.;
Chris@25:     *sterm = 1.;
Chris@25:     return 0;
Chris@25: } /* g1_ */
Chris@25: 
Chris@25: /*     SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) */
Chris@25: 
Chris@25: /*  CONSTRUCTION AND/OR APPLICATION OF A SINGLE */
Chris@25: /*  HOUSEHOLDER TRANSFORMATION..     Q = I + U*(U**T)/B */
Chris@25: 
Chris@25: /*  The original version of this code was developed by */
Chris@25: /*  Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
Chris@25: /*  1973 JUN 12, and published in the book */
Chris@25: /*  "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
Chris@25: /*  Revised FEB 1995 to accompany reprinting of the book by SIAM. */
Chris@25: /*     ------------------------------------------------------------------ */
Chris@25: /*                     Subroutine Arguments */
Chris@25: 
Chris@25: /*     MODE   = 1 OR 2   Selects Algorithm H1 to construct and apply a */
Chris@25: /*            Householder transformation, or Algorithm H2 to apply a */
Chris@25: /*            previously constructed transformation. */
Chris@25: /*     LPIVOT IS THE INDEX OF THE PIVOT ELEMENT. */
Chris@25: /*     L1,M   IF L1 .LE. M   THE TRANSFORMATION WILL BE CONSTRUCTED TO */
Chris@25: /*            ZERO ELEMENTS INDEXED FROM L1 THROUGH M.   IF L1 GT. M */
Chris@25: /*            THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION. */
Chris@25: /*     U(),IUE,UP    On entry with MODE = 1, U() contains the pivot */
Chris@25: /*            vector.  IUE is the storage increment between elements. */
Chris@25: /*            On exit when MODE = 1, U() and UP contain quantities */
Chris@25: /*            defining the vector U of the Householder transformation. */
Chris@25: /*            on entry with MODE = 2, U() and UP should contain */
Chris@25: /*            quantities previously computed with MODE = 1.  These will */
Chris@25: /*            not be modified during the entry with MODE = 2. */
Chris@25: /*     C()    ON ENTRY with MODE = 1 or 2, C() CONTAINS A MATRIX WHICH */
Chris@25: /*            WILL BE REGARDED AS A SET OF VECTORS TO WHICH THE */
Chris@25: /*            HOUSEHOLDER TRANSFORMATION IS TO BE APPLIED. */
Chris@25: /*            ON EXIT C() CONTAINS THE SET OF TRANSFORMED VECTORS. */
Chris@25: /*     ICE    STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C(). */
Chris@25: /*     ICV    STORAGE INCREMENT BETWEEN VECTORS IN C(). */
Chris@25: /*     NCV    NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0 */
Chris@25: /*            NO OPERATIONS WILL BE DONE ON C(). */
Chris@25: /*     ------------------------------------------------------------------ */
Chris@25: /* Subroutine */ int h12_(integer *mode, integer *lpivot, integer *l1, 
Chris@25: 	integer *m, doublereal *u, integer *iue, doublereal *up, doublereal *
Chris@25: 	c__, integer *ice, integer *icv, integer *ncv)
Chris@25: {
Chris@25:     /* System generated locals */
Chris@25:     integer u_dim1, u_offset, i__1, i__2;
Chris@25:     doublereal d__1, d__2;
Chris@25: 
Chris@25:     /* Builtin functions */
Chris@25:     double sqrt(doublereal);
Chris@25: 
Chris@25:     /* Local variables */
Chris@25:     static doublereal b;
Chris@25:     static integer i__, j, i2, i3, i4;
Chris@25:     static doublereal cl, sm;
Chris@25:     static integer incr;
Chris@25:     static doublereal clinv;
Chris@25: 
Chris@25: /*     ------------------------------------------------------------------ */
Chris@25: /*     double precision U(IUE,M) */
Chris@25: /*     ------------------------------------------------------------------ */
Chris@25:     /* Parameter adjustments */
Chris@25:     u_dim1 = *iue;
Chris@25:     u_offset = 1 + u_dim1;
Chris@25:     u -= u_offset;
Chris@25:     --c__;
Chris@25: 
Chris@25:     /* Function Body */
Chris@25:     if (0 >= *lpivot || *lpivot >= *l1 || *l1 > *m) {
Chris@25: 	return 0;
Chris@25:     }
Chris@25:     cl = (d__1 = u[*lpivot * u_dim1 + 1], abs(d__1));
Chris@25:     if (*mode == 2) {
Chris@25: 	goto L60;
Chris@25:     }
Chris@25: /*                            ****** CONSTRUCT THE TRANSFORMATION. ****** */
Chris@25:     i__1 = *m;
Chris@25:     for (j = *l1; j <= i__1; ++j) {
Chris@25: /* L10: */
Chris@25: /* Computing MAX */
Chris@25: 	d__2 = (d__1 = u[j * u_dim1 + 1], abs(d__1));
Chris@25: 	cl = max(d__2,cl);
Chris@25:     }
Chris@25:     if (cl <= 0.) {
Chris@25: 	goto L130;
Chris@25:     } else {
Chris@25: 	goto L20;
Chris@25:     }
Chris@25: L20:
Chris@25:     clinv = 1. / cl;
Chris@25: /* Computing 2nd power */
Chris@25:     d__1 = u[*lpivot * u_dim1 + 1] * clinv;
Chris@25:     sm = d__1 * d__1;
Chris@25:     i__1 = *m;
Chris@25:     for (j = *l1; j <= i__1; ++j) {
Chris@25: /* L30: */
Chris@25: /* Computing 2nd power */
Chris@25: 	d__1 = u[j * u_dim1 + 1] * clinv;
Chris@25: 	sm += d__1 * d__1;
Chris@25:     }
Chris@25:     cl *= sqrt(sm);
Chris@25:     if (u[*lpivot * u_dim1 + 1] <= 0.) {
Chris@25: 	goto L50;
Chris@25:     } else {
Chris@25: 	goto L40;
Chris@25:     }
Chris@25: L40:
Chris@25:     cl = -cl;
Chris@25: L50:
Chris@25:     *up = u[*lpivot * u_dim1 + 1] - cl;
Chris@25:     u[*lpivot * u_dim1 + 1] = cl;
Chris@25:     goto L70;
Chris@25: /*            ****** APPLY THE TRANSFORMATION  I+U*(U**T)/B  TO C. ****** */
Chris@25: 
Chris@25: L60:
Chris@25:     if (cl <= 0.) {
Chris@25: 	goto L130;
Chris@25:     } else {
Chris@25: 	goto L70;
Chris@25:     }
Chris@25: L70:
Chris@25:     if (*ncv <= 0) {
Chris@25: 	return 0;
Chris@25:     }
Chris@25:     b = *up * u[*lpivot * u_dim1 + 1];
Chris@25: /*                       B  MUST BE NONPOSITIVE HERE.  IF B = 0., RETURN. */
Chris@25: 
Chris@25:     if (b >= 0.) {
Chris@25: 	goto L130;
Chris@25:     } else {
Chris@25: 	goto L80;
Chris@25:     }
Chris@25: L80:
Chris@25:     b = 1. / b;
Chris@25:     i2 = 1 - *icv + *ice * (*lpivot - 1);
Chris@25:     incr = *ice * (*l1 - *lpivot);
Chris@25:     i__1 = *ncv;
Chris@25:     for (j = 1; j <= i__1; ++j) {
Chris@25: 	i2 += *icv;
Chris@25: 	i3 = i2 + incr;
Chris@25: 	i4 = i3;
Chris@25: 	sm = c__[i2] * *up;
Chris@25: 	i__2 = *m;
Chris@25: 	for (i__ = *l1; i__ <= i__2; ++i__) {
Chris@25: 	    sm += c__[i3] * u[i__ * u_dim1 + 1];
Chris@25: /* L90: */
Chris@25: 	    i3 += *ice;
Chris@25: 	}
Chris@25: 	if (sm != 0.) {
Chris@25: 	    goto L100;
Chris@25: 	} else {
Chris@25: 	    goto L120;
Chris@25: 	}
Chris@25: L100:
Chris@25: 	sm *= b;
Chris@25: 	c__[i2] += sm * *up;
Chris@25: 	i__2 = *m;
Chris@25: 	for (i__ = *l1; i__ <= i__2; ++i__) {
Chris@25: 	    c__[i4] += sm * u[i__ * u_dim1 + 1];
Chris@25: /* L110: */
Chris@25: 	    i4 += *ice;
Chris@25: 	}
Chris@25: L120:
Chris@25: 	;
Chris@25:     }
Chris@25: L130:
Chris@25:     return 0;
Chris@25: } /* h12_ */
Chris@25: 
Chris@25: /*     SUBROUTINE NNLS  (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) */
Chris@25: 
Chris@25: /*  Algorithm NNLS: NONNEGATIVE LEAST SQUARES */
Chris@25: 
Chris@25: /*  The original version of this code was developed by */
Chris@25: /*  Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
Chris@25: /*  1973 JUN 15, and published in the book */
Chris@25: /*  "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
Chris@25: /*  Revised FEB 1995 to accompany reprinting of the book by SIAM. */
Chris@25: 
Chris@25: /*     GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B,  COMPUTE AN */
Chris@25: /*     N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM */
Chris@25: 
Chris@25: /*                      A * X = B  SUBJECT TO X .GE. 0 */
Chris@25: /*     ------------------------------------------------------------------ */
Chris@25: /*                     Subroutine Arguments */
Chris@25: 
Chris@25: /*     A(),MDA,M,N     MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE */
Chris@25: /*                     ARRAY, A().   ON ENTRY A() CONTAINS THE M BY N */
Chris@25: /*                     MATRIX, A.           ON EXIT A() CONTAINS */
Chris@25: /*                     THE PRODUCT MATRIX, Q*A , WHERE Q IS AN */
Chris@25: /*                     M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY */
Chris@25: /*                     THIS SUBROUTINE. */
Chris@25: /*     B()     ON ENTRY B() CONTAINS THE M-VECTOR, B.   ON EXIT B() CON- */
Chris@25: /*             TAINS Q*B. */
Chris@25: /*     X()     ON ENTRY X() NEED NOT BE INITIALIZED.  ON EXIT X() WILL */
Chris@25: /*             CONTAIN THE SOLUTION VECTOR. */
Chris@25: /*     RNORM   ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE */
Chris@25: /*             RESIDUAL VECTOR. */
Chris@25: /*     W()     AN N-ARRAY OF WORKING SPACE.  ON EXIT W() WILL CONTAIN */
Chris@25: /*             THE DUAL SOLUTION VECTOR.   W WILL SATISFY W(I) = 0. */
Chris@25: /*             FOR ALL I IN SET P  AND W(I) .LE. 0. FOR ALL I IN SET Z */
Chris@25: /*     ZZ()     AN M-ARRAY OF WORKING SPACE. */
Chris@25: /*     INDEX()     AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N. */
Chris@25: /*                 ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS */
Chris@25: /*                 P AND Z AS FOLLOWS.. */
Chris@25: 
Chris@25: /*                 INDEX(1)   THRU INDEX(NSETP) = SET P. */
Chris@25: /*                 INDEX(IZ1) THRU INDEX(IZ2)   = SET Z. */
Chris@25: /*                 IZ1 = NSETP + 1 = NPP1 */
Chris@25: /*                 IZ2 = N */
Chris@25: /*     MODE    THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING */
Chris@25: /*             MEANINGS. */
Chris@25: /*             1     THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY. */
Chris@25: /*             2     THE DIMENSIONS OF THE PROBLEM ARE BAD. */
Chris@25: /*                   EITHER M .LE. 0 OR N .LE. 0. */
Chris@25: /*             3    ITERATION COUNT EXCEEDED.  MORE THAN 3*N ITERATIONS. */
Chris@25: 
Chris@25: /*     ------------------------------------------------------------------ */
Chris@25: /* Subroutine */ int nnls_(doublereal *a, integer *mda, integer *m, integer *
Chris@25: 	n, doublereal *b, doublereal *x, doublereal *rnorm, doublereal *w, 
Chris@25: 	doublereal *zz, integer *index, integer *mode)
Chris@25: {
Chris@25:     /* System generated locals */
Chris@25:     integer a_dim1, a_offset, i__1, i__2;
Chris@25:     doublereal d__1, d__2;
Chris@25: 
Chris@25:     /* Builtin functions */
Chris@25:     double sqrt(doublereal);
Chris@25: 
Chris@25:     /* Local variables */
Chris@25:     static integer i__, j, l;
Chris@25:     static doublereal t;
Chris@25:     extern /* Subroutine */ int g1_(doublereal *, doublereal *, doublereal *, 
Chris@25: 	    doublereal *, doublereal *);
Chris@25:     static doublereal cc;
Chris@25:     extern /* Subroutine */ int h12_(integer *, integer *, integer *, integer 
Chris@25: 	    *, doublereal *, integer *, doublereal *, doublereal *, integer *,
Chris@25: 	     integer *, integer *);
Chris@25:     static integer ii, jj, ip;
Chris@25:     static doublereal sm;
Chris@25:     static integer iz, jz;
Chris@25:     static doublereal up, ss;
Chris@25:     static integer iz1, iz2, npp1;
Chris@25:     extern doublereal diff_(doublereal *, doublereal *);
Chris@25:     static integer iter;
Chris@25:     static doublereal temp, wmax, alpha, asave;
Chris@25:     static integer itmax, izmax, nsetp;
Chris@25:     static doublereal dummy, unorm, ztest;
Chris@25:     static integer rtnkey;
Chris@25: 
Chris@25: /*     ------------------------------------------------------------------ */
Chris@25: /*     integer INDEX(N) */
Chris@25: /*     double precision A(MDA,N), B(M), W(N), X(N), ZZ(M) */
Chris@25: /*     ------------------------------------------------------------------ */
Chris@25:     /* Parameter adjustments */
Chris@25:     a_dim1 = *mda;
Chris@25:     a_offset = 1 + a_dim1;
Chris@25:     a -= a_offset;
Chris@25:     --b;
Chris@25:     --x;
Chris@25:     --w;
Chris@25:     --zz;
Chris@25:     --index;
Chris@25: 
Chris@25:     /* Function Body */
Chris@25:     *mode = 1;
Chris@25:     if (*m <= 0 || *n <= 0) {
Chris@25: 	*mode = 2;
Chris@25: 	return 0;
Chris@25:     }
Chris@25:     iter = 0;
Chris@25:     itmax = *n * 3;
Chris@25: 
Chris@25: /*                    INITIALIZE THE ARRAYS INDEX() AND X(). */
Chris@25: 
Chris@25:     i__1 = *n;
Chris@25:     for (i__ = 1; i__ <= i__1; ++i__) {
Chris@25: 	x[i__] = 0.;
Chris@25: /* L20: */
Chris@25: 	index[i__] = i__;
Chris@25:     }
Chris@25: 
Chris@25:     iz2 = *n;
Chris@25:     iz1 = 1;
Chris@25:     nsetp = 0;
Chris@25:     npp1 = 1;
Chris@25: /*                             ******  MAIN LOOP BEGINS HERE  ****** */
Chris@25: L30:
Chris@25: /*                  QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION. */
Chris@25: /*                        OR IF M COLS OF A HAVE BEEN TRIANGULARIZED. */
Chris@25: 
Chris@25:     if (iz1 > iz2 || nsetp >= *m) {
Chris@25: 	goto L350;
Chris@25:     }
Chris@25: 
Chris@25: /*         COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W(). */
Chris@25: 
Chris@25:     i__1 = iz2;
Chris@25:     for (iz = iz1; iz <= i__1; ++iz) {
Chris@25: 	j = index[iz];
Chris@25: 	sm = 0.;
Chris@25: 	i__2 = *m;
Chris@25: 	for (l = npp1; l <= i__2; ++l) {
Chris@25: /* L40: */
Chris@25: 	    sm += a[l + j * a_dim1] * b[l];
Chris@25: 	}
Chris@25: 	w[j] = sm;
Chris@25: /* L50: */
Chris@25:     }
Chris@25: /*                                   FIND LARGEST POSITIVE W(J). */
Chris@25: L60:
Chris@25:     wmax = 0.;
Chris@25:     i__1 = iz2;
Chris@25:     for (iz = iz1; iz <= i__1; ++iz) {
Chris@25: 	j = index[iz];
Chris@25: 	if (w[j] > wmax) {
Chris@25: 	    wmax = w[j];
Chris@25: 	    izmax = iz;
Chris@25: 	}
Chris@25: /* L70: */
Chris@25:     }
Chris@25: 
Chris@25: /*             IF WMAX .LE. 0. GO TO TERMINATION. */
Chris@25: /*             THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS. */
Chris@25: 
Chris@25:     if (wmax <= 0.) {
Chris@25: 	goto L350;
Chris@25:     }
Chris@25:     iz = izmax;
Chris@25:     j = index[iz];
Chris@25: 
Chris@25: /*     THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P. */
Chris@25: /*     BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID */
Chris@25: /*     NEAR LINEAR DEPENDENCE. */
Chris@25: 
Chris@25:     asave = a[npp1 + j * a_dim1];
Chris@25:     i__1 = npp1 + 1;
Chris@25:     h12_(&c__1, &npp1, &i__1, m, &a[j * a_dim1 + 1], &c__1, &up, &dummy, &
Chris@25: 	    c__1, &c__1, &c__0);
Chris@25:     unorm = 0.;
Chris@25:     if (nsetp != 0) {
Chris@25: 	i__1 = nsetp;
Chris@25: 	for (l = 1; l <= i__1; ++l) {
Chris@25: /* L90: */
Chris@25: /* Computing 2nd power */
Chris@25: 	    d__1 = a[l + j * a_dim1];
Chris@25: 	    unorm += d__1 * d__1;
Chris@25: 	}
Chris@25:     }
Chris@25:     unorm = sqrt(unorm);
Chris@25:     d__2 = unorm + (d__1 = a[npp1 + j * a_dim1], abs(d__1)) * .01;
Chris@25:     if (diff_(&d__2, &unorm) > 0.) {
Chris@25: 
Chris@25: /*        COL J IS SUFFICIENTLY INDEPENDENT.  COPY B INTO ZZ, UPDATE ZZ */
Chris@25: /*        AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ). */
Chris@25: 
Chris@25: 	i__1 = *m;
Chris@25: 	for (l = 1; l <= i__1; ++l) {
Chris@25: /* L120: */
Chris@25: 	    zz[l] = b[l];
Chris@25: 	}
Chris@25: 	i__1 = npp1 + 1;
Chris@25: 	h12_(&c__2, &npp1, &i__1, m, &a[j * a_dim1 + 1], &c__1, &up, &zz[1], &
Chris@25: 		c__1, &c__1, &c__1);
Chris@25: 	ztest = zz[npp1] / a[npp1 + j * a_dim1];
Chris@25: 
Chris@25: /*                                     SEE IF ZTEST IS POSITIVE */
Chris@25: 
Chris@25: 	if (ztest > 0.) {
Chris@25: 	    goto L140;
Chris@25: 	}
Chris@25:     }
Chris@25: 
Chris@25: /*     REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P. */
Chris@25: /*     RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL */
Chris@25: /*     COEFFS AGAIN. */
Chris@25: 
Chris@25:     a[npp1 + j * a_dim1] = asave;
Chris@25:     w[j] = 0.;
Chris@25:     goto L60;
Chris@25: 
Chris@25: /*     THE INDEX  J=INDEX(IZ)  HAS BEEN SELECTED TO BE MOVED FROM */
Chris@25: /*     SET Z TO SET P.    UPDATE B,  UPDATE INDICES,  APPLY HOUSEHOLDER */
Chris@25: /*     TRANSFORMATIONS TO COLS IN NEW SET Z,  ZERO SUBDIAGONAL ELTS IN */
Chris@25: /*     COL J,  SET W(J)=0. */
Chris@25: 
Chris@25: L140:
Chris@25:     i__1 = *m;
Chris@25:     for (l = 1; l <= i__1; ++l) {
Chris@25: /* L150: */
Chris@25: 	b[l] = zz[l];
Chris@25:     }
Chris@25: 
Chris@25:     index[iz] = index[iz1];
Chris@25:     index[iz1] = j;
Chris@25:     ++iz1;
Chris@25:     nsetp = npp1;
Chris@25:     ++npp1;
Chris@25: 
Chris@25:     if (iz1 <= iz2) {
Chris@25: 	i__1 = iz2;
Chris@25: 	for (jz = iz1; jz <= i__1; ++jz) {
Chris@25: 	    jj = index[jz];
Chris@25: 	    h12_(&c__2, &nsetp, &npp1, m, &a[j * a_dim1 + 1], &c__1, &up, &a[
Chris@25: 		    jj * a_dim1 + 1], &c__1, mda, &c__1);
Chris@25: /* L160: */
Chris@25: 	}
Chris@25:     }
Chris@25: 
Chris@25:     if (nsetp != *m) {
Chris@25: 	i__1 = *m;
Chris@25: 	for (l = npp1; l <= i__1; ++l) {
Chris@25: /* L180: */
Chris@25: 	    a[l + j * a_dim1] = 0.;
Chris@25: 	}
Chris@25:     }
Chris@25: 
Chris@25:     w[j] = 0.;
Chris@25: /*                                SOLVE THE TRIANGULAR SYSTEM. */
Chris@25: /*                                STORE THE SOLUTION TEMPORARILY IN ZZ(). */
Chris@25:     rtnkey = 1;
Chris@25:     goto L400;
Chris@25: L200:
Chris@25: 
Chris@25: /*                       ******  SECONDARY LOOP BEGINS HERE ****** */
Chris@25: 
Chris@25: /*                          ITERATION COUNTER. */
Chris@25: 
Chris@25: L210:
Chris@25:     ++iter;
Chris@25:     if (iter > itmax) {
Chris@25: 	*mode = 3;
Chris@25: /*         write (*,'(/a)') ' NNLS quitting on iteration count.' */
Chris@25: 	goto L350;
Chris@25:     }
Chris@25: 
Chris@25: /*                    SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE. */
Chris@25: /*                                  IF NOT COMPUTE ALPHA. */
Chris@25: 
Chris@25:     alpha = 2.;
Chris@25:     i__1 = nsetp;
Chris@25:     for (ip = 1; ip <= i__1; ++ip) {
Chris@25: 	l = index[ip];
Chris@25: 	if (zz[ip] <= 0.) {
Chris@25: 	    t = -x[l] / (zz[ip] - x[l]);
Chris@25: 	    if (alpha > t) {
Chris@25: 		alpha = t;
Chris@25: 		jj = ip;
Chris@25: 	    }
Chris@25: 	}
Chris@25: /* L240: */
Chris@25:     }
Chris@25: 
Chris@25: /*          IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL */
Chris@25: /*          STILL = 2.    IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP. */
Chris@25: 
Chris@25:     if (alpha == 2.) {
Chris@25: 	goto L330;
Chris@25:     }
Chris@25: 
Chris@25: /*          OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO */
Chris@25: /*          INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ. */
Chris@25: 
Chris@25:     i__1 = nsetp;
Chris@25:     for (ip = 1; ip <= i__1; ++ip) {
Chris@25: 	l = index[ip];
Chris@25: 	x[l] += alpha * (zz[ip] - x[l]);
Chris@25: /* L250: */
Chris@25:     }
Chris@25: 
Chris@25: /*        MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I */
Chris@25: /*        FROM SET P TO SET Z. */
Chris@25: 
Chris@25:     i__ = index[jj];
Chris@25: L260:
Chris@25:     x[i__] = 0.;
Chris@25: 
Chris@25:     if (jj != nsetp) {
Chris@25: 	++jj;
Chris@25: 	i__1 = nsetp;
Chris@25: 	for (j = jj; j <= i__1; ++j) {
Chris@25: 	    ii = index[j];
Chris@25: 	    index[j - 1] = ii;
Chris@25: 	    g1_(&a[j - 1 + ii * a_dim1], &a[j + ii * a_dim1], &cc, &ss, &a[j 
Chris@25: 		    - 1 + ii * a_dim1]);
Chris@25: 	    a[j + ii * a_dim1] = 0.;
Chris@25: 	    i__2 = *n;
Chris@25: 	    for (l = 1; l <= i__2; ++l) {
Chris@25: 		if (l != ii) {
Chris@25: 
Chris@25: /*                 Apply procedure G2 (CC,SS,A(J-1,L),A(J,L)) */
Chris@25: 
Chris@25: 		    temp = a[j - 1 + l * a_dim1];
Chris@25: 		    a[j - 1 + l * a_dim1] = cc * temp + ss * a[j + l * a_dim1]
Chris@25: 			    ;
Chris@25: 		    a[j + l * a_dim1] = -ss * temp + cc * a[j + l * a_dim1];
Chris@25: 		}
Chris@25: /* L270: */
Chris@25: 	    }
Chris@25: 
Chris@25: /*                 Apply procedure G2 (CC,SS,B(J-1),B(J)) */
Chris@25: 
Chris@25: 	    temp = b[j - 1];
Chris@25: 	    b[j - 1] = cc * temp + ss * b[j];
Chris@25: 	    b[j] = -ss * temp + cc * b[j];
Chris@25: /* L280: */
Chris@25: 	}
Chris@25:     }
Chris@25: 
Chris@25:     npp1 = nsetp;
Chris@25:     --nsetp;
Chris@25:     --iz1;
Chris@25:     index[iz1] = i__;
Chris@25: 
Chris@25: /*        SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE.  THEY SHOULD */
Chris@25: /*        BE BECAUSE OF THE WAY ALPHA WAS DETERMINED. */
Chris@25: /*        IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR.  ANY */
Chris@25: /*        THAT ARE NONPOSITIVE WILL BE SET TO ZERO */
Chris@25: /*        AND MOVED FROM SET P TO SET Z. */
Chris@25: 
Chris@25:     i__1 = nsetp;
Chris@25:     for (jj = 1; jj <= i__1; ++jj) {
Chris@25: 	i__ = index[jj];
Chris@25: 	if (x[i__] <= 0.) {
Chris@25: 	    goto L260;
Chris@25: 	}
Chris@25: /* L300: */
Chris@25:     }
Chris@25: 
Chris@25: /*         COPY B( ) INTO ZZ( ).  THEN SOLVE AGAIN AND LOOP BACK. */
Chris@25: 
Chris@25:     i__1 = *m;
Chris@25:     for (i__ = 1; i__ <= i__1; ++i__) {
Chris@25: /* L310: */
Chris@25: 	zz[i__] = b[i__];
Chris@25:     }
Chris@25:     rtnkey = 2;
Chris@25:     goto L400;
Chris@25: L320:
Chris@25:     goto L210;
Chris@25: /*                      ******  END OF SECONDARY LOOP  ****** */
Chris@25: 
Chris@25: L330:
Chris@25:     i__1 = nsetp;
Chris@25:     for (ip = 1; ip <= i__1; ++ip) {
Chris@25: 	i__ = index[ip];
Chris@25: /* L340: */
Chris@25: 	x[i__] = zz[ip];
Chris@25:     }
Chris@25: /*        ALL NEW COEFFS ARE POSITIVE.  LOOP BACK TO BEGINNING. */
Chris@25:     goto L30;
Chris@25: 
Chris@25: /*                        ******  END OF MAIN LOOP  ****** */
Chris@25: 
Chris@25: /*                        COME TO HERE FOR TERMINATION. */
Chris@25: /*                     COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR. */
Chris@25: 
Chris@25: L350:
Chris@25:     sm = 0.;
Chris@25:     if (npp1 <= *m) {
Chris@25: 	i__1 = *m;
Chris@25: 	for (i__ = npp1; i__ <= i__1; ++i__) {
Chris@25: /* L360: */
Chris@25: /* Computing 2nd power */
Chris@25: 	    d__1 = b[i__];
Chris@25: 	    sm += d__1 * d__1;
Chris@25: 	}
Chris@25:     } else {
Chris@25: 	i__1 = *n;
Chris@25: 	for (j = 1; j <= i__1; ++j) {
Chris@25: /* L380: */
Chris@25: 	    w[j] = 0.;
Chris@25: 	}
Chris@25:     }
Chris@25:     *rnorm = sqrt(sm);
Chris@25:     return 0;
Chris@25: 
Chris@25: /*     THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE */
Chris@25: /*     TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ(). */
Chris@25: 
Chris@25: L400:
Chris@25:     i__1 = nsetp;
Chris@25:     for (l = 1; l <= i__1; ++l) {
Chris@25: 	ip = nsetp + 1 - l;
Chris@25: 	if (l != 1) {
Chris@25: 	    i__2 = ip;
Chris@25: 	    for (ii = 1; ii <= i__2; ++ii) {
Chris@25: 		zz[ii] -= a[ii + jj * a_dim1] * zz[ip + 1];
Chris@25: /* L410: */
Chris@25: 	    }
Chris@25: 	}
Chris@25: 	jj = index[ip];
Chris@25: 	zz[ip] /= a[ip + jj * a_dim1];
Chris@25: /* L430: */
Chris@25:     }
Chris@25:     switch (rtnkey) {
Chris@25: 	case 1:  goto L200;
Chris@25: 	case 2:  goto L320;
Chris@25:     }
Chris@25:     return 0;
Chris@25: } /* nnls_ */
Chris@25: