newmat.h

矩阵计算库

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

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

#ifndef MATRIX_LIB
#define MATRIX_LIB 0

#include "include.h"

#ifdef NO_LONG_NAMES
#define UpperTriangularMatrix UTMatrix
#define LowerTriangularMatrix LTMatrix
#define SymmetricMatrix SMatrix
#define DiagonalMatrix DMatrix
#define BandMatrix BMatrix
#define UpperBandMatrix UBMatrix
#define LowerBandMatrix LBMatrix
#define SymmetricBandMatrix SBMatrix
#define BandLUMatrix BLUMatrix
#endif

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

#include "boolean.h"
#include "myexcept.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 to show whether to check for loss of data ************/

class MatrixConversionCheck : public Janitor
{
   static Boolean DoCheck;
public:
   MatrixConversionCheck() { DoCheck=TRUE; }          // turn check on
   ~MatrixConversionCheck() { DoCheck=FALSE; }        // turn check off
	void CleanUp() { DoCheck=FALSE; }
	static Boolean IsOn() { return DoCheck; }
	static void DataLoss();
	friend class SkipConversionCheck;
};

class SkipConversionCheck : public Janitor
// to turn check off in current block
{
	Boolean LastValue;
public:
	SkipConversionCheck() : LastValue(MatrixConversionCheck::DoCheck)
		{ MatrixConversionCheck::DoCheck = FALSE; }
	~SkipConversionCheck() { MatrixConversionCheck::DoCheck = LastValue; }
	void CleanUp() { MatrixConversionCheck::DoCheck = LastValue; }
};

/**************************** 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 MatrixInput;                              // defined later

class MatrixType
{
public:
   enum Attribute {  Valid     = 1,
                     Diagonal  = 2,             // order of these is important
                     Symmetric = 4,
                     Band      = 8,
                     Lower     = 16,
                     Upper     = 32,
                     LUDeco    = 64 };

   enum            { US = 0,
                     UT = Valid + Upper,
                     LT = Valid + Lower,
                     Rt = Valid,
                     Sm = Valid + Symmetric,
                     Dg = Valid + Diagonal + Band + Lower + Upper + Symmetric,
		     RV = Valid,     //   don't separate out
		     CV = Valid,     //   vectors
		     BM = Valid + Band,
		     UB = Valid + Band + Upper,
		     LB = Valid + Band + Lower,
		     SB = Valid + Band + Symmetric,
		     Ct = Valid + LUDeco,
		     BC = Valid + Band + LUDeco
                   };


   static nTypes() { return 9; }               // number of different types
					       // exclude Ct, US, BC
public:
   int attribute;
public:
   MatrixType () {}
   MatrixType (int i) : attribute(i) {}
   int operator+() const { return attribute; }
   MatrixType operator+(MatrixType mt) const
      { return MatrixType(attribute & mt.attribute); }
   MatrixType operator*(const MatrixType&) const;
   MatrixType SP(const MatrixType&) const;
   MatrixType operator|(const MatrixType& mt) const
      { return MatrixType(attribute & mt.attribute & Valid); }
   MatrixType operator&(const MatrixType& mt) const
      { return MatrixType(attribute & mt.attribute & Valid); }
   Boolean operator>=(MatrixType mt) const
      { return ( attribute & mt.attribute ) == attribute; }
   Boolean operator==(MatrixType t) const
      { return (attribute == t.attribute); }
   Boolean operator!=(MatrixType t) const
      { return (attribute != t.attribute); }
   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
      { return MatrixType(attribute & (Valid + Symmetric)); }
   MatrixType MultRHS() const;                 // type for rhs of multiply
   MatrixType sub() const                      // type of submatrix
      { return MatrixType(attribute & Valid); }
   MatrixType ssub() const                     // type of sym submatrix
      { return MatrixType(attribute); }        // not for selection matrix
   GeneralMatrix* New() const;                 // new matrix of given type
   GeneralMatrix* New(int,int,BaseMatrix*) const;
                                               // new matrix of given type
   char* Value() const;                        // to print type
   friend Boolean Rectangular(MatrixType a, MatrixType b, MatrixType c);
   friend Boolean Compare(const MatrixType&, MatrixType&);
                                               // compare and check conv.
   Boolean IsBand() const { return (attribute & Band) != 0; }
   FREE_CHECK(MatrixType)
};

void TestTypeAdd();                            // test +
void TestTypeMult();                           // test *
void TestTypeSP();                             // test SP
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 minimum(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 SPMatrix;
class ConcatenatedMatrix;
class StackedMatrix;
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;
class GenericMatrix;



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

class BaseMatrix : public Janitor               // base of all matrix classes
{
protected:
   virtual int search(const BaseMatrix*) const = 0;
						// count number of times matrix
						// is referred to
public:
   virtual GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp) = 0;
						// evaluate temporary
   // for old version of G++
   //   virtual GeneralMatrix* Evaluate(MatrixType mt) = 0;
   //   GeneralMatrix* Evaluate() { return Evaluate(MatrixTypeUnSp); }
#ifndef TEMPS_DESTROYED_QUICKLY
   AddedMatrix operator+(const BaseMatrix&) const;    // results of operations
   MultipliedMatrix operator*(const BaseMatrix&) const;
   SubtractedMatrix operator-(const BaseMatrix&) const;
   ConcatenatedMatrix operator|(const BaseMatrix&) const;
   StackedMatrix 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;
   ConcatenatedMatrix& operator|(const BaseMatrix&) const;
   StackedMatrix& 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 Sum() 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 class GeneralMatrix;
   friend class Matrix;
   friend class nricMatrix;
   friend class RowVector;
   friend class ColumnVector;
   friend class SymmetricMatrix;
   friend class UpperTriangularMatrix;
   friend class LowerTriangularMatrix;
   friend class DiagonalMatrix;
   friend class CroutMatrix;
   friend class BandMatrix;
   friend class LowerBandMatrix;
   friend class UpperBandMatrix;
   friend class SymmetricBandMatrix;
   friend class AddedMatrix;
   friend class MultipliedMatrix;
   friend class SubtractedMatrix;
   friend class SPMatrix;
   friend class ConcatenatedMatrix;
   friend class StackedMatrix;
   friend class SolvedMatrix;
   friend class ShiftedMatrix;
   friend class ScaledMatrix;
   friend class TransposedMatrix;
   friend class NegatedMatrix;
   friend class InvertedMatrix;
   friend class RowedMatrix;
   friend class ColedMatrix;
   friend class DiagedMatrix;
   friend class MatedMatrix;