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
*++++++++++++++++++++++++ End of file hfti.f +++++++++++++++++++++++++++