home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / FORTRAN / LLSQ.ZIP / PROG2.FOR < prev    next >
Encoding:
Text File  |  1984-02-24  |  10.2 KB  |  130 lines
  1. C     PROG2                                                             PR200100
  2. C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 PR200200
  3. C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974PR200300
  4. C        DEMONSTRATE ALGORITHM  HFTI  FOR SOLVING LEAST SQUARES PROBLEMSPR200400
  5. C     AND ALGORITHM  COV  FOR COMPUTING THE ASSOCIATED UNSCALED         PR200500
  6. C     COVARIANCE MATRIX.                                                PR200600
  7. C                                                                       PR200700
  8.       DIMENSION A(8,8),H(8),B(8),G(8)                                   PR200800
  9.       REAL GEN,ANOISE                                                   PR200900
  10.       INTEGER IP(8)                                                     PR201000
  11.       DOUBLE PRECISION SM                                               PR201100
  12.       DATA MDA /8/                                                      PR201200
  13. C                                                                       PR201300
  14. C Set math to execute single precision in 23-bit mantissa format.               
  15.       CALL MPBRQQ
  16. C
  17.           DO 180 NOISE=1,2                                              PR201400
  18.           ANORM=500.                                                    PR201500
  19.           ANOISE=0.                                                     PR201600
  20.           TAU=.5                                                        PR201700
  21.           IF (NOISE.EQ.1) GO TO 10                                      PR201800
  22.           ANOISE=1.E-4                                                  PR201900
  23.           TAU=ANORM*ANOISE*10.                                          PR202000
  24.    10     CONTINUE                                                      PR202100
  25. C     INITIALIZE THE DATA GENERATION FUNCTION                           PR202200
  26. C    ..                                                                 PR202300
  27.           DUMMY=GEN(-1.)                                                PR202400
  28.           WRITE (6,230)                                                 PR202500
  29.           WRITE (6,240) ANOISE,ANORM,TAU                                PR202600
  30. C                                                                       PR202700
  31.               DO 180 MN1=1,6,5                                          PR202800
  32.               MN2=MN1+2                                                 PR202900
  33.                   DO 180 M=MN1,MN2                                      PR203000
  34.                       DO 180 N=MN1,MN2                                  PR203100
  35.                       WRITE (6,250) M,N                                 PR203200
  36. C     GENERATE DATA                                                     PR203300
  37. C    ..                                                                 PR203400
  38.                           DO 20 I=1,M                                   PR203500
  39.                               DO 20 J=1,N                               PR203600
  40.    20                         A(I,J)=GEN(ANOISE)                        PR203700
  41.                           DO 30 I=1,M                                   PR203800
  42.    30                     B(I)=GEN(ANOISE)                              PR203900
  43. C                                                                       PR204000
  44. C     ****** CALL HFTI   ******                                         PR204100
  45. C                                                                       PR204200
  46.                       CALL HFTI(A,MDA,M,N,B,1,1,TAU,KRANK,SRSMSQ,H,G,IP)PR204300
  47. C                                                                       PR204400
  48. C                                                                       PR204500
  49.                       WRITE (6,260) KRANK                               PR204600
  50.                       WRITE (6,200) (I,B(I),I=1,N)                      PR204700
  51.                       WRITE (6,190) SRSMSQ                              PR204800
  52.                       IF (KRANK.LT.N) GO TO 180                         PR204900
  53. C     ******  ALGORITHM COV BIGINS HERE  ******                         PR205000
  54. C    ..                                                                 PR205100
  55.                           DO 40 J=1,N                                   PR205200
  56.    40                     A(J,J)=1./A(J,J)                              PR205300
  57.                       IF (N.EQ.1) GO TO 70                              PR205400
  58.                       NM1=N-1                                           PR205500
  59.                           DO 60 I=1,NM1                                 PR205600
  60.                           IP1=I+1                                       PR205700
  61.                               DO 60 J=IP1,N                             PR205800
  62.                               JM1=J-1                                   PR205900
  63.                               SM=0.D0                                   PR206000
  64.                                   DO 50 L=I,JM1                         PR206100
  65.    50                             SM=SM+A(I,L)*DBLE(A(L,J))             PR206200
  66.    60                         A(I,J)=-SM*A(J,J)                         PR206300
  67. C    ..                                                                 PR206400
  68. C     THE UPPER TRIANGLE OF A HAS BEEN INVERTED                         PR206500
  69. C     UPON ITSELF.                                                      PR206600
  70.    70                     DO 90 I=1,N                                   PR206700
  71.                               DO 90 J=I,N                               PR206800
  72.                               SM=0.D0                                   PR206900
  73.                                   DO 80 L=J,N                           PR207000
  74.    80                             SM=SM+A(I,L)*DBLE(A(J,L))             PR207100
  75.    90                         A(I,J)=SM                                 PR207200
  76.                       IF (N.LT.2) GO TO 160                             PR207300
  77.                           DO 150 II=2,N                                 PR207400
  78.                           I=N+1-II                                      PR207500
  79.                           IF (IP(I).EQ.I) GO TO 150                     PR207600
  80.                           K=IP(I)                                       PR207700
  81.                           TMP=A(I,I)                                    PR207800
  82.                           A(I,I)=A(K,K)                                 PR207900
  83.                           A(K,K)=TMP                                    PR208000
  84.                           IF (I.EQ.1) GO TO 110                         PR208100
  85.                               DO 100 L=2,I                              PR208200
  86.                               TMP=A(L-1,I)                              PR208300
  87.                               A(L-1,I)=A(L-1,K)                         PR208400
  88.   100                         A(L-1,K)=TMP                              PR208500
  89.   110                     IP1=I+1                                       PR208600
  90.                           KM1=K-1                                       PR208700
  91.                           IF (IP1.GT.KM1) GO TO 130                     PR208800
  92.                               DO 120 L=IP1,KM1                          PR208900
  93.                               TMP=A(I,L)                                PR209000
  94.                               A(I,L)=A(L,K)                             PR209100
  95.   120                         A(L,K)=TMP                                PR209200
  96.   130                     IF (K.EQ.N) GO TO 150                         PR209300
  97.                           KP1=K+1                                       PR209400
  98.                               DO 140 L=KP1,N                            PR209500
  99.                               TMP=A(I,L)                                PR209600
  100.                               A(I,L)=A(K,L)                             PR209700
  101.   140                         A(K,L)=TMP                                PR209800
  102.   150                     CONTINUE                                      PR209900
  103.   160                 CONTINUE                                          PR210000
  104. C    ..                                                                 PR210100
  105. C     COVARIANCE HAS BEEN COMPUTED AND REPERMUTED.                      PR210200
  106. C     THE UPPER TRIANGULAR PART OF THE                                  PR210300
  107. C     SYMMETRIC MATRIX (A**T*A)**(-1) HAS                               PR210400
  108. C     REPLACED THE UPPER TRIANGULAR PART OF                             PR210500
  109. C     THE A ARRAY.                                                      PR210600
  110.                       WRITE (6,210)                                     PR210700
  111.                           DO 170 I=1,N                                  PR210800
  112.   170                     WRITE (6,220) (I,J,A(I,J),J=I,N)              PR210900
  113.   180                 CONTINUE                                          PR211000
  114.       STOP                                                              PR211100
  115.   190 FORMAT (1H0,8X,17HRESIDUAL LENGTH =,E12.4)                        PR211200
  116.   200 FORMAT (1H0,8X,34HESTIMATED PARAMETERS,  X=A**(+)*B,,22H COMPUTED PR211300
  117.      1BY 'HFTI'   //(9X,I6,E16.8,I6,E16.8,I6,E16.8,I6,E16.8,I6,E16.8))  PR211400
  118.   210 FORMAT (1H0,8X,31HCOVARIANCE MATRIX (UNSCALED) OF,22H ESTIMATED PAPR211500
  119.      1RAMETERS.,19H COMPUTED BY 'COV'./1X)                              PR211600
  120.   220 FORMAT (9X,2I3,E16.8,2I3,E16.8,2I3,E16.8,2I3,E16.8,2I3,E16.8)     PR211700
  121.   230 FORMAT (52H1 PROG2.     THIS PROGRAM DEMONSTATES THE ALGORITHMS,16PR211800
  122.      1H HFTI  AND  COV.)                                                PR211900
  123.   240 FORMAT (1H0,54HTHE RELATIVE NOISE LEVEL OF THE GENERATED DATA WILLPR212000
  124.      1 BE,E16.4/33H0THE MATRIX NORM IS APPROXIMATELY,E12.4/43H0THE ABSOLPR212100
  125.      2UTE PSEUDORANK TOLERANCE, TAU, IS,E12.4)                          PR212200
  126.   250 FORMAT (1H0////9H0   M   N/1X,2I4)                                PR212300
  127.   260 FORMAT (1H0,8X,12HPSEUDORANK =,I4)                                PR212400
  128.       END                                                               PR212500
  129. C     PROG3                                                             PR300100
  130.