home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / FORTRAN / LLSQ.ZIP / LDP.FOR < prev    next >
Encoding:
Text File  |  1984-02-23  |  6.3 KB  |  80 lines
  1.       SUBROUTINE LDP (G,MDG,M,N,H,X,XNORM,W,INDEX,MODE)                 LDP00100
  2. C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1974 MAR 1  LDP00200
  3. C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974LDP00300
  4. C                                                                       LDP00400
  5. C            **********  LEAST DISTANCE PROGRAMMING  **********         LDP00500
  6. C                                                                       LDP00600
  7.       INTEGER INDEX(M)                                                  LDP00700
  8.       DIMENSION G(MDG,N), H(M), X(N), W(1)                              LDP00800
  9.       ZERO=0.                                                           LDP00900
  10.       ONE=1.                                                            LDP01000
  11.       IF (N.LE.0) GO TO 120                                             LDP01100
  12.           DO 10 J=1,N                                                   LDP01200
  13.    10     X(J)=ZERO                                                     LDP01300
  14.       XNORM=ZERO                                                        LDP01400
  15.       IF (M.LE.0) GO TO 110                                             LDP01500
  16. C                                                                       LDP01600
  17. C     THE DECLARED DIMENSION OF W() MUST BE AT LEAST (N+1)*(M+2)+2*M.   LDP01700
  18. C                                                                       LDP01800
  19. C      FIRST (N+1)*M LOCS OF W()   =  MATRIX E FOR PROBLEM NNLS.        LDP01900
  20. C       NEXT     N+1 LOCS OF W()   =  VECTOR F FOR PROBLEM NNLS.        LDP02000
  21. C       NEXT     N+1 LOCS OF W()   =  VECTOR Z FOR PROBLEM NNLS.        LDP02100
  22. C       NEXT       M LOCS OF W()   =  VECTOR Y FOR PROBLEM NNLS.        LDP02200
  23. C       NEXT       M LOCS OF W()   =  VECTOR WDUAL FOR PROBLEM NNLS.    LDP02300
  24. C     COPY G**T INTO FIRST N ROWS AND M COLUMNS OF E.                   LDP02400
  25. C     COPY H**T INTO ROW N+1 OF E.                                      LDP02500
  26. C                                                                       LDP02600
  27.       IW=0                                                              LDP02700
  28.           DO 30 J=1,M                                                   LDP02800
  29.               DO 20 I=1,N                                               LDP02900
  30.               IW=IW+1                                                   LDP03000
  31.    20         W(IW)=G(J,I)                                              LDP03100
  32.           IW=IW+1                                                       LDP03200
  33.    30     W(IW)=H(J)                                                    LDP03300
  34.       IF=IW+1                                                           LDP03400
  35. C                                STORE N ZEROS FOLLOWED BY A ONE INTO F.LDP03500
  36.           DO 40 I=1,N                                                   LDP03600
  37.           IW=IW+1                                                       LDP03700
  38.    40     W(IW)=ZERO                                                    LDP03800
  39.       W(IW+1)=ONE                                                       LDP03900
  40. C                                                                       LDP04000
  41.       NP1=N+1                                                           LDP04100
  42.       IZ=IW+2                                                           LDP04200
  43.       IY=IZ+NP1                                                         LDP04300
  44.       IWDUAL=IY+M                                                       LDP04400
  45. C                                                                       LDP04500
  46.       CALL NNLS (W,NP1,NP1,M,W(IF),W(IY),RNORM,W(IWDUAL),W(IZ),INDEX,   LDP04600
  47.      *           MODE)                                                  LDP04700
  48. C                      USE THE FOLLOWING RETURN IF UNSUCCESSFUL IN NNLS.LDP04800
  49.       IF (MODE.NE.1) RETURN                                             LDP04900
  50.       IF (RNORM) 130,130,50                                             LDP05000
  51.    50 FAC=ONE                                                           LDP05100
  52.       IW=IY-1                                                           LDP05200
  53.           DO 60 I=1,M                                                   LDP05300
  54.           IW=IW+1                                                       LDP05400
  55. C                               HERE WE ARE USING THE SOLUTION VECTOR Y.LDP05500
  56.    60     FAC=FAC-H(I)*W(IW)                                            LDP05600
  57. C                                                                       LDP05700
  58.       IF (DIFF(ONE+FAC,ONE)) 130,130,70                                 LDP05800
  59.    70 FAC=ONE/FAC                                                       LDP05900
  60.           DO 90 J=1,N                                                   LDP06000
  61.           IW=IY-1                                                       LDP06100
  62.               DO 80 I=1,M                                               LDP06200
  63.               IW=IW+1                                                   LDP06300
  64. C                               HERE WE ARE USING THE SOLUTION VECTOR Y.LDP06400
  65.    80         X(J)=X(J)+G(I,J)*W(IW)                                    LDP06500
  66.    90     X(J)=X(J)*FAC                                                 LDP06600
  67.           DO 100 J=1,N                                                  LDP06700
  68.   100     XNORM=XNORM+X(J)**2                                           LDP06800
  69.       XNORM=SQRT(XNORM)                                                 LDP06900
  70. C                             SUCCESSFUL RETURN.                        LDP07000
  71.   110 MODE=1                                                            LDP07100
  72.       RETURN                                                            LDP07200
  73. C                             ERROR RETURN.       N .LE. 0.             LDP07300
  74.   120 MODE=2                                                            LDP07400
  75.       RETURN                                                            LDP07500
  76. C                             RETURNING WITH CONSTRAINTS NOT COMPATIBLE.LDP07600
  77.   130 MODE=4                                                            LDP07700
  78.       RETURN                                                            LDP07800
  79.       END                                                               LDP07900
  80.