C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 15 PR600200
C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974PR600300
C          DEMONSTRATE THE USE OF THE SUBROUTINE   LDP  FOR LEAST       PR600400
C     DISTANCE PROGRAMMING BY SOLVING THE CONSTRAINED LINE DATA FITTING PR600500
C     PROBLEM OF CHAPTER 23.                                            PR600600
C                                                                       PR600700
      DIMENSION E(4,2), F(4), G(3,2), H(3), G2(3,2), H2(3), X(2), Z(2), PR600800
     1W(4)                                                              PR600900
      DIMENSION WLDP(21), S(6), T(4)                                    PR601000
      INTEGER INDEX(3)                                                  PR601100
C                                                                       PR601200
C Set math to execute single precision in 23-bit mantissa format.               
      CALL MPBRQQ
C
      WRITE (6,110)                                                     PR601300
      MDE=4                                                             PR601400
      MDGH=3                                                            PR601500
C                                                                       PR601600
      ME=4                                                              PR601700
      MG=3                                                              PR601800
      N=2                                                               PR601900
C                DEFINE THE LEAST SQUARES AND CONSTRAINT MATRICES.      PR602000
      T(1)=0.25                                                         PR602100
      T(2)=0.5                                                          PR602200
      T(3)=0.5                                                          PR602300
      T(4)=0.8                                                          PR602400
C                                                                       PR602500
      W(1)=0.5                                                          PR602600
      W(2)=0.6                                                          PR602700
      W(3)=0.7                                                          PR602800
      W(4)=1.2                                                          PR602900
C                                                                       PR603000
          DO 10 I=1,ME                                                  PR603100
          E(I,1)=T(I)                                                   PR603200
          E(I,2)=1.                                                     PR603300
   10     F(I)=W(I)                                                     PR603400
C                                                                       PR603500
      G(1,1)=1.                                                         PR603600
      G(1,2)=0.                                                         PR603700
      G(2,1)=0.                                                         PR603800
      G(2,2)=1.                                                         PR603900
      G(3,1)=-1.                                                        PR604000
      G(3,2)=-1.                                                        PR604100
C                                                                       PR604200
      H(1)=0.                                                           PR604300
      H(2)=0.                                                           PR604400
      H(3)=-1.                                                          PR604500
C                                                                       PR604600
C     COMPUTE THE SINGULAR VALUE DECOMPOSITION OF THE MATRIX, E.        PR604700
C                                                                       PR604800
      CALL SVDRS (E,MDE,ME,N,F,1,1,S)                                   PR604900
C                                                                       PR605000
      WRITE (6,120) ((E(I,J),J=1,N),I=1,N)                              PR605100
      WRITE (6,130) F,(S(J),J=1,N)                                      PR605200
C                                                                       PR605300
C     GENERALLY RANK DETERMINATION AND LEVENBERG-MARQUARDT              PR605400
C     STABILIZATION COULD BE INSERTED HERE.                             PR605500
C                                                                       PR605600
C        DEFINE THE CONSTRAINT MATRIX FOR THE Z COORDINATE SYSTEM.      PR605700
          DO 30 I=1,MG                                                  PR605800
              DO 30 J=1,N                                               PR605900
              SM=0.                                                     PR606000
                  DO 20 L=1,N                                           PR606100
   20             SM=SM+G(I,L)*E(L,J)                                   PR606200
   30         G2(I,J)=SM/S(J)                                           PR606300
C         DEFINE CONSTRAINT RT SIDE FOR THE Z COORDINATE SYSTEM.        PR606400
          DO 50 I=1,MG                                                  PR606500
          SM=0.                                                         PR606600
              DO 40 J=1,N                                               PR606700
   40         SM=SM+G2(I,J)*F(J)                                        PR606800
   50     H2(I)=H(I)-SM                                                 PR606900
C                                                                       PR607000
      WRITE (6,140) ((G2(I,J),J=1,N),I=1,MG)                            PR607100
      WRITE (6,150) H2                                                  PR607200
C                                                                       PR607300
C                        SOLVE THE CONSTRAINED PROBLEM IN Z-COORDINATES.PR607400
C                                                                       PR607500
      CALL LDP (G2,MDGH,MG,N,H2,Z,ZNORM,WLDP,INDEX,MODE)                PR607600
C                                                                       PR607700
      WRITE (6,200) MODE,ZNORM                                          PR607800
      WRITE (6,160) Z                                                   PR607900
C                                                                       PR608000
C                    TRANSFORM BACK FROM Z-COORDINATES TO X-COORDINATES.PR608100
          DO 60 J=1,N                                                   PR608200
   60     Z(J)=(Z(J)+F(J))/S(J)                                         PR608300
          DO 80 I=1,N                                                   PR608400
          SM=0.                                                         PR608500
              DO 70 J=1,N                                               PR608600
   70         SM=SM+E(I,J)*Z(J)                                         PR608700
   80     X(I)=SM                                                       PR608800
      RES=ZNORM**2                                                      PR608900
      NP1=N+1                                                           PR609000
          DO 90 I=NP1,ME                                                PR609100
   90     RES=RES+F(I)**2                                               PR609200
      RES=SQRT(RES)                                                     PR609300
C                                    COMPUTE THE RESIDUALS.             PR609400
          DO 100 I=1,ME                                                 PR609500
  100     F(I)=W(I)-X(1)*T(I)-X(2)                                      PR609600
      WRITE (6,170) (X(J),J=1,N)                                        PR609700
      WRITE (6,180) (I,F(I),I=1,ME)                                     PR609800
      WRITE (6,190) RES                                                 PR609900
      STOP                                                              PR610000
C                                                                       PR610100
  110 FORMAT (46H0PROG6..  EXAMPLE OF CONSTRAINED CURVE FITTING,26H USINPR610200
     1G THE SUBROUTINE LDP.,/43H0RELATED INTERMEDIATE QUANTITIES ARE GIVPR610300
     2EN.)                                                              PR610400
  120 FORMAT (10H0V =      ,2F10.5/(10X,2F10.5))                        PR610500
  130 FORMAT (10H0F TILDA =,4F10.5/10H0S =      ,2F10.5)                PR610600
  140 FORMAT (10H0G TILDA =,2F10.5/(10X,2F10.5))                        PR610700
  150 FORMAT (10H0H TILDA =,3F10.5)                                     PR610800
  160 FORMAT (10H0Z =      ,2F10.5)                                     PR610900
  170 FORMAT (52H0THE COEFICIENTS OF THE FITTED LINE F(T)=X(1)*T+X(2),12PR611000
     1H ARE :           /
     2'     X(1) = ',F10.5,14H   AND X(2) = ,F10.5)                     PR611100
  180 FORMAT (30H0THE CONSECUTIVE RESIDUALS ARE/1X,4(I5,F10.5))         PR611200
  190 FORMAT (23H0THE RESIDUALS NORM IS ,F10.5)                         PR611300
  200 FORMAT (18H0MODE (FROM LDP) =,I3,2X,7HZNORM =,F10.5)              PR611400
C                                                                       PR611500
      END                                                               PR611600
