      SUBROUTINE HFTI (A,MDA,M,N,B,MDB,NB,TAU,K,RNORM,H,G,IP)
C-----------------------------------------------------------------------
C     DIMENSION A(MDA,N),(B(MDB,NB) OR B(M)),RNORM(NB),H(N),G(N),IP(N)
C
C     Written by C.L. Lawson and R.J. Hanson.
C     From the book Solving Least Squares Problems, Prentice-Hall, 1974.
C     For algorithmic details see algorithm HFTI in chapter 14.
C
C     ABSTRACT
C
C     This subroutine solves a linear least squares problem or a set of
C     linear least squares problems having the same matrix but different
C     Right-side vectors.  The problem data consists of an M by N matrix
C     A, an M by NB matrix B, and an absolute tolerance parameter TAU
C     whose usage is described below.  The NB column vectors of B
C     represent right-side vectors for NB distinct linear least squares
C     problems.
C
C     This set of problems can also be written as the matrix least
C     squares problem
C
C                       AX = B,
C
C     Where X is the N by NB solution matrix.
C
C     Note that if B is the M by M identity matrix, then X will be the
C     pseudo-inverse of A.
C
C     This subroutine first transforms the augmented matrix (A B) to a
C     matrix (R C) using premultiplying Householder transformations with
C     column interchanges.  All subdiagonal elements in the matrix R are
C     zero and its diagonal elements satisfy
C
C                       ABS(R(I,I)).GE.ABS(R(I+1,I+1)),
C
C                       I = 1,...,L-1, where
C
C                       L = MIN(M,N).
C
C     The subroutine sets K to be the number of diagonal elements
C     of R that exceed TAU in magnitude. Then the solution of minimum
C     Euclidean length is computed using the first K rows of (R C).
C
C     To be specific we suggest that the user consider an easily
C     computable matrix norm, such as, the maximum of all column sums of
C     magnitudes.
C
C     Now if the relative uncertainty of B is EPS, (norm of uncertainty/
C     norm of B), it is suggested that TAU be set approximately equal to
C     EPS*(norm of A).
C
C     The entire set of parameters for HFTI are
C
C     INPUT..
C
C     A(*,*),MDA,M,N    The array A(*,*) initially contains the M by N
C                       matrix A of the least squares problem AX = B.
C                       The first dimensioning parameter of the array
C                       A(*,*) is MDA, which must satisfy MDA.GE.M
C                       either M.GE.N or M.LT.N is permitted.  There
C                       is no restriction on the rank of A.  The
C                       condition MDA.LT.M is considered an error.
C
C     B(*),MDB,NB       If NB = 0 the subroutine will perform the
C                       orthogonal decomposition but will make no
C                       references to the array B(*).  If NB.GT.0
C                       the array B(*) must initially contain the M by
C                       NB matrix B of the least squares problem AX =
C                       B.  If NB.GE.2 the array B(*) must be doubly
C                       subscripted with first dimensioning parameter
C                       MDB.GE.MAX(M,N).  If NB = 1 the array B(*) may
C                       be either doubly or singly subscripted.  In
C                       the latter case the value of MDB is arbitrary
C                       but it should be set to some valid integer
C                       value such as MDB = M.
C
C                       The condition of NB.GT.1.AND.MDB.LT. MAX(M,N)
C                       is considered an error.
C
C     TAU               Absolute tolerance parameter provided by user
C                       for pseudorank determination.
C
C     H(*),G(*),IP(*)   Arrays of working space used by HFTI.
C
C     OUTPUT..
C
C     A(*,*)            The contents of the array A(*,*) will be
C                       modified by the subroutine.  These contents
C                       are not generally required by the user.
C
C     B(*)              On return the array B(*) will contain the N by
C                       NB solution matrix X.
C
C     K                 Set by the subroutine to indicate the
C                       pseudorank of A.
C
C     RNORM(*)          On return, RNORM(J) will contain the Euclidean
C                       norm of the residual vector for the problem
C                       defined by the J-th column vector of the array
C                       B(*,*) for J = 1,...,NB.
C
C     H(*),G(*)         On return these arrays respectively contain
C                       elements of the pre- and post-multiplying
C                       Householder transformations used to compute
C                       the minimum Euclidean length solution.
C
C     IP(*)             Array in which the subroutine records indices
C                       describing the permutation of column vectors.
C                       The contents of arrays H(*),G(*) and IP(*)
C                       are not generally required by the user.
C-----------------------------------------------------------------------
      DIMENSION A(MDA,N),B(MDB,*),RNORM(*),H(N),G(N)
      INTEGER IP(N)
      DOUBLE PRECISION SM
