You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
3560 lines
174 KiB
HTML
3560 lines
174 KiB
HTML
<!DOCTYPE HTML PUBLIC "-//SoftQuad Software//DTD HoTMetaL PRO 5.0::19981217::extensions to HTML 4.0//EN" "hmpro5.dtd">
|
|
|
|
<HTML>
|
|
|
|
<HEAD>
|
|
<meta http-equiv="Content-Language" content="en-nz">
|
|
<meta name="GENERATOR" content="Microsoft FrontPage 6.0">
|
|
<meta name="ProgId" content="FrontPage.Editor.Document">
|
|
|
|
<TITLE>Newmat10 documentation</TITLE>
|
|
<link rel="stylesheet" type="text/css" href="rbd.css">
|
|
<link REL="SHORTCUT ICON" HREF="favicon.ico">
|
|
</HEAD>
|
|
|
|
<BODY>
|
|
<H1><A NAME="top"></A>Documentation for newmat10D, a matrix library in C++</H1>
|
|
<P CLASS="small"><A HREF="#intro">next</A> - <A HREF="#top">skip</A> -
|
|
<A HREF="ol_doc.htm">up</A> - <A HREF="#top">start</A><BR>
|
|
<A HREF="ol_doc.htm">return to online documentation page</A></P>
|
|
<P><B>Copyright (C) 2006: R B Davies</B></P>
|
|
<P><i>2 April, 2006</i></P>
|
|
<TABLE WIDTH=100%>
|
|
<TR>
|
|
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A HREF="#intro">1.
|
|
Introduction</A><BR>
|
|
<A HREF="#starting">2. Getting started</A><BR>
|
|
<A HREF="#refer">3. Reference manual</A></TD>
|
|
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A HREF="#error">4.
|
|
Error messages</A><BR>
|
|
<A HREF="#design">5. Design of the library</A></TD>
|
|
</TR>
|
|
</TABLE>
|
|
<P>This is the <I>how to use</I> documentation for <I>newmat</I> plus some
|
|
background information on its design. </P>
|
|
<P>There is additional support material on my <A HREF="#where">web site</A>.
|
|
</P>
|
|
<P CLASS="small">Navigation: This page is arranged in sections,
|
|
sub-sections and sub-sub-sections; four cross-references are given at the top
|
|
of these. <I>Next</I> takes you through the sections, sub-sections and
|
|
sub-sub-sections in order. <I>Skip</I> goes to the next section, sub-section or
|
|
sub-sub-section at the same level in the hierarchy as the section, sub-section
|
|
or sub-sub-section that you are currently reading. <I>Up</I> takes you up one
|
|
level in the hierarchy and <I>start</I> gets you back here.</P>
|
|
<P><B>Please read the sections on <A HREF="#custom">customising</A> and
|
|
<a href="#compiler">compilers</a> before
|
|
attempting to compile <i>newmat</i>.</B> </P>
|
|
<H2><A NAME="intro"></A>1. Introduction</H2>
|
|
<P CLASS="small"><A HREF="#use">next</A> - <A HREF="#starting">skip</A> -
|
|
<A HREF="#top">up</A> - <A HREF="#top">start</A></P>
|
|
<TABLE WIDTH=100%>
|
|
<TR>
|
|
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A HREF="#use">1.1
|
|
Conditions of use</A><BR>
|
|
<A HREF="#descript">1.2 Description</A><BR>
|
|
<A HREF="#which">1.3 Is this the library for you?</A><BR>
|
|
<A HREF="#other">1.4 Other matrix libraries</A></TD>
|
|
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A HREF="#where">1.5
|
|
Where to find this library</A><BR>
|
|
<A HREF="#author">1.6 How to contact the author</A><BR>
|
|
<A HREF="#changes">1.7 Change history</A><BR>
|
|
<A HREF="#sources">1.8 References</A></TD>
|
|
</TR>
|
|
</TABLE>
|
|
<H2><A NAME="use"></A>1.1 Conditions of use</H2>
|
|
<P CLASS="small"><A HREF="#descript">next</A> - <A HREF="#descript">skip</A> -
|
|
<A HREF="#intro">up</A> - <A HREF="#top">start</A></P>
|
|
<hr>
|
|
<p>I place no restrictions on the use of <i>newmat</i> except that I take no liability for any problems that may arise from its use,
|
|
distribution or other dealings with it.<p>
|
|
<p>
|
|
You can use it in your commercial projects.<p>
|
|
<p>
|
|
You can make and distribute modified or merged versions. You can
|
|
include parts of it in your own software.<p>
|
|
<p>
|
|
If you distribute modified or merged versions, please make it clear
|
|
which parts are mine and which parts are modified.<p>
|
|
<p>
|
|
For a substantially modified version, simply note that it is, in
|
|
part, derived from my software. A comment in the code will be
|
|
sufficient.<p>
|
|
<p>
|
|
The software is provided "as is", without warranty of any kind.<p>
|
|
<p>
|
|
Please understand that there may still be bugs and errors. Use at
|
|
your own risk. I (Robert Davies) take no responsibility for any errors
|
|
or omissions in this package or for any misfortune that may befall you
|
|
or others as a result of your use, distribution or other dealings with it.<HR>
|
|
<P>Please report bugs to me at <B>robert (at) statsresearch.co.nz</B> </P>
|
|
<P>When reporting a bug please tell me which C++ compiler you are using, and
|
|
what version. Also give me details of your computer. And tell me which version
|
|
of <I>newmat</I> (e.g. newmat03 or newmat04) you are using. Note any changes
|
|
you have made to my code. If at all possible give me a piece of code
|
|
illustrating the bug. See the <A HREF="#problem">problem report form</A>. </P>
|
|
<P><I>Please do report bugs to me.</I> </P>
|
|
<H2><A NAME="descript"></A>1.2 General description
|
|
</H2>
|
|
<P CLASS="small"><A HREF="#which">next</A> - <A HREF="#which">skip</A> -
|
|
<A HREF="#intro">up</A> - <A HREF="#top">start</A></P>
|
|
<P>The package is intended for scientists and engineers who need to manipulate
|
|
a variety of types of matrices using standard matrix operations. Emphasis is on
|
|
the kind of operations needed in statistical calculations such as least
|
|
squares, linear equation solve and eigenvalues. </P>
|
|
<P>It supports matrix types </P>
|
|
<BLOCKQUOTE>
|
|
<TABLE WIDTH="80%" cellspacing="0" style="border-collapse: collapse" bordercolor="#111111" cellpadding="0">
|
|
<TR>
|
|
<TD>Matrix</TD>
|
|
<TD>rectangular matrix</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>nricMatrix</TD>
|
|
<TD>for use with <I>Numerical Recipes in C</I> programs</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>UpperTriangularMatrix</TD>
|
|
<TD> </TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>LowerTriangularMatrix</TD>
|
|
<TD> </TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>DiagonalMatrix</TD>
|
|
<TD> </TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>SymmetricMatrix</TD>
|
|
<TD> </TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>BandMatrix</TD>
|
|
<TD> </TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>UpperBandMatrix</TD>
|
|
<TD>upper triangular band matrix</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>LowerBandMatrix</TD>
|
|
<TD>lower triangular band matrix</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>SymmetricBandMatrix</TD>
|
|
<TD> </TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>RowVector</TD>
|
|
<TD>derived from Matrix</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>ColumnVector</TD>
|
|
<TD>derived from Matrix</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>IdentityMatrix</TD>
|
|
<TD>diagonal matrix, elements have same value</TD>
|
|
</TR>
|
|
</TABLE>
|
|
</BLOCKQUOTE>
|
|
<P>Only one element type (float or double) is supported. </P>
|
|
<P>The package includes the operations <TT>*</TT>, <TT>+</TT>, <TT>-</TT>,
|
|
Kronecker product, Schur product, concatenation, inverse, transpose, conversion between types, submatrix,
|
|
determinant, Cholesky decomposition, QR triangularisation, singular value
|
|
decomposition, eigenvalues of a symmetric matrix, sorting, fast Fourier
|
|
transform, printing and an interface with <I>Numerical Recipes in C</I>. </P>
|
|
<P>It is intended for matrices in the range 10 x 10 to the maximum size your
|
|
machine will accommodate in a single array. The number of elements in an array
|
|
cannot exceed the maximum size of an <I>int</I>. The package will work for very
|
|
small matrices but becomes rather inefficient. Some of the factorisation
|
|
functions are not (yet) optimised for paged memory and so become inefficient
|
|
when used with very large matrices. </P>
|
|
<P>A <I>lazy evaluation</I> approach to evaluating matrix expressions is used
|
|
to improve efficiency and reduce the use of temporary storage. </P>
|
|
<P>I have tested versions of the package on variety of compilers and platforms
|
|
including Borland, Gnu, Microsoft, Sun and Watcom. For more details
|
|
see the section on <A HREF="#compiler">compilers</A>. </P>
|
|
<H2><A NAME="which"></A>1.3 Is this the library for
|
|
you?</H2>
|
|
<P CLASS="small"><A HREF="#other">next</A> - <A HREF="#other">skip</A> -
|
|
<A HREF="#intro">up</A> - <A HREF="#top">start</A></P>
|
|
<P>Do you</P>
|
|
<UL>
|
|
<LI>understand <TT>*</TT> to mean matrix multiply and not element by element
|
|
multiply </LI>
|
|
<LI>need matrix operators such as <TT>*</TT> and <TT>+</TT> defined as
|
|
operators so you can write things like <TT> X = A * (B + C);</TT></LI>
|
|
<LI>need a variety of types of matrices (but not sparse)</LI>
|
|
<LI>need only one element type (float or double)</LI>
|
|
<LI>work with matrices in the range 10 x 10 up to what can be stored in memory
|
|
</LI>
|
|
<LI>tolerate a moderately large but not huge package</LI>
|
|
<LI>need high quality but not necessarily the latest numerical methods. </LI>
|
|
</UL>
|
|
<P>Then <I>newmat</I> may be the right matrix library for you. </P>
|
|
<H2><A NAME="other"></A>1.4 Other matrix libraries
|
|
</H2>
|
|
<P CLASS="small"><A HREF="#where">next</A> - <A HREF="#where">skip</A> -
|
|
<A HREF="#intro">up</A> - <A HREF="#top">start</A></P>
|
|
<P>For details of other C++ matrix libraries look at
|
|
<A
|
|
HREF="http://www.robertnz.net/cpp_site.html">http://www.robertnz.net/cpp_site.html</A>.
|
|
Look at the section <I>lists of libraries</I> which gives the locations of
|
|
several very comprehensive lists of matrix and other C++ libraries and at the
|
|
section <I>source code</I>. </P>
|
|
<H2><A NAME="where"></A>1.5 Where to find this
|
|
library</H2>
|
|
<P CLASS="small"><A HREF="#author">next</A> - <A HREF="#author">skip</A> -
|
|
<A HREF="#intro">up</A> - <A HREF="#top">start</A> </P>
|
|
<UL>
|
|
<LI><a href="http://www.robertnz.net">http://www.robertnz.net</a></LI>
|
|
</UL>
|
|
<H2><A NAME="author"></A>1.6 How to contact the
|
|
author</H2>
|
|
<P CLASS="small"><A HREF="#changes">next</A> - <A HREF="#changes">skip</A> -
|
|
<A HREF="#intro">up</A> - <A HREF="#top">start</A></P>
|
|
<PRE> Robert Davies
|
|
16 Gloucester Street
|
|
Wilton
|
|
Wellington
|
|
New Zealand
|
|
|
|
<I>email:</I> robert<b> at </b>statsresearch.co.nz
|
|
</PRE>
|
|
|
|
<H2><A NAME="changes"></A>1.7 Change history</H2>
|
|
<P CLASS="small"><A HREF="#sources">next</A> - <A HREF="#sources">skip</A> -
|
|
<A HREF="#intro">up</A> - <A HREF="#top">start</A></P>
|
|
<P><b>Newmat10D - April, 2006:</b></P>
|
|
<P class="small">Compatibility fix for Gnu G++ 4.1.</P>
|
|
<P><b>Newmat10C - March, 2006:</b></P>
|
|
<P class="small">Update conditions of use to be acceptable for Debian
|
|
distribution, a few minor fixes; additional <i>make</i> files; update include.h, myexpect.h,
|
|
myexcept.cpp, precisio.h. </P>
|
|
<P><b>Newmat10B - January, 2005:</b></P>
|
|
<P class="small">Fix compatibility problems with Gnu G++ 3.4 and Intel 8.1
|
|
compilers; update include.h, myexpect.h,
|
|
myexcept.cpp, precisio.h.</P>
|
|
<P><b>Newmat10A - October, 2002:</b></P>
|
|
<P CLASS="small">Fix error in Kronecker product; fixes for Intel and GCC3
|
|
compilers.</P>
|
|
<P><B>Newmat10 - January, 2002:</B> </P>
|
|
<P CLASS="small">Improve compatibility with GCC, fix errors in FFT and
|
|
GenericMatrix, update simulated exceptions, maxima, minima, determinant, dot
|
|
product and Frobenius norm functions, update make files for CC and GCC, faster
|
|
FFT, <TT>A.ReSize(B)</TT>, fix pointer arithmetic, <TT><<</TT> for loading
|
|
data into rows, IdentityMatrix, Kronecker product, sort singular values.</P>
|
|
<P><B>Newmat09 - September, 1997:</B> </P>
|
|
<P CLASS="small">Operator <TT>==</TT>, <TT>!=</TT>, <TT>+=</TT>, <TT>-=</TT>,
|
|
<TT>*=</TT>, <TT>/=</TT>, <TT>|=</TT>, <TT>&=</TT>. Follow new rules for
|
|
<I>for (int i; ... )</I> construct. Change Boolean, TRUE, FALSE to bool, true,
|
|
false. Change ReDimension to ReSize. SubMatrix allows zero rows and columns.
|
|
Scalar <TT>+</TT>, <TT>-</TT> or <TT>*</TT> matrix is OK. Simplify simulated
|
|
exceptions. Fix non-linear programs for AT&T compilers. Dummy inequality
|
|
operators. Improve internal row/column operations. Improve matrix LU
|
|
decomposition. Improve sort. Reverse function. IsSingular function. Fast trig
|
|
transforms. Namespace definitions. </P>
|
|
<P><B>Newmat08A - July, 1995:</B> </P>
|
|
<P CLASS="small">Fix error in SVD. </P>
|
|
<P><B>Newmat08 - January, 1995:</B> </P>
|
|
<P CLASS="small">Corrections to improve compatibility with Gnu, Watcom.
|
|
Concatenation of matrices. Elementwise products. Option to use compilers
|
|
supporting exceptions. Correction to exception module to allow global
|
|
declarations of matrices. Fix problem with inverse of symmetric matrices. Fix
|
|
divide-by-zero problem in SVD. Include new QR routines. Sum function.
|
|
Non-linear optimisation. GenericMatrices. </P>
|
|
<P><B>Newmat07 - January, 1993</B> </P>
|
|
<P CLASS="small">Minor corrections to improve compatibility with Zortech,
|
|
Microsoft and Gnu. Correction to exception module. Additional FFT functions.
|
|
Some minor increases in efficiency. Submatrices can now be used on RHS of =.
|
|
Option for allowing C type subscripts. Method for loading short lists of
|
|
numbers. </P>
|
|
<P><B>Newmat06 - December 1992:</B> </P>
|
|
<P CLASS="small">Added band matrices; 'real' changed to 'Real' (to avoid
|
|
potential conflict in complex class); Inject doesn't check for no loss of
|
|
information; fixes for AT&T C++ version 3.0; real(A) becomes A.AsScalar();
|
|
CopyToMatrix becomes AsMatrix, etc; .c() is no longer required (to be deleted
|
|
in next version); option for version 2.1 or later. Suffix for include files
|
|
changed to .h; BOOL changed to Boolean (BOOL doesn't work in g++ v 2.0);
|
|
modifications to allow for compilers that destroy temporaries very quickly;
|
|
(Gnu users - see the section of compilers). Added CleanUp,
|
|
LinearEquationSolver, primitive version of exceptions. </P>
|
|
<P><B>Newmat05 - June 1992:</B> </P>
|
|
<P CLASS="small">For private release only </P>
|
|
<P><B>Newmat04 - December 1991:</B> </P>
|
|
<P CLASS="small">Fix problem with G++1.40, some extra documentation
|
|
</P>
|
|
<P><B>Newmat03 - November 1991:</B> </P>
|
|
<P CLASS="small">Col and Cols become Column and Columns. Added Sort, SVD,
|
|
Jacobi, Eigenvalues, FFT, real conversion of 1x1 matrix, <I>Numerical Recipes
|
|
in C</I> interface, output operations, various scalar functions. Improved
|
|
return from functions. Reorganised setting options in "include.hxx".
|
|
</P>
|
|
<P><B>Newmat02 - July 1991:</B> </P>
|
|
<P CLASS="small">Version with matrix row/column operations and numerous
|
|
additional functions. </P>
|
|
<P><B>Matrix - October 1990:</B> </P>
|
|
<P CLASS="small">Early version of package. </P>
|
|
<H2><A NAME="sources"></A>1.8 References</H2>
|
|
<P CLASS="small"><A HREF="#starting">next</A> - <A HREF="#starting">skip</A> -
|
|
<A HREF="#intro">up</A> - <A HREF="#top">start</A></P>
|
|
<UL>
|
|
<LI>The matrix LU decomposition is from Golub, G.H. & Van Loan, C.F.
|
|
(1996), <I>Matrix Computations</I>, published by Johns Hopkins University
|
|
Press. </LI>
|
|
<LI>Part of the matrix inverse/solve routine is adapted from Press, Flannery,
|
|
Teukolsky, Vetterling (1988), <I>Numerical Recipes in C</I>, published by the
|
|
Cambridge University Press. </LI>
|
|
<LI>Many of the advanced matrix routines are adapted from routines in Wilkinson
|
|
and Reinsch (1971), <I>Handbook for Automatic Computation, Vol II, Linear
|
|
Algebra</I> published by Springer Verlag. </LI>
|
|
<LI>The fast Fourier transform is adapted from Carl de Boor (1980), <I>Siam J
|
|
Sci Stat Comput</I>, pp173-8 and the fast trigonometric transforms from Charles
|
|
Van Loan (1992) in <I>Computational frameworks for the fast Fourier
|
|
transform</I> published by SIAM. </LI>
|
|
<LI>The sort function is derived from Sedgewick, Robert (1992), <I>Algorithms
|
|
in C++</I> published by Addison Wesley. </LI>
|
|
</UL>
|
|
<P>For references about <I>Newmat</I> see </P>
|
|
<UL>
|
|
<LI>Davies, R.B. (1994) Writing a matrix package in C++. In OON-SKI'94: The
|
|
second annual object-oriented numerics conference, pp 207-213. Rogue Wave
|
|
Software, Corvallis. </LI>
|
|
<LI>Eddelbuttel, Dirk (1996) Object-oriented econometrics: matrix programming
|
|
in C++ using GCC and Newmat. Journal of Applied Econometrics, Vol 11, No 2, pp
|
|
199-209. </LI>
|
|
</UL>
|
|
<H2><A NAME="starting"></A>2. Getting started</H2>
|
|
<P CLASS="small"><A HREF="#overview">next</A> - <A HREF="#refer">skip</A> -
|
|
<A HREF="#top">up</A> - <A HREF="#top">start</A></P>
|
|
<TABLE WIDTH="100%">
|
|
<TR>
|
|
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A
|
|
HREF="#overview">2.1 Overview</A><BR>
|
|
<A HREF="#make">2.2 Make files</A><BR>
|
|
<A HREF="#custom">2.3 Customising</A><BR>
|
|
<A HREF="#compiler">2.4 Compilers</A><BR>
|
|
<a HREF="#update">2.5 Updating from previous versions</a><br>
|
|
<a HREF="#except_1">2.6 Catching exceptions</a></TD>
|
|
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A
|
|
HREF="#example">2.7 Example</A><BR>
|
|
<A HREF="#testing">2.8 Testing</A><BR>
|
|
<A HREF="#bugs">2.9 Bugs</A><BR>
|
|
<A HREF="#files">2.10 Files in newmat10</A><BR>
|
|
<A HREF="#problem">2.11 Problem report form</A></TD>
|
|
</TR>
|
|
</TABLE>
|
|
<H2><A NAME="overview"></A>2.1 Overview</H2>
|
|
<P CLASS="small"><A HREF="#make">next</A> - <A HREF="#make">skip</A> -
|
|
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
|
|
<P>I use .h as the suffix of definition files and .cpp as the suffix of C++
|
|
source files. </P>
|
|
<P>You will need to compile all the *.cpp files listed as program files in the
|
|
<A HREF="#files">files section</A> to get the complete package. Ideally you
|
|
should store the resulting object files as a library. The tmt*.cpp files are
|
|
used for <A HREF="#testing">testing</A>, example.cpp is an <A
|
|
HREF="#example">example</A> and sl_ex.cpp, nl_ex.cpp and garch.cpp are examples
|
|
of the <A HREF="#nonlin">non-linear</A> solve and optimisation routines. A
|
|
demonstration and test of the exception mechanism is in test_exc.cpp. </P>
|
|
<P>I include a number of <I>make</I> files for compiling the example and the
|
|
test package. See the section on <A HREF="#make">make files</A> for details.
|
|
But with the PC compilers, its pretty quick just to load all the files in the
|
|
interactive environments by pointing and clicking. </P>
|
|
<P>Use the large or win32 console model when you are using a PC. Do not
|
|
<I>outline</I> inline functions. You may need to increase the stack size. </P>
|
|
<P>Your source files that access the newmat will need to #include one or more
|
|
of the following files. </P>
|
|
<TABLE WIDTH="100%" CELLPADDING="3" cellspacing="0" style="border-collapse: collapse" bordercolor="#111111">
|
|
<TR>
|
|
<TD VALIGN="TOP">include.h</TD>
|
|
<TD>to access just the compiler options</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD VALIGN="TOP">newmat.h</TD>
|
|
<TD>to access just the main matrix library (includes include.h)</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD VALIGN="TOP">newmatap.h</TD>
|
|
<TD>to access the advanced matrix routines such as Cholesky decomposition, QR
|
|
triangularisation etc (includes newmat.h)</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD VALIGN="TOP">newmatio.h</TD>
|
|
<TD>to access the <A HREF="#output">output</A> routines (includes newmat.h) You
|
|
can use this only with compilers that support the standard input/output
|
|
routines including manipulators</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD VALIGN="TOP">newmatnl.h</TD>
|
|
<TD>to access the non-linear optimisation routines (includes newmat.h)</TD>
|
|
</TR>
|
|
</TABLE>
|
|
<P>See the section on <A HREF="#custom">customising</A> to see how to edit
|
|
include.h for your environment and the section on <A
|
|
HREF="#compiler">compilers</A> for any special problems with the compiler you
|
|
are using. </P>
|
|
<H2><A NAME="make"></A>2.2 Make files</H2>
|
|
<P CLASS="small"><A HREF="#custom">next</A> - <A HREF="#custom">skip</A> -
|
|
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
|
|
<P>I have included <I>make</I> files for CC, Microsoft, Intel, Borland 5.5 and Gnu compilers
|
|
for compiling the examples. You
|
|
can generate make files for a number of other compilers with my
|
|
<a href="genmake.htm">genmake</a> utility. Make files provide
|
|
a way of compiling your programs without using the IDE that comes with
|
|
PC compilers. See the <A HREF="#files">files section</A> for details. See the
|
|
<A HREF="#example">example</A> for how to use them. Leave out the target name
|
|
to compile and link all my examples and test files. For more information on how
|
|
to use these files see the documentation for my
|
|
<a href="genmake.htm">genmake</a> utility.</P>
|
|
<H3>PC</H3>
|
|
<P>I include make files for Microsoft, Intel, Borland 5.5. For Borland you will need to edit it to show where
|
|
you have stored your Borland compiler. For make files for other compilers use my
|
|
<a href="genmake.htm">genmake</a> utility.</P>
|
|
<H3>Unix</H3>
|
|
<P>The <I>make</I> file for the Unix CC compilers link a .cxx file to each .cpp
|
|
file since some of these compilers do not recognise .cpp as a legitimate
|
|
extension for a C++ file. I suggest you delete this part of the <I>make</I>
|
|
file and, if necessary, rename the .cpp files to something your compiler
|
|
recognises. </P>
|
|
<P>My <I>make</I> file for Gnu GCC on Unix systems is for use with
|
|
<TT>gmake</TT> rather than <TT>make</TT>. I assume your compiler recognises the
|
|
.cpp extension. Ordinary <TT>make</TT> works with it on the Sun but not the
|
|
Silicon Graphics or HP machines. On Linux use <TT>make</TT>. </P>
|
|
<P>My make file for the CC compilers works with the ordinary make. </P>
|
|
<P>To compile everything with the CC compiler use </P>
|
|
<PRE> make -f nm_cc.mak
|
|
</PRE>
|
|
|
|
<P>or for the gnu compiler use </P>
|
|
<PRE> gmake -f nm_gnu.mak
|
|
</PRE>
|
|
|
|
<P>There is a line in the make file for CC <TT>rm -f $*.cxx</TT>. Some systems
|
|
won't accept this line and you will need to delete it. In this case, if you
|
|
have a bad compile and you are using my scheme for linking .cxx files, you will
|
|
need to delete the .cxx file link generated by that compile before you can do
|
|
the next one.</P>
|
|
<P>There is also a <i>make</i> file for the Intel compiler for Linux. </P>
|
|
<H2><A NAME="custom"></A>2.3 Customising</H2>
|
|
<P CLASS="small"><A HREF="#compiler">next</A> - <A HREF="#compiler">skip</A> -
|
|
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
|
|
<P>The file <I>include.h</I> sets a variety of options including several
|
|
compiler dependent options. You may need to edit include.h to get the options
|
|
you require. If you are using a compiler different from one I have worked with
|
|
you may have to set up a new section in include.h appropriate for your
|
|
compiler. </P>
|
|
<P>Borland, Turbo, Gnu, Microsoft and Watcom are recognised automatically. If
|
|
none of these are recognised a default set of options is used. These are fine
|
|
for AT&T, HPUX and Sun C++. If you using a compiler I don't know about you
|
|
may have to write a new set of options. </P>
|
|
<P>There is an option in include.h for selecting whether you use compiler
|
|
supported exceptions, simulated exceptions, or disable exceptions. I now set
|
|
<I> compiler supported exceptions</I> as the default. Use the option for
|
|
compiler supported exceptions <I>if and only if</I> you have set the option on
|
|
your compiler to recognise exceptions. Disabling exceptions sometimes helps
|
|
with compilers that are incompatible with my exception simulation scheme. </P>
|
|
<P>If you are using an older compiler that does <I>not</I> recognises<B>
|
|
bool</B> as required by the standard then <I>de-activate</I> the statement
|
|
<TT>#define bool_LIB</TT>. This will turn <I>on</I> my Boolean class. </P>
|
|
<P>Activate the appropriate statement to make the element type float or double.
|
|
</P>
|
|
<P>I suggest you leave the options TEMPS_DESTROYED_QUICKLY,
|
|
TEMPS_DESTROYED_QUICKLY_R <I>de-activated</I>, unless you are using a very old
|
|
version of Gnu compiler (<2.6). This stores the <I>trees</I> describing
|
|
matrix expressions on the stack rather than the heap. See the discussion on
|
|
<A HREF="#destr">destruction of temporaries</A> for more explanation.</P>
|
|
<P>The option <A HREF="#testing">DO_FREE_CHECK</A> is used for tracking memory
|
|
leaks and normally should not be activated. </P>
|
|
<P>Activate SETUP_C_SUBSCRIPTS if you want to use traditional C style
|
|
<A HREF="#elements">element access</A>. Note that this does <I>not</I> change
|
|
the starting point for indices when you are using round brackets for accessing
|
|
elements or selecting submatrices. It does enable you to use C style square
|
|
brackets.</P>
|
|
<P>Activate <TT>#define use_namespace</TT> if you want to use <A
|
|
HREF="#namesp">namespaces</A>. Do this only if you are sure your compiler
|
|
supports namespaces. If you do turn this option on, be prepared to turn it off
|
|
again if the compiler reports inaccessible variables or the linker reports
|
|
missing links. </P>
|
|
<P>Activate <TT>#define _STANDARD_</TT> to use the standard names for the
|
|
included files and to find the floating point precision data using the floating
|
|
point standard. This will work only with the most recent compilers. This is
|
|
automatically turned on for the Gnu compiler version 3 and the Intel compiler
|
|
for Linux. </P>
|
|
<P>If you haven't defined <TT>_STANDARD_</TT> and are using a compiler that
|
|
<I>include.h</I> does not recognise and you want to pick up the floating point
|
|
precision data from <I>float.h</I> then activate <TT>#define use_float_h</TT>.
|
|
Otherwise the floating point precision data will be accessed from
|
|
<I>values.h</I>. You may need to do this with computers from Digital, in
|
|
particular. </P>
|
|
<H2><A NAME="compiler"></A>2.4 Compilers</H2>
|
|
<P CLASS="small"><A HREF="#atandt">next</A> - <A HREF="#update">skip</A> -
|
|
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
|
|
<TABLE WIDTH="100%">
|
|
<TR>
|
|
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%">
|
|
<A HREF="#atandt">2.4.1 AT&T </A><BR>
|
|
<A HREF="#borland">2.4.2 Borland </A><BR>
|
|
<A HREF="#gcc">2.4.3 Gnu G++ </A><BR>
|
|
<A HREF="#hpux">2.4.4 HPUX </A> </TD>
|
|
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%">
|
|
<A HREF="#intel">2.4.5 Intel </A><BR>
|
|
<A HREF="#microso">2.4.6 Microsoft </A><BR>
|
|
<A HREF="#sun">2.4.7 Sun </A><BR>
|
|
<A HREF="#watcom">2.4.8 Watcom </A></TD>
|
|
</TR>
|
|
</TABLE>
|
|
<P>I have tested this library on a number of compilers. Here are the levels of
|
|
success and any special considerations. In most cases I have chosen code that
|
|
works under all the compilers I have access to, but I have had to include some
|
|
specific work-arounds for some compilers. For the PC versions, I use Pentium 3
|
|
& 4 computers running windows 2000 or XP or various varieties of Linux
|
|
(Red Hat or Fedora). The Unix versions are on a Sun Sparc
|
|
station. Thanks to Victoria University for access to the Sparc. </P>
|
|
<P>I have set up a block of code for each of the compilers in include.h. Turbo,
|
|
Borland, Gnu, Microsoft and Watcom are recognised automatically. There is a
|
|
default option that works for AT&T, Sun C++ and HPUX. So you don't
|
|
have to make any changes for these compilers. Otherwise you may have to build
|
|
your own set of options in include.h. </P>
|
|
<H2><A NAME="atandt"></A>2.4.1 AT&T</H2>
|
|
<P CLASS="small"><A HREF="#borland">next</A> - <A HREF="#borland">skip</A> -
|
|
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>
|
|
<P>The AT&T compiler used to be available on a wide variety of Unix
|
|
workstations. I don't know if anyone still uses it. However the AT&T options are
|
|
the default if your compiler is not recognised.</P>
|
|
<P>AT&T C++ 2.1; 3.0.1 on a Sun: Previous versions worked on these
|
|
compilers, which I no longer have access to. </P>
|
|
<P>In AT&T 2.1 you may get an error when you use an expression for the
|
|
single argument when constructing a Vector or DiagonalMatrix or one of the
|
|
Triangular Matrices. You need to evaluate the expression separately. </P>
|
|
<H2><A NAME="borland"></A>2.4.2 Borland</H2>
|
|
<P CLASS="small"><A HREF="#gcc">next</A> - <A HREF="#gcc">skip</A> -
|
|
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>
|
|
<H3>Newer compilers</H3>
|
|
<P><b>Borland Builder version 6:</b> This is <b>not</b> compatible with <i>
|
|
newmat10</i>. Use <i>newmat11</i> instead.</P>
|
|
<P><B>Borland Builder version 5:</B> This works fine in console mode and no
|
|
special editing of the source codes is required. I haven't tested it in GUI
|
|
mode. You can set the <i>newmat10</i> options to use namespace and the standard
|
|
library. <b> You
|
|
should turn <I>off</I> the Borland option to use pre-compiled headers.</b> There
|
|
are notes on compiling with the IDE on my <a href="#where">website</a>.
|
|
Alternatively you can use the <i>nm_b55.mak</i> make file<i>.</i></P>
|
|
<P><B>Borland Builder version 4</B>: I have successfully used this on older
|
|
versions of newmat using the
|
|
console wizard (menu item file/new - select new tab). Use compiler
|
|
exceptions. Suppose you are compiling my test program <I>tmt</I>. Rename my
|
|
<I>main()</I> function in <I>tmt.cpp</I> to <I>my_main()</I>. Rename
|
|
<I>tmt.cpp</I> to <I>tmt_main.cpp</I>. Borland will generate a new file
|
|
<I>tmt.cpp</I> containing their <I>main()</I> function. Put the line <TT>int
|
|
my_main();</TT> above this function and put <TT>return my_main();</TT> into the
|
|
body of <I>main()</I>.</P>
|
|
<P><B>Borland compiler version 5.5</B>: this is the free C++ compiler available
|
|
from Borland's web site. I suggest you use
|
|
the compiler supported exceptions and turn on <I>standard</I> in include.h. You
|
|
can use the make file <i>nm_b55.mak</i> after editing to correct the file locations for
|
|
your system.</P>
|
|
<H3>Older compilers</H3>
|
|
<P><B>Borland C++ 3.1, 5.02</B>: Use the simulated exceptions with these. Then
|
|
version 5.02 works OK. You will need to use
|
|
the large or 32 bit flat model. If you are not debugging, turn off the options
|
|
that collect debugging information. It compiles with version 3.1 but you can't
|
|
run the tmt test program.</P>
|
|
<P>If you are using versions earlier than 5 remember to edit include.h to
|
|
activate my Boolean class. </P>
|
|
<P>When running my test program under ms-dos you may run out of memory. Either
|
|
compile the test routine to run under <I>easywin</I> or use simulated exceptions
|
|
rather than the built in exceptions. </P>
|
|
<P>If you can, upgrade to windows 9X or window NT and use the 32 bit console
|
|
model. </P>
|
|
<P>If you are using the 16 bit large model, don't forget to keep all matrices
|
|
less than 64K bytes in length (90x90 for a rectangular matrix if you are using
|
|
<TT>double</TT> as your element type). Otherwise your program will crash
|
|
without warning or explanation. You will need to break the <A
|
|
HREF="#testing">tmt</A> set of test files into several parts to get the program
|
|
to fit into your computer and run without stack overflow. </P>
|
|
<P>One version of Borland had DBL_MIN incorrectly defined. If you are using an
|
|
older version of Borland and are getting strange numerical errors in the test
|
|
programs reinstate the commented out statements in precision.h. </P>
|
|
<P>You can generate make files for versions 5 or 3.1 with my <a href="genmake.htm">
|
|
genmake</a> utility.</P>
|
|
<H2><A NAME="gcc"></A>2.4.3 Gnu G++</H2>
|
|
<P CLASS="small"><A HREF="#hpux">next</A> - <A HREF="#hpux">skip</A> -
|
|
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>
|
|
<P><B>Gnu G++ 3.3, 4.0, 4.1:</B> These work OK. If you are using a much
|
|
earlier version see if you can upgrade. Standard is automatically turned on with 3.X.</P>
|
|
<P>If you are using 2.6 or earlier remember to edit include.h to activate my
|
|
Boolean class. In 2.6.?, <TT>fabs(*X++)</TT> causes a problem. You may need to write you
|
|
own non-inlined version. </P>
|
|
<P>For versions earlier than 2.6.0 you must enable the options
|
|
TEMPS_DESTROYED_QUICKLY and TEMPS_DESTROYED_QUICKLY_R. You can't use
|
|
expressions like <TT>Matrix(X*Y)</TT> in the middle of an expression and
|
|
<TT>(Matrix)(X*Y)</TT> is unreliable. If you write a function returning a
|
|
matrix, you MUST use the <A HREF="#memory">ReturnMatrix</A> method described in
|
|
this documentation. This is because g++ destroys temporaries occurring in an
|
|
expression too soon for the two stage way of evaluating expressions that newmat
|
|
uses. You will have problems with versions of Gnu earlier than 2.3.1. </P>
|
|
<H2><A NAME="hpux"></A>2.4.4 HP-UX</H2>
|
|
<P CLASS="small"><A HREF="#intel">next</A> - <A HREF="#intel">skip</A> -
|
|
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>
|
|
<P>HP 9000 series HP-UX. I no longer have access to this compiler. Newmat09
|
|
worked without problems with the simulated exceptions; haven't tried the
|
|
built-in exceptions. </P>
|
|
<P>With recent versions of the compiler you may get warning messages like
|
|
<TT>Unsafe cast between pointers/references to incomplete classes</TT>. At
|
|
present, I think these can be ignored. </P>
|
|
<P>Here are comments I made in 1997. </P>
|
|
<P>I have tried the library on two versions of HP-UX. (I don't know the version
|
|
numbers, the older is a clone of AT&T 3, the newer is HP's version with
|
|
exceptions). Both worked after the modifications described in this section.
|
|
</P>
|
|
<P>With the older version of the compiler I needed to edit the math.h library
|
|
file to remove a duplicate definition of abs. </P>
|
|
<P>With the newer version you can set the +eh option to enable exceptions and
|
|
activate the UseExceptions option in include.h. If you are using my make file,
|
|
you will need to replace CC with CC +eh where ever CC occurs. I recommend that
|
|
you do not do this and either disable exceptions or use my simulated
|
|
exceptions. I get core dumps when I use the built-in exceptions and suspect
|
|
they are not sufficiently debugged as yet. </P>
|
|
<P>If you are using my simulated exceptions you may get a mass of error
|
|
messages from the linker about __EH_JMPBUF_TEMP. In this case get file setjmp.h
|
|
(in directory /usr/include/CC ?) and put extern in front of the line </P>
|
|
<PRE> jmp_buf * __EH_JMPBUF_TEMP;
|
|
</PRE>
|
|
|
|
<P>The file setjmp.h is accessed in my file myexcept.h. You may want to change
|
|
the #include statement to access your edited copy of setjmp.h. </P>
|
|
|
|
<h2><a name="intel"></a>2.4.5 Intel</h2>
|
|
<P CLASS="small"><A HREF="#microso">next</A> - <A HREF="#microso">skip</A> -
|
|
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>
|
|
|
|
<P>Newmat works correctly with the Intel 9 C++ compilers for Windows and for
|
|
Linux. Standard is automatically switched on for the Linux version and for the
|
|
Windows version if you are emulating VC++ 7 or higher.</P>
|
|
|
|
<H2><A NAME="microso"></A>2.4.6 Microsoft</H2>
|
|
<P CLASS="small"><A HREF="#sun">next</A> - <A HREF="#sun">skip</A> -
|
|
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>
|
|
<H3>Newer versions</H3>
|
|
<P>See my <A HREF="#where">web site</A> for instructions how to work
|
|
Microsoft's IDE.</P>
|
|
<P><b>Microsoft Visual C++ 8:</b> I have tested the <i>express</i> version using
|
|
my make file run the compiler.</P>
|
|
<P><b>Microsoft Visual C++ 7, 7.1:</b> This works OK. Note that all my tests have
|
|
been in console mode. The standard option is on by default but I am still a bit
|
|
wary about the namespace option.</P>
|
|
<P><B>Microsoft Visual C++ 6:</B> <b>Get the latest service pack</b>. I have tried this
|
|
in console mode and it seems to work satisfactorily. Use the compiler supported exceptions. You
|
|
<i>may</i> be able to
|
|
use the namespace and standard options but I suggest not using namespace. If you want to work under MFC
|
|
you may need to <TT>#include "stdafx.h"</TT> at the beginning of each .cpp file. </P>
|
|
<P><B>Microsoft Visual C++ 5</B>: I have tried this in console mode and it
|
|
seems to work satisfactorily. There may be a problem with <A
|
|
HREF="#namesp">namespace</A> (fixed by Service Pack 3?). <B>Turn optimisation
|
|
off</B>. Use the compiler supported exceptions. If
|
|
you want to work under MFC
|
|
you may need to <TT>#include "stdafx.h"</TT> at the
|
|
beginning of each .cpp file.</P>
|
|
<H3>Older versions</H3>
|
|
<P>Microsoft Visual C++ 2.0: This used to work OK. I haven't tried it with
|
|
recent versions of newmat.</P>
|
|
<P>You must <TT>#define TEMPS_DESTROYED_QUICKLY</TT> owing to a bug in version
|
|
7 (at least) of MSC. There are some notes in the file <TT>include.h</TT> on
|
|
changes to run under version 7. I haven't tried newmat10 on version 7. </P>
|
|
<P>Microsoft Visual C++ 1.51. Disable exceptions, comment out the line in
|
|
include.h <TT>#define TEMPS_DESTROYED_QUICKLY_R</TT>. In <TT>tmt.cpp</TT>,
|
|
comment out the <TT>Try</TT> and <TT>CatchAll</TT> lines at the beginning of
|
|
<TT>main()</TT> and the line <TT>trymati()</TT>. You can use the makefile
|
|
<TT>ms.mak</TT>. You will probably need to break the <A HREF="#testing">tmt</A>
|
|
test files into two parts to get the program to link. </P>
|
|
<P>If you can, upgrade to windows 95, 98 or window NT and use the 32 bit
|
|
console model. </P>
|
|
<P>If you are using the 16 bit large model, don't forget to keep all matrices
|
|
less than 64K bytes in length (90x90 for a rectangular matrix if you are using
|
|
<TT>double</TT> as your element type). Otherwise your program will crash
|
|
without warning or explanation. You may need to break the <A
|
|
HREF="#testing">tmt</A> set of test files into two parts to get the program to
|
|
fit into your computer. </P>
|
|
<P>Microsoft Visual C++ 4: I haven't tried this - a correspondent reports: I
|
|
use Microsoft Visual C++ Version 4. there is only one minor problem. In all
|
|
files you must include <TT>#include "stdafx.h"</TT> (presumably if
|
|
you are using MFC). This file contains essential information for VC++. Leave it
|
|
out and you get <I>Unexpected end of file</I>. </P>
|
|
<H2><A NAME="sun"></A>2.4.7 Sun</H2>
|
|
<P CLASS="small"><A HREF="#watcom">next</A> - <A HREF="#watcom">skip</A> -
|
|
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>
|
|
<P><B>Sun C++:</B> The current version works fine with
|
|
compiler supported exceptions. <B>Sun C++ (version
|
|
5):</B> There seems to be a problem with exceptions. If you use my simulated
|
|
exceptions the non-linear optimisation programs hang. If you use the compiler
|
|
supported exceptions my tmt and test_exc programs crash. You should
|
|
<I>disable</I> exceptions.</P>
|
|
<H2><A NAME="watcom"></A>2.4.8 Watcom</H2>
|
|
<P CLASS="small"><A HREF="#update">next</A> - <A HREF="#update">skip</A> -
|
|
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>
|
|
<P><B>Open Watcom C++:</B> this works fine.</P>
|
|
<H2><A NAME="update"></A>2.5 Updating from previous
|
|
versions</H2>
|
|
<P CLASS="small"><A HREF="#example">next</A> - <A HREF="#example">skip</A> -
|
|
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
|
|
<P><B>Newmat10</B> includes new <A HREF="#scalar2">maxima, minima</A>,
|
|
<A HREF="#scalar3">determinant, dot product and Frobenius norm</A> functions, a
|
|
faster <A HREF="#fft">FFT</A>, revised <A HREF="#make">make</A> files for GCC
|
|
and CC compilers, several corrections, new <A HREF="#dimen">ReSize</A> function, <a href="#constr">IdentityMatrix</a>
|
|
and <a href="#binary">Kronecker Product</a>. Singular values from <a href="#svd">
|
|
SVD</a> are sorted. The program files include a new file, <TT>newfft.cpp</TT>, so you
|
|
will need to include this in the list of files in your IDE and make files. There
|
|
is also a new test file tmtm.cpp.
|
|
<A HREF="#pointer">Pointer arithmetic</A> now mostly meets requirements of
|
|
standard. You can use <A HREF="#entering"><<</A> to load data into rows
|
|
of a matrix. The <A HREF="#custom">default options</A> in include.h have been
|
|
changed. If you are updating from a beta version of newmat09 look through the
|
|
next section as there were some late changes to newmat09. </P>
|
|
<P>If you are upgrading from <B>newmat08</B> note the following:
|
|
</P>
|
|
<UL>
|
|
<LI>Boolean, TRUE, FALSE are now bool, true, false. See
|
|
<A HREF="#custom">customising</A> if your compiler supports the bool class.
|
|
</LI>
|
|
<LI>ReDimension is now <A HREF="#dimen">ReSize</A>.
|
|
</LI>
|
|
<LI>The <A HREF="#except">simulated exception</A> package has
|
|
been updated. </LI>
|
|
<LI>Operators <TT>==</TT>, <TT>!=</TT>, <TT>+=</TT>,
|
|
<TT>-=</TT>, <TT>*=</TT>, <TT>|=</TT>, <TT>&=</TT> are now supported as
|
|
<A HREF="#binary">binary</A> matrix operators. </LI>
|
|
<LI><TT>A+=f</TT>, <TT>A-=f</TT>, <TT>A*=f</TT>, <TT>A/=f</TT>,
|
|
<TT>f+A</TT>, <TT>f-A</TT>, <TT>f*A</TT> are supported for <A HREF="#matscal">A
|
|
matrix, f scalar</A>. </LI>
|
|
<LI><A HREF="#trigtran">Fast trigonometric transforms</A>.
|
|
</LI>
|
|
<LI><A HREF="#unary">Reverse</A> function for reversing order of
|
|
elements in a vector or matrix. </LI>
|
|
<LI><A HREF="#scalar3">IsSingular</A> function. </LI>
|
|
<LI>An option is included for defining <A
|
|
HREF="#namesp">namespaces</A>. </LI>
|
|
<LI>Dummy inequality operators are defined for compatibility
|
|
with the STL. </LI>
|
|
<LI>The row/column classes in newmat3.cpp have been modified to
|
|
improve efficiency and correct an invalid use of pointer arithmetic. Most users
|
|
won't be using these classes explicitly; if you are, please contact me for
|
|
details of the changes. </LI>
|
|
<LI>Matrix LU decomposition rewritten (faster for large arrays).
|
|
</LI>
|
|
<LI>The sort function rewritten (faster). </LI>
|
|
<LI>The documentation files newmata.txt and newmatb.txt have
|
|
been amalgamated and both are included in the hypertext version. </LI>
|
|
<LI>Some of the <A HREF="#make">make</A> files reorganised
|
|
again. </LI>
|
|
</UL>
|
|
<P>If you are upgrading from <B>newmat07</B> note the following:
|
|
</P>
|
|
<UL>
|
|
<LI>.cxx files are now .cpp files. Some versions of won't accept
|
|
.cpp. The <I>make</I> files for Gnu and AT&T link the .cpp files to .cxx
|
|
files before compilation and delete the links after compilation. </LI>
|
|
<LI>An <A HREF="#custom">option</A> in include.h allows you to
|
|
use compiler supported exceptions, simulated exceptions or disable exceptions.
|
|
Edit the file include.h to select one of these three options. Don't simulate
|
|
exceptions if you have set your compiler's option to implement exceptions.
|
|
</LI>
|
|
<LI>New <A HREF="#qr">QR decomposition</A> functions.
|
|
</LI>
|
|
<LI>A <A HREF="#nonlin">non-linear least squares</A> class.
|
|
</LI>
|
|
<LI>No need to explicitly set the AT&T option in include.h.
|
|
</LI>
|
|
<LI><A HREF="#binary">Concatenation and elementwise
|
|
multiplication</A>. </LI>
|
|
<LI>A new <A HREF="#unspec">GenericMatrix</A> class.
|
|
</LI>
|
|
<LI><A HREF="#scalar3">Sum</A> function. </LI>
|
|
<LI>Some of the <A HREF="#make">make</A> files reorganised.
|
|
</LI>
|
|
</UL>
|
|
<P>If you are upgrading from <B>newmat06</B> note the following:
|
|
</P>
|
|
<UL>
|
|
<LI>If you are using << to load a Real into a submatrix
|
|
change this to =. </LI>
|
|
</UL>
|
|
<P>If you are upgrading from <B>newmat03</B> or <B>newmat04</B>
|
|
note the following </P>
|
|
<UL>
|
|
<LI>.hxx files are now .h files </LI>
|
|
<LI>real changed to Real </LI>
|
|
<LI>BOOL changed to Boolean </LI>
|
|
<LI>CopyToMatrix changed to AsMatrix, etc </LI>
|
|
<LI>real(A) changed to A.AsScalar() </LI>
|
|
</UL>
|
|
<P>The current version is quite a bit longer that newmat04, so
|
|
if you are almost out of space with newmat04, don't throw newmat04 away until
|
|
you have checked your program will work under this version. </P>
|
|
<P>See the <A HREF="#changes">change history</A> for other changes. </P>
|
|
<h2><a name="except_1"></a>2.6 Catching exceptions</h2>
|
|
<P>This section applies particularly to people using <i>compiler supported</i>
|
|
exceptions rather than my <i>simulated</i> exceptions. </P>
|
|
<P>If newmat detects an error it will throw an exception. It is important that
|
|
you catch this exception and print the error message. Otherwise you will get an
|
|
unhelpful message like <i>abnormal termination</i>. I suggest you set up your
|
|
main program like </P>
|
|
<pre>#define WANT_STREAM // or #include <iostream>
|
|
#include "include.h" // or #include "newmat.h"
|
|
#include "myexcept.h"
|
|
|
|
main()
|
|
{
|
|
try
|
|
{
|
|
... your program here
|
|
}
|
|
// catch exceptions thrown by my programs
|
|
catch(Exception) { cout << Exception::what() << endl; }
|
|
// catch exceptions thrown by other people's programs
|
|
catch(...) { cout << "exception caught in main program" << endl; }
|
|
return 0;
|
|
}</pre>
|
|
<P>If you are using a GUI version rather a console version of the program you
|
|
will need to replace the <i>cout </i>statements by windows pop-up messages.</P>
|
|
<P>If you are using my simulated exceptions or have set the disable exceptions
|
|
option in <i>include.h</i> then uncaught exceptions automatically print the
|
|
error message generated by the exception so you can ignore this section.
|
|
Alternatively use <i>Try</i>, <i>Catch</i> and <i>CatchAll</i> in place of <i>try</i>, <i>catch</i>
|
|
and <i>catch(...)</i> in the preceding code. It is probably a good idea to do
|
|
this if you are using a GUI version of the program as opposed to a console
|
|
version as the <i>cout </i>statement used in newmat's <i>Terminate</i> function
|
|
may be ignored in a GUI version.</P>
|
|
<P>See the <a href="#except">section on exceptions</a> for more information on
|
|
the exception structure. </P>
|
|
<H2><A NAME="example"></A>2.7 Example</H2>
|
|
<P CLASS="small"><A HREF="#testing">next</A> - <A HREF="#testing">skip</A> -
|
|
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
|
|
<P>An example is given in <TT>example.cpp</TT>. This gives a simple linear
|
|
regression example using five different algorithms. The correct output is given
|
|
in <TT>example.txt</TT>. The program carries out a rough check that no memory
|
|
is left allocated on the heap when it terminates. See the section on
|
|
<A HREF="#testing">testing</A> for a comment on the reliability of this check
|
|
(generally it doesn't work with the newer compilers) and the use of the
|
|
DO_FREE_CHECK option. </P>
|
|
<P>I include a variety of make files. To compile the example use a command like
|
|
</P>
|
|
<pre> gmake -f nm_gnu.mak example (Gnu G++)
|
|
gmake -f nm_cc.mak example (AT&T, HPUX, Sun)
|
|
make -f nm_b55.mak example.exe (Borland C++ 5.5)</pre>
|
|
<p>You can generate make make files for a number of other compilers with my
|
|
<a href="genmake.htm">genmake</a> utility.</p>
|
|
|
|
<P>To compile all the example and test files use a command like </P>
|
|
<PRE> gmake -f nm_gnu.mak (Gnu G++)
|
|
</PRE>
|
|
|
|
<P> The example uses io manipulators. It will not work with a
|
|
compiler that does not support the standard io manipulators. </P>
|
|
<P>Other example files are <TT>nl_ex.cpp</TT> and <TT>garch.cpp</TT> for
|
|
demonstrating the non-linear fitting routines, <TT>sl_ex</TT> for demonstrating
|
|
the solve function and <TT>test_exc</TT> for demonstrating the exceptions. </P>
|
|
<H2><A NAME="testing"></A>2.8 Testing</H2>
|
|
<P CLASS="small"><A HREF="#bugs">next</A> - <A HREF="#bugs">skip</A> -
|
|
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
|
|
<P>The library package contains a comprehensive test program in the form of a
|
|
series of files with names of the form tmt?.cxx. The files consist of a large
|
|
number of matrix formulae all of which evaluate to zero (except the first one
|
|
which is used to check that we are detecting non-zero matrices). The printout
|
|
should state that it has found just one non-zero matrix. </P>
|
|
<P>The test program should be run with <I>Real</I> typedefed to <I>double</I>
|
|
rather than <I>float</I> in <A HREF="#custom">include.h</A>. </P>
|
|
<P>Make sure the <A HREF="#elements">C subscripts</A> are enabled if you want
|
|
to test these. </P>
|
|
<P>If you are carrying out some form of bounds checking, for example, with
|
|
Borland's <I>CodeGuard</I>, then disable the testing of the <A
|
|
HREF="#nric">Numerical Recipes in C</A> interface. Activate the statement
|
|
<TT>#define DONT_DO_NRIC</TT> in tmt.h. </P>
|
|
<P>Various versions of the make file (extension .mak) are included with the
|
|
package. See the section on <A HREF="#make">make files</A>. </P>
|
|
<P>The program also allocates and deletes a large block and small block of
|
|
memory before it starts the main testing and then at the end of the test. It
|
|
then checks that the blocks of memory were allocated in the same place. If not,
|
|
then one suspects that there has been a memory leak. i.e. a piece of memory has
|
|
been allocated and not deleted.</P>
|
|
<P>This is not completely foolproof. Programs may allocate extra print buffers
|
|
while the program is running. I have tried to overcome this by doing a print
|
|
before I allocate the first memory block. Programs may allocate memory for
|
|
different sized items in different places, or might not allocate items
|
|
consecutively. Or they might mix the items with memory blocks from other
|
|
programs. Nevertheless, I seem to get consistent answers from <i>some</i> of the
|
|
compilers I work with, so I think this is a worthwhile test. The compilers that
|
|
the test seems to work for include the Borland compilers, Microsoft VC++ 6 ,
|
|
Watcom 10a, and Gnu 2.96 for Linux.</P>
|
|
<P>If the <A HREF="#custom">DO_FREE_CHECK</A> option in include.h is activated,
|
|
the program checks that each <TT>new</TT> is balanced with exactly one
|
|
<TT>delete</TT>. This provides a more definitive test of no memory leaks. There
|
|
are additional statements in myexcept.cpp which can be activated to print out
|
|
details of the memory being allocated and released. </P>
|
|
<P>I have included a facility for checking that each piece of code in the
|
|
library is really exercised by the test routines. Each block of code in the
|
|
main part of the library contains a word <TT>REPORT</TT>. <TT>newmat.h</TT> has
|
|
a line defining <TT>REPORT</TT> that can be activated (deactivate the dummy
|
|
version). This gives a printout of the number of times each of the
|
|
<TT>REPORT</TT> statements in the <TT>.cpp</TT> files is accessed. Use a grep
|
|
with line numbers to locate the lines on which <TT>REPORT</TT> occurs and
|
|
compare these with the lines that the printout shows were actually accessed.
|
|
One can then see which lines of code were not accessed. </P>
|
|
<H2><A NAME="bugs"></A>2.9 Bugs</H2>
|
|
<P CLASS="small"><A HREF="#files">next</A> - <A HREF="#files">skip</A> -
|
|
<A HREF="#starting">up</A> - <A HREF="#top">start</A> </P>
|
|
<UL>
|
|
<LI>Small memory leaks may occur when an exception is thrown and caught. </LI>
|
|
<LI>My exception scheme may not be not properly linked in with the standard
|
|
library exceptions. In particular, my scheme may fail to catch out-of-memory
|
|
exceptions. </LI>
|
|
</UL>
|
|
<H2><A NAME="files"></A>2.10 List of files</H2>
|
|
<P CLASS="small"><A HREF="#problem">next</A> - <A HREF="#problem">skip</A> -
|
|
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
|
|
<TABLE WIDTH="100%" BORDER="0" cellspacing="0" style="border-collapse: collapse" bordercolor="#111111" cellpadding="0">
|
|
<TR>
|
|
<TD width="25%"><B>Documentation</B></TD>
|
|
<TD width="25%">README(.txt)</TD>
|
|
<TD width="50%">readme file</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">COPYING.(.txt)</TD>
|
|
<TD width="50%">permission to copy, modify etc</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">AUTHORS(.txt)</TD>
|
|
<TD width="50%">author details</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">nm10.htm</TD>
|
|
<TD width="50%">documentation file</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">add_time.pgn</TD>
|
|
<TD width="50%">image used by nm10.htm - move to subdirectory "images"</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">rbd.css</TD>
|
|
<TD width="50%">style sheet for nm10.htm</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"><B>Definition files</B></TD>
|
|
<TD width="25%">boolean.h</TD>
|
|
<TD width="50%">boolean class definition</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">controlw.h</TD>
|
|
<TD width="50%">control word definition file</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">include.h</TD>
|
|
<TD width="50%">details of include files and options</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">myexcept.h</TD>
|
|
<TD width="50%">general exception handler definitions</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmat.h</TD>
|
|
<TD width="50%">main matrix class definition file</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmatap.h</TD>
|
|
<TD width="50%">applications definition file</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmatio.h</TD>
|
|
<TD width="50%">input/output definition file</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmatnl.h</TD>
|
|
<TD width="50%">non-linear optimisation definition file</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmatrc.h</TD>
|
|
<TD width="50%">row/column functions definition files</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmatrm.h</TD>
|
|
<TD width="50%">rectangular matrix access definition files</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">precisio.h</TD>
|
|
<TD width="50%">numerical precision constants</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">solution.h</TD>
|
|
<TD width="50%">one dimensional solve definition file</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"><B>Program files</B></TD>
|
|
<TD width="25%">bandmat.cpp</TD>
|
|
<TD width="50%">band matrix routines</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">cholesky.cpp</TD>
|
|
<TD width="50%">Cholesky decomposition</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">evalue.cpp</TD>
|
|
<TD width="50%">eigenvalues and eigenvector calculation</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">fft.cpp</TD>
|
|
<TD width="50%">fast Fourier, trig. transforms</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">hholder.cpp</TD>
|
|
<TD width="50%">QR routines</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">jacobi.cpp</TD>
|
|
<TD width="50%">eigenvalues by the Jacobi method</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">myexcept.cpp</TD>
|
|
<TD width="50%">general error and exception handler</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newfft.cpp</TD>
|
|
<TD width="50%">new fast Fourier transform</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmat1.cpp</TD>
|
|
<TD width="50%">type manipulation routines</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmat2.cpp</TD>
|
|
<TD width="50%">row and column manipulation functions</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmat3.cpp</TD>
|
|
<TD width="50%">row and column access functions</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmat4.cpp</TD>
|
|
<TD width="50%">constructors, resize, utilities</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmat5.cpp</TD>
|
|
<TD width="50%">transpose, evaluate, matrix functions</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmat6.cpp</TD>
|
|
<TD width="50%">operators, element access</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmat7.cpp</TD>
|
|
<TD width="50%">invert, solve, binary operations</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmat8.cpp</TD>
|
|
<TD width="50%">LU decomposition, scalar functions</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmat9.cpp</TD>
|
|
<TD width="50%">output routines</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmatex.cpp</TD>
|
|
<TD width="50%">matrix exception handler</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmatnl.cpp</TD>
|
|
<TD width="50%">non-linear optimisation</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">newmatrm.cpp</TD>
|
|
<TD width="50%">rectangular matrix access functions</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">sort.cpp</TD>
|
|
<TD width="50%">sorting functions</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">solution.cpp</TD>
|
|
<TD width="50%">one dimensional solve</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">submat.cpp</TD>
|
|
<TD width="50%">submatrix functions</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">svd.cpp</TD>
|
|
<TD width="50%">singular value decomposition</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"><B>Example files</B></TD>
|
|
<TD width="25%">example.cpp</TD>
|
|
<TD width="50%">example of use of package</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">example.txt</TD>
|
|
<TD width="50%">output from example</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">sl_ex.cpp</TD>
|
|
<TD width="50%">example of OneDimSolve routine</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">sl_ex.txt</TD>
|
|
<TD width="50%">output from example</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">nl_ex.cpp</TD>
|
|
<TD width="50%">example of non-linear least squares</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">nl_ex.txt</TD>
|
|
<TD width="50%">output from example</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">garch.cpp</TD>
|
|
<TD width="50%">example of maximum-likelihood fit</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">garch.dat</TD>
|
|
<TD width="50%">data file for garch.cpp</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">garch.txt</TD>
|
|
<TD width="50%">output from example</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">test_exc.cpp</TD>
|
|
<TD width="50%">demonstration exceptions</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">test_exc.txt</TD>
|
|
<TD width="50%">output from test_exc.cpp</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"><B>Test files</B></TD>
|
|
<TD width="25%">tmt.h</TD>
|
|
<TD width="50%">header file for test files</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">tmt*.cpp</TD>
|
|
<TD width="50%">test files (see the section on <A
|
|
HREF="#testing">testing</A>)</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">tmt.txt</TD>
|
|
<TD width="50%">output from test files</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"><B>Make files</B></TD>
|
|
<TD width="25%">nm_gnu.mak</TD>
|
|
<TD width="50%">make file for Gnu G++</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">nm_cc.mak</TD>
|
|
<TD width="50%">make file for AT&T, Sun and HPUX</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">nm_b55.mak</TD>
|
|
<TD width="50%">make file for Borland C++ 5.5</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">nm_m6.mak</TD>
|
|
<TD width="50%">make file for Microsoft Visual C++ 6, 7, 7.1</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">nm_m8.mak</TD>
|
|
<TD width="50%">make file for Microsoft Visual C++ 8</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">nm_i8.mak</TD>
|
|
<TD width="50%">make file for Intel C++ 8,9 for Windows</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">nm_il8.mak</TD>
|
|
<TD width="50%">make file for Intel C++ 8,9 for Linux</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">nm_ow.mak</TD>
|
|
<TD width="50%">make file for Open Watcom</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%" height="22"> </TD>
|
|
<TD width="25%" height="22">newmat.lfl</TD>
|
|
<TD width="50%" height="22">library file list for use with <a href="genmake.htm">genmake</a></TD>
|
|
</TR>
|
|
<TR>
|
|
<TD width="25%"> </TD>
|
|
<TD width="25%">nm_targ.txt</TD>
|
|
<TD width="50%">target file list for use with <a href="genmake.htm">genmake</a></TD>
|
|
</TR>
|
|
</TABLE>
|
|
<H2><A NAME="problem"></A>2.11 Problem report form
|
|
</H2>
|
|
<P CLASS="small"><A HREF="#refer">next</A> - <A HREF="#refer">skip</A> -
|
|
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
|
|
<P>Copy and paste this to your editor; fill it out and email to
|
|
<B>robert at statsresearch.co.nz</B> </P>
|
|
<P>But first look in my web page <a href="http://www.robertnz.net/bugs">http://www.robertnz.net/bugs.htm</a> to see if the bug has
|
|
already been reported. </P>
|
|
<PRE> Version: ............... newmat10D (2 April 2006)
|
|
Your email address: ....
|
|
Today's date: ..........
|
|
Your machine: ..........
|
|
Operating system: ......
|
|
Compiler & version: ....
|
|
Compiler options
|
|
(eg GUI or console)...
|
|
Describe the problem - attach examples if possible:
|
|
|
|
|
|
|
|
|
|
|
|
-----------------------------------------------------------
|
|
</PRE>
|
|
|
|
<H2><A NAME="refer"></A>3. Reference manual</H2>
|
|
<P CLASS="small"><A HREF="#constr">next</A> - <A HREF="#error">skip</A> -
|
|
<A HREF="#top">up</A> - <A HREF="#top">start</A></P>
|
|
<TABLE WIDTH="100%">
|
|
<TR>
|
|
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A HREF="#constr">3.1
|
|
Constructors </A><BR>
|
|
<A HREF="#elements">3.2 Accessing elements </A><BR>
|
|
<A HREF="#copy">3.3 Assignment and copying </A><BR>
|
|
<A HREF="#entering">3.4 Entering values </A><BR>
|
|
<A HREF="#unary">3.5 Unary operations </A><BR>
|
|
<A HREF="#binary">3.6 Binary operations </A><BR>
|
|
<A HREF="#matscal">3.7 Matrix and scalar ops </A><BR>
|
|
<A HREF="#scalar1">3.8 Scalar functions - size & shape </A><BR>
|
|
<A HREF="#scalar2">3.9 Scalar functions - maximum & minimum </A><BR>
|
|
<A HREF="#scalar3">3.10 Scalar functions - numerical </A><BR>
|
|
<A HREF="#submat">3.11 Submatrices </A><BR>
|
|
<A HREF="#dimen">3.12 Change dimension </A><BR>
|
|
<A HREF="#ch_type">3.13 Change type </A><BR>
|
|
<A HREF="#solve">3.14 Multiple matrix solve </A><BR>
|
|
<A HREF="#memory">3.15 Memory management </A><BR>
|
|
<A HREF="#efficien">3.16 Efficiency </A></TD>
|
|
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A
|
|
HREF="#output">3.17 Output </A><BR>
|
|
<A HREF="#unspec">3.18 Accessing unspecified type </A><BR>
|
|
<A HREF="#cholesky">3.19 Cholesky decomposition </A><BR>
|
|
<A HREF="#qr">3.20 QR decomposition </A><BR>
|
|
<A HREF="#svd">3.21 Singular value decomposition </A><BR>
|
|
<A HREF="#evalues">3.22 Eigenvalue decomposition </A><BR>
|
|
<A HREF="#sorting">3.23 Sorting </A><BR>
|
|
<A HREF="#fft">3.24 Fast Fourier transform </A><BR>
|
|
<A HREF="#trigtran">3.25 Fast trigonometric transforms </A><BR>
|
|
<A HREF="#nric">3.26 Numerical recipes in C </A><BR>
|
|
<A HREF="#except">3.27 Exceptions </A><BR>
|
|
<A HREF="#cleanup">3.28 Cleanup following exception </A><BR>
|
|
<A HREF="#nonlin">3.29 Non-linear applications </A><BR>
|
|
<A HREF="#stl">3.30 Standard template library </A><BR>
|
|
<A HREF="#namesp">3.31 Namespace </A></TD>
|
|
</TR>
|
|
</TABLE>
|
|
<H2><A NAME="constr"></A>3.1 Constructors</H2>
|
|
<P CLASS="small"><A HREF="#elements">next</A> - <A HREF="#elements">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>To construct an <I>m</I> x <I>n</I> matrix, <TT>A</TT>, (<I>m</I> and
|
|
<I>n</I> are integers) use </P>
|
|
<PRE> Matrix A(m,n);
|
|
</PRE>
|
|
|
|
<P>The UpperTriangularMatrix, LowerTriangularMatrix, SymmetricMatrix and
|
|
DiagonalMatrix types are square. To construct an <I>n</I> x <I>n</I> matrix
|
|
use, for example </P>
|
|
<PRE> UpperTriangularMatrix UT(n);
|
|
LowerTriangularMatrix LT(n);
|
|
SymmetricMatrix S(n);
|
|
DiagonalMatrix D(n);
|
|
</PRE>
|
|
|
|
<P>Band matrices need to include bandwidth information in their constructors.
|
|
</P>
|
|
<PRE> BandMatrix BM(n, lower, upper);
|
|
UpperBandMatrix UB(n, upper);
|
|
LowerBandMatrix LB(n, lower);
|
|
SymmetricBandMatrix SB(n, lower);
|
|
</PRE>
|
|
|
|
<P>The integers <I>upper</I> and <I>lower</I> are the number of non-zero
|
|
diagonals above and below the diagonal (<I>excluding</I> the diagonal)
|
|
respectively. </P>
|
|
<P>The RowVector and ColumnVector types take just one argument in their
|
|
constructors: </P>
|
|
<PRE> RowVector RV(n);
|
|
ColumnVector CV(n);
|
|
</PRE>
|
|
|
|
<P><b>These constructors do <EM>not</EM> initialise the elements of the matrices.
|
|
</b>
|
|
To set all the elements to zero use, for example, </P>
|
|
<PRE> Matrix A(m, n); A = 0.0;
|
|
</PRE>
|
|
|
|
<P>The IdentityMatrix takes one argument in its constructor specifying its
|
|
dimension.</P>
|
|
<pre> IdentityMatrix I(n);</pre>
|
|
<p>The value of the diagonal elements <b>is</b> set to 1 by default, but you can
|
|
change this value as with other matrix types. </p>
|
|
|
|
<P>You can also construct vectors and matrices without specifying the
|
|
dimension. For example </P>
|
|
<PRE> Matrix A;
|
|
</PRE>
|
|
|
|
<P>In this case the dimension must be set by an <A HREF="#copy">assignment
|
|
statement</A> or a <A HREF="#dimen">re-size statement</A>. </P>
|
|
<P>You can also use a constructor to set a matrix equal to another matrix or
|
|
matrix expression. </P>
|
|
<PRE> Matrix A = UT;
|
|
Matrix A = UT * LT;
|
|
</PRE>
|
|
|
|
<P>Only conversions that don't lose information are supported - eg you cannot
|
|
convert an upper triangular matrix into a diagonal matrix using =. </P>
|
|
<H2><A NAME="elements"></A>3.2 Accessing elements
|
|
</H2>
|
|
<P CLASS="small"><A HREF="#copy">next</A> - <A HREF="#copy">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>Elements are accessed by expressions of the form <TT>A(i,j)</TT> where
|
|
<I>i</I> and <I>j</I> run from 1 to the appropriate dimension. Access elements
|
|
of vectors with just one argument. Diagonal matrices can accept one or two
|
|
subscripts. </P>
|
|
<P>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
|
|
<TT>A.element(i,j)</TT> if you want this earlier convention. </P>
|
|
<P><TT>A(i,j)</TT> and <TT>A.element(i,j)</TT> can appear on either side of an
|
|
= sign. </P>
|
|
<P>If you activate the <TT>#define SETUP_C_SUBSCRIPTS</TT> in
|
|
<TT>include.h</TT> you can also access elements using the traditional C style
|
|
notation. That is <TT>A[i][j]</TT> for matrices (except diagonal) and
|
|
<TT>V[i]</TT> for vectors and diagonal matrices. The subscripts start at zero
|
|
(i.e. like element) and there is <I>no</I> range checking. Because of the
|
|
possibility of confusing <TT>V(i)</TT> and <TT>V[i]</TT>, I suggest you do
|
|
<I>not</I> activate this option unless you really want to use it. </P>
|
|
<P>Symmetric matrices are stored as lower triangular matrices. It is important
|
|
to remember this if you are using the <TT>A[i][j]</TT> method of accessing
|
|
elements. Make sure the first subscript is greater than or equal to the second
|
|
subscript. However, if you are using the <TT>A(i,j)</TT> method the program
|
|
will swap <TT>i</TT> and <TT>j</TT> if necessary; so it doesn't matter if you
|
|
think of the storage as being in the upper triangle (but it <I>does</I> matter
|
|
in some other situations such as when <A HREF="#entering">entering</A> data).</P>
|
|
<P>The IdentityMatrix type does not support element access.
|
|
</P>
|
|
<H2><A NAME="copy"></A>3.3 Assignment and copying</H2>
|
|
<P CLASS="small"><A HREF="#entering">next</A> - <A HREF="#entering">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>The operator <TT>=</TT> is used for copying matrices, converting matrices,
|
|
or evaluating expressions. For example </P>
|
|
<PRE> A = B; A = L; A = L * U;
|
|
</PRE>
|
|
|
|
<P>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. </P>
|
|
<P>The operator <TT><<</TT> can be used in place of <TT>=</TT> where it
|
|
is permissible for information to be lost. </P>
|
|
<P>For example </P>
|
|
<PRE> SymmetricMatrix S; Matrix A;
|
|
......
|
|
S << A.t() * A;
|
|
</PRE>
|
|
|
|
<P>is acceptable whereas </P>
|
|
<PRE> S = A.t() * A; // error
|
|
</PRE>
|
|
|
|
<P>will cause a runtime error since the package does not (yet?) recognise
|
|
<TT>A.t()*A</TT> as symmetric. </P>
|
|
<P>Note that you can <I>not</I> use <TT><<</TT> with constructors. For
|
|
example </P>
|
|
<PRE> SymmetricMatrix S << A.t() * A; // error
|
|
</PRE>
|
|
|
|
<P>does <I>not</I> work. </P>
|
|
<P>Also note that <TT><<</TT> 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. </P>
|
|
<P>A third copy routine is used in a similar role to <TT>=</TT>. Use </P>
|
|
<PRE> A.Inject(D);
|
|
</PRE>
|
|
|
|
<P>to copy the elements of <TT>D</TT> to the corresponding elements of
|
|
<TT>A</TT> but leave the elements of <TT>A</TT> unchanged if there is no
|
|
corresponding element of <TT>D</TT> (the <TT>=</TT> 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 <TT>=</TT> and
|
|
<TT><<</TT>, Inject does not reset the dimensions of <TT>A</TT>, which
|
|
must match those of <TT>D</TT>. Inject does <I>not</I> test for no loss of
|
|
information. </P>
|
|
<P>You cannot replace <TT>D</TT> by a matrix expression. The effect of
|
|
<TT>Inject(D)</TT> depends on the type of <TT>D</TT>. If <TT>D</TT> is an
|
|
expression it might not be obvious to the user what type it would have. So I
|
|
thought it best to disallow expressions. </P>
|
|
<P>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). </P>
|
|
<P>Both <TT><<</TT> and Inject can be used with submatrix expressions on
|
|
the left hand side. See the section on <A HREF="#submat">submatrices</A>. </P>
|
|
<P>To set the elements of a matrix to a scalar use operator <TT>=</TT> </P>
|
|
<PRE> Real r; int m,n;
|
|
......
|
|
Matrix A(m,n); A = r;</PRE>
|
|
|
|
<p>Notes:</p>
|
|
|
|
<ul>
|
|
<li>When you do a matrix assignment to another matrix or matrix expression with
|
|
either <tt>=</tt> or <tt><<</tt> the original data array associated with the matrix
|
|
being assigned to is destroyed
|
|
even if there is no change in length. See the section on <a href="#stor">storage</a>.
|
|
This means, in particular, that pointers to matrix elements - e.g.
|
|
<tt>Real* a; a = &(A(1,1));</tt> become invalid. If you want avoid this you
|
|
can use
|
|
<tt>Inject</tt> rather than <tt>=</tt>. But remember that you may need to zero
|
|
the matrix first.<br>
|
|
</li>
|
|
</ul>
|
|
|
|
<H2><A NAME="entering"></A>3.4 Entering values</H2>
|
|
<P CLASS="small"><A HREF="#unary">next</A> - <A HREF="#unary">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>You can load the elements of a matrix from an array: </P>
|
|
<PRE> Matrix A(3,2);
|
|
Real a[] = { 11,12,21,22,31,33 };
|
|
A << a;
|
|
</PRE>
|
|
|
|
<P>This construction does <I>not</I> check that the numbers of elements match
|
|
correctly. This version of <TT><<</TT> can be used with submatrices on
|
|
the left hand side. It is not defined for band matrices. </P>
|
|
<P>Alternatively you can enter short lists using a sequence of numbers
|
|
separated by <TT><<</TT> . </P>
|
|
<PRE> Matrix A(3,2);
|
|
A << 11 << 12
|
|
<< 21 << 22
|
|
<< 31 << 32;
|
|
</PRE>
|
|
|
|
<P>This does check for the correct total number of entries, although the
|
|
message for there being insufficient numbers in the list may be delayed until
|
|
the end of the block or the next use of this construction. This does <I>not</I>
|
|
work for band matrices or for long lists. It does work for submatrices if the
|
|
submatrix is a single complete row. For example </P>
|
|
<PRE> Matrix A(3,2);
|
|
A.Row(1) << 11 << 12;
|
|
A.Row(2) << 21 << 22;
|
|
A.Row(3) << 31 << 32;
|
|
</PRE>
|
|
|
|
<P>Load only values that are actually stored in the matrix. For example </P>
|
|
<PRE> SymmetricMatrix S(2);
|
|
S.Row(1) << 11;
|
|
S.Row(2) << 21 << 22;
|
|
</PRE>
|
|
|
|
<P>Try to restrict this way of loading data to numbers. You can include
|
|
expressions, but these must not call a function which includes the same
|
|
construction. </P>
|
|
|
|
<p>Remember that matrices are stored by rows and that symmetric matrices are
|
|
stored as <I> lower</I> triangular matrices when using these methods to enter
|
|
data. </p>
|
|
<H2><A NAME="unary"></A>3.5 Unary operators</H2>
|
|
<P CLASS="small"><A HREF="#binary">next</A> - <A HREF="#binary">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>The package supports unary operations </P>
|
|
<PRE> X = -A; // change sign of elements
|
|
X = A.t(); // transpose
|
|
X = A.i(); // inverse (of square matrix A)
|
|
X = A.Reverse(); // reverse order of elements of vector
|
|
// or matrix (not band matrix)
|
|
</PRE>
|
|
|
|
<H2><A NAME="binary"></A>3.6 Binary operators</H2>
|
|
<P CLASS="small"><A HREF="#matscal">next</A> - <A HREF="#matscal">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>The package supports binary operations </P>
|
|
<PRE> X = A + B; // matrix addition
|
|
X = A - B; // matrix subtraction
|
|
X = A * B; // matrix multiplication
|
|
X = A.i() * B; // equation solve (square matrix A)
|
|
X = A | B; // concatenate horizontally (concatenate the rows)
|
|
X = A & B; // concatenate vertically (concatenate the columns)
|
|
X = SP(A, B); // elementwise product of A and B (Schur product)
|
|
X = KP(A, B); // Kronecker product of A and B
|
|
bool b = A == B; // test whether A and B are equal
|
|
bool b = A != B; // ! (A == B)
|
|
A += B; // A = A + B;
|
|
A -= B; // A = A - B;
|
|
A *= B; // A = A * B;
|
|
A |= B; // A = A | B;
|
|
A &= B; // A = A & B;
|
|
<, >, <=, >= // included for compatibility with STL - see notes
|
|
</PRE>
|
|
|
|
<P>Notes: </P>
|
|
<UL>
|
|
<LI>If you are doing repeated multiplication. For example <TT>A*B*C</TT>, use
|
|
brackets to force the order of evaluation to minimise the number of operations.
|
|
If <TT>C</TT> is a column vector and <TT>A</TT> is not a vector, then it will
|
|
usually reduce the number of operations to use <TT>A*(B*C)</TT>. </LI>
|
|
<LI>In the equation solve example case the inverse is not explicitly
|
|
calculated. An LU decomposition of <TT>A</TT> is performed and this is applied
|
|
to <TT>B</TT>. This is more efficient than calculating the inverse and then
|
|
multiplying. See also <A HREF="#solve">multiple matrix solving</A>. </LI>
|
|
<LI>The package does not (yet?) recognise <TT>B*A.i()</TT> as an equation solve
|
|
and the inverse of <TT>A</TT> would be calculated. It is probably better to use
|
|
<TT>(A.t().i()*B.t()).t()</TT>. </LI>
|
|
<LI>Horizontal or vertical concatenation returns a result of type Matrix,
|
|
RowVector or ColumnVector. </LI>
|
|
<LI>If <TT>A</TT> is <I> m</I> x <I>p</I>, <TT>B</TT> is <I> m</I> x <I>q</I>,
|
|
then <TT>A | B</TT> is <I> m</I> x (<I>p</I>+<I>q</I>) with the <I>k</I>-th row
|
|
being the elements of the <I>k</I>-th row of <TT>A</TT> followed by the
|
|
elements of the <I>k</I>-th row of <TT>B</TT>. </LI>
|
|
<LI>If <TT>A</TT> is <I> p</I> x <I>n</I>, <TT>B</TT> is <I> q</I> x <I>n</I>,
|
|
then <TT>A & B</TT> is (<I>p</I>+<I>q</I>) x <I> n</I> with the <I>k</I>-th
|
|
column being the elements of the <I>k</I>-th column of <TT>A</TT> followed by
|
|
the elements of the <I>k</I>-th column of <TT>B</TT>. </LI>
|
|
<LI>For complicated concatenations of matrices, consider instead using
|
|
<A HREF="#submat">submatrices</A>. </LI>
|
|
<LI>See the section on <A HREF="#submat">submatrices</A> on using a submatrix
|
|
on the RHS of an expression. </LI>
|
|
<LI>Two matrices are equal if their difference is zero. They may be of
|
|
different types. For the CroutMatrix or BandLUMatrix they must be of the same
|
|
type and have all their elements equal. This is not a very useful operator and
|
|
is included for compatibility with some container templates. </LI>
|
|
<LI>The inequality operators are included for compatibility with the
|
|
<A HREF="#stl">standard template library</A>. If actually called, they will
|
|
throw an exception. So don't try to sort a <I>list</I> of matrices. </LI>
|
|
<LI> A row vector multiplied by a column vector yields a 1x1 matrix, <I>not</I>
|
|
a Real. To get a Real result use either <A HREF="#scalar3"> AsScalar or
|
|
DotProduct</A>.</LI>
|
|
<LI> The result from Kronecker product, KP(A, B), possesses an attribute such as
|
|
upper triangular, lower triangular, band, symmetric, diagonal if both of the
|
|
matrices A and B have the attribute. (This differs slightly from the way the
|
|
January 2002 version of newmat10 worked).</LI>
|
|
</UL>
|
|
<P CLASS="small">Remember that the product of symmetric matrices is not
|
|
necessarily symmetric so the following code will not run:</P>
|
|
<PRE> SymmetricMatrix A, B;
|
|
.... put values in A, B ....
|
|
SymmetricMatrix C = A * B; // run time error</PRE>
|
|
|
|
<P CLASS="small">Use instead</P>
|
|
<PRE> Matrix C = A * B;</PRE>
|
|
|
|
<P CLASS="small">or, if you <EM>know</EM> the product will be symmetric,
|
|
use</P>
|
|
<PRE> SymmetricMatrix C; C << A * B;
|
|
</PRE>
|
|
|
|
<H2><A NAME="matscal"></A>3.7 Matrix and scalar</H2>
|
|
<P CLASS="small"><A HREF="#scalar1">next</A> - <A HREF="#scalar1">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>The following expressions multiply the elements of a matrix <TT>A</TT> by a
|
|
scalar <TT>f</TT>: <TT>A * f</TT> or <TT>f * A</TT> . Likewise one can divide
|
|
the elements of a matrix <TT>A</TT> by a scalar <TT>f</TT>: <TT>A / f</TT> .
|
|
</P>
|
|
<P>The expressions <TT>A + f</TT> and <TT>A - f</TT> add or subtract a
|
|
rectangular matrix of the same dimension as <TT>A</TT> with elements equal to
|
|
<TT>f</TT> to or from the matrix <TT>A</TT> . </P>
|
|
<P>The expression <TT>f + A</TT> is an alternative to <TT>A + f</TT>. The
|
|
expression <TT>f - A</TT> subtracts matrix <TT>A</TT> from a rectangular matrix
|
|
of the same dimension as <TT>A</TT> and with elements equal to <TT>f</TT> .
|
|
</P>
|
|
<P>The expression <TT>A += f</TT> replaces <TT>A</TT> by <TT>A + f</TT>.
|
|
Operators <TT>-=</TT>, <TT>*=</TT>, <TT>/=</TT> are similarly defined. </P>
|
|
<H2><A NAME="scalar1"></A>3.8 Scalar functions of a
|
|
matrix - size & shape</H2>
|
|
<P CLASS="small"><A HREF="#scalar2">next</A> - <A HREF="#scalar2">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>This page describes functions returning the values associated with the size
|
|
and shape of matrices. The following pages describe other scalar matrix
|
|
functions. </P>
|
|
<PRE> int m = A.Nrows(); // number of rows
|
|
int n = A.Ncols(); // number of columns
|
|
MatrixType mt = A.Type(); // type of matrix
|
|
Real* s = A.Store(); // pointer to array of elements
|
|
int l = A.Storage(); // length of array of elements
|
|
MatrixBandWidth mbw = A.BandWidth(); // upper and lower bandwidths
|
|
</PRE>
|
|
|
|
<P><TT>MatrixType mt = A.Type()</TT> returns the type of a matrix. Use
|
|
<TT>mt.Value()</TT> to get a string (UT, LT, Rect, Sym, Diag, Band, UB, LB,
|
|
Crout, BndLU) showing the type (Vector types are returned as Rect). </P>
|
|
<P><TT>MatrixBandWidth</TT> has member functions <TT>Upper()</TT> and
|
|
<TT>Lower()</TT> for finding the upper and lower bandwidths (number of
|
|
diagonals above and below the diagonal, both zero for a diagonal matrix). For
|
|
non-band matrices -1 is returned for both these values. </P>
|
|
<H2><A NAME="scalar2"></A>3.9 Scalar functions of a
|
|
matrix - maximum & minimum</H2>
|
|
<P CLASS="small"><A HREF="#scalar3">next</A> - <A HREF="#scalar3">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>This page describes functions for finding the maximum and minimum elements
|
|
of a matrix. </P>
|
|
<PRE> int i, j;
|
|
Real mv = A.MaximumAbsoluteValue(); // maximum of absolute values
|
|
Real mv = A.MinimumAbsoluteValue(); // minimum of absolute values
|
|
Real mv = A.Maximum(); // maximum value
|
|
Real mv = A.Minimum(); // minimum value
|
|
Real mv = A.MaximumAbsoluteValue1(i); // maximum of absolute values
|
|
Real mv = A.MinimumAbsoluteValue1(i); // minimum of absolute values
|
|
Real mv = A.Maximum1(i); // maximum value
|
|
Real mv = A.Minimum1(i); // minimum value
|
|
Real mv = A.MaximumAbsoluteValue2(i,j);// maximum of absolute values
|
|
Real mv = A.MinimumAbsoluteValue2(i,j);// minimum of absolute values
|
|
Real mv = A.Maximum2(i,j); // maximum value
|
|
Real mv = A.Minimum2(i,j); // minimum value
|
|
</PRE>
|
|
|
|
<P>All these functions throw an exception if <TT>A</TT> has no rows or no
|
|
columns. </P>
|
|
<P>The versions <TT>A.MaximumAbsoluteValue1(i)</TT>, etc return the location of
|
|
the extreme element in a RowVector, ColumnVector or DiagonalMatrix. The
|
|
versions <TT>A.MaximumAbsoluteValue2(i,j)</TT>, etc return the row and column
|
|
numbers of the extreme element. If the extreme value occurs more than once the
|
|
location of the last one is given. </P>
|
|
<P>The versions MaximumAbsoluteValue(A), MinimumAbsoluteValue(A), Maximum(A),
|
|
Minimum(A) can be used in place of A.MaximumAbsoluteValue(),
|
|
A.MinimumAbsoluteValue(), A.Maximum(), A.Minimum(). </P>
|
|
<H2><A NAME="scalar3"></A>3.10 Scalar functions of a
|
|
matrix - numerical</H2>
|
|
<P CLASS="small"><A HREF="#submat">next</A> - <A HREF="#submat">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<PRE> Real r = A.AsScalar(); // value of 1x1 matrix
|
|
Real ssq = A.SumSquare(); // sum of squares of elements
|
|
Real sav = A.SumAbsoluteValue(); // sum of absolute values
|
|
Real s = A.Sum(); // sum of 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 norm = A.NormFrobenius(); // square root of sum of squares
|
|
// of the elements
|
|
Real t = A.Trace(); // trace
|
|
Real d = A.Determinant(); // determinant
|
|
LogAndSign ld = A.LogDeterminant(); // log of determinant
|
|
bool z = A.IsZero(); // test all elements zero
|
|
bool s = A.IsSingular(); // A is a CroutMatrix or
|
|
BandLUMatrix
|
|
Real s = DotProduct(A, B); // dot product of A and B
|
|
// interpreted as vectors
|
|
</PRE>
|
|
|
|
<P><TT>A.LogDeterminant()</TT> returns a value of type LogAndSign. If ld is of
|
|
type LogAndSign use </P>
|
|
<PRE> 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.
|
|
</PRE>
|
|
|
|
<P>Note that the direct use of the function <TT>Determinant()</TT> will often
|
|
cause a floating point overflow exception. </P>
|
|
<P><TT>A.IsZero()</TT> returns Boolean value <TT>true</TT> if the matrix
|
|
<TT>A</TT> has all elements equal to 0.0. </P>
|
|
<P><TT>IsSingular</TT> is defined only for CroutMatrix and BandLUMatrix. It
|
|
returns <TT>true</TT> if one of the diagonal elements of the LU decomposition
|
|
is exactly zero. </P>
|
|
<P><TT>DotProduct(const Matrix& A,const Matrix& B)</TT> converts both
|
|
of the arguments to rectangular matrices, checks that they have the same number
|
|
of elements and then calculates the first element of <TT>A * </TT>first element
|
|
of <TT>B + </TT>second element of <TT>A * </TT>second element of <TT>B +
|
|
...</TT> ignoring the row/column structure of <TT>A</TT> and <TT>B</TT>. It is
|
|
primarily intended for the situation where <TT>A</TT> and <TT>B</TT> are row or
|
|
column vectors. </P>
|
|
<P>The versions Sum(A), SumSquare(A), SumAbsoluteValue(A), Trace(A),
|
|
LogDeterminant(A), Determinant(A), Norm1(A), NormInfinity(A), NormFrobenius(A)
|
|
can be used in place of A.Sum(), A.SumSquare(), A.SumAbsoluteValue(),
|
|
A.Trace(), A.LogDeterminant(), A.Norm1(), A.NormInfinity(), A.NormFrobenius().
|
|
</P>
|
|
<H2><A NAME="submat"></A>3.11 Submatrices</H2>
|
|
<P CLASS="small"><A HREF="#dimen">next</A> - <A HREF="#dimen">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<PRE> A.SubMatrix(fr,lr,fc,lc)
|
|
</PRE>
|
|
|
|
<P>This selects a submatrix from <TT>A</TT>. The arguments
|
|
<TT>fr</TT>,<TT>lr</TT>,<TT>fc</TT>,<TT>lc</TT> are the first row, last row,
|
|
first column, last column of the submatrix with the numbering beginning at 1.
|
|
</P>
|
|
<P>I allow <TT>lr = fr-1</TT> or <TT>lc = fc-1</TT> or to indicate that a
|
|
matrix of zero rows or columns is to be returned. </P>
|
|
<P>A submatrix command may be used in any matrix expression or on the left hand
|
|
side of <TT>=</TT>, <TT><<</TT> or Inject. Inject does <I>not</I> check
|
|
no information loss. You can also use the construction </P>
|
|
<PRE> Real c; .... A.SubMatrix(fr,lr,fc,lc) = c;
|
|
</PRE>
|
|
|
|
<P>to set a submatrix equal to a constant. </P>
|
|
<P>The following are variants of SubMatrix: </P>
|
|
<PRE> 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
|
|
</PRE>
|
|
|
|
<P>In each case <TT>f</TT> and <TT>l</TT> mean the first and last row or column
|
|
to be selected (starting at 1). </P>
|
|
<P>I allow <TT>l = f-1</TT> to indicate that a matrix of zero rows or columns
|
|
is to be returned. </P>
|
|
<P>If SubMatrix or its variant occurs on the right hand side of an <TT>=</TT>
|
|
or <TT><<</TT> or within an expression think of its type as follows </P>
|
|
<PRE> 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
|
|
</PRE>
|
|
|
|
<P>If SubMatrix or its variant appears on the left hand side of <TT>=</TT> or
|
|
<TT><<</TT> , think of its type being Matrix. Thus <TT>L.Row(1)</TT>
|
|
where <TT>L</TT> is LowerTriangularMatrix expects <TT>L.Ncols()</TT> elements
|
|
even though it will use only one of them. If you are using <TT>=</TT> the
|
|
program will check for no loss of data.</P>
|
|
|
|
<P>A SubMatrix can appear on the left-hand side of <TT>+=</TT> or <TT>-=</TT>
|
|
with a matrix expression on the right-hand side. It can also appear on the
|
|
left-hand side of <TT>+=</TT>, <TT>-=</TT>, <TT>*=</TT> or <TT>/=</TT> with a
|
|
Real on the right-hand side. In each case there must be no loss of information.
|
|
</P>
|
|
<P>The <TT>Row</TT> version can appear on the left hand side of
|
|
<TT><<</TT> for <A HREF="#entering">loading literal data</A> into a row.
|
|
Load only the number of elements that are actually going to be stored in
|
|
memory. </P>
|
|
<P>Do not use the <TT>+=</TT> and <TT>-=</TT> operations with a submatrix of a
|
|
SymmetricMatrix or BandSymmetricMatrix on the LHS and a Real on the RHS.</P>
|
|
<P>You can't pass a submatrix (or any of its variants) as a reference
|
|
non-constant matrix in a function argument. For example, the following will not
|
|
work:</P>
|
|
<pre> void YourFunction(Matrix& A);
|
|
...
|
|
Matrix B(10,10);
|
|
YourFunction(B.SubMatrix(1,5,1,5)) // won't compile</pre>
|
|
<P>If you are are using the submatrix facility to build a matrix from a small
|
|
number of components, consider instead using the <A
|
|
HREF="#binary">concatenation operators</A>. </P>
|
|
<H2><A NAME="dimen"></A>3.12 Change dimensions</H2>
|
|
<P CLASS="small"><A HREF="#ch_type">next</A> - <A HREF="#ch_type">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>The following operations change the dimensions of a matrix. The values of
|
|
the elements are lost. </P>
|
|
<PRE> A.ReSize(nrows,ncols); // for type Matrix or nricMatrix
|
|
A.ReSize(n); // for all other types, except Band
|
|
A.ReSize(n,lower,upper); // for BandMatrix
|
|
A.ReSize(n,lower); // for LowerBandMatrix
|
|
A.ReSize(n,upper); // for UpperBandMatrix
|
|
A.ReSize(n,lower); // for SymmetricBandMatrix
|
|
A.ReSize(B); // set dims to those of B
|
|
</PRE>
|
|
|
|
<P>Use <TT>A.CleanUp()</TT> to set the dimensions of <TT>A</TT> to zero and
|
|
release all the heap memory. </P>
|
|
<P><TT>A.ReSize(B)</TT> sets the dimensions of <TT>A</TT> to those of a matrix
|
|
<TT>B</TT>. This includes the band-width in the case of a band matrix. It is an
|
|
error for <TT>A</TT> to be a band matrix and <TT>B</TT> not a band matrix (or
|
|
diagonal matrix). </P>
|
|
<P>Remember that <TT>ReSize</TT> destroys values. If you want to
|
|
<TT>ReSize</TT>, but keep the values in the bit that is left use something like
|
|
</P>
|
|
<PRE> ColumnVector V(100);
|
|
... // load values
|
|
V = V.Rows(1,50); // to get first 50 values.
|
|
</PRE>
|
|
|
|
<P>If you want to extend a matrix or vector use something like </P>
|
|
<PRE> ColumnVector V(50);
|
|
... // load values
|
|
{ V.Release(); ColumnVector X=V; V.ReSize(100); V.Rows(1,50)=X; }
|
|
// V now length 100
|
|
</PRE>
|
|
|
|
<H2><A NAME="ch_type"></A>3.13 Change type</H2>
|
|
<P CLASS="small"><A HREF="#solve">next</A> - <A HREF="#solve">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>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. </P>
|
|
<PRE> A.AsRow()
|
|
A.AsColumn()
|
|
A.AsDiagonal()
|
|
A.AsMatrix(nrows,ncols)
|
|
A.AsScalar()
|
|
</PRE>
|
|
|
|
<P>The expression <TT>A.AsScalar()</TT> is used to convert a 1 x 1 matrix to a
|
|
scalar. </P>
|
|
<H2><A NAME="solve"></A>3.14 Multiple matrix solve</H2>
|
|
<P CLASS="small"><A HREF="#memory">next</A> - <A HREF="#memory">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>To solve the matrix equation <TT>Ay = b</TT> where <TT>A</TT> is a square
|
|
matrix of equation coefficients, <TT>y</TT> is a column vector of values to be
|
|
solved for, and <TT>b</TT> is a column vector, use the code </P>
|
|
<PRE> int n = something
|
|
Matrix A(n,n); ColumnVector b(n);
|
|
... put values in A and b
|
|
ColumnVector y = A.i() * b; // solves matrix equation
|
|
</PRE>
|
|
|
|
<P>The following notes are for the case where you want to solve more than one
|
|
matrix equation with different values of <TT>b</TT> but the same <TT>A</TT>. Or
|
|
where you want to solve a matrix equation and also find the determinant of
|
|
<TT>A</TT>. In these cases you probably want to avoid repeating the LU
|
|
decomposition of <TT>A</TT> for each solve or determinant calculation. </P>
|
|
<P>If <TT>A</TT> is a square or symmetric matrix use </P>
|
|
<PRE> CroutMatrix X = A; // carries out LU decomposition
|
|
Matrix AP = X.i()*P; Matrix AQ = X.i()*Q;
|
|
LogAndSign ld = X.LogDeterminant();
|
|
</PRE>
|
|
|
|
<P>rather than </P>
|
|
<PRE> Matrix AP = A.i()*P; Matrix AQ = A.i()*Q;
|
|
LogAndSign ld = A.LogDeterminant();
|
|
</PRE>
|
|
|
|
<P>since each operation will repeat the LU decomposition. </P>
|
|
<P>If <TT>A</TT> is a BandMatrix or a SymmetricBandMatrix begin with </P>
|
|
<PRE> BandLUMatrix X = A; // carries out LU decomposition
|
|
</PRE>
|
|
|
|
<P>A CroutMatrix or a BandLUMatrix can't be manipulated or copied. Use
|
|
references as an alternative to copying. </P>
|
|
<P>Alternatively use </P>
|
|
<PRE> LinearEquationSolver X = A;
|
|
</PRE>
|
|
|
|
<P>This will choose the most appropriate decomposition of <TT>A</TT>. That is,
|
|
the band form if <TT>A</TT> is banded; the Crout decomposition if <TT>A</TT> is
|
|
square or symmetric and no decomposition if <TT>A</TT> is triangular or
|
|
diagonal.</P>
|
|
<H2><A NAME="memory"></A>3.15 Memory management</H2>
|
|
<P CLASS="small"><A HREF="#efficien">next</A> - <A HREF="#efficien">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>The package does not support delayed copy. Several strategies are required
|
|
to prevent unnecessary matrix copies. </P>
|
|
<P>Where a matrix is called as a function argument use a constant reference.
|
|
For example </P>
|
|
<PRE> YourFunction(const Matrix& A)
|
|
</PRE>
|
|
|
|
<P>rather than </P>
|
|
<PRE> YourFunction(Matrix A)
|
|
</PRE>
|
|
|
|
<P>Skip the rest of this section on your first reading. </P>
|
|
<P CLASS="small">Gnu g++ (< 2.6) users please read on: if you are returning
|
|
matrix values from a function, then you must use the ReturnMatrix
|
|
construct.</P>
|
|
<P>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. </P>
|
|
<P>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 </P>
|
|
<PRE> ReturnMatrix MakeAMatrix()
|
|
{
|
|
Matrix A; // or any other matrix type
|
|
......
|
|
A.Release(); return A;
|
|
}
|
|
</PRE>
|
|
|
|
<P>or </P>
|
|
<PRE> ReturnMatrix MakeAMatrix()
|
|
{
|
|
Matrix* m = new Matrix;
|
|
......
|
|
m->ReleaseAndDelete(); return *m;
|
|
}
|
|
</PRE>
|
|
|
|
<P>If your compiler objects to this code, replace the return statements with
|
|
</P>
|
|
<PRE> return A.ForReturn();
|
|
</PRE>
|
|
|
|
<P>or </P>
|
|
<PRE> return m->ForReturn();
|
|
</PRE>
|
|
|
|
<P>If you are using AT&T C++ you may wish to replace <TT>return A;</TT> by
|
|
return <TT>(ReturnMatrix)A;</TT> to avoid a warning message; but this will give
|
|
a runtime error with Gnu. (You can't please everyone.) </P>
|
|
<HR>
|
|
<P><B>Do not forget to make the function of type ReturnMatrix; otherwise you
|
|
may get incomprehensible run-time errors.</B> </P>
|
|
<HR>
|
|
<P>You can also use <TT>.Release()</TT> or <TT>->ReleaseAndDelete()</TT> to
|
|
allow a matrix expression to recycle space. Suppose you call </P>
|
|
<PRE> A.Release();
|
|
</PRE>
|
|
|
|
<P>just before <TT>A</TT> is used just once in an expression. Then the memory
|
|
used by <TT>A</TT> is either returned to the system or reused in the
|
|
expression. In either case, <TT>A</TT>'s memory is destroyed. This procedure
|
|
can be used to improve efficiency and reduce the use of memory. </P>
|
|
<P>Use <TT>->ReleaseAndDelete</TT> for matrices created by new if you want
|
|
to completely delete the matrix after it is accessed. </P>
|
|
<H2><A NAME="efficien"></A>3.16 Efficiency</H2>
|
|
<P CLASS="small"><A HREF="#output">next</A> - <A HREF="#output">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>The package tends to be not very efficient for dealing with matrices with
|
|
short rows. This is because some administration is required for accessing rows
|
|
for a variety of types of matrices. To reduce the administration a special
|
|
multiply routine is used for rectangular matrices in place of the generic one.
|
|
Where operations can be done without reference to the individual rows (such as
|
|
adding matrices of the same type) appropriate routines are used. </P>
|
|
<P>When you are using small matrices (say smaller than 10 x 10) you may find it
|
|
faster to use rectangular matrices rather than the triangular or symmetric
|
|
ones. </P>
|
|
<H2><A NAME="output"></A>3.17 Output</H2>
|
|
<P CLASS="small"><A HREF="#unspec">next</A> - <A HREF="#unspec">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>To print a matrix use an expression like </P>
|
|
<PRE> Matrix A;
|
|
......
|
|
cout << setw(10) << setprecision(5) << A;
|
|
</PRE>
|
|
|
|
<P>This will work only with systems that support the standard input/output
|
|
routines including manipulators. You need to #include the files iostream.h,
|
|
iomanip.h, newmatio.h in your C++ source files that use this facility. The
|
|
files iostream.h, iomanip.h will be included automatically if you include the
|
|
statement <TT>#define WANT_STREAM</TT> at the beginning of your source file. So
|
|
you can begin your file with either </P>
|
|
<PRE> #define WANT_STREAM
|
|
#include "newmatio.h"
|
|
</PRE>
|
|
|
|
<P>or </P>
|
|
<PRE> #include <iostream.h>
|
|
#include <iomanip.h>
|
|
#include "newmatio.h"
|
|
</PRE>
|
|
|
|
<P>The present version of this routine is useful only for matrices small enough
|
|
to fit within a page or screen width. </P>
|
|
<P>To print several vectors or matrices in columns use a <A
|
|
HREF="#binary">concatenation operator</A>: </P>
|
|
<PRE> ColumnVector A, B;
|
|
.....
|
|
cout << setw(10) << setprecision(5) << (A | B);
|
|
</PRE>
|
|
|
|
<H2><A NAME="unspec"></A>3.18 Unspecified type</H2>
|
|
<P CLASS="small"><A HREF="#cholesky">next</A> - <A HREF="#cholesky">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>Skip this section on your first reading. </P>
|
|
<P>If you want to work with a matrix of unknown type, say in a function. You
|
|
can construct a matrix of type <TT>GenericMatrix</TT>. Eg </P>
|
|
<PRE> Matrix A;
|
|
..... // put some values in A
|
|
GenericMatrix GM = A;
|
|
</PRE>
|
|
|
|
<P>A GenericMatrix matrix can be used anywhere where a matrix expression can be
|
|
used and also on the left hand side of an <TT>=</TT>. You can pass any type of
|
|
matrix (excluding the Crout and BandLUMatrix types) to a <TT>const
|
|
GenericMatrix&</TT> argument in a function. However most scalar functions
|
|
including Nrows(), Ncols(), Type() and element access do not work with it. Nor
|
|
does the ReturnMatrix construct. See also the paragraph on <A
|
|
HREF="#solve">LinearEquationSolver</A>. </P>
|
|
<P>An alternative and less flexible approach is to use BaseMatrix or
|
|
GeneralMatrix. </P>
|
|
<P>Suppose you wish to write a function which accesses a matrix of unknown type
|
|
including expressions (eg <TT>A*B</TT>). Then use a layout similar to the
|
|
following: </P>
|
|
<PRE> void YourFunction(BaseMatrix& X)
|
|
{
|
|
GeneralMatrix* gm = X.Evaluate(); // evaluate an expression
|
|
// if necessary
|
|
........ // operations on *gm
|
|
gm->tDelete(); // delete *gm if a temporary
|
|
}
|
|
</PRE>
|
|
|
|
<P>See, as an example, the definitions of <TT>operator<<</TT> in
|
|
newmat9.cpp. </P>
|
|
<P>Under certain circumstances; particularly where <TT>X</TT> is to be used
|
|
just once in an expression you can leave out the <TT>Evaluate()</TT> statement
|
|
and the corresponding <TT>tDelete()</TT>. Just use <TT>X</TT> in the
|
|
expression. </P>
|
|
<P>If you know YourFunction will never have to handle a formula as its argument
|
|
you could also use </P>
|
|
<PRE> void YourFunction(const GeneralMatrix& X)
|
|
{
|
|
........ // operations on X
|
|
}
|
|
</PRE>
|
|
|
|
<P>Do not try to construct a GeneralMatrix or BaseMatrix. </P>
|
|
<H2><A NAME="cholesky"></A>3.19 Cholesky
|
|
decomposition</H2>
|
|
<P CLASS="small"><A HREF="#qr">next</A> - <A HREF="#qr">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>Suppose <TT>S</TT> is symmetric and positive definite. Then there exists a
|
|
unique lower triangular matrix <TT>L</TT> such that <TT>L * L.t() = S</TT>. To
|
|
calculate this use </P>
|
|
<PRE> SymmetricMatrix S;
|
|
......
|
|
LowerTriangularMatrix L = Cholesky(S);
|
|
</PRE>
|
|
|
|
<P>If <TT>S</TT> is a symmetric band matrix then <TT>L</TT> is a band matrix
|
|
and an alternative procedure is provided for carrying out the decomposition:
|
|
</P>
|
|
<PRE> SymmetricBandMatrix S;
|
|
......
|
|
LowerBandMatrix L = Cholesky(S);
|
|
</PRE>
|
|
|
|
<H2><A NAME="qr"></A>3.20 QR decomposition</H2>
|
|
<P CLASS="small"><A HREF="#svd">next</A> - <A HREF="#svd">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>This is a variant on the usual QR transformation. </P>
|
|
<P>Start with matrix (dimensions shown to left and below the matrix)</P>
|
|
<PRE> / 0 0 \ s
|
|
\ X Y / n
|
|
|
|
s t
|
|
</PRE>
|
|
|
|
<P>Our version of the QR decomposition multiplies this matrix by an orthogonal
|
|
matrix Q to get </P>
|
|
<PRE> / U M \ s
|
|
\ 0 Z / n
|
|
|
|
s t
|
|
</PRE>
|
|
|
|
<P>where <TT>U</TT> is upper triangular (the R of the QR transform). That is
|
|
</P>
|
|
<PRE> Q / 0 0 \ = / U M \
|
|
\ X Y / \ 0 Z / </PRE>
|
|
|
|
<P>This is good for solving least squares problems: choose b (matrix or column
|
|
vector) to minimise the sum of the squares of the elements of </P>
|
|
<PRE> Y - X*b
|
|
</PRE>
|
|
|
|
<P>Then choose <TT>b = U.i()*M;</TT> The residuals <TT>Y - X*b</TT> are in
|
|
<TT>Z</TT>. </P>
|
|
<P>This is the usual QR transformation applied to the matrix <TT>X</TT> with
|
|
the square zero matrix concatenated on top of it. It gives the same triangular
|
|
matrix as the QR transform applied directly to <TT>X</TT> and generally seems
|
|
to work in the same way as the usual QR transform. However it fits into the
|
|
matrix package better and also gives us the residuals directly. It turns out to
|
|
be essentially a modified Gram-Schmidt decomposition. </P>
|
|
<P>Two routines are provided in <I>newmat</I>: </P>
|
|
<PRE> QRZ(X, U);
|
|
</PRE>
|
|
|
|
<P>replaces <TT>X</TT> by orthogonal columns and forms <TT>U</TT>. </P>
|
|
<PRE> QRZ(X, Y, M);
|
|
</PRE>
|
|
|
|
<P>uses <TT>X</TT> from the first routine, replaces <TT>Y</TT> by <TT>Z</TT>
|
|
and forms <TT>M</TT>. </P>
|
|
<P>The are also two routines <TT>QRZT(X, L)</TT> and <TT>QRZT(X, Y, M)</TT>
|
|
which do the same decomposition on the transposes of all these matrices. QRZT
|
|
replaces the routines HHDecompose in earlier versions of newmat. HHDecompose is
|
|
still defined but just calls QRZT.</P>
|
|
<P>For an example of the use of this decomposition see the file
|
|
<A HREF="#example">example.cpp</A>. </P>
|
|
<H2><A NAME="svd"></A>3.21 Singular value
|
|
decomposition</H2>
|
|
<P CLASS="small"><A HREF="#evalues">next</A> - <A HREF="#evalues">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>The singular value decomposition of an <I>m</I> x <I>n</I> <TT>Matrix</TT>
|
|
<TT>A</TT> (where <I>m</I> >= <I>n</I>) is a decomposition </P>
|
|
<PRE> A = U * D * V.t()
|
|
</PRE>
|
|
|
|
<P>where <TT>U</TT> is <I>m</I> x <I>n</I> with <TT>U.t() * U</TT> equalling
|
|
the identity, <TT>D</TT> is an <I>n</I> x <I>n </I><TT>DiagonalMatrix</TT> and
|
|
<TT>V</TT> is an <I>n</I> x <I>n</I> orthogonal matrix (type <TT>Matrix</TT> in
|
|
<I>Newmat</I>). </P>
|
|
<P>Singular value decompositions are useful for understanding the structure of
|
|
ill-conditioned matrices, solving least squares problems, and for finding the
|
|
eigenvalues of <TT>A.t() * A</TT>. </P>
|
|
<P>To calculate the singular value decomposition of <TT>A</TT> (with <I>m</I>
|
|
>= <I>n</I>) use one of </P>
|
|
<PRE> 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
|
|
</PRE>
|
|
|
|
<P>where <TT>A</TT>, <TT>U</TT> and <TT>V</TT> are of type <TT>Matrix</TT> and
|
|
<TT>D</TT> is a <TT>DiagonalMatrix</TT>. The values of <TT>A</TT> are not
|
|
changed unless <TT>A</TT> is also inserted as the third argument.</P>
|
|
|
|
<P>The elements of <tt>D</tt> are sorted in <i>descending</i> order.</P>
|
|
<P CLASS="small">Remember that the SVD decomposition is not completely unique. The signs of the elements in a column of <TT>U</TT> may be reversed
|
|
if the signs in the corresponding column in <TT>V</TT> are reversed. If a
|
|
number of the singular values are identical one can apply an orthogonal
|
|
transformation to the corresponding columns of <TT>U</TT> and the corresponding
|
|
columns of <TT>V</TT>.</P>
|
|
<H2><A NAME="evalues"></A>3.22 Eigenvalue
|
|
decomposition</H2>
|
|
<P CLASS="small"><A HREF="#sorting">next</A> - <A HREF="#sorting">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>An eigenvalue decomposition of a SymmetricMatrix <TT>A</TT> is a
|
|
decomposition </P>
|
|
<PRE> A = V * D * V.t()
|
|
</PRE>
|
|
|
|
<P>where <TT>V</TT> is an orthogonal matrix (type <TT>Matrix</TT> in
|
|
<I>Newmat</I>) and <TT>D</TT> is a DiagonalMatrix. </P>
|
|
<P>Eigenvalue analyses are used in a wide variety of engineering, statistical
|
|
and other mathematical analyses. </P>
|
|
<P>The package includes two algorithms: Jacobi and Householder. The first is
|
|
extremely reliable but much slower than the second. </P>
|
|
<P>The code is adapted from routines in <I>Handbook for Automatic Computation,
|
|
Vol II, Linear Algebra</I> by Wilkinson and Reinsch, published by Springer
|
|
Verlag. </P>
|
|
<PRE> Jacobi(A,D,S,V); // A, S symmetric; S is workspace,
|
|
// S = A is OK; V is a matrix
|
|
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; V is a matrix
|
|
|
|
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; V is a matrix
|
|
</PRE>
|
|
|
|
<P>where <TT>A</TT>, <TT>S</TT> are of type <TT>SymmetricMatrix</TT>,
|
|
<TT>D</TT> is of type <TT>DiagonalMatrix</TT> and <TT>V</TT> is of type
|
|
<TT>Matrix</TT>. The values of <TT>A</TT> are not changed unless <TT>A</TT> is
|
|
also inserted as the third argument. If you need eigenvectors use one of the
|
|
forms with matrix <TT>V</TT>. The eigenvectors are returned as the columns of
|
|
<TT>V</TT>.</P>
|
|
|
|
<P>The elements of <tt>D</tt> are sorted in <i>ascending</i> order.</P>
|
|
<P CLASS="small">Remember that an eigenvalue decomposition is not completely
|
|
unique - see the comments about the <A HREF="#svd">SVD</A>
|
|
decomposition. </P>
|
|
<H2><A NAME="sorting"></A>3.23 Sorting</H2>
|
|
<P CLASS="small"><A HREF="#fft">next</A> - <A HREF="#fft">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>To sort the values in a matrix or vector, <TT>A</TT>, (in general this
|
|
operation makes sense only for vectors and diagonal matrices) use one of</P>
|
|
<PRE> SortAscending(A);
|
|
|
|
SortDescending(A);
|
|
</PRE>
|
|
|
|
<P>I use the quicksort algorithm. The algorithm is similar to that in
|
|
Sedgewick's algorithms in C++. If the sort seems to be failing (as quicksort
|
|
can do) an exception is thrown. </P>
|
|
<P>You will get incorrect results if you try to sort a band matrix - but why
|
|
would you want to sort a band matrix? </P>
|
|
<H2><A NAME="fft"></A>3.24 Fast Fourier transform
|
|
</H2>
|
|
<P CLASS="small"><A HREF="#trigtran">next</A> - <A HREF="#trigtran">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<PRE> FFT(X, Y, F, G); // F=X and G=Y are OK
|
|
</PRE>
|
|
|
|
<P>where <TT>X</TT>, <TT>Y</TT>, <TT>F</TT>, <TT>G</TT> are column vectors.
|
|
<TT>X</TT> and <TT>Y</TT> are the real and imaginary input vectors; <TT>F</TT>
|
|
and <TT>G</TT> are the real and imaginary output vectors. The lengths of
|
|
<TT>X</TT> and <TT>Y</TT> must be equal and should be the product of numbers
|
|
less than about 10 for fast execution. </P>
|
|
<P>The formula is </P>
|
|
<PRE> n-1
|
|
h[k] = SUM z[j] exp (-2 pi i jk/n)
|
|
j=0
|
|
</PRE>
|
|
|
|
<P>where <TT>z[j]</TT> is stored complex and stored in <TT>X(j+1)</TT> and
|
|
<TT>Y(j+1)</TT>. Likewise <TT>h[k]</TT> is complex and stored in
|
|
<TT>F(k+1)</TT> and <TT>G(k+1)</TT>. The fast Fourier algorithm takes order <I>
|
|
n</I> log(<I>n</I>) operations (for <I>good</I> values of <I>n</I>) rather than
|
|
<I>n</I>**2 that straight evaluation (see the file <TT>tmtf.cpp</TT>) takes.
|
|
</P>
|
|
<P>I use one of two methods: </P>
|
|
<UL>
|
|
<LI>A program originally written by Sande and Gentleman. This requires that
|
|
<I>n</I> can be expressed as a product of small numbers.</LI>
|
|
<LI>A method of Carl de Boor (1980), <I>Siam J Sci Stat Comput</I>, pp 173-8.
|
|
The sines and cosines are calculated explicitly. This gives better accuracy, at
|
|
an expense of being a little slower than is otherwise possible. This is slower
|
|
than the Sande-Gentleman program but will work for all <I>n</I> --- although it
|
|
will be very slow for <I>bad</I> values of <I>n</I>.</LI>
|
|
</UL>
|
|
<P>Related functions </P>
|
|
<PRE> FFTI(F, G, X, Y); // X=F and Y=G are OK
|
|
RealFFT(X, F, G);
|
|
RealFFTI(F, G, X);
|
|
</PRE>
|
|
|
|
<P><TT>FFTI</TT> is the inverse transform for <TT>FFT</TT>. <TT>RealFFT</TT> is
|
|
for the case when the input vector is real, that is <TT>Y = 0</TT>. I assume
|
|
the length of <TT>X</TT>, denoted by <I>n</I>, is <I>even</I>. That is,
|
|
<I>n</I> must be divisible by 2. The program sets the lengths of <TT>F</TT> and
|
|
<TT>G</TT> to <I>n</I>/2 + 1. <TT>RealFFTI</TT> is the inverse of
|
|
<TT>RealFFT</TT>. </P>
|
|
<P>See also the section on fast <A HREF="#trigtran">trigonometric
|
|
transforms</A>. </P>
|
|
<H2><A NAME="trigtran"></A>3.25 Fast trigonometric
|
|
transforms</H2>
|
|
<P CLASS="small"><A HREF="#nric">next</A> - <A HREF="#nric">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>These are the sin and cosine transforms as defined by Charles Van Loan
|
|
(1992) in <I>Computational frameworks for the fast Fourier transform</I>
|
|
published by SIAM. See page 229. Some other authors use slightly different
|
|
conventions. All the functions call the <A HREF="#fft">fast Fourier
|
|
transforms</A> and require an <I> even</I> transform length, denoted by
|
|
<I>m</I> in these notes. That is, <I>m</I> must be divisible by 2. As with the
|
|
FFT <I>m</I> should be the product of numbers less than about 10 for fast
|
|
execution. </P>
|
|
<P>The functions I define are </P>
|
|
<PRE> DCT(U,V); // U, V are ColumnVectors, length <I>m+1</I>
|
|
DCT_inverse(V,U); // inverse of DCT
|
|
DST(U,V); // U, V are ColumnVectors, length <I>m+1</I>
|
|
DST_inverse(V,U); // inverse of DST
|
|
DCT_II(U,V); // U, V are ColumnVectors, length <I>m</I>
|
|
DCT_II_inverse(V,U); // inverse of DCT_II
|
|
DST_II(U,V); // U, V are ColumnVectors, length <I>m</I>
|
|
DST_II_inverse(V,U); // inverse of DST_II
|
|
</PRE>
|
|
|
|
<P>where the first argument is the input and the second argument is the output.
|
|
<TT>V = U</TT> is OK. The length of the output ColumnVector is set by the
|
|
functions. </P>
|
|
<P>Here are the formulae: </P>
|
|
<H3>DCT</H3>
|
|
<PRE> m-1 k
|
|
v[k] = u[0]/2 + SUM { u[j] cos (pi jk/m) } + (-) u[m]/2
|
|
j=1
|
|
</PRE>
|
|
|
|
<P>for <TT>k = 0...m</TT>, where <TT>u[j]</TT> and <TT>v[k]</TT> are stored in
|
|
<TT>U(j+1)</TT> and <TT>V(k+1)</TT>. </P>
|
|
<H3>DST</H3>
|
|
<PRE> m-1
|
|
v[k] = SUM { u[j] sin (pi jk/m) }
|
|
j=1
|
|
</PRE>
|
|
|
|
<P>for <TT>k = 1...(m-1)</TT>, where <TT>u[j]</TT> and <TT>v[k]</TT> are stored
|
|
in <TT>U(j+1)</TT> and <TT>V(k+1)</TT>and where <TT>u[0]</TT> and <TT>u[m]</TT>
|
|
are ignored and <TT>v[0]</TT> and <TT>v[m]</TT> are set to zero. For the
|
|
inverse function <TT>v[0]</TT> and <TT>v[m]</TT> are ignored and <TT>u[0]</TT>
|
|
and <TT>u[m]</TT> are set to zero. </P>
|
|
<H3>DCT_II</H3>
|
|
<PRE> m-1
|
|
v[k] = SUM { u[j] cos (pi (j+1/2)k/m) }
|
|
j=0
|
|
</PRE>
|
|
|
|
<P>for <TT>k = 0...(m-1)</TT>, where <TT>u[j]</TT> and <TT>v[k]</TT> are stored
|
|
in <TT>U(j+1)</TT> and <TT>V(k+1)</TT>. </P>
|
|
<H3>DST_II</H3>
|
|
<PRE> m
|
|
v[k] = SUM { u[j] sin (pi (j-1/2)k/m) }
|
|
j=1
|
|
</PRE>
|
|
|
|
<P>for <TT>k = 1...m</TT>, where <TT>u[j]</TT> and <TT>v[k]</TT> are stored in
|
|
<TT>U(j)</TT> and <TT>V(k)</TT>. </P>
|
|
<P>Note that the relationship between the subscripts in the formulae and those
|
|
used in <I>newmat</I> is different for DST_II (and DST_II_inverse). </P>
|
|
<H2><A NAME="nric"></A>3.26 Interface to Numerical
|
|
Recipes in C</H2>
|
|
<P CLASS="small"><A HREF="#except">next</A> - <A HREF="#except">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>This package can be used with the vectors and matrices defined in
|
|
<I>Numerical Recipes in C</I>. 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 may
|
|
need to edit the function definitions to use the version acceptable to your
|
|
compiler (if you are using the first edition of NRIC). You may need to enclose
|
|
the code from Numerical Recipes in <TT>extern "C" { ... }</TT>. You
|
|
will also need to include the matrix and vector utility routines. </P>
|
|
<P>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
|
|
<TT>X</TT> in an function from Numerical Recipes, use <TT>X.nric()</TT> in the
|
|
function call. </P>
|
|
<P>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. </P>
|
|
<P>For example </P>
|
|
<PRE> 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
|
|
</PRE>
|
|
|
|
<H2><A NAME="except"></A>3.27 Exceptions</H2>
|
|
<P CLASS="small"><A HREF="#cleanup">next</A> - <A HREF="#cleanup">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>Here is the class structure for exceptions: </P>
|
|
<PRE>Exception
|
|
Logic_error
|
|
ProgramException miscellaneous matrix error
|
|
IndexException index out of bounds
|
|
VectorException unable to convert matrix to vector
|
|
NotSquareException matrix is not square (invert, solve)
|
|
SubMatrixDimensionException out of bounds index of submatrix
|
|
IncompatibleDimensionsException (multiply, add etc)
|
|
NotDefinedException operation not defined (eg <)
|
|
CannotBuildException copying a matrix where copy is undefined
|
|
InternalException probably an error in newmat
|
|
Runtime_error
|
|
NPDException matrix not positive definite (Cholesky)
|
|
ConvergenceException no convergence (e-values, non-linear, sort)
|
|
SingularException matrix is singular (invert, solve)
|
|
SolutionException no convergence in solution routine
|
|
OverflowException floating point overflow
|
|
Bad_alloc out of space (new fails)
|
|
</PRE>
|
|
|
|
<P>I have attempted to mimic the exception class structure in the C++ standard
|
|
library, by defining the Logic_error and Runtime_error classes. </P>
|
|
<P>Suppose you have edited <tt>include.h</tt> to use my <i>simulated</i>
|
|
exceptions or to <i>disable</i> exceptions. If there is no catch statement or exceptions are disabled then my
|
|
<TT>Terminate()</TT> function in <TT>myexcept.h</TT> is called when you throw an
|
|
exception. This prints out
|
|
an error message, the dimensions and types of the matrices involved, the name
|
|
of the routine detecting the exception, and any other information set by the
|
|
<A HREF="#error">Tracer</A> class. Also see the section on <A
|
|
HREF="#error">error messages</A> for additional notes on the messages generated
|
|
by the exceptions. </P>
|
|
<P>You can also print this information in a <i>catch</i> clause by printing <TT>Exception::what()</TT>.
|
|
</P>
|
|
<P>If you are using <i> compiler supported</i> exceptions then see the section
|
|
on <a href="#except_1">catching exceptions</a>.
|
|
</P>
|
|
<P>See the file <TT>test_exc.cpp</TT> as an example of catching an exception
|
|
and printing the error message. </P>
|
|
<P>The 08 version of newmat defined a member function <TT>void
|
|
SetAction(int)</TT> to help customise the action when an exception is called.
|
|
This has been deleted in the 09 and 10 versions. Now include an instruction
|
|
such as <TT>cout << Exception::what() << endl;</TT> in the
|
|
<TT>Catch</TT> or <TT>CatchAll</TT> block to determine the action. </P>
|
|
<P>The library includes the alternatives of using the inbuilt exceptions
|
|
provided by a compiler, simulating exceptions, or disabling exceptions. See
|
|
<A HREF="#custom">customising</A> for selecting the correct exception option.
|
|
</P>
|
|
<P>The rest of this section describes my partial simulation of exceptions for
|
|
compilers which do not support C++ exceptions. I use Carlos Vidal's article in
|
|
the September 1992 <I>C Users Journal</I> as a starting point. </P>
|
|
<P>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
|
|
major problem, but it is something that needs to be sorted out. </P>
|
|
<P>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. </P>
|
|
<P>I have added another macro ReThrow to mean a rethrow, Throw(). This was
|
|
necessary to enable the package to be compatible with both my exception package
|
|
and C++ exceptions. </P>
|
|
<P>If you want to throw an exception, use a statement like </P>
|
|
<PRE> Throw(Exception("Error message\n"));
|
|
</PRE>
|
|
|
|
<P>It is important to have the exception declaration in the Throw statement,
|
|
rather than as a separate statement. </P>
|
|
<P>All exception classes must be derived from the class, Exception, defined in
|
|
newmat and can contain only static variables. See the examples in newmat if you
|
|
want to define additional exceptions.</P>
|
|
<P>Note that the simulation exception mechanism does not work if you define
|
|
arrays of matrices. </P>
|
|
<H2><A NAME="cleanup"></A>3.28 Cleanup after an
|
|
exception</H2>
|
|
<P CLASS="small"><A HREF="#nonlin">next</A> - <A HREF="#nonlin">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>This section is about the simulated exceptions used in newmat. It is
|
|
irrelevant if you are using the exceptions built into a compiler or have set
|
|
the disable-exceptions option. </P>
|
|
<P>The simulated 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 <I>new</I> as
|
|
heap memory). For example, when you call </P>
|
|
<PRE> Matrix A(20,30);
|
|
</PRE>
|
|
|
|
<P>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. </P>
|
|
<P>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. </P>
|
|
<P>However it will not return heap memory from objects from other packages.
|
|
</P>
|
|
<P>If you want the mechanism to work with another class you will have to do
|
|
four things: </P>
|
|
<OL>
|
|
<LI>derive your class from class Janitor defined in except.h; </LI>
|
|
<LI>define a function <TT>void CleanUp()</TT> in that class to return all heap
|
|
memory; </LI>
|
|
<LI>include the following lines in the class definition <PRE> public:
|
|
void* operator new(size_t size)
|
|
{ do_not_link=true; void* t = ::operator new(size); return t; }
|
|
void operator delete(void* t) { ::operator delete(t); }
|
|
</PRE>
|
|
|
|
</LI>
|
|
<LI>be sure to include a copy constructor in you class definition, that is,
|
|
something like <PRE> X(const X&);
|
|
</PRE>
|
|
|
|
</LI>
|
|
</OL>
|
|
<P>Note that the function <TT>CleanUp()</TT> does somewhat the same duties as
|
|
the destructor. However <TT>CleanUp()</TT> has to do the <I>cleaning</I> for
|
|
the class you are working with and also the classes it is derived from. So it
|
|
will often be wrong to use exactly the same code for both <TT>CleanUp()</TT>
|
|
and the destructor or to define your destructor as a call to
|
|
<TT>CleanUp()</TT>. </P>
|
|
<H2><A NAME="nonlin"></A>3.29 Non-linear
|
|
applications</H2>
|
|
<P CLASS="small"><A HREF="#stl">next</A> - <A HREF="#stl">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>Files solution.h, solution.cpp contain a class for solving for <I>x</I> in
|
|
<I>y</I> = <I>f</I>(<I>x)</I> where <I>x</I> is a one-dimensional continuous
|
|
monotonic function. This is not a matrix thing at all but is included because
|
|
it is a useful thing and because it is a simpler version of the technique used
|
|
in the non-linear least squares. </P>
|
|
<P>Files newmatnl.h, newmatnl.cpp contain a series of classes for non-linear
|
|
least squares and maximum likelihood. These classes work on very well-behaved
|
|
functions but need upgrading for less well-behaved functions.</P>
|
|
<P>Documentation for both of these is in the definition files. Simple examples
|
|
are in sl_ex.cpp, nl_ex.cpp and garch.cpp.</P>
|
|
<H2><A NAME="stl"></A>3.30 Standard template
|
|
library</H2>
|
|
<P CLASS="small"><A HREF="#namesp">next</A> - <A HREF="#namesp">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P>The standard template library (STL) is the set of <I>container templates</I>
|
|
(vector, deque, list etc) defined by the C++ standards committee. Newmat is
|
|
intended to be compatible with the STL in the sense that you can store matrices
|
|
in the standard containers. I have defined <A HREF="#binary"><TT>==</TT> and
|
|
inequality </A>operators which seem to be required by some versions of the STL. Probably there will have
|
|
to be some other changes. My experiments with the Rogue Wave STL that comes
|
|
with Borland C++ 5.0 showed that some things worked and some things
|
|
unexpectedly didn't work. </P>
|
|
<P>If you want to use the container classes with Newmat please note </P>
|
|
<UL>
|
|
<LI>Don't use simulated exceptions. </LI>
|
|
<LI>Make sure the option <A HREF="#custom">DO_FREE_CHECK</A> is <EM>not</EM>
|
|
turned on. </LI>
|
|
<LI>You can store only one type of matrix in a container. If you want to use a
|
|
variety of types use the GenericMatrix type or store pointers to the matrices.
|
|
</LI>
|
|
<LI>The vector and deque container templates like to copy their elements. For
|
|
the vector container this happens when you insert an element anywhere except at
|
|
the end or when you append an element and the current vector storage overflows.
|
|
Since Newmat does not have <I>copy-on-write</I> this could get very
|
|
inefficient. (Later versions may have <I>copy-on-write</I> for the
|
|
GenericMatrix type). </LI>
|
|
<LI>You won't be able to sort the container or do anything that would call an
|
|
inequality operator. </LI>
|
|
</UL>
|
|
<P>I doubt whether the STL container will be used often for matrices. So I
|
|
don't think these limitations are very critical. If you think otherwise, please
|
|
tell me. </P>
|
|
<H2><A NAME="namesp"></A>3.31 Namespace</H2>
|
|
<P CLASS="small"><A HREF="#error">next</A> - <A HREF="#error">skip</A> -
|
|
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
|
|
<P><I>Namespace</I> is a new facility in C++. Its purpose is to avoid name
|
|
clashes between different libraries. I have included the namespace capability.
|
|
Activate the line <TT>#define use_namespace</TT> in <TT>include.h</TT>. Then
|
|
include either the statement </P>
|
|
<PRE> using namespace NEWMAT;
|
|
</PRE>
|
|
|
|
<P>at the beginning of any file that needs to access the newmat library or </P>
|
|
<PRE> using namespace RBD_LIBRARIES;
|
|
</PRE>
|
|
|
|
<P>at the beginning of any file that needs to access all my libraries. </P>
|
|
<P>This works correctly with <A HREF="#borland">Borland</A> C++ version 5. </P>
|
|
<P><A HREF="#microso">Microsoft</A> Visual C++ version 5 works in my example
|
|
and test files, but fails with apparently insignificant changes (it may be more
|
|
reliable if you have applied service pack 3). If you #include
|
|
"newmatap.h", but no other newmat include file, then also #include
|
|
"newmatio.h". It seems to work with <A HREF="#microso">Microsoft</A>
|
|
Visual C++ version 6 <I>if</I> you have applied at least service pack 2.</P>
|
|
<P>My use of namespace works with <A HREF="#gcc">Gnu g++</A> version
|
|
3.3 </P>
|
|
<P>I have defined the following namespaces: </P>
|
|
<UL>
|
|
<LI>RBD_COMMON for functions and classes used by most of my libraries </LI>
|
|
<LI>NEWMAT for the newmat library </LI>
|
|
<LI>RBD_LIBRARIES for all my libraries </LI>
|
|
</UL>
|
|
<H2><A NAME="error"></A>4. Error messages</H2>
|
|
<P CLASS="small"><A HREF="#design">next</A> - <A HREF="#design">skip</A> -
|
|
<A HREF="#top">up</A> - <A HREF="#top">start</A></P>
|
|
<P>Most error messages are self-explanatory. The message gives the size of the
|
|
matrices involved. Matrix types are referred to by the following codes: </P>
|
|
<PRE> Matrix or vector Rect
|
|
Symmetric matrix Sym
|
|
Band matrix Band
|
|
Symmetric band matrix SmBnd
|
|
Lower triangular matrix LT
|
|
Lower triangular band matrix LwBnd
|
|
Upper triangular matrix UT
|
|
Upper triangular band matrix UpBnd
|
|
Diagonal matrix Diag
|
|
Crout matrix (LU matrix) Crout
|
|
Band LU matrix BndLU
|
|
</PRE>
|
|
|
|
<P>Other codes should not occur. </P>
|
|
<P>See the section on <A HREF="#except">exceptions</A> for more details on the
|
|
structure of the exception classes. </P>
|
|
<P>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 </P>
|
|
<PRE> Tracer tr("name");
|
|
</PRE>
|
|
|
|
<P>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 </P>
|
|
<PRE> tr.ReName("new name");
|
|
</PRE>
|
|
|
|
<P>if, for example, you want to track progress through the function. </P>
|
|
<H2><A NAME="design"></A>5. Notes on the design of the
|
|
library</H2>
|
|
<P CLASS="small"><A HREF="#sue">next</A> - <A HREF="#top">skip</A> -
|
|
<A HREF="#top">up</A> - <A HREF="#top">start</A></P>
|
|
<TABLE WIDTH="100%">
|
|
<TR>
|
|
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A HREF="#sue">5.1
|
|
Safety, usability, efficiency </A><BR>
|
|
<A HREF="#mat_arr">5.2 Matrix vs array library </A><BR>
|
|
<A HREF="#question">5.3 Design questions </A><BR>
|
|
<A HREF="#stor">5.4 Data storage </A><BR>
|
|
<A HREF="#mem_man">5.5 Memory management - 1 </A><BR>
|
|
<A HREF="#mem_man2">5.6 Memory management - 2 </A><BR>
|
|
<A HREF="#evalx">5.7 Evaluation of expressions </A></TD>
|
|
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A
|
|
HREF="#explode">5.8 Explosion in the number of operations </A><BR>
|
|
<A HREF="#destr">5.9 Destruction of temporaries </A><BR>
|
|
<A HREF="#calc">5.10 A calculus of matrix types </A><BR>
|
|
<A HREF="#pointer">5.11 Pointer arithmetic </A><BR>
|
|
<A HREF="#err_hand">5.12 Error handling </A><BR>
|
|
<A HREF="#sparse">5.13 Sparse matrices </A><BR>
|
|
<A HREF="#comp_mat">5.14 Complex matrices </A></TD>
|
|
</TR>
|
|
</TABLE>
|
|
<P>I describe some of the ideas behind this package, some of the decisions that
|
|
I needed to make and give some details about the way it works. You don't need
|
|
to read this part of the documentation in order to use the package. </P>
|
|
<P>It isn't obvious what is the best way of going about structuring a matrix
|
|
package. I don't think you can figure this out with <I>thought</I> 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 <I>ultimate</I> package, but simply a trial of a number of ideas.
|
|
The following pages give some background on these ideas. </P>
|
|
<H2><A NAME="sue"></A>5.1 Safety, usability,
|
|
efficiency</H2>
|
|
<P CLASS="small"><A HREF="#mat_arr">next</A> - <A HREF="#mat_arr">skip</A> -
|
|
<A HREF="#design">up</A> - <A HREF="#top">start</A></P>
|
|
<H3>Some general comments</H3>
|
|
<P>A library like <I>newmat</I> needs to balance <I>safety</I>,
|
|
<I>usability</I> and <I>efficiency</I>. </P>
|
|
<P>By <B>safety</B>, I mean getting the right answer, and not causing crashes
|
|
or damage to the computer system. </P>
|
|
<P>By <B>usability</B>, I mean being easy to learn and use, including not being
|
|
too complicated, being intuitive, saving the users' time, being nice to use.
|
|
</P>
|
|
<P><B>Efficiency</B> means minimising the use of computer memory and time. </P>
|
|
<P>In the early days of computers the emphasis was on efficiency. But computer
|
|
power gets cheaper and cheaper, halving in price every 18 months. On the other
|
|
hand the unaided human brain is probably not a lot better than it was 100,000
|
|
years ago! So we should expect the balance to shift to put more emphasis on
|
|
safety and usability and a little less on efficiency. So I don't mind if my
|
|
programs are a little less efficient than programs written in pure C (or
|
|
Fortran) if I gain substantially in safety and usability. But I would mind if
|
|
they were a lot less efficient. </P>
|
|
<H3>Type of use</H3>
|
|
<P>Second reason for putting extra emphasis on safety and usability is the way
|
|
I and, I suspect, most other users actually use <I>newmat</I>. Most completed
|
|
programs are used only a few times. Some result is required for a client, paper
|
|
or thesis. The program is developed and tested, the result is obtained, and the
|
|
program archived. Of course bits of the program will be recycled for the next
|
|
project. But it may be less usual for the same program to be run over and over
|
|
again. So the cost, computer time + people time, is in the development time and
|
|
often, much less in the actual time to run the final program. So good use of
|
|
people time, especially during development is really important. This means you
|
|
need highly usable libraries. </P>
|
|
<P>So if you are dealing with matrices, you want the good interface that I have
|
|
tried to provide in <I>newmat</I>, and, of course, reliable methods underneath
|
|
it. </P>
|
|
<P>Of course, efficiency is still important. We often want to run the biggest
|
|
problem our computer will handle and often a little bigger. The C++ language
|
|
almost lets us have both worlds. We can define a reasonably good interface, and
|
|
get good efficiency in the use of the computer. </P>
|
|
<H3>Levels of access</H3>
|
|
<P>We can imagine the <I>black box</I> model of a <I>newmat</I>. Suppose the
|
|
inside is hidden but can be accessed by the methods described in the
|
|
<A HREF="#refer">reference</A> section. Then the interface is reasonably
|
|
consistent and intuitive. Matrices can be accessed and manipulated in much the
|
|
same way as doubles or ints in regular C. All accesses are checked. It is most
|
|
unlikely that an incorrect index will crash the system. In general, users do
|
|
not need to use pointers, so one shouldn't get pointers pointing into space.
|
|
And, hopefully, you will get simpler code and so less errors. </P>
|
|
<P>There are some exceptions to this. In particular, the <A
|
|
HREF="#elements">C-like subscripts</A> are not checked for validity. They give
|
|
faster access but with a lower level of safety. </P>
|
|
<P>Then there is the <A HREF="#scalar1">Store()</A> function which takes you to
|
|
the data array within a matrix. This takes you right inside the <I>black
|
|
box</I>. But this is what you have to use if you are writing, for example, a
|
|
new matrix factorisation, and require fast access to the data array. I have
|
|
tried to write code to simplify access to the interior of a rectangular matrix,
|
|
see file newmatrm.cpp, but I don't regard this as very successful, as yet, and
|
|
have not included it in the documentation. Ideally we should have improved
|
|
versions of this code for each of the major types of matrix. But, in reality,
|
|
most of my matrix factorisations are written in what is basically the C
|
|
language with very little C++. </P>
|
|
<P>So our <I>box</I> is not very <I>black</I>. You have a choice of how far you
|
|
penetrate. On the outside you have a good level of safety, but in some cases
|
|
efficiency is compromised a little. If you penetrate inside the <I>box</I>
|
|
safety is reduced but you can get better efficiency. </P>
|
|
<H3>Some performance data</H3>
|
|
<P>This section looks at the performance on <I>newmat</I> for simple sums,
|
|
comparing it with C code and with a simple array program.
|
|
</P>
|
|
<P>The following table lists the time (in seconds) for carrying out the
|
|
operations <TT>X=A+B;</TT>, <TT>X=A+B+C;</TT>, <TT>X=A+B+C+D;</TT>, <TT>X=A+B+C+D+E;</TT> where
|
|
<TT>X,A,B,C,D,E</TT> are of type ColumnVector, with a variety of programs. I am
|
|
using Microsoft VC++, version 6 in console mode under windows 2000 on a PC with
|
|
a 1 ghz Pentium III and 512 mbytes of memory. </P>
|
|
<pre> length iters. newmat C C-res. subs. array
|
|
<b>X = A + B</b>
|
|
2 5000000 27.8 0.3 8.8 1.9 9.5
|
|
20 500000 3.0 0.3 1.1 1.9 1.2
|
|
200 50000 0.5 0.3 0.4 1.9 0.3
|
|
2000 5000 0.4 0.3 0.4 2.0 1.0
|
|
20000 500 4.5 4.5 4.5 6.7 4.4
|
|
200000 50 5.2 4.7 5.5 5.8 5.2
|
|
|
|
<b>X = A + B + C</b>
|
|
2 5000000 36.6 0.4 8.9 2.5 12.2
|
|
20 500000 4.0 0.4 1.2 2.5 1.6
|
|
200 50000 0.8 0.3 0.5 2.5 0.5
|
|
2000 5000 3.6 4.4 4.6 9.0 4.4
|
|
20000 500 6.8 5.4 5.4 9.6 6.8
|
|
200000 50 8.6 6.0 6.7 7.1 8.6
|
|
|
|
<b>X = A + B + C + D</b>
|
|
2 5000000 44.0 0.7 9.3 3.1 14.6
|
|
20 500000 4.9 0.6 1.5 3.1 1.9
|
|
200 50000 1.0 0.6 0.8 3.2 0.8
|
|
2000 5000 5.6 6.6 6.8 11.5 5.9
|
|
20000 500 9.0 6.7 6.8 11.0 8.5
|
|
200000 50 11.9 7.1 7.9 9.5 12.0
|
|
|
|
<b>X = A + B + C + D + E</b>
|
|
2 5000000 50.6 1.0 9.5 3.8 17.1
|
|
20 500000 5.7 0.8 1.7 3.9 2.4
|
|
200 50000 1.3 0.9 1.0 3.9 1.0
|
|
2000 5000 7.0 8.3 8.2 13.8 7.1
|
|
20000 500 11.5 8.1 8.4 13.2 11.0
|
|
200000 50 15.2 8.7 9.5 12.4 15.4
|
|
</pre>
|
|
<P>I have <a target="_blank" href="add_time.png">graphed</a> the results
|
|
and included rather more array lengths.</P>
|
|
|
|
<P>The first column gives the lengths of the arrays, the second the number of
|
|
iterations and the remaining columns the total time required in seconds. If the
|
|
only thing that consumed time was the double precision addition then the
|
|
numbers within each block of the table would be the same. The summation is
|
|
repeated 5 times within each loop, for example:</P>
|
|
<pre> for (i=1; i<=m; ++i)
|
|
{
|
|
X1 = A1+B1+C1; X2 = A2+B2+C2; X3 = A3+B3+C3;
|
|
X4 = A4+B4+C4; X5 = A5+B5+C5;
|
|
}</pre>
|
|
<P>The column labelled <I>newmat</I> is using the standard <I>newmat</I> add. The column labelled <I>C</I> uses the usual C method: <TT>while
|
|
(j1--) *x1++ = *a1++ + *b1++;</TT> . The following column also includes an
|
|
<TT>X.ReSize()</TT> in the outer loop to correspond to the reassignment of
|
|
memory that <I>newmat</I> would do. In the next column the calculation is using
|
|
the usual C style <I>for</I> loop
|
|
and accessing the elements using <I>newmat</I> subscripts such as
|
|
<TT>A(i)</TT>. The final column is the time taken by a
|
|
simple array package. This uses an alternative method for avoiding temporaries
|
|
and unnecessary copies that does not involve runtime tests. It does its sums in blocks of 4 and copies in blocks of
|
|
8 in the same way that <I>newmat</I> does. </P>
|
|
<P>Here are my conclusions. </P>
|
|
<UL>
|
|
<LI><I>Newmat</I> does very badly for length 2 and doesn't do well for
|
|
length 20. There is a lot of code in <I>newmat</I> for
|
|
determining which sum algorithm to use and it is not surprising that this
|
|
impacts on performance for small lengths.
|
|
However the <I>array</I> program is also having difficulty with length 2 so it
|
|
is unlikely that the problem could be completely eliminated. </LI>
|
|
<LI>For arrays of length 2000 or longer <I>newmat</I> is doing about as well as
|
|
C and slightly better than C with resize in the <TT>X=A+B</TT> table. For the
|
|
other two tables it tends to be slower, but not dramatically so. </LI>
|
|
<LI>It is really important for fast processing with the Pentium III to stay
|
|
within the Pentium cache.</LI>
|
|
<LI>Addition using the <I>newmat</I> subscripts, while considerably slower than
|
|
the others, is still surprisingly good for the longer arrays. </LI>
|
|
<LI>The <I>array</I> program and <I>newmat</I> are similar for
|
|
lengths 2000 or higher (the longer times for the array program for the longest
|
|
arrays shown on the graph are probably a quirk of the timing program). </LI>
|
|
</UL>
|
|
<P>In summary: for the situation considered here, <I>newmat</I> is doing very
|
|
well for large ColumnVectors, even for sums with several terms, but not so well
|
|
for shorter ColumnVectors. </P>
|
|
<H2><A NAME="mat_arr"></A>5.2 Matrix vs array
|
|
library</H2>
|
|
<P CLASS="small"><A HREF="#question">next</A> - <A HREF="#question">skip</A> -
|
|
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
|
|
<P>The <I>newmat</I> library is for the manipulation of matrices, including the
|
|
standard operations such as multiplication as understood by numerical analysts,
|
|
engineers and mathematicians. </P>
|
|
<P>A matrix is a two dimensional array of numbers. However, very special
|
|
operations such as matrix multiplication are defined specifically for matrices.
|
|
This means that a <I>matrix</I> library, as I understand the term, is different
|
|
from a general <I>array</I> library. Here are some contrasting properties.</P>
|
|
<TABLE WIDTH="100%" BORDER="1">
|
|
<TR>
|
|
<TH VALIGN="TOP" WIDTH="20%" ALIGN="LEFT">Feature</TH>
|
|
<TH VALIGN="TOP" WIDTH="40%" ALIGN="LEFT">Matrix
|
|
library</TH>
|
|
<TH VALIGN="TOP" WIDTH="40%" ALIGN="LEFT">Array
|
|
library</TH>
|
|
</TR>
|
|
<TR>
|
|
<TD VALIGN="TOP" WIDTH="20%">Expressions</TD>
|
|
<TD VALIGN="TOP" WIDTH="40%">Matrix expressions<TT>;</TT>
|
|
<TT>*</TT> means matrix multiply; inverse function</TD>
|
|
<TD VALIGN="TOP" WIDTH="40%">Arithmetic operations, if
|
|
supported, mean elementwise combination of arrays</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD VALIGN="TOP" WIDTH="20%">Element access</TD>
|
|
<TD VALIGN="TOP" WIDTH="40%">Access to the elements of a
|
|
matrix</TD>
|
|
<TD VALIGN="TOP" WIDTH="40%">High speed access to elements
|
|
directly and perhaps with iterators</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD VALIGN="TOP" WIDTH="20%">Elementary functions</TD>
|
|
<TD VALIGN="TOP" WIDTH="40%">For example: determinant,
|
|
trace</TD>
|
|
<TD VALIGN="TOP" WIDTH="40%">Matrix multiplication as a
|
|
function</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD VALIGN="TOP" WIDTH="20%">Advanced functions</TD>
|
|
<TD VALIGN="TOP" WIDTH="40%">For example: eigenvalue
|
|
analysis</TD>
|
|
<TD VALIGN="TOP" WIDTH="40%"> </TD>
|
|
</TR>
|
|
<TR>
|
|
<TD VALIGN="TOP" WIDTH="20%">Element types</TD>
|
|
<TD VALIGN="TOP" WIDTH="40%">Real and possibly
|
|
complex</TD>
|
|
<TD VALIGN="TOP" WIDTH="40%">Wide range - real, integer, string
|
|
etc</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD VALIGN="TOP" WIDTH="20%">Types</TD>
|
|
<TD VALIGN="TOP" WIDTH="40%">Rectangular, symmetric, diagonal,
|
|
etc</TD>
|
|
<TD VALIGN="TOP" WIDTH="40%">One, two and three dimensional
|
|
arrays, at least</TD>
|
|
</TR>
|
|
</TABLE>
|
|
<P>Both types of library need to support access to sub-matrices or sub-arrays,
|
|
have good efficiency and storage management, and graceful exit for errors. In
|
|
both cases, we probably need two versions, one optimised for large matrices or
|
|
arrays and one for small matrices or arrays.</P>
|
|
<P>It may be possible to amalgamate the two sets of requirements to some
|
|
extent. However <I>newmat</I> is definitely oriented towards the matrix library
|
|
set.</P>
|
|
<H2><A NAME="question"></A>5.3 Design questions</H2>
|
|
<P CLASS="small"><A HREF="#stor">next</A> - <A HREF="#stor">skip</A> -
|
|
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
|
|
<P>Even within the bounds set by the requirements of a matrix library there is
|
|
a substantial opportunity for variation between what different matrix packages
|
|
might provide. It is not possible to 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 in 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. This package is intended to be one of these packages;
|
|
but not all of them. </P>
|
|
<P>Since my background is in statistical methods, this package is oriented
|
|
towards the kinds things you need for statistical analyses. </P>
|
|
<P>Now looking at some specific questions. </P>
|
|
<H3>What size of matrices?</H3>
|
|
<P>A matrix library may target small matrices (say 3 x 3), or medium sized
|
|
matrices, or very large matrices. </P>
|
|
<P>A library targeting very small matrices will seek to minimise
|
|
administration. A library for medium sized or very large matrices can spend
|
|
more time on administration in order to conserve space or optimise the
|
|
evaluation of expressions. A library for very large matrices will need to pay
|
|
special attention to storage and numerical properties. This library is designed
|
|
for medium sized matrices. This means it is worth introducing some
|
|
optimisations, but I don't have to worry about setting up some form of virtual
|
|
memory. </P>
|
|
<H3>Which matrix types?</H3>
|
|
<P>As well as the usual rectangular matrices, matrices occurring repeatedly in
|
|
numerical calculations are upper and lower triangular matrices, symmetric
|
|
matrices and diagonal matrices. This is particularly the case in calculations
|
|
involving least squares and eigenvalue calculations. So as a first stage these
|
|
were the types I decided to include. </P>
|
|
<P>It is also necessary to have types row vector and column vector. In a
|
|
<I>matrix</I> package, in contrast to an <I>array</I> package, it is necessary
|
|
to have both these types since they behave differently in matrix expressions.
|
|
The vector types can be derived for the rectangular matrix type, so having them
|
|
does not greatly increase the complexity of the package. </P>
|
|
<P>The problem with having several matrix types is the number of versions of
|
|
the binary operators one needs. If one has 5 distinct matrix types then a
|
|
simple library will need 25 versions of each of the binary operators. In fact,
|
|
we can evade this problem, but at the cost of some complexity. </P>
|
|
<H3>What element types?</H3>
|
|
<P>Ideally we would allow element types double, float, complex and int, at
|
|
least. It might be reasonably easy, using templates or equivalent, to provide a
|
|
library which could handle a variety of element types. However, as soon as one
|
|
starts implementing the binary operators between matrices with different
|
|
element types, again one gets an explosion in the number of operations one
|
|
needs to consider. At the present time the compilers I deal with are not up to
|
|
handling this problem with templates. (Of course, when I started writing
|
|
<I>newmat</I> there were no templates). But even when the compilers do meet the
|
|
specifications of the draft standard, writing a matrix package that allows for
|
|
a variety of element types using the template mechanism is going to be very
|
|
difficult. I am inclined to use templates in an <I>array</I> library but not in
|
|
a <I>matrix</I> library. </P>
|
|
<P>Hence I decided to implement only one element type. But the user can decide
|
|
whether this is float or double. The package assumes elements are of type Real.
|
|
The user typedefs Real to float or double. </P>
|
|
<P>It might also be worth including symmetric and triangular matrices with
|
|
extra precision elements (double or long double) to be used for storage only
|
|
and with a minimum of operations defined. These would be used for accumulating
|
|
the results of sums of squares and product matrices or multi-stage QR
|
|
triangularisations. </P>
|
|
<H3>Allow matrix expressions</H3>
|
|
<P>I want to be able to write matrix expressions the way I would on paper. So
|
|
if I want to multiply two matrices and then add the transpose of a third one I
|
|
can write something like <TT>X = A * B + C.t();</TT>. I want this expression to
|
|
be evaluated with close to the same efficiency as a hand-coded version. This is
|
|
not so much of a problem with expressions including a multiply since the
|
|
multiply will dominate the time. However, it is not so easy to achieve with
|
|
expressions with just <TT>+</TT> and <TT>-</TT>. </P>
|
|
<P>A second requirement is that temporary matrices generated during the
|
|
evaluation of an expression are destroyed as quickly as possible. </P>
|
|
<P>A desirable feature is that a certain amount of <I>intelligence</I> be
|
|
displayed in the evaluation of an expression. For example, in the expression
|
|
<TT>X = A.i() * B;</TT> where <TT>i()</TT> denotes inverse, it would be
|
|
desirable if the inverse wasn't explicitly calculated. </P>
|
|
<H3>Naming convention</H3>
|
|
<P>How are classes and public member functions to be named? As a general rule I
|
|
have spelt identifiers out in full with individual words being capitalised. For
|
|
example <I>UpperTriangularMatrix</I>. If you don't like this you can #define or
|
|
typedef shorter names. This convention means you can select an abbreviation
|
|
scheme that makes sense to you. </P>
|
|
<P>Exceptions to the general rule are the functions for transpose and inverse.
|
|
To make matrix expressions more like the corresponding mathematical formulae, I
|
|
have used the single letter abbreviations, <TT>t()</TT> and <TT>i()</TT>. </P>
|
|
<H3>Row and column index ranges</H3>
|
|
<P>In mathematical work matrix subscripts usually start at one. In C, array
|
|
subscripts start at zero. In Fortran, they start at one. Possibilities for this
|
|
package were to make them start at 0 or 1 or be arbitrary. </P>
|
|
<P>Alternatively one could specify an <I>index set</I> for indexing the rows
|
|
and columns of a matrix. One would be able to add or multiply matrices only if
|
|
the appropriate row and column index sets were identical. </P>
|
|
<P>In fact, I adopted the simpler convention of making the rows and columns of
|
|
a matrix be indexed by an integer starting at one, following the traditional
|
|
convention. In an earlier version of the package I had them starting at zero,
|
|
but even I was getting mixed up when trying to use this earlier package. So I
|
|
reverted to the more usual notation and started at 1. </P>
|
|
<H3>Element access - method and checking</H3>
|
|
<P>We want to be able to use the notation <TT>A(i,j)</TT> to specify the
|
|
<TT>(i,j)</TT>-th element of a matrix. This is the way mathematicians expect to
|
|
address the elements of matrices. I consider the notation <TT>A[i][j]</TT>
|
|
totally alien. However I include this as an option to help people converting
|
|
from C. </P>
|
|
<P>There are two ways of working out the address of <TT>A(i,j)</TT>. One is
|
|
using a <I>dope</I> vector which contains the first address of each row.
|
|
Alternatively you can calculate the address using the formula appropriate for
|
|
the structure of <TT>A</TT>. I use this second approach. It is probably slower,
|
|
but saves worrying about an extra bit of storage. </P>
|
|
<P>The other question is whether to check for <TT>i</TT> and <TT>j</TT> 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. </P>
|
|
<H3>Use iterators</H3>
|
|
<P>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 internally for some row
|
|
and column functions. </P>
|
|
<H2><A NAME="stor"></A>5.4 Data storage</H2>
|
|
<P CLASS="small"><A HREF="#mem_man">next</A> - <A HREF="#mem_man">skip</A> -
|
|
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
|
|
<H3>The stack and heap</H3>
|
|
<p>To understand how <i>newmat</i> stores matrices you need to know a little bit
|
|
about the <i>heap</i> and <i>stack</i>.</p>
|
|
<p>The data values of variables or objects in a C++ program are stored in either
|
|
of two sections of memory called the <i>stack</i> and the <i>heap</i>. Sometimes
|
|
there is more than one <i>heap</i> to cater for different sized variables.</p>
|
|
<p>If you declare an <i>automatic</i> variable</p>
|
|
<pre> int x;</pre>
|
|
<p>then the value of <i>x</i> is stored on the <i>stack</i>. As you declare more
|
|
variables the stack gets bigger. When you exit a block (i.e a section of code
|
|
delimited by curly brackets <i>{...}</i>) the memory used by the automatic
|
|
variables declared in the block is released and the <i>stack</i> shrinks.</p>
|
|
<p>When you declare a variable with <i>new</i>, for example,</p>
|
|
<pre> int* y = new int;</pre>
|
|
<p>the pointer <i>y</i> is stored on the <i>stack</i> but the value it is
|
|
pointing to is stored on the <i>heap</i>. Memory on the <i>heap</i> is not
|
|
released until the program explicitly does this with a <i>delete</i> statement</p>
|
|
<pre> delete *y;</pre>
|
|
<p>or the program exits.</p>
|
|
<p>On the <i>stack</i>, variables and objects are is always added to the end of
|
|
the <i>stack</i> and are removed in the reverse order to that in which they are
|
|
added - that is the last on will be the first off. This is not the case with the <i>
|
|
heap</i>, where the variables and objects can be removed in any order. So one
|
|
can get alternating pieces of used and unused memory. When a new variable or
|
|
object is declared on the <i>heap</i> the system needs to search for piece of
|
|
unused memory large enough to hold it. This means that storing on the <i>heap</i>
|
|
will usually be a slower process than storing on the <i>stack</i>. There is also
|
|
likely to be waste space on the <i>heap</i> because of gaps between the used
|
|
blocks of memory that are too small for the next object you want to store on the
|
|
<i>heap</i>. There is also the possibility of wasting space if you forget to
|
|
remove a variable or object on the <i>heap</i> even though you have finished
|
|
using it. However, the <i>stack</i> is usually limited to holding small objects
|
|
with size known at compile time. Large objects, objects whose size you don't
|
|
know at compile time, and objects that you want to persist after the end of the
|
|
block need to be stored on the <i>heap</i>.</p>
|
|
<p>In C++, the <i>constructor</i>/<i>destructor</i> system enables one to build
|
|
complicated objects such as matrices that behave as automatic variables stored
|
|
on the <i>stack</i>, so the programmer doesn't have to worry about deleting them
|
|
at the end of the block, but which really utilise the <i>heap</i> for storing
|
|
their data.</p>
|
|
<H3>Structure of matrix objects</H3>
|
|
<P>Each matrix object contains the basic information such as the number of rows
|
|
and columns, the amount of memory used, a status variable and a pointer to the data array which is on
|
|
the heap. So if you declare a matrix</P>
|
|
<pre> Matrix A(1000,1000);</pre>
|
|
<p>there is an small amount of memory used on the stack for storing the numbers
|
|
of rows and columns, the amount of memory used, the status variable and
|
|
the pointer together with 1,000,000 <i>Real</i> locations stored on the heap.
|
|
When you exit the block in which <i>A</i> is declared, the heap memory used by
|
|
<i>A</i> is automatically returned to the system, as well as the memory used on
|
|
the stack.</p>
|
|
<p>Of course, if you use new to declare a matrix</p>
|
|
<pre> Matrix* B = new Matrix(1000,1000);</pre>
|
|
<p>both the information about the size and the actual data are stored on heap
|
|
and not deleted until the program exits or you do an explicit delete:</p>
|
|
<pre> delete *B;</pre>
|
|
<p>If you carry out an assignment with <tt>=</tt> or <tt><<</tt> or do a
|
|
<tt>resize()</tt> the data array currently associated with a matrix is destroyed and
|
|
a new array generated. For example</p>
|
|
<pre> Matrix A(1000,1000);
|
|
Matrix B(50, 50);
|
|
... put some values in B
|
|
A = B;</pre>
|
|
<p>At the last step the heap memory associated with <i>A</i> is returned to the
|
|
system and a new block of heap memory is assigned to contain the new values.
|
|
This happens even if there is no change in the amount of memory required. </p>
|
|
<H3>One block or several</H3>
|
|
<P>The elements of the matrix are stored as a single array. Alternatives would
|
|
have been to store each row as a separate array or a set of adjacent rows as a
|
|
separate array. The present solution simplifies the program but limits the size
|
|
of matrices in 16 bit PCs that have a 64k byte limit on the size of arrays (I
|
|
don't use the <TT>huge</TT> keyword). The large arrays may also cause problems
|
|
for memory management in smaller machines. [The 16 bit PC problem has largely
|
|
gone away but it was a problem when much of <I>newmat</I> was written. Now,
|
|
occasionally I run into the 32 bit PC problem.]</P>
|
|
<H3>By row or by column or other</H3>
|
|
<P>In Fortran two dimensional arrays are stored by column. In most other
|
|
systems they are stored by row. I have followed this later convention. This
|
|
makes it easier to interface with other packages written in C but harder to
|
|
interface with those written in Fortran. This may have been a wrong decision.
|
|
Most work on the efficient manipulation of large matrices is being done in
|
|
Fortran. It would have been easier to use this work if I had adopted the
|
|
Fortran convention. </P>
|
|
<P>An alternative would be to store the elements by mid-sized rectangular
|
|
blocks. This might impose less strain on memory management when one needs to
|
|
access both rows and columns. </P>
|
|
<H3>Storage of symmetric matrices</H3>
|
|
<P>Symmetric matrices are stored as lower triangular matrices. The decision was
|
|
pretty arbitrary, but it does slightly simplify the Cholesky decomposition
|
|
program. </P>
|
|
<H2><A NAME="mem_man"></A>5.5 Memory management -
|
|
reference counting or status variable?</H2>
|
|
<P CLASS="small"><A HREF="#mem_man2">next</A> - <A HREF="#mem_man2">skip</A> -
|
|
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
|
|
<P>Consider the instruction </P>
|
|
<PRE> X = A + B + C;
|
|
</PRE>
|
|
|
|
<P>To evaluate this a simple program will add <TT>A</TT> to <TT>B</TT> putting
|
|
the total in a temporary <TT>T1</TT>. Then it will add <TT>T1</TT> to
|
|
<TT>C</TT> creating another temporary <TT>T2</TT> which will be copied into
|
|
<TT>X</TT>. <TT>T1</TT> and <TT>T2</TT> will sit around till the end of the
|
|
execution of the statement and perhaps of the block. It would be faster if the
|
|
program recognised that <TT>T1</TT> was temporary and stored the sum of
|
|
<TT>T1</TT> and <TT>C</TT> back into <TT>T1</TT> instead of creating
|
|
<TT>T2</TT> and then avoided the final copy by just assigning the contents of
|
|
<TT>T1</TT> to <TT>X</TT> 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). </P>
|
|
<P>For an instruction like </P>
|
|
<PRE> X = (A * B) + (C * D);
|
|
</PRE>
|
|
|
|
<P>we can't easily avoid one temporary being left over, so we would like this
|
|
temporary deleted as quickly as possible. </P>
|
|
<P>I provide the functionality for doing all 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. </P>
|
|
<P>An alternative or additional approach would be to use <I>reference counting
|
|
and delayed copying</I> - also known as <I>copy on write</I>. 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.
|
|
If the original matrix is deleted before either matrix is modified, in effect,
|
|
the values of the original matrix are transferred to the copy without any
|
|
actual copying taking place. This solves the difficult problem of returning an
|
|
object from a function without copying and saves the unnecessary copying in the
|
|
previous examples. </P>
|
|
<P>There are downsides to the delayed copying approach. Typically, for delayed
|
|
copying one uses a structure like the following: </P>
|
|
<PRE> Matrix
|
|
|
|
|
+------> Array Object
|
|
| |
|
|
| +------> Data array
|
|
| |
|
|
| +------- Counter
|
|
|
|
|
+------ Dimension information
|
|
|
|
</PRE>
|
|
|
|
<P>where the arrows denote a pointer to a data structure. If one wants to
|
|
access the <I>Data array</I> one will need to track through two pointers. If
|
|
one is going to write, one will have to check whether one needs to copy first.
|
|
This is not important when one is going to access the whole array, say, for a
|
|
add operation. But if one wants to access just a single element, then it
|
|
imposes a significant additional overhead on that operation. Any subscript
|
|
operation would need to check whether an update was required - even read since
|
|
it is hard for the compiler to tell whether a subscript access is a read or
|
|
write. </P>
|
|
<P>Some matrix libraries don't bother to do this. So if you write <TT>A =
|
|
B;</TT> and then modify an element of one of <TT>A</TT> or <TT>B</TT>, then the
|
|
same element of the other is also modified. I don't think this is acceptable
|
|
behaviour. </P>
|
|
<P>Delayed copy does not provide the additional functionality of my approach
|
|
but I suppose it would be possible to have both delayed copy and tagging
|
|
temporaries. </P>
|
|
<P>My approach does not automatically avoid all copying. In particular, you
|
|
need use a special technique to return a matrix from a function without
|
|
copying. </P>
|
|
<H2><A NAME="mem_man2"></A>5.6 Memory management -
|
|
accessing contiguous locations</H2>
|
|
<P CLASS="small"><A HREF="#evalx">next</A> - <A HREF="#evalx">skip</A> -
|
|
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
|
|
<P>Modern computers work faster if one accesses memory by running through
|
|
contiguous locations rather than by jumping around all over the place. Newmat
|
|
stores matrices <A HREF="#stor">by rows</A> so that algorithms that access
|
|
memory by running along rows will tend to work faster than one that runs down
|
|
columns. A number of the algorithms used in Newmat were developed before this
|
|
was an issue and so are not as efficient as possible. </P>
|
|
<P>I have gradually upgrading the algorithms to access memory by rows. The
|
|
following table shows the current status of this process.</P>
|
|
<TABLE WIDTH="100%" BORDER="1">
|
|
<TR>
|
|
<TH ALIGN="LEFT">Function</TH>
|
|
<TH ALIGN="LEFT">Contiguous memory access</TH>
|
|
<TH ALIGN="LEFT">Comment</TH>
|
|
</TR>
|
|
<TR>
|
|
<TD>Add, subtract</TD>
|
|
<TD>Yes</TD>
|
|
<TD> </TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>Multiply</TD>
|
|
<TD>Yes</TD>
|
|
<TD> </TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>Concatenate</TD>
|
|
<TD>Yes</TD>
|
|
<TD> </TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>Transpose</TD>
|
|
<TD>No</TD>
|
|
<TD> </TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>Invert and solve</TD>
|
|
<TD>Yes</TD>
|
|
<TD>Mostly</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>Cholesky</TD>
|
|
<TD>Yes</TD>
|
|
<TD> </TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>QRZ, QRZT</TD>
|
|
<TD>Yes</TD>
|
|
<TD> </TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>SVD</TD>
|
|
<TD>No</TD>
|
|
<TD> </TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>Jacobi</TD>
|
|
<TD>No</TD>
|
|
<TD>Not an issue; used only for smaller matrices</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>Eigenvalues</TD>
|
|
<TD>No</TD>
|
|
<TD> </TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>Sort</TD>
|
|
<TD>Yes</TD>
|
|
<TD>Quick-sort is naturally good</TD>
|
|
</TR>
|
|
<TR>
|
|
<TD>FFT</TD>
|
|
<TD>?</TD>
|
|
<TD>Could be improved?</TD>
|
|
</TR>
|
|
</TABLE>
|
|
<H2><A NAME="evalx"></A>5.7 Evaluation of expressions -
|
|
lazy evaluation</H2>
|
|
<P CLASS="small"><A HREF="#explode">next</A> - <A HREF="#explode">skip</A> -
|
|
<A HREF="#design">up</A> - <A HREF="#top">start</A></P>
|
|
<P>Consider the instruction </P>
|
|
<PRE> X = B - X;
|
|
</PRE>
|
|
|
|
<P>A simple program will subtract <TT>X</TT> from <TT>B</TT>, store the result
|
|
in a temporary <TT>T1</TT> and copy <TT>T1</TT> into <TT>X</TT>. It would be
|
|
faster if the program recognised that the result could be stored directly into
|
|
<TT>X</TT>. This would happen automatically if the program could look at the
|
|
instruction first and mark <TT>X</TT> as temporary. </P>
|
|
<P>C programmers would expect to avoid the same problem with </P>
|
|
<PRE> X = X - B;
|
|
</PRE>
|
|
|
|
<P>by using an operator <TT>-=</TT> </P>
|
|
<PRE> X -= B;
|
|
</PRE>
|
|
|
|
<P>However this is an unnatural notation for non C users and it may be nicer to
|
|
write <TT>X = X - B</TT>; and know that the program will carry out the
|
|
simplification. </P>
|
|
<P>Another example where this intelligent analysis of an instruction is helpful
|
|
is in </P>
|
|
<PRE> X = A.i() * B;
|
|
</PRE>
|
|
|
|
<P>where <TT>i()</TT> 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. </P>
|
|
<P>I regard this interpretation of <TT>A.i() * B</TT> as just providing a
|
|
convenient notation. The objective is not primarily to correct the errors of
|
|
people who are unaware of the inefficiency of <TT>A.i() * B</TT> if interpreted
|
|
literally. </P>
|
|
<P>There is a third reason for the two-stage evaluation of expressions and this
|
|
is probably the most important one. In C++ it is quite hard to return an
|
|
expression from a function such as (<TT>*</TT>, <TT>+</TT> etc) without a copy.
|
|
This is particularly the case when an assignment (<TT>=</TT>) is involved. The
|
|
mechanism described here provides one way for avoiding this in matrix
|
|
expressions. </P>
|
|
<P>The C++ standard (section 12.8/15) allows the compiler to optimise away the
|
|
copy when returning an object from a function (but there will still be one copy
|
|
is an <span lang="en-nz">assignment </span>(=) is involved). This means special handling of returns from a
|
|
function is less important when a modern optimising compiler is being
|
|
used. </P>
|
|
<P>To carry out this <I>intelligent</I> 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 <TT>(A+B)*C</TT> is
|
|
represented by a tree </P>
|
|
<PRE><TT>
|
|
*
|
|
/ \
|
|
+ C
|
|
/ \
|
|
A B
|
|
</TT></PRE>
|
|
|
|
<P>Rather than adding <TT>A</TT> and <TT>B</TT> the <TT>+</TT> operator yields
|
|
an object of a class <I>AddedMatrix</I> which is just a pair of pointers to
|
|
<TT>A</TT> and <TT>B</TT>. Then the <TT>*</TT> operator yields a
|
|
<I>MultipliedMatrix</I> which is a pair of pointers to the <I>AddedMatrix</I>
|
|
and <TT>C</TT>. The tree is examined for any simplifications and then evaluated
|
|
recursively. </P>
|
|
<P>Further possibilities not yet included are to recognise <TT>A.t()*A</TT> and
|
|
<TT>A.t()+A</TT> as symmetric or to improve the efficiency of evaluation of
|
|
expressions like <TT>A+B+C</TT>, <TT>A*B*C</TT>, <TT>A*B.t()</TT> (<TT>t()</TT>
|
|
denotes transpose). </P>
|
|
<P>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 </P>
|
|
<PRE> Matrix M;
|
|
DiagonalMatrix D;
|
|
....;
|
|
D = M;
|
|
</PRE>
|
|
|
|
<P>We don't allow conversions using <TT>=</TT> when information would be lost.
|
|
Such errors will be detected when the statement is executed. </P>
|
|
<H2><A NAME="explode"></A>5.8 How to overcome an
|
|
explosion in number of operations</H2>
|
|
<P CLASS="small"><A HREF="#destr">next</A> - <A HREF="#destr">skip</A> -
|
|
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
|
|
<P>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. </P>
|
|
<P>With <I>n</I> types of matrices the binary operations will each require
|
|
<I>n</I>-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. </P>
|
|
<P>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. </P>
|
|
<P>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. </P>
|
|
<P>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. </P>
|
|
<P>Generally the method works well. However symmetric matrices are not always
|
|
handled very efficiently (yet) since complete rows are not stored explicitly.
|
|
</P>
|
|
<P>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.
|
|
</P>
|
|
<P>However it must also be admitted that there is a substantial overhead in the
|
|
approach adopted in the present package for small matrices. The test program
|
|
developed for the original version of the package takes 30 to 50% longer to run
|
|
with the current version (though there may be some other reasons for this).
|
|
This is for matrices in the range 6x6 to 10x10. </P>
|
|
<P>To try to improve the situation a little I do provide an ordinary matrix
|
|
multiplication routine for the case when all the matrices involved are
|
|
rectangular. </P>
|
|
<H2><A NAME="destr"></A>5.9 Destruction of
|
|
temporaries</H2>
|
|
<P CLASS="small"><A HREF="#calc">next</A> - <A HREF="#calc">skip</A> -
|
|
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
|
|
<P>Versions before version 5 of newmat did not work correctly with Gnu C++
|
|
(version 5 or earlier). 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++. </P>
|
|
<P>However early version Gnu C++ destroys temporary structures as soon as the
|
|
function that accesses them finishes. The other compilers wait until the end of
|
|
the current expression or current block. 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. </P>
|
|
<P>Now that the C++ standards committee has said that temporary structures
|
|
should not be destroyed before a statement finishes, I suggest using the stack,
|
|
because of the difficulty of managing exceptions with the heap version. </P>
|
|
<H2><A NAME="calc"></A>5.10 A calculus of matrix
|
|
types</H2>
|
|
<P CLASS="small"><A HREF="#pointer">next</A> - <A HREF="#pointer">skip</A> -
|
|
<A HREF="#design">up</A> - <A HREF="#top">start</A></P>
|
|
<P>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 an intermediate result. To assist with this, a class MatrixType is
|
|
defined. Operators <TT>+</TT>, <TT>-</TT>, <TT>*</TT>, <TT>>=</TT> are
|
|
defined to calculate the types of the results of expressions or to check that
|
|
conversions are legal. </P>
|
|
<P>Early versions of <I>newmat</I> stored the types of the results of
|
|
operations in a table. So, for example, if you multiplied an
|
|
UpperTriangularMatrix by a LowerTriangularMatrix, <I>newmat</I> would look up
|
|
the table and see that the result was of type Matrix. With this approach the
|
|
<A HREF="#explode">exploding</A> number of operations problem recurred although
|
|
not as seriously as when code had to be written for each pair of types. But
|
|
there was always the suspicion that somewhere, there was an error in one of
|
|
those 9x9 tables, that would be very hard to find. And the problem would get
|
|
worse as additional matrix types or operators were included. </P>
|
|
<P>The present version of <I>newmat</I> solves the problem by assigning
|
|
<I>attributes</I> such as <I>diagonal</I> or <I>band</I> or <I>upper
|
|
triangular</I> to each matrix type. Which attributes a matrix type has, is
|
|
stored as bits in an integer. As an example, the DiagonalMatrix type has the
|
|
bits corresponding to <I>diagonal</I>, <I>symmetric</I> and <I>band</I> equal
|
|
to 1. By looking at the attributes of each of the operands of a binary
|
|
operator, the program can work out the attributes of the result of the
|
|
operation with simple bitwise operations. Hence it can deduce an appropriate
|
|
type. The <I>symmetric</I> attribute is a minor problem because
|
|
<I>symmetric</I> * <I>symmetric</I> does not yield <I>symmetric</I> unless both
|
|
operands are <I>diagonal</I>. But otherwise very simple code can be used to
|
|
deduce the attributes of the result of a binary operation. </P>
|
|
<P>Tables of the types resulting from the binary operators are output at the
|
|
beginning of the <A HREF="#testing">test</A> program. </P>
|
|
<H2><A NAME="pointer"></A>5.11 Pointer arithmetic</H2>
|
|
<P CLASS="small"><A HREF="#err_hand">next</A> - <A HREF="#err_hand">skip</A> -
|
|
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
|
|
<P>Suppose you do something like </P>
|
|
<PRE> int* y = new int[100];
|
|
y += 200; // y points to something outside the array
|
|
// y is <B>never</B> accessed
|
|
</PRE>
|
|
|
|
<P>Then the standard says that the behaviour of the program is <I>undefined</I>
|
|
even if <TT>y</TT> is never accessed. (You are allowed to calculate a pointer
|
|
value one location beyond the end of the array). In practice, a program like
|
|
this does not cause any problems with any compiler I have come across and
|
|
no-one has reported any such problems to me. </P>
|
|
<P>However, this <I>error</I> is detected by Borland's <I>Code Guard</I>
|
|
bound's checker and this makes it very difficult to use this to use <I>Code
|
|
Guard</I> to detect other problems since the output is swamped by reports of
|
|
this <I>error</I>. </P>
|
|
<P>Now consider </P>
|
|
<PRE> int* y = new int[100];
|
|
y += 200; // y points to something outside the array
|
|
y -= 150; // y points to something inside the array
|
|
// y <B>is</B> accessed
|
|
</PRE>
|
|
|
|
<P>Again this is not strictly correct but does not seem to cause a problem. But
|
|
it is much more doubtful than the previous example. </P>
|
|
<P>I removed most instances of the second version of the problem from Newmat09.
|
|
Hopefully the remainder of these instances were removed from the current
|
|
version of Newmat10. In addition, most instances of the first version of the
|
|
problem have also been fixed. </P>
|
|
<P>There is one exception. The interface to the <A HREF="#nric">Numerical
|
|
Recipes in C</A> does still contain the second version of the problem. This is
|
|
inevitable because of the way Numerical Recipes in C stores vectors and
|
|
matrices. If you are running the <A HREF="#testing">test program</A> with a
|
|
bounds checking program, edit <TT>tmt.h</TT> to disable the testing of the NRIC
|
|
interface. </P>
|
|
<P>The rule does does cause a problem for authors of matrix and
|
|
multidimensional array packages. If we want to run down a column of a matrix we
|
|
would like to do something like </P>
|
|
<PRE> // set values of column 1
|
|
Matrix A;
|
|
... set dimensions and put values in A
|
|
Real* a = A.Store(); // points to first element
|
|
int nr = A.Nrows(); // number of rows
|
|
int nc = A.Ncols(); // number of columns
|
|
while (nr--)
|
|
{
|
|
*a = something to put in first element of row
|
|
a += nc; // jump to next element of column
|
|
}
|
|
</PRE>
|
|
|
|
<P>If the matrix has more than one column the last execution of <TT>a +=
|
|
nc;</TT> will run off the end of the space allocated to the matrix and we'll
|
|
get a bounds error report. </P>
|
|
<P>Instead we have to use a program like </P>
|
|
<PRE> // set values of column 1
|
|
Matrix A;
|
|
... set dimensions and put values in A
|
|
Real* a = A.Store(); // points to first element
|
|
int nr = A.Nrows(); // number of rows
|
|
int nc = A.Ncols(); // number of columns
|
|
if (nr != 0)
|
|
{
|
|
for(;;)
|
|
{
|
|
*a = something to put in first element of row
|
|
if (!(--nr)) break;
|
|
a += nc; // jump to next element of column
|
|
}
|
|
}
|
|
</PRE>
|
|
|
|
<P>which is more complicated and consequently introduces more chance of error.
|
|
</P>
|
|
<H2><A NAME="err_hand"></A>5.12 Error handling
|
|
</H2>
|
|
<P CLASS="small"><A HREF="#sparse">next</A> - <A HREF="#sparse">skip</A> -
|
|
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
|
|
<P>The library now does have a moderately graceful exit from errors. One can
|
|
use either the simulated exceptions or the compiler supported exceptions. When
|
|
newmat08 was released (in 1995), compiler exception handling in the compilers I
|
|
had access to was unreliable. I recommended you used my simulated exceptions.
|
|
In 1997 compiler supported exceptions seemed to work on a variety of
|
|
compilers - but not all compilers. This is still true in 2001. Try using the
|
|
compiler supported exceptions if you have a recent compiler, but if you are
|
|
getting strange crashes or errors try going back to my simulated exceptions.
|
|
</P>
|
|
<P>The approach in the present library, attempting to simulate C++ exceptions,
|
|
is not completely satisfactory, but seems a good interim solution for those who
|
|
cannot use compiler supported exceptions. People who don't want exceptions in
|
|
any shape or form, can set the option to exit the program if an exception is
|
|
thrown. </P>
|
|
<P>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 completely 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. </P>
|
|
<H2><A NAME="sparse"></A>5.13 Sparse matrices</H2>
|
|
<P CLASS="small"><A HREF="#comp_mat">next</A> - <A HREF="#comp_mat">skip</A> -
|
|
<A HREF="#design">up</A> - <A HREF="#top">start</A></P>
|
|
<P>The library does not support sparse matrices. </P>
|
|
<P>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. </P>
|
|
<P>Sparse matrices are important for people solving large sets of differential
|
|
equations as well as being important for statistical and operational research
|
|
applications. </P>
|
|
<P>But there are packages being developed specifically for sparse matrices and
|
|
these might present the best approach, at least where sparse matrices are the
|
|
main interest. </P>
|
|
<H2><A NAME="comp_mat"></A>5.14 Complex matrices</H2>
|
|
<P CLASS="small"><A HREF="#top">next</A> - <A HREF="#top">skip</A> -
|
|
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
|
|
<P>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. </P>
|
|
<P>This probably means making new versions of the basic row and column
|
|
operations for Real*Complex, Complex*Complex, Complex*Real and similarly for
|
|
<TT>+</TT> and <TT>-</TT>. This would be OK, except that if I also want to do
|
|
this for sparse matrices, then when you put these together, the whole thing
|
|
will get out of hand. </P>
|
|
<P>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. But there is going to be a problem
|
|
with accessing elements and it does not seem possible to solve this in an
|
|
entirely satisfactory way. </P>
|
|
<P>Complex matrices are used extensively by electrical engineers and physicists
|
|
and really should be fully supported in a comprehensive package. </P>
|
|
<P>You can simulate most complex operations by representing <TT>Z = X + iY</TT>
|
|
by </P>
|
|
<PRE> / X Y \
|
|
\ -Y X /
|
|
</PRE>
|
|
|
|
<P>Most matrix operations will simulate the corresponding complex operation,
|
|
when applied to this matrix. But, of course, this matrix is essentially twice as big as you
|
|
would need with a genuine complex matrix library.</P>
|
|
<P CLASS="small"><A HREF="#top">return to top</A><BR>
|
|
<A HREF="ol_doc.htm">return to online documentation page</A></P>
|
|
<P> </P>
|
|
<P> </P>
|
|
<P> </P>
|
|
<P> </P>
|
|
<P> </P>
|
|
<P> </P>
|
|
<P> </P>
|
|
</BODY>
|
|
</HTML> |