      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
