newmat.h

科学和工程计算中使用统计功能开发工具包

//$$ newmat.h           definition file for new version of matrix package

// Copyright (C) 1991,2: R B Davies

#ifndef MATRIX_LIB
#define MATRIX_LIB 0

#ifdef NO_LONG_NAMES
#define UpperTriangularMatrix UTMatrix
#define LowerTriangularMatrix LTMatrix
#define SymmetricMatrix SMatrix
#define DiagonalMatrix DMatrix
#endif

#ifndef TEMPS_DESTROYED_QUICKLY
#define ReturnMatrix ReturnMatrixX
#else
#define ReturnMatrix ReturnMatrixX&
#endif

#include "boolean.h"
#include "except.h"

/**************************** general utilities ****************************/

class GeneralMatrix;
void MatrixErrorNoSpace(void*);                 // no space handler

class LogAndSign
// Return from LogDeterminant function
//    - value of the log plus the sign (+, - or 0)
{
   Real log_value;
   int sign;
public:
   LogAndSign() { log_value=0.0; sign=1; }
   LogAndSign(Real);
   void operator*=(Real);
   void ChangeSign() { sign = -sign; }
   Real LogValue() const { return log_value; }
   int Sign() const { return sign; }
   Real Value() const;
   FREE_CHECK(LogAndSign)
};

// the following class is for counting the number of times a piece of code
// is executed. It is used for locating any code not executed by test
// routines. Use turbo GREP locate all places this code is called and
// check which ones are not accessed.
// Somewhat implementation dependent as it relies on "cout" still being
// present when ExeCounter objects are destructed.

class ExeCounter
{
   int line;                                    // code line number
   int fileid;                                  // file identifier
   long nexe;                                   // number of executions
   static int nreports;                         // number of reports
public:
   ExeCounter(int,int);
   void operator++() { nexe++; }
   ~ExeCounter();                               // prints out reports
};


/**************************** class MatrixType *****************************/

// Is used for finding the type of a matrix resulting from the binary operations
// +, -, * and identifying what conversions are permissible.
// This class must be updated when new matrix types are added.

class GeneralMatrix;                            // defined later
class BaseMatrix;                               // defined later

class MatrixType
{
public:
   enum Attribute {  Valid     = 1,
                     Symmetric = 2,
                     Band      = 4,
                     Upper     = 8,
                     Lower     = 16,
                     LUDeco    = 32,
                     OneRow    = 64,
                     OneCol    = 128 };

   enum            { US = 0,
                     UT = Valid + Upper,
                     LT = Valid + Lower,
                     Rt = Valid,
                     Sm = Valid + Symmetric,
                     Dg = Valid + Band + Lower + Upper + Symmetric,
		     RV = Valid + OneRow,
		     CV = Valid + OneCol,
		     BM = Valid + Band,
		     UB = Valid + Band + Upper,
		     LB = Valid + Band + Lower,
		     SB = Valid + Band + Symmetric,
		     Ct = Valid + LUDeco,
		     BC = Valid + Band + LUDeco,
                   };

   enum { USX,UTX,LTX,RtX,SmX,DgX,RVX,CVX,BMX,UBX,LBX,SBX,CtX,BCX };

   static nTypes() { return 11; }              // number of different types
					       // exclude Ct, US, BC
public:
   int attribute;
public:
   MatrixType operator+(const MatrixType& t) const;
   MatrixType operator*(const MatrixType&) const;
   Boolean operator>=(const MatrixType& t) const;
   Boolean operator==(const MatrixType& t) const
      { return (attribute == t.attribute); }
   Boolean operator!=(const MatrixType& t) const { return !operator==(t); }
   Boolean operator!() const { return (attribute & Valid) == 0; }
   MatrixType i() const;                       // type of inverse
   MatrixType t() const;                       // type of transpose
   MatrixType AddEqualEl() const;              // Add constant to matrix
   MatrixType MultRHS() const;                 // type for rhs of multiply
   MatrixType sub() const;                     // type of submatrix
   MatrixType ssub() const;                    // type of sym submatrix
//   MatrixType (Type tx) : attribute((int)tx) {}
//					       // (& doesn't work with AT&T)
   MatrixType () {}
   MatrixType (int i) : attribute(i) {}
   GeneralMatrix* New() const;                 // new matrix of given type
   GeneralMatrix* New(int,int,BaseMatrix*) const;
                                               // new matrix of given type
   operator char*() const;                     // for printing type
   int BestFit() const;
   FREE_CHECK(MatrixType)
};

