diff nnls.f @ 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)

--- /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   
+