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