Source Code 1994 March/Source_Code_CD-ROM_Walnut_Creek_March_1994.iso/compsrcs/misc/volume34/newmat06/part03

Text File | 1992-12-06 | 54.7 KB | 1,699 lines

Newsgroups: comp.sources.misc
From: [e-mail] (Robert Davies)
Subject:  v34i009:  newmat06 - A matrix package in C++, Part03/07
Message-ID: <[e-mail]>
X-Md4-Signature: 326a18926a36885cc99e1edc2bf7391f
Date: Sun, 6 Dec 1992 04:52:59 GMT
Approved: [e-mail]

Submitted-by: [e-mail] (Robert Davies)
Posting-number: Volume 34, Issue 9
Archive-name: newmat06/part03
Environment: C++
Supersedes: newmat04: Volume 26, Issue 87-91

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#        "End of archive 3 (of 7)."
# Contents:  newmat.h newmatc.txt tmt.cxx tmt.mak
# Wrapped by robert@kea on Thu Dec  3 23:02:17 1992
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'newmat.h' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'newmat.h'\"
else
echo shar: Extracting \"'newmat.h'\" \(42428 characters\)
sed "s/^X//" >'newmat.h' <<'END_OF_FILE'
X//$$ newmat.h           definition file for new version of matrix package
X
X// Copyright (C) 1991,2: R B Davies
X
X#ifndef MATRIX_LIB
X#define MATRIX_LIB 0
X
X#ifdef NO_LONG_NAMES
X#define UpperTriangularMatrix UTMatrix
X#define LowerTriangularMatrix LTMatrix
X#define SymmetricMatrix SMatrix
X#define DiagonalMatrix DMatrix
X#endif
X
X#ifndef TEMPS_DESTROYED_QUICKLY
X#define ReturnMatrix ReturnMatrixX
X#else
X#define ReturnMatrix ReturnMatrixX&
X#endif
X
X#include "boolean.h"
X#include "except.h"
X
X/**************************** general utilities ****************************/
X
Xclass GeneralMatrix;
Xvoid MatrixErrorNoSpace(void*);                 // no space handler
X
Xclass LogAndSign
X// Return from LogDeterminant function
X//    - value of the log plus the sign (+, - or 0)
X{
X   Real log_value;
X   int sign;
Xpublic:
X   LogAndSign() { log_value=0.0; sign=1; }
X   LogAndSign(Real);
X   void operator*=(Real);
X   void ChangeSign() { sign = -sign; }
X   Real LogValue() const { return log_value; }
X   int Sign() const { return sign; }
X   Real Value() const;
X   FREE_CHECK(LogAndSign)
X};
X
X// the following class is for counting the number of times a piece of code
X// is executed. It is used for locating any code not executed by test
X// routines. Use turbo GREP locate all places this code is called and
X// check which ones are not accessed.
X// Somewhat implementation dependent as it relies on "cout" still being
X// present when ExeCounter objects are destructed.
X
Xclass ExeCounter
X{
X   int line;                                    // code line number
X   int fileid;                                  // file identifier
X   long nexe;                                   // number of executions
X   static int nreports;                         // number of reports
Xpublic:
X   ExeCounter(int,int);
X   void operator++() { nexe++; }
X   ~ExeCounter();                               // prints out reports
X};
X
X
X/**************************** class MatrixType *****************************/
X
X// Is used for finding the type of a matrix resulting from the binary operations
X// +, -, * and identifying what conversions are permissible.
X// This class must be updated when new matrix types are added.
X
Xclass GeneralMatrix;                            // defined later
Xclass BaseMatrix;                               // defined later
X
Xclass MatrixType
X{
Xpublic:
X   enum Attribute {  Valid     = 1,
X                     Symmetric = 2,
X                     Band      = 4,
X                     Upper     = 8,
X                     Lower     = 16,
X                     LUDeco    = 32,
X                     OneRow    = 64,
X                     OneCol    = 128 };
X
X   enum            { US = 0,
X                     UT = Valid + Upper,
X                     LT = Valid + Lower,
X                     Rt = Valid,
X                     Sm = Valid + Symmetric,
X                     Dg = Valid + Band + Lower + Upper + Symmetric,
X             RV = Valid + OneRow,
X             CV = Valid + OneCol,
X             BM = Valid + Band,
X             UB = Valid + Band + Upper,
X             LB = Valid + Band + Lower,
X             SB = Valid + Band + Symmetric,
X             Ct = Valid + LUDeco,
X             BC = Valid + Band + LUDeco,
X                   };
X
X   enum { USX,UTX,LTX,RtX,SmX,DgX,RVX,CVX,BMX,UBX,LBX,SBX,CtX,BCX };
X
X   static nTypes() { return 11; }              // number of different types
X                           // exclude Ct, US, BC
Xpublic:
X   int attribute;
Xpublic:
X   MatrixType operator+(const MatrixType& t) const;
X   MatrixType operator*(const MatrixType&) const;
X   Boolean operator>=(const MatrixType& t) const;
X   Boolean operator==(const MatrixType& t) const
X      { return (attribute == t.attribute); }
X   Boolean operator!=(const MatrixType& t) const { return !operator==(t); }
X   Boolean operator!() const { return (attribute & Valid) == 0; }
X   MatrixType i() const;                       // type of inverse
X   MatrixType t() const;                       // type of transpose
X   MatrixType AddEqualEl() const;              // Add constant to matrix
X   MatrixType MultRHS() const;                 // type for rhs of multiply
X   MatrixType sub() const;                     // type of submatrix
X   MatrixType ssub() const;                    // type of sym submatrix
X//   MatrixType (Type tx) : attribute((int)tx) {}
X//                           // (& doesn't work with AT&T)
X   MatrixType () {}
X   MatrixType (int i) : attribute(i) {}
X   GeneralMatrix* New() const;                 // new matrix of given type
X   GeneralMatrix* New(int,int,BaseMatrix*) const;
X                                               // new matrix of given type
X   operator char*() const;                     // for printing type
X   int BestFit() const;
X   FREE_CHECK(MatrixType)
X};
X
Xvoid TestTypeAdd();                            // test +
Xvoid TestTypeMult();                           // test *
Xvoid TestTypeOrder();                          // test >=
X
X
X/************************* class MatrixBandWidth ***********************/
X
Xclass MatrixBandWidth
X{
Xpublic:
X   int lower;
X   int upper;
X   MatrixBandWidth(const int l, const int u) : lower(l), upper (u) {}
X   MatrixBandWidth(const int i) : lower(i), upper(i) {}
X   MatrixBandWidth operator+(const MatrixBandWidth&) const;
X   MatrixBandWidth operator*(const MatrixBandWidth&) const;
X   MatrixBandWidth t() const { return MatrixBandWidth(upper,lower); }
X   Boolean operator==(const MatrixBandWidth& bw) const
X      { return (lower == bw.lower) && (upper == bw.upper); }
X   FREE_CHECK(MatrixBandWidth)
X};
X
X
X/*********************** Array length specifier ************************/
X
X// This class is introduced to avoid constructors such as
X//   ColumnVector(int)
X// being used for conversions
X
Xclass ArrayLengthSpecifier
X{
X   int value;
Xpublic:
X   int Value() const { return value; }
X   ArrayLengthSpecifier(int l) : value(l) {}
X   FREE_CHECK(ArrayLengthSpecifier)
X};
X
X
X/*************************** Matrix routines ***************************/
X
X
Xclass MatrixRowCol;                             // defined later
Xclass MatrixRow;
Xclass MatrixCol;
X
Xclass GeneralMatrix;                            // defined later
Xclass AddedMatrix;
Xclass MultipliedMatrix;
Xclass SubtractedMatrix;
Xclass SolvedMatrix;
Xclass ShiftedMatrix;
Xclass ScaledMatrix;
Xclass TransposedMatrix;
Xclass NegatedMatrix;
Xclass InvertedMatrix;
Xclass RowedMatrix;
Xclass ColedMatrix;
Xclass DiagedMatrix;
Xclass MatedMatrix;
Xclass GetSubMatrix;
Xclass ConstMatrix;
Xclass ReturnMatrixX;
Xclass Matrix;
Xclass nricMatrix;
Xclass RowVector;
Xclass ColumnVector;
Xclass SymmetricMatrix;
Xclass UpperTriangularMatrix;
Xclass LowerTriangularMatrix;
Xclass DiagonalMatrix;
Xclass CroutMatrix;
Xclass BandMatrix;
Xclass LowerBandMatrix;
Xclass UpperBandMatrix;
Xclass SymmetricBandMatrix;
Xclass LinearEquationSolver;
X
Xstatic MatrixType MatrixTypeUnSp(MatrixType::US);
X                        // AT&T needs this
X
Xclass BaseMatrix : public Janitor               // base of all matrix classes
X{
Xprotected:
X//   BaseMatrix() {}
X   virtual ~BaseMatrix() {}
X   virtual int search(const BaseMatrix*) const = 0;
X                        // count number of times matrix
X                        // is referred to
X   virtual MatrixType Type() const = 0;         // type of a matrix
X//   virtual int NrowsV()  const = 0;
X//   virtual int NcolsV()  const = 0;
Xpublic:
X#ifndef GXX
X   virtual GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp) = 0;
X                        // evaluate temporary
X#else
X   virtual GeneralMatrix* Evaluate(MatrixType mt) = 0;
X   GeneralMatrix* Evaluate() { return Evaluate(MatrixTypeUnSp); }
X#endif
X//   void MatrixParameters(int& nr, int& nc, MatrixType& mt)
X//      { nr = NrowsV(); nc = NcolsV(); mt = Type(); }
X#ifndef TEMPS_DESTROYED_QUICKLY
X   AddedMatrix operator+(const BaseMatrix&) const;    // results of operations
X   MultipliedMatrix operator*(const BaseMatrix&) const;
X   SubtractedMatrix operator-(const BaseMatrix&) const;
X   ShiftedMatrix operator+(Real) const;
X   ScaledMatrix operator*(Real) const;
X   ScaledMatrix operator/(Real) const;
X   ShiftedMatrix operator-(Real) const;
X   TransposedMatrix t() const;
X//   TransposedMatrix t;
X   NegatedMatrix operator-() const;                   // change sign of elements
X   InvertedMatrix i() const;
X//   InvertedMatrix i;
X   RowedMatrix AsRow() const;
X   ColedMatrix AsColumn() const;
X   DiagedMatrix AsDiagonal() const;
X   MatedMatrix AsMatrix(int,int) const;
X   GetSubMatrix SubMatrix(int,int,int,int) const;
X   GetSubMatrix SymSubMatrix(int,int) const;
X   GetSubMatrix Row(int) const;
X   GetSubMatrix Rows(int,int) const;
X   GetSubMatrix Column(int) const;
X   GetSubMatrix Columns(int,int) const;
X#else
X   AddedMatrix& operator+(const BaseMatrix&) const;    // results of operations
X   MultipliedMatrix& operator*(const BaseMatrix&) const;
X   SubtractedMatrix& operator-(const BaseMatrix&) const;
X   ShiftedMatrix& operator+(Real) const;
X   ScaledMatrix& operator*(Real) const;
X   ScaledMatrix& operator/(Real) const;
X   ShiftedMatrix& operator-(Real) const;
X   TransposedMatrix& t() const;
X//   TransposedMatrix& t;
X   NegatedMatrix& operator-() const;                   // change sign of elements
X   InvertedMatrix& i() const;
X//   InvertedMatrix& i;
X   RowedMatrix& AsRow() const;
X   ColedMatrix& AsColumn() const;
X   DiagedMatrix& AsDiagonal() const;
X   MatedMatrix& AsMatrix(int,int) const;
X   GetSubMatrix& SubMatrix(int,int,int,int) const;
X   GetSubMatrix& SymSubMatrix(int,int) const;
X   GetSubMatrix& Row(int) const;
X   GetSubMatrix& Rows(int,int) const;
X   GetSubMatrix& Column(int) const;
X   GetSubMatrix& Columns(int,int) const;
X#endif
X   Real AsScalar() const;                      // conversion of 1 x 1 matrix
X   virtual LogAndSign LogDeterminant() const;
X   virtual Real SumSquare() const;
X   virtual Real SumAbsoluteValue() const;
X   virtual Real MaximumAbsoluteValue() const;
X   virtual Real Trace() const;
X   Real Norm1() const;
X   Real NormInfinity() const;
X   virtual MatrixBandWidth BandWidth() const;  // bandwidths of band matrix
X   virtual void CleanUp() {}                     // to clear store
X//protected:
X//   BaseMatrix() : t(this), i(this) {}
X
X   friend GeneralMatrix;
X   friend Matrix;
X   friend nricMatrix;
X   friend RowVector;
X   friend ColumnVector;
X   friend SymmetricMatrix;
X   friend UpperTriangularMatrix;
X   friend LowerTriangularMatrix;
X   friend DiagonalMatrix;
X   friend CroutMatrix;
X   friend BandMatrix;
X   friend LowerBandMatrix;
X   friend UpperBandMatrix;
X   friend SymmetricBandMatrix;
X   friend AddedMatrix;
X   friend MultipliedMatrix;
X   friend SubtractedMatrix;
X   friend SolvedMatrix;
X   friend ShiftedMatrix;
X   friend ScaledMatrix;
X   friend TransposedMatrix;
X   friend NegatedMatrix;
X   friend InvertedMatrix;
X   friend RowedMatrix;
X   friend ColedMatrix;
X   friend DiagedMatrix;
X   friend MatedMatrix;
X   friend GetSubMatrix;
X   friend ConstMatrix;
X   friend ReturnMatrixX;
X   friend LinearEquationSolver;
X   NEW_DELETE(BaseMatrix)
X};
X
X
X/******************************* working classes **************************/
X
Xclass GeneralMatrix : public BaseMatrix         // declarable matrix types
X{
Xprotected:
X   int tag;                                     // shows whether can reuse
X   int nrows, ncols;                            // dimensions
X   int storage;                                 // total store required
X   Real* store;                                 // point to store (0=not set)
X   GeneralMatrix();                             // initialise with no store
X   GeneralMatrix(ArrayLengthSpecifier);         // constructor getting store
X   void Add(GeneralMatrix*, Real);              // sum of GM and Real
X   void Add(Real);                              // add Real to this
X   void Multiply(GeneralMatrix*, Real);         // product of GM and Real
X   void Multiply(Real);                         // multiply this by Real
X   void Negate(GeneralMatrix*);                 // change sign
X   void Negate();                               // change sign
X   void operator=(Real);                        // set matrix to constant
X   Real* GetStore();                            // get store or copy
X   GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
X                                                // temporarily access store
X   void GetMatrix(const GeneralMatrix*);        // used by = and initialise
X   void Eq(const BaseMatrix&, MatrixType);      // used by =
X   int search(const BaseMatrix*) const;
X   virtual GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
X   void CheckConversion(const BaseMatrix&);     // check conversion OK
X   void ReDimension(int, int, int);             // change dimensions
X//   int NrowsV() const;                          // get dimensions
X//   int NcolsV() const;                          // virtual version
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   MatrixType Type() const = 0;                 // type of a matrix
X   int Nrows() const { return nrows; }          // get dimensions
X   int Ncols() const { return ncols; }
X   int Storage() const { return storage; }
X   Real* Store() const { return store; }
X   virtual ~GeneralMatrix();                    // delete store if set
X   void tDelete();                              // delete if tag permits
X   Boolean reuse();                                // TRUE if tag allows reuse
X   void Protect() { tag=-1; }                   // can't delete or reuse
X   int Tag() const { return tag; }
X   Boolean IsZero() const;                         // test matrix has all zeros
X   void Release() { tag=1; }                    // del store after next use
X   void Release(int t) { tag=t; }               // del store after t accesses
X   void ReleaseAndDelete() { tag=0; }           // delete matrix after use
X   void operator<<(const Real*);                // assignment from an array
X   void operator<<(const BaseMatrix& X) { Eq(X,this->Type()); }
X                                                // = without checking type
X   void Inject(const GeneralMatrix&);           // copy stored els only
X   virtual GeneralMatrix* MakeSolver();         // for solving
X   virtual void Solver(MatrixRowCol&, const MatrixRowCol&) {}
X   virtual void GetRow(MatrixRowCol&) = 0;      // Get matrix row
X   virtual void RestoreRow(MatrixRowCol&) {}    // Restore matrix row
X   virtual void NextRow(MatrixRowCol&);         // Go to next row
X   virtual void GetCol(MatrixRowCol&) = 0;      // Get matrix col
X   virtual void RestoreCol(MatrixRowCol&) {}    // Restore matrix col
X   virtual void NextCol(MatrixRowCol&);         // Go to next col
X   Real SumSquare() const;
X   Real SumAbsoluteValue() const;
X   Real MaximumAbsoluteValue() const;
X   LogAndSign LogDeterminant() const;
X#ifndef TEMPS_DESTROYED_QUICKLY
X   ConstMatrix c() const;                       // to access constant matrices
X#else
X   ConstMatrix& c() const;                      // to access constant matrices
X#endif
X   void CheckStore() const;                     // check store is non-zero
X   virtual void SetParameters(const GeneralMatrix*) {}
X                                                // set parameters in GetMatrix
X#ifndef TEMPS_DESTROYED_QUICKLY
X   operator ReturnMatrixX() const;              // for building a ReturnMatrix
X#else
X   operator ReturnMatrixX&() const;             // for building a ReturnMatrix
X#endif
X   void CleanUp();                                // to clear store
X
X   friend Matrix;
X   friend nricMatrix;
X   friend SymmetricMatrix;
X   friend UpperTriangularMatrix;
X   friend LowerTriangularMatrix;
X   friend DiagonalMatrix;
X   friend CroutMatrix;
X   friend RowVector;
X   friend ColumnVector;
X   friend BandMatrix;
X   friend LowerBandMatrix;
X   friend UpperBandMatrix;
X   friend SymmetricBandMatrix;
X   friend BaseMatrix;
X   friend AddedMatrix;
X   friend MultipliedMatrix;
X   friend SubtractedMatrix;
X   friend SolvedMatrix;
X   friend ShiftedMatrix;
X   friend ScaledMatrix;
X   friend TransposedMatrix;
X   friend NegatedMatrix;
X   friend InvertedMatrix;
X   friend RowedMatrix;
X   friend ColedMatrix;
X   friend DiagedMatrix;
X   friend MatedMatrix;
X   friend GetSubMatrix;
X   friend ConstMatrix;
X   friend ReturnMatrixX;
X   friend LinearEquationSolver;
X   NEW_DELETE(GeneralMatrix)
X};
X
Xclass Matrix : public GeneralMatrix             // usual rectangular matrix
X{
Xpublic:
X   Matrix() {}
X   Matrix(int, int);                            // standard declaration
X   Matrix(const BaseMatrix&);                   // evaluate BaseMatrix
X   void operator=(const BaseMatrix&);
X   void operator=(Real f) { GeneralMatrix::operator=(f); }
X   MatrixType Type() const;
X   Real& operator()(int, int);                  // access element
X   Real& element(int, int);                     // access element
X   Matrix(const Matrix& gm) { GetMatrix(&gm); }
X#ifndef __ZTC__
X   Real operator()(int, int) const;             // access element
X#endif
X   GeneralMatrix* MakeSolver();
X   Real Trace() const;
X   void GetRow(MatrixRowCol&);
X   void GetCol(MatrixRowCol&);
X   void RestoreCol(MatrixRowCol&);
X   void NextRow(MatrixRowCol&);
X   void NextCol(MatrixRowCol&);
X   void ReDimension(int,int);                   // change dimensions
X   NEW_DELETE(Matrix)
X};
X
Xclass nricMatrix : public Matrix                // for use with Numerical
X                                                // Recipes in C
X{
X   Real** row_pointer;                          // points to rows
X   void MakeRowPointer();                       // build rowpointer
X   void DeleteRowPointer();
Xpublic:
X   nricMatrix() {}
X   nricMatrix(int m, int n)                     // standard declaration
X      :  Matrix(m,n) { MakeRowPointer(); }
X   nricMatrix(const BaseMatrix& bm)             // evaluate BaseMatrix
X      :  Matrix(bm) { MakeRowPointer(); }
X   void operator=(const BaseMatrix& bm)
X      { DeleteRowPointer(); Matrix::operator=(bm); MakeRowPointer(); }
X   void operator=(Real f) { GeneralMatrix::operator=(f); }
X   void operator<<(const BaseMatrix& X)
X      { DeleteRowPointer(); Eq(X,this->Type()); MakeRowPointer(); }
X   nricMatrix(const nricMatrix& gm) { GetMatrix(&gm); MakeRowPointer(); }
X   void ReDimension(int m, int n)               // change dimensions
X      { DeleteRowPointer(); Matrix::ReDimension(m,n); MakeRowPointer(); }
X   ~nricMatrix() { DeleteRowPointer(); }
X#ifndef __ZTC__
X   Real** nric() const { CheckStore(); return row_pointer-1; }
X#endif
X   void CleanUp();                                // to clear store
X   NEW_DELETE(nricMatrix)
X};
X
Xclass SymmetricMatrix : public GeneralMatrix
X{
Xpublic:
X   SymmetricMatrix() {}
X   SymmetricMatrix(ArrayLengthSpecifier);
X   SymmetricMatrix(const BaseMatrix&);
X   void operator=(const BaseMatrix&);
X   void operator=(Real f) { GeneralMatrix::operator=(f); }
X   Real& operator()(int, int);                  // access element
X   Real& element(int, int);                     // access element
X   MatrixType Type() const;
X   SymmetricMatrix(const SymmetricMatrix& gm) { GetMatrix(&gm); }
X#ifndef __ZTC__
X   Real operator()(int, int) const;             // access element
X#endif
X   Real SumSquare() const;
X   Real SumAbsoluteValue() const;
X   Real Trace() const;
X   void GetRow(MatrixRowCol&);
X   void GetCol(MatrixRowCol&);
X   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
X   void ReDimension(int);                       // change dimensions
X   NEW_DELETE(SymmetricMatrix)
X};
X
Xclass UpperTriangularMatrix : public GeneralMatrix
X{
Xpublic:
X   UpperTriangularMatrix() {}
X   UpperTriangularMatrix(ArrayLengthSpecifier);
X   void operator=(const BaseMatrix&);
X   UpperTriangularMatrix(const BaseMatrix&);
X   UpperTriangularMatrix(const UpperTriangularMatrix& gm) { GetMatrix(&gm); }
X#ifndef __ZTC__
X   Real operator()(int, int) const;             // access element
X#endif
X   void operator=(Real f) { GeneralMatrix::operator=(f); }
X   Real& operator()(int, int);                  // access element
X   Real& element(int, int);                     // access element
X   MatrixType Type() const;
X   GeneralMatrix* MakeSolver() { return this; } // for solving
X   void Solver(MatrixRowCol&, const MatrixRowCol&);
X   LogAndSign LogDeterminant() const;
X   Real Trace() const;
X   void GetRow(MatrixRowCol&);
X   void GetCol(MatrixRowCol&);
X   void RestoreCol(MatrixRowCol&);
X   void NextRow(MatrixRowCol&);
X   void ReDimension(int);                       // change dimensions
X   NEW_DELETE(UpperTriangularMatrix)
X};
X
Xclass LowerTriangularMatrix : public GeneralMatrix
X{
Xpublic:
X   LowerTriangularMatrix() {}
X   LowerTriangularMatrix(ArrayLengthSpecifier);
X   LowerTriangularMatrix(const LowerTriangularMatrix& gm) { GetMatrix(&gm); }
X#ifndef __ZTC__
X   Real operator()(int, int) const;             // access element
X#endif
X   LowerTriangularMatrix(const BaseMatrix& M);
X   void operator=(const BaseMatrix&);
X   void operator=(Real f) { GeneralMatrix::operator=(f); }
X   Real& operator()(int, int);                  // access element
X   Real& element(int, int);                     // access element
X   MatrixType Type() const;
X   GeneralMatrix* MakeSolver() { return this; } // for solving
X   void Solver(MatrixRowCol&, const MatrixRowCol&);
X   LogAndSign LogDeterminant() const;
X   Real Trace() const;
X   void GetRow(MatrixRowCol&);
X   void GetCol(MatrixRowCol&);
X   void RestoreCol(MatrixRowCol&);
X   void NextRow(MatrixRowCol&);
X   void ReDimension(int);                       // change dimensions
X   NEW_DELETE(LowerTriangularMatrix)
X};
X
Xclass DiagonalMatrix : public GeneralMatrix
X{
Xpublic:
X   DiagonalMatrix() {}
X   DiagonalMatrix(ArrayLengthSpecifier);
X   DiagonalMatrix(const BaseMatrix&);
X   DiagonalMatrix(const DiagonalMatrix& gm) { GetMatrix(&gm); }
X#ifndef __ZTC__
X   Real operator()(int, int) const;             // access element
X   Real operator()(int) const;
X#endif
X   void operator=(const BaseMatrix&);
X   void operator=(Real f) { GeneralMatrix::operator=(f); }
X   Real& operator()(int, int);                  // access element
X   Real& operator()(int);                       // access element
X   Real& element(int, int);                     // access element
X   Real& element(int);                          // access element
X   MatrixType Type() const;
X
X   LogAndSign LogDeterminant() const;
X   Real Trace() const;
X   void GetRow(MatrixRowCol&);
X   void GetCol(MatrixRowCol&);
X   void NextRow(MatrixRowCol&);
X   void NextCol(MatrixRowCol&);
X   GeneralMatrix* MakeSolver() { return this; } // for solving
X   void Solver(MatrixRowCol&, const MatrixRowCol&);
X   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
X   void ReDimension(int);                       // change dimensions
X#ifndef __ZTC__
X   Real* nric() const
X      { CheckStore(); return store-1; }         // for use by NRIC
X#endif
X   MatrixBandWidth BandWidth() const;
X   NEW_DELETE(DiagonalMatrix)
X};
X
Xclass RowVector : public Matrix
X{
Xpublic:
X   RowVector() {}
X   RowVector(ArrayLengthSpecifier n) : Matrix(1,n.Value()) {}
X   RowVector(const BaseMatrix&);
X   RowVector(const RowVector& gm) { GetMatrix(&gm); }
X#ifndef __ZTC__
X   Real operator()(int) const;                  // access element
X#endif
X   void operator=(const BaseMatrix&);
X   void operator=(Real f) { GeneralMatrix::operator=(f); }
X   Real& operator()(int);                       // access element
X   Real& element(int);                          // access element
X   MatrixType Type() const;
X   void GetCol(MatrixRowCol&);
X   void NextCol(MatrixRowCol&);
X   void RestoreCol(MatrixRowCol&);
X   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
X   void ReDimension(int);                       // change dimensions
X#ifndef __ZTC__
X   Real* nric() const
X      { CheckStore(); return store-1; }         // for use by NRIC
X#endif
X   void CleanUp();                                // to clear store
X   NEW_DELETE(RowVector)
X};
X
Xclass ColumnVector : public Matrix
X{
Xpublic:
X   ColumnVector() {}
X   ColumnVector(ArrayLengthSpecifier n) : Matrix(n.Value(),1) {}
X   ColumnVector(const BaseMatrix&);
X   ColumnVector(const ColumnVector& gm) { GetMatrix(&gm); }
X#ifndef __ZTC__
X   Real operator()(int) const;                  // access element
X#endif
X   void operator=(const BaseMatrix&);
X   void operator=(Real f) { GeneralMatrix::operator=(f); }
X   Real& operator()(int);                       // access element
X   Real& element(int);                          // access element
X   MatrixType Type() const;
X   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
X   void ReDimension(int);                       // change dimensions
X#ifndef __ZTC__
X   Real* nric() const
X      { CheckStore(); return store-1; }         // for use by NRIC
X#endif
X   void CleanUp();                                // to clear store
X   NEW_DELETE(ColumnVector)
X};
X
Xclass CroutMatrix : public GeneralMatrix        // for LU decomposition
X{
X   int* indx;
X   Boolean d;
X   Boolean sing;
X   void ludcmp();
Xpublic:
X   CroutMatrix(const BaseMatrix&);
X   MatrixType Type() const;
X   void lubksb(Real*, int=0);
X   ~CroutMatrix();
X   GeneralMatrix* MakeSolver() { return this; } // for solving
X   LogAndSign LogDeterminant() const;
X   void Solver(MatrixRowCol&, const MatrixRowCol&);
X   void GetRow(MatrixRowCol&);
X   void GetCol(MatrixRowCol&);
X   void operator=(const BaseMatrix&);
X   void CleanUp();                                // to clear store
X   NEW_DELETE(CroutMatrix)
X};
X
X/******************************* band matrices ***************************/
X
Xclass BandMatrix : public GeneralMatrix         // band matrix
X{
Xprotected:
X   void CornerClear() const;                    // set unused elements to zero
Xpublic:
X   int lower, upper;                            // band widths
X   BandMatrix() { lower=0; upper=0; CornerClear(); }
X   BandMatrix(int n,int lb,int ub) { ReDimension(n,lb,ub); CornerClear(); }
X                                                // standard declaration
X   BandMatrix(const BaseMatrix&);               // evaluate BaseMatrix
X   void operator=(const BaseMatrix&);
X   void operator=(Real f) { GeneralMatrix::operator=(f); }
X   MatrixType Type() const;
X   Real& operator()(int, int);                  // access element
X   Real& element(int, int);                     // access element
X   BandMatrix(const BandMatrix& gm) { GetMatrix(&gm); }
X#ifndef __ZTC__
X   Real operator()(int, int) const;             // access element
X#endif
X   LogAndSign LogDeterminant() const;
X   GeneralMatrix* MakeSolver();
X   Real Trace() const;
X   Real SumSquare() const { CornerClear(); return GeneralMatrix::SumSquare(); }
X   Real SumAbsoluteValue() const
X      { CornerClear(); return GeneralMatrix::SumAbsoluteValue(); }
X   Real MaximumAbsoluteValue() const
X      { CornerClear(); return GeneralMatrix::MaximumAbsoluteValue(); }
X   void GetRow(MatrixRowCol&);
X   void GetCol(MatrixRowCol&);
X   void RestoreCol(MatrixRowCol&);
X   void NextRow(MatrixRowCol&);
X   void ReDimension(int, int, int);             // change dimensions
X   MatrixBandWidth BandWidth() const;
X   void SetParameters(const GeneralMatrix*);
X   NEW_DELETE(BandMatrix)
X};
X
Xclass UpperBandMatrix : public BandMatrix       // upper band matrix
X{
Xpublic:
X   UpperBandMatrix() {}
X   UpperBandMatrix(int n, int ubw)              // standard declaration
X      : BandMatrix(n, 0, ubw) {}
X   UpperBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
X   void operator=(const BaseMatrix&);
X   void operator=(Real f) { GeneralMatrix::operator=(f); }
X   MatrixType Type() const;
X   UpperBandMatrix(const UpperBandMatrix& gm) { GetMatrix(&gm); }
X   GeneralMatrix* MakeSolver() { return this; }
X   void Solver(MatrixRowCol&, const MatrixRowCol&);
X   LogAndSign LogDeterminant() const;
X   void ReDimension(int n,int ubw)              // change dimensions
X      { BandMatrix::ReDimension(n,0,ubw); }
X   Real& operator()(int, int);
X   Real& element(int, int);
X   NEW_DELETE(UpperBandMatrix)
X};
X
Xclass LowerBandMatrix : public BandMatrix       // upper band matrix
X{
Xpublic:
X   LowerBandMatrix() {}
X   LowerBandMatrix(int n, int lbw)              // standard declaration
X      : BandMatrix(n, lbw, 0) {}
X   LowerBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
X   void operator=(const BaseMatrix&);
X   void operator=(Real f) { GeneralMatrix::operator=(f); }
X   MatrixType Type() const;
X   LowerBandMatrix(const LowerBandMatrix& gm) { GetMatrix(&gm); }
X   GeneralMatrix* MakeSolver() { return this; }
X   void Solver(MatrixRowCol&, const MatrixRowCol&);
X   LogAndSign LogDeterminant() const;
X   void ReDimension(int n,int lbw)             // change dimensions
X      { BandMatrix::ReDimension(n,lbw,0); }
X   Real& operator()(int, int);
X   Real& element(int, int);
X   NEW_DELETE(LowerBandMatrix)
X};
X
Xclass SymmetricBandMatrix : public GeneralMatrix
X{
X   void CornerClear() const;                    // set unused elements to zero
Xpublic:
X   int lower;                                   // lower band width
X   SymmetricBandMatrix() { lower=0; CornerClear(); }
X   SymmetricBandMatrix(int n, int lb) { ReDimension(n,lb); CornerClear(); }
X   SymmetricBandMatrix(const BaseMatrix&);
X   void operator=(const BaseMatrix&);
X   void operator=(Real f) { GeneralMatrix::operator=(f); }
X   Real& operator()(int, int);                  // access element
X   Real& element(int, int);                     // access element
X   MatrixType Type() const;
X   SymmetricBandMatrix(const SymmetricBandMatrix& gm) { GetMatrix(&gm); }
X#ifndef __ZTC__
X   Real operator()(int, int) const;             // access element
X#endif
X   GeneralMatrix* MakeSolver();
X   Real SumSquare() const;
X   Real SumAbsoluteValue() const;
X   Real MaximumAbsoluteValue() const
X      { CornerClear(); return GeneralMatrix::MaximumAbsoluteValue(); }
X   Real Trace() const;
X   LogAndSign LogDeterminant() const;
X   void GetRow(MatrixRowCol&);
X   void GetCol(MatrixRowCol&);
X   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
X   void ReDimension(int,int);                       // change dimensions
X   MatrixBandWidth BandWidth() const;
X   void SetParameters(const GeneralMatrix*);
X   NEW_DELETE(SymmetricBandMatrix)
X};
X
Xclass BandLUMatrix : public GeneralMatrix
X// for LU decomposition of band matrix
X{
X   int* indx;
X   Boolean d;
X   Boolean sing;                                   // TRUE if singular
X   Real* store2;
X   int storage2;
X   void ludcmp();
X   int m1,m2;                                   // lower and upper
Xpublic:
X   BandLUMatrix(const BaseMatrix&);
X   MatrixType Type() const;
X   void lubksb(Real*, int=0);
X   ~BandLUMatrix();
X   GeneralMatrix* MakeSolver() { return this; } // for solving
X   LogAndSign LogDeterminant() const;
X   void Solver(MatrixRowCol&, const MatrixRowCol&);
X   void GetRow(MatrixRowCol&);
X   void GetCol(MatrixRowCol&);
X   void operator=(const BaseMatrix&);
X   NEW_DELETE(BandLUMatrix)
X   void CleanUp();                                // to clear store
X};
X
X
X/***************************** temporary classes *************************/
X
Xclass MultipliedMatrix : public BaseMatrix
X{
Xprotected:
X   union { const BaseMatrix* bm1; GeneralMatrix* gm1; };
X                          // pointers to summands
X   union { const BaseMatrix* bm2; GeneralMatrix* gm2; };
X   MultipliedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
X      : bm1(bm1x),bm2(bm2x) {}
X   int search(const BaseMatrix*) const;
X//   int NrowsV() const;
X//   int NcolsV() const;
X   MatrixType Type() const;
X   friend BaseMatrix;
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   MatrixBandWidth BandWidth() const;
X   NEW_DELETE(MultipliedMatrix)
X};
X
Xclass AddedMatrix : public MultipliedMatrix
X{
Xprotected:
X   AddedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
X      : MultipliedMatrix(bm1x,bm2x) {}
X
Xprivate:
X   MatrixType Type() const;
X   friend BaseMatrix;
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   MatrixBandWidth BandWidth() const;
X#ifdef GXX
X   void SelectVersion(MatrixType, int&, int&) const;
X#else
X   void SelectVersion(MatrixType, Boolean&, Boolean&) const;
X#endif
X   NEW_DELETE(AddedMatrix)
X};
X
Xclass SolvedMatrix : public MultipliedMatrix
X{
X   SolvedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
X      : MultipliedMatrix(bm1x,bm2x) {}
X   MatrixType Type() const;
X   friend BaseMatrix;
X   friend InvertedMatrix;                        // for operator*
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   MatrixBandWidth BandWidth() const;
X   NEW_DELETE(SolvedMatrix)
X};
X
Xclass SubtractedMatrix : public AddedMatrix
X{
X   SubtractedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
X      : AddedMatrix(bm1x,bm2x) {}
X   MatrixType Type() const;
X   friend BaseMatrix;
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   NEW_DELETE(SubtractedMatrix)
X};
X
Xclass ShiftedMatrix : public BaseMatrix
X{
Xprotected:
X   Real f;
X   union { const BaseMatrix* bm; GeneralMatrix* gm; };
X   ShiftedMatrix(const BaseMatrix* bmx, Real fx) : bm(bmx),f(fx) {}
X   int search(const BaseMatrix*) const;
X//   int NrowsV() const;
X//   int NcolsV() const;
Xprivate:
X   MatrixType Type() const;
X   friend BaseMatrix;
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   NEW_DELETE(ShiftedMatrix)
X};
X
Xclass ScaledMatrix : public ShiftedMatrix
X{
X   ScaledMatrix(const BaseMatrix* bmx, Real fx) : ShiftedMatrix(bmx,fx) {}
X   MatrixType Type() const;
X   friend BaseMatrix;
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   MatrixBandWidth BandWidth() const;
X   NEW_DELETE(ScaledMatrix)
X};
X
Xclass NegatedMatrix : public BaseMatrix
X{
Xprotected:
X   union { const BaseMatrix* bm; GeneralMatrix* gm; };
X   NegatedMatrix(const BaseMatrix* bmx) : bm(bmx) {}
X   int search(const BaseMatrix*) const;
X//   int NrowsV() const;
X//   int NcolsV() const;
Xprivate:
X   MatrixType Type() const;
X   friend BaseMatrix;
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   MatrixBandWidth BandWidth() const;
X   NEW_DELETE(NegatedMatrix)
X};
X
Xclass TransposedMatrix : public NegatedMatrix
X{
X   TransposedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
X//   int NrowsV() const;
X//   int NcolsV() const;
X   MatrixType Type() const;
X   friend BaseMatrix;
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   MatrixBandWidth BandWidth() const;
X   NEW_DELETE(TransposedMatrix)
X};
X
Xclass InvertedMatrix : public NegatedMatrix
X{
X   InvertedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
X   MatrixType Type() const;
Xpublic:
X#ifndef TEMPS_DESTROYED_QUICKLY
X   SolvedMatrix operator*(const BaseMatrix&) const;       // inverse(A) * B
X#else
X   SolvedMatrix& operator*(const BaseMatrix&) const;       // inverse(A) * B
X#endif
X   friend BaseMatrix;
X   GeneralMatrix* Evaluate(MatrixType);
X   MatrixBandWidth BandWidth() const;
X   NEW_DELETE(InvertedMatrix)
X};
X
Xclass RowedMatrix : public NegatedMatrix
X{
X   MatrixType Type() const;
X//   int NrowsV() const;
X//   int NcolsV() const;
X   RowedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
X   friend BaseMatrix;
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   MatrixBandWidth BandWidth() const;
X   NEW_DELETE(RowedMatrix)
X};
X
Xclass ColedMatrix : public NegatedMatrix
X{
X   MatrixType Type() const;
X//   int NrowsV() const;
X//   int NcolsV() const;
X   ColedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
X   friend BaseMatrix;
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   MatrixBandWidth BandWidth() const;
X   NEW_DELETE(ColedMatrix)
X};
X
Xclass DiagedMatrix : public NegatedMatrix
X{
X   MatrixType Type() const;
X//   int NrowsV() const;
X//   int NcolsV() const;
X   DiagedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
X   friend BaseMatrix;
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   MatrixBandWidth BandWidth() const;
X   NEW_DELETE(DiagedMatrix)
X};
X
Xclass MatedMatrix : public NegatedMatrix
X{
X   int nr, nc;
X   MatrixType Type() const;
X//   int NrowsV() const;
X//   int NcolsV() const;
X   MatedMatrix(const BaseMatrix* bmx, int nrx, int ncx)
X      : NegatedMatrix(bmx), nr(nrx), nc(ncx) {}
X   friend BaseMatrix;
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   MatrixBandWidth BandWidth() const;
X   NEW_DELETE(MatedMatrix)
X};
X
Xclass ConstMatrix : public BaseMatrix
X{
X   const GeneralMatrix* cgm;
X   MatrixType Type() const;
X//   int NrowsV() const;
X//   int NcolsV() const;
X   int search(const BaseMatrix*) const;
X   ConstMatrix(const GeneralMatrix* cgmx) : cgm(cgmx) {}
X   friend BaseMatrix;
X   friend GeneralMatrix;
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   MatrixBandWidth BandWidth() const;
X   NEW_DELETE(ConstMatrix)
X};
X
Xclass ReturnMatrixX : public BaseMatrix    // for matrix return
X{
X   GeneralMatrix* gm;
X   MatrixType Type() const;
X//   int NrowsV() const;
X//   int NcolsV() const;
X   int search(const BaseMatrix*) const;
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   friend BaseMatrix;
X#ifdef TEMPS_DESTROYED_QUICKLY
X   ReturnMatrixX(const ReturnMatrixX& tm);
X#else
X   ReturnMatrixX(const ReturnMatrixX& tm) : gm(tm.gm) {}
X#endif
X   ReturnMatrixX(const GeneralMatrix* gmx) : gm((GeneralMatrix*&)gmx) {}
X//   ReturnMatrixX(GeneralMatrix&);
X   MatrixBandWidth BandWidth() const;
X   NEW_DELETE(ReturnMatrixX)
X};
X
X
X/**************************** submatrices ******************************/
X
Xclass GetSubMatrix : public NegatedMatrix
X{
X   int row_skip;
X   int row_number;
X   int col_skip;
X   int col_number;
X   MatrixType mt;
X
X   GetSubMatrix
X      (const BaseMatrix* bmx, int rs, int rn, int cs, int cn, MatrixType mtx)
X      : NegatedMatrix(bmx),
X      row_skip(rs), row_number(rn), col_skip(cs), col_number(cn), mt(mtx) {}
X   GetSubMatrix(const GetSubMatrix& g)
X      : NegatedMatrix(g.bm), row_skip(g.row_skip), row_number(g.row_number),
X      col_skip(g.col_skip), col_number(g.col_number), mt(g.mt) {}
X   void SetUpLHS();
X   MatrixType Type() const;
X//   int NrowsV() const;
X//   int NcolsV() const;
X   friend BaseMatrix;
Xpublic:
X   GeneralMatrix* Evaluate(MatrixType);
X   void operator=(const BaseMatrix&);
X   void operator<<(const BaseMatrix&);
X   void operator<<(const Real*);                // copy from array
X   void operator<<(Real);                       // copy from constant
X   void Inject(const GeneralMatrix&);           // copy stored els only
X   MatrixBandWidth BandWidth() const;
X   NEW_DELETE(GetSubMatrix)
X};
X
X
X/***************************** exceptions ********************************/
X
Xclass MatrixDetails
X{
X   MatrixType type;
X   int nrows;
X   int ncols;
X   int ubw;
X   int lbw;
Xpublic:
X   MatrixDetails(const GeneralMatrix& A);
X   void PrintOut();
X};
X
X
X
Xclass SpaceException : public Exception
X{
Xpublic:
X   static long st_type() { return 2; }
X   long type() const { return 2; }
X   static int action;
X   SpaceException();
X   static void SetAction(int a) { action=a; }
X};
X
X
Xclass MatrixException : public Exception
X{
Xpublic:
X   static long st_type() { return 3; }
X   long type() const { return 3; }
X   MatrixException(int);
X   MatrixException(int, const GeneralMatrix&);
X   MatrixException(int, const GeneralMatrix&, const GeneralMatrix&);
X};
X
Xclass DataException : public MatrixException
X{
Xpublic:
X   static long st_type() { return 3*53; }
X   long type() const { return 3*53; }
X   static int action;
X   DataException(const GeneralMatrix& A);
X   static void SetAction(int a) { action=a; }
X};
X
Xclass SingularException : public DataException
X{
Xpublic:
X   static long st_type() { return 3*53*109; }
X   long type() const { return 3*53*109; }
X   SingularException(const GeneralMatrix& A);
X};
X
Xclass NPDException : public DataException     // Not positive definite
X{
Xpublic:
X   static long st_type() { return 3*53*113; }
X   long type() const { return 3*53*113; }
X   NPDException(const GeneralMatrix&);
X};
X
Xclass ConvergenceException : public MatrixException
X{
Xpublic:
X   static long st_type() { return 3*59; }
X   long type() const { return 3*59; }
X   static int action;
X   ConvergenceException(const GeneralMatrix& A);
X   static void SetAction(int a) { action=a; }
X};
X
Xclass ProgramException : public MatrixException
X{
Xpublic:
X   static long st_type() { return 3*61; }
X   long type() const { return 3*61; }
X   static int action;
X   ProgramException(char* c);
X   ProgramException(char* c, const GeneralMatrix&);
X   ProgramException(char* c, const GeneralMatrix&, const GeneralMatrix&);
X   ProgramException();
X   ProgramException(const GeneralMatrix&);
X   static void SetAction(int a) { action=a; }
X};
X
Xclass IndexException : public ProgramException
X{
Xpublic:
X   static long st_type() { return 3*61*101; }
X   long type() const { return 3*61*101; }
X   IndexException(int i, const GeneralMatrix& A);
X   IndexException(int i, int j, const GeneralMatrix& A);
X   // next two are for access via element function
X   IndexException(int i, const GeneralMatrix& A, Boolean);
X   IndexException(int i, int j, const GeneralMatrix& A, Boolean);
X};
X
Xclass VectorException : public ProgramException    // can't convert to vector
X{
Xpublic:
X   static long st_type() { return 3*61*107; }
X   long type() const { return 3*61*107; }
X   VectorException();
X   VectorException(const GeneralMatrix& A);
X};
X
Xclass NotSquareException : public ProgramException
X{
Xpublic:
X   static long st_type() { return 3*61*109; }
X   long type() const { return 3*61*109; }
X   NotSquareException(const GeneralMatrix& A);
X};
X
Xclass SubMatrixDimensionException : public ProgramException
X{
Xpublic:
X   static long st_type() { return 3*61*113; }
X   long type() const { return 3*61*113; }
X   SubMatrixDimensionException();
X};
X
Xclass IncompatibleDimensionsException : public ProgramException
X{
Xpublic:
X   static long st_type() { return 3*61*127; }
X   long type() const { return 3*61*127; }
X   IncompatibleDimensionsException();
X};
X
Xclass NotDefinedException : public ProgramException
X{
Xpublic:
X   static long st_type() { return 3*61*131; }
X   long type() const { return 3*61*131; }
X   NotDefinedException(char* op, char* matrix);
X};
X
Xclass CannotBuildException : public ProgramException
X{
Xpublic:
X   static long st_type() { return 3*61*137; }
X   long type() const { return 3*61*137; }
X   CannotBuildException(char* matrix);
X};
X
X
Xclass InternalException : public MatrixException
X{
Xpublic:
X   static long st_type() { return 3*67; }
X   long type() const { return 3*67; }
X   static int action;
X   InternalException(char* c);
X   static void SetAction(int a) { action=a; }
X};
X
X
X/***************************** functions ***********************************/
X
X
Xinline LogAndSign LogDeterminant(const BaseMatrix& B)
X   { return B.LogDeterminant(); }
Xinline Real SumSquare(const BaseMatrix& B) { return B.SumSquare(); }
Xinline Real Trace(const BaseMatrix& B) { return B.Trace(); }
Xinline Real SumAbsoluteValue(const BaseMatrix& B)
X   { return B.SumAbsoluteValue(); }
Xinline Real MaximumAbsoluteValue(const BaseMatrix& B)
X   { return B.MaximumAbsoluteValue(); }
Xinline Real Norm1(const BaseMatrix& B) { return B.Norm1(); }
Xinline Real Norm1(RowVector& RV) { return RV.MaximumAbsoluteValue(); }
Xinline Real NormInfinity(const BaseMatrix& B) { return B.NormInfinity(); }
Xinline Real NormInfinity(ColumnVector& CV)
X   { return CV.MaximumAbsoluteValue(); }
X
X#endif
END_OF_FILE
if test 42428 -ne `wc -c <'newmat.h'`; then
    echo shar: \"'newmat.h'\" unpacked with wrong size!
