matthiasm@6
|
1 #ifndef NNLS_H
|
matthiasm@6
|
2 #define NNLS_H
|
matthiasm@6
|
3
|
matthiasm@6
|
4 #include <stdio.h>
|
matthiasm@6
|
5 #include <math.h>
|
matthiasm@6
|
6 #define nnls_max(a,b) ((a) >= (b) ? (a) : (b))
|
matthiasm@6
|
7 #define nnls_abs(x) ((x) >= 0 ? (x) : -(x))
|
matthiasm@6
|
8
|
matthiasm@6
|
9 typedef int integer;
|
matthiasm@6
|
10 typedef float floatreal;
|
matthiasm@6
|
11
|
matthiasm@6
|
12 /* SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) */
|
matthiasm@6
|
13
|
matthiasm@6
|
14 /* Algorithm NNLS: NONNEGATIVE LEAST SQUARES */
|
matthiasm@6
|
15
|
matthiasm@6
|
16 /* The original version of this code was developed by */
|
matthiasm@6
|
17 /* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
|
matthiasm@6
|
18 /* 1973 JUN 15, and published in the book */
|
matthiasm@6
|
19 /* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
|
matthiasm@6
|
20 /* Revised FEB 1995 to accompany reprinting of the book by SIAM. */
|
matthiasm@6
|
21
|
matthiasm@6
|
22 /* GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN */
|
matthiasm@6
|
23 /* N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM */
|
matthiasm@6
|
24
|
matthiasm@6
|
25 /* A * X = B SUBJECT TO X .GE. 0 */
|
matthiasm@6
|
26 /* ------------------------------------------------------------------ */
|
matthiasm@6
|
27 /* Subroutine Arguments */
|
matthiasm@6
|
28
|
matthiasm@6
|
29 /* A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE */
|
matthiasm@6
|
30 /* ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N */
|
matthiasm@6
|
31 /* MATRIX, A. ON EXIT A() CONTAINS */
|
matthiasm@6
|
32 /* THE PRODUCT MATRIX, Q*A , WHERE Q IS AN */
|
matthiasm@6
|
33 /* M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY */
|
matthiasm@6
|
34 /* THIS SUBROUTINE. */
|
matthiasm@6
|
35 /* B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() CON- */
|
matthiasm@6
|
36 /* TAINS Q*B. */
|
matthiasm@6
|
37 /* X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL */
|
matthiasm@6
|
38 /* CONTAIN THE SOLUTION VECTOR. */
|
matthiasm@6
|
39 /* RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE */
|
matthiasm@6
|
40 /* RESIDUAL VECTOR. */
|
matthiasm@6
|
41 /* W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN */
|
matthiasm@6
|
42 /* THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0. */
|
matthiasm@6
|
43 /* FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z */
|
matthiasm@6
|
44 /* ZZ() AN M-ARRAY OF WORKING SPACE. */
|
matthiasm@6
|
45 /* INDEX() AN INT WORKING ARRAY OF LENGTH AT LEAST N. */
|
matthiasm@6
|
46 /* ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS */
|
matthiasm@6
|
47 /* P AND Z AS FOLLOWS.. */
|
matthiasm@6
|
48
|
matthiasm@6
|
49 /* INDEX(1) THRU INDEX(NSETP) = SET P. */
|
matthiasm@6
|
50 /* INDEX(IZ1) THRU INDEX(IZ2) = SET Z. */
|
matthiasm@6
|
51 /* IZ1 = NSETP + 1 = NPP1 */
|
matthiasm@6
|
52 /* IZ2 = N */
|
matthiasm@6
|
53 /* MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING */
|
matthiasm@6
|
54 /* MEANINGS. */
|
matthiasm@6
|
55 /* 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY. */
|
matthiasm@6
|
56 /* 2 THE DIMENSIONS OF THE PROBLEM ARE BAD. */
|
matthiasm@6
|
57 /* EITHER M .LE. 0 OR N .LE. 0. */
|
matthiasm@6
|
58 /* 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS. */
|
matthiasm@6
|
59
|
matthiasm@6
|
60 /* ------------------------------------------------------------------ */
|
matthiasm@6
|
61 /* Subroutine */
|
matthiasm@6
|
62 int nnls(float* a, int mda, int m, int n,
|
matthiasm@6
|
63 float* b, float* x, float* rnorm,
|
matthiasm@6
|
64 float* w, float* zz, int* index, int* mode);
|
matthiasm@6
|
65
|
matthiasm@6
|
66
|
matthiasm@6
|
67
|
matthiasm@6
|
68 /* SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) */
|
matthiasm@6
|
69
|
matthiasm@6
|
70 /* CONSTRUCTION AND/OR APPLICATION OF A SINGLE */
|
matthiasm@6
|
71 /* HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B */
|
matthiasm@6
|
72
|
matthiasm@6
|
73 /* The original version of this code was developed by */
|
matthiasm@6
|
74 /* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */
|
matthiasm@6
|
75 /* 1973 JUN 12, and published in the book */
|
matthiasm@6
|
76 /* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
|
matthiasm@6
|
77 /* Revised FEB 1995 to accompany reprinting of the book by SIAM. */
|
matthiasm@6
|
78 /* ------------------------------------------------------------------ */
|
matthiasm@6
|
79 /* Subroutine Arguments */
|
matthiasm@6
|
80
|
matthiasm@6
|
81 /* MODE = 1 OR 2 Selects Algorithm H1 to construct and apply a */
|
matthiasm@6
|
82 /* Householder transformation, or Algorithm H2 to apply a */
|
matthiasm@6
|
83 /* previously constructed transformation. */
|
matthiasm@6
|
84 /* LPIVOT IS THE INDEX OF THE PIVOT ELEMENT. */
|
matthiasm@6
|
85 /* L1,M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO */
|
matthiasm@6
|
86 /* ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M */
|
matthiasm@6
|
87 /* THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION. */
|
matthiasm@6
|
88 /* U(),IUE,UP On entry with MODE = 1, U() contains the pivot */
|
matthiasm@6
|
89 /* vector. IUE is the storage increment between elements. */
|
matthiasm@6
|
90 /* On exit when MODE = 1, U() and UP contain quantities */
|
matthiasm@6
|
91 /* defining the vector U of the Householder transformation. */
|
matthiasm@6
|
92 /* on entry with MODE = 2, U() and UP should contain */
|
matthiasm@6
|
93 /* quantities previously computed with MODE = 1. These will */
|
matthiasm@6
|
94 /* not be modified during the entry with MODE = 2. */
|
matthiasm@6
|
95 /* C() ON ENTRY with MODE = 1 or 2, C() CONTAINS A MATRIX WHICH */
|
matthiasm@6
|
96 /* WILL BE REGARDED AS A SET OF VECTORS TO WHICH THE */
|
matthiasm@6
|
97 /* HOUSEHOLDER TRANSFORMATION IS TO BE APPLIED. */
|
matthiasm@6
|
98 /* ON EXIT C() CONTAINS THE SET OF TRANSFORMED VECTORS. */
|
matthiasm@6
|
99 /* ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C(). */
|
matthiasm@6
|
100 /* ICV STORAGE INCREMENT BETWEEN VECTORS IN C(). */
|
matthiasm@6
|
101 /* NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0 */
|
matthiasm@6
|
102 /* NO OPERATIONS WILL BE DONE ON C(). */
|
matthiasm@6
|
103 /* ------------------------------------------------------------------ */
|
matthiasm@6
|
104 /* Subroutine */
|
matthiasm@6
|
105 int h12(int mode, int* lpivot, int* l1,
|
matthiasm@6
|
106 int m, float* u, int* iue, float* up, float* c__,
|
matthiasm@6
|
107 int* ice, int* icv, int* ncv);
|
matthiasm@6
|
108
|
matthiasm@6
|
109
|
matthiasm@6
|
110 /* COMPUTE ORTHOGONAL ROTATION MATRIX.. */
|
matthiasm@6
|
111
|
matthiasm@6
|
112 /* The original version of this code was developed by */
|
matthiasm@6
|
113 /* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
|
matthiasm@6
|
114 */
|
matthiasm@6
|
115 /* 1973 JUN 12, and published in the book */
|
matthiasm@6
|
116 /* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */
|
matthiasm@6
|
117 /* Revised FEB 1995 to accompany reprinting of the book by SIAM. */
|
matthiasm@6
|
118
|
matthiasm@6
|
119 /* COMPUTE.. MATRIX (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2)) */
|
matthiasm@6
|
120 /* (-S,C) (-S,C)(B) ( 0 ) */
|
matthiasm@6
|
121 /* COMPUTE SIG = SQRT(A**2+B**2) */
|
matthiasm@6
|
122 /* SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT */
|
matthiasm@6
|
123 /* SIG MAY BE IN THE SAME LOCATION AS A OR B . */
|
matthiasm@6
|
124 int g1(float* a, float* b, float* cterm, float* sterm, float* sig);
|
matthiasm@6
|
125 #endif
|
matthiasm@6
|
126
|