C---------------------
      DATA FACTOR /1.E-3/
C
      K = 0
      LDIAG = MIN0(M,N)
      IF (LDIAG .LE. 0) GO TO 270
          DO 80 J = 1,LDIAG
          IF (J .EQ. 1) GO TO 20
C
C     Update squared column lengths and find LMAX
C    ..
          LMAX = J
              DO 10 L = J,N
              H(L) = H(L) - A(J-1,L)**2
              IF (H(L) .GT. H(LMAX)) LMAX = L
   10         CONTINUE
          Z = HMAX + FACTOR*H(LMAX)
          IF (Z .GT. HMAX) GO TO 50
C
C     Compute squared column lengths and find LMAX
C    ..
   20     LMAX=J
              DO 40 L = J,N
              H(L) = 0.0
                  DO 30 I = J,M
   30             H(L) = H(L) + A(I,L)**2
              IF (H(L) .GT. H(LMAX)) LMAX = L
   40         CONTINUE
          HMAX = H(LMAX)
C    ..
C     LMAX has been determined
C
C     Do column interchanges if needed.
C    ..
   50     IP(J) = LMAX
          IF (IP(J) .EQ. J) GO TO 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
C     Compute the J-th transformation and apply it to A and B.
C    ..
   70     CALL H12 (1,J,J+1,M,A(1,J),1,H(J),A(1,J+1),1,MDA,N-J)
   80     CALL H12 (2,J,J+1,M,A(1,J),1,H(J),B,1,MDB,NB)
C
C     Determine the pseudorank, K, using the tolerance, TAU.
C    ..
          DO 90 J = 1,LDIAG
          IF (ABS(A(J,J)) .LE. TAU) GO TO 100
   90     CONTINUE
      K = LDIAG
      GO TO 110
  100 K = J - 1
  110 KP1 = K + 1
C
C     Compute the norms of the residual vectors.
C
      IF (NB .LE. 0) GO TO 140
          DO 130 JB = 1,NB
          TMP = 0.0
          IF (KP1 .GT. M) GO TO 130
              DO 120 I = KP1,M
  120         TMP = TMP + B(I,JB)**2
  130     RNORM(JB) = SQRT(TMP)
  140 CONTINUE
C                                           Special for pseudorank = 0
      IF (K .GT. 0) GO TO 160
      IF (NB .LE. 0) GO TO 270
          DO 151 JB = 1,NB
              DO 150 I = 1,N
  150         B(I,JB) = 0.0
  151     CONTINUE
      GO TO 270
C
C     If the pseudorank is less than N compute Householder
C     decomposition of first K rows.
C    ..
  160 IF (K .EQ. N) GO TO 180
          DO 170 II = 1,K
          I = KP1 - II
  170     CALL H12 (1,I,KP1,N,A(I,1),MDA,G(I),A,MDA,1,I-1)
  180 CONTINUE
C
C
      IF (NB .LE. 0) GO TO 270
          DO 260 JB = 1,NB
C
C     Solve the K by K triangular system.
C    ..
              DO 210 L = 1,K
              SM = 0.D0
              I = KP1 - L
              IF (I .EQ. K) GO TO 200
              IP1 = I + 1
                  DO 190 J = IP1,K
  190             SM = SM + DBLE(A(I,J))*DBLE(B(J,JB))
  200         SM1 = SNGL(DBLE(B(I,JB)) - SM)
  210         B(I,JB) = SM1/A(I,I)
C
C     Complete computation of solution vector.
C    ..
          IF (K .EQ. N) GO TO 240
              DO 220 J = KP1,N
  220         B(J,JB) = 0.0
              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
