newmatb.txt

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

Alternatively you can calculate the address using the formula
appropriate for the structure of A. I use this second approach. It is
probably slower, but saves worrying about an extra bit of storage. The
other question is whether to check for i and j being in range. I do
carry out this check following years of experience with both systems
that do and systems that don't do this check.

I would hope that the routines I supply with this package will reduce
your need to access elements of matrices so speed of access is not a
high priority.


Use iterators?
--------------

Iterators are an alternative way of providing fast access to the
elements of an array or matrix when they are to be accessed
sequentially. They need to be customised for each type of matrix. I have
not implemented iterators in this package, although some iterator like
functions are used for some row and column functions.


Memory management - reference counting or status variable?
----------------------------------------------------------

Consider the instruction

   X = A + B + C;

To evaluate this a simple program will add A to B putting the total in a
temporary T1. Then it will add T1 to C creating another temporary T2
which will be copied into X. T1 and T2 will sit around till the end of
the block. It would be faster if the program recognised that T1 was
temporary and stored the sum of T1 and C back into T1 instead of
creating T2 and then avoided the final copy by just assigning the
contents of T1 to X rather than copying. In this case there will be no
temporaries requiring deletion. (More precisely there will be a header
to be deleted but no contents).

For an instruction like

   X = (A * B) + (C * D);

we can't easily avoid one temporary being left over, so we would like
this temporary deleted as quickly as possible.

I provide the functionality for doing this by attaching a status
variable to each matrix. This indicates if the matrix is temporary so
that its memory is available for recycling or deleting. Any matrix
operation checks the status variables of the matrices it is working with
and recycles or deletes any temporary memory.

An alternative or additional approach would be to use delayed copying.
If a program requests a matrix to be copied, the copy is delayed until
an instruction is executed which modifies the memory of either the
original matrix or the copy. This saves the unnecessary copying in the
previous examples. However, it does not provide the additional
functionality of my approach.

It would be possible to have both delayed copy and tagging temporaries,
but this seemed an unnecessary complexity. In particular delayed copy
mechanisms seem to require two calls to the heap manager, using extra
time and making building a memory compacting mechanism more difficult.


Evaluation of expressions - use two stage method?
-------------------------------------------------

Consider the instruction

   X = B - X;

The simple program will subtract X from B, store the result in a
temporary T1 and copy T1 into X. It would be faster if the program
recognised that the result could be stored directly into X. This would
happen automatically if the program could look at the instruction first
and mark X as temporary.

C programmers would expect to avoid the same problem with

   X = X - B;

