From pa.dec.com!decwrl!uunet!sparky!kent Sun Aug 11 15:26:26 PDT 1991 Article: 2559 of comp.sources.misc Newsgroups: comp.sources.misc Path: pa.dec.com!decwrl!uunet!sparky!kent From: Robert Davies Subject: v21i049: newmat02 - A matrix package in C++, Part00/03 Message-ID: Sender: kent@sparky.IMD.Sterling.COM (Kent Landfield) Organization: Sterling Software, IMD Date: Fri, 2 Aug 1991 01:28:32 GMT Approved: kent@sparky.imd.sterling.com Lines: 238 X-Md4-Signature: 33d0792923f597317ee2364260446a92 Submitted-by: Robert Davies 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 robert@am.dsir.govt.nz 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 (robert@am.dsir.govt.nz) exit 0 # Just in case...