C     PROG2                                                             PR200100
C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 PR200200
C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974PR200300
C        DEMONSTRATE ALGORITHM  HFTI  FOR SOLVING LEAST SQUARES PROBLEMSPR200400
C     AND ALGORITHM  COV  FOR COMPUTING THE ASSOCIATED UNSCALED         PR200500
C     COVARIANCE MATRIX.                                                PR200600
C                                                                       PR200700
      DIMENSION A(8,8),H(8),B(8),G(8)                                   PR200800
      REAL GEN,ANOISE                                                   PR200900
      INTEGER IP(8)                                                     PR201000
      DOUBLE PRECISION SM                                               PR201100
      DATA MDA /8/                                                      PR201200
C                                                                       PR201300
C Set math to execute single precision in 23-bit mantissa format.               
      CALL MPBRQQ
C
          DO 180 NOISE=1,2                                              PR201400
          ANORM=500.                                                    PR201500
          ANOISE=0.                                                     PR201600
          TAU=.5                                                        PR201700
          IF (NOISE.EQ.1) GO TO 10                                      PR201800
          ANOISE=1.E-4                                                  PR201900
          TAU=ANORM*ANOISE*10.                                          PR202000
   10     CONTINUE                                                      PR202100
C     INITIALIZE THE DATA GENERATION FUNCTION                           PR202200
C    ..                                                                 PR202300
          DUMMY=GEN(-1.)                                                PR202400
          WRITE (6,230)                                                 PR202500
          WRITE (6,240) ANOISE,ANORM,TAU                                PR202600
C                                                                       PR202700
              DO 180 MN1=1,6,5                                          PR202800
              MN2=MN1+2                                                 PR202900
                  DO 180 M=MN1,MN2                                      PR203000
                      DO 180 N=MN1,MN2                                  PR203100
                      WRITE (6,250) M,N                                 PR203200
C     GENERATE DATA                                                     PR203300
C    ..                                                                 PR203400
                          DO 20 I=1,M                                   PR203500
                              DO 20 J=1,N                               PR203600
   20                         A(I,J)=GEN(ANOISE)                        PR203700
                          DO 30 I=1,M                                   PR203800
   30                     B(I)=GEN(ANOISE)                              PR203900
C                                                                       PR204000
C     ****** CALL HFTI   ******                                         PR204100
C                                                                       PR204200
                      CALL HFTI(A,MDA,M,N,B,1,1,TAU,KRANK,SRSMSQ,H,G,IP)PR204300
C                                                                       PR204400
C                                                                       PR204500
                      WRITE (6,260) KRANK                               PR204600
                      WRITE (6,200) (I,B(I),I=1,N)                      PR204700
                      WRITE (6,190) SRSMSQ                              PR204800
                      IF (KRANK.LT.N) GO TO 180                         PR204900
C     ******  ALGORITHM COV BIGINS HERE  ******                         PR205000
C    ..                                                                 PR205100
                          DO 40 J=1,N                                   PR205200
   40                     A(J,J)=1./A(J,J)                              PR205300
                      IF (N.EQ.1) GO TO 70                              PR205400
                      NM1=N-1                                           PR205500
                          DO 60 I=1,NM1                                 PR205600
                          IP1=I+1                                       PR205700
                              DO 60 J=IP1,N                             PR205800
                              JM1=J-1                                   PR205900
                              SM=0.D0                                   PR206000
                                  DO 50 L=I,JM1                         PR206100
   50                             SM=SM+A(I,L)*DBLE(A(L,J))             PR206200
   60                         A(I,J)=-SM*A(J,J)                         PR206300
