SUBROUTINE HFTI(A,MDA,M,N,B,MDB,NB,TAU,KRANK,RNORM,H,G,IP) C RANK-DEFICIENT LEAST SQUARES ALGORITHM AS DESCRIBED IN: C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 C TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974 C A(*,*),MDA,M,N THE ARRAY A INITIALLY CONTAINS THE M x N MATRIX A C OF THE LEAST SQUARES PROBLEM AX = B. C THE FIRST DIMENSIONING PARAMETER MDA MUST SATISFY C MDA >= M. EITHER M >= N OR M < N IS PERMITTED. C THERE IS NO RESTRICTION ON THE RANK OF A. C THE MATRIX A WILL BE MODIFIED BY THE SUBROUTINE. C B(*,*),MDB,NB IF NB = 0 THE SUBROUTINE WILL MAKE NO REFERENCE C TO THE ARRAY B. IF NB > 0 THE ARRAY B() MUST C INITIALLY CONTAIN THE M x NB MATRIX B OF THE C THE LEAST SQUARES PROBLEM AX = B AND ON RETURN C THE ARRAY B() WILL CONTAIN THE N x NB SOLUTION X. C IF NB>1 THE ARRAY B() MUST BE DOUBLE SUBSCRIPTED C WITH FIRST DIMENSIONING PARAMETER MDB>=MAX(M,N), C IF NB=1 THE ARRAY B() MAY BE EITHER SINGLE OR C DOUBLE SUBSCRIPTED. C TAU ABSOLUTE TOLERANCE PARAMETER FOR PSEUDORANK C DETERMINATION, PROVIDED BY THE USER. C KRANK PSEUDORANK OF A, SET BY THE SUBROUTINE. C RNORM ON EXIT, RNORM(J) WILL CONTAIN THE EUCLIDIAN C NORM OF THE RESIDUAL VECTOR FOR THE PROBLEM C DEFINED BY THE J-TH COLUMN VECTOR OF THE ARRAY B. C H(), G() ARRAYS OF WORKING SPACE OF LENGTH >= N. C IP() INTEGER ARRAY OF WORKING SPACE OF LENGTH >= N C RECORDING PERMUTATION INDICES OF COLUMN VECTORS INTEGER I,J,JB,K,KP1,KRANK,L,LDIAG,LMAX,M, . MDA,MDB,N,NB,IP(N) DOUBLE PRECISION A(MDA,N),B(MDB,NB),H(N),G(N),RNORM(NB),FACTOR, . TAU,ZERO,HMAX,DIFF,TMP,DDOT,DNRM2,U,V DIFF(U,V)= U-V DATA ZERO/0.0D0/, FACTOR/1.0D-3/ K=0 LDIAG=MIN(M,N) IF(LDIAG.LE.0) GOTO 270 C COMPUTE LMAX DO 80 J=1,LDIAG IF(J.EQ.1) GOTO 20 LMAX=J DO 10 L=J,N H(L)=H(L)-A(J-1,L)**2 10 IF(H(L).GT.H(LMAX)) LMAX=L IF(DIFF(HMAX+FACTOR*H(LMAX),HMAX).GT.ZERO) . GOTO 50 20 LMAX=J DO 40 L=J,N H(L)=ZERO DO 30 I=J,M 30 H(L)=H(L)+A(I,L)**2 40 IF(H(L).GT.H(LMAX)) LMAX=L HMAX=H(LMAX) C COLUMN INTERCHANGES IF NEEDED 50 IP(J)=LMAX IF(IP(J).EQ.J) GOTO 70 DO 60 I=1,M TMP=A(I,J) A(I,J)=A(I,LMAX) 60 A(I,LMAX)=TMP H(LMAX)=H(J) C J-TH TRANSFORMATION AND APPLICATION TO A AND B 70 I=MIN(J+1,N) CALL H12(1,J,J+1,M,A(1,J),1,H(J),A(1,I),1,MDA,N-J) 80 CALL H12(2,J,J+1,M,A(1,J),1,H(J),B,1,MDB,NB) C DETERMINE PSEUDORANK DO 90 J=1,LDIAG 90 IF(ABS(A(J,J)).LE.TAU) GOTO 100 K=LDIAG GOTO 110 100 K=J-1 110 KP1=K+1 C NORM OF RESIDUALS DO 130 JB=1,NB 130 RNORM(JB)=DNRM2(M-K,B(KP1,JB),1) IF(K.GT.0) GOTO 160 DO 150 JB=1,NB DO 150 I=1,N 150 B(I,JB)=ZERO GOTO 270 160 IF(K.EQ.N) GOTO 180 C HOUSEHOLDER DECOMPOSITION OF FIRST K ROWS DO 170 I=K,1,-1 170 CALL H12(1,I,KP1,N,A(I,1),MDA,G(I),A,MDA,1,I-1) 180 DO 250 JB=1,NB C SOLVE K*K TRIANGULAR SYSTEM DO 210 I=K,1,-1 J=MIN(I+1,N) 210 B(I,JB)=(B(I,JB)-DDOT(K-I,A(I,J),MDA,B(J,JB),1))/A(I,I) C COMPLETE SOLUTION VECTOR IF(K.EQ.N) GOTO 240 DO 220 J=KP1,N 220 B(J,JB)=ZERO DO 230 I=1,K 230 CALL H12(2,I,KP1,N,A(I,1),MDA,G(I),B(1,JB),1,MDB,1) C REORDER SOLUTION ACCORDING TO PREVIOUS COLUMN INTERCHANGES 240 DO 250 J=LDIAG,1,-1 IF(IP(J).EQ.J) GOTO 250 L=IP(J) TMP=B(L,JB) B(L,JB)=B(J,JB) B(J,JB)=TMP 250 CONTINUE 270 KRANK=K END