123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288 |
- *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
|