void TestTypeAdd();                            // test +
void TestTypeMult();                           // test *
void TestTypeOrder();                          // test >=


/************************* class MatrixBandWidth ***********************/

class MatrixBandWidth
{
public:
   int lower;
   int upper;
   MatrixBandWidth(const int l, const int u) : lower(l), upper (u) {}
   MatrixBandWidth(const int i) : lower(i), upper(i) {}
   MatrixBandWidth operator+(const MatrixBandWidth&) const;
   MatrixBandWidth operator*(const MatrixBandWidth&) const;
   MatrixBandWidth t() const { return MatrixBandWidth(upper,lower); }
   Boolean operator==(const MatrixBandWidth& bw) const
      { return (lower == bw.lower) && (upper == bw.upper); }
   FREE_CHECK(MatrixBandWidth)
};


/*********************** Array length specifier ************************/

// This class is introduced to avoid constructors such as
//   ColumnVector(int)
// being used for conversions

class ArrayLengthSpecifier
{
   int value;
public:
   int Value() const { return value; }
   ArrayLengthSpecifier(int l) : value(l) {}
   FREE_CHECK(ArrayLengthSpecifier)
};


/*************************** Matrix routines ***************************/


class MatrixRowCol;                             // defined later
class MatrixRow;
class MatrixCol;

class GeneralMatrix;                            // defined later
class AddedMatrix;
class MultipliedMatrix;
class SubtractedMatrix;
class SolvedMatrix;
class ShiftedMatrix;
class ScaledMatrix;
class TransposedMatrix;
class NegatedMatrix;
class InvertedMatrix;
class RowedMatrix;
class ColedMatrix;
class DiagedMatrix;
class MatedMatrix;
class GetSubMatrix;
class ConstMatrix;
class ReturnMatrixX;
class Matrix;
class nricMatrix;
class RowVector;
class ColumnVector;
class SymmetricMatrix;
class UpperTriangularMatrix;
class LowerTriangularMatrix;
class DiagonalMatrix;
class CroutMatrix;
class BandMatrix;
class LowerBandMatrix;
class UpperBandMatrix;
class SymmetricBandMatrix;
class LinearEquationSolver;

static MatrixType MatrixTypeUnSp(MatrixType::US);
						// AT&T needs this

class BaseMatrix : public Janitor               // base of all matrix classes
{
protected:
//   BaseMatrix() {}
   virtual ~BaseMatrix() {}
   virtual int search(const BaseMatrix*) const = 0;
						// count number of times matrix
						// is referred to
   virtual MatrixType Type() const = 0;         // type of a matrix
//   virtual int NrowsV()  const = 0;
//   virtual int NcolsV()  const = 0;
public:
#ifndef GXX
   virtual GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp) = 0;
						// evaluate temporary
#else
   virtual GeneralMatrix* Evaluate(MatrixType mt) = 0;
   GeneralMatrix* Evaluate() { return Evaluate(MatrixTypeUnSp); }
#endif
//   void MatrixParameters(int& nr, int& nc, MatrixType& mt)
//      { nr = NrowsV(); nc = NcolsV(); mt = Type(); }
#ifndef TEMPS_DESTROYED_QUICKLY
   AddedMatrix operator+(const BaseMatrix&) const;    // results of operations
   MultipliedMatrix operator*(const BaseMatrix&) const;
   SubtractedMatrix operator-(const BaseMatrix&) const;
   ShiftedMatrix operator+(Real) const;
   ScaledMatrix operator*(Real) const;
   ScaledMatrix operator/(Real) const;
   ShiftedMatrix operator-(Real) const;
   TransposedMatrix t() const;
