[Scipy-svn] r4757 - in trunk/scipy/optimize: . nnls
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Sep 30 03:00:13 EDT 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