C     PROG1                                                             PR100100
C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 PR100200
C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974PR100300
C        DEMONSTRATE ALGORITHMS HFT AND HS1 FOR SOLVING LEAST SQUARES   PR100400
C     PROBLEMS AND ALGORITHM COV FOR OMPUTING THE ASSOCIATED COVARIANCE PR100500
C     MATRICES.                                                         PR100600
C                                                                       PR100700
      DIMENSION A(8,8),H(8),B(8)                                        PR100800
      REAL GEN,ANOISE                                                   PR100900
      DOUBLE PRECISION SM                                               PR101000
      DATA MDA/8/                                                       PR101100
C                                                                       PR101200
C Set math to execute single precision in 23-bit mantissa format.               
      CALL MPBRQQ
C
           DO 180 NOISE=1,2                                             PR101300
           ANOISE=0.                                                    PR101400
           IF (NOISE.EQ.2) ANOISE=1.E-4                                 PR101500
           WRITE (6,230)                                                PR101600
           WRITE (6,240) ANOISE                                         PR101700
C     INITIALIZE THE DATA GENERATION FUNCTION                           PR101800
C    ..                                                                 PR101900
           DUMMY=GEN(-1.)                                               PR102000
                DO 180 MN1=1,6,5                                        PR102100
                MN2=MN1+2                                               PR102200
                     DO 180 M=MN1,MN2                                   PR102300
                     DO 180 N=MN1,MN2                                   PR102400
                     NP1=N+1                                            PR102500
                     WRITE (6,250) M,N                                  PR102600
C     GENERATE DATA                                                     PR102700
C    ..                                                                 PR102800
                          DO 10 I=1,M                                   PR102900
                               DO 10 J=1,N                              PR103000
   10                          A(I,J)=GEN(ANOISE)                       PR103100
                          DO 20 I=1,M                                   PR103200
   20                     B(I)=GEN(ANOISE)                              PR103300
                     IF(M .LT. N) GO TO 180                             PR103400
C                                                                       PR103500
C     ******  BEGIN ALGORITHM HFT  ******                               PR103600
C    ..                                                                 PR103700
                          DO 30 J=1,N                                   PR103800
   30              CALL H12 (1,J,J+1,M,A(1,J),1,H(J),A(1,J+1),1,MDA,N-J)PR103900
C    ..                                                                 PR104000
C     THE ALGORITHM 'HFT' IS COMPLETED.                                 PR104100
C                                                                       PR104200
C     ******  BEGIN ALGORITHM HS1  ******                               PR104300
C     APPLY THE TRANSFORMATIONS  Q(N)...Q(1)=Q TO B                     PR104400
C     REPLACING THE PREVIOUS CONTENTS OF THE ARRAY, B .                 PR104500
C    ..                                                                 PR104600
                          DO 40 J=1,N                                   PR104700
   40                         CALL H12 (2,J,J+1,M,A(1,J),1,H(J),B,1,1,1)PR104800
C     SOLVE THE TRIANGULAR SYSTEM FOR THE SOLUTION X.                   PR104900
C     STORE X IN THE ARRAY, B .                                         PR105000
C    ..                                                                 PR105100
                          DO 80 K=1,N                                   PR105200
                          I=NP1-K                                       PR105300
                          SM=0.D0                                       PR105400
                          IF (I.EQ.N) GO TO 60                          PR105500
                          IP1=I+1                                       PR105600
                               DO 50 J=IP1,N                            PR105700
   50                          SM=SM+A(I,J)*DBLE(B(J))                  PR105800
   60                     IF (A(I,I)) 80,70,80                          PR105900
   70                     WRITE (6,260)                                 PR106000
                          GO TO 180                                     PR106100
   80                     B(I)=(B(I)-SM)/A(I,I)                         PR106200
C      COMPUTE LENGTH OF RESIDUAL VECTOR.                               PR106300
C     ..                                                                PR106400
                     SRSMSQ=0.                                          PR106500
                     IF (N.EQ.M) GO TO 100                              PR106600
                     MMN=M-N                                            PR106700
                          DO 90 J=1,MMN                                 PR106800
                          NPJ=N+J                                       PR106900
   90                     SRSMSQ=SRSMSQ+B(NPJ)**2                       PR107000
                     SRSMSQ=SQRT(SRSMSQ)                                PR107100
