      SUBROUTINE LDP (G,MDG,M,N,H,X,XNORM,W,INDEX,MODE)                 LDP00100
C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1974 MAR 1  LDP00200
C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974LDP00300
C                                                                       LDP00400
C            **********  LEAST DISTANCE PROGRAMMING  **********         LDP00500
C                                                                       LDP00600
      INTEGER INDEX(M)                                                  LDP00700
      DIMENSION G(MDG,N), H(M), X(N), W(1)                              LDP00800
      ZERO=0.                                                           LDP00900
      ONE=1.                                                            LDP01000
      IF (N.LE.0) GO TO 120                                             LDP01100
          DO 10 J=1,N                                                   LDP01200
   10     X(J)=ZERO                                                     LDP01300
      XNORM=ZERO                                                        LDP01400
      IF (M.LE.0) GO TO 110                                             LDP01500
C                                                                       LDP01600
C     THE DECLARED DIMENSION OF W() MUST BE AT LEAST (N+1)*(M+2)+2*M.   LDP01700
C                                                                       LDP01800
C      FIRST (N+1)*M LOCS OF W()   =  MATRIX E FOR PROBLEM NNLS.        LDP01900
C       NEXT     N+1 LOCS OF W()   =  VECTOR F FOR PROBLEM NNLS.        LDP02000
C       NEXT     N+1 LOCS OF W()   =  VECTOR Z FOR PROBLEM NNLS.        LDP02100
C       NEXT       M LOCS OF W()   =  VECTOR Y FOR PROBLEM NNLS.        LDP02200
C       NEXT       M LOCS OF W()   =  VECTOR WDUAL FOR PROBLEM NNLS.    LDP02300
C     COPY G**T INTO FIRST N ROWS AND M COLUMNS OF E.                   LDP02400
C     COPY H**T INTO ROW N+1 OF E.                                      LDP02500
C                                                                       LDP02600
      IW=0                                                              LDP02700
          DO 30 J=1,M                                                   LDP02800
              DO 20 I=1,N                                               LDP02900
              IW=IW+1                                                   LDP03000
   20         W(IW)=G(J,I)                                              LDP03100
          IW=IW+1                                                       LDP03200
   30     W(IW)=H(J)                                                    LDP03300
      IF=IW+1                                                           LDP03400
C                                STORE N ZEROS FOLLOWED BY A ONE INTO F.LDP03500
          DO 40 I=1,N                                                   LDP03600
          IW=IW+1                                                       LDP03700
   40     W(IW)=ZERO                                                    LDP03800
      W(IW+1)=ONE                                                       LDP03900
C                                                                       LDP04000
      NP1=N+1                                                           LDP04100
      IZ=IW+2                                                           LDP04200
      IY=IZ+NP1                                                         LDP04300
      IWDUAL=IY+M                                                       LDP04400
C                                                                       LDP04500
      CALL NNLS (W,NP1,NP1,M,W(IF),W(IY),RNORM,W(IWDUAL),W(IZ),INDEX,   LDP04600
     *           MODE)                                                  LDP04700
C                      USE THE FOLLOWING RETURN IF UNSUCCESSFUL IN NNLS.LDP04800
      IF (MODE.NE.1) RETURN                                             LDP04900
      IF (RNORM) 130,130,50                                             LDP05000
   50 FAC=ONE                                                           LDP05100
      IW=IY-1                                                           LDP05200
          DO 60 I=1,M                                                   LDP05300
          IW=IW+1                                                       LDP05400
C                               HERE WE ARE USING THE SOLUTION VECTOR Y.LDP05500
   60     FAC=FAC-H(I)*W(IW)                                            LDP05600
C                                                                       LDP05700
      IF (DIFF(ONE+FAC,ONE)) 130,130,70                                 LDP05800
   70 FAC=ONE/FAC                                                       LDP05900
          DO 90 J=1,N                                                   LDP06000
          IW=IY-1                                                       LDP06100
              DO 80 I=1,M                                               LDP06200
              IW=IW+1                                                   LDP06300
C                               HERE WE ARE USING THE SOLUTION VECTOR Y.LDP06400
   80         X(J)=X(J)+G(I,J)*W(IW)                                    LDP06500
   90     X(J)=X(J)*FAC                                                 LDP06600
          DO 100 J=1,N                                                  LDP06700
  100     XNORM=XNORM+X(J)**2                                           LDP06800
      XNORM=SQRT(XNORM)                                                 LDP06900
C                             SUCCESSFUL RETURN.                        LDP07000
  110 MODE=1                                                            LDP07100
      RETURN                                                            LDP07200
C                             ERROR RETURN.       N .LE. 0.             LDP07300
  120 MODE=2                                                            LDP07400
      RETURN                                                            LDP07500
C                             RETURNING WITH CONSTRAINTS NOT COMPATIBLE.LDP07600
  130 MODE=4                                                            LDP07700
      RETURN                                                            LDP07800
      END                                                               LDP07900
