newmata.txt

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

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.


Accessing matrices of unspecified type
--------------------------------------

Skip this section on your first reading.

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
   }


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


Householder triangularisation
-----------------------------

Start with matrix

       / X    0 \      s
       \ Y    0 /      t

         n    s

The Householder triangularisation post multiplies by an orthogonal
matrix Q such that the matrix becomes

       / 0    L \      s
       \ Z    M /      t

         n    s

where L is lower triangular. Note that X is the transpose of the matrix
sometimes considered in this context.

This is good for solving least squares problems: choose b (matrix or row
vector) to minimize the sum of the squares of the elements of

         Y - b*X

Then choose b = M * L.i();

Two routines are provided:

    HHDecompose(X, L);

replaces X by orthogonal columns and forms L.

    HHDecompose(X, Y, M);

uses X from the first routine, replaces Y by Z and forms M.


Singular Value Decomposition
----------------------------

The singular value decomposition of an m x n matrix A ( where m >= n) is
a decomposition

    A  = U * D * V.t()

where U is m x n with  U.t() * U  equalling the identity, D is an n x n
diagonal matrix and V is an n x n orthogonal matrix.

Singular value decompositions are useful for understanding the structure
of ill-conditioned matrices, solving least squares problems, and for
finding the eigenvalues of A.t() * A.

To calculate the singular value decomposition of A (with m >= n) use one
of

    SVD(A, D, U, V);                  // U (= A is OK)
    SVD(A, D);
    SVD(A, D, U);                     // U (= A is OK)
    SVD(A, D, U, FALSE);              // U (can = A) for workspace only
    SVD(A, D, U, V, FALSE);           // U (can = A) for workspace only

The values of A are not changed unless A is also inserted as the third
argument.


Eigenvalues
-----------

An eigenvalue decomposition of a symmetric matrix A is a decomposition

    A  = V * D * V.t()

where V is an orthogonal matrix and D is a diagonal matrix.

Eigenvalue analyses are used in a wide variety of engineering,
statistical and other mathematical analyses.

The package includes two algorithms: Jacobi and Householder. The first
is extremely reliable but much slower than the second.

The code is adapted from routines in "Handbook for Automatic
Computation, Vol II, Linear Algebra" by Wilkinson and Reinsch, published
by Springer Verlag. 


    Jacobi(A,D,S,V);                  // A, S symmetric; S is workspace,
                                      //    S = A is OK
    Jacobi(A,D);                      // A symmetric
    Jacobi(A,D,S);                    // A, S symmetric; S is workspace,
                                      //    S = A is OK
    Jacobi(A,D,V);                    // A symmetric

    EigenValues(A,D);                 // A symmetric
    EigenValues(A,D,S);               // A, S symmetric; S is for back
                                      //    transforming, S = A is OK
    EigenValues(A,D,V);               // A symmetric


Sorting
-------

To sort the values in a matrix or vector, A, (in general this operation
makes sense only for vectors and diagonal matrices) use

    SortAscending(A);

or

    SortDescending(A);

I use the Shell-sort algorithm. This is a medium speed algorithm, you
might want to replace it with something faster if speed is critical and
your matrices are large.


Fast Fourier Transform
----------------------

FFT(CV1, CV2, CV3, CV4);       // CV3=CV1 and CV4=CV2 is OK

where CV1, CV2, CV3, CV4 are column vectors. CV1 and CV2 are the real
and imaginary input vectors; CV3 and CV4 are the real and imaginary
output vectors. The lengths of CV1 and CV2 must be equal and should be
the product of numbers less than about 10 for fast execution.


Interface to Numerical Recipes in C
-----------------------------------

This package can be used with the vectors and matrices defined in
"Numerical Recipes in C". You need to edit the routines in Numerical
Recipes so that the elements are of the same type as used in this
package. Eg replace float by double, vector by dvector and matrix by
dmatrix, etc. You will also need to edit the function definitions to use
the version acceptable to your compiler. Then enclose the code from
Numerical Recipes in  extern "C" { ... }. You will also need to include
the matrix and vector utility routines.

Then any vector in Numerical Recipes with subscripts starting from 1 in
a function call can be accessed by a RowVector, ColumnVector or
DiagonalMatrix in the present package. Similarly any matrix with
subscripts starting from 1 can be accessed by an  nricMatrix  in the
present package. The class nricMatrix is derived from Matrix and can be
used in place of Matrix. In each case, if you wish to refer to a
RowVector, ColumnVector, DiagonalMatrix or nricMatrix X in an function
from Numerical Recipes, use X.nric() in the function call.

Numerical Recipes cannot change the dimensions of a matrix or vector. So
matrices or vectors must be correctly dimensioned before a Numerical
Recipes routine is called.

For example

   SymmetricMatrix B(44);
   .....                             // load values into B
   nricMatrix BX = B;                // copy values to an nricMatrix
   DiagonalMatrix D(44);             // Matrices for output
   nricMatrix V(44,44);              //    correctly dimensioned
   int nrot;
   jacobi(BX.nric(),44,D.nric(),V.nric(),&nrot);
                                     // jacobi from NRIC
   cout << D;                        // print eigenvalues


Exceptions
----------

This package includes a partial implementation of exceptions. I used
Carlos Vidal's article in the September 1992 C Users Journal as a
starting point.

Newmat does a partial clean up of memory following throwing an exception
- see the next section. However, the present version will leave a little
heap memory unrecovered under some circumstances. I would not expect
this to be a problem in most circumstances, but it is something that
needs to be sorted out.