C     ******  BEGIN ALGORITHM  COV  ******                              PR107200
C     COMPUTE UNSCALED COVARIANCE MATRIX   ((A**T)*A)**(-1)             PR107300
C    ..                                                                 PR107400
  100                     DO 110 J=1,N                                  PR107500
  110                     A(J,J)=1./A(J,J)                              PR107600
                     IF (N.EQ.1) GO TO 140                              PR107700
                     NM1=N-1                                            PR107800
                          DO 130 I=1,NM1                                PR107900
                          IP1=I+1                                       PR108000
                               DO 130 J=IP1,N                           PR108100
                               JM1=J-1                                  PR108200
                               SM=0.D0                                  PR108300
                                    DO 120 L=I,JM1                      PR108400
  120                               SM=SM+A(I,L)*DBLE(A(L,J))           PR108500
  130                          A(I,J)=-SM*A(J,J)                        PR108600
C    ..                                                                 PR108700
C     THE UPPER TRIANGLE OF A HAS BEEN INVERTED                         PR108800
C     UPON ITSELF.                                                      PR108900
  140                     DO 160 I=1,N                                  PR109000
                               DO 160 J=I,N                             PR109100
                               SM=0.D0                                  PR109200
                                    DO 150 L=J,N                        PR109300
  150                               SM=SM+A(I,L)*DBLE(A(J,L))           PR109400
  160                          A(I,J)=SM                                PR109500
C    ..                                                                 PR109600
C     THE UPPER TRIANGULAR PART OF THE                                  PR109700
C     SYMMETRIC MATRIX (A**T*A)**(-1) HAS                               PR109800
C     REPLACED THE UPPER TRIANGULAR PART OF                             PR109900
C     THE A ARRAY.                                                      PR110000
                     WRITE (6,200) (I,B(I),I=1,N)                       PR110100
                     WRITE (6,190) SRSMSQ                               PR110200
                     WRITE (6,210)                                      PR110300
                          DO 170 I=1,N                                  PR110400
  170                     WRITE (6,220) (I,J,A(I,J),J=I,N)              PR110500
  180                CONTINUE                                           PR110600
      STOP                                                              PR110700
  190 FORMAT (1H0,8X,17HRESIDUAL LENGTH =,E12.4)                        PR110800
  200 FORMAT (1H0,8X,34HESTIMATED PARAMETERS,  X=A**(+)*B,,22H COMPUTED PR110900
     1BY 'HFT,HS1'//(9X,I6,E16.8,I6,E16.8,I6,E16.8,I6,E16.8,I6,E16.8))  PR111000
  210 FORMAT (1H0,8X,31HCOVARIANCE MATRIX (UNSCALED) OF,22H ESTIMATED PAPR111100
     1RAMETERS.,19H COMPUTED BY 'COV'./1X)                              PR111200
  220 FORMAT (9X,2I3,E16.8,2I3,E16.8,2I3,E16.8,2I3,E16.8,2I3,E16.8)     PR111300
  230 FORMAT (52H1 PROG1.    THIS PROGRAM DEMONSTRATES THE ALGORITHMS,19PR111400
     1H HFT, HS1, AND COV.//40H CAUTION.. 'PROG1' DOES NO CHECKING FOR ,PR111500
     252HNEAR RANK DEFICIENT MATRICES.  RESULTS IN SUCH CASES,20H MAY BEPR111600
     3 MEANINGLESS.,/34H SUCH CASES ARE HANDLED BY 'PROG2',11H OR 'PROG3PR111700
     4')                                                                PR111800
  240 FORMAT (1H0,54HTHE RELATIVE NOISE LEVEL OF THE GENERATED DATA WILLPR111900
     1 BE,E16.4)                                                        PR112000
  250 FORMAT (1H0////9H0   M   N/1X,2I4)                                PR112100
  260 FORMAT (1H0,8X,36H******  TERMINATING THIS CASE DUE TO,37H A DIVISPR112200
     1OR BEING EXACTLY ZERO. ******)                                    PR112300
      END                                                               PR112400
