home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / FORTRAN / LLSQ.ZIP / PROG5.FOR < prev    next >
Encoding:
Text File  |  1984-03-01  |  11.7 KB  |  148 lines
  1. C     PROG5                                                             PR500100
  2. C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 PR500200
  3. C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974PR500300
  4. C          EXAMPLE OF THE USE OF SUBROUTINES BNDACC AND BNDSOL TO SOLVE PR500400
  5. C     SEQUENTIALLY THE BANDED LEAST SQUARES PROBLEM WHICH ARISES IN     PR500500
  6. C     SPLINE CURVE FITTING.                                             PR500600
  7. C                                                                       PR500700
  8.       DIMENSION X(12),Y(12),B(10),G(12,5),C(12),Q(4),COV(12)            PR500800
  9. C                                                                       PR500900
  10. C        DEFINE FUNCTIONS P1 AND P2 TO BE USED IN CONSTRUCTING          PR501000
  11. C        CUBIC SPLINES OVER UNIFORMLY SPACED BREAKPOINTS.               PR501100
  12. C                                                                       PR501200
  13.       P1(T)=.25*T**2*T                                                  PR501300
  14.       P2(T)=-(1.-T)**2*(1.+T)*.75+1.                                    PR501400
  15.       CALL MPBRQQ
  16.       ZERO=0.                                                           PR501500
  17. C                                                                       PR501600
  18.       WRITE (6,230)                                                     PR501700
  19.       MDG=12                                                            PR501800
  20.       NBAND=4                                                           PR501900
  21.       M=12                                                              PR502000
  22. C                                  SET ORDINATES OF DATA                PR502100
  23.       Y(1)=2.2                                                          PR502200
  24.       Y(2)=4.0                                                          PR502300
  25.       Y(3)=5.0                                                          PR502400
  26.       Y(4)=4.6                                                          PR502500
  27.       Y(5)=2.8                                                          PR502600
  28.       Y(6)=2.7                                                          PR502700
  29.       Y(7)=3.8                                                          PR502800
  30.       Y(8)=5.1                                                          PR502900
  31.       Y(9)=6.1                                                          PR503000
  32.       Y(10)=6.3                                                         PR503100
  33.       Y(11)=5.0                                                         PR503200
  34.       Y(12)=2.0                                                         PR503300
  35. C                                  SET ABCISSAS OF DATA                 PR503400
  36.            DO 10 I=1,M                                                  PR503500
  37.    10      X(I)=2*I                                                     PR503600
  38. C                                                                       PR503700
  39. C     BEGIN LOOP THRU CASES USING INCREASING NOS OF BREAKPOINTS.        PR503800
  40. C                                                                       PR503900
  41.            DO 150 NBP=5,10                                              PR504000
  42.            NC=NBP+2                                                     PR504100
  43. C                                  SET BREAKPOINTS                      PR504200
  44.            B(1)=X(1)                                                    PR504300
  45.            B(NBP)=X(M)                                                  PR504400
  46.            H=(B(NBP)-B(1))/FLOAT(NBP-1)                                 PR504500
  47.            IF (NBP.LE.2) GO TO 30                                       PR504600
  48.                 DO 20 I=3,NBP                                           PR504700
  49.    20           B(I-1)=B(I-2)+H                                         PR504800
  50.    30      CONTINUE                                                     PR504900
  51.            WRITE (6,240) NBP,(B(I),I=1,NBP)                             PR505000
  52. C                                                                       PR505100
  53. C     INITIALIZE IR AND IP BEFORE FIRST CALL TO BNDACC.                 PR505200
  54. C                                                                       PR505300
  55.            IR=1                                                         PR505400
  56.            IP=1                                                         PR505500
  57.            I=1                                                          PR505600
  58.            JT=1                                                         PR505700
  59.    40      MT=0                                                         PR505800
  60.    50      CONTINUE                                                     PR505900
  61.            IF (X(I).GT.B(JT+1)) GO TO 60                                PR506000
  62. C                                                                       PR506100
  63. C                        SET ROW  FOR ITH DATA POINT                    PR506200
  64. C                                                                       PR506300
  65.            U=(X(I)-B(JT))/H                                             PR506400
  66.            IG=IR+MT                                                     PR506500
  67.            G(IG,1)=P1(1.-U)                                             PR506600
  68.            G(IG,2)=P2(1.-U)                                             PR506700
  69.            G(IG,3)=P2(U)                                                PR506800
  70.            G(IG,4)=P1(U)                                                PR506900
  71.            G(IG,5)=Y(I)                                                 PR507000
  72.            MT=MT+1                                                      PR507100
  73.            IF (I.EQ.M) GO TO 60                                         PR507200
  74.            I=I+1                                                        PR507300
  75.            GO TO 50                                                     PR507400
  76. C                                                                       PR507500
  77. C                   SEND BLOCK OF DATA TO PROCESSOR                     PR507600
  78. C                                                                       PR507700
  79.    60      CONTINUE                                                     PR507800
  80.            CALL BNDACC (G,MDG,NBAND,IP,IR,MT,JT)                        PR507900
  81.            IF (I.EQ.M) GO TO 70                                         PR508000
  82.            JT=JT+1                                                      PR508100
  83.            GO TO 40                                                     PR508200
  84. C                   COMPUTE SOLUTION C()                                PR508300
  85.    70      CONTINUE                                                     PR508400
  86.            CALL BNDSOL (1,G,MDG,NBAND,IP,IR,C,NC,RNORM)                 PR508500
  87. C                  WRITE SOLUTION COEFFICIENTS                          PR508600
  88.            WRITE (6,180) (C(L),L=1,NC)                                  PR508700
  89.            WRITE (6,210) RNORM                                          PR508800
  90. C                                                                       PR508900
  91. C              COMPUTE AND PRINT X,Y,YFIT,R=Y-YFIT                      PR509000
  92. C                                                                       PR509100
  93.            WRITE (6,160)                                                PR509200
  94.            JT=1                                                         PR509300
  95.                 DO 110 I=1,M                                            PR509400
  96.    80           IF (X(I).LE.B(JT+1)) GO TO 90                           PR509500
  97.                 JT=JT+1                                                 PR509600
  98.                 GO TO 80                                                PR509700
  99. C                                                                       PR509800
  100.    90           U=(X(I)-B(JT))/H                                        PR509900
  101.                 Q(1)=P1(1.-U)                                           PR510000
  102.                 Q(2)=P2(1.-U)                                           PR510100
  103.                 Q(3)=P2(U)                                              PR510200
  104.                 Q(4)=P1(U)                                              PR510300
  105.                 YFIT=ZERO                                               PR510400
  106.                      DO 100 L=1,4                                       PR510500
  107.                      IC=JT-1+L                                          PR510600
  108.   100                YFIT=YFIT+C(IC)*Q(L)                               PR510700
  109.                 R=Y(I)-YFIT                                             PR510800
  110.                 WRITE (6,170) I,X(I),Y(I),YFIT,R                        PR510900
  111.   110           CONTINUE                                                PR511000
  112. C                                                                       PR511100
  113. C     COMPUTE RESIDUAL VECTOR NORM.                                     PR511200
  114. C                                                                       PR511300
  115.            IF (M.LE.NC) GO TO 150                                       PR511400
  116.            SIGSQ=(RNORM**2)/FLOAT(M-NC)                                 PR511500
  117.            SIGFAC=SQRT(SIGSQ)                                           PR511600
  118.            WRITE (6,220) SIGFAC                                         PR511700
  119.            WRITE (6,200)                                                PR511800
  120. C                                                                       PR511900
  121. C     COMPUTE AND PRINT COLS. OF COVARIANCE.                            PR512000
  122. C                                                                       PR512100
  123.                 DO 140 J=1,NC                                           PR512200
  124.                      DO 120 I=1,NC                                      PR512300
  125.   120                COV(I)=ZERO                                        PR512400
  126.                 COV(J)=1.                                               PR512500
  127.                 CALL BNDSOL (2,G,MDG,NBAND,IP,IR,COV,NC,RDUMMY)         PR512600
  128.                 CALL BNDSOL (3,G,MDG,NBAND,IP,IR,COV,NC,RDUMMY)         PR512700
  129. C                                                                       PR512800
  130. C    COMPUTE THE JTH COL. OF THE COVARIANCE MATRIX.                     PR512900
  131.                      DO 130 I=1,NC                                      PR513000
  132.   130                COV(I)=COV(I)*SIGSQ                                PR513100
  133.   140           WRITE (6,190) (L,J,COV(L),L=1,NC)                       PR513200
  134.   150      CONTINUE                                                     PR513300
  135.       STOP                                                              PR513400
  136.   160 FORMAT (4H0  I,8X,1HX,10X,1HY,6X,4HYFIT,4X,8HR=Y-YFIT/1X)         PR513500
  137.   170 FORMAT (1X,I3,4X,F6.0,4X,F6.2,4X,F6.2,4X,F8.4)                    PR513600
  138.   180 FORMAT (4H0C =,10F10.5/(4X,10F10.5))                              PR513700
  139.   190 FORMAT (3(02X,2I4,E15.7))                                         PR513800
  140.   200 FORMAT (46H0COVARIANCE MATRIX OF THE SPLINE COEFFICIENTS.)        PR513900
  141.   210 FORMAT (9H0RNORM  =,E15.8)                                        PR514000
  142.   220 FORMAT (9H0SIGFAC =,E15.8)                                        PR514100
  143.   230 FORMAT (50H1PROG5.    EXECUTE A SEQUENCE OF CUBIC SPLINE FITS,27H PR514200
  144.      1TO A DISCRETE SET OF DATA.)                                       PR514300
  145.   240 FORMAT (1X////,11H0NEW CASE..,/47H0NUMBER OF BREAKPOINTS, INCLUDINPR514400
  146.      1G ENDPOINTS, IS,I5/14H0BREAKPOINTS =,10F10.5,/(14X,10F10.5))      PR514500
  147.       END                                                               PR514600
  148.