home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / FORTRAN / LLSQ.ZIP / PROG1.FOR < prev    next >
Encoding:
Text File  |  1984-02-23  |  10.0 KB  |  128 lines
  1. C     PROG1                                                             PR100100
  2. C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 PR100200
  3. C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974PR100300
  4. C        DEMONSTRATE ALGORITHMS HFT AND HS1 FOR SOLVING LEAST SQUARES   PR100400
  5. C     PROBLEMS AND ALGORITHM COV FOR OMPUTING THE ASSOCIATED COVARIANCE PR100500
  6. C     MATRICES.                                                         PR100600
  7. C                                                                       PR100700
  8.       DIMENSION A(8,8),H(8),B(8)                                        PR100800
  9.       REAL GEN,ANOISE                                                   PR100900
  10.       DOUBLE PRECISION SM                                               PR101000
  11.       DATA MDA/8/                                                       PR101100
  12. C                                                                       PR101200
  13. C Set math to execute single precision in 23-bit mantissa format.               
  14.       CALL MPBRQQ
  15. C
  16.            DO 180 NOISE=1,2                                             PR101300
  17.            ANOISE=0.                                                    PR101400
  18.            IF (NOISE.EQ.2) ANOISE=1.E-4                                 PR101500
  19.            WRITE (6,230)                                                PR101600
  20.            WRITE (6,240) ANOISE                                         PR101700
  21. C     INITIALIZE THE DATA GENERATION FUNCTION                           PR101800
  22. C    ..                                                                 PR101900
  23.            DUMMY=GEN(-1.)                                               PR102000
  24.                 DO 180 MN1=1,6,5                                        PR102100
  25.                 MN2=MN1+2                                               PR102200
  26.                      DO 180 M=MN1,MN2                                   PR102300
  27.                      DO 180 N=MN1,MN2                                   PR102400
  28.                      NP1=N+1                                            PR102500
  29.                      WRITE (6,250) M,N                                  PR102600
  30. C     GENERATE DATA                                                     PR102700
  31. C    ..                                                                 PR102800
  32.                           DO 10 I=1,M                                   PR102900
  33.                                DO 10 J=1,N                              PR103000
  34.    10                          A(I,J)=GEN(ANOISE)                       PR103100
  35.                           DO 20 I=1,M                                   PR103200
  36.    20                     B(I)=GEN(ANOISE)                              PR103300
  37.                      IF(M .LT. N) GO TO 180                             PR103400
  38. C                                                                       PR103500
  39. C     ******  BEGIN ALGORITHM HFT  ******                               PR103600
  40. C    ..                                                                 PR103700
  41.                           DO 30 J=1,N                                   PR103800
  42.    30              CALL H12 (1,J,J+1,M,A(1,J),1,H(J),A(1,J+1),1,MDA,N-J)PR103900
  43. C    ..                                                                 PR104000
  44. C     THE ALGORITHM 'HFT' IS COMPLETED.                                 PR104100
  45. C                                                                       PR104200
  46. C     ******  BEGIN ALGORITHM HS1  ******                               PR104300
  47. C     APPLY THE TRANSFORMATIONS  Q(N)...Q(1)=Q TO B                     PR104400
  48. C     REPLACING THE PREVIOUS CONTENTS OF THE ARRAY, B .                 PR104500
  49. C    ..                                                                 PR104600
  50.                           DO 40 J=1,N                                   PR104700
  51.    40                         CALL H12 (2,J,J+1,M,A(1,J),1,H(J),B,1,1,1)PR104800
  52. C     SOLVE THE TRIANGULAR SYSTEM FOR THE SOLUTION X.                   PR104900
  53. C     STORE X IN THE ARRAY, B .                                         PR105000
  54. C    ..                                                                 PR105100
  55.                           DO 80 K=1,N                                   PR105200
  56.                           I=NP1-K                                       PR105300
  57.                           SM=0.D0                                       PR105400
  58.                           IF (I.EQ.N) GO TO 60                          PR105500
  59.                           IP1=I+1                                       PR105600
  60.                                DO 50 J=IP1,N                            PR105700
  61.    50                          SM=SM+A(I,J)*DBLE(B(J))                  PR105800
  62.    60                     IF (A(I,I)) 80,70,80                          PR105900
  63.    70                     WRITE (6,260)                                 PR106000
  64.                           GO TO 180                                     PR106100
  65.    80                     B(I)=(B(I)-SM)/A(I,I)                         PR106200
  66. C      COMPUTE LENGTH OF RESIDUAL VECTOR.                               PR106300
  67. C     ..                                                                PR106400
  68.                      SRSMSQ=0.                                          PR106500
  69.                      IF (N.EQ.M) GO TO 100                              PR106600
  70.                      MMN=M-N                                            PR106700
  71.                           DO 90 J=1,MMN                                 PR106800
  72.                           NPJ=N+J                                       PR106900
  73.    90                     SRSMSQ=SRSMSQ+B(NPJ)**2                       PR107000
  74.                      SRSMSQ=SQRT(SRSMSQ)                                PR107100
  75. C     ******  BEGIN ALGORITHM  COV  ******                              PR107200
  76. C     COMPUTE UNSCALED COVARIANCE MATRIX   ((A**T)*A)**(-1)             PR107300
  77. C    ..                                                                 PR107400
  78.   100                     DO 110 J=1,N                                  PR107500
  79.   110                     A(J,J)=1./A(J,J)                              PR107600
  80.                      IF (N.EQ.1) GO TO 140                              PR107700
  81.                      NM1=N-1                                            PR107800
  82.                           DO 130 I=1,NM1                                PR107900
  83.                           IP1=I+1                                       PR108000
  84.                                DO 130 J=IP1,N                           PR108100
  85.                                JM1=J-1                                  PR108200
  86.                                SM=0.D0                                  PR108300
  87.                                     DO 120 L=I,JM1                      PR108400
  88.   120                               SM=SM+A(I,L)*DBLE(A(L,J))           PR108500
  89.   130                          A(I,J)=-SM*A(J,J)                        PR108600
  90. C    ..                                                                 PR108700
  91. C     THE UPPER TRIANGLE OF A HAS BEEN INVERTED                         PR108800
  92. C     UPON ITSELF.                                                      PR108900
  93.   140                     DO 160 I=1,N                                  PR109000
  94.                                DO 160 J=I,N                             PR109100
  95.                                SM=0.D0                                  PR109200
  96.                                     DO 150 L=J,N                        PR109300
  97.   150                               SM=SM+A(I,L)*DBLE(A(J,L))           PR109400
  98.   160                          A(I,J)=SM                                PR109500
  99. C    ..                                                                 PR109600
  100. C     THE UPPER TRIANGULAR PART OF THE                                  PR109700
  101. C     SYMMETRIC MATRIX (A**T*A)**(-1) HAS                               PR109800
  102. C     REPLACED THE UPPER TRIANGULAR PART OF                             PR109900
  103. C     THE A ARRAY.                                                      PR110000
  104.                      WRITE (6,200) (I,B(I),I=1,N)                       PR110100
  105.                      WRITE (6,190) SRSMSQ                               PR110200
  106.                      WRITE (6,210)                                      PR110300
  107.                           DO 170 I=1,N                                  PR110400
  108.   170                     WRITE (6,220) (I,J,A(I,J),J=I,N)              PR110500
  109.   180                CONTINUE                                           PR110600
  110.       STOP                                                              PR110700
  111.   190 FORMAT (1H0,8X,17HRESIDUAL LENGTH =,E12.4)                        PR110800
  112.   200 FORMAT (1H0,8X,34HESTIMATED PARAMETERS,  X=A**(+)*B,,22H COMPUTED PR110900
  113.      1BY 'HFT,HS1'//(9X,I6,E16.8,I6,E16.8,I6,E16.8,I6,E16.8,I6,E16.8))  PR111000
  114.   210 FORMAT (1H0,8X,31HCOVARIANCE MATRIX (UNSCALED) OF,22H ESTIMATED PAPR111100
  115.      1RAMETERS.,19H COMPUTED BY 'COV'./1X)                              PR111200
  116.   220 FORMAT (9X,2I3,E16.8,2I3,E16.8,2I3,E16.8,2I3,E16.8,2I3,E16.8)     PR111300
  117.   230 FORMAT (52H1 PROG1.    THIS PROGRAM DEMONSTRATES THE ALGORITHMS,19PR111400
  118.      1H HFT, HS1, AND COV.//40H CAUTION.. 'PROG1' DOES NO CHECKING FOR ,PR111500
  119.      252HNEAR RANK DEFICIENT MATRICES.  RESULTS IN SUCH CASES,20H MAY BEPR111600
  120.      3 MEANINGLESS.,/34H SUCH CASES ARE HANDLED BY 'PROG2',11H OR 'PROG3PR111700
  121.      4')                                                                PR111800
  122.   240 FORMAT (1H0,54HTHE RELATIVE NOISE LEVEL OF THE GENERATED DATA WILLPR111900
  123.      1 BE,E16.4)                                                        PR112000
  124.   250 FORMAT (1H0////9H0   M   N/1X,2I4)                                PR112100
  125.   260 FORMAT (1H0,8X,36H******  TERMINATING THIS CASE DUE TO,37H A DIVISPR112200
  126.      1OR BEING EXACTLY ZERO. ******)                                    PR112300
  127.       END                                                               PR112400
  128.