matthiasm@2: // 	$Id: nnls.cc,v 1.1.2.1 2003/03/06 16:28:07 suvrit Exp $	
matthiasm@2: // File: nnls.cc
matthiasm@2: // Implements the Lawson-Hanson NNLS algorithm
matthiasm@2: // Copied over from nnls.c so i don't ahve copyright on this
matthiasm@2: //...somebody else has.....don't know if this should be under GPL or LGPL or
matthiasm@2: // whatever, but i'll put a lil note here anyways:
matthiasm@2: 
matthiasm@2: // Copyright (C) 2004 Lawson-Hanson
matthiasm@2: // Modifications to adapt to c++ by Suvrit Sra.
matthiasm@2: 
matthiasm@2: // This program is free software; you can redistribute it and/or
matthiasm@2: // modify it under the terms of the GNU General Public License
matthiasm@2: // as published by the Free Software Foundation; either version 2
matthiasm@2: // of the License, or (at your option) any later version.
matthiasm@2: 
matthiasm@2: // This program is distributed in the hope that it will be useful,
matthiasm@2: // but WITHOUT ANY WARRANTY; without even the implied warranty of
matthiasm@2: // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
matthiasm@2: // GNU General Public License for more details.
matthiasm@2: 
matthiasm@2: // You should have received a copy of the GNU General Public License
matthiasm@2: // along with this program; if not, write to the Free Software
matthiasm@2: // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
matthiasm@2: 
matthiasm@2: 
matthiasm@2: #include "nnls.h"
matthiasm@3: #include <iostream>
matthiasm@2: 
matthiasm@3: using namespace std;
matthiasm@3: 
matthiasm@3: float d_sign(float& a, float& b)
matthiasm@2: {
matthiasm@3:   float x;
matthiasm@2:   x = (a >= 0 ? a : - a);
matthiasm@2:   return (b >= 0 ? x : -x);
matthiasm@2: }
matthiasm@2:  
matthiasm@2: /* Table of constant values */
matthiasm@2:  
matthiasm@2: int c__1 = 1;
matthiasm@2: int c__0 = 0;
matthiasm@2: int c__2 = 2;
matthiasm@2: 
matthiasm@2: 
matthiasm@3: int nnls(float* a,  int mda,  int m,  int n, float* b, 
matthiasm@3: 	 float* x, float* rnorm, float* w, float* zz, int* index, 
matthiasm@2: 	 int* mode)
matthiasm@2: {
matthiasm@2:   /* System generated locals */
matthiasm@2:   int a_dim1, a_offset, idx1, idx2;
matthiasm@3:   float d1, d2;
matthiasm@2:  
matthiasm@2:  
matthiasm@2:   /* Local variables */
matthiasm@2:   static int iter;
matthiasm@3:   static float temp, wmax;
matthiasm@2:   static int i__, j, l;
matthiasm@3:   static float t, alpha, asave;
matthiasm@2:   static int itmax, izmax, nsetp;
matthiasm@3:   static float unorm, ztest, cc;
matthiasm@3:   float dummy[2];
matthiasm@2:   static int ii, jj, ip;
matthiasm@3:   static float sm;
matthiasm@2:   static int iz, jz;
matthiasm@3:   static float up, ss;
matthiasm@2:   static int rtnkey, iz1, iz2, npp1;
matthiasm@2:  
matthiasm@2:   /*     ------------------------------------------------------------------ 
matthiasm@2:    */
matthiasm@2:   /*     int INDEX(N) */
matthiasm@3:   /*     float precision A(MDA,N), B(M), W(N), X(N), ZZ(M) */
matthiasm@2:   /*     ------------------------------------------------------------------ 
matthiasm@2:    */
matthiasm@2:   /* Parameter adjustments */
matthiasm@2:   a_dim1 = mda;
matthiasm@2:   a_offset = a_dim1 + 1;
matthiasm@2:   a -= a_offset;
matthiasm@2:   --b;
matthiasm@2:   --x;
matthiasm@2:   --w;
matthiasm@2:   --zz;
matthiasm@2:   --index;
matthiasm@2:  
matthiasm@2:   /* Function Body */
matthiasm@2:   *mode = 1;
matthiasm@2:   if (m <= 0 || n <= 0) {
matthiasm@2:     *mode = 2;
matthiasm@2:     return 0;
matthiasm@2:   }
matthiasm@2:   iter = 0;
matthiasm@3:   itmax = n * 30;
matthiasm@2:  
matthiasm@2:   /*                    INITIALIZE THE ARRAYS INDEX() AND X(). */
matthiasm@2:  
matthiasm@2:   idx1 = n;
matthiasm@2:   for (i__ = 1; i__ <= idx1; ++i__) {
matthiasm@2:     x[i__] = 0.;
matthiasm@2:     /* L20: */
matthiasm@2:     index[i__] = i__;
matthiasm@2:   }
matthiasm@2:  
matthiasm@2:   iz2 = n;
matthiasm@2:   iz1 = 1;
matthiasm@2:   nsetp = 0;
matthiasm@2:   npp1 = 1;
matthiasm@2:   /*                             ******  MAIN LOOP BEGINS HERE  ****** */
matthiasm@2:  L30:
matthiasm@2:   /*                  QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION. 
matthiasm@2:    */
matthiasm@2:   /*                        OR IF M COLS OF A HAVE BEEN TRIANGULARIZED. */
matthiasm@2:  
matthiasm@2:   if (iz1 > iz2 || nsetp >= m) {
matthiasm@2:     goto L350;
matthiasm@2:   }
matthiasm@2:  
matthiasm@2:   /*         COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W(). 
matthiasm@2:    */
matthiasm@2:  
matthiasm@2:   idx1 = iz2;
matthiasm@2:   for (iz = iz1; iz <= idx1; ++iz) {
matthiasm@2:     j = index[iz];
matthiasm@2:     sm = 0.;
matthiasm@2:     idx2 = m;
matthiasm@2:     for (l = npp1; l <= idx2; ++l) {
matthiasm@2:       /* L40: */
matthiasm@2:       sm += a[l + j * a_dim1] * b[l];
matthiasm@2:     }
matthiasm@2:     w[j] = sm;
matthiasm@2:     /* L50: */
matthiasm@2:   }
matthiasm@2:   /*                                   FIND LARGEST POSITIVE W(J). */
matthiasm@2:  L60:
matthiasm@2:   wmax = 0.;
matthiasm@2:   idx1 = iz2;
matthiasm@2:   for (iz = iz1; iz <= idx1; ++iz) {
matthiasm@2:     j = index[iz];
matthiasm@2:     if (w[j] > wmax) {
matthiasm@2:       wmax = w[j];
matthiasm@2:       izmax = iz;
matthiasm@2:     }
matthiasm@2:     /* L70: */
matthiasm@2:   }
matthiasm@2:  
matthiasm@2:   /*             IF WMAX .LE. 0. GO TO TERMINATION. */
matthiasm@2:   /*             THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS. 
matthiasm@2:    */
matthiasm@2:  
matthiasm@2:   if (wmax <= 0.) {
matthiasm@2:     goto L350;
matthiasm@2:   }
matthiasm@2:   iz = izmax;
matthiasm@2:   j = index[iz];
matthiasm@2:  
matthiasm@2:   /*     THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P. */
matthiasm@2:   /*     BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID */
matthiasm@2:   /*     NEAR LINEAR DEPENDENCE. */
matthiasm@2:  
matthiasm@2:   asave = a[npp1 + j * a_dim1];
matthiasm@2:   idx1 = npp1 + 1;
matthiasm@2:   h12(c__1, &npp1, &idx1, m, &a[j * a_dim1 + 1], &c__1, &up, dummy, &
matthiasm@2:       c__1, &c__1, &c__0);
matthiasm@2:   unorm = 0.;
matthiasm@2:   if (nsetp != 0) {
matthiasm@2:     idx1 = nsetp;
matthiasm@2:     for (l = 1; l <= idx1; ++l) {
matthiasm@2:       /* L90: */
matthiasm@2:       /* Computing 2nd power */
matthiasm@2:       d1 = a[l + j * a_dim1];
matthiasm@2:       unorm += d1 * d1;
matthiasm@2:     }
matthiasm@2:   }
matthiasm@2:   unorm = sqrt(unorm);
matthiasm@2:   d2 = unorm + (d1 = a[npp1 + j * a_dim1], nnls_abs(d1)) * .01;
matthiasm@2:   if ((d2- unorm) > 0.) {
matthiasm@2:  
matthiasm@2:     /*        COL J IS SUFFICIENTLY INDEPENDENT.  COPY B INTO ZZ, UPDATE Z
matthiasm@2: 	      Z */
matthiasm@2:     /*        AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ). */
matthiasm@2:  
matthiasm@2:     idx1 = m;
matthiasm@2:     for (l = 1; l <= idx1; ++l) {
matthiasm@2:       /* L120: */
matthiasm@2:       zz[l] = b[l];
matthiasm@2:     }
matthiasm@2:     idx1 = npp1 + 1;
matthiasm@2:     h12(c__2, &npp1, &idx1, m, &a[j * a_dim1 + 1], &c__1, &up, (zz+1), &
matthiasm@2: 	c__1, &c__1, &c__1);
matthiasm@2:     ztest = zz[npp1] / a[npp1 + j * a_dim1];
matthiasm@2:  
matthiasm@2:     /*                                     SEE IF ZTEST IS POSITIVE */
matthiasm@2:  
matthiasm@2:     if (ztest > 0.) {
matthiasm@2:       goto L140;
matthiasm@2:     }
matthiasm@2:   }
matthiasm@2:  
matthiasm@2:   /*     REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P. */
matthiasm@2:   /*     RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL */
matthiasm@2:   /*     COEFFS AGAIN. */
matthiasm@2:  
matthiasm@2:   a[npp1 + j * a_dim1] = asave;
matthiasm@2:   w[j] = 0.;
matthiasm@2:   goto L60;
matthiasm@2:  
matthiasm@2:   /*     THE INDEX  J=INDEX(IZ)  HAS BEEN SELECTED TO BE MOVED FROM */
matthiasm@2:   /*     SET Z TO SET P.    UPDATE B,  UPDATE INDICES,  APPLY HOUSEHOLDER */
matthiasm@2:   /*     TRANSFORMATIONS TO COLS IN NEW SET Z,  ZERO SUBDIAGONAL ELTS IN */
matthiasm@2:   /*     COL J,  SET W(J)=0. */
matthiasm@2:  
matthiasm@2:  L140:
matthiasm@2:   idx1 = m;
matthiasm@2:   for (l = 1; l <= idx1; ++l) {
matthiasm@2:     /* L150: */
matthiasm@2:     b[l] = zz[l];
matthiasm@2:   }
matthiasm@2:  
matthiasm@2:   index[iz] = index[iz1];
matthiasm@2:   index[iz1] = j;
matthiasm@2:   ++iz1;
matthiasm@2:   nsetp = npp1;
matthiasm@2:   ++npp1;
matthiasm@2:  
matthiasm@2:   if (iz1 <= iz2) {
matthiasm@2:     idx1 = iz2;
matthiasm@2:     for (jz = iz1; jz <= idx1; ++jz) {
matthiasm@2:       jj = index[jz];
matthiasm@2:       h12(c__2, &nsetp, &npp1, m, 
matthiasm@2: 	  &a[j * a_dim1 + 1], &c__1, &up, 
matthiasm@2: 	  &a[jj * a_dim1 + 1], &c__1, &mda, &c__1);
matthiasm@2:       /* L160: */
matthiasm@2:     }
matthiasm@2:   }
matthiasm@2:  
matthiasm@2:   if (nsetp != m) {
matthiasm@2:     idx1 = m;
matthiasm@2:     for (l = npp1; l <= idx1; ++l) {
matthiasm@2:       /* L180: */
matthiasm@2:       // SS: CHECK THIS DAMAGE....
matthiasm@2:       a[l + j * a_dim1] = 0.;
matthiasm@2:     }
matthiasm@2:   }
matthiasm@2:  
matthiasm@2:   w[j] = 0.;
matthiasm@2:   /*                                SOLVE THE TRIANGULAR SYSTEM. */
matthiasm@2:   /*                                STORE THE SOLUTION TEMPORARILY IN ZZ(). 
matthiasm@2:    */
matthiasm@2:   rtnkey = 1;
matthiasm@2:   goto L400;
matthiasm@2:  L200:
matthiasm@2:  
matthiasm@2:   /*                       ******  SECONDARY LOOP BEGINS HERE ****** */
matthiasm@2:  
matthiasm@2:   /*                          ITERATION COUNTER. */
matthiasm@2:  
matthiasm@2:  L210:
matthiasm@2:   ++iter;
matthiasm@3: // cerr << iter << endl;
matthiasm@2:   if (iter > itmax) {
matthiasm@2:     *mode = 3;
matthiasm@2:     /* The following lines were replaced after the f2c translation */
matthiasm@2:     /* s_wsfe(&io___22); */
matthiasm@2:     /* do_fio(&c__1, " NNLS quitting on iteration count.", 34L); */
matthiasm@2:     /* e_wsfe(); */
matthiasm@2:     fprintf(stdout, "\n NNLS quitting on iteration count.\n");
matthiasm@2:     fflush(stdout);
matthiasm@2:     goto L350;
matthiasm@2:   }
matthiasm@2:  
matthiasm@2:   /*                    SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE. */
matthiasm@2:   /*                                  IF NOT COMPUTE ALPHA. */
matthiasm@2:  
matthiasm@2:   alpha = 2.;
matthiasm@2:   idx1 = nsetp;
matthiasm@2:   for (ip = 1; ip <= idx1; ++ip) {
matthiasm@2:     l = index[ip];
matthiasm@2:     if (zz[ip] <= 0.) {
matthiasm@2:       t = -x[l] / (zz[ip] - x[l]);
matthiasm@2:       if (alpha > t) {
matthiasm@2: 	alpha = t;
matthiasm@2: 	jj = ip;
matthiasm@2:       }
matthiasm@2:     }
matthiasm@2:     /* L240: */
matthiasm@2:   }
matthiasm@2:  
matthiasm@2:   /*          IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL */
matthiasm@2:   /*          STILL = 2.    IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP. */
matthiasm@2:  
matthiasm@2:   if (alpha == 2.) {
matthiasm@2:     goto L330;
matthiasm@2:   }
matthiasm@2:  
matthiasm@2:   /*          OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO */
matthiasm@2:   /*          INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ. */
matthiasm@2:  
matthiasm@2:   idx1 = nsetp;
matthiasm@2:   for (ip = 1; ip <= idx1; ++ip) {
matthiasm@2:     l = index[ip];
matthiasm@2:     x[l] += alpha * (zz[ip] - x[l]);
matthiasm@2:     /* L250: */
matthiasm@2:   }
matthiasm@2:  
matthiasm@2:   /*        MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I */
matthiasm@2:   /*        FROM SET P TO SET Z. */
matthiasm@2:  
matthiasm@2:   i__ = index[jj];
matthiasm@2:  L260:
matthiasm@2:   x[i__] = 0.;
matthiasm@2:  
matthiasm@2:   if (jj != nsetp) {
matthiasm@2:     ++jj;
matthiasm@2:     idx1 = nsetp;
matthiasm@2:     for (j = jj; j <= idx1; ++j) {
matthiasm@2:       ii = index[j];
matthiasm@2:       index[j - 1] = ii;
matthiasm@2:       g1(&a[j - 1 + ii * a_dim1], &a[j + ii * a_dim1], 
matthiasm@2: 	 &cc, &ss, &a[j - 1 + ii * a_dim1]);
matthiasm@2:       // SS: CHECK THIS DAMAGE...
matthiasm@2:       a[j + ii * a_dim1] = 0.;
matthiasm@2:       idx2 = n;
matthiasm@2:       for (l = 1; l <= idx2; ++l) {
matthiasm@2: 	if (l != ii) {
matthiasm@2:  
matthiasm@2: 	  /*                 Apply procedure G2 (CC,SS,A(J-1,L),A(J,
matthiasm@2: 			     L)) */
matthiasm@2:  
matthiasm@2: 	  temp = a[j - 1 + l * a_dim1];
matthiasm@2: 	  // SS: CHECK THIS DAMAGE
matthiasm@2: 	  a[j - 1 + l * a_dim1] = cc * temp + ss * a[j + l * a_dim1];
matthiasm@2: 	  a[j + l * a_dim1] = -ss * temp + cc * a[j + l * a_dim1];
matthiasm@2: 	}
matthiasm@2: 	/* L270: */
matthiasm@2:       }
matthiasm@2:  
matthiasm@2:       /*                 Apply procedure G2 (CC,SS,B(J-1),B(J)) */
matthiasm@2:  
matthiasm@2:       temp = b[j - 1];
matthiasm@2:       b[j - 1] = cc * temp + ss * b[j];
matthiasm@2:       b[j] = -ss * temp + cc * b[j];
matthiasm@2:       /* L280: */
matthiasm@2:     }
matthiasm@2:   }
matthiasm@2:  
matthiasm@2:   npp1 = nsetp;
matthiasm@2:   --nsetp;
matthiasm@2:   --iz1;
matthiasm@2:   index[iz1] = i__;
matthiasm@3: // cerr << "index[" << iz1 << "]  is " << i__ << endl;
matthiasm@2:  
matthiasm@2:   /*        SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE.  THEY SHOULD 
matthiasm@2:    */
matthiasm@2:   /*        BE BECAUSE OF THE WAY ALPHA WAS DETERMINED. */
matthiasm@2:   /*        IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR.  ANY */
matthiasm@2:   /*        THAT ARE NONPOSITIVE WILL BE SET TO ZERO */
matthiasm@2:   /*        AND MOVED FROM SET P TO SET Z. */
matthiasm@2:  
matthiasm@2:   idx1 = nsetp;
matthiasm@2:   for (jj = 1; jj <= idx1; ++jj) {
matthiasm@2:     i__ = index[jj];
matthiasm@2:     if (x[i__] <= 0.) {
matthiasm@2:       goto L260;
matthiasm@2:     }
matthiasm@2:     /* L300: */
matthiasm@2:   }
matthiasm@2:  
matthiasm@2:   /*         COPY B( ) INTO ZZ( ).  THEN SOLVE AGAIN AND LOOP BACK. */
matthiasm@2:  
matthiasm@2:   idx1 = m;
matthiasm@2:   for (i__ = 1; i__ <= idx1; ++i__) {
matthiasm@2:     /* L310: */
matthiasm@2:     zz[i__] = b[i__];
matthiasm@2:   }
matthiasm@2:   rtnkey = 2;
matthiasm@2:   goto L400;
matthiasm@2:  L320:
matthiasm@2:   goto L210;
matthiasm@2:   /*                      ******  END OF SECONDARY LOOP  ****** */
matthiasm@2:  
matthiasm@2:  L330:
matthiasm@2:   idx1 = nsetp;
matthiasm@2:   for (ip = 1; ip <= idx1; ++ip) {
matthiasm@2:     i__ = index[ip];
matthiasm@2:     /* L340: */
matthiasm@2:     x[i__] = zz[ip];
matthiasm@2:   }
matthiasm@2:   /*        ALL NEW COEFFS ARE POSITIVE.  LOOP BACK TO BEGINNING. */
matthiasm@2:   goto L30;
matthiasm@2:  
matthiasm@2:   /*                        ******  END OF MAIN LOOP  ****** */
matthiasm@2:  
matthiasm@2:   /*                        COME TO HERE FOR TERMINATION. */
matthiasm@2:   /*                     COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR. */
matthiasm@2:  
matthiasm@2:  L350:
matthiasm@2:   sm = 0.;
matthiasm@2:   if (npp1 <= m) {
matthiasm@2:     idx1 = m;
matthiasm@2:     for (i__ = npp1; i__ <= idx1; ++i__) {
matthiasm@2:       /* L360: */
matthiasm@2:       /* Computing 2nd power */
matthiasm@2:       d1 = b[i__];
matthiasm@2:       sm += d1 * d1;
matthiasm@2:     }
matthiasm@2:   } else {
matthiasm@2:     idx1 = n;
matthiasm@2:     for (j = 1; j <= idx1; ++j) {
matthiasm@2:       /* L380: */
matthiasm@2:       w[j] = 0.;
matthiasm@2:     }
matthiasm@2:   }
matthiasm@2:   *rnorm = sqrt(sm);
matthiasm@2:   return 0;
matthiasm@2:  
matthiasm@2:   /*     THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE */
matthiasm@2:   /*     TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ(). */
matthiasm@2:  
matthiasm@2:  L400:
matthiasm@2:   idx1 = nsetp;
matthiasm@2:   for (l = 1; l <= idx1; ++l) {
matthiasm@2:     ip = nsetp + 1 - l;
matthiasm@2:     if (l != 1) {
matthiasm@2:       idx2 = ip;
matthiasm@2:       for (ii = 1; ii <= idx2; ++ii) {
matthiasm@2: 	zz[ii] -= a[ii + jj * a_dim1] * zz[ip + 1];
matthiasm@2: 	/* L410: */
matthiasm@2:       }
matthiasm@2:     }
matthiasm@2:     jj = index[ip];
matthiasm@3: 	// cerr << ip << " " << jj << " " << a_dim1 << endl;
matthiasm@2:     zz[ip] /= a[ip + jj * a_dim1];
matthiasm@2:     /* L430: */
matthiasm@2:   }
matthiasm@2:   switch ((int)rtnkey) {
matthiasm@2:   case 1:  goto L200;
matthiasm@2:   case 2:  goto L320;
matthiasm@2:   }
matthiasm@2:  
matthiasm@2:   /* The next line was added after the f2c translation to keep
matthiasm@2:      compilers from complaining about a void return from a non-void
matthiasm@2:      function. */
matthiasm@2:   return 0;
matthiasm@2:  
matthiasm@2: } /* nnls_ */
matthiasm@2: 
matthiasm@2: 
matthiasm@3: int g1(float* a, float* b, float* cterm, float* sterm, float* sig)
matthiasm@2: {
matthiasm@2:   /* System generated locals */
matthiasm@3:   float d;
matthiasm@2:  
matthiasm@3:   static float xr, yr;
matthiasm@2:  
matthiasm@2:  
matthiasm@2:   if (nnls_abs(*a) > nnls_abs(*b)) {
matthiasm@2:     xr = *b / *a;
matthiasm@2:     /* Computing 2nd power */
matthiasm@2:     d = xr;
matthiasm@2:     yr = sqrt(d * d + 1.);
matthiasm@2:     d = 1. / yr;
matthiasm@2:     *cterm = d_sign(d, *a);
matthiasm@2:     *sterm = *cterm * xr;
matthiasm@2:     *sig = nnls_abs(*a) * yr;
matthiasm@2:     return 0;
matthiasm@2:   }
matthiasm@2:   if (*b != 0.) {
matthiasm@2:     xr = *a / *b;
matthiasm@2:     /* Computing 2nd power */
matthiasm@2:     d = xr;
matthiasm@2:     yr = sqrt(d * d + 1.);
matthiasm@2:     d = 1. / yr;
matthiasm@2:     *sterm = d_sign(d, *b);
matthiasm@2:     *cterm = *sterm * xr;
matthiasm@2:     *sig = nnls_abs(*b) * yr;
matthiasm@2:     return 0;
matthiasm@2:   }
matthiasm@2:   *sig = 0.;
matthiasm@2:   *cterm = 0.;
matthiasm@2:   *sterm = 1.;
matthiasm@2:   return 0;
matthiasm@2: } /* g1_ */
matthiasm@2:  
matthiasm@2: 
matthiasm@2: /* See nnls.h for explanation */
matthiasm@2: int h12(int mode, int* lpivot, int* l1, 
matthiasm@3: 	int m, float* u, int* iue, float* up, float* c__, 
matthiasm@2: 	int* ice, int* icv, int* ncv)
matthiasm@2: {
matthiasm@2:   /* System generated locals */
matthiasm@2:   int u_dim1, u_offset, idx1, idx2;
matthiasm@3:   float d, d2;
matthiasm@2:  
matthiasm@2:   /* Builtin functions */
matthiasm@2:   /* The following line was commented out after the f2c translation */
matthiasm@3:   /* float sqrt(); */
matthiasm@2:  
matthiasm@2:   /* Local variables */
matthiasm@2:   static int incr;
matthiasm@3:   static float b;
matthiasm@2:   static int i__, j;
matthiasm@3:   static float clinv;
matthiasm@2:   static int i2, i3, i4;
matthiasm@3:   static float cl, sm;
matthiasm@2:  
matthiasm@2:   /*     ------------------------------------------------------------------ 
matthiasm@2:    */
matthiasm@3:   /*     float precision U(IUE,M) */
matthiasm@2:   /*     ------------------------------------------------------------------ 
matthiasm@2:    */
matthiasm@2:   /* Parameter adjustments */
matthiasm@2:   u_dim1 = *iue;
matthiasm@2:   u_offset = u_dim1 + 1;
matthiasm@2:   u -= u_offset;
matthiasm@2:   --c__;
matthiasm@2:  
matthiasm@2:   /* Function Body */
matthiasm@2:   if (0 >= *lpivot || *lpivot >= *l1 || *l1 > m) {
matthiasm@2:     return 0;
matthiasm@2:   }
matthiasm@2:   cl = (d = u[*lpivot * u_dim1 + 1], nnls_abs(d));
matthiasm@2:   if (mode == 2) {
matthiasm@2:     goto L60;
matthiasm@2:   }
matthiasm@2:   /*                            ****** CONSTRUCT THE TRANSFORMATION. ****** 
matthiasm@2:    */
matthiasm@2:   idx1 = m;
matthiasm@2:   for (j = *l1; j <= idx1; ++j) {
matthiasm@2:     /* L10: */
matthiasm@2:     /* Computing MAX */
matthiasm@2:     d2 = (d = u[j * u_dim1 + 1], nnls_abs(d));
matthiasm@2:     cl = nnls_max(d2,cl);
matthiasm@2:   }
matthiasm@2:   if (cl <= 0.) {
matthiasm@2:     goto L130;
matthiasm@2:   } else {
matthiasm@2:     goto L20;
matthiasm@2:   }
matthiasm@2:  L20:
matthiasm@2:   clinv = 1. / cl;
matthiasm@2:   /* Computing 2nd power */
matthiasm@2:   d = u[*lpivot * u_dim1 + 1] * clinv;
matthiasm@2:   sm = d * d;
matthiasm@2:   idx1 = m;
matthiasm@2:   for (j = *l1; j <= idx1; ++j) {
matthiasm@2:     /* L30: */
matthiasm@2:     /* Computing 2nd power */
matthiasm@2:     d = u[j * u_dim1 + 1] * clinv;
matthiasm@2:     sm += d * d;
matthiasm@2:   }
matthiasm@2:   cl *= sqrt(sm);
matthiasm@2:   if (u[*lpivot * u_dim1 + 1] <= 0.) {
matthiasm@2:     goto L50;
matthiasm@2:   } else {
matthiasm@2:     goto L40;
matthiasm@2:   }
matthiasm@2:  L40:
matthiasm@2:   cl = -cl;
matthiasm@2:  L50:
matthiasm@2:   *up = u[*lpivot * u_dim1 + 1] - cl;
matthiasm@2:   u[*lpivot * u_dim1 + 1] = cl;
matthiasm@2:   goto L70;
matthiasm@2:   /*            ****** APPLY THE TRANSFORMATION  I+U*(U**T)/B  TO C. ****** 
matthiasm@2:    */
matthiasm@2:  
matthiasm@2:  L60:
matthiasm@2:   if (cl <= 0.) {
matthiasm@2:     goto L130;
matthiasm@2:   } else {
matthiasm@2:     goto L70;
matthiasm@2:   }
matthiasm@2:  L70:
matthiasm@2:   if (*ncv <= 0) {
matthiasm@2:     return 0;
matthiasm@2:   }
matthiasm@2:   b = *up * u[*lpivot * u_dim1 + 1];
matthiasm@2:   /*                       B  MUST BE NONPOSITIVE HERE.  IF B = 0., RETURN. 
matthiasm@2:    */
matthiasm@2:  
matthiasm@2:   if (b >= 0.) {
matthiasm@2:     goto L130;
matthiasm@2:   } else {
matthiasm@2:     goto L80;
matthiasm@2:   }
matthiasm@2:  L80:
matthiasm@2:   b = 1. / b;
matthiasm@2:   i2 = 1 - *icv + *ice * (*lpivot - 1);
matthiasm@2:   incr = *ice * (*l1 - *lpivot);
matthiasm@2:   idx1 = *ncv;
matthiasm@2:   for (j = 1; j <= idx1; ++j) {
matthiasm@2:     i2 += *icv;
matthiasm@2:     i3 = i2 + incr;
matthiasm@2:     i4 = i3;
matthiasm@2:     sm = c__[i2] * *up;
matthiasm@2:     idx2 = m;
matthiasm@2:     for (i__ = *l1; i__ <= idx2; ++i__) {
matthiasm@2:       sm += c__[i3] * u[i__ * u_dim1 + 1];
matthiasm@2:       /* L90: */
matthiasm@2:       i3 += *ice;
matthiasm@2:     }
matthiasm@2:     if (sm != 0.) {
matthiasm@2:       goto L100;
matthiasm@2:     } else {
matthiasm@2:       goto L120;
matthiasm@2:     }
matthiasm@2:   L100:
matthiasm@2:     sm *= b;
matthiasm@2:     c__[i2] += sm * *up;
matthiasm@2:     idx2 = m;
matthiasm@2:     for (i__ = *l1; i__ <= idx2; ++i__) {
matthiasm@2:       c__[i4] += sm * u[i__ * u_dim1 + 1];
matthiasm@2:       /* L110: */
matthiasm@2:       i4 += *ice;
matthiasm@2:     }
matthiasm@2:   L120:
matthiasm@2:     ;
matthiasm@2:   }
matthiasm@2:  L130:
matthiasm@2:   return 0;
matthiasm@2: } /* h12 */
matthiasm@2: