program main c*********************************************************************72 c cc MAIN is the main program for LAWSON_PRB6. c C DEMONSTRATE THE USE OF THE SUBROUTINE LDP FOR LEAST C DISTANCE PROGRAMMING BY SOLVING THE CONSTRAINED LINE DATA FITTING C PROBLEM OF CHAPTER 23. C c Modified: c c 19 October 2008 c c Author: c c Charles Lawson, Richard Hanson c c Reference: c c Charles Lawson, Richard Hanson, c Solving Least Squares Problems, c SIAM, 1995, c ISBN: 0898713560, c LC: QA275.L38. c integer MX parameter(MX = 2) integer I, INDEX(3), J, L, MDE, MDGH, ME, MG, MODE, N, NP1 double precision E(4,MX), F(4), G(3,MX), G2(3,MX) double precision H(3), H2(3), ONE, RES, S(MX), SM, T(4) double precision W(4), WLDP(21), WORK(2*MX), X(MX) double precision Z(MX), ZERO, ZNORM parameter(ONE = 1.0d0, ZERO = 0.0d0) data T / 0.25d0, 0.5d0, 0.5d0, 0.8d0 / data W / 0.50d0, 0.6d0, 0.7d0, 1.2d0 / 110 format (' PROG6.. EXAMPLE OF CONSTRAINED CURVE FITTING'/ * 10x,'USING THE SUBROUTINE LDP.'// * 10x,'RELATED INTERMEDIATE QUANTITIES ARE PRINTED.') 120 format (/' V = ',2F10.5/(10X,2F10.5)) 130 format (/' F TILDA =',4F10.5/' S = ',2F10.5) 140 format (/' G TILDA =',2F10.5/(10X,2F10.5)) 150 format (/' H TILDA =',3F10.5) 160 format (/' Z = ',2F10.5) 170 format (/' THE COEFICIENTS OF THE FITTED LINE F(T)=X(1)*T+X(2)'/ * ' ARE X(1) = ',F10.5,' AND X(2) = ',F10.5) 180 format (/' THE CONSECUTIVE RESIDUALS ARE'/1X,4(I4,F10.5)) 190 format (/' THE RESIDUALS NORM IS ',F10.5) 200 format (/' MODE (FROM LDP) = ',I3,', ZNORM = ',F10.5) call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LAWSON_PRB6' write (*,110) MDE=4 MDGH=3 ME=4 MG=3 N=2 c C DEFINE THE LEAST SQUARES AND CONSTRAINT MATRICES. C DO 10 I=1,ME E(I,1)=T(I) E(I,2)= ONE 10 F(I)=W(I) G(1,1)= ONE G(1,2)= ZERO G(2,1)= ZERO G(2,2)= ONE G(3,1)=-ONE G(3,2)=-ONE H(1)= ZERO H(2)= ZERO H(3)=-ONE C C COMPUTE THE SINGULAR VALUE DECOMPOSITION OF THE MATRIX, E. C CALL SVDRS (E, MDE, ME, N, F, 1, 1, S, WORK) write (*,120) ((E(I,J),J=1,N),I=1,N) write (*,130) F,(S(J),J=1,N) C C GENERALLY RANK DETERMINATION AND LEVENBERG-MARQUARDT C STABILIZATION COULD BE INSERTED HERE. C C DEFINE THE CONSTRAINT MATRIX FOR THE Z COORDINATE SYSTEM. c DO 30 I=1,MG DO 30 J=1,N SM= ZERO DO 20 L=1,N 20 SM=SM+G(I,L)*E(L,J) 30 G2(I,J)=SM/S(J) c C DEFINE CONSTRAINT RT SIDE FOR THE Z COORDINATE SYSTEM. c DO 50 I=1,MG SM= ZERO DO 40 J=1,N 40 SM=SM+G2(I,J)*F(J) 50 H2(I)=H(I)-SM write (*,140) ((G2(I,J),J=1,N),I=1,MG) write (*,150) H2 C C SOLVE THE CONSTRAINED PROBLEM IN Z-COORDINATES. C CALL LDP (G2,MDGH,MG,N,H2,Z,ZNORM,WLDP,INDEX,MODE) write (*,200) MODE,ZNORM write (*,160) Z C C TRANSFORM BACK FROM Z-COORDINATES TO X-COORDINATES. c DO 60 J=1,N 60 Z(J)=(Z(J)+F(J))/S(J) DO 80 I=1,N SM= ZERO DO 70 J=1,N 70 SM=SM+E(I,J)*Z(J) 80 X(I)=SM RES=ZNORM**2 NP1=N+1 DO 90 I=NP1,ME 90 RES=RES+F(I)**2 RES=sqrt(RES) c C COMPUTE THE RESIDUALS. c DO 100 I=1,ME 100 F(I)=W(I)-X(1)*T(I)-X(2) write (*,170) (X(J),J=1,N) write (*,180) (I,F(I),I=1,ME) write (*,190) RES write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LAWSON_PRB6:' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) STOP END