by using an operator -= (which I haven't provided, yet)

   X -= B;

However this is an unnatural notation for non C users and it is much
nicer to write

   X = X - B;

and know that the program will carry out the simplification.

Another example where this intelligent analysis of an instruction is
helpful is in

   X = A.i() * B;

where i() denotes inverse. Numerical analysts know it is inefficient to
evaluate this expression by carrying out the inverse operation and then
the multiply. Yet it is a convenient way of writing the instruction. It
would be helpful if the program recognised this expression and carried
out the more appropriate approach.

To carry out this "intelligent" analysis of an instruction  matrix
expressions are evaluated in two stages. In the the first stage a tree
representation of the expression is formed.

For example (A+B)*C is represented by a tree

                    *
                   / \
                  +   C
                 / \
                A   B

Rather than adding A and B the + operator yields an object of a class
"AddedMatrix" which is just a pair of pointers to A and B. Then the *
operator yields a "MultipliedMatrix" which is a pair of pointers to the
"AddedMatrix" and C. The tree is examined for any simplifications and
then evaluated recursively.

Further possibilities not yet included are to recognise A.t()*A and
A.t()+A as symmetric or to improve the efficiency of evaluation of
expressions like A+B+C, A*B*C, A*B.t()  [t() denotes transpose].

One of the disadvantages of the two-stage approach is that the types of
matrix expressions are determined at run-time. So the compiler will not
detect errors of the type

   Matrix M; DiagonalMatrix D; ....; D = M;

We don't allow conversions using = when information would be lost. Such
errors will be detected when the statement is executed.


How to overcome an explosion in number of operations
----------------------------------------------------

The package attempts to solve the problem of the large number of
versions of the binary operations required when one has a variety of
types. With n types of matrices the binary operations will each require
n-squared separate algorithms. Some reduction in the number may be
possible by carrying out conversions. However the situation rapidly
becomes impossible with more than 4 or 5 types.

Doug Lea told me that it was possible to avoid this problem. I don't
know what his solution is. Here's mine.

Each matrix type includes routines for extracting individual rows or
columns. I assume a row or column consists of a sequence of zeros, a
sequence of stored values and then another sequence of zeros. Only a
single algorithm is then required for each binary operation. The rows
can be located very quickly since most of the matrices are stored row by
row. Columns must be copied and so the access is somewhat slower. As far
as possible my algorithms access the matrices by row.

An alternative approach of using iterators will be slower since the
iterators will involve a virtual function access at each step.

There is another approach. Each of the matrix types defined in this
package can be set up so both rows and columns have their elements at
equal intervals provided we are prepared to store the rows and columns
in up to three chunks. With such an approach one could write a single
"generic" algorithm for each of multiply and add. This would be a
reasonable alternative to my approach.

I provide several algorithms for operations like + . If one is adding
two matrices of the same type then there is no need to access the
individual rows or columns and a faster general algorithm is
appropriate.

Generally the method works well. However symmetric matrices are not
always handled very efficiently (yet) since complete rows are not stored
explicitly.

The original version of the package did not use this access by row or
column method and provided the multitude of algorithms for the
combination of different matrix types. The code file length turned out
to be just a little longer than the present one when providing the same
facilities with 5 distinct types of matrices. It would have been very
difficult to increase the number of matrix types in the original
version. Apparently 4 to 5 types is about the break even point for
switching to the approach adopted in the present package.


Destruction of temporaries
--------------------------

Versions before version 5 of newmat did not work correctly with Gnu C++.
This was because the tree structure used to represent a matrix
expression was set up on the stack. This was fine for AT&T, Borland and
Zortech C++. However Gnu C++ destroys temporary structures as soon as
the function that accesses them finishes. The other compilers wait until
the end of the current block. (It would be good enough for them to wait
till the end of the expression.) To overcome this problem, there is now
an option to store the temporaries forming the tree structure on the
heap (created with new) and to delete them explicitly. Activate the
definition of TEMPS_DESTROYED_QUICKLY to set this option. In fact, I
suggest this be the default option as, with it, the package uses less
storage and runs faster.

There still exist situations Gnu C++ will go wrong. These include
statements like

A = X * Matrix(P * Q); Real r = (A*B)(3,4);

Neither of these kinds of statements will occur often in practice.


A calculus of matrix types
--------------------------

The program needs to be able to work out the class of the result of a
matrix expression. This is to check that a conversion is legal or to
determine the class of a temporary. To assist with this, a class
MatrixType is defined. Operators +, -, *, >= are defined to calculate
the types of the results of expressions or to check that conversions are
legal.


Error handling
--------------

The package now does have a moderately graceful exit from errors.
Originally I thought I would wait until exceptions became available in
C++. Trying to set up an error handling scheme that did not involve an
exception-like facility was just too complicated. Version 5 of this
package included the beginnings of such a facility. The approach in the
present package, attempting to implement C++ exceptions, is not
completely satisfactory, but seems a good interim solution.

The exception mechanism cannot clean-up objects explicitly created by
new. This must be explicitly carried out by the package writer or the
package user. I have not yet done this with the present package so
occasionally a little garbage may be left behind after an exception. I
don't think this is a big problem, but it is one that needs fixing.

I identify five categories of errors:

   Programming error - eg illegal conversion of a matrix, subscript out
   of bounds, matrix dimensions don't match;

   Illegal data error - eg Cholesky of a non-positive definite matrix;

   Out of space error - "new" returns a null pointer;

   Convergence failure - an iterative operation fails to converge;

   Internal error - failure of an internal check.

Some of these have subcategories, which are nicely represented by
derived exception classes. But I don't think it is a good idea for
package users to refer to these as they may change in the next release
of the package.


Sparse matrices
---------------

The package does not yet support sparse matrices.

For sparse matrices there is going to be some kind of structure vector.
It is going to have to be calculated for the results of expressions in
much the same way that types are calculated. In addition, a whole new
set of row and column operations would have to be written.

Sparse matrices are important for people solving large sets of
differential equations as well as being important for statistical and
operational research applications. I think comprehensive matrix package
needs to implement sparse matrices.


Complex matrices
----------------

The package does not yet support matrices with complex elements. There
are at least two approaches to including these. One is to have matrices
with complex elements. This probably means making new versions of the
basic row and column operations for Real*Complex, Complex*Complex,
Complex*Real and similarly for + and -. This would be OK, except that I
am also going to have to also do this for sparse and when you put these
together, the whole thing will get out of hand.

The alternative is to represent a Complex matrix by a pair of Real
matrices. One probably needs another level of decoding expressions but I
think it might still be simpler than the first approach. There is going
to be a problem with accessing elements and the final package may be a
little less efficient that one using the previous representation.
Neverthless, this is the approach I favour.

Complex matrices are used extensively by electrical engineers and really
should be fully supported in a comprehensive package.