The functions/macros I define are Try, Throw, Catch, CatchAll and
CatchAndThrow. Try, Throw, Catch and CatchAll correspond to try, throw,
catch and catch(...) in the C++ standard. A list of Catch clauses must
be terminated by either CatchAll or CatchAndThrow but not both. Throw
takes an Exception as an argument or takes no argument (for passing on
an exception). I do not have a version of Throw for specifying which
exceptions a function might throw. Catch takes an exception class name
as an argument; CatchAll and CatchAndThrow don't have any arguments.
Try, Catch and CatchAll must be followed by blocks enclosed in curly
brackets.

All Exceptions must be derived from a class, Exception, defined in
newmat and can contain only static variables. See the examples in newmat
if you want to define additional exceptions.

I have defined 5 clases of exceptions for users (there are others but I
suggest you stick to these ones):

   SpaceException                 Insufficient space on the heap
   ProgramException               Errors such as out of range index or
                                  incompatible matrix types or
                                  dimensions
   ConvergenceException           Iterative process does not converge
   DataException                  Errors such as attempting to invert a
                                  singular matrix
   InternalException              Probably a programming error in newmat

For each of these exception classes, I have defined a member function
void SetAction(int). If you call SetAction(1), and a corresponding
exception occurs, you will get an error message. If there is a Catch
clause for that exception, execution will be passed to that clause,
otherwise the program will exit. If you call SetAction(0) you will get
the same response, except that there will be no error message. If you
call SetAction(-1), you will get the error message but the program will
always exit.

I have defined a class Tracer that is intended to help locate the place
where an error has occurred. At the beginning of a function I suggest
you include a statement like

   Tracer tr("name");

where name is the name of the function. This name will be printed as
part of the error message, if an exception occurs in that function, or
in a function called from that function. You can change the name as you
proceed through a function with the ReName function

   tr.ReName("new name");

if, for example, you want to track progress through the function.


Clean up following an exception
-------------------------------

The exception mechanisms in newmat are based on the C functions setjmp
and longjmp. These functions do not call destructors so can lead to
garbage being left on the heap. (I refer to memory allocated by "new" as
heap memory). For example, when you call

Matrix A(20,30);

a small amount of space is used on the stack containing the row and
column dimensions of the matrix and 600 doubles are allocated on the
heap for the actual values of the matrix. At the end of the block in
which A is declared, the destructor for A is called and the 600
doubles are freed. The locations on the stack are freed as part of the
normal operations of the stack. If you leave the block using a longjmp
command those 600 doubles will not be freed and will occupy space until
the program terminates.

To overcome this problem newmat keeps a list of all the currently
declared matrices and its exception mechanism will return heap memory
when you do a Throw and Catch.

However it will not return heap memory from objects from other packages.
If you want the mechanism to work with another class you will have to do
three things:

1: derive your class from class Janitor defined in except.h;

2: define a function void CleanUp() in that class to return all heap
   memory;

3: include the following lines in the class definition

      public:
         void* operator new(size_t size)
         { do_not_link=TRUE; void* t = ::operator new(size); return t; }






---------------------------------------------------------------------------


List of files
=============

README          readme file
NEWMATA  TXT    documentation file
NEWMATB  TXT    notes on the package design
NEWMATC  TXT    notes on testing the package

BOOLEAN  H      boolean class definition
CONTROLW H      control word definition file
EXCEPT   H      general exception handler definitions
INCLUDE  H      details of include files and options
NEWMAT   H      main matrix class definition file
NEWMATAP H      applications definition file
NEWMATIO H      input/output definition file
NEWMATRC H      row/column functions definition files
NEWMATRM H      rectangular matrix access definition files
PRECISIO H      numerical precision constants

BANDMAT  CXX    band matrix routines
CHOLESKY CXX    Cholesky decomposition
ERROR    CXX    general error handler
EVALUE   CXX    eigenvalues and eigenvector calculation
FFT      CXX    fast Fourier transform
HHOLDER  CXX    Householder triangularisation
JACOBI   CXX    eigenvalues by the Jacobi method
NEWMAT1  CXX    type manipulation routines
NEWMAT2  CXX    row and column manipulation functions
NEWMAT3  CXX    row and column access functions
NEWMAT4  CXX    constructors, redimension, utilities
NEWMAT5  CXX    transpose, evaluate, matrix functions
NEWMAT6  CXX    operators, element access
NEWMAT7  CXX    invert, solve, binary operations
NEWMAT8  CXX    LU decomposition, scalar functions
NEWMAT9  CXX    output routines
NEWMATEX CXX    matrix exception handler
NEWMATRM CXX    rectangular matrix access functions
SORT     CXX    sorting functions
SUBMAT   CXX    submatrix functions
SVD      CXX    singular value decomposition

EXAMPLE  CXX    example of use of package
EXAMPLE  TXT    output from example
EXAMPLE  MAK    make file for example

---------------------------------------------------------------------------

                   Matrix package problem report form
                   ----------------------------------

Version: ...............newmat06
Date of release: .......December 1st, 1992
Primary site: ..........
Downloaded from: .......
Your email address: ....
Today's date: ..........
Your machine: ..........
Operating system: ......
Compiler & version: ....
Describe the problem - attach examples if possible:









Email to  [e-mail]  or  Compuserve 72777,656 

-------------------------------------------------------------------------------