newmata.txt

矩阵计算库

        2.6  Testing
        ===  =======

The library package contains a comprehensive test program in the form of a
series of files with names of the form tmt?.cxx. The files consist of a large
number of matrix formulae all of which evaluate to zero (except the first one
which is used to check that we are detecting non-zero matrices). The printout
should state that it has found just one non-zero matrix.

Various versions of the make file (extension .mak) are included with the
package. See the section on files [6].

The program also allocates and deletes a large block and small block of memory
before it starts the main testing and then at the end of the test. It then
checks that the blocks of memory were allocated in the same place. If not then
one suspects that there has been a memory leak. i.e. a piece of memory has
been allocated and not deleted.

This is not completely foolproof. Programs may allocate extra print buffers
while the program is running. I have tried to overcome this by doing a print
before I allocate the first memory block. Programs may allocate memory for
different sized items in different places, or might not allocate items
consecutively. Or they might mix the items with memory blocks from other
programs. Nevertheless, I seem to get consistent answers from many of the
compilers I am working with, so I think this is a worthwhile test.

If the DO_FREE_CHECK [2.2] option in include.h is activated, the program
checks that each 'new' is balanced with exactly one 'delete'. This provides a
more definitive test of no memory leaks. There are additional statements in
myexcept.cpp which can be activated to print out details of the memory being
allocated and released.

Several of the files have a line defining 'REPORT' that can be activated
(deactivate the dummy version). This gives a printout of the number of times
each of the 'REPORT' statements in the file is accessed. Use a grep with line
numbers to locate the lines on which 'REPORT' occurs and compare these with
the lines that the printout shows were actually accessed.


        3  Reference manual
        =  ========= ======

        Constructors                                 3.1
        Accessing elements                           3.2
        Assignment and copying                       3.3
        Entering values                              3.4
        Unary operations                             3.5
        Binary operations                            3.6
        Matrix and scalar ops                        3.7
        Scalar functions                             3.8
        Submatrices                                  3.9 
        Change dimension                             3.10 
        Change type                                  3.11 
        Multiple matrix solve                        3.12 
        Memory management                            3.13 
        Efficiency                                   3.14 
        Output                                       3.15 
        Accessing unspecified type                   3.16 
        Cholesky decomposition                       3.17 
        QR decomposition                             3.18 
        Singular value decomposition                 3.19
        Eigenvalue decomposition                     3.20 
        Sorting                                      3.21
        Fast Fourier transform                       3.22 
        Numerical recipes in C                       3.23
        Exceptions                                   3.24 
        Cleanup following exception                  3.25
        Non-linear applications                      3.26 


        3.1  Constructors
        ===  ============

To construct an m x n matrix, 'A', (m and n are integers) use

    Matrix A(m,n);

The UpperTriangularMatrix, LowerTriangularMatrix, SymmetricMatrix and
DiagonalMatrix types are square. To construct an n x n matrix use, for example

    UpperTriangularMatrix UT(n);
    LowerTriangularMatrix LT(n);
    SymmetricMatrix S(n);
    DiagonalMatrix D(n);

Band matrices need to include bandwidth information in their constructors.

    BandMatrix BM(n, lower, upper);
    UpperBandMatrix UB(n, upper);
    LowerBandMatrix LB(n, lower);
    SymmetrixBandMatrix SB(n, lower);

The integers upper and lower are the number of non-zero diagonals above and
below the diagonal (excluding the diagonal) respectively.

The RowVector and ColumnVector types take just one argument in their
constructors:

    RowVector RV(n);
    ColumnVector CV(n);

You can also construct vectors and matrices without specifying the dimension.
For example

    Matrix A;

In this case the dimension must be set by an assignment statement [3.3] or a
re-dimension statement [3.10].

You can also use a constructor to set a matrix equal to another matrix or
matrix expression.

    Matrix A = UT;
    Matrix A = UT * LT;

Only conversions that don't lose information are supported - eg you cannot
convert an upper triangular matrix into a diagonal matrix using =.



        3.2  Accessing elements
        ===  ========= ========

Elements are accessed by expressions of the form 'A(i,j)' where i and j run
from 1 to the appropriate dimension. Access elements of vectors with just one
argument. Diagonal matrices can accept one or two subscripts.

This is different from the earliest version of the package in which the
subscripts ran from 0 to one less than the appropriate dimension. Use
'A.element(i,j)' if you want this earlier convention.

'A(i,j)' and 'A.element(i,j)' can appear on either side of an = sign.

If you activate the '#define SETUP_C_SUBSCRIPTS' in 'include.h' you can also
access elements using the tradition C style notation. That is 'A[i][j]' for
matrices (except diagonal) and 'V[i]' for vectors and diagonal matrices. The
subscripts start at zero (ie like element) and there is no range checking.
Because of the possibility of confusing 'V(i)' and 'V[i]', I suggest you do
not activate this option unless you really want to use it.


        3.3  Assignment and copying
        ===  ========== === =======

The operator = is used for copying matrices, converting matrices, or
evaluating expressions. For example

    A = B;  A = L;  A = L * U;

Only conversions that don't lose information are supported. The dimensions of
the matrix on the left hand side are adjusted to those of the matrix or
expression on the right hand side. Elements on the right hand side which are
not present on the left hand side are set to zero.

The operator << can be used in place of = where it is permissible for
information to be lost.

For example

    SymmetricMatrix S; Matrix A;
    ......
    S << A.t() * A;

is acceptable whereas

    S = A.t() * A;                            // error

will cause a runtime error since the package does not (yet?) recognise
'A.t()*A' as symmetric.

