Usenet 1994 October/usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso/misc/volume21/newmat02/part00

Text File | 1991-08-01 | 11KB | 247 lines

Newsgroups: comp.sources.misc
From: Robert Davies <[e-mail]>
Subject:  v21i049:  newmat02 - A matrix package in C++, Part00/03
Message-ID: <csm-v21i049=[e-mail]>
Date: Fri, 2 Aug 1991 01:28:32 GMT
Approved: [e-mail]
X-Md4-Signature: 33d0792923f597317ee2364260446a92

Submitted-by: Robert Davies <[e-mail]>
Posting-number: Volume 21, Issue 49
Archive-name: newmat02/part00
Environment: C++

MATRIX PACKAGE                                          30 July, 1991

Introduction:

  This is a description of an experimental matrix package. The package is
  an upgrade of one I distributed to a number of people in 1990.

  The package supports matrix types

  Matrix                       (rectangular matrix)
  UpperTriangularMatrix
  LowerTriangularMatrix
  DiagonalMatrix
  SymmetricMatrix
  RowVector                    (derived from Matrix)
  ColumnVector                 (derived from Matrix).

  Only one element type is supported.

  The package includes the operations *, +, -, inverse, transpose,
  conversion between types, submatrix, determinant, Cholesky decomposition
  and Householder triangularisation.

  This description begins with some general notes on matrix packages, then
  a general description of my package and a list of features and problems.

What we want from a matrix package:

  A matrix package needs to provide

  1.   Evaluation of matrix expressions in a form familiar to
       scientists and engineers. For example  A = B * (C + D.t());
  2.   Access to the elements of a matrix;
  3.   Access to submatrices;
  4.   Common elementary matrix functions such as determinant and trace;
  5.   A platform for developing advanced matrix functions such as eigen-
       value analysis;
  6.   Good efficiency and storage management.
  7.   Graceful exit from errors.

  It may also provide

  8.   A variety of types of elements (eg real and complex);
  9.   A variety of types of matrices (eg rectangular and symmetric);

  It may target small matrices (say 3 x 3), or medium sized matrices, or
  very large matrices.

  I do not believe one can build a matrix package that will meet
  everyone's requirements. In many cases if you put in one facility, you
  impose overheads on everyone using the package. This both is storage
  required for the program and in efficiency. Likewise a package that is
  optimised towards handling large matrices is likely to become less
  efficient for very small matrices where the administration time for the
  matrix may become significant compared with the time to carry out the
  operations.

  It is better to provide a variety of packages (hopefully compatible) so
  that most users can find one that meets their requirements. My package
  is intended to be one of these packages; but not all of them.

  I don't think it is obvious what is the best way of going about
  structuring a matrix package. And I don't think you can figure this
  out with "thought" experiments. Different people have to try out
  different approaches. And someone else may have to figure out which is
  best. Or, more likely, the ultimate packages will lift some ideas from
  each of a variety of trial packages. So, I don't claim my package is an
  "ultimate" package, but simply a trial of a number of ideas.

  But I do hope it will meet the immediate requirements of some people who
  need to carry out matrix calculations.

The present package:

  The package described here is intended to meet the requirements of
  someone carrying out statistical analysis and related calculations.

  Specifically it caters for rectangular matrices, upper and lower
  triangular matrices, diagonal matrices, symmetric matrices, and row and
  column vectors. Only one element type (either float or double) is
  allowed. A later version may include band matrices and/or a limited
  sparse matrix capability.

  It is intended for matrices in the range 4 x 4 to about 80 x 80. The
  upper limit is imposed by the maximum number of elements that can be
  contained in a single array (8000 doubles in some machines).

  Operations provided include +, -, *, inverse, transpose, scaling by a
  scalar, conversion between types, submatrix, determinant, trace, sum of
  squares of elements, Householder triangularisation, Cholesky
  decomposition. Singular value decomposition is to come but is not in the
  present package.

  The package is designed for version 2.0 of C++. It works with Sun C++,
  Turbo C++, Borland C++, Glockenspiel C++ (2.00a) on a PC. There are
  problems with Zortech C++ (version 2.12) and JPI C++.

Features:

  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

                      *
                     / \
                    +   C
                   / \
                  A   B

  The expression is evaluated recursively during the second stage. The two
  stage approach means that an expression can be scanned before evaluation
  and the best method of evaluation found.

  The package recognises combinations of the form  A.i()*B  where i() is
  the inverse operation and avoids evaluating the inverse explicitly.

  The package recognises and takes advantage of a matrix appearing on both
  sides of an '=' as in  A=B-A; resulting in faster execution and less
  temporary storage.  There is no need to provide an operator += since
  A=A+B  is automatically recognised as an "add into" operation.

  Matrices generated temporarily when an expression is evaluated are
  destroyed or recycled as soon as their values are no longer required.
  Each matrix has a tag which indicates if it is temporary so that its
  memory is available for recycling or deleting.

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

  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.

  In this package each matrix type includes routines for extracting
  individual rows or columns. 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. The alternative approach of using iterators
  will be slower since the iterators will involve virtual functions.

  In fact, I provide two 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.

  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.

Problems and limitations:

  The package does not have graceful exit from errors. There seems little
  point in developing this until exceptions become available in C++.

  The package does not have delayed copy. In most situations this is not a
  significant disadvantage. It is a problem with functions returning a
  matrix; see the detailed documentation.

  The package does not readily handle matrices declared constant. Every
  operation must have the option of recycling the memory of matrices in
  its arguments. So it is not possible to declare these arguments as
  constant. It is too complicated to declare versions of each operation
  for constant and non-constant arguments and there doesn't seem to be a
  way of doing automatic conversions. In the present package one needs to
  carry out an explicit conversion.

  The package is long with a large .h file. I do not recommended it if you
  are short of storage and require only matrices and vectors and do not
  require triangular or symmetric matrices. If you need to deal with
  several types of matrices then I am pretty sure you need the row and
  column operations used in this package. I am not so sure about the two
  stage approach to evaluating matrix expressions.

  It is also not recommended when you are dealing with very small matrices
  when it becomes inefficient or very large matrices when memory
  management becomes a problem (I may introduce a HugeMatrix type or alter
  the memory storage method in a later version).

  There is a problem with overloading a function to access a variety of
  matrix types. This works fine if you access only matrices with the type
  known to the compiler. But it does not work with matrix expressions
  since the compiler cannot determine the type.

  Matrix types of expressions are determined at run-time so the compiler
  will not detect errors of the type Matrix M; DiagonalMatrix D; .... D=M;
  Such errors will be detected when the statement is executed.

  Symmetric matrices are not handled very efficiently (yet) since complete
  rows are not stored explicitly.

  It will be non-trivial to include band matrices or sparse matrices in
  this package.

How to get a copy of this package:

  I am putting copies on at least some of Compuserve (Borland library, zip
  format), SIMTEL20 (zip format), comp.sources.misc on Internet (shar
  format), and on one of the program libraries at Victoria University,
  Wellington.

  Users must understand that this is a test version; there may still be
  bugs and errors. Use at your own risk. Neither I nor DSIR take any
  responsibility for any errors or omissions in this package or for any
  misfortune that may befall you or others as a result of its use.

  Please report bugs to me at

      [e-mail]

  or

      Compuserve 72777,656

  The matrix inverse routine is adapted from "Numerical Recipies in C" by
  Press, Flannery, Teukolsky, Vetterling, published by the Cambridge
  University Press.

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


Robert Davies  ([e-mail])
exit 0 # Just in case...