Return to standard_floating_least_squares.ads CVS log

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     ! 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;