fi
# end of 'newmat.h'
fi
if test -f 'newmatc.txt' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'newmatc.txt'\"
else
echo shar: Extracting \"'newmatc.txt'\" \(1887 characters\)
sed "s/^X//" >'newmatc.txt' <<'END_OF_FILE'
X//$$ newmatc.txt                                      Testing
X
X                    Testing newmat on your compiler
X                    ===============================
X
XThere are a series of files of the form tmt?.cxx included with some
Xversions of the package. These are used to check that the package is
Xperforming correctly. They consist of a large number of matrix formulae
Xall of which evaluate to zero (except the first one which is used to
Xcheck that we are detecting non-zero matrices). The printout should
Xstate that it has found one non-zero matrix.
X
XThere is an optional #define DO_FREE_CHECK in include.h. If this is
Xactivated the program will check that the new-s and delete-s are
Xbalanced. There should be no unbalanced new-s or delete-s. You need the
XANSI version of the preprocessor to run DO_FREE_CHECK.
X
XThe program also allocates and deletes a large block and small block of
Xmemory before it starts the main testing and then at the end of the
Xtest. It then checks that the blocks of memory were allocated in the
Xsame place. If not then one suspects that there has been a memory leak.
Xi.e. a piece of memory has been allocated and not deleted.
X
XThis is not completely foolproof. Programs may allocate extra print
Xbuffers while the program is running. I have tried to overcome this by
Xdoing a print before I allocate the first memory block. Programs may
Xallocate memory for different sized items in different places, or might
Xnot allocate items consecutively. Or they might mix the items with memory
Xblocks from other programs. Nevertheless, I seem to get consistent
Xanswers from most of the compilers I am working with, so I think this is
Xa worthwhile test.
X
XGnu 2.2 and Zortech do not pass this test. There may be good reasons for
Xthis, but I think it would be a good idea if the authors of these
Xpackages made sure they knew what was happening.
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
END_OF_FILE
if test 1887 -ne `wc -c <'newmatc.txt'`; then
    echo shar: \"'newmatc.txt'\" unpacked with wrong size!
