C SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE,MAXITER) 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 MAXITER ITERATIONS. C MAXITER THE MAXIMUM ALLOWED NUMBER OF ITERATIONS. C IF NEGATIVE, THE LIMIT IS TAKEN AS THE DEFAULT, 3*N. C C ------------------------------------------------------------------ SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE,MAXITER) C ------------------------------------------------------------------ integer I, II, IP, ITER, ITMAX, IZ, IZ1, IZ2, IZMAX, J, JJ, JZ, L integer M, MDA, MODE,N, NPP1, NSETP, RTNKEY integer MAXITER 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(1), 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 IF (MAXITER .le. 0) then ITMAX=3*N ELSE ITMAX = MAXITER ENDIF 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