C    ..                                                                 PR206400
C     THE UPPER TRIANGLE OF A HAS BEEN INVERTED                         PR206500
C     UPON ITSELF.                                                      PR206600
   70                     DO 90 I=1,N                                   PR206700
                              DO 90 J=I,N                               PR206800
                              SM=0.D0                                   PR206900
                                  DO 80 L=J,N                           PR207000
   80                             SM=SM+A(I,L)*DBLE(A(J,L))             PR207100
   90                         A(I,J)=SM                                 PR207200
                      IF (N.LT.2) GO TO 160                             PR207300
                          DO 150 II=2,N                                 PR207400
                          I=N+1-II                                      PR207500
                          IF (IP(I).EQ.I) GO TO 150                     PR207600
                          K=IP(I)                                       PR207700
                          TMP=A(I,I)                                    PR207800
                          A(I,I)=A(K,K)                                 PR207900
                          A(K,K)=TMP                                    PR208000
                          IF (I.EQ.1) GO TO 110                         PR208100
                              DO 100 L=2,I                              PR208200
                              TMP=A(L-1,I)                              PR208300
                              A(L-1,I)=A(L-1,K)                         PR208400
  100                         A(L-1,K)=TMP                              PR208500
  110                     IP1=I+1                                       PR208600
                          KM1=K-1                                       PR208700
                          IF (IP1.GT.KM1) GO TO 130                     PR208800
                              DO 120 L=IP1,KM1                          PR208900
                              TMP=A(I,L)                                PR209000
                              A(I,L)=A(L,K)                             PR209100
  120                         A(L,K)=TMP                                PR209200
  130                     IF (K.EQ.N) GO TO 150                         PR209300
                          KP1=K+1                                       PR209400
                              DO 140 L=KP1,N                            PR209500
                              TMP=A(I,L)                                PR209600
                              A(I,L)=A(K,L)                             PR209700
  140                         A(K,L)=TMP                                PR209800
  150                     CONTINUE                                      PR209900
  160                 CONTINUE                                          PR210000
C    ..                                                                 PR210100
C     COVARIANCE HAS BEEN COMPUTED AND REPERMUTED.                      PR210200
C     THE UPPER TRIANGULAR PART OF THE                                  PR210300
C     SYMMETRIC MATRIX (A**T*A)**(-1) HAS                               PR210400
C     REPLACED THE UPPER TRIANGULAR PART OF                             PR210500
C     THE A ARRAY.                                                      PR210600
                      WRITE (6,210)                                     PR210700
                          DO 170 I=1,N                                  PR210800
  170                     WRITE (6,220) (I,J,A(I,J),J=I,N)              PR210900
  180                 CONTINUE                                          PR211000
      STOP                                                              PR211100
  190 FORMAT (1H0,8X,17HRESIDUAL LENGTH =,E12.4)                        PR211200
  200 FORMAT (1H0,8X,34HESTIMATED PARAMETERS,  X=A**(+)*B,,22H COMPUTED PR211300
     1BY 'HFTI'   //(9X,I6,E16.8,I6,E16.8,I6,E16.8,I6,E16.8,I6,E16.8))  PR211400
  210 FORMAT (1H0,8X,31HCOVARIANCE MATRIX (UNSCALED) OF,22H ESTIMATED PAPR211500
     1RAMETERS.,19H COMPUTED BY 'COV'./1X)                              PR211600
  220 FORMAT (9X,2I3,E16.8,2I3,E16.8,2I3,E16.8,2I3,E16.8,2I3,E16.8)     PR211700
  230 FORMAT (52H1 PROG2.     THIS PROGRAM DEMONSTATES THE ALGORITHMS,16PR211800
     1H HFTI  AND  COV.)                                                PR211900
  240 FORMAT (1H0,54HTHE RELATIVE NOISE LEVEL OF THE GENERATED DATA WILLPR212000
     1 BE,E16.4/33H0THE MATRIX NORM IS APPROXIMATELY,E12.4/43H0THE ABSOLPR212100
     2UTE PSEUDORANK TOLERANCE, TAU, IS,E12.4)                          PR212200
  250 FORMAT (1H0////9H0   M   N/1X,2I4)                                PR212300
  260 FORMAT (1H0,8X,12HPSEUDORANK =,I4)                                PR212400
      END                                                               PR212500
C     PROG3                                                             PR300100
