      SUBROUTINE HFTI2(A,MDA,M,N,B,MDB,NB,D,TAU,K,RNORM,H,G,IP,IERR)
C     ------------------------------------------------------------------
C     DIMENSION A(MDA,N),(B(MDB,NB) OR B(M))
C     DIMENSION D(L) WHERE L = MIN(M,N)
C     DIMENSION RNORM(NB),H(N),G(N),IP(N)
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     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 ARRAY D WILL CONTAIN THE DIAGONAL ELEMENTS R(1,1),...,R(L,L).
C     THE SUBROUTINE SETS K TO BE THE NUMBER OF DIAGONAL ELEMENTS THAT
C     EXCEED TAU IN MAGNITUDE.  THEN THE SOLUTION OF MINIMUM EUCLIDEAN
C     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 USER MUST DIMENSION ALL ARRAYS APPEARING IN THE CALL LIST..
C     A(MDA,N),(B(MDB,NB) OR B(M)),D(L),RNORM(NB),H(N),G(N),IP(N). THIS
C     PERMITS THE SOLUTION OF A RANGE OF PROBLEMS IN THE SAME ARRAY
C     SPACE.
C
C     THE ENTIRE SET OF PARAMETERS FOR HFTI2 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     D(*)              THE ARRAY OF DIAGONAL ELEMENTS OF THE
C                       TRIANGULAR MATRIX R
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
C     IERR              ERROR INDICATOR. IF NO INPUT ERRORS ARE
C                       DETECTED THEN IERR IS SET TO 0. OTHERWISE
C                          IERR = 1  IF MDA.LT.M
C                          IERR = 2  IF NB.GT.1 AND MDB.LT.MAX(M,N)
C                       THESE ERRORS ARE FATAL.
C     ------------------------------------------------------------------
      DIMENSION A(MDA,N),B(MDB,NB),D(1),H(N),G(N),RNORM(NB)
      INTEGER IP(N)
      DOUBLE PRECISION SM,DZERO
      SZERO=0.
      DZERO=0.D0
      FACTOR=0.001
C
      K=0
      LDIAG=MIN0(M,N)
      IF (LDIAG.LE.0) GO TO 270
      IF (MDA.GE.M) GO TO 5
      IERR=1
      RETURN
    5 CONTINUE
C
      IF (.NOT.(NB.GT.1.AND.MAX0(M,N).GT.MDB)) GO TO 6
      IERR=2
      RETURN
    6 CONTINUE
C
          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)=SZERO
                  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     CONTINUE
          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     ALSO STORE THE DIAGONAL ELEMENTS IN THE ARRAY D.
C    ..
          DO 90 J=1,LDIAG
          IF (ABS(A(J,J)).LE.TAU) GO TO 100
   90     D(J)=A(J,J)
      K=LDIAG
      KP1=K+1
      GO TO 110
  100 K=J-1
      KP1=J
          DO 105 J=KP1,LDIAG
  105     D(J)=A(J,J)
C
C     COMPUTE THE NORMS OF THE RESIDUAL VECTORS.
C
  110 IF (NB.LE.0) GO TO 140
          DO 130 JB=1,NB
          TMP=SZERO
          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 150 JB=1,NB
              DO 150 I=1,N
  150         B(I,JB)=SZERO
      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=DZERO
              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=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)=SZERO
              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 IERR=0
      RETURN
      END
