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