Note that you can not use << with constructors. For example

    SymmetricMatrix S << A.t() * A;           // error

does not work.

Also note that << cannot be used to load values from a full matrix into a band
matrix, since it will be unable to determine the bandwidth of the band matrix.

A third copy routine is used in a similar role to =. Use

    A.Inject(D);

to copy the elements of 'D' to the corresponding elements of 'A' but leave the
elements of 'A' unchanged if there is no corresponding element of 'D' (the =
operator would set them to 0). This is useful, for example, for setting the
diagonal elements of a matrix without disturbing the rest of the matrix.
Unlike = and <<, Inject does not reset the dimensions of 'A', which must match
those of 'D'. Inject does not test for no loss of information.

You cannot replace 'D' by a matrix expression. The effect of 'Inject(D)'
depends on the type of 'D'. If 'D' is an expression it might not be obvious to
the user what type it would have. So I thought it best to disallow
expressions.

Inject can be used for loading values from a regular matrix into a band
matrix. (Don't forget to zero any elements of the left hand side that will not
be set by the loading operation).

Both << and Inject can be used with submatrix expressions on the left hand
side. See the section on submatrices [3.9].

To set the elements of a matrix to a scalar use operator =

    Real r; Matrix A(m,n);
    ......
    Matrix A(m,n); A = r;




        3.4  Entering values
        ===  ======== ======

You can load the elements of a matrix from an array:

    Matrix A(3,2);
    Real a[] = { 11,12,21,22,31,33 };
    A << a;

This construction does not check that the numbers of elements match correctly.
This version of << can be used with submatrices on the left hand side. It is
not defined for band matrices.

Alternatively you can enter short lists using a sequence of numbers separated
by << .

    Matrix A(3,2);
    A << 11 << 12
      << 21 << 22
      << 31 << 32;

This does check for the correct total number of entries, although the message
for there being insufficient numbers in the list may be delayed until the end
of the block or the next use of this construction. This does not work for band
matrices or submatrices, or for long lists. Also try to restrict its use to
numbers. You can include expressions, but these must not call a function which
includes the same construction.



        3.5  Unary operators
        ===  ===== =========

The package supports unary operations

    X = -A      // change sign of elements
    X = A.t()   // transpose
    X = A.i()   // inverse (of square matrix A)




        3.6  Binary operators
        ===  ====== =========

The package supports binary operations

    X = A + B      // matrix addition
    X = A - B      // matrix subtraction
    X = A * B      // matrix multiplication
    X = A.i() * B  // equation solve (square matrix A)
    X = A | B      // concatenate horizontally (concatenate the rows)
    X = A & B      // concatenate vertically (concatenate the columns)
    X = SP(A, B)   // elementwise product of A and B (Schur product)


Notes:

*  If you are doing repeated multiplication. For example 'A*B*C', use brackets
to force the order of evaluation to minimize the number of operations. If 'C'
is a column vector and 'A' is not a vector, then it will usually reduce the
number of operations to use 'A*(B*C)'.

*  In the equation solve example case the inverse is not explicitly
calculated. An LU decomposition of 'A' is performed and this is applied to
'B'. This is more efficient than calculating the inverse and then multiplying.
See also multiple matrix solving [3.12].

*  The package does not (yet?) recognise 'B*A.i()' as an equation solve and
the inverse of 'A' would be calculated. It is probably better to use
'(A.t().i()*B.t()).t()'.

*  Horizontal or vertical concatenation returns a result of type Matrix,
RowVector or ColumnVector.

*  If 'A' is m x p, 'B' is m x q, then 'A | B' is m x (p+q) with the k-th row
being the elements of the k-th row of 'A' followed by the elements of the k-th
row of 'B'.

*  If 'A' is p x n, 'B' is q x n, then 'A & B' is (p+q) x n with the k-th
column being the elements of the k-th column of 'A' followed by the elements
of the k-th column of 'B'.

*  For complicated concatenations of matrices, consider instead using
submatrices [3.9].




        3.7  Matrix and scalar
        ===  ====== === ======

The following expression multiplies the elements of a matrix 'A' by a scalar
f: 'A * f;' Likewise one can divide the elements of a matrix 'A' by a scalar
f:' A / f;'

The expressions  'A + f' and 'A - f' add or subtract a rectangular matrix of
the same dimension as 'A' with elements equal to 'f' to or from the matrix
'A'.

In each case the matrix must be the first term in the expression. Expressions
such  'f + A'  or  'f * A'  are not recognised (yet?).




        3.8  Scalar functions of a matrix
        ===  ====== ========= == = ======


    int m = A.Nrows();                    // number of rows
    int n = A.Ncols();                    // number of columns
    Real ssq = A.SumSquare();             // sum of squares of elements
    Real sav = A.SumAbsoluteValue();      // sum of absolute values
    Real s = A.Sum();                     // sum of values
    Real mav = A.MaximumAbsoluteValue();  // maximum of absolute values
    Real norm = A.Norm1();                // maximum of sum of absolute
                                             values of elements of a column
    Real norm = A.NormInfinity();         // maximum of sum of absolute
                                             values of elements of a row
    Real t = A.Trace();                   // trace
    LogandSign ld = A.LogDeterminant();   // log of determinant
    Boolean z = A.IsZero();               // test all elements zero
    MatrixType mt = A.Type();             // type of matrix
    Real* s = Store();                    // pointer to array of elements
    int l = Storage();                    // length of array of elements

'A.LogDeterminant()' returns a value of type LogandSign. If ld is of type
LogAndSign  use

    ld.Value()    to get the value of the determinant
    ld.Sign()     to get the sign of the determinant (values 1, 0, -1)