//   TransposedMatrix t;
   NegatedMatrix operator-() const;                   // change sign of elements
   InvertedMatrix i() const;
//   InvertedMatrix i;
   RowedMatrix AsRow() const;
   ColedMatrix AsColumn() const;
   DiagedMatrix AsDiagonal() const;
   MatedMatrix AsMatrix(int,int) const;
   GetSubMatrix SubMatrix(int,int,int,int) const;
   GetSubMatrix SymSubMatrix(int,int) const;
   GetSubMatrix Row(int) const;
   GetSubMatrix Rows(int,int) const;
   GetSubMatrix Column(int) const;
   GetSubMatrix Columns(int,int) const;
#else
   AddedMatrix& operator+(const BaseMatrix&) const;    // results of operations
   MultipliedMatrix& operator*(const BaseMatrix&) const;
   SubtractedMatrix& operator-(const BaseMatrix&) const;
   ShiftedMatrix& operator+(Real) const;
   ScaledMatrix& operator*(Real) const;
   ScaledMatrix& operator/(Real) const;
   ShiftedMatrix& operator-(Real) const;
   TransposedMatrix& t() const;
//   TransposedMatrix& t;
   NegatedMatrix& operator-() const;                   // change sign of elements
   InvertedMatrix& i() const;
//   InvertedMatrix& i;
   RowedMatrix& AsRow() const;
   ColedMatrix& AsColumn() const;
   DiagedMatrix& AsDiagonal() const;
   MatedMatrix& AsMatrix(int,int) const;
   GetSubMatrix& SubMatrix(int,int,int,int) const;
   GetSubMatrix& SymSubMatrix(int,int) const;
   GetSubMatrix& Row(int) const;
   GetSubMatrix& Rows(int,int) const;
   GetSubMatrix& Column(int) const;
   GetSubMatrix& Columns(int,int) const;
#endif
   Real AsScalar() const;                      // conversion of 1 x 1 matrix
   virtual LogAndSign LogDeterminant() const;
   virtual Real SumSquare() const;
   virtual Real SumAbsoluteValue() const;
   virtual Real MaximumAbsoluteValue() const;
   virtual Real Trace() const;
   Real Norm1() const;
   Real NormInfinity() const;
   virtual MatrixBandWidth BandWidth() const;  // bandwidths of band matrix
   virtual void CleanUp() {}                     // to clear store
//protected:
//   BaseMatrix() : t(this), i(this) {}

   friend GeneralMatrix;
   friend Matrix;
   friend nricMatrix;
   friend RowVector;
   friend ColumnVector;
   friend SymmetricMatrix;
   friend UpperTriangularMatrix;
   friend LowerTriangularMatrix;
   friend DiagonalMatrix;
   friend CroutMatrix;
   friend BandMatrix;
   friend LowerBandMatrix;
   friend UpperBandMatrix;
   friend SymmetricBandMatrix;
   friend AddedMatrix;
   friend MultipliedMatrix;
   friend SubtractedMatrix;
   friend SolvedMatrix;
   friend ShiftedMatrix;
   friend ScaledMatrix;
   friend TransposedMatrix;
   friend NegatedMatrix;
   friend InvertedMatrix;
   friend RowedMatrix;
   friend ColedMatrix;
   friend DiagedMatrix;
   friend MatedMatrix;
   friend GetSubMatrix;
   friend ConstMatrix;
   friend ReturnMatrixX;
   friend LinearEquationSolver;
   NEW_DELETE(BaseMatrix)
};


/******************************* working classes **************************/

