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