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