newmat.h

C++矩阵算法库

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

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

#ifndef NEWMAT_LIB
#define NEWMAT_LIB 0

#include "include.h"

#include "boolean.h"
#include "myexcept.h"


#ifdef use_namespace
namespace NEWMAT { using namespace RBD_COMMON; }
namespace RBD_LIBRARIES { using namespace NEWMAT; }
namespace NEWMAT {
#endif

//#define DO_REPORT                     // to activate REPORT

#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

// ************************** 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 PowEq(int k);  // raise to power of k
   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.

#ifdef DO_REPORT

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
};

#endif


// ************************** 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,
                     Ones      = 128 };

   enum            { US = 0,
                     UT = Valid + Upper,
                     LT = Valid + Lower,
                     Rt = Valid,
                     Sm = Valid + Symmetric,
                     Dg = Valid + Diagonal + Band + Lower + Upper + Symmetric,
                     Id = Valid + Diagonal + Band + Lower + Upper + Symmetric
                        + Ones,
                     RV = Valid,     //   do not 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 int nTypes() { return 10; }           // number of different types
					       // exclude Ct, US, BC
public:
   int attribute;
   bool DataLossOK;                            // true if data loss is OK when
                                               // this represents a destination
public:
   MatrixType () : DataLossOK(false) {}
   MatrixType (int i) : attribute(i), DataLossOK(false) {}
   MatrixType (int i, bool dlok) : attribute(i), DataLossOK(dlok) {}
   MatrixType (const MatrixType& mt)
      : attribute(mt.attribute), DataLossOK(mt.DataLossOK) {}
   void operator=(const MatrixType& mt)
      { attribute = mt.attribute; DataLossOK = mt.DataLossOK; }
   void SetDataLossOK() { DataLossOK = true; }
   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 KP(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); }
   bool operator>=(MatrixType mt) const
      { return ( attribute & mt.attribute ) == attribute; }
   bool operator<(MatrixType mt) const         // for MS Visual C++ 4
      { return ( attribute & mt.attribute ) != attribute; }
   bool operator==(MatrixType t) const
      { return (attribute == t.attribute); }
   bool operator!=(MatrixType t) const
      { return (attribute != t.attribute); }
   bool 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
   const char* Value() const;                  // to print type
   friend bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
   friend bool Compare(const MatrixType&, MatrixType&);
                                               // compare and check conv.
   bool IsBand() const { return (attribute & Band) != 0; }
   bool IsDiagonal() const { return (attribute & Diagonal) != 0; }
   bool IsSymmetric() const { return (attribute & Symmetric) != 0; }
   bool CannotConvert() const { return (attribute & LUDeco) != 0; }
                                               // used by operator== 
   FREE_CHECK(MatrixType)
};


// *********************** 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); }
   bool operator==(const MatrixBandWidth& bw) const
      { return (lower == bw.lower) && (upper == bw.upper); }
   bool operator!=(const MatrixBandWidth& bw) const { return !operator==(bw); }
   int Upper() const { return upper; }
   int Lower() const { return lower; }
   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) {}
};

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


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

class GeneralMatrix;                            // defined later
class AddedMatrix;
class MultipliedMatrix;
class SubtractedMatrix;
class SPMatrix;
class KPMatrix;
class ConcatenatedMatrix;
class StackedMatrix;
class SolvedMatrix;
class ShiftedMatrix;
class NegShiftedMatrix;
class ScaledMatrix;
class TransposedMatrix;
class ReversedMatrix;
class NegatedMatrix;
class InvertedMatrix;
class RowedMatrix;
class ColedMatrix;
class DiagedMatrix;
class MatedMatrix;
class GetSubMatrix;
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;


#define MatrixTypeUnSp 0
//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
   ReversedMatrix Reverse() const;
   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
   ReversedMatrix& Reverse() const;
   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;
   Real Determinant() const;
   virtual Real SumSquare() const;