C PROG5 PR500100 C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 PR500200 C TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974PR500300 C EXAMPLE OF THE USE OF SUBROUTINES BNDACC AND BNDSOL TO SOLVE PR500400 C SEQUENTIALLY THE BANDED LEAST SQUARES PROBLEM WHICH ARISES IN PR500500 C SPLINE CURVE FITTING. PR500600 C PR500700 DIMENSION X(12),Y(12),B(10),G(12,5),C(12),Q(4),COV(12) PR500800 C PR500900 C DEFINE FUNCTIONS P1 AND P2 TO BE USED IN CONSTRUCTING PR501000 C CUBIC SPLINES OVER UNIFORMLY SPACED BREAKPOINTS. PR501100 C PR501200 P1(T)=.25*T**2*T PR501300 P2(T)=-(1.-T)**2*(1.+T)*.75+1. PR501400 CALL MPBRQQ ZERO=0. PR501500 C PR501600 WRITE (6,230) PR501700 MDG=12 PR501800 NBAND=4 PR501900 M=12 PR502000 C SET ORDINATES OF DATA PR502100 Y(1)=2.2 PR502200 Y(2)=4.0 PR502300 Y(3)=5.0 PR502400 Y(4)=4.6 PR502500 Y(5)=2.8 PR502600 Y(6)=2.7 PR502700 Y(7)=3.8 PR502800 Y(8)=5.1 PR502900 Y(9)=6.1 PR503000 Y(10)=6.3 PR503100 Y(11)=5.0 PR503200 Y(12)=2.0 PR503300 C SET ABCISSAS OF DATA PR503400 DO 10 I=1,M PR503500 10 X(I)=2*I PR503600 C PR503700 C BEGIN LOOP THRU CASES USING INCREASING NOS OF BREAKPOINTS. PR503800 C PR503900 DO 150 NBP=5,10 PR504000 NC=NBP+2 PR504100 C SET BREAKPOINTS PR504200 B(1)=X(1) PR504300 B(NBP)=X(M) PR504400 H=(B(NBP)-B(1))/FLOAT(NBP-1) PR504500 IF (NBP.LE.2) GO TO 30 PR504600 DO 20 I=3,NBP PR504700 20 B(I-1)=B(I-2)+H PR504800 30 CONTINUE PR504900 WRITE (6,240) NBP,(B(I),I=1,NBP) PR505000 C PR505100 C INITIALIZE IR AND IP BEFORE FIRST CALL TO BNDACC. PR505200 C PR505300 IR=1 PR505400 IP=1 PR505500 I=1 PR505600 JT=1 PR505700 40 MT=0 PR505800 50 CONTINUE PR505900 IF (X(I).GT.B(JT+1)) GO TO 60 PR506000 C PR506100 C SET ROW FOR ITH DATA POINT PR506200 C PR506300 U=(X(I)-B(JT))/H PR506400 IG=IR+MT PR506500 G(IG,1)=P1(1.-U) PR506600 G(IG,2)=P2(1.-U) PR506700 G(IG,3)=P2(U) PR506800 G(IG,4)=P1(U) PR506900 G(IG,5)=Y(I) PR507000 MT=MT+1 PR507100 IF (I.EQ.M) GO TO 60 PR507200 I=I+1 PR507300 GO TO 50 PR507400 C PR507500 C SEND BLOCK OF DATA TO PROCESSOR PR507600 C PR507700 60 CONTINUE PR507800 CALL BNDACC (G,MDG,NBAND,IP,IR,MT,JT) PR507900 IF (I.EQ.M) GO TO 70 PR508000 JT=JT+1 PR508100 GO TO 40 PR508200 C COMPUTE SOLUTION C() PR508300 70 CONTINUE PR508400 CALL BNDSOL (1,G,MDG,NBAND,IP,IR,C,NC,RNORM) PR508500 C WRITE SOLUTION COEFFICIENTS PR508600 WRITE (6,180) (C(L),L=1,NC) PR508700 WRITE (6,210) RNORM PR508800 C PR508900 C COMPUTE AND PRINT X,Y,YFIT,R=Y-YFIT PR509000 C PR509100 WRITE (6,160) PR509200 JT=1 PR509300 DO 110 I=1,M PR509400 80 IF (X(I).LE.B(JT+1)) GO TO 90 PR509500 JT=JT+1 PR509600 GO TO 80 PR509700 C PR509800 90 U=(X(I)-B(JT))/H PR509900 Q(1)=P1(1.-U) PR510000 Q(2)=P2(1.-U) PR510100 Q(3)=P2(U) PR510200 Q(4)=P1(U) PR510300 YFIT=ZERO PR510400 DO 100 L=1,4 PR510500 IC=JT-1+L PR510600 100 YFIT=YFIT+C(IC)*Q(L) PR510700 R=Y(I)-YFIT PR510800 WRITE (6,170) I,X(I),Y(I),YFIT,R PR510900 110 CONTINUE PR511000 C PR511100 C COMPUTE RESIDUAL VECTOR NORM. PR511200 C PR511300 IF (M.LE.NC) GO TO 150 PR511400 SIGSQ=(RNORM**2)/FLOAT(M-NC) PR511500 SIGFAC=SQRT(SIGSQ) PR511600 WRITE (6,220) SIGFAC PR511700 WRITE (6,200) PR511800 C PR511900 C COMPUTE AND PRINT COLS. OF COVARIANCE. PR512000 C PR512100 DO 140 J=1,NC PR512200 DO 120 I=1,NC PR512300 120 COV(I)=ZERO PR512400 COV(J)=1. PR512500 CALL BNDSOL (2,G,MDG,NBAND,IP,IR,COV,NC,RDUMMY) PR512600 CALL BNDSOL (3,G,MDG,NBAND,IP,IR,COV,NC,RDUMMY) PR512700 C PR512800 C COMPUTE THE JTH COL. OF THE COVARIANCE MATRIX. PR512900 DO 130 I=1,NC PR513000 130 COV(I)=COV(I)*SIGSQ PR513100 140 WRITE (6,190) (L,J,COV(L),L=1,NC) PR513200 150 CONTINUE PR513300 STOP PR513400 160 FORMAT (4H0 I,8X,1HX,10X,1HY,6X,4HYFIT,4X,8HR=Y-YFIT/1X) PR513500 170 FORMAT (1X,I3,4X,F6.0,4X,F6.2,4X,F6.2,4X,F8.4) PR513600 180 FORMAT (4H0C =,10F10.5/(4X,10F10.5)) PR513700 190 FORMAT (3(02X,2I4,E15.7)) PR513800 200 FORMAT (46H0COVARIANCE MATRIX OF THE SPLINE COEFFICIENTS.) PR513900 210 FORMAT (9H0RNORM =,E15.8) PR514000 220 FORMAT (9H0SIGFAC =,E15.8) PR514100 230 FORMAT (50H1PROG5. EXECUTE A SEQUENCE OF CUBIC SPLINE FITS,27H PR514200 1TO A DISCRETE SET OF DATA.) PR514300 240 FORMAT (1X////,11H0NEW CASE..,/47H0NUMBER OF BREAKPOINTS, INCLUDINPR514400 1G ENDPOINTS, IS,I5/14H0BREAKPOINTS =,10F10.5,/(14X,10F10.5)) PR514500 END PR514600