From: robertd@kauri.vuw.ac.nz (Robert Davies)
Subject: Newmat06, a matrix package in C++
Organization: Victoria University of Wellington
Date: Wed, 16 Dec 1992 12:02:44 GMT

I have released a new version of my matrix package in C++. This version
is called newmat06.

The package supports classes Matrix, UpperTriangularMatrix,
LowerTriangularMatrix, DiagonalMatrix, SymmetricMatrix, BandMatrix,
UpperBandMatrix, LowerBandMatrix, SymmetricBandMatrix, RowVector,
ColumnVector.

Only one element type (float or double) is supported.

The package includes the operations *, +, -, inverse, transpose,
conversion between types, submatrix, determinant, Cholesky
decomposition, Householder triangularisation, singular value
decomposition, symmetric eigenvalue analysis, fast Fourier transform,
sort, print, an interface to "Numerical Recipes in C" programs, and an
emulation of exceptions.

It is suitable for matrices in the size range 15 by 15 up to the maximum
size that can conveniently be stored in a single block of memory.

The package works with AT&T C++ (versions 2.1 and 3), Gnu G++ (version
2.2), Borland C++ (version 3.1), Microsoft C++ (7.0) (with changes - see
below), and mostly with Zortech (version 3.04).

It is available on Simtel20 (192.88.110.20) in the pd1:<msdos.cplusplus>
directory in .zip format, and on the archives of comp.sources.misc,
volume 34 issues 7 to 13 (shar format). Look for it on
wuarchive.wustl.edu (128.252.135.4). It is also on Compuserve on the
Borland bulletin board. And it is being distributed by the C users'
group.

Here are some notes on the package.

The documentation says the package is good for matrices in the range 4 x
4 up to the maximum size that can be stored in a single block of memory.
In fact the package becomes inefficient for matrices smaller than about
15 x 15 (the actual figure will depend on what you are doing and the
speed of your floating point operations). The upper bound is set by the
maximum size of an array and the maximum value of an "int".


Gnu 2.3.x
---------

The latest version of Gnu picked up a couple of errors in the package.
You should probably make these changes whatever compiler you are using.

In except.h change the destructor for class Janitor to

   virtual ~Janitor();

In newmat.h, class InvertedMatrix, delete the final const in the
second definition of operator*

#ifndef TEMPS_DESTROYED_QUICKLY
   SolvedMatrix operator*(const BaseMatrix&) const;  // inverse(A) * B
#else
   SolvedMatrix& operator*(const BaseMatrix&);       // inverse(A) * B
#endif


and make the corresponding changes in newmat6.cxx


SolvedMatrix& InvertedMatrix::operator*(const BaseMatrix& bmx)
{
   REPORT
   SolvedMatrix* x;
   Try { x = new SolvedMatrix(bm, &bmx); MatrixErrorNoSpace(x); }
   CatchAll { delete this; Throw(); }
   delete this;                // since we are using bm rather than this
   return *x;
}


I have not checked that G++ 2.3.x works as it is not yet loaded here.
Please report any further problems to me. One correspondent reported a
compiler error in tmte.cxx. Has anyone else had this error?


Zortech 3.04
------------

1. Replace controlw.h by the following


//$$ controlw.h                Control word class

#ifndef CONTROL_WORD_LIB
#define CONTROL_WORD_LIB 0

// for organising an int as a series of bits which indicate whether an
// option is on or off.

class ControlWord
{
protected:
   int cw;                                      // the control word
public:
   ControlWord() : cw(0) {}                     // do nothing
   ControlWord(int i) : cw(i) {}                // load an integer

      // select specific bits (for testing at least one set)
   ControlWord operator*(const ControlWord& i) const
      { return ControlWord(cw & i.cw); }
   void operator*=(const ControlWord& i)  { cw &= i.cw; }

      // set bits
   ControlWord operator+(const ControlWord& i) const
      { return ControlWord(cw | i.cw); }
   void operator+=(const ControlWord& i)  { cw |= i.cw; }

      // reset bits
   ControlWord operator-(const ControlWord& i) const
      { return ControlWord(cw - (cw & i.cw)); }
   void operator-=(const ControlWord& i) { cw -= (cw & i.cw); }

      // check if all of selected bits set or reset
   Boolean operator>=(const ControlWord& i) const { return (cw & i.cw) == i.cw; }
   Boolean operator<=(const ControlWord& i) const { return (cw & i.cw) == cw; }

      // flip selected bits
   ControlWord operator^(const ControlWord& i) const
      { return ControlWord(cw ^ i.cw); }
   ControlWord operator~() const { return ControlWord(~cw); }

      // convert to integer
   int operator+() const { return cw; }
   int operator!() const { return cw==0; }
   FREE_CHECK(ControlWord)
};


#endif


2. Do not #define TEMPS_DESTROYED_QUICKLY; do not #define Version21

3. Do not use exceptions or the NRIC interface.

4. Use the large model


Microsoft C/C++, 7.0
--------------------

Include the following block in include.h and delete the "const"s as for
G++.

#define __MSC__


#ifdef __MSC__                          // Microsoft
   #include <stdlib.h>

   typedef int jmp_buf[9];
   extern "C"
   {
      int __cdecl setjmp(jmp_buf);
      void __cdecl longjmp(jmp_buf, int);
   }

   #ifdef WANT_STREAM
      #include <iostream.h>
      #include <iomanip.h>
   #endif
   #ifdef WANT_MATH
      #include <math.h>
      #include <float.h>
   #endif
#endif


In newmat.h you must include to following statement

void operator=(const Matrix& m) { operator=((const BaseMatrix&)m); }

in every class where operator=(const BaseMatrix&) is defined. "Matrix"
is to be replaced by the class name of the class you are amending.

If you don't make this change Microsoft uses the default assignment
operator. (What should it do, ARM doesn't make it clear?)

In jacobi.cxx, function void Jacobi(const SymmetricMatrix& X,
DiagonalMatrix& D, SymmetricMatrix& A, Matrix& V, Boolean eivec), make
the following changes:

Immediately after the initial { insert the line

   Real epsilon = FloatingPointPrecision::Epsilon();

change
	    if (i>4 && adp+g==adp && adq+g==adq) apq = 0.0;
to
	    if (i>4 && g < epsilon*adp && g < epsilon*adq) apq = 0.0;
and change
	       if (ah+g==ah) t = apq / h;
to
	       if (g < epsilon*ah) t = apq / h;

Do #define TEMPS_DESTROYED_QUICKLY as the cleanup after an exception
seems to work better.

Microsoft warns that setjmp and longjmp, which are used in the emulation
of the exceptions, are not compatible with C++. So, be warned.


AT&T (and others?)
------------------

If you are using this with interviews you may get a conflict with Catch.
Either #undefine Catch or replace Catch with CATCH throughout my
package.

I have assumed that math.h contains a function "abs". I understand I
should be using iabs. So, if abs causes a problem, replace it with iabs.
(In the HP unix workstation, the problem was the opposite, abs was
defined in both math.h and stdlib.h. You need to delete it from one of
these.)