fi
# end of 'newmatc.txt'
fi
if test -f 'tmt.cxx' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'tmt.cxx'\"
else
echo shar: Extracting \"'tmt.cxx'\" \(4634 characters\)
sed "s/^X//" >'tmt.cxx' <<'END_OF_FILE'
X#define WANT_STREAM
X
X#include "include.h"
X
X#include "newmat.h"
X
X/**************************** test program ******************************/
X
Xclass PrintCounter
X{
X   int count;
X   char* s;
Xpublic:
X   ~PrintCounter();
X   PrintCounter(char * sx) : count(0), s(sx) {}
X   void operator++() { count++; }
X};
X
X   PrintCounter PCZ("Number of non-zero matrices = ");
X   PrintCounter PCN("Number of matrices tested   = ");
X
XPrintCounter::~PrintCounter()
X{ cout << s << count << "\n"; }
X
X
Xvoid Print(const Matrix& X)
X{
X   ++PCN;
X   cout << "\nMatrix type: " << (char*)X.Type() << " (";
X   cout << X.Nrows() << ", ";
X   cout << X.Ncols() << ")\n\n";
X   if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
X   int nr=X.Nrows(); int nc=X.Ncols();
X   for (int i=1; i<=nr; i++)
X   {
X      for (int j=1; j<=nc; j++)  cout << X(i,j) << "\t";
X      cout << "\n";
X   }
X   cout << flush; ++PCZ;
X}
X
Xvoid Print(const UpperTriangularMatrix& X)
X{
X   ++PCN;
X   cout << "\nMatrix type: " << (char*)X.Type() << " (";
X   cout << X.Nrows() << ", ";
X   cout << X.Ncols() << ")\n\n";
X   if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
X   int nr=X.Nrows(); int nc=X.Ncols();
X   for (int i=1; i<=nr; i++)
X   {
X      for (int j=1; j<i; j++) cout << "\t";
X      for (j=i; j<=nc; j++)  cout << X(i,j) << "\t";
X      cout << "\n";
X   }
X   cout << flush; ++PCZ;
X}
X
Xvoid Print(const DiagonalMatrix& X)
X{
X   ++PCN;
X   cout << "\nMatrix type: " << (char*)X.Type() << " (";
X   cout << X.Nrows() << ", ";
X   cout << X.Ncols() << ")\n\n";
X   if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
X   int nr=X.Nrows(); int nc=X.Ncols();
X   for (int i=1; i<=nr; i++)
X   {
X      for (int j=1; j<i; j++) cout << "\t";
X      if (i<=nc) cout << X(i,i) << "\t";
X      cout << "\n";
X   }
X   cout << flush; ++PCZ;
X}
X
Xvoid Print(const SymmetricMatrix& X)
X{
X   ++PCN;
X   cout << "\nMatrix type: " << (char*)X.Type() << " (";
X   cout << X.Nrows() << ", ";
X   cout << X.Ncols() << ")\n\n";
X   if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
X   int nr=X.Nrows(); int nc=X.Ncols();
X   for (int i=1; i<=nr; i++)
X   {
X      for (int j=1; j<i; j++) cout << X(j,i) << "\t";
X      for (j=i; j<=nc; j++)  cout << X(i,j) << "\t";
X      cout << "\n";
X   }
X   cout << flush; ++PCZ;
X}
X
Xvoid Print(const LowerTriangularMatrix& X)
X{
X   ++PCN;
X   cout << "\nMatrix type: " << (char*)X.Type() << " (";
X   cout << X.Nrows() << ", ";
X   cout << X.Ncols() << ")\n\n";
X   if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
X   int nr=X.Nrows();
X   for (int i=1; i<=nr; i++)
X   {
X      for (int j=1; j<=i; j++) cout << X(i,j) << "\t";
X      cout << "\n";
X   }
X   cout << flush; ++PCZ;
X}
X
X
Xvoid Clean(Matrix& A, Real c)
X{
X   int nr = A.Nrows(); int nc = A.Ncols();
X   for (int i=1; i<=nr; i++)
X   {
X      for ( int j=1; j<=nc; j++)
X      { Real a = A(i,j); if ((a < c) && (a > -c)) A(i,j) = 0.0; }
X   }
X}
X
Xvoid Clean(DiagonalMatrix& A, Real c)
X{
X   int nr = A.Nrows();
X   for (int i=1; i<=nr; i++)
X   { Real a = A(i,i); if ((a < c) && (a > -c)) A(i,i) = 0.0; }
X}
X
X
Xvoid trymat1(); void trymat2(); void trymat3();
Xvoid trymat4(); void trymat5(); void trymat6();
Xvoid trymat7(); void trymat8(); void trymat9();
Xvoid trymata(); void trymatb(); void trymatc();
Xvoid trymatd(); void trymate(); void trymatf();
Xvoid trymatg(); void trymath(); void trymati();
X
Xmain()
X{
X   Real* s1; Real* s2; Real* s3; Real* s4;
X   cout << "\nBegin test\n";   // Forces cout to allocate memory at beginning
X   { Matrix A1(25,200); s1 = A1.Store(); }
X   { Matrix A1(1,1); s3 = A1.Store(); }
X   {
X      Tracer et("Matrix test program");
X
X      Matrix A(25,150);
X      {
X     int i;
X     RowVector A(8);
X     for (i=1;i<=7;i++) A(i)=0.0; A(8)=1.0;
X     Print(A);
X      }
X      cout << "\n";
X
X      TestTypeAdd(); TestTypeMult(); TestTypeOrder();
X
X      trymat1();
X      trymat2();
X      trymat3();
X      trymat4();
X      trymat5();
X      trymat6();
X      trymat7();
X      trymat8();
X      trymat9();
X      trymata();
X      trymatb();
X      trymatc();
X      trymatd();
X      trymate();
X      trymatf();
X      trymatg();
X      trymath();
X      trymati();
X   }
X   { Matrix A1(25,200); s2 = A1.Store(); }
X   cout << "\nChecking for lost memory: "
X      << (unsigned long)s1 << " " << (unsigned long)s2 << " ";
X   if (s1 != s2) cout << " - error\n"; else cout << " - ok\n\n";
X   { Matrix A1(1,1); s4 = A1.Store(); }
X   cout << "\nChecking for lost memory: "
X      << (unsigned long)s3 << " " << (unsigned long)s4 << " ";
X   if (s3 != s4) cout << " - error\n"; else cout << " - ok\n\n";
X#ifdef DO_FREE_CHECK
X   FreeCheck::Status();
X#endif
X   return 0;
X}
X
END_OF_FILE
if test 4634 -ne `wc -c <'tmt.cxx'`; then
    echo shar: \"'tmt.cxx'\" unpacked with wrong size!
