matthiasm@6: #ifndef NNLS_H matthiasm@6: #define NNLS_H matthiasm@6: matthiasm@6: #include matthiasm@6: #include matthiasm@6: #define nnls_max(a,b) ((a) >= (b) ? (a) : (b)) matthiasm@6: #define nnls_abs(x) ((x) >= 0 ? (x) : -(x)) matthiasm@6: matthiasm@6: typedef int integer; matthiasm@6: typedef float floatreal; matthiasm@6: matthiasm@6: /* SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) */ matthiasm@6: matthiasm@6: /* Algorithm NNLS: NONNEGATIVE LEAST SQUARES */ matthiasm@6: matthiasm@6: /* The original version of this code was developed by */ matthiasm@6: /* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */ matthiasm@6: /* 1973 JUN 15, and published in the book */ matthiasm@6: /* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */ matthiasm@6: /* Revised FEB 1995 to accompany reprinting of the book by SIAM. */ matthiasm@6: matthiasm@6: /* GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN */ matthiasm@6: /* N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM */ matthiasm@6: matthiasm@6: /* A * X = B SUBJECT TO X .GE. 0 */ matthiasm@6: /* ------------------------------------------------------------------ */ matthiasm@6: /* Subroutine Arguments */ matthiasm@6: matthiasm@6: /* A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE */ matthiasm@6: /* ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N */ matthiasm@6: /* MATRIX, A. ON EXIT A() CONTAINS */ matthiasm@6: /* THE PRODUCT MATRIX, Q*A , WHERE Q IS AN */ matthiasm@6: /* M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY */ matthiasm@6: /* THIS SUBROUTINE. */ matthiasm@6: /* B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() CON- */ matthiasm@6: /* TAINS Q*B. */ matthiasm@6: /* X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL */ matthiasm@6: /* CONTAIN THE SOLUTION VECTOR. */ matthiasm@6: /* RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE */ matthiasm@6: /* RESIDUAL VECTOR. */ matthiasm@6: /* W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN */ matthiasm@6: /* THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0. */ matthiasm@6: /* FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z */ matthiasm@6: /* ZZ() AN M-ARRAY OF WORKING SPACE. */ matthiasm@6: /* INDEX() AN INT WORKING ARRAY OF LENGTH AT LEAST N. */ matthiasm@6: /* ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS */ matthiasm@6: /* P AND Z AS FOLLOWS.. */ matthiasm@6: matthiasm@6: /* INDEX(1) THRU INDEX(NSETP) = SET P. */ matthiasm@6: /* INDEX(IZ1) THRU INDEX(IZ2) = SET Z. */ matthiasm@6: /* IZ1 = NSETP + 1 = NPP1 */ matthiasm@6: /* IZ2 = N */ matthiasm@6: /* MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING */ matthiasm@6: /* MEANINGS. */ matthiasm@6: /* 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY. */ matthiasm@6: /* 2 THE DIMENSIONS OF THE PROBLEM ARE BAD. */ matthiasm@6: /* EITHER M .LE. 0 OR N .LE. 0. */ matthiasm@6: /* 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS. */ matthiasm@6: matthiasm@6: /* ------------------------------------------------------------------ */ matthiasm@6: /* Subroutine */ matthiasm@6: int nnls(float* a, int mda, int m, int n, matthiasm@6: float* b, float* x, float* rnorm, matthiasm@6: float* w, float* zz, int* index, int* mode); matthiasm@6: matthiasm@6: matthiasm@6: matthiasm@6: /* SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) */ matthiasm@6: matthiasm@6: /* CONSTRUCTION AND/OR APPLICATION OF A SINGLE */ matthiasm@6: /* HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B */ matthiasm@6: matthiasm@6: /* The original version of this code was developed by */ matthiasm@6: /* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */ matthiasm@6: /* 1973 JUN 12, and published in the book */ matthiasm@6: /* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */ matthiasm@6: /* Revised FEB 1995 to accompany reprinting of the book by SIAM. */ matthiasm@6: /* ------------------------------------------------------------------ */ matthiasm@6: /* Subroutine Arguments */ matthiasm@6: matthiasm@6: /* MODE = 1 OR 2 Selects Algorithm H1 to construct and apply a */ matthiasm@6: /* Householder transformation, or Algorithm H2 to apply a */ matthiasm@6: /* previously constructed transformation. */ matthiasm@6: /* LPIVOT IS THE INDEX OF THE PIVOT ELEMENT. */ matthiasm@6: /* L1,M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO */ matthiasm@6: /* ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M */ matthiasm@6: /* THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION. */ matthiasm@6: /* U(),IUE,UP On entry with MODE = 1, U() contains the pivot */ matthiasm@6: /* vector. IUE is the storage increment between elements. */ matthiasm@6: /* On exit when MODE = 1, U() and UP contain quantities */ matthiasm@6: /* defining the vector U of the Householder transformation. */ matthiasm@6: /* on entry with MODE = 2, U() and UP should contain */ matthiasm@6: /* quantities previously computed with MODE = 1. These will */ matthiasm@6: /* not be modified during the entry with MODE = 2. */ matthiasm@6: /* C() ON ENTRY with MODE = 1 or 2, C() CONTAINS A MATRIX WHICH */ matthiasm@6: /* WILL BE REGARDED AS A SET OF VECTORS TO WHICH THE */ matthiasm@6: /* HOUSEHOLDER TRANSFORMATION IS TO BE APPLIED. */ matthiasm@6: /* ON EXIT C() CONTAINS THE SET OF TRANSFORMED VECTORS. */ matthiasm@6: /* ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C(). */ matthiasm@6: /* ICV STORAGE INCREMENT BETWEEN VECTORS IN C(). */ matthiasm@6: /* NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0 */ matthiasm@6: /* NO OPERATIONS WILL BE DONE ON C(). */ matthiasm@6: /* ------------------------------------------------------------------ */ matthiasm@6: /* Subroutine */ matthiasm@6: int h12(int mode, int* lpivot, int* l1, matthiasm@6: int m, float* u, int* iue, float* up, float* c__, matthiasm@6: int* ice, int* icv, int* ncv); matthiasm@6: matthiasm@6: matthiasm@6: /* COMPUTE ORTHOGONAL ROTATION MATRIX.. */ matthiasm@6: matthiasm@6: /* The original version of this code was developed by */ matthiasm@6: /* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory matthiasm@6: */ matthiasm@6: /* 1973 JUN 12, and published in the book */ matthiasm@6: /* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */ matthiasm@6: /* Revised FEB 1995 to accompany reprinting of the book by SIAM. */ matthiasm@6: matthiasm@6: /* COMPUTE.. MATRIX (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2)) */ matthiasm@6: /* (-S,C) (-S,C)(B) ( 0 ) */ matthiasm@6: /* COMPUTE SIG = SQRT(A**2+B**2) */ matthiasm@6: /* SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT */ matthiasm@6: /* SIG MAY BE IN THE SAME LOCATION AS A OR B . */ matthiasm@6: int g1(float* a, float* b, float* cterm, float* sterm, float* sig); matthiasm@6: #endif matthiasm@6: