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
