File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices / standard_floating_least_squares.ads (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:24 2000 UTC (23 years, 8 months ago) by maekawa
Branch: PHC, MAIN
CVS Tags: v2, maekawa-ipv6, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, HEAD Changes since 1.1: +0 -0
lines
Import the second public release of PHCpack.
OKed by Jan Verschelde.
|
with Standard_Floating_Vectors;
with Standard_Floating_Matrices;
package Standard_Floating_Least_Squares is
-- DESCRIPTION :
-- This package provides an implementation of Least Squares algorithm
-- for matrices of standard floating-point numbers.
procedure QRLS ( x : in out Standard_Floating_Matrices.Matrix;
ldx,n,k : in natural;
qraux,y : in Standard_Floating_Vectors.Vector;
qy,qty,b,rsd,xb : out Standard_Floating_Vectors.Vector;
job : in natural; info : out natural );
-- DESCRIPTION :
-- Applies the output of Standard_Floating_QR_Decomposition.QRD to
-- compute coordinate transformations, projections and least quares
-- solutions. For k <= min(n,p), let xk be the matrix
--
-- xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k)))
--
-- formed from columnns jpvt(1), ... ,jpvt(k) of the original
-- n x p matrix x that was input to dqrdc (if no pivoting was done,
-- xk consists of the first k columns of x in their original order).
-- QRD produces a factored orthogonal matrix q and an upper triangular
-- matrix r such that
-- xk = q * (r)
-- (0)
-- this information is contained in coded form in x and qraux.
-- ON ENTRY :
-- x contains the output of QRD, of size ldx times p;
-- ldx leading dimension of x;
-- n number of rows in the matrix xk, must be same as in QRD
-- k number of columns of the matrix xk, k <= min(n,p)
-- qraux contains p entries, auxiliary output from QRD;
-- y n-vector to be manipulated by QRLS;
-- job specifies what is to be computed, job has the decimal
-- expansion abcde, with the following meaning :
-- if a /= 0, compute qy,
-- if b,c,d, or e /= 0, compute qty,
-- if c /= 0, compute b,
-- if d /= 0, compute rsd,
-- if e /= 0, compute xb.
-- Note that a request to compute b, rsd, or xb
-- automatically triggers the computation of qty, for
-- which an array must be provided in the calling sequence.
-- ON RETURN :
-- x may be altered, used as work space;
-- qy qy'length = n,
-- qy conntains q*y, if its computation has been requested;
-- qty qty'length = n,
-- qty contains trans(q)*y, if its computation has been
-- requested; here trans(q) is the transpose of the matrix q;
-- b b'length = k,
-- b contains the solution of the least squares problem
--
-- minimize norm2(y - xk*b),
--
-- if its computation has been requested. (Note that
-- if pivoting was requested in dqrdc, the j-th
-- component of b will be associated with column jpvt(j)
-- of the original matrix x that was input into dqrdc.)
-- rsd rsd'length = n,
-- rsd contains the least squares residual y - xk*b,
-- if its computation has been requested; rsd is also the
-- orthogonal projection of y onto the orthogonal complement
-- of the column space of xk;
-- xb x'length = n, xb contains the least squares approximation
-- xk*b, if its computation has been requested; xb is also
-- the orthogonal projection of y onto the column space of x;
-- info is zero unless the computation of b has been requested
-- and r is exactly singular.
-- In this case, info is the index of the first zero
-- diagonal element of r and b is left unaltered.
-- The parameters qy, qty, b, rsd, and xb are not referenced
-- if their computation is not requested and in this case can be
-- replaced by dummy variables in the calling program.
-- To save storage, the user may in some cases use the same array
-- for different parameters in the calling sequence.
-- A frequently occuring example is when one wishes to compute
-- any of b, rsd, or xb and does not need y or qty. In this
-- case one may identify y, qty, and one of b, rsd, or xb, while
-- providing separate arrays for anything else that is to be computed.
-- Thus the calling sequence
--
-- QRLS(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info)
--
-- will result in the computation of b and rsd, with rsd overwriting y.
-- More generally, each item in the following list contains groups of
-- permissible identifications for a single callinng sequence.
--
-- 1. (y,qty,b) (rsd) (xb) (qy)
-- 2. (y,qty,rsd) (b) (xb) (qy)
-- 3. (y,qty,xb) (b) (rsd) (qy)
-- 4. (y,qy) (qty,b) (rsd) (xb)
-- 5. (y,qy) (qty,rsd) (b) (xb)
-- 6. (y,qy) (qty,xb) (b) (rsd)
--
-- In any group the value returned in the array allocated to the group
-- corresponds to the last member of the group.
-- ACKNOWLEDGMENT :
-- This Ada version is a translation of the LINPACK version, datet 08/14/78,
-- written by G.W. Stewart, University of Maryland, Argonne National Lab.
end Standard_Floating_Least_Squares;