[Scipy-svn] r4757 - in trunk/scipy/optimize: . nnls

scipy-svn@scip... scipy-svn@scip...
Tue Sep 30 02:00:13 CDT 2008


Author: uwe.schmitt
Date: 2008-09-30 02:00:07 -0500 (Tue, 30 Sep 2008)
New Revision: 4757

Added:
   trunk/scipy/optimize/nnls/nnls.f
Removed:
   trunk/scipy/optimize/nnls/NNLS.F
Modified:
   trunk/scipy/optimize/SConscript
   trunk/scipy/optimize/setup.py
Log:
renamed NNLS.F to nnls.f

Modified: trunk/scipy/optimize/SConscript
===================================================================
--- trunk/scipy/optimize/SConscript	2008-09-30 06:48:15 UTC (rev 4756)
+++ trunk/scipy/optimize/SConscript	2008-09-30 07:00:07 UTC (rev 4757)
@@ -74,7 +74,7 @@
 env.NumpyPythonExtension('minpack2', source = src)
 
 # _nnls pyextension
-src = [pjoin('nnls', i) for i in ['NNLS.F', 'nnls.pyf']]
+src = [pjoin('nnls', i) for i in ['nnls.f', 'nnls.pyf']]
 env.NumpyPythonExtension('_nnls', source = src)
 
 # moduleTNC pyextension

Deleted: trunk/scipy/optimize/nnls/NNLS.F
===================================================================
--- trunk/scipy/optimize/nnls/NNLS.F	2008-09-30 06:48:15 UTC (rev 4756)
+++ trunk/scipy/optimize/nnls/NNLS.F	2008-09-30 07:00:07 UTC (rev 4757)
@@ -1,477 +0,0 @@
-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   
-
-
-       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   
-
-
-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   
-
-
-
-      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   

Added: trunk/scipy/optimize/nnls/nnls.f
===================================================================
--- trunk/scipy/optimize/nnls/nnls.f	2008-09-30 06:48:15 UTC (rev 4756)
+++ trunk/scipy/optimize/nnls/nnls.f	2008-09-30 07:00:07 UTC (rev 4757)
@@ -0,0 +1,477 @@
+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   
+
+
+       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   
+
+
+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   
+
+
+
+      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   

Modified: trunk/scipy/optimize/setup.py
===================================================================
--- trunk/scipy/optimize/setup.py	2008-09-30 06:48:15 UTC (rev 4756)
+++ trunk/scipy/optimize/setup.py	2008-09-30 07:00:07 UTC (rev 4757)
@@ -42,7 +42,7 @@
     config.add_extension('_slsqp', sources=[join('slsqp', x) for x in sources])
 
     config.add_extension('_nnls', sources=[join('nnls', x) \
-                                          for x in ["NNLS.F","nnls.pyf"]])
+                                          for x in ["nnls.f","nnls.pyf"]])
 
     config.add_data_dir('tests')
     config.add_data_dir('benchmarks')



More information about the Scipy-svn mailing list