newmata.txt

矩阵计算库

    ld.LogValue() to get the log of the absolute value.

'A.IsZero()' returns Boolean value TRUE if the matrix 'A' has all elements
equal to 0.0.

'MatrixType mt = A.Type()' returns the type of a matrix. Use '(char*)mt' to
get a string  (UT, LT, Rect, Sym, Diag, Band, UB, LB, Crout, BndLU) showing
the type (Vector types are returned as Rect).

The versions Sum(A), SumSquare(A), SumAbsoluteValue(A),
MaximumAbsoluteValue(A), Trace(A), LogDeterminant(A), Norm1(A),
NormInfinity(A)  can be used in place of A.Sum(), A.SumSquare(),
A.SumAbsoluteValue(), A.MaximumAbsoluteValue(), A.Trace(), A.LogDeterminant(),
A.Norm1(), A.NormInfinity().



        3.9  Submatrices
        ===  ===========


    A.SubMatrix(fr,lr,fc,lc)

This selects a submatrix from 'A'. The arguments  fr,lr,fc,lc are the first
row, last row, first column, last column of the submatrix with the numbering
beginning at 1. This may be used in any matrix expression or on the left hand
side of =, << or Inject. Inject does not check no information loss. You can
also use the construction

    Real c; .... A.SubMatrix(fr,lr,fc,lc) = c;

to set a submatrix equal to a constant.

The follwing are variants of SubMatrix:

    A.SymSubMatrix(f,l)             //   This assumes fr=fc and lr=lc.
    A.Rows(f,l)                     //   select rows
    A.Row(f)                        //   select single row
    A.Columns(f,l)                  //   select columns
    A.Column(f)                     //   select single column

In each case f and l mean the first and last row or column to be selected
(starting at 1).

If SubMatrix or its variant occurs on the right hand side of an = or << or
within an expression its type is as follows

    A.Submatrix(fr,lr,fc,lc):           If A is RowVector or
                                        ColumnVector then same type
                                        otherwise type Matrix
    A.SymSubMatrix(f,l):                Same type as A
    A.Rows(f,l):                        Type Matrix
    A.Row(f):                           Type RowVector
    A.Columns(f,l):                     Type Matrix
    A.Column(f):                        Type ColumnVector

If SubMatrix or its variant appears on the left hand side of  = or << , think
of its type being Matrix. Thus 'L.Row(1)' where 'L' is LowerTriangularMatrix
expects 'L.Ncols()' elements even though it will use only one of them. If you
are using = the program will check for no loss of data.

If you are are using the submatrix facility to build a matrix from a small
number of components, consider instead using the concatenation operators
[3.6].



        3.10  Change dimensions
        ====  ====== ==========

The following operations change the dimensions of a matrix. The values of the
elements are lost.

    A.ReDimension(nrows,ncols);     // for type Matrix or nricMatrix
    A.ReDimension(n);               // for all other types, except Band
    A.ReDimension(n,lower,upper);   // for BandMatrix
    A.ReDimension(n,lower);         // for LowerBandMatrix
    A.ReDimension(n,upper);         // for UpperBandMatrix
    A.ReDimension(n,lower);         // for SymmetricBandMatrix

Use   'A.CleanUp()'  to set the dimensions of 'A' to zero and release all the
heap memory.

Remember that 'ReDimension' destroys values. If you want to 'ReDimension', but
keep the values in the bit that is left use something like

   ColumnVector V(100);
   ...                            // load values
   V = V.Rows(1,50);              // to get first 50 values.

If you want to extend a matrix or vector use something like

   ColumnVector V(50);
   ...                            // load values
   { V.Release(); ColumnVector X=V; V.ReDimension(100); V.Rows(1,50)=X; }
                                  // V now length 100




        3.11  Change type
        ====  ====== ====

The following functions interpret the elements of a matrix (stored row by row)
to be a vector or matrix of a different type. Actual copying is usually
avoided where these occur as part of a more complicated expression.

    A.AsRow()
    A.AsColumn()
    A.AsDiagonal()
    A.AsMatrix(nrows,ncols)
    A.AsScalar()

The expression 'A.AsScalar()' is used to convert a 1 x 1 matrix to a scalar.



        3.12  Multiple matrix solve
        ====  ======== ====== =====

If 'A' is a square or symmetric matrix use

    CroutMatrix X = A;                // carries out LU decomposition
    Matrix AP = X.i()*P; Matrix AQ = X.i()*Q;
    LogAndSign ld = X.LogDeterminant();

rather than

    Matrix AP = A.i()*P; Matrix AQ = A.i()*Q;
    LogAndSign ld = A.LogDeterminant();

since each operation will repeat the LU decomposition.

If 'A' is a BandMatrix or a SymmetricBandMatrix begin with

    BandLUMatrix X = A;               // carries out LU decomposition

A CroutMatrix or a BandLUMatrix can't be manipulated or copied. Use references
as an alternative to copying.

Alternatively use

    LinearEquationSolver X = A;

This will choose the most appropiate decomposition of 'A'. That is, the band
form if 'A' is banded; the Crout decomposition if 'A' is square or symmetric
and no decomposition if 'A' is triangular or diagonal. If you want to use the
LinearEquationSolver '#include newmatap.h'.



        3.13  Memory management
        ====  ====== ==========