fi
# end of 'tmt.cxx'
fi
if test -f 'tmt.mak' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'tmt.mak'\"
else
echo shar: Extracting \"'tmt.mak'\" \(2558 characters\)
sed "s/^X//" >'tmt.mak' <<'END_OF_FILE'
XOBJ = fft.o evalue.o submat.o cholesky.o hholder.o sort.o newmatrm.o     \
X  jacobi.o tmtf.o svd.o tmte.o tmtd.o newmat8.o tmtc.o tmtb.o            \
X  newmat7.o newmat6.o newmat5.o newmat3.o newmat4.o newmat2.o newmat1.o  \
X  tmt.o tmt1.o tmt2.o tmt3.o tmt4.o tmt5.o tmt6.o tmt7.o tmt8.o          \
X  tmt9.o tmta.o tmtg.o tmth.o tmti.o bandmat.o except.o newmatex.o
X
Xtmt:          $(OBJ)
X          g++ -o $@ $(OBJ) -lm
X
X%.o:          %.cxx
X          g++ -c $*.cxx
X
Xnewmatxx:     include.h newmat.h boolean.h except.h
X          rm -f newmatxx
X          echo "main .h files uptodate?" > newmatxx
X
Xexcept.o:     except.h except.cxx
X
Xnewmatex.o:   newmatxx newmatex.cxx
X
Xexample.o:    newmatxx newmatap.h example.cxx
X
Xcholesky.o:   newmatxx cholesky.cxx
X
Xevalue.o:     newmatxx newmatrm.h precisio.h evalue.cxx
X
Xfft.o:        newmatxx newmatap.h fft.cxx
X
Xhholder.o:    newmatxx newmatap.h hholder.cxx
X
Xjacobi.o:     newmatxx precisio.h newmatrm.h jacobi.cxx
X
Xbandmat.o:    newmatxx newmatrc.h controlw.h bandmat.cxx
X
Xnewmat1.o:    newmatxx newmat1.cxx
X
Xnewmat2.o:    newmatxx newmatrc.h controlw.h newmat2.cxx
X
Xnewmat3.o:    newmatxx newmatrc.h controlw.h newmat3.cxx
X
Xnewmat4.o:    newmatxx newmatrc.h controlw.h newmat4.cxx
X
Xnewmat5.o:    newmatxx newmatrc.h controlw.h newmat5.cxx
X
Xnewmat6.o:    newmatxx newmatrc.h controlw.h newmat6.cxx
X
Xnewmat7.o:    newmatxx newmatrc.h controlw.h newmat7.cxx
X
Xnewmat8.o:    newmatxx newmatap.h newmat8.cxx
X
Xnewmat9.o:    newmatxx newmatrc.h controlw.h newmatio.h newmat9.cxx
X
Xnewmatrm.o:   newmatxx newmatrm.h newmatrm.cxx
X
Xsort.o:       newmatxx newmatap.h sort.cxx
X
Xsubmat.o:     newmatxx newmatrc.h controlw.h submat.cxx
X
Xsvd.o:        newmatxx newmatrm.h precisio.h svd.cxx
X
Xtmt.o:        newmatxx newmatap.h tmt.cxx
X
Xtmt1.o:       newmatxx newmatap.h tmt1.cxx
X
Xtmt2.o:       newmatxx newmatap.h tmt2.cxx
X
Xtmt3.o:       newmatxx newmatap.h tmt3.cxx
X
Xtmt4.o:       newmatxx newmatap.h tmt4.cxx
X
Xtmt5.o:       newmatxx newmatap.h tmt5.cxx
X
Xtmt6.o:       newmatxx newmatap.h tmt6.cxx
X
Xtmt7.o:       newmatxx newmatap.h tmt7.cxx
X
Xtmt8.o:       newmatxx newmatap.h tmt8.cxx
X
Xtmt9.o:       newmatxx newmatap.h tmt9.cxx
X
Xtmta.o:       newmatxx newmatap.h tmta.cxx
X
Xtmtb.o:       newmatxx newmatap.h tmtb.cxx
X
Xtmtc.o:       newmatxx newmatap.h tmtc.cxx
X
Xtmtd.o:       newmatxx newmatap.h tmtd.cxx
X
Xtmte.o:       newmatxx newmatap.h tmte.cxx
X
Xtmtf.o:       newmatxx newmatap.h tmtf.cxx
X
Xtmtg.o:       newmatxx newmatap.h tmtg.cxx
X
Xtmth.o:       newmatxx newmatap.h tmth.cxx
X
Xtmti.o:       newmatxx newmatap.h tmti.cxx
X
X
END_OF_FILE
if test 2558 -ne `wc -c <'tmt.mak'`; then
    echo shar: \"'tmt.mak'\" unpacked with wrong size!
fi
# end of 'tmt.mak'
fi
echo shar: End of archive 3 \(of 7\).
cp /dev/null ark3isdone
MISSING=""
for I in 1 2 3 4 5 6 7 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 7 archives.
    rm -f ark[1-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0

exit 0 # Just in case...