program main c*********************************************************************72 c cc MAIN is the main program for LAWSON_PRB05. c C EXAMPLE OF THE USE OF SUBROUTINES BNDACC AND BNDSOL TO SOLVE C SEQUENTIALLY THE BANDED LEAST SQUARES PROBLEM THAT ARISES IN C SPLINE CURVE FITTING. 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 mdg parameter ( mdg = 12 ) integer mxy parameter ( mxy = 12 ) integer nband parameter ( nband = 4 ) double precision b(10), c(mxy), cov(mxy), g(mdg,5), h integer i integer ic integer ig integer ip integer ir integer j integer jt integer l integer m integer mt integer nbp integer nc double precision one, p1, p2, q(4), r, rdummy, rnorm double precision sigfac, sigsq, u, x(mxy), y(mxy), yfit, zero parameter(one = 1.0d0, zero = 0.0d0) data y /2.2d0, 4.0d0, 5.0d0, 4.6d0, 2.8d0, 2.7d0, * 3.8d0, 5.1d0, 6.1d0, 6.3d0, 5.0d0, 2.0d0/ 160 format (/' I',8X,'X',10X,'Y',6X,'YFIT',4X,'R=Y-YFIT/1X') 170 format (1X,I3,4X,F6.0,4X,F6.2,4X,F6.2,4X,F8.4) 180 format (/' C =',6F10.5/(4X,6F10.5)) 190 format (3(2X,2I4,g15.7)) 200 format (/' COVARIANCE MATRIX OF THE SPLINE COEFFICIENTS.') 210 format (/' RNORM =',g15.8) 220 format (/' SIGFAC =',g15.8) call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LAWSON_PRB5' write ( *, '(a)' ) ' BNDACC accumulates a banded matrix.' write ( *, '(a)' ) ' BNDSOL solves an associated banded least' write ( *, '(a)' ) ' squares problem.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Execute a sequence of cubic spline fits' write ( *, '(a)' ) ' to a discrete set of data.' M=MDG c c Set the data points. c do i = 1, m x(i) = dble ( 2 * i ) end do C C BEGIN LOOP THRU CASES USING INCREASING NOS OF BREAKPOINTS. C DO 150 NBP = 5, 10 write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The number of breakpoints is ', nbp NC = NBP + 2 c C SET BREAKPOINTS c B(1)=X(1) B(NBP)=X(M) H=(B(NBP)-B(1))/(NBP-1) DO I=3,NBP B(I-1)=B(I-2)+H end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Breakpoints:' write ( *, '(a)' ) ' ' write (*,'(5(2x,g14.6))' ) ( b(i), i = 1, nbp ) C C INITIALIZE IR AND IP BEFORE FIRST CALL TO BNDACC. C IR=1 IP=1 I=1 JT=1 40 continue MT=0 50 CONTINUE IF (X(I).GT.B(JT+1)) GO TO 60 C C SET ROW FOR ITH DATA POINT C U=(X(I)-B(JT))/H IG=IR+MT G(IG,1)=P1(ONE - U) G(IG,2)=P2(ONE - U) G(IG,3)=P2(U) G(IG,4)=P1(U) G(IG,5)=Y(I) MT=MT+1 IF (I.EQ.M) GO TO 60 I=I+1 GO TO 50 C C SEND BLOCK OF DATA TO PROCESSOR C 60 CONTINUE CALL BNDACC (G,MDG,NBAND,IP,IR,MT,JT) IF (I.EQ.M) GO TO 70 JT=JT+1 GO TO 40 c C COMPUTE SOLUTION C() c 70 CONTINUE CALL BNDSOL (1,G,MDG,NBAND,IP,IR,C,NC,RNORM) c C WRITE SOLUTION COEFFICIENTS c write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'C:' write ( *, '(a)' ) ' ' write ( *, '(5(2x,g14.6))' ) ( c(l), l = 1, nc ) write (*,210) RNORM C C COMPUTE AND PRINT X,Y,YFIT,R=Y-YFIT C write (*,160) JT=1 DO 110 I=1,M 80 IF (X(I).LE.B(JT+1)) GO TO 90 JT=JT+1 GO TO 80 90 U=(X(I)-B(JT))/H Q(1)=P1(ONE - U) Q(2)=P2(ONE - U) Q(3)=P2(U) Q(4)=P1(U) YFIT=ZERO DO 100 L=1,4 IC=JT-1+L 100 YFIT=YFIT+C(IC)*Q(L) R=Y(I)-YFIT write (*,170) I,X(I),Y(I),YFIT,R 110 CONTINUE C C COMPUTE RESIDUAL VECTOR NORM. C IF (M.LE.NC) GO TO 150 SIGSQ=(RNORM**2)/(M-NC) SIGFAC=sqrt(SIGSQ) write (*,220) SIGFAC write (*,200) C C COMPUTE AND PRINT COLS. OF COVARIANCE. C DO 140 J=1,NC DO 120 I=1,NC 120 COV(I)=ZERO COV(J)= ONE CALL BNDSOL (2,G,MDG,NBAND,IP,IR,COV,NC,RDUMMY) CALL BNDSOL (3,G,MDG,NBAND,IP,IR,COV,NC,RDUMMY) C C COMPUTE THE JTH COL. OF THE COVARIANCE MATRIX. c DO 130 I=1,NC 130 COV(I)=COV(I)*SIGSQ 140 write (*,190) (L,J,COV(L),L=1,NC) 150 CONTINUE write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LAWSON_PRB5:' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) STOP END double precision function P1(T) c*********************************************************************72 c cc P1 c double precision T P1 = 0.25d0 * T**2 * T return end double precision function P2(T) c*********************************************************************72 c cc P2 c double precision T P2 = -(1.0d0 - T)**2 * (1.0d0 + T) * 0.75d0 + 1.0d0 return end