# dhfti.f

          SUBROUTINE DHFTI (A, MDA, M, N, B, MDB, NB, TAU, KRANK, RNORM, H,
         +   G, IP)
    C***BEGIN PROLOGUE  DHFTI
    C***PURPOSE  Solve a least squares problem for banded matrices using
    C            sequential accumulation of rows of the data matrix.
    C            Exactly one right-hand side vector is permitted.
    C***LIBRARY   SLATEC
    C***CATEGORY  D9
    C***TYPE      DOUBLE PRECISION (HFTI-S, DHFTI-D)
    C***KEYWORDS  CURVE FITTING, LEAST SQUARES
    C***AUTHOR  Lawson, C. L., (JPL)
    C           Hanson, R. J., (SNLA)
    C***DESCRIPTION
    C
    C     DIMENSION A(MDA,N),(B(MDB,NB) or B(M)),RNORM(NB),H(N),G(N),IP(N)
    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 will compute an integer, KRANK, equal to the number
    C     of diagonal terms of R that exceed TAU in magnitude. Then a
    C     solution of minimum Euclidean length is computed using the first
    C     KRANK 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)),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 DHFTI are
    C
    C     INPUT.. All TYPE REAL variables are DOUBLE PRECISION
    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 DHFTI.
    C
    C     OUTPUT.. All TYPE REAL variables are DOUBLE PRECISION
    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     KRANK             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***REFERENCES  C. L. Lawson and R. J. Hanson, Solving Least Squares
    C                 Problems, Prentice-Hall, Inc., 1974, Chapter 14.
    C***ROUTINES CALLED  D1MACH, DH12, XERMSG
    C***REVISION HISTORY  (YYMMDD)
    C   790101  DATE WRITTEN
    C   890531  Changed all specific intrinsics to generic.  (WRB)
    C   891006  Cosmetic changes to prologue.  (WRB)
    C   891006  REVISION DATE from Version 3.2
    C   891214  Prologue converted to Version 4.0 format.  (BAB)
    C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    C   901005  Replace usage of DDIFF with usage of D1MACH.  (RWC)
    C   920501  Reformatted the REFERENCES section.  (WRB)
    C***END PROLOGUE  DHFTI