class GeneralMatrix : public BaseMatrix         // declarable matrix types
{
protected:
   int tag;                                     // shows whether can reuse
   int nrows, ncols;                            // dimensions
   int storage;                                 // total store required
   Real* store;                                 // point to store (0=not set)
   GeneralMatrix();                             // initialise with no store
   GeneralMatrix(ArrayLengthSpecifier);         // constructor getting store
   void Add(GeneralMatrix*, Real);              // sum of GM and Real
   void Add(Real);                              // add Real to this
   void Multiply(GeneralMatrix*, Real);         // product of GM and Real
   void Multiply(Real);                         // multiply this by Real
   void Negate(GeneralMatrix*);                 // change sign
   void Negate();                               // change sign
   void operator=(Real);                        // set matrix to constant
   Real* GetStore();                            // get store or copy
   GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
                                                // temporarily access store
   void GetMatrix(const GeneralMatrix*);        // used by = and initialise
   void Eq(const BaseMatrix&, MatrixType);      // used by =
   int search(const BaseMatrix*) const;
   virtual GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
   void CheckConversion(const BaseMatrix&);     // check conversion OK
   void ReDimension(int, int, int);             // change dimensions
//   int NrowsV() const;                          // get dimensions
//   int NcolsV() const;                          // virtual version
public:
   GeneralMatrix* Evaluate(MatrixType);
   MatrixType Type() const = 0;                 // type of a matrix
   int Nrows() const { return nrows; }          // get dimensions
   int Ncols() const { return ncols; }
   int Storage() const { return storage; }
   Real* Store() const { return store; }
   virtual ~GeneralMatrix();                    // delete store if set
   void tDelete();                              // delete if tag permits
   Boolean reuse();                                // TRUE if tag allows reuse
   void Protect() { tag=-1; }                   // can't delete or reuse
   int Tag() const { return tag; }
   Boolean IsZero() const;                         // test matrix has all zeros
   void Release() { tag=1; }                    // del store after next use
   void Release(int t) { tag=t; }               // del store after t accesses
   void ReleaseAndDelete() { tag=0; }           // delete matrix after use
   void operator<<(const Real*);                // assignment from an array
   void operator<<(const BaseMatrix& X) { Eq(X,this->Type()); }
                                                // = without checking type
   void Inject(const GeneralMatrix&);           // copy stored els only
   virtual GeneralMatrix* MakeSolver();         // for solving
   virtual void Solver(MatrixRowCol&, const MatrixRowCol&) {}
   virtual void GetRow(MatrixRowCol&) = 0;      // Get matrix row
   virtual void RestoreRow(MatrixRowCol&) {}    // Restore matrix row
   virtual void NextRow(MatrixRowCol&);         // Go to next row
   virtual void GetCol(MatrixRowCol&) = 0;      // Get matrix col
   virtual void RestoreCol(MatrixRowCol&) {}    // Restore matrix col
   virtual void NextCol(MatrixRowCol&);         // Go to next col
   Real SumSquare() const;
   Real SumAbsoluteValue() const;
   Real MaximumAbsoluteValue() const;
   LogAndSign LogDeterminant() const;
#ifndef TEMPS_DESTROYED_QUICKLY
   ConstMatrix c() const;                       // to access constant matrices
#else
   ConstMatrix& c() const;                      // to access constant matrices
#endif
   void CheckStore() const;                     // check store is non-zero
   virtual void SetParameters(const GeneralMatrix*) {}
                                                // set parameters in GetMatrix
#ifndef TEMPS_DESTROYED_QUICKLY
   operator ReturnMatrixX() const;              // for building a ReturnMatrix
#else
   operator ReturnMatrixX&() const;             // for building a ReturnMatrix
#endif
   void CleanUp();                                // to clear store

   friend Matrix;
   friend nricMatrix;
   friend SymmetricMatrix;
   friend UpperTriangularMatrix;
   friend LowerTriangularMatrix;
   friend DiagonalMatrix;
   friend CroutMatrix;
   friend RowVector;
   friend ColumnVector;
   friend BandMatrix;
   friend LowerBandMatrix;
   friend UpperBandMatrix;
   friend SymmetricBandMatrix;
   friend BaseMatrix;
   friend AddedMatrix;
   friend MultipliedMatrix;