*DECK HFTI SUBROUTINE HFTI (A, MDA, M, N, B, MDB, NB, TAU, KRANK, RNORM, H, + G, IP) C***BEGIN PROLOGUE HFTI C***PURPOSE Solve a linear least squares problems by performing a QR C factorization of the matrix using Householder C transformations. C***LIBRARY SLATEC C***CATEGORY D9 C***TYPE SINGLE PRECISION (HFTI-S, DHFTI-D) C***KEYWORDS CURVE FITTING, LINEAR LEAST SQUARES, QR FACTORIZATION 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 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 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 H12, R1MACH, 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 DIFF with usage of R1MACH. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE HFTI DIMENSION A(MDA,*),B(MDB,*),H(*),G(*),RNORM(*) INTEGER IP(*) DOUBLE PRECISION SM,DZERO SAVE RELEPS DATA RELEPS /0.E0/ C***FIRST EXECUTABLE STATEMENT HFTI IF (RELEPS.EQ.0) RELEPS = R1MACH(4) SZERO=0. DZERO=0.D0 FACTOR=0.001 C K=0 LDIAG=MIN(M,N) IF (LDIAG.LE.0) GO TO 270 IF (.NOT.MDA.LT.M) GO TO 5 NERR=1 IOPT=2 CALL XERMSG ('SLATEC', 'HFTI', 'MDA.LT.M, PROBABLE ERROR.', + NERR, IOPT) RETURN 5 CONTINUE C IF (.NOT.(NB.GT.1.AND.MAX(M,N).GT.MDB)) GO TO 6 NERR=2 IOPT=2 CALL XERMSG ('SLATEC', 'HFTI', + 'MDB.LT.MAX(M,N).AND.NB.GT.1. PROBABLE ERROR.', NERR, IOPT) 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 IF (FACTOR*H(LMAX) .GT. HMAX*RELEPS) GO TO 50 C C COMPUTE SQUARED COLUMN LENGTHS AND FIND LMAX C .. 20 LMAX=J DO 40 L=J,N H(L)=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 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 .. 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=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+A(I,J)*DBLE(B(J,JB)) 200 SM1=SM 210 B(I,JB)=(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 KRANK=K RETURN END