[BACK]Return to standard_floating_least_squares.ads CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices

Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Matrices/standard_floating_least_squares.ads, Revision 1.1.1.1

1.1       maekawa     1: with Standard_Floating_Vectors;
                      2: with Standard_Floating_Matrices;
                      3:
                      4: package Standard_Floating_Least_Squares is
                      5:
                      6: -- DESCRIPTION :
                      7: --   This package provides an implementation of Least Squares algorithm
                      8: --   for matrices of standard floating-point numbers.
                      9:
                     10:   procedure QRLS ( x : in out Standard_Floating_Matrices.Matrix;
                     11:                    ldx,n,k : in natural;
                     12:                    qraux,y : in Standard_Floating_Vectors.Vector;
                     13:                    qy,qty,b,rsd,xb : out Standard_Floating_Vectors.Vector;
                     14:                    job : in natural; info : out natural );
                     15:
                     16:   -- DESCRIPTION :
                     17:   --   Applies the output of Standard_Floating_QR_Decomposition.QRD to
                     18:   --   compute coordinate transformations, projections and least quares
                     19:   --   solutions.  For k <= min(n,p), let xk be the matrix
                     20:   --
                     21:   --             xk = (x(jpvt(1)),x(jpvt(2)), ... ,x(jpvt(k)))
                     22:   --
                     23:   --   formed from columnns jpvt(1), ... ,jpvt(k) of the original
                     24:   --   n x p matrix x that was input to dqrdc (if no pivoting was done,
                     25:   --   xk consists of the first k columns of x in their original order).
                     26:   --   QRD produces a factored orthogonal matrix q and an upper triangular
                     27:   --   matrix r such that
                     28:   --                       xk = q * (r)
                     29:   --                                (0)
                     30:   --   this information is contained in coded form in x and qraux.
                     31:
                     32:   -- ON ENTRY :
                     33:   --   x         contains the output of QRD, of size ldx times p;
                     34:   --   ldx       leading dimension of x;
                     35:   --   n         number of rows in the matrix xk, must be same as in QRD
                     36:   --   k         number of columns of the matrix xk, k <= min(n,p)
                     37:   --   qraux     contains p entries, auxiliary output from QRD;
                     38:   --   y         n-vector to be manipulated by QRLS;
                     39:   --   job       specifies what is to be computed, job has the decimal
                     40:   --             expansion abcde, with the following meaning :
                     41:   --                  if a /= 0, compute qy,
                     42:   --                  if b,c,d, or e /= 0, compute qty,
                     43:   --                  if c /= 0, compute b,
                     44:   --                  if d /= 0, compute rsd,
                     45:   --                  if e /= 0, compute xb.
                     46:   --             Note that a request to compute b, rsd, or xb
                     47:   --             automatically triggers the computation of qty, for
                     48:   --             which an array must be provided in the calling sequence.
                     49:
                     50:   -- ON RETURN :
                     51:   --   x         may be altered, used as work space;
                     52:   --   qy        qy'length = n,
                     53:   --             qy conntains q*y, if its computation has been requested;
                     54:   --   qty       qty'length = n,
                     55:   --             qty contains trans(q)*y, if its computation has been
                     56:   --             requested;  here trans(q) is the transpose of the matrix q;
                     57:   --   b         b'length = k,
                     58:   --             b contains the solution of the least squares problem
                     59:   --
                     60:   --                    minimize norm2(y - xk*b),
                     61:   --
                     62:   --             if its computation has been requested.  (Note that
                     63:   --             if pivoting was requested in dqrdc, the j-th
                     64:   --             component of b will be associated with column jpvt(j)
                     65:   --             of the original matrix x that was input into dqrdc.)
                     66:   --   rsd       rsd'length = n,
                     67:   --             rsd contains the least squares residual y - xk*b,
                     68:   --             if its computation has been requested;  rsd is also the
                     69:   --             orthogonal projection of y onto the orthogonal complement
                     70:   --             of the column space of xk;
                     71:   --   xb        x'length = n, xb contains the least squares approximation
                     72:   --             xk*b, if its computation has been requested;  xb is also
                     73:   --             the orthogonal projection of y onto the column space of x;
                     74:   --   info      is zero unless the computation of b has been requested
                     75:   --             and r is exactly singular.
                     76:   --             In this case, info is the index of the first zero
                     77:   --             diagonal element of r and b is left unaltered.
                     78:
                     79:   -- The parameters qy, qty, b, rsd, and xb are not referenced
                     80:   -- if their computation is not requested and in this case can be
                     81:   -- replaced by dummy variables in the calling program.
                     82:   -- To save storage, the user may in some cases use the same array
                     83:   -- for different parameters in the calling sequence.
                     84:   -- A frequently occuring example is when one wishes to compute
                     85:   -- any of b, rsd, or xb and does not need y or qty.  In this
                     86:   -- case one may identify y, qty, and one of b, rsd, or xb, while
                     87:   -- providing separate arrays for anything else that is to be computed.
                     88:   -- Thus the calling sequence
                     89:   --
                     90:   --     QRLS(x,ldx,n,k,qraux,y,dum,y,b,y,dum,110,info)
                     91:   --
                     92:   -- will result in the computation of b and rsd, with rsd overwriting y.
                     93:   -- More generally, each item in the following list contains groups of
                     94:   -- permissible identifications for a single callinng sequence.
                     95:   --
                     96:   --     1. (y,qty,b) (rsd) (xb) (qy)
                     97:   --     2. (y,qty,rsd) (b) (xb) (qy)
                     98:   --     3. (y,qty,xb) (b) (rsd) (qy)
                     99:   --     4. (y,qy) (qty,b) (rsd) (xb)
                    100:   --     5. (y,qy) (qty,rsd) (b) (xb)
                    101:   --     6. (y,qy) (qty,xb) (b) (rsd)
                    102:   --
                    103:   -- In any group the value returned in the array allocated to the group
                    104:   -- corresponds to the last member of the group.
                    105:
                    106:   -- ACKNOWLEDGMENT :
                    107:   --   This Ada version is a translation of the LINPACK version, datet 08/14/78,
                    108:   --   written by G.W. Stewart, University of Maryland, Argonne National Lab.
                    109:
                    110: end Standard_Floating_Least_Squares;

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>