SUBROUTINE NNLS(A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) C C C.L.LAWSON AND R.J.HANSON,JET PROPULSION LABORATORY, 1973 JUNE 15 C TO APPEAR IN ^SOLVING LEAST SQUARES PROBLEMS^, PRENTICE HALL,1974 C C *****NON-NEGATIVE LEAST SQUARES***** C C GIVEN AN M BY N MATRIX, A, AND AN M VECTOR, B, COMPUTE AN N C VECTOR, X, WHICH SOLVES THE LEAST SQUARES PROBLEM C C A*X = B SUBJECT TO X >= 0, AND 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 THE PRODUCT C MATRIX, Q*A, WHERE Q IS A M BY M ORTHOGONAL C MATRIX GENERATED IMPLICITLY BY THIS SUBROUTINE. C B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() C CONTAINS 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 EUCLIDIAN NORN OF THE RESIDUAL C VECTOR C W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN THE C DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0. FOR ALL C I IN THE SET P AND W(I) .LE. 0 FOR ALL I IN THE 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) = SETP 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 IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C DIMENSION A(MDA,N),B(M),X(N),W(N),ZZ(M),INDEX(N) ZERO = 0.0 ONE = 1.0 TWO = 2.0 FACTOR = 0.01 C MODE = 1 IF(M.GT.0.AND.N.GT.0) GO TO 10 MODE = 2 RETURN 10 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 C *****MAIN LOOP BEGINS HERE***** C 30 CONTINUE C C QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION OF IF M C COLUMNS 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) 50 W(J) = SM C C FIND LARGEST POSITIVE W(J) C 60 WMAX = ZERO DO 70 IZ=IZ1,IZ2 J = INDEX(IZ) IF(W(J).LE.WMAX) GO TO 70 WMAX = W(J) IZMAX = IZ 70 CONTINUE C C IF WMAX .LE. 0 GO TO TERMINATION C THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS C IF(WMAX) 350,350,80 80 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.EQ.0) GO TO 100 DO 90 L=1,NSETP 90 UNORM = UNORM + A(L,J)**2 100 UNORM = SQRT(UNORM) IF(DIFF(UNORM+ABS(A(NPP1,J))*FACTOR,UNORM)) 130,130,110 C C SOL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE ZZ AND C SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J)). C 110 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) 130,130,140 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 130 A(NPP1,J) = ASAVE W(J) = ZERO GO TO 60 C C THE INDEX J = INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM THE 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 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.GT.IZ2) GO TO 170 DO 160 JZ=IZ1,IZ2 JJ = INDEX(JZ) 160 CALL H12(2,NSETP,NPP1,M,A(1,J),1,UP,A(1,JJ),1,MDA,1) 170 CONTINUE C IF(NSETP.EQ.M) GO TO 190 DO 180 L=NPP1,M 180 A(L,J) = ZERO 190 CONTINUE C W(J) = ZERO C C SOLVE THE TRIANGULAR SYSTEM C STORE THE SOLUTION TEMPORARILY IN ZZ() C ASSIGN 200 TO NEXT GO TO 400 200 CONTINUE C C *****SECONDARY LOOP BEGINS HERE***** C C ITERATION COUNTER C 210 ITER = ITER + 1 IF(ITER.LE.ITMAX) GO TO 220 MODE = 3 WRITE(6,500) GO TO 350 220 CONTINUE 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)) 230,230,240 C 230 T = -X(L)/(ZZ(IP)-X(L)) IF(ALPHA.LE.T) GO TO 240 ALPHA = T JJ = IP 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) 250 X(L) = X(L) + ALPHA*(ZZ(IP)-X(L)) 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 X(I) = ZERO C IF(JJ.EQ.NSETP) GO TO 290 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) CALL G2(CC,SS,A(J-1,L),A(J,L)) 270 CONTINUE 280 CALL G2(CC,SS,B(J-1),B(J)) 290 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 NON-POSITIVE 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)) 260,260,300 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) ASSIGN 320 TO NEXT GO TO 400 320 CONTINUE GO TO 210 C C *****END OF SECONDARY LOOP***** C 330 DO 340 IP=1,NSETP I = INDEX(IP) 340 X(I) = ZZ(IP) C C ALL NEW COEFFICIENTS ARE POSITIVE. LOOP BACK TO BEGINNING. C 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 SM = ZERO IF(NPP1.GT.M) GO TO 370 DO 360 I=NPP1,M 360 SM = SM + B(I)**2 GO TO 390 370 DO 380 J=1,N 380 W(J) = ZERO 390 RNORM = SQRT(SM) RETURN C 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 DO 430 L=1,NSETP IP = NSETP + 1 - L IF(L.EQ.1) GO TO 420 DO 410 II=1,IP 410 ZZ(II) = ZZ(II) - A(II,JJ)*ZZ(IP+1) 420 JJ = INDEX(IP) 430 ZZ(IP) = ZZ(IP)/A(IP,JJ) GO TO NEXT,(200,320) 500 FORMAT(35H0 NNLS QUITTING ON ITERATION COUNT ) END