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