C SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) NLS00100 C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUNE 15NLS00200 C TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974NLS00300 C NLS00400 C ********** NONNEGATIVE LEAST SQUARES ********** NLS00500 C NLS00600 C GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN NLS00700 C N-VECTOR, X, WHICH SOLVES THE LEAST SQUARES PROBLEM NLS00800 C NLS00900 C A * X = B SUBJECT TO X .GE. 0 NLS01000 C NLS01100 C A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE NLS01200 C ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N NLS01300 C MATRIX, A. ON EXIT A() CONTAINS NLS01400 C THE PRODUCT MATRIX, Q*A , WHERE Q IS AN NLS01500 C M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY NLS01600 C THIS SUBROUTINE. NLS01700 C B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() CON- NLS01800 C TAINS Q*B. NLS01900 C X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL NLS02000 C CONTAIN THE SOLUTION VECTOR. NLS02100 C RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE NLS02200 C RESIDUAL VECTOR. NLS02300 C W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN NLS02400 C THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0. NLS02500 C FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z NLS02600 C ZZ() AN M-ARRAY OF WORKING SPACE. NLS02700 C INDEX() AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N. NLS02800 C ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS NLS02900 C P AND Z AS FOLLOWS.. NLS03000 C NLS03100 C INDEX(1) THRU INDEX(NSETP) = SET P. NLS03200 C INDEX(IZ1) THRU INDEX(IZ2) = SET Z. NLS03300 C IZ1 = NSETP + 1 = NPP1 NLS03400 C IZ2 = N NLS03500 C MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING NLS03600 C MEANINGS. NLS03700 C 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY. NLS03800 C 2 THE DIMENSIONS OF THE PROBLEM ARE BAD. NLS03900 C EITHER M .LE. 0 OR N .LE. 0. NLS04000 C 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS. NLS04100 C NLS04200 SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) NLS04300 DIMENSION A(MDA,N), B(M), X(N), W(N), ZZ(M) NLS04400 INTEGER INDEX(N) NLS04500 ZERO=0. NLS04600 ONE=1. NLS04700 TWO=2. NLS04800 FACTOR=0.01 NLS04900 C NLS05000 MODE=1 NLS05100 IF (M.GT.0.AND.N.GT.0) GO TO 10 NLS05200 MODE=2 NLS05300 RETURN NLS05400 10 ITER=0 NLS05500 ITMAX=3*N NLS05600 C NLS05700 C INITIALIZE THE ARRAYS INDEX() AND X(). NLS05800 C NLS05900 DO 20 I=1,N NLS06000 X(I)=ZERO NLS06100 20 INDEX(I)=I NLS06200 C NLS06300 IZ2=N NLS06400 IZ1=1 NLS06500 NSETP=0 NLS06600 NPP1=1 NLS06700 C ****** MAIN LOOP BEGINS HERE ****** NLS06800 30 CONTINUE NLS06900 C QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION.NLS07000 C OR IF M COLS OF A HAVE BEEN TRIANGULARIZED. NLS07100 C NLS07200 IF (IZ1.GT.IZ2.OR.NSETP.GE.M) GO TO 350 NLS07300 C NLS07400 C COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W().NLS07500 C NLS07600 DO 50 IZ=IZ1,IZ2 NLS07700 J=INDEX(IZ) NLS07800 SM=ZERO NLS07900 DO 40 L=NPP1,M NLS08000 40 SM=SM+A(L,J)*B(L) NLS08100 50 W(J)=SM NLS08200 C FIND LARGEST POSITIVE W(J). NLS08300 60 WMAX=ZERO NLS08400 DO 70 IZ=IZ1,IZ2 NLS08500 J=INDEX(IZ) NLS08600 IF (W(J).LE.WMAX) GO TO 70 NLS08700 WMAX=W(J) NLS08800 IZMAX=IZ NLS08900 70 CONTINUE NLS09000 C NLS09100 C IF WMAX .LE. 0. GO TO TERMINATION. NLS09200 C THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS.NLS09300 C NLS09400 IF (WMAX) 350,350,80 NLS09500 80 IZ=IZMAX NLS09600 J=INDEX(IZ) NLS09700 C NLS09800 C THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P. NLS09900 C BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID NLS10000 C NEAR LINEAR DEPENDENCE. NLS10100 C NLS10200 ASAVE=A(NPP1,J) NLS10300 CALL H12 (1,NPP1,NPP1+1,M,A(1,J),1,UP,DUMMY,1,1,0) NLS10400 UNORM=ZERO NLS10500 IF (NSETP.EQ.0) GO TO 100 NLS10600 DO 90 L=1,NSETP NLS10700 90 UNORM=UNORM+A(L,J)**2 NLS10800 100 UNORM=SQRT(UNORM) NLS10900 IF (DIFF(UNORM+ABS(A(NPP1,J))*FACTOR,UNORM)) 130,130,110 NLS11000 C NLS11100 C COL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE ZZ AND NLS11200 C > SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ). NLS11300 C NLS11400 110 DO 120 L=1,M NLS11500 120 ZZ(L)=B(L) NLS11600 CALL H12 (2,NPP1,NPP1+1,M,A(1,J),1,UP,ZZ,1,1,1) NLS11700 ZTEST=ZZ(NPP1)/A(NPP1,J) NLS11800 C NLS11900 C SEE IF ZTEST IS POSITIVE NLS12000 C REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P. NLS12400 C RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL NLS12500 C NLS12100 IF (ZTEST) 130,130,140 NLS12200 C NLS12300 C COEFFS AGAIN. NLS12600 C NLS12700 130 A(NPP1,J)=ASAVE NLS12800 W(J)=ZERO NLS12900 GO TO 60 NLS13000 C NLS13100 C THE INDEX J=INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM NLS13200 C SET Z TO SET P. UPDATE B, UPDATE INDICES, APPLY HOUSEHOLDER NLS13300 C TRANSFORMATIONS TO COLS IN NEW SET Z, ZERO SUBDIAGONAL ELTS IN NLS13400 C COL J, SET W(J)=0. NLS13500 C NLS13600 140 DO 150 L=1,M NLS13700 150 B(L)=ZZ(L) NLS13800 C NLS13900 INDEX(IZ)=INDEX(IZ1) NLS14000 INDEX(IZ1)=J NLS14100 IZ1=IZ1+1 NLS14200 NSETP=NPP1 NLS14300 NPP1=NPP1+1 NLS14400 C NLS14500 IF (IZ1.GT.IZ2) GO TO 170 NLS14600 DO 160 JZ=IZ1,IZ2 NLS14700 JJ=INDEX(JZ) NLS14800 160 CALL H12 (2,NSETP,NPP1,M,A(1,J),1,UP,A(1,JJ),1,MDA,1) NLS14900 170 CONTINUE NLS15000 C NLS15100 IF (NSETP.EQ.M) GO TO 190 NLS15200 DO 180 L=NPP1,M NLS15300 180 A(L,J)=ZERO NLS15400 190 CONTINUE NLS15500 C NLS15600 W(J)=ZERO NLS15700 C SOLVE THE TRIANGULAR SYSTEM. NLS15800 C STORE THE SOLUTION TEMPORARILY IN ZZ().NLS15900 ASSIGN 200 TO NEXT NLS16000 GO TO 400 NLS16100 200 CONTINUE NLS16200 C NLS16300 C ****** SECONDARY LOOP BEGINS HERE ****** NLS16400 C NLS16500 C ITERATION COUNTER. NLS16600 C NLS16700 210 ITER=ITER+1 NLS16800 IF (ITER.LE.ITMAX) GO TO 220 NLS16900 MODE=3 NLS17000 WRITE (6,440) NLS17100 GO TO 350 NLS17200 220 CONTINUE NLS17300 C NLS17400 C SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE. NLS17500 C IF NOT COMPUTE ALPHA. NLS17600 C NLS17700 ALPHA=TWO NLS17800 DO 240 IP=1,NSETP NLS17900 L=INDEX(IP) NLS18000 IF (ZZ(IP)) 230,230,240 NLS18100 C NLS18200 230 T=-X(L)/(ZZ(IP)-X(L)) NLS18300 IF (ALPHA.LE.T) GO TO 240 NLS18400 ALPHA=T NLS18500 JJ=IP NLS18600 240 CONTINUE NLS18700 C NLS18800 C IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL NLS18900 C STILL = 2. IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP. NLS19000 C NLS19100 IF (ALPHA.EQ.TWO) GO TO 330 NLS19200 C NLS19300 C OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO NLS19400 C INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ. NLS19500 C NLS19600 DO 250 IP=1,NSETP NLS19700 L=INDEX(IP) NLS19800 250 X(L)=X(L)+ALPHA*(ZZ(IP)-X(L)) NLS19900 C NLS20000 C MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I NLS20100 C FROM SET P TO SET Z. NLS20200 C NLS20300 I=INDEX(JJ) NLS20400 260 X(I)=ZERO NLS20500 C NLS20600 IF (JJ.EQ.NSETP) GO TO 290 NLS20700 JJ=JJ+1 NLS20800 DO 280 J=JJ,NSETP NLS20900 II=INDEX(J) NLS21000 INDEX(J-1)=II NLS21100 CALL G1 (A(J-1,II),A(J,II),CC,SS,A(J-1,II)) NLS21200 A(J,II)=ZERO NLS21300 DO 270 L=1,N NLS21400 IF (L.NE.II) CALL G2 (CC,SS,A(J-1,L),A(J,L)) NLS21500 270 CONTINUE NLS21600 280 CALL G2 (CC,SS,B(J-1),B(J)) NLS21700 290 NPP1=NSETP NLS21800 NSETP=NSETP-1 NLS21900 IZ1=IZ1-1 NLS22000 INDEX(IZ1)=I NLS22100 C NLS22200 C SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE. THEY SHOULDNLS22300 C BE BECAUSE OF THE WAY ALPHA WAS DETERMINED. NLS22400 C IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR. ANY NLS22500 C THAT ARE NONPOSITIVE WILL BE SET TO ZERO NLS22600 C AND MOVED FROM SET P TO SET Z. NLS22700 C NLS22800 DO 300 JJ=1,NSETP NLS22900 I=INDEX(JJ) NLS23000 IF (X(I)) 260,260,300 NLS23100 300 CONTINUE NLS23200 C NLS23300 C COPY B( ) INTO ZZ( ). THEN SOLVE AGAIN AND LOOP BACK. NLS23400 C NLS23500 DO 310 I=1,M NLS23600 310 ZZ(I)=B(I) NLS23700 ASSIGN 320 TO NEXT NLS23800 GO TO 400 NLS23900 320 CONTINUE NLS24000 GO TO 210 NLS24100 C ****** END OF SECONDARY LOOP ****** NLS24200 C NLS24300 330 DO 340 IP=1,NSETP NLS24400 I=INDEX(IP) NLS24500 340 X(I)=ZZ(IP) NLS24600 C ALL NEW COEFFS ARE POSITIVE. LOOP BACK TO BEGINNING. NLS24700 GO TO 30 NLS24800 C NLS24900 C ****** END OF MAIN LOOP ****** NLS25000 C NLS25100 C COME TO HERE FOR TERMINATION. NLS25200 C COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR. NLS25300 C NLS25400 350 SM=ZERO NLS25500 IF (NPP1.GT.M) GO TO 370 NLS25600 DO 360 I=NPP1,M NLS25700 360 SM=SM+B(I)**2 NLS25800 GO TO 390 NLS25900 370 DO 380 J=1,N NLS26000 380 W(J)=ZERO NLS26100 390 RNORM=SQRT(SM) NLS26200 RETURN NLS26300 C NLS26400 C THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE NLS26500 C TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ(). NLS26600 C NLS26700 400 DO 430 L=1,NSETP NLS26800 IP=NSETP+1-L NLS26900 IF (L.EQ.1) GO TO 420 NLS27000 DO 410 II=1,IP NLS27100 410 ZZ(II)=ZZ(II)-A(II,JJ)*ZZ(IP+1) NLS27200 420 JJ=INDEX(IP) NLS27300 430 ZZ(IP)=ZZ(IP)/A(IP,JJ) NLS27400 GO TO NEXT, (200,320) NLS27500 440 FORMAT (35H0 NNLS QUITTING ON ITERATION COUNT.) NLS27600 END NLS27700