Mercurial > hg > nnls-chroma
changeset 22:444c344681f3 matthiasm-plugin
* Rather than worry about provenance of C versions, why not just use the FORTRAN? (We can back this out if it doesn't go well with the build scripts)
| author | Chris Cannam |
|---|---|
| date | Thu, 21 Oct 2010 11:58:28 +0100 |
| parents | 1e2d0db15920 |
| children | 93c836cfb8c5 |
| files | Makefile.cc-linux NNLSChroma.cpp nnls.cc nnls.cpp nnls.f nnls.h |
| diffstat | 6 files changed, 497 insertions(+), 1335 deletions(-) [+] |
line wrap: on
line diff
--- a/Makefile.cc-linux Mon Oct 18 09:35:26 2010 +0000 +++ b/Makefile.cc-linux Thu Oct 21 11:58:28 2010 +0100 @@ -31,7 +31,6 @@ NNLS_DIR = ../tsnnls/tsnnls - ## Uncomment these for an OS/X native build using command-line tools: #CXXFLAGS = -I$(VAMP_SDK_DIR) -I$(LAPACK_DIR) -I$(FFT_DIR) -I$(NNLS_DIR) -Wall -fPIC -g #PLUGIN_EXT = .dylib @@ -80,6 +79,9 @@ $(PLUGIN): $(PLUGIN_CODE_OBJECTS) $(CXX) -o $@ $^ $(LDFLAGS) +nnls.o: nnls.f + gfortran -std=legacy -fPIC -Wall -c $< + clean: rm -f *.o
--- a/NNLSChroma.cpp Mon Oct 18 09:35:26 2010 +0000 +++ b/NNLSChroma.cpp Thu Oct 21 11:58:28 2010 +0100 @@ -1199,7 +1199,7 @@ f6.hasTimestamp = true; f6.timestamp = f2.timestamp; - float b[256]; + double b[256]; bool some_b_greater_zero = false; float sumb = 0; @@ -1233,7 +1233,7 @@ } } else { - float x[84+1000]; + double x[84+1000]; for (int i = 1; i < 1084; ++i) x[i] = 1.0; vector<int> signifIndex; int index=0; @@ -1247,20 +1247,22 @@ f3.values.push_back(0); // fill the values, change later index++; } - float rnorm; - float w[84+1000]; - float zz[84+1000]; + double rnorm; + double w[84+1000]; + double zz[84+1000]; int indx[84+1000]; int mode; int dictsize = 256*signifIndex.size(); // cerr << "dictsize is " << dictsize << "and values size" << f3.values.size()<< endl; - float *curr_dict = new float[dictsize]; + double *curr_dict = new double[dictsize]; for (unsigned iNote = 0; iNote < signifIndex.size(); ++iNote) { for (unsigned iBin = 0; iBin < 256; iBin++) { curr_dict[iNote * 256 + iBin] = 1.0 * m_dict[signifIndex[iNote] * 256 + iBin]; } } - nnls(curr_dict, nNote, nNote, signifIndex.size(), b, x, &rnorm, w, zz, indx, &mode); + int sz = signifIndex.size(); + int nn = nNote; + NNLS(curr_dict, &nn, &nn, &sz, b, x, &rnorm, w, zz, indx, &mode); delete [] curr_dict; for (unsigned iNote = 0; iNote < signifIndex.size(); ++iNote) { f3.values[signifIndex[iNote]] = x[iNote];
--- a/nnls.cc Mon Oct 18 09:35:26 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,591 +0,0 @@ -// $Id: nnls.cc,v 1.1.1.1 2003/04/02 22:06:19 suvrit Exp $ -// File: nnls.cc -// Implements the Lawson-Hanson NNLS algorithm - -#include "nnls.h" - -float d_sign(float& a, float& b) -{ - float x; - x = (a >= 0 ? a : - a); - return (b >= 0 ? x : -x); -} - -/* Table of constant values */ - -int c__1 = 1; -int c__0 = 0; -int c__2 = 2; - - -int nnls(float* a, int mda, int m, int n, float* b, - float* x, float* rnorm, float* w, float* zz, int* index, - int* mode) -{ - /* System generated locals */ - int a_dim1, a_offset, idx1, idx2; - float d1, d2; - - - /* Local variables */ - static int iter; - static float temp, wmax; - static int i__, j, l; - static float t, alpha, asave; - static int itmax, izmax, nsetp; - static float unorm, ztest, cc; - float dummy[2]; - static int ii, jj, ip; - static float sm; - static int iz, jz; - static float up, ss; - static int rtnkey, iz1, iz2, npp1; - - /* ------------------------------------------------------------------ - */ - /* int INDEX(N) */ - /* float precision A(MDA,N), B(M), W(N), X(N), ZZ(M) */ - /* ------------------------------------------------------------------ - */ - /* Parameter adjustments */ - a_dim1 = mda; - a_offset = a_dim1 + 1; - a -= a_offset; - --b; - --x; - --w; - --zz; - --index; - - /* Function Body */ - *mode = 1; - if (m <= 0 || n <= 0) { - *mode = 2; - return 0; - } - iter = 0; - itmax = n * 3; - - /* INITIALIZE THE ARRAYS INDEX() AND X(). */ - - idx1 = n; - for (i__ = 1; i__ <= idx1; ++i__) { - x[i__] = 0.; - /* L20: */ - index[i__] = i__; - } - - iz2 = n; - iz1 = 1; - nsetp = 0; - npp1 = 1; - /* ****** MAIN LOOP BEGINS HERE ****** */ - L30: - /* QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION. - */ - /* OR IF M COLS OF A HAVE BEEN TRIANGULARIZED. */ - - if (iz1 > iz2 || nsetp >= m) { - goto L350; - } - - /* COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W(). - */ - - idx1 = iz2; - for (iz = iz1; iz <= idx1; ++iz) { - j = index[iz]; - sm = 0.; - idx2 = m; - for (l = npp1; l <= idx2; ++l) { - /* L40: */ - sm += a[l + j * a_dim1] * b[l]; - } - w[j] = sm; - /* L50: */ - } - /* FIND LARGEST POSITIVE W(J). */ - L60: - wmax = 0.; - idx1 = iz2; - for (iz = iz1; iz <= idx1; ++iz) { - j = index[iz]; - if (w[j] > wmax) { - wmax = w[j]; - izmax = iz; - } - /* L70: */ - } - - /* IF WMAX .LE. 0. GO TO TERMINATION. */ - /* THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS. - */ - - if (wmax <= 0.) { - goto L350; - } - iz = izmax; - j = index[iz]; - - /* THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P. */ - /* BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID */ - /* NEAR LINEAR DEPENDENCE. */ - - asave = a[npp1 + j * a_dim1]; - idx1 = npp1 + 1; - h12(c__1, &npp1, &idx1, m, &a[j * a_dim1 + 1], &c__1, &up, dummy, & - c__1, &c__1, &c__0); - unorm = 0.; - if (nsetp != 0) { - idx1 = nsetp; - for (l = 1; l <= idx1; ++l) { - /* L90: */ - /* Computing 2nd power */ - d1 = a[l + j * a_dim1]; - unorm += d1 * d1; - } - } - unorm = sqrt(unorm); - d2 = unorm + (d1 = a[npp1 + j * a_dim1], nnls_abs(d1)) * .01; - if ((d2- unorm) > 0.) { - - /* COL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE Z - Z */ - /* AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ). */ - - idx1 = m; - for (l = 1; l <= idx1; ++l) { - /* L120: */ - zz[l] = b[l]; - } - idx1 = npp1 + 1; - h12(c__2, &npp1, &idx1, m, &a[j * a_dim1 + 1], &c__1, &up, (zz+1), & - c__1, &c__1, &c__1); - ztest = zz[npp1] / a[npp1 + j * a_dim1]; - - /* SEE IF ZTEST IS POSITIVE */ - - if (ztest > 0.) { - goto L140; - } - } - - /* REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P. */ - /* RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL */ - /* COEFFS AGAIN. */ - - a[npp1 + j * a_dim1] = asave; - w[j] = 0.; - goto L60; - - /* THE INDEX J=INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM */ - /* SET Z TO SET P. UPDATE B, UPDATE INDICES, APPLY HOUSEHOLDER */ - /* TRANSFORMATIONS TO COLS IN NEW SET Z, ZERO SUBDIAGONAL ELTS IN */ - /* COL J, SET W(J)=0. */ - - L140: - idx1 = m; - for (l = 1; l <= idx1; ++l) { - /* L150: */ - b[l] = zz[l]; - } - - index[iz] = index[iz1]; - index[iz1] = j; - ++iz1; - nsetp = npp1; - ++npp1; - - if (iz1 <= iz2) { - idx1 = iz2; - for (jz = iz1; jz <= idx1; ++jz) { - jj = index[jz]; - h12(c__2, &nsetp, &npp1, m, - &a[j * a_dim1 + 1], &c__1, &up, - &a[jj * a_dim1 + 1], &c__1, &mda, &c__1); - /* L160: */ - } - } - - if (nsetp != m) { - idx1 = m; - for (l = npp1; l <= idx1; ++l) { - /* L180: */ - // SS: CHECK THIS DAMAGE.... - a[l + j * a_dim1] = 0.; - } - } - - w[j] = 0.; - /* SOLVE THE TRIANGULAR SYSTEM. */ - /* STORE THE SOLUTION TEMPORARILY IN ZZ(). - */ - rtnkey = 1; - goto L400; - L200: - - /* ****** SECONDARY LOOP BEGINS HERE ****** */ - - /* ITERATION COUNTER. */ - - L210: - ++iter; - if (iter > itmax) { - *mode = 3; - /* The following lines were replaced after the f2c translation */ - /* s_wsfe(&io___22); */ - /* do_fio(&c__1, " NNLS quitting on iteration count.", 34L); */ - /* e_wsfe(); */ - fprintf(stdout, "\n NNLS quitting on iteration count.\n"); - fflush(stdout); - goto L350; - } - - /* SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE. */ - /* IF NOT COMPUTE ALPHA. */ - - alpha = 2.; - idx1 = nsetp; - for (ip = 1; ip <= idx1; ++ip) { - l = index[ip]; - if (zz[ip] <= 0.) { - t = -x[l] / (zz[ip] - x[l]); - if (alpha > t) { - alpha = t; - jj = ip; - } - } - /* L240: */ - } - - /* IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL */ - /* STILL = 2. IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP. */ - - if (alpha == 2.) { - goto L330; - } - - /* OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO */ - /* INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ. */ - - idx1 = nsetp; - for (ip = 1; ip <= idx1; ++ip) { - l = index[ip]; - x[l] += alpha * (zz[ip] - x[l]); - /* L250: */ - } - - /* MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I */ - /* FROM SET P TO SET Z. */ - - i__ = index[jj]; - L260: - x[i__] = 0.; - - if (jj != nsetp) { - ++jj; - idx1 = nsetp; - for (j = jj; j <= idx1; ++j) { - ii = index[j]; - index[j - 1] = ii; - g1(&a[j - 1 + ii * a_dim1], &a[j + ii * a_dim1], - &cc, &ss, &a[j - 1 + ii * a_dim1]); - // SS: CHECK THIS DAMAGE... - a[j + ii * a_dim1] = 0.; - idx2 = n; - for (l = 1; l <= idx2; ++l) { - if (l != ii) { - - /* Apply procedure G2 (CC,SS,A(J-1,L),A(J, - L)) */ - - temp = a[j - 1 + l * a_dim1]; - // SS: CHECK THIS DAMAGE - a[j - 1 + l * a_dim1] = cc * temp + ss * a[j + l * a_dim1]; - a[j + l * a_dim1] = -ss * temp + cc * a[j + l * a_dim1]; - } - /* L270: */ - } - - /* Apply procedure G2 (CC,SS,B(J-1),B(J)) */ - - temp = b[j - 1]; - b[j - 1] = cc * temp + ss * b[j]; - b[j] = -ss * temp + cc * b[j]; - /* L280: */ - } - } - - npp1 = nsetp; - --nsetp; - --iz1; - index[iz1] = i__; - - /* SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE. THEY SHOULD - */ - /* BE BECAUSE OF THE WAY ALPHA WAS DETERMINED. */ - /* IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR. ANY */ - /* THAT ARE NONPOSITIVE WILL BE SET TO ZERO */ - /* AND MOVED FROM SET P TO SET Z. */ - - idx1 = nsetp; - for (jj = 1; jj <= idx1; ++jj) { - i__ = index[jj]; - if (x[i__] <= 0.) { - goto L260; - } - /* L300: */ - } - - /* COPY B( ) INTO ZZ( ). THEN SOLVE AGAIN AND LOOP BACK. */ - - idx1 = m; - for (i__ = 1; i__ <= idx1; ++i__) { - /* L310: */ - zz[i__] = b[i__]; - } - rtnkey = 2; - goto L400; - L320: - goto L210; - /* ****** END OF SECONDARY LOOP ****** */ - - L330: - idx1 = nsetp; - for (ip = 1; ip <= idx1; ++ip) { - i__ = index[ip]; - /* L340: */ - x[i__] = zz[ip]; - } - /* ALL NEW COEFFS ARE POSITIVE. LOOP BACK TO BEGINNING. */ - goto L30; - - /* ****** END OF MAIN LOOP ****** */ - - /* COME TO HERE FOR TERMINATION. */ - /* COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR. */ - - L350: - sm = 0.; - if (npp1 <= m) { - idx1 = m; - for (i__ = npp1; i__ <= idx1; ++i__) { - /* L360: */ - /* Computing 2nd power */ - d1 = b[i__]; - sm += d1 * d1; - } - } else { - idx1 = n; - for (j = 1; j <= idx1; ++j) { - /* L380: */ - w[j] = 0.; - } - } - *rnorm = sqrt(sm); - return 0; - - /* THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE */ - /* TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ(). */ - - L400: - idx1 = nsetp; - for (l = 1; l <= idx1; ++l) { - ip = nsetp + 1 - l; - if (l != 1) { - idx2 = ip; - for (ii = 1; ii <= idx2; ++ii) { - zz[ii] -= a[ii + jj * a_dim1] * zz[ip + 1]; - /* L410: */ - } - } - jj = index[ip]; - zz[ip] /= a[ip + jj * a_dim1]; - /* L430: */ - } - switch ((int)rtnkey) { - case 1: goto L200; - case 2: goto L320; - } - - /* The next line was added after the f2c translation to keep - compilers from complaining about a void return from a non-void - function. */ - return 0; - -} /* nnls_ */ - - -int g1(float* a, float* b, float* cterm, float* sterm, float* sig) -{ - /* System generated locals */ - float d; - - static float xr, yr; - - - if (nnls_abs(*a) > nnls_abs(*b)) { - xr = *b / *a; - /* Computing 2nd power */ - d = xr; - yr = sqrt(d * d + 1.); - d = 1. / yr; - *cterm = d_sign(d, *a); - *sterm = *cterm * xr; - *sig = nnls_abs(*a) * yr; - return 0; - } - if (*b != 0.) { - xr = *a / *b; - /* Computing 2nd power */ - d = xr; - yr = sqrt(d * d + 1.); - d = 1. / yr; - *sterm = d_sign(d, *b); - *cterm = *sterm * xr; - *sig = nnls_abs(*b) * yr; - return 0; - } - *sig = 0.; - *cterm = 0.; - *sterm = 1.; - return 0; -} /* g1_ */ - - -/* See nnls.h for explanation */ -int h12(int mode, int* lpivot, int* l1, - int m, float* u, int* iue, float* up, float* c__, - int* ice, int* icv, int* ncv) -{ - /* System generated locals */ - int u_dim1, u_offset, idx1, idx2; - float d, d2; - - /* Builtin functions */ - /* The following line was commented out after the f2c translation */ - /* float sqrt(); */ - - /* Local variables */ - static int incr; - static float b; - static int i__, j; - static float clinv; - static int i2, i3, i4; - static float cl, sm; - - /* ------------------------------------------------------------------ - */ - /* float precision U(IUE,M) */ - /* ------------------------------------------------------------------ - */ - /* Parameter adjustments */ - u_dim1 = *iue; - u_offset = u_dim1 + 1; - u -= u_offset; - --c__; - - /* Function Body */ - if (0 >= *lpivot || *lpivot >= *l1 || *l1 > m) { - return 0; - } - cl = (d = u[*lpivot * u_dim1 + 1], nnls_abs(d)); - if (mode == 2) { - goto L60; - } - /* ****** CONSTRUCT THE TRANSFORMATION. ****** - */ - idx1 = m; - for (j = *l1; j <= idx1; ++j) { - /* L10: */ - /* Computing MAX */ - d2 = (d = u[j * u_dim1 + 1], nnls_abs(d)); - cl = nnls_max(d2,cl); - } - if (cl <= 0.) { - goto L130; - } else { - goto L20; - } - L20: - clinv = 1. / cl; - /* Computing 2nd power */ - d = u[*lpivot * u_dim1 + 1] * clinv; - sm = d * d; - idx1 = m; - for (j = *l1; j <= idx1; ++j) { - /* L30: */ - /* Computing 2nd power */ - d = u[j * u_dim1 + 1] * clinv; - sm += d * d; - } - cl *= sqrt(sm); - if (u[*lpivot * u_dim1 + 1] <= 0.) { - goto L50; - } else { - goto L40; - } - L40: - cl = -cl; - L50: - *up = u[*lpivot * u_dim1 + 1] - cl; - u[*lpivot * u_dim1 + 1] = cl; - goto L70; - /* ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C. ****** - */ - - L60: - if (cl <= 0.) { - goto L130; - } else { - goto L70; - } - L70: - if (*ncv <= 0) { - return 0; - } - b = *up * u[*lpivot * u_dim1 + 1]; - /* B MUST BE NONPOSITIVE HERE. IF B = 0., RETURN. - */ - - if (b >= 0.) { - goto L130; - } else { - goto L80; - } - L80: - b = 1. / b; - i2 = 1 - *icv + *ice * (*lpivot - 1); - incr = *ice * (*l1 - *lpivot); - idx1 = *ncv; - for (j = 1; j <= idx1; ++j) { - i2 += *icv; - i3 = i2 + incr; - i4 = i3; - sm = c__[i2] * *up; - idx2 = m; - for (i__ = *l1; i__ <= idx2; ++i__) { - sm += c__[i3] * u[i__ * u_dim1 + 1]; - /* L90: */ - i3 += *ice; - } - if (sm != 0.) { - goto L100; - } else { - goto L120; - } - L100: - sm *= b; - c__[i2] += sm * *up; - idx2 = m; - for (i__ = *l1; i__ <= idx2; ++i__) { - c__[i4] += sm * u[i__ * u_dim1 + 1]; - /* L110: */ - i4 += *ice; - } - L120: - ; - } - L130: - return 0; -} /* h12 */
--- a/nnls.cpp Mon Oct 18 09:35:26 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,619 +0,0 @@ -// $Id: nnls.cc,v 1.1.2.1 2003/03/06 16:28:07 suvrit Exp $ -// File: nnls.cc -// Implements the Lawson-Hanson NNLS algorithm -// Copied over from nnls.c so i don't ahve copyright on this -//...somebody else has.....don't know if this should be under GPL or LGPL or -// whatever, but i'll put a lil note here anyways: - -// Copyright (C) 2004 Lawson-Hanson -// Modifications to adapt to c++ by Suvrit Sra. - -// This program is free software; you can redistribute it and/or -// modify it under the terms of the GNU General Public License -// as published by the Free Software Foundation; either version 2 -// of the License, or (at your option) any later version. - -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. - -// You should have received a copy of the GNU General Public License -// along with this program; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. - - -#include "nnls.h" -#include <iostream> - -using namespace std; - -float d_sign(float& a, float& b) -{ - float x; - x = (a >= 0 ? a : - a); - return (b >= 0 ? x : -x); -} - -/* Table of constant values */ - -int c__1 = 1; -int c__0 = 0; -int c__2 = 2; - - -int nnls(float* a, int mda, int m, int n, float* b, - float* x, float* rnorm, float* w, float* zz, int* index, - int* mode) -{ - /* System generated locals */ - int a_dim1, a_offset, idx1, idx2; - float d1, d2; - - - /* Local variables */ - static int iter; - static float temp, wmax; - static int i__, j, l; - static float t, alpha, asave; - static int itmax, izmax, nsetp; - static float unorm, ztest, cc; - float dummy[2]; - static int ii, jj, ip; - static float sm; - static int iz, jz; - static float up, ss; - static int rtnkey, iz1, iz2, npp1; - - /* ------------------------------------------------------------------ - */ - /* int INDEX(N) */ - /* float precision A(MDA,N), B(M), W(N), X(N), ZZ(M) */ - /* ------------------------------------------------------------------ - */ - /* Parameter adjustments */ - a_dim1 = mda; - a_offset = a_dim1 + 1; - a -= a_offset; - --b; - --x; - --w; - --zz; - --index; - - /* Function Body */ - *mode = 1; - if (m <= 0 || n <= 0) { - *mode = 2; - return 0; - } - iter = 0; - itmax = n * 30; - - /* INITIALIZE THE ARRAYS INDEX() AND X(). */ - - idx1 = n; - for (i__ = 1; i__ <= idx1; ++i__) { - x[i__] = 0.; - /* L20: */ - index[i__] = i__; - } - - iz2 = n; - iz1 = 1; - nsetp = 0; - npp1 = 1; - /* ****** MAIN LOOP BEGINS HERE ****** */ - L30: - /* QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION. - */ - /* OR IF M COLS OF A HAVE BEEN TRIANGULARIZED. */ - - if (iz1 > iz2 || nsetp >= m) { - goto L350; - } - - /* COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W(). - */ - - idx1 = iz2; - for (iz = iz1; iz <= idx1; ++iz) { - j = index[iz]; - sm = 0.; - idx2 = m; - for (l = npp1; l <= idx2; ++l) { - /* L40: */ - sm += a[l + j * a_dim1] * b[l]; - } - w[j] = sm; - /* L50: */ - } - /* FIND LARGEST POSITIVE W(J). */ - L60: - wmax = 0.; - idx1 = iz2; - for (iz = iz1; iz <= idx1; ++iz) { - j = index[iz]; - if (w[j] > wmax) { - wmax = w[j]; - izmax = iz; - } - /* L70: */ - } - - /* IF WMAX .LE. 0. GO TO TERMINATION. */ - /* THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS. - */ - - if (wmax <= 0.) { - goto L350; - } - iz = izmax; - j = index[iz]; - - /* THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P. */ - /* BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID */ - /* NEAR LINEAR DEPENDENCE. */ - - asave = a[npp1 + j * a_dim1]; - idx1 = npp1 + 1; - h12(c__1, &npp1, &idx1, m, &a[j * a_dim1 + 1], &c__1, &up, dummy, & - c__1, &c__1, &c__0); - unorm = 0.; - if (nsetp != 0) { - idx1 = nsetp; - for (l = 1; l <= idx1; ++l) { - /* L90: */ - /* Computing 2nd power */ - d1 = a[l + j * a_dim1]; - unorm += d1 * d1; - } - } - unorm = sqrt(unorm); - d2 = unorm + (d1 = a[npp1 + j * a_dim1], nnls_abs(d1)) * .01; - if ((d2- unorm) > 0.) { - - /* COL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE Z - Z */ - /* AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ). */ - - idx1 = m; - for (l = 1; l <= idx1; ++l) { - /* L120: */ - zz[l] = b[l]; - } - idx1 = npp1 + 1; - h12(c__2, &npp1, &idx1, m, &a[j * a_dim1 + 1], &c__1, &up, (zz+1), & - c__1, &c__1, &c__1); - ztest = zz[npp1] / a[npp1 + j * a_dim1]; - - /* SEE IF ZTEST IS POSITIVE */ - - if (ztest > 0.) { - goto L140; - } - } - - /* REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P. */ - /* RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL */ - /* COEFFS AGAIN. */ - - a[npp1 + j * a_dim1] = asave; - w[j] = 0.; - goto L60; - - /* THE INDEX J=INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM */ - /* SET Z TO SET P. UPDATE B, UPDATE INDICES, APPLY HOUSEHOLDER */ - /* TRANSFORMATIONS TO COLS IN NEW SET Z, ZERO SUBDIAGONAL ELTS IN */ - /* COL J, SET W(J)=0. */ - - L140: - idx1 = m; - for (l = 1; l <= idx1; ++l) { - /* L150: */ - b[l] = zz[l]; - } - - index[iz] = index[iz1]; - index[iz1] = j; - ++iz1; - nsetp = npp1; - ++npp1; - - if (iz1 <= iz2) { - idx1 = iz2; - for (jz = iz1; jz <= idx1; ++jz) { - jj = index[jz]; - h12(c__2, &nsetp, &npp1, m, - &a[j * a_dim1 + 1], &c__1, &up, - &a[jj * a_dim1 + 1], &c__1, &mda, &c__1); - /* L160: */ - } - } - - if (nsetp != m) { - idx1 = m; - for (l = npp1; l <= idx1; ++l) { - /* L180: */ - // SS: CHECK THIS DAMAGE.... - a[l + j * a_dim1] = 0.; - } - } - - w[j] = 0.; - /* SOLVE THE TRIANGULAR SYSTEM. */ - /* STORE THE SOLUTION TEMPORARILY IN ZZ(). - */ - rtnkey = 1; - goto L400; - L200: - - /* ****** SECONDARY LOOP BEGINS HERE ****** */ - - /* ITERATION COUNTER. */ - - L210: - ++iter; -// cerr << iter << endl; - if (iter > itmax) { - *mode = 3; - /* The following lines were replaced after the f2c translation */ - /* s_wsfe(&io___22); */ - /* do_fio(&c__1, " NNLS quitting on iteration count.", 34L); */ - /* e_wsfe(); */ - fprintf(stdout, "\n NNLS quitting on iteration count.\n"); - fflush(stdout); - goto L350; - } - - /* SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE. */ - /* IF NOT COMPUTE ALPHA. */ - - alpha = 2.; - idx1 = nsetp; - for (ip = 1; ip <= idx1; ++ip) { - l = index[ip]; - if (zz[ip] <= 0.) { - t = -x[l] / (zz[ip] - x[l]); - if (alpha > t) { - alpha = t; - jj = ip; - } - } - /* L240: */ - } - - /* IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL */ - /* STILL = 2. IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP. */ - - if (alpha == 2.) { - goto L330; - } - - /* OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO */ - /* INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ. */ - - idx1 = nsetp; - for (ip = 1; ip <= idx1; ++ip) { - l = index[ip]; - x[l] += alpha * (zz[ip] - x[l]); - /* L250: */ - } - - /* MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I */ - /* FROM SET P TO SET Z. */ - - i__ = index[jj]; - L260: - x[i__] = 0.; - - if (jj != nsetp) { - ++jj; - idx1 = nsetp; - for (j = jj; j <= idx1; ++j) { - ii = index[j]; - index[j - 1] = ii; - g1(&a[j - 1 + ii * a_dim1], &a[j + ii * a_dim1], - &cc, &ss, &a[j - 1 + ii * a_dim1]); - // SS: CHECK THIS DAMAGE... - a[j + ii * a_dim1] = 0.; - idx2 = n; - for (l = 1; l <= idx2; ++l) { - if (l != ii) { - - /* Apply procedure G2 (CC,SS,A(J-1,L),A(J, - L)) */ - - temp = a[j - 1 + l * a_dim1]; - // SS: CHECK THIS DAMAGE - a[j - 1 + l * a_dim1] = cc * temp + ss * a[j + l * a_dim1]; - a[j + l * a_dim1] = -ss * temp + cc * a[j + l * a_dim1]; - } - /* L270: */ - } - - /* Apply procedure G2 (CC,SS,B(J-1),B(J)) */ - - temp = b[j - 1]; - b[j - 1] = cc * temp + ss * b[j]; - b[j] = -ss * temp + cc * b[j]; - /* L280: */ - } - } - - npp1 = nsetp; - --nsetp; - --iz1; - index[iz1] = i__; -// cerr << "index[" << iz1 << "] is " << i__ << endl; - - /* SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE. THEY SHOULD - */ - /* BE BECAUSE OF THE WAY ALPHA WAS DETERMINED. */ - /* IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR. ANY */ - /* THAT ARE NONPOSITIVE WILL BE SET TO ZERO */ - /* AND MOVED FROM SET P TO SET Z. */ - - idx1 = nsetp; - for (jj = 1; jj <= idx1; ++jj) { - i__ = index[jj]; - if (x[i__] <= 0.) { - goto L260; - } - /* L300: */ - } - - /* COPY B( ) INTO ZZ( ). THEN SOLVE AGAIN AND LOOP BACK. */ - - idx1 = m; - for (i__ = 1; i__ <= idx1; ++i__) { - /* L310: */ - zz[i__] = b[i__]; - } - rtnkey = 2; - goto L400; - L320: - goto L210; - /* ****** END OF SECONDARY LOOP ****** */ - - L330: - idx1 = nsetp; - for (ip = 1; ip <= idx1; ++ip) { - i__ = index[ip]; - /* L340: */ - x[i__] = zz[ip]; - } - /* ALL NEW COEFFS ARE POSITIVE. LOOP BACK TO BEGINNING. */ - goto L30; - - /* ****** END OF MAIN LOOP ****** */ - - /* COME TO HERE FOR TERMINATION. */ - /* COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR. */ - - L350: - sm = 0.; - if (npp1 <= m) { - idx1 = m; - for (i__ = npp1; i__ <= idx1; ++i__) { - /* L360: */ - /* Computing 2nd power */ - d1 = b[i__]; - sm += d1 * d1; - } - } else { - idx1 = n; - for (j = 1; j <= idx1; ++j) { - /* L380: */ - w[j] = 0.; - } - } - *rnorm = sqrt(sm); - return 0; - - /* THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE */ - /* TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ(). */ - - L400: - idx1 = nsetp; - for (l = 1; l <= idx1; ++l) { - ip = nsetp + 1 - l; - if (l != 1) { - idx2 = ip; - for (ii = 1; ii <= idx2; ++ii) { - zz[ii] -= a[ii + jj * a_dim1] * zz[ip + 1]; - /* L410: */ - } - } - jj = index[ip]; - // cerr << ip << " " << jj << " " << a_dim1 << endl; - zz[ip] /= a[ip + jj * a_dim1]; - /* L430: */ - } - switch ((int)rtnkey) { - case 1: goto L200; - case 2: goto L320; - } - - /* The next line was added after the f2c translation to keep - compilers from complaining about a void return from a non-void - function. */ - return 0; - -} /* nnls_ */ - - -int g1(float* a, float* b, float* cterm, float* sterm, float* sig) -{ - /* System generated locals */ - float d; - - static float xr, yr; - - - if (nnls_abs(*a) > nnls_abs(*b)) { - xr = *b / *a; - /* Computing 2nd power */ - d = xr; - yr = sqrt(d * d + 1.); - d = 1. / yr; - *cterm = d_sign(d, *a); - *sterm = *cterm * xr; - *sig = nnls_abs(*a) * yr; - return 0; - } - if (*b != 0.) { - xr = *a / *b; - /* Computing 2nd power */ - d = xr; - yr = sqrt(d * d + 1.); - d = 1. / yr; - *sterm = d_sign(d, *b); - *cterm = *sterm * xr; - *sig = nnls_abs(*b) * yr; - return 0; - } - *sig = 0.; - *cterm = 0.; - *sterm = 1.; - return 0; -} /* g1_ */ - - -/* See nnls.h for explanation */ -int h12(int mode, int* lpivot, int* l1, - int m, float* u, int* iue, float* up, float* c__, - int* ice, int* icv, int* ncv) -{ - /* System generated locals */ - int u_dim1, u_offset, idx1, idx2; - float d, d2; - - /* Builtin functions */ - /* The following line was commented out after the f2c translation */ - /* float sqrt(); */ - - /* Local variables */ - static int incr; - static float b; - static int i__, j; - static float clinv; - static int i2, i3, i4; - static float cl, sm; - - /* ------------------------------------------------------------------ - */ - /* float precision U(IUE,M) */ - /* ------------------------------------------------------------------ - */ - /* Parameter adjustments */ - u_dim1 = *iue; - u_offset = u_dim1 + 1; - u -= u_offset; - --c__; - - /* Function Body */ - if (0 >= *lpivot || *lpivot >= *l1 || *l1 > m) { - return 0; - } - cl = (d = u[*lpivot * u_dim1 + 1], nnls_abs(d)); - if (mode == 2) { - goto L60; - } - /* ****** CONSTRUCT THE TRANSFORMATION. ****** - */ - idx1 = m; - for (j = *l1; j <= idx1; ++j) { - /* L10: */ - /* Computing MAX */ - d2 = (d = u[j * u_dim1 + 1], nnls_abs(d)); - cl = nnls_max(d2,cl); - } - if (cl <= 0.) { - goto L130; - } else { - goto L20; - } - L20: - clinv = 1. / cl; - /* Computing 2nd power */ - d = u[*lpivot * u_dim1 + 1] * clinv; - sm = d * d; - idx1 = m; - for (j = *l1; j <= idx1; ++j) { - /* L30: */ - /* Computing 2nd power */ - d = u[j * u_dim1 + 1] * clinv; - sm += d * d; - } - cl *= sqrt(sm); - if (u[*lpivot * u_dim1 + 1] <= 0.) { - goto L50; - } else { - goto L40; - } - L40: - cl = -cl; - L50: - *up = u[*lpivot * u_dim1 + 1] - cl; - u[*lpivot * u_dim1 + 1] = cl; - goto L70; - /* ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C. ****** - */ - - L60: - if (cl <= 0.) { - goto L130; - } else { - goto L70; - } - L70: - if (*ncv <= 0) { - return 0; - } - b = *up * u[*lpivot * u_dim1 + 1]; - /* B MUST BE NONPOSITIVE HERE. IF B = 0., RETURN. - */ - - if (b >= 0.) { - goto L130; - } else { - goto L80; - } - L80: - b = 1. / b; - i2 = 1 - *icv + *ice * (*lpivot - 1); - incr = *ice * (*l1 - *lpivot); - idx1 = *ncv; - for (j = 1; j <= idx1; ++j) { - i2 += *icv; - i3 = i2 + incr; - i4 = i3; - sm = c__[i2] * *up; - idx2 = m; - for (i__ = *l1; i__ <= idx2; ++i__) { - sm += c__[i3] * u[i__ * u_dim1 + 1]; - /* L90: */ - i3 += *ice; - } - if (sm != 0.) { - goto L100; - } else { - goto L120; - } - L100: - sm *= b; - c__[i2] += sm * *up; - idx2 = m; - for (i__ = *l1; i__ <= idx2; ++i__) { - c__[i4] += sm * u[i__ * u_dim1 + 1]; - /* L110: */ - i4 += *ice; - } - L120: - ; - } - L130: - return 0; -} /* h12 */ -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nnls.f Thu Oct 21 11:58:28 2010 +0100 @@ -0,0 +1,479 @@ + + double precision FUNCTION DIFF(X,Y) +c +c Function used in tests that depend on machine precision. +c +c The original version of this code was developed by +c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory +c 1973 JUN 7, and published in the book +c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. +c Revised FEB 1995 to accompany reprinting of the book by SIAM. +C + double precision X, Y + DIFF=X-Y + RETURN + END + + + SUBROUTINE G1 (A,B,CTERM,STERM,SIG) +c +C COMPUTE ORTHOGONAL ROTATION MATRIX.. +c +c The original version of this code was developed by +c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory +c 1973 JUN 12, and published in the book +c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. +c Revised FEB 1995 to accompany reprinting of the book by SIAM. +C +C COMPUTE.. MATRIX (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2)) +C (-S,C) (-S,C)(B) ( 0 ) +C COMPUTE SIG = SQRT(A**2+B**2) +C SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT +C SIG MAY BE IN THE SAME LOCATION AS A OR B . +C ------------------------------------------------------------------ + double precision A, B, CTERM, ONE, SIG, STERM, XR, YR, ZERO + parameter(ONE = 1.0d0, ZERO = 0.0d0) +C ------------------------------------------------------------------ + if (abs(A) .gt. abs(B)) then + XR=B/A + YR=sqrt(ONE+XR**2) + CTERM=sign(ONE/YR,A) + STERM=CTERM*XR + SIG=abs(A)*YR + RETURN + endif + + if (B .ne. ZERO) then + XR=A/B + YR=sqrt(ONE+XR**2) + STERM=sign(ONE/YR,B) + CTERM=STERM*XR + SIG=abs(B)*YR + RETURN + endif + + SIG=ZERO + CTERM=ZERO + STERM=ONE + RETURN + END + + +C SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) +C +C CONSTRUCTION AND/OR APPLICATION OF A SINGLE +C HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B +C +c The original version of this code was developed by +c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory +c 1973 JUN 12, and published in the book +c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. +c Revised FEB 1995 to accompany reprinting of the book by SIAM. +C ------------------------------------------------------------------ +c Subroutine Arguments +c +C MODE = 1 OR 2 Selects Algorithm H1 to construct and apply a +c Householder transformation, or Algorithm H2 to apply a +c previously constructed transformation. +C LPIVOT IS THE INDEX OF THE PIVOT ELEMENT. +C L1,M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO +C ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M +C THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION. +C U(),IUE,UP On entry with MODE = 1, U() contains the pivot +c vector. IUE is the storage increment between elements. +c On exit when MODE = 1, U() and UP contain quantities +c defining the vector U of the Householder transformation. +c on entry with MODE = 2, U() and UP should contain +c quantities previously computed with MODE = 1. These will +c not be modified during the entry with MODE = 2. +C C() ON ENTRY with MODE = 1 or 2, C() CONTAINS A MATRIX WHICH +c WILL BE REGARDED AS A SET OF VECTORS TO WHICH THE +c HOUSEHOLDER TRANSFORMATION IS TO BE APPLIED. +c ON EXIT C() CONTAINS THE SET OF TRANSFORMED VECTORS. +C ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C(). +C ICV STORAGE INCREMENT BETWEEN VECTORS IN C(). +C NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0 +C NO OPERATIONS WILL BE DONE ON C(). +C ------------------------------------------------------------------ + SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) +C ------------------------------------------------------------------ + integer I, I2, I3, I4, ICE, ICV, INCR, IUE, J + integer L1, LPIVOT, M, MODE, NCV + double precision B, C(*), CL, CLINV, ONE, SM +c double precision U(IUE,M) + double precision U(IUE,*) + double precision UP + parameter(ONE = 1.0d0) +C ------------------------------------------------------------------ + IF (0.GE.LPIVOT.OR.LPIVOT.GE.L1.OR.L1.GT.M) RETURN + CL=abs(U(1,LPIVOT)) + IF (MODE.EQ.2) GO TO 60 +C ****** CONSTRUCT THE TRANSFORMATION. ****** + DO 10 J=L1,M + 10 CL=MAX(abs(U(1,J)),CL) + IF (CL) 130,130,20 + 20 CLINV=ONE/CL + SM=(U(1,LPIVOT)*CLINV)**2 + DO 30 J=L1,M + 30 SM=SM+(U(1,J)*CLINV)**2 + CL=CL*SQRT(SM) + IF (U(1,LPIVOT)) 50,50,40 + 40 CL=-CL + 50 UP=U(1,LPIVOT)-CL + U(1,LPIVOT)=CL + GO TO 70 +C ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C. ****** +C + 60 IF (CL) 130,130,70 + 70 IF (NCV.LE.0) RETURN + B= UP*U(1,LPIVOT) +C B MUST BE NONPOSITIVE HERE. IF B = 0., RETURN. +C + IF (B) 80,130,130 + 80 B=ONE/B + I2=1-ICV+ICE*(LPIVOT-1) + INCR=ICE*(L1-LPIVOT) + DO 120 J=1,NCV + I2=I2+ICV + I3=I2+INCR + I4=I3 + SM=C(I2)*UP + DO 90 I=L1,M + SM=SM+C(I3)*U(1,I) + 90 I3=I3+ICE + IF (SM) 100,120,100 + 100 SM=SM*B + C(I2)=C(I2)+SM*UP + DO 110 I=L1,M + C(I4)=C(I4)+SM*U(1,I) + 110 I4=I4+ICE + 120 CONTINUE + 130 RETURN + END + + + +C SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) +C +C Algorithm NNLS: NONNEGATIVE LEAST SQUARES +C +c The original version of this code was developed by +c Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory +c 1973 JUN 15, and published in the book +c "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. +c Revised FEB 1995 to accompany reprinting of the book by SIAM. +c +C GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN +C N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM +C +C A * X = B SUBJECT TO X .GE. 0 +C ------------------------------------------------------------------ +c Subroutine Arguments +c +C A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE +C ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N +C MATRIX, A. ON EXIT A() CONTAINS +C THE PRODUCT MATRIX, Q*A , WHERE Q IS AN +C M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY +C THIS SUBROUTINE. +C B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() CON- +C TAINS Q*B. +C X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL +C CONTAIN THE SOLUTION VECTOR. +C RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE +C RESIDUAL VECTOR. +C W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN +C THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0. +C FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z +C ZZ() AN M-ARRAY OF WORKING SPACE. +C INDEX() AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N. +C ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS +C P AND Z AS FOLLOWS.. +C +C INDEX(1) THRU INDEX(NSETP) = SET P. +C INDEX(IZ1) THRU INDEX(IZ2) = SET Z. +C IZ1 = NSETP + 1 = NPP1 +C IZ2 = N +C MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING +C MEANINGS. +C 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY. +C 2 THE DIMENSIONS OF THE PROBLEM ARE BAD. +C EITHER M .LE. 0 OR N .LE. 0. +C 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS. +C +C ------------------------------------------------------------------ + SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) +C ------------------------------------------------------------------ + integer I, II, IP, ITER, ITMAX, IZ, IZ1, IZ2, IZMAX, J, JJ, JZ, L + integer M, MDA, MODE,N, NPP1, NSETP, RTNKEY +c integer INDEX(N) +c double precision A(MDA,N), B(M), W(N), X(N), ZZ(M) + integer INDEX(*) + double precision A(MDA,*), B(*), W(*), X(*), ZZ(*) + double precision ALPHA, ASAVE, CC, DIFF, DUMMY, FACTOR, RNORM + double precision SM, SS, T, TEMP, TWO, UNORM, UP, WMAX + double precision ZERO, ZTEST + parameter(FACTOR = 0.01d0) + parameter(TWO = 2.0d0, ZERO = 0.0d0) +C ------------------------------------------------------------------ + MODE=1 + IF (M .le. 0 .or. N .le. 0) then + MODE=2 + RETURN + endif + ITER=0 + ITMAX=3*N +C +C INITIALIZE THE ARRAYS INDEX() AND X(). +C + DO 20 I=1,N + X(I)=ZERO + 20 INDEX(I)=I +C + IZ2=N + IZ1=1 + NSETP=0 + NPP1=1 +C ****** MAIN LOOP BEGINS HERE ****** + 30 CONTINUE +C QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION. +C OR IF M COLS OF A HAVE BEEN TRIANGULARIZED. +C + IF (IZ1 .GT.IZ2.OR.NSETP.GE.M) GO TO 350 +C +C COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W(). +C + DO 50 IZ=IZ1,IZ2 + J=INDEX(IZ) + SM=ZERO + DO 40 L=NPP1,M + 40 SM=SM+A(L,J)*B(L) + W(J)=SM + 50 continue +C FIND LARGEST POSITIVE W(J). + 60 continue + WMAX=ZERO + DO 70 IZ=IZ1,IZ2 + J=INDEX(IZ) + IF (W(J) .gt. WMAX) then + WMAX=W(J) + IZMAX=IZ + endif + 70 CONTINUE +C +C IF WMAX .LE. 0. GO TO TERMINATION. +C THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS. +C + IF (WMAX .le. ZERO) go to 350 + IZ=IZMAX + J=INDEX(IZ) +C +C THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P. +C BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID +C NEAR LINEAR DEPENDENCE. +C + ASAVE=A(NPP1,J) + CALL H12 (1,NPP1,NPP1+1,M,A(1,J),1,UP,DUMMY,1,1,0) + UNORM=ZERO + IF (NSETP .ne. 0) then + DO 90 L=1,NSETP + 90 UNORM=UNORM+A(L,J)**2 + endif + UNORM=sqrt(UNORM) + IF (DIFF(UNORM+ABS(A(NPP1,J))*FACTOR,UNORM) .gt. ZERO) then +C +C COL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE ZZ +C AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ). +C + DO 120 L=1,M + 120 ZZ(L)=B(L) + CALL H12 (2,NPP1,NPP1+1,M,A(1,J),1,UP,ZZ,1,1,1) + ZTEST=ZZ(NPP1)/A(NPP1,J) +C +C SEE IF ZTEST IS POSITIVE +C + IF (ZTEST .gt. ZERO) go to 140 + endif +C +C REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P. +C RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL +C COEFFS AGAIN. +C + A(NPP1,J)=ASAVE + W(J)=ZERO + GO TO 60 +C +C THE INDEX J=INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM +C SET Z TO SET P. UPDATE B, UPDATE INDICES, APPLY HOUSEHOLDER +C TRANSFORMATIONS TO COLS IN NEW SET Z, ZERO SUBDIAGONAL ELTS IN +C COL J, SET W(J)=0. +C + 140 continue + DO 150 L=1,M + 150 B(L)=ZZ(L) +C + INDEX(IZ)=INDEX(IZ1) + INDEX(IZ1)=J + IZ1=IZ1+1 + NSETP=NPP1 + NPP1=NPP1+1 +C + IF (IZ1 .le. IZ2) then + DO 160 JZ=IZ1,IZ2 + JJ=INDEX(JZ) + CALL H12 (2,NSETP,NPP1,M,A(1,J),1,UP,A(1,JJ),1,MDA,1) + 160 continue + endif +C + IF (NSETP .ne. M) then + DO 180 L=NPP1,M + 180 A(L,J)=ZERO + endif +C + W(J)=ZERO +C SOLVE THE TRIANGULAR SYSTEM. +C STORE THE SOLUTION TEMPORARILY IN ZZ(). + RTNKEY = 1 + GO TO 400 + 200 CONTINUE +C +C ****** SECONDARY LOOP BEGINS HERE ****** +C +C ITERATION COUNTER. +C + 210 continue + ITER=ITER+1 + IF (ITER .gt. ITMAX) then + MODE=3 + write (*,'(/a)') ' NNLS quitting on iteration count.' + GO TO 350 + endif +C +C SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE. +C IF NOT COMPUTE ALPHA. +C + ALPHA=TWO + DO 240 IP=1,NSETP + L=INDEX(IP) + IF (ZZ(IP) .le. ZERO) then + T=-X(L)/(ZZ(IP)-X(L)) + IF (ALPHA .gt. T) then + ALPHA=T + JJ=IP + endif + endif + 240 CONTINUE +C +C IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL +C STILL = 2. IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP. +C + IF (ALPHA.EQ.TWO) GO TO 330 +C +C OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO +C INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ. +C + DO 250 IP=1,NSETP + L=INDEX(IP) + X(L)=X(L)+ALPHA*(ZZ(IP)-X(L)) + 250 continue +C +C MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I +C FROM SET P TO SET Z. +C + I=INDEX(JJ) + 260 continue + X(I)=ZERO +C + IF (JJ .ne. NSETP) then + JJ=JJ+1 + DO 280 J=JJ,NSETP + II=INDEX(J) + INDEX(J-1)=II + CALL G1 (A(J-1,II),A(J,II),CC,SS,A(J-1,II)) + A(J,II)=ZERO + DO 270 L=1,N + IF (L.NE.II) then +c +c Apply procedure G2 (CC,SS,A(J-1,L),A(J,L)) +c + TEMP = A(J-1,L) + A(J-1,L) = CC*TEMP + SS*A(J,L) + A(J,L) =-SS*TEMP + CC*A(J,L) + endif + 270 CONTINUE +c +c Apply procedure G2 (CC,SS,B(J-1),B(J)) +c + TEMP = B(J-1) + B(J-1) = CC*TEMP + SS*B(J) + B(J) =-SS*TEMP + CC*B(J) + 280 continue + endif +c + NPP1=NSETP + NSETP=NSETP-1 + IZ1=IZ1-1 + INDEX(IZ1)=I +C +C SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE. THEY SHOULD +C BE BECAUSE OF THE WAY ALPHA WAS DETERMINED. +C IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR. ANY +C THAT ARE NONPOSITIVE WILL BE SET TO ZERO +C AND MOVED FROM SET P TO SET Z. +C + DO 300 JJ=1,NSETP + I=INDEX(JJ) + IF (X(I) .le. ZERO) go to 260 + 300 CONTINUE +C +C COPY B( ) INTO ZZ( ). THEN SOLVE AGAIN AND LOOP BACK. +C + DO 310 I=1,M + 310 ZZ(I)=B(I) + RTNKEY = 2 + GO TO 400 + 320 CONTINUE + GO TO 210 +C ****** END OF SECONDARY LOOP ****** +C + 330 continue + DO 340 IP=1,NSETP + I=INDEX(IP) + 340 X(I)=ZZ(IP) +C ALL NEW COEFFS ARE POSITIVE. LOOP BACK TO BEGINNING. + GO TO 30 +C +C ****** END OF MAIN LOOP ****** +C +C COME TO HERE FOR TERMINATION. +C COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR. +C + 350 continue + SM=ZERO + IF (NPP1 .le. M) then + DO 360 I=NPP1,M + 360 SM=SM+B(I)**2 + else + DO 380 J=1,N + 380 W(J)=ZERO + endif + RNORM=sqrt(SM) + RETURN +C +C THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE +C TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ(). +C + 400 continue + DO 430 L=1,NSETP + IP=NSETP+1-L + IF (L .ne. 1) then + DO 410 II=1,IP + ZZ(II)=ZZ(II)-A(II,JJ)*ZZ(IP+1) + 410 continue + endif + JJ=INDEX(IP) + ZZ(IP)=ZZ(IP)/A(IP,JJ) + 430 continue + go to (200, 320), RTNKEY + END +
--- a/nnls.h Mon Oct 18 09:35:26 2010 +0000 +++ b/nnls.h Thu Oct 21 11:58:28 2010 +0100 @@ -1,126 +1,15 @@ #ifndef NNLS_H #define NNLS_H -#include <stdio.h> -#include <math.h> -#define nnls_max(a,b) ((a) >= (b) ? (a) : (b)) -#define nnls_abs(x) ((x) >= 0 ? (x) : -(x)) +extern "C" { -typedef int integer; -typedef float floatreal; +int nnls_(double *a, int *mda, int *m, int *n, + double *b, double *x, double *rnorm, + double *w, double *zz, int *index, int *mode); -/* SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) */ - -/* Algorithm NNLS: NONNEGATIVE LEAST SQUARES */ - -/* The original version of this code was developed by */ -/* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */ -/* 1973 JUN 15, and published in the book */ -/* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */ -/* Revised FEB 1995 to accompany reprinting of the book by SIAM. */ - -/* GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN */ -/* N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM */ - -/* A * X = B SUBJECT TO X .GE. 0 */ -/* ------------------------------------------------------------------ */ -/* Subroutine Arguments */ - -/* A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE */ -/* ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N */ -/* MATRIX, A. ON EXIT A() CONTAINS */ -/* THE PRODUCT MATRIX, Q*A , WHERE Q IS AN */ -/* M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY */ -/* THIS SUBROUTINE. */ -/* B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() CON- */ -/* TAINS Q*B. */ -/* X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL */ -/* CONTAIN THE SOLUTION VECTOR. */ -/* RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE */ -/* RESIDUAL VECTOR. */ -/* W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN */ -/* THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0. */ -/* FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z */ -/* ZZ() AN M-ARRAY OF WORKING SPACE. */ -/* INDEX() AN INT WORKING ARRAY OF LENGTH AT LEAST N. */ -/* ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS */ -/* P AND Z AS FOLLOWS.. */ - -/* INDEX(1) THRU INDEX(NSETP) = SET P. */ -/* INDEX(IZ1) THRU INDEX(IZ2) = SET Z. */ -/* IZ1 = NSETP + 1 = NPP1 */ -/* IZ2 = N */ -/* MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING */ -/* MEANINGS. */ -/* 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY. */ -/* 2 THE DIMENSIONS OF THE PROBLEM ARE BAD. */ -/* EITHER M .LE. 0 OR N .LE. 0. */ -/* 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS. */ - -/* ------------------------------------------------------------------ */ -/* Subroutine */ -int nnls(float* a, int mda, int m, int n, - float* b, float* x, float* rnorm, - float* w, float* zz, int* index, int* mode); +#define NNLS nnls_ +} - -/* SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) */ - -/* CONSTRUCTION AND/OR APPLICATION OF A SINGLE */ -/* HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B */ - -/* The original version of this code was developed by */ -/* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */ -/* 1973 JUN 12, and published in the book */ -/* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */ -/* Revised FEB 1995 to accompany reprinting of the book by SIAM. */ -/* ------------------------------------------------------------------ */ -/* Subroutine Arguments */ - -/* MODE = 1 OR 2 Selects Algorithm H1 to construct and apply a */ -/* Householder transformation, or Algorithm H2 to apply a */ -/* previously constructed transformation. */ -/* LPIVOT IS THE INDEX OF THE PIVOT ELEMENT. */ -/* L1,M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO */ -/* ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M */ -/* THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION. */ -/* U(),IUE,UP On entry with MODE = 1, U() contains the pivot */ -/* vector. IUE is the storage increment between elements. */ -/* On exit when MODE = 1, U() and UP contain quantities */ -/* defining the vector U of the Householder transformation. */ -/* on entry with MODE = 2, U() and UP should contain */ -/* quantities previously computed with MODE = 1. These will */ -/* not be modified during the entry with MODE = 2. */ -/* C() ON ENTRY with MODE = 1 or 2, C() CONTAINS A MATRIX WHICH */ -/* WILL BE REGARDED AS A SET OF VECTORS TO WHICH THE */ -/* HOUSEHOLDER TRANSFORMATION IS TO BE APPLIED. */ -/* ON EXIT C() CONTAINS THE SET OF TRANSFORMED VECTORS. */ -/* ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C(). */ -/* ICV STORAGE INCREMENT BETWEEN VECTORS IN C(). */ -/* NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0 */ -/* NO OPERATIONS WILL BE DONE ON C(). */ -/* ------------------------------------------------------------------ */ -/* Subroutine */ -int h12(int mode, int* lpivot, int* l1, - int m, float* u, int* iue, float* up, float* c__, - int* ice, int* icv, int* ncv); - - - /* COMPUTE ORTHOGONAL ROTATION MATRIX.. */ - - /* The original version of this code was developed by */ - /* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory - */ - /* 1973 JUN 12, and published in the book */ - /* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */ - /* Revised FEB 1995 to accompany reprinting of the book by SIAM. */ - - /* COMPUTE.. MATRIX (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2)) */ - /* (-S,C) (-S,C)(B) ( 0 ) */ - /* COMPUTE SIG = SQRT(A**2+B**2) */ - /* SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT */ - /* SIG MAY BE IN THE SAME LOCATION AS A OR B . */ -int g1(float* a, float* b, float* cterm, float* sterm, float* sig); #endif