C      Re-order the solution vector to compensate for the
C      column interchanges.
C    ..
  240         DO 250 JJ = 1,LDIAG
              J = LDIAG + 1 - JJ
              IF (IP(J) .EQ. J) GO TO 250
              L = IP(J)
              TMP = B(L,JB)
              B(L,JB) = B(J,JB)
              B(J,JB) = TMP
  250         CONTINUE
  260     CONTINUE
C    ..
C     The solution vectors, X, are now
C     in the first  N  rows of the array B(,).
C
  270 RETURN
      END ! of hfti


      SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)
C-----------------------------------------------------------------------
C     Written by C.L. Lawson and R.J. Hanson.  Modified by A.H. Morris.
C     From the book Solving Least Squares Problems, Prentice-Hall, 1974.
C
C     Construction and/or application of a single
C     Householder transformation..     Q = I + U*(U**T)/B
C
C     MODE    = 1 or 2   to select algorithm  H1  or  H2 .
C     LPIVOT is the index of the pivot element.
C     L1,M   If L1 .LE. M   the transformation will be constructed to
C            zero elements indexed from L1 through M.   If L1 GT. M
C            the subroutine does an identity transformation.
C     U(),IUE,UP    On entry to H1 U() contains the pivot vector.
C                   IUE is the storage increment between elements.
C                   On exit from H1 U() and UP contain quantities
C                   defining the vector U of the Householder 
C                   transformation.   On entry to H2 U() and UP should 
C                   contain quantities previously computed by H1.
C                   These will not be modified by H2.
C     C()    On entry to H1 or H2 C() contains a matrix which will be
C            regarded as a set of vectors to which the Householder
C            transformation is to be applied.  On exit C() contains the
C            set of transformed vectors.
C     ICE    Storage increment between elements of vectors in C().
C     ICV    Storage increment between vectors in C().
C     NCV    Number of vectors in C() to be transformed. If NCV .LE. 0
C            no operations will be done on C().
C-----------------------------------------------------------------------
      DIMENSION U(IUE,M), C(*)
      DOUBLE PRECISION SM,B
C
      IF (0 .GE. LPIVOT .OR. LPIVOT .GE. L1 .OR. L1 .GT. M) RETURN
      CL = ABS(U(1,LPIVOT))
      IF (MODE .EQ. 2) GO TO 60
C
C                            ****** Construct the transformation. ******
C
          DO 10 J = L1,M
   10     CL = AMAX1(ABS(U(1,J)),CL)
      IF (CL .LE. 0.0) GO TO 130
      D = U(1,LPIVOT)/CL
      SM = D*D
          DO 20 J = L1,M
          D = U(1,J)/CL
   20     SM = SM + DBLE(D*D)
C
      SM1 = SNGL(SM)
      CL = CL*SQRT(SM1)
      IF (U(1,LPIVOT) .GT. 0.0) CL = -CL
      UP = U(1,LPIVOT) - CL
      U(1,LPIVOT) = CL
      GO TO 70
C
C            ****** Apply the transformation  I+U*(U**T)/B  to C. ******
C
   60 IF (CL) 130,130,70
   70 IF (NCV .LE. 0) RETURN
      B = DBLE(UP)*DBLE(U(1,LPIVOT))
C
C                       B  must be nonpositive here.  If B = 0., return.
C
      IF (B .GE. 0.D0) GO TO 130
      B = 1.D0/B
      I2 = 1 - ICV + ICE*(LPIVOT - 1)
      INCR = ICE*(L1 - LPIVOT)
          DO 120 J = 1,NCV
          I2 = I2 + ICV
          I3 = I2 + INCR
          I4 = I3
          SM = DBLE(C(I2))*DBLE(UP)
              DO 90 I = L1,M
              SM = SM + DBLE(C(I3))*DBLE(U(1,I))
   90         I3 = I3 + ICE
          IF (SM .EQ. 0.D0) GO TO 120
          SM = SM*B
          C(I2) = C(I2) + SNGL(SM*DBLE(UP))
              DO 110 I = L1,M
              C(I4) = C(I4) + SNGL(SM*DBLE(U(1,I)))
  110         I4 = I4 + ICE
  120     CONTINUE
  130 RETURN
      END ! of H12
