1
// -*- mode: C++; tab-width: 2; -*-
4
// --------------------------------------------------------------------------
5
// OpenMS Mass Spectrometry Framework
6
// --------------------------------------------------------------------------
7
// Copyright (C) 2003-2011 -- Oliver Kohlbacher, Knut Reinert
9
// This library is free software; you can redistribute it and/or
10
// modify it under the terms of the GNU Lesser General Public
11
// License as published by the Free Software Foundation; either
12
// version 2.1 of the License, or (at your option) any later version.
14
// This library is distributed in the hope that it will be useful,
15
// but WITHOUT ANY WARRANTY; without even the implied warranty of
16
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17
// Lesser General Public License for more details.
19
// You should have received a copy of the GNU Lesser General Public
20
// License along with this library; if not, write to the Free Software
21
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23
// --------------------------------------------------------------------------
24
// $Maintainer: Chris Bielow $
25
// $Authors: Chris Bielow $
26
// --------------------------------------------------------------------------
31
#include <OpenMS/MATH/MISC/NNLS/NNLS.h>
35
The code below was converted from FORTRAN using f2c from http://www.netlib.org/lawson-hanson/all
36
Some modifications were made, in order for it to run properly (search for "--removed", "-- added" and "--changed" in the code below)
46
/* start of original code (with modification as described above) */
48
/* nnls.F -- translated by f2c (version 20100827).
49
You must link the resulting object file with libf2c:
50
on Microsoft Windows system, link with libf2c.lib;
51
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
52
or, if you install libf2c.a in a standard place, with -lf2c -lm
53
-- in that order, at the end of the command line, as in
55
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
57
http://www.netlib.org/f2c/libf2c.zip
60
/* #include "f2c.h" -- removed */
62
/* Table of constant values */
64
static integer c__1 = 1;
65
static integer c__0 = 0;
66
static integer c__2 = 2;
68
/* SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) */
70
/* Algorithm NNLS: NONNEGATIVE LEAST SQUARES */
72
/* The original version of this code was developed by */
73
/* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
74
/* 1973 JUN 15, and published in the book */
75
/* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
76
/* Revised FEB 1995 to accompany reprinting of the book by SIAM. */
78
/* GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN */
79
/* N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM */
81
/* A * X = B SUBJECT TO X .GE. 0 */
82
/* ------------------------------------------------------------------ */
83
/* Subroutine Arguments */
85
/* A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE */
86
/* ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N */
87
/* MATRIX, A. ON EXIT A() CONTAINS */
88
/* THE PRODUCT MATRIX, Q*A , WHERE Q IS AN */
89
/* M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY */
90
/* THIS SUBROUTINE. */
91
/* B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() CON- */
93
/* X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL */
94
/* CONTAIN THE SOLUTION VECTOR. */
95
/* RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE */
96
/* RESIDUAL VECTOR. */
97
/* W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN */
98
/* THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0. */
99
/* FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z */
100
/* ZZ() AN M-ARRAY OF WORKING SPACE. */
101
/* INDEX() AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N. */
102
/* ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS */
103
/* P AND Z AS FOLLOWS.. */
105
/* INDEX(1) THRU INDEX(NSETP) = SET P. */
106
/* INDEX(IZ1) THRU INDEX(IZ2) = SET Z. */
107
/* IZ1 = NSETP + 1 = NPP1 */
109
/* MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING */
111
/* 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY. */
112
/* 2 THE DIMENSIONS OF THE PROBLEM ARE BAD. */
113
/* EITHER M .LE. 0 OR N .LE. 0. */
114
/* 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS. */
116
/* ------------------------------------------------------------------ */
117
/* Subroutine */ int nnls_(doublereal *a, integer *mda, integer *m, integer *
118
n, doublereal *b, doublereal *x, doublereal *rnorm, doublereal *w,
119
doublereal *zz, integer *index, integer *mode)
121
/* System generated locals */
122
integer a_dim1, a_offset, i__1, i__2;
123
doublereal d__1, d__2;
125
/* Builtin functions */
126
/* double sqrt(doublereal); --removed */
127
/* integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); -- removed */
129
/* Local variables */
130
static integer i__, j, l;
132
/* Subroutine */ int g1_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *);
133
static doublereal cc;
134
/* Subroutine */ int h12_(integer *, integer *, integer *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *, integer *, integer *);
135
static integer ii, jj, ip;
136
static doublereal sm;
137
static integer iz, jz;
138
static doublereal up, ss;
139
static integer iz1, iz2, npp1;
140
doublereal diff_(doublereal *, doublereal *);
142
static doublereal temp, wmax, alpha, asave;
143
static integer itmax, izmax, nsetp;
144
static doublereal dummy, unorm, ztest;
145
static integer rtnkey;
147
/* Fortran I/O blocks */
148
/* static cilist io___22 = { 0, 6, 0, "(/a)", 0 }; --removed */
151
/* ------------------------------------------------------------------ */
152
/* integer INDEX(N) */
153
/* double precision A(MDA,N), B(M), W(N), X(N), ZZ(M) */
154
/* ------------------------------------------------------------------ */
155
/* Parameter adjustments */
157
a_offset = 1 + a_dim1;
167
if (*m <= 0 || *n <= 0) {
174
/* INITIALIZE THE ARRAYS INDEX() AND X(). */
177
for (i__ = 1; i__ <= i__1; ++i__) {
187
/* ****** MAIN LOOP BEGINS HERE ****** */
189
/* QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION. */
190
/* OR IF M COLS OF A HAVE BEEN TRIANGULARIZED. */
192
if (iz1 > iz2 || nsetp >= *m) {
196
/* COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W(). */
199
for (iz = iz1; iz <= i__1; ++iz) {
203
for (l = npp1; l <= i__2; ++l) {
205
sm += a[l + j * a_dim1] * b[l];
210
/* FIND LARGEST POSITIVE W(J). */
214
for (iz = iz1; iz <= i__1; ++iz) {
223
/* IF WMAX .LE. 0. GO TO TERMINATION. */
224
/* THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS. */
232
/* THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P. */
233
/* BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID */
234
/* NEAR LINEAR DEPENDENCE. */
236
asave = a[npp1 + j * a_dim1];
238
h12_(&c__1, &npp1, &i__1, m, &a[j * a_dim1 + 1], &c__1, &up, &dummy, &
243
for (l = 1; l <= i__1; ++l) {
245
/* Computing 2nd power */
246
d__1 = a[l + j * a_dim1];
247
unorm += d__1 * d__1;
251
d__2 = unorm + (d__1 = a[npp1 + j * a_dim1], fabs(d__1)) * .01; /* --changed */
252
if (diff_(&d__2, &unorm) > 0.) {
254
/* COL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE ZZ */
255
/* AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ). */
258
for (l = 1; l <= i__1; ++l) {
263
h12_(&c__2, &npp1, &i__1, m, &a[j * a_dim1 + 1], &c__1, &up, &zz[1], &
265
ztest = zz[npp1] / a[npp1 + j * a_dim1];
267
/* SEE IF ZTEST IS POSITIVE */
274
/* REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P. */
275
/* RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL */
278
a[npp1 + j * a_dim1] = asave;
282
/* THE INDEX J=INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM */
283
/* SET Z TO SET P. UPDATE B, UPDATE INDICES, APPLY HOUSEHOLDER */
284
/* TRANSFORMATIONS TO COLS IN NEW SET Z, ZERO SUBDIAGONAL ELTS IN */
285
/* COL J, SET W(J)=0. */
289
for (l = 1; l <= i__1; ++l) {
294
index[iz] = index[iz1];
302
for (jz = iz1; jz <= i__1; ++jz) {
304
h12_(&c__2, &nsetp, &npp1, m, &a[j * a_dim1 + 1], &c__1, &up, &a[
305
jj * a_dim1 + 1], &c__1, mda, &c__1);
312
for (l = npp1; l <= i__1; ++l) {
314
a[l + j * a_dim1] = 0.;
319
/* SOLVE THE TRIANGULAR SYSTEM. */
320
/* STORE THE SOLUTION TEMPORARILY IN ZZ(). */
325
/* ****** SECONDARY LOOP BEGINS HERE ****** */
327
/* ITERATION COUNTER. */
334
do_fio(&c__1, " NNLS quitting on iteration count.", (ftnlen)34);
335
e_wsfe(); --removed */
339
/* SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE. */
340
/* IF NOT COMPUTE ALPHA. */
344
for (ip = 1; ip <= i__1; ++ip) {
347
t = -x[l] / (zz[ip] - x[l]);
356
/* IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL */
357
/* STILL = 2. IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP. */
363
/* OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO */
364
/* INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ. */
367
for (ip = 1; ip <= i__1; ++ip) {
369
x[l] += alpha * (zz[ip] - x[l]);
373
/* MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I */
374
/* FROM SET P TO SET Z. */
383
for (j = jj; j <= i__1; ++j) {
386
g1_(&a[j - 1 + ii * a_dim1], &a[j + ii * a_dim1], &cc, &ss, &a[j
388
a[j + ii * a_dim1] = 0.;
390
for (l = 1; l <= i__2; ++l) {
393
/* Apply procedure G2 (CC,SS,A(J-1,L),A(J,L)) */
395
temp = a[j - 1 + l * a_dim1];
396
a[j - 1 + l * a_dim1] = cc * temp + ss * a[j + l * a_dim1]
398
a[j + l * a_dim1] = -ss * temp + cc * a[j + l * a_dim1];
403
/* Apply procedure G2 (CC,SS,B(J-1),B(J)) */
406
b[j - 1] = cc * temp + ss * b[j];
407
b[j] = -ss * temp + cc * b[j];
417
/* SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE. THEY SHOULD */
418
/* BE BECAUSE OF THE WAY ALPHA WAS DETERMINED. */
419
/* IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR. ANY */
420
/* THAT ARE NONPOSITIVE WILL BE SET TO ZERO */
421
/* AND MOVED FROM SET P TO SET Z. */
424
for (jj = 1; jj <= i__1; ++jj) {
432
/* COPY B( ) INTO ZZ( ). THEN SOLVE AGAIN AND LOOP BACK. */
435
for (i__ = 1; i__ <= i__1; ++i__) {
443
/* ****** END OF SECONDARY LOOP ****** */
447
for (ip = 1; ip <= i__1; ++ip) {
452
/* ALL NEW COEFFS ARE POSITIVE. LOOP BACK TO BEGINNING. */
455
/* ****** END OF MAIN LOOP ****** */
457
/* COME TO HERE FOR TERMINATION. */
458
/* COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR. */
464
for (i__ = npp1; i__ <= i__1; ++i__) {
466
/* Computing 2nd power */
472
for (j = 1; j <= i__1; ++j) {
480
/* THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE */
481
/* TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ(). */
485
for (l = 1; l <= i__1; ++l) {
489
for (ii = 1; ii <= i__2; ++ii) {
490
zz[ii] -= a[ii + jj * a_dim1] * zz[ip + 1];
495
zz[ip] /= a[ip + jj * a_dim1];
505
/* Subroutine */ int g1_(doublereal *a, doublereal *b, doublereal *cterm,
506
doublereal *sterm, doublereal *sig)
508
/* System generated locals */
511
/* Builtin functions */
512
/* double sqrt(doublereal), d_sign(doublereal *, doublereal *); --removed */
514
/* Local variables */
515
static doublereal xr, yr;
518
/* COMPUTE ORTHOGONAL ROTATION MATRIX.. */
520
/* The original version of this code was developed by */
521
/* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
522
/* 1973 JUN 12, and published in the book */
523
/* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
524
/* Revised FEB 1995 to accompany reprinting of the book by SIAM. */
526
/* COMPUTE.. MATRIX (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2)) */
527
/* (-S,C) (-S,C)(B) ( 0 ) */
528
/* COMPUTE SIG = SQRT(A**2+B**2) */
529
/* SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT */
530
/* SIG MAY BE IN THE SAME LOCATION AS A OR B . */
531
/* ------------------------------------------------------------------ */
532
/* ------------------------------------------------------------------ */
533
if (fabs(*a) > fabs(*b)) {
535
/* Computing 2nd power */
537
yr = sqrt(d__1 * d__1 + 1.);
539
*cterm = d_sign_(d__1, *a); /* --changed */
540
*sterm = *cterm * xr;
541
*sig = fabs(*a) * yr;
546
/* Computing 2nd power */
548
yr = sqrt(d__1 * d__1 + 1.);
550
*sterm = d_sign_(d__1, *b); /* --changed */
551
*cterm = *sterm * xr;
552
*sig = fabs(*b) * yr;
561
/* SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) */
563
/* CONSTRUCTION AND/OR APPLICATION OF A SINGLE */
564
/* HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B */
566
/* The original version of this code was developed by */
567
/* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
568
/* 1973 JUN 12, and published in the book */
569
/* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
570
/* Revised FEB 1995 to accompany reprinting of the book by SIAM. */
571
/* ------------------------------------------------------------------ */
572
/* Subroutine Arguments */
574
/* MODE = 1 OR 2 Selects Algorithm H1 to construct and apply a */
575
/* Householder transformation, or Algorithm H2 to apply a */
576
/* previously constructed transformation. */
577
/* LPIVOT IS THE INDEX OF THE PIVOT ELEMENT. */
578
/* L1,M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO */
579
/* ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M */
580
/* THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION. */
581
/* U(),IUE,UP On entry with MODE = 1, U() contains the pivot */
582
/* vector. IUE is the storage increment between elements. */
583
/* On exit when MODE = 1, U() and UP contain quantities */
584
/* defining the vector U of the Householder transformation. */
585
/* on entry with MODE = 2, U() and UP should contain */
586
/* quantities previously computed with MODE = 1. These will */
587
/* not be modified during the entry with MODE = 2. */
588
/* C() ON ENTRY with MODE = 1 or 2, C() CONTAINS A MATRIX WHICH */
589
/* WILL BE REGARDED AS A SET OF VECTORS TO WHICH THE */
590
/* HOUSEHOLDER TRANSFORMATION IS TO BE APPLIED. */
591
/* ON EXIT C() CONTAINS THE SET OF TRANSFORMED VECTORS. */
592
/* ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C(). */
593
/* ICV STORAGE INCREMENT BETWEEN VECTORS IN C(). */
594
/* NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0 */
595
/* NO OPERATIONS WILL BE DONE ON C(). */
596
/* ------------------------------------------------------------------ */
597
/* Subroutine */ int h12_(integer *mode, integer *lpivot, integer *l1,
598
integer *m, doublereal *u, integer *iue, doublereal *up, doublereal *
599
c__, integer *ice, integer *icv, integer *ncv)
601
/* System generated locals */
602
integer u_dim1, u_offset, i__1, i__2;
603
doublereal d__1, d__2;
605
/* Builtin functions */
606
/* double sqrt(doublereal); --removed */
608
/* Local variables */
610
static integer i__, j, i2, i3, i4;
611
static doublereal cl, sm;
613
static doublereal clinv;
615
/* ------------------------------------------------------------------ */
616
/* double precision U(IUE,M) */
617
/* ------------------------------------------------------------------ */
618
/* Parameter adjustments */
620
u_offset = 1 + u_dim1;
625
if (0 >= *lpivot || *lpivot >= *l1 || *l1 > *m) {
628
cl = (d__1 = u[*lpivot * u_dim1 + 1], fabs(d__1));
632
/* ****** CONSTRUCT THE TRANSFORMATION. ****** */
634
for (j = *l1; j <= i__1; ++j) {
637
d__2 = (d__1 = u[j * u_dim1 + 1], fabs(d__1));
638
cl = std::max(d__2,cl); /* --changed */
647
/* Computing 2nd power */
648
d__1 = u[*lpivot * u_dim1 + 1] * clinv;
651
for (j = *l1; j <= i__1; ++j) {
653
/* Computing 2nd power */
654
d__1 = u[j * u_dim1 + 1] * clinv;
658
if (u[*lpivot * u_dim1 + 1] <= 0.) {
666
*up = u[*lpivot * u_dim1 + 1] - cl;
667
u[*lpivot * u_dim1 + 1] = cl;
669
/* ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C. ****** */
681
b = *up * u[*lpivot * u_dim1 + 1];
682
/* B MUST BE NONPOSITIVE HERE. IF B = 0., RETURN. */
691
i2 = 1 - *icv + *ice * (*lpivot - 1);
692
incr = *ice * (*l1 - *lpivot);
694
for (j = 1; j <= i__1; ++j) {
700
for (i__ = *l1; i__ <= i__2; ++i__) {
701
sm += c__[i3] * u[i__ * u_dim1 + 1];
714
for (i__ = *l1; i__ <= i__2; ++i__) {
715
c__[i4] += sm * u[i__ * u_dim1 + 1];
726
doublereal diff_(doublereal *x, doublereal *y)
728
/* System generated locals */
732
/* Function used in tests that depend on machine precision. */
734
/* The original version of this code was developed by */
735
/* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
736
/* 1973 JUN 7, and published in the book */
737
/* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
738
/* Revised FEB 1995 to accompany reprinting of the book by SIAM. */
744
/* -- added manually */
745
double d_sign_(double& a, double& b)
747
double x = (a >= 0 ? a : - a);
748
return (b >= 0 ? x : -x);
752
} // namespace OpenMS
added
removed
source/MATH/MISC/NNLS/NNLS.C
