newmata.txt

科学和工程计算中使用统计功能开发工具包

files in the interactive environment by pointing and clicking.


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 U(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);

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 or a
re-dimension statement.

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

    Matrix A = U;

    Matrix A = U * L;

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


Elements of matrices
--------------------

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.


Matrix copy
-----------

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.

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

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

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 cannot check that the numbers of elements match
correctly. This version of << can be used with submatrices on the left
hand side.


Unary operators
---------------

The package supports unary operations

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


Binary operations
-----------------

The package supports binary operations

    matrix addition                    A+B
    matrix subtraction                 A-B
    matrix multiplication              A*B
    equation solve (square matrix A)   A.i()*B

In the last case the inverse is not calculated.

Notes:

If you are doing repeated multiplication. For example A*B*C, use
brackets to force the order 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) .

The package does not recognise B*A.i() as an equation solve. It is
probably better to use (A.t().i()*B.t()).t() .


Combination of a 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.


Scalar functions of matrices
----------------------------
            
    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 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)
    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, RowV, ColV, Crout, BndLU)
showing the type.

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


Submatrix operations
--------------------

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  << , 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.


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.


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.


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.


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++ 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 you are using AT&T C++ you may wish to replace  return A; by
return (ReturnMatrix)A;  to avoid a warning message.

 
 --------------------------------------------------------------------- 
| 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.