SUBROUTINE HFTI (A,MDA,M,N,B,MDB,NB,TAU,KRANK,RNORM,H,G,IP) HFT00100 C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 HFT00200 C TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974HFT00300 C SOLVE LEAST SQUARES PROBLEM USING ALGORITHM, HFTI. HFT00400 C HFT00500 DIMENSION A(MDA,N),B(MDB,NB),H(N),G(N),RNORM(NB) HFT00600 INTEGER IP(N) HFT00700 DOUBLE PRECISION SM,DZERO HFT00800 SZERO=0. HFT00900 DZERO=0.D0 HFT01000 FACTOR=0.001 HFT01100 C HFT01200 K=0 HFT01300 LDIAG=MIN0(M,N) HFT01400 IF (LDIAG.LE.0) GO TO 270 HFT01500 DO 80 J=1,LDIAG HFT01600 IF (J.EQ.1) GO TO 20 HFT01700 C HFT01800 C UPDATE SQUARED COLUMN LENGTHS AND FIND LMAX HFT01900 C .. HFT02000 LMAX=J HFT02100 DO 10 L=J,N HFT02200 H(L)=H(L)-A(J-1,L)**2 HFT02300 IF (H(L).GT.H(LMAX)) LMAX=L HFT02400 10 CONTINUE HFT02500 IF(DIFF(HMAX+FACTOR*H(LMAX),HMAX)) 20,20,50 HFT02600 C HFT02700 C COMPUTE SQUARED COLUMN LENGTHS AND FIND LMAX HFT02800 C .. HFT02900 20 LMAX=J HFT03000 DO 40 L=J,N HFT03100 H(L)=0. HFT03200 DO 30 I=J,M HFT03300 30 H(L)=H(L)+A(I,L)**2 HFT03400 IF (H(L).GT.H(LMAX)) LMAX=L HFT03500 40 CONTINUE HFT03600 HMAX=H(LMAX) HFT03700 C .. HFT03800 C LMAX HAS BEEN DETERMINED HFT03900 C HFT04000 C DO COLUMN INTERCHANGES IF NEEDED. HFT04100 C .. HFT04200 50 CONTINUE HFT04300 IP(J)=LMAX HFT04400 IF (IP(J).EQ.J) GO TO 70 HFT04500 DO 60 I=1,M HFT04600 TMP=A(I,J) HFT04700 A(I,J)=A(I,LMAX) HFT04800 60 A(I,LMAX)=TMP HFT04900 H(LMAX)=H(J) HFT05000 C HFT05100 C COMPUTE THE J-TH TRANSFORMATION AND APPLY IT TO A AND B. HFT05200 C .. HFT05300 70 CALL H12 (1,J,J+1,M,A(1,J),1,H(J),A(1,J+1),1,MDA,N-J) HFT05400 80 CALL H12 (2,J,J+1,M,A(1,J),1,H(J),B,1,MDB,NB) HFT05500 C HFT05600 C DETERMINE THE PSEUDORANK, K, USING THE TOLERANCE, TAU. HFT05700 C .. HFT05800 DO 90 J=1,LDIAG HFT05900 IF (ABS(A(J,J)).LE.TAU) GO TO 100 HFT06000 90 CONTINUE HFT06100 K=LDIAG HFT06200 GO TO 110 HFT06300 100 K=J-1 HFT06400 110 KP1=K+1 HFT06500 C HFT06600 C COMPUTE THE NORMS OF THE RESIDUAL VECTORS. HFT06700 C HFT06800 IF (NB.LE.0) GO TO 140 HFT06900 DO 130 JB=1,NB HFT07000 TMP=SZERO HFT07100 IF (KP1.GT.M) GO TO 130 HFT07200 DO 120 I=KP1,M HFT07300 120 TMP=TMP+B(I,JB)**2 HFT07400 130 RNORM(JB)=SQRT(TMP) HFT07500 140 CONTINUE HFT07600 C SPECIAL FOR PSEUDORANK = 0 HFT07700 IF (K.GT.0) GO TO 160 HFT07800 IF (NB.LE.0) GO TO 270 HFT07900 DO 150 JB=1,NB HFT08000 DO 150 I=1,N HFT08100 150 B(I,JB)=SZERO HFT08200 GO TO 270 HFT08300 C HFT08400 C IF THE PSEUDORANK IS LESS THAN N COMPUTE HOUSEHOLDER HFT08500 C DECOMPOSITION OF FIRST K ROWS. HFT08600 C .. HFT08700 160 IF (K.EQ.N) GO TO 180 HFT08800 DO 170 II=1,K HFT08900 I=KP1-II HFT09000 170 CALL H12 (1,I,KP1,N,A(I,1),MDA,G(I),A,MDA,1,I-1) HFT09100 180 CONTINUE HFT09200 C HFT09300 C HFT09400 IF (NB.LE.0) GO TO 270 HFT09500 DO 260 JB=1,NB HFT09600 C HFT09700 C SOLVE THE K BY K TRIANGULAR SYSTEM. HFT09800 C .. HFT09900 DO 210 L=1,K HFT10000 SM=DZERO HFT10100 I=KP1-L HFT10200 IF (I.EQ.K) GO TO 200 HFT10300 IP1=I+1 HFT10400 DO 190 J=IP1,K HFT10500 190 SM=SM+A(I,J)*DBLE(B(J,JB)) HFT10600 200 SM1=SM HFT10700 210 B(I,JB)=(B(I,JB)-SM1)/A(I,I) HFT10800 C HFT10900 C COMPLETE COMPUTATION OF SOLUTION VECTOR. HFT11000 C .. HFT11100 IF (K.EQ.N) GO TO 240 HFT11200 DO 220 J=KP1,N HFT11300 220 B(J,JB)=SZERO HFT11400 DO 230 I=1,K HFT11500 230 CALL H12 (2,I,KP1,N,A(I,1),MDA,G(I),B(1,JB),1,MDB,1) HFT11600 C HFT11700 C RE-ORDER THE SOLUTION VECTOR TO COMPENSATE FOR THE HFT11800 C COLUMN INTERCHANGES. HFT11900 C .. HFT12000 240 DO 250 JJ=1,LDIAG HFT12100 J=LDIAG+1-JJ HFT12200 IF (IP(J).EQ.J) GO TO 250 HFT12300 L=IP(J) HFT12400 TMP=B(L,JB) HFT12500 B(L,JB)=B(J,JB) HFT12600 B(J,JB)=TMP HFT12700 250 CONTINUE HFT12800 260 CONTINUE HFT12900 C .. HFT13000 C THE SOLUTION VECTORS, X, ARE NOW HFT13100 C IN THE FIRST N ROWS OF THE ARRAY B(,). HFT13200 C HFT13300 270 KRANK=K HFT13400 RETURN HFT13500 END HFT13600