The package does not support delayed copy. Several strategies are required to
prevent unnecessary matrix copies.

Where a matrix is called as a function argument use a constant reference. For
example

    YourFunction(const Matrix& A)

rather than

    YourFunction(Matrix A)


Skip the rest of this section on your first reading.
 ------------------------------------------------------------------------------
Gnu g++ (< 2.6) users please read on; if you are returning matrix values from
a function, then you must use the ReturnMatrix construct.     
 ------------------------------------------------------------------------------

A second place where it is desirable to avoid unnecessary copies is when a
function is returning a matrix. Matrices can be returned from a function with
the return command as you would expect. However these may incur one and
possibly two copyings of the matrix. To avoid this use the following
instructions.

Make your function of type  ReturnMatrix . Then precede the return statement
with a Release statement (or a ReleaseAndDelete statement if the matrix was
created with new). For example


    ReturnMatrix MakeAMatrix()
    {
       Matrix A;
       ......
       A.Release(); return A;
    }

or

    ReturnMatrix MakeAMatrix()
    {
       Matrix* m = new Matrix;
       ......
       m->ReleaseAndDelete(); return *m;
    }

If your compiler objects to this code, replace the return statements with

    return A.ForReturn();

or

    return m->ForReturn();


If you are using AT&T C++ you may wish to replace  'return A;' by return
'(ReturnMatrix)A;'  to avoid a warning message; but this will give a runtime
error with Gnu. (You can't please everyone.)
 ------------------------------------------------------------------------------
Do not forget to make the function of type ReturnMatrix; otherwise you may get
incomprehensible run-time errors.                       
 ------------------------------------------------------------------------------
You can also use '.Release()' or '->ReleaseAndDelete()' to allow a matrix
expression to recycle space. Suppose you call

    A.Release();

just before 'A' is used just once in an expression. Then the memory used by
'A' is either returned to the system or reused in the expression. In either
case, 'A''s memory is destroyed. This procedure can be used to improve
efficiency and reduce the use of memory.

Use '->ReleaseAndDelete' for matrices created by new if you want to completely
delete the matrix after it is accessed.



        3.14  Efficiency
        ====  ==========

The package tends to be not very efficient for dealing with matrices with
short rows. This is because some administration is required for accessing rows
for a variety of types of matrices. To reduce the administration a special
multiply routine is used for rectangular matrices in place of the generic one.
Where operations can be done without reference to the individual rows (such as
adding matrices of the same type) appropriate routines are used.

When you are using small matrices (say smaller than 10 x 10) you may find it
faster to use rectangular matrices rather than the triangular or symmetric
ones.



        3.15  Output
        ====  ======

To print a matrix use an expression like

    Matrix A;
    ......
    cout << setw(10) << setprecision(5) << A;

This will work only with systems that support the AT&T input/output routines
including manipulators. You need to #include the file newmat9.h.

The present version of this routine is useful only for matrices small enough
to fit within a page or screen width.

To print several vectors or matrices in columns use a concatenation operator
[3.6]:

   Matrix A, B;
   .....
   cout << setw(10) << setprecision(5) << (A | B);




        3.16  Unspecified type
        ====  =========== ====

Skip this section on your first reading.

If you want to work with a matrix of unknown type, say in a function. You can
construct a matrix of type 'GenericMatrix'. Eg

   Matrix A;
   .....                                  // put some values in A
   GenericMatrix GM = A;

A GenericMatrix matrix can be used anywhere where a matrix expression can be
used and also on the left hand side of an '='. You can pass any type of matrix
(excluding the Crout and BandLUMatrix types) to a 'const GenericMatrix&'
argument in a function. However most scalar functions including Nrows(),
Ncols(), Type() and element access do not work with it. Nor does the
ReturnMatrix construct. See also the paragraph on LinearEquationSolver [3.12].

An alternative and less flexible approach is to use Basematrix or
GeneralMatrix.

Suppose you wish to write a function which accesses a matrix of unknown type
including expressions (eg 'A*B'). Then use a layout similar to the following:

   void YourFunction(BaseMatrix& X)
   {
      GeneralMatrix* gm = X.Evaluate();   // evaluate an expression
                                          // if necessary
      ........                            // operations on *gm
      gm->tDelete();                      // delete *gm if a temporary
   }

See, as an example, the definitions of 'operator<<' in newmat9.cxx.

Under certain circumstances; particularly where 'X' is to be used just once in
an expression you can leave out the 'Evaluate()' statement and the
corresponding 'tDelete()'. Just use 'X' in the expression.

If you know YourFunction will never have to handle a formula as its argument
you could also use

   void YourFunction(const GeneralMatrix& X)
   {
      ........                            // operations on X
   }

Do not try to construct a GeneralMatrix or BaseMatrix.



        3.17  Cholesky decomposition
        ====  ======== =============

Suppose S is symmetric and positive definite. Then there exists a unique lower
triangular matrix L such that L * L.t() = S. To calculate this use

    SymmetricMatrix S;
    ......
    LowerTriangularMatrix L = Cholesky(S);

If S is a symmetric band matrix then L is a band matrix and an alternative
procedure is provied for carrying out the decomposition:

    SymmetricBandMatrix S;
    ......
    LowerBandMatrix L = Cholesky(S);




        3.18  QR decomposition