;+ ; NAME: ; QRSOLV ; ; AUTHOR: ; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770 ; craigm@lheamail.gsfc.nasa.gov ; ; PURPOSE: ; Solve a linear equation after performing QR factorization ; ; MAJOR TOPICS: ; Linear Systems ; ; CALLING SEQUENCE: ; X = QRSOLV(A, R, B, PIVOTS=IPVT) ; ; DESCRIPTION: ; ; The procedure QRSOLV completes the solution of a linear equation, ; ; A ## x = B ; ; after the MxN matrix has been factorized by QR decomposition. ; After being factorized once using QRFAC, the matrices can be used ; for multiple righthand sides (i.e., different B's). ; ; The solution technique is to first compute the factorization using ; QRFAC, which yields the orthogonal matrix Q and the upper ; triangular matrix R. [ Actually, Q is represented by its ; Householder reflectors. ] Then the solution vector, X, is computed ; using QRSOLV. ; ; If pivoting was performed in the factorization, the permutation ; vector IPVT returned by QRFAC must also be passed to QRSOLV. ; ; ; PARAMETERS: ; ; A - upon input, the factorized matrix A, returned by QRFAC. ; ; R - upon input, the upper diagonal matrix R, returned by QRFAC. ; ; B - upon input, the righthand vector B, which fits into the ; equation, A ## x = B ; ; X - upon ouptut, the solution vector X, to the above linear ; equation. For an overdetermined system, X is the least ; squares solution which minimizes TOTAL( (A ## X - B)^2 ). ; ; ; KEYWORD PARAMETERS: ; ; PIVOTS - upon input, the permutation matrix IPVT returned by ; QRFAC, if pivoting is to be performed. ; ; ; EXAMPLE: ; ; Solve the equation A ## X = B, in the least squares sense, where: ; ; A = [[1.0,1.0,1.0,1.0,1.0,1.0],$ ; [0.6,0.8,0.5,0.8,0.7,0.9],$ ; [0.2,0.3,0.1,0.4,0.3,0.4]] ; ; and B = [0.57E,0.69,0.5,0.7,0.6,0.8] ; ; qrfac, a, r, ipvt, /PIVOT ; x = qrsolv(a, r, b, PIVOTS=ipvt) ; ; print, x ; 0.0834092 0.852273 -0.179545 ; ; REFERENCES: ; ; More', Jorge J., "The Levenberg-Marquardt Algorithm: ; Implementation and Theory," in *Numerical Analysis*, ed. Watson, ; G. A., Lecture Notes in Mathematics 630, Springer-Verlag, 1977. ; ; MODIFICATION HISTORY: ; Written (taken from MPFIT), CM, Feb 2002 ; Usage message, error checking, CM, 15 Mar 2002 ; Error checking is fixed, CM, 10 May 2002 ; Found error in return of permuted results, CM, 21 May 2004 ; ; $Id: qrsolv.pro,v 1.4 2004/05/22 02:16:02 craigm Exp $ ; ;- ; Copyright (C) 2002, 2004, Craig Markwardt ; This software is provided as is without any warranty whatsoever. ; Permission to use, copy, modify, and distribute modified or ; unmodified copies is granted, provided this copyright and disclaimer ; are included unchanged. ;- function qrsolv, q, r, b, pivots=ipvt0 if n_params() EQ 0 then begin USAGE: message, 'USAGE:', /info message, 'X = QRSOLV(Q, R, B, [PIVOTS=IPVT])', /info message, ' Q and R are factorization from QRFAC', /info message, ' B is right hand side of equation Q # X = B', /info return, 0 endif sz = size(q) m = sz(1) n = sz(2) if sz(0) NE 2 OR m LT n then $ goto, USAGE szr = size(r) if szr(0) NE 2 OR szr(1) NE szr(2) OR szr(1) NE n then $ goto, USAGE if n_elements(b) NE m then $ goto, USAGE delm = lindgen(n) * (n+1) ;; Diagonal elements of r ;; Default pivoting if not specified if n_elements(ipvt0) EQ 0 then ipvt = lindgen(n) $ else ipvt = ipvt0 ;; Multiply QT * B to get QTB. Because Q is stored as reflectors, ;; they must be expanded explicitly. wa4 = b for j=0L, n-1 do begin lj = ipvt(j) temp3 = q(j,lj) if temp3 NE 0 then begin fj = q(j:*,lj) wj = wa4(j:*) ;; *** optimization wa4(j:*) wa4(j) = wj - fj * total(fj*wj) / temp3 endif endfor qtb = wa4(0:n-1) ;; Get some initial parameters diag = r(delm) x = qtb ;; Solve the triangular system for x. If the system is singular ;; then obtain a least squares solution nsing = n wh = where(diag EQ 0, ct) if ct GT 0 then begin nsing = wh(0) x(nsing:*) = 0 endif if nsing GE 1 then begin ;; Could use LUSOL here, but it turns out that the largest ;; expense for most reasonable problems is in computing QT*B, ;; rather than the solution here. For huge problems where that ;; is not true, you are going to run into other problems. x(nsing-1) = x(nsing-1)/diag(nsing-1) ;; Degenerate case ;; *** Reverse loop *** for j=nsing-2,0,-1 do begin sum = total(r(j+1:nsing-1,j)*x(j+1:nsing-1)) x(j) = (x(j)-sum)/diag(j) endfor endif ;; Permute the components of x back to original xp = x & xp(ipvt) = x return, xp end