C DQRLS SOURCE PV090527 23/02/07 11:58:37 11592
subroutine dqrls ( a, lda, m, n, tol, kr, b, x,r,jpvt,qraux,work,
# itask, ind )
c*****************************************************************************80
c
c! DQRLS factors and solves a linear system in the least squares sense.
c
c Discussion:
c
c The linear system may be overdetermined, underdetermined or singular.
c The solution is obtained using a QR factorization of the
c coefficient matrix.
c
c DQRLS can be efficiently used to solve several least squares
c problems with the same matrix A. The first system is solved
c with ITASK = 1. The subsequent systems are solved with
c ITASK = 2, to avoid the recomputation of the matrix factors.
c The parameters KR, JPVT, and QRAUX must not be modified
c between calls to DQRLS.
c
c DQRLS is used to solve in a least squares sense
c overdetermined, underdetermined and singular linear systems.
c The system is A*X approximates B where A is M by N.
c B is a given M-vector, and X is the N-vector to be computed.
c A solution X is found which minimimzes the sum of squares (2-norm)
c of the residual, A*X - B.
c
c The numerical rank of A is determined using the tolerance TOL.
c
c DQRLS uses the LINPACK subroutine DQRDC to compute the QR
c factorization, with column pivoting, of an M by N matrix A.
c
c Modified:
c
c 15 April 2012
c
c Reference:
c
c David Kahaner, Cleve Moler, Steven Nash,
c Numerical Methods and Software,
c Prentice Hall, 1989,
c ISBN: 0-13-627258-4,
c LC: TA345.K34.
c
c Parameters:
c
c Input/output, real*8 A(LDA,N), an M by N matrix.
c On input, the matrix whose decomposition is to be computed.
c In a least squares data fitting problem, A(I,J) is the
c value of the J-th basis (model) function at the I-th data point.
c On output, A contains the output from DQRDC. The triangular matrix R
c of the QR factorization is contained in the upper triangle and
c information needed to recover the orthogonal matrix Q is stored
c below the diagonal in A and in the vector QRAUX.
c
c Input, integer LDA, the leading dimension of A.
c
c Input, integer M, the number of rows of A.
c
c Input, integer N, the number of columns of A.
c
c Input, real*8 TOL, a relative tolerance used to determine the
c numerical rank. The problem should be scaled so that all the elements
c of A have roughly the same absolute accuracy EPS. Then a reasonable
c value for TOL is roughly EPS divided by the magnitude of the largest
c element.
c
c Output, integer KR, the numerical rank.
c
c Input, real*8 B(M), the right hand side of the linear system.
c
c Output, real*8 X(N), a least squares solution to the linear
c system.
c
c Output, real*8 R(M), the residual, B - A*X. R may
c overwrite B.
c
c Workspace, integer JPVT(N), required if ITASK = 1.
c Columns JPVT(1), ..., JPVT(KR) of the original matrix are linearly
c independent to within the tolerance TOL and the remaining columns
c are linearly dependent. ABS ( A(1,1) ) / ABS ( A(KR,KR) ) is an estimate
c of the condition number of the matrix of independent columns,
c and of R. This estimate will be .le. 1/TOL.
c
c Workspace, real*8 QRAUX(N), required if ITASK = 1.
c
c Workspace, real*8 WORK(N), required if ITASK = 1.
c
c Input, integer ITASK.
c 1, DQRLS factors the matrix A and solves the least squares problem.
c 2, DQRLS assumes that the matrix A was factored with an earlier
c call to DQRLS, and only solves the least squares problem.
c
c Output, integer IND, error code.
c 0: no error
c -1: LDA .lt. M (fatal error)
c -2: N .lt. 1 (fatal error)
c -3: ITASK .lt. 1 (fatal error)
c
implicit real*8 (a-h,o-z)
implicit integer (i-n)
integer lda
integer m
integer n
real*8 a(lda,n)
real*8 b(m)
integer ind
integer itask
integer jpvt(n)
integer kr
real*8 qraux(n)
real*8 r(m)
real*8 tol
real*8 work(n)
real*8 x(n)
if ( lda .lt. m ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'DQRLS - Fatal error!'
write ( *, '(a)' ) ' LDA .lt. M.'
stop
end if
if ( n .le. 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'DQRLS - Fatal error!'
write ( *, '(a)' ) ' N .le. 0.'
stop
end if
if ( itask .lt. 1 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'DQRLS - Fatal error!'
write ( *, '(a)' ) ' ITASK .lt. 1.'
stop
end if
ind = 0
c Factor the matrix.
if ( itask .eq. 1 ) then
call dqrank ( a, lda, m, n, tol, kr, jpvt, qraux, work )
end if
c Solve the least-squares problem.
call dqrlss ( a, lda, m, n, kr, b, x, r, jpvt, qraux )
return
end