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