newmat.h

矩阵计算库

#endif
   void operator=(const BaseMatrix&);
   void operator=(Real f) { GeneralMatrix::operator=(f); }
   void operator=(const ColumnVector& m) { operator=((const BaseMatrix&)m); }
   Real& operator()(int);                       // access element
   Real& element(int);                          // access element
#ifdef SETUP_C_SUBSCRIPTS
   Real& operator[](int m) { return store[m]; }
   const Real& operator[](int m) const { return store[m]; }
#endif
   MatrixType Type() const;
   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
   void ReDimension(int);                       // change dimensions
   void ReDimension(int,int);                   // in case access is matrix
#ifndef __ZTC__
   Real* nric() const
      { CheckStore(); return store-1; }         // for use by NRIC
#endif
   void CleanUp();                              // to clear store
   NEW_DELETE(ColumnVector)
};

class CroutMatrix : public GeneralMatrix        // for LU decomposition
{
   int* indx;
   Boolean d;
   Boolean sing;
   void ludcmp();
public:
   CroutMatrix(const BaseMatrix&);
   MatrixType Type() const;
   void lubksb(Real*, int=0);
   ~CroutMatrix();
   GeneralMatrix* MakeSolver() { return this; } // for solving
   LogAndSign LogDeterminant() const;
   void Solver(MatrixRowCol&, const MatrixRowCol&);
   void GetRow(MatrixRowCol&);
   void GetCol(MatrixRowCol&);
   void operator=(const BaseMatrix&);
   void operator=(const CroutMatrix& m) { operator=((const BaseMatrix&)m); }
   void CleanUp();                                // to clear store
   NEW_DELETE(CroutMatrix)
};

/******************************* band matrices ***************************/

class BandMatrix : public GeneralMatrix         // band matrix
{
   GeneralMatrix* Image() const;                // copy of matrix
protected:
   void CornerClear() const;                    // set unused elements to zero
public:
   int lower, upper;                            // band widths
   BandMatrix() { lower=0; upper=0; CornerClear(); }
   BandMatrix(int n,int lb,int ub) { ReDimension(n,lb,ub); CornerClear(); }
                                                // standard declaration
   BandMatrix(const BaseMatrix&);               // evaluate BaseMatrix
   void operator=(const BaseMatrix&);
   void operator=(Real f) { GeneralMatrix::operator=(f); }
   void operator=(const BandMatrix& m) { operator=((const BaseMatrix&)m); }
   MatrixType Type() const;
   Real& operator()(int, int);                  // access element
   Real& element(int, int);                     // access element
#ifdef SETUP_C_SUBSCRIPTS
   Real* operator[](int m) { return store+(upper+lower)*m+lower; }
   const Real* operator[](int m) const { return store+(upper+lower)*m+lower; }
#endif
   BandMatrix(const BandMatrix& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   Real operator()(int, int) const;             // access element
#endif
   LogAndSign LogDeterminant() const;
   GeneralMatrix* MakeSolver();
   Real Trace() const;
   Real SumSquare() const { CornerClear(); return GeneralMatrix::SumSquare(); }
   Real SumAbsoluteValue() const
      { CornerClear(); return GeneralMatrix::SumAbsoluteValue(); }
   Real Sum() const
      { CornerClear(); return GeneralMatrix::Sum(); }
   Real MaximumAbsoluteValue() const
      { CornerClear(); return GeneralMatrix::MaximumAbsoluteValue(); }
   void GetRow(MatrixRowCol&);
   void GetCol(MatrixRowCol&);
   void RestoreCol(MatrixRowCol&);
   void NextRow(MatrixRowCol&);
   void ReDimension(int, int, int);             // change dimensions
   MatrixBandWidth BandWidth() const;
   void SetParameters(const GeneralMatrix*);
   MatrixInput operator<<(Real);                // will give error
   void operator<<(const Real* r);              // will give error
      // the next is included because Zortech and Borland
      // can't find the copy in GeneralMatrix
   void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
   NEW_DELETE(BandMatrix)
};

class UpperBandMatrix : public BandMatrix       // upper band matrix
{
   GeneralMatrix* Image() const;                // copy of matrix
public:
   UpperBandMatrix() {}
   UpperBandMatrix(int n, int ubw)              // standard declaration
      : BandMatrix(n, 0, ubw) {}
   UpperBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
   void operator=(const BaseMatrix&);
   void operator=(Real f) { GeneralMatrix::operator=(f); }
   void operator=(const UpperBandMatrix& m)
      { operator=((const BaseMatrix&)m); }
   MatrixType Type() const;
   UpperBandMatrix(const UpperBandMatrix& gm) { GetMatrix(&gm); }
   GeneralMatrix* MakeSolver() { return this; }
   void Solver(MatrixRowCol&, const MatrixRowCol&);
   LogAndSign LogDeterminant() const;
   void ReDimension(int n,int ubw)              // change dimensions
      { BandMatrix::ReDimension(n,0,ubw); }
   Real& operator()(int, int);
   Real& element(int, int);
#ifdef SETUP_C_SUBSCRIPTS
   Real* operator[](int m) { return store+upper*m; }
   const Real* operator[](int m) const { return store+upper*m; }
#endif
   NEW_DELETE(UpperBandMatrix)
};

class LowerBandMatrix : public BandMatrix       // upper band matrix
{
   GeneralMatrix* Image() const;                // copy of matrix
public:
   LowerBandMatrix() {}
   LowerBandMatrix(int n, int lbw)              // standard declaration
      : BandMatrix(n, lbw, 0) {}
   LowerBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
   void operator=(const BaseMatrix&);
   void operator=(Real f) { GeneralMatrix::operator=(f); }
   void operator=(const LowerBandMatrix& m)
      { operator=((const BaseMatrix&)m); }
   MatrixType Type() const;
   LowerBandMatrix(const LowerBandMatrix& gm) { GetMatrix(&gm); }
   GeneralMatrix* MakeSolver() { return this; }
   void Solver(MatrixRowCol&, const MatrixRowCol&);
   LogAndSign LogDeterminant() const;
   void ReDimension(int n,int lbw)             // change dimensions
      { BandMatrix::ReDimension(n,lbw,0); }
   Real& operator()(int, int);
   Real& element(int, int);
#ifdef SETUP_C_SUBSCRIPTS
   Real* operator[](int m) { return store+lower*(m+1); }
   const Real* operator[](int m) const { return store+lower*(m+1); }
#endif
   NEW_DELETE(LowerBandMatrix)
};

class SymmetricBandMatrix : public GeneralMatrix
{
   GeneralMatrix* Image() const;                // copy of matrix
   void CornerClear() const;                    // set unused elements to zero
public:
   int lower;                                   // lower band width
   SymmetricBandMatrix() { lower=0; CornerClear(); }
   SymmetricBandMatrix(int n, int lb) { ReDimension(n,lb); CornerClear(); }
   SymmetricBandMatrix(const BaseMatrix&);
   void operator=(const BaseMatrix&);
   void operator=(Real f) { GeneralMatrix::operator=(f); }
   void operator=(const SymmetricBandMatrix& m)
      { operator=((const BaseMatrix&)m); }
   Real& operator()(int, int);                  // access element
   Real& element(int, int);                     // access element
#ifdef SETUP_C_SUBSCRIPTS
   Real* operator[](int m) { return store+lower*(m+1); }
   const Real* operator[](int m) const { return store+lower*(m+1); }
#endif
   MatrixType Type() const;
   SymmetricBandMatrix(const SymmetricBandMatrix& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   Real operator()(int, int) const;             // access element
#endif
   GeneralMatrix* MakeSolver();
   Real SumSquare() const;
   Real SumAbsoluteValue() const;
   Real Sum() const;
   Real MaximumAbsoluteValue() const
      { CornerClear(); return GeneralMatrix::MaximumAbsoluteValue(); }
   Real Trace() const;
   LogAndSign LogDeterminant() const;
   void GetRow(MatrixRowCol&);
   void GetCol(MatrixRowCol&);
   void RestoreCol(MatrixRowCol&);
   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
   void ReDimension(int,int);                       // change dimensions
   MatrixBandWidth BandWidth() const;
   void SetParameters(const GeneralMatrix*);
   NEW_DELETE(SymmetricBandMatrix)
};

class BandLUMatrix : public GeneralMatrix
// for LU decomposition of band matrix
{
   int* indx;
   Boolean d;
   Boolean sing;                                   // TRUE if singular
   Real* store2;
   int storage2;
   void ludcmp();
   int m1,m2;                                   // lower and upper
public:
   BandLUMatrix(const BaseMatrix&);
   MatrixType Type() const;
   void lubksb(Real*, int=0);
   ~BandLUMatrix();
   GeneralMatrix* MakeSolver() { return this; } // for solving
   LogAndSign LogDeterminant() const;
   void Solver(MatrixRowCol&, const MatrixRowCol&);
   void GetRow(MatrixRowCol&);
   void GetCol(MatrixRowCol&);
   void operator=(const BaseMatrix&);
   void operator=(const BandLUMatrix& m) { operator=((const BaseMatrix&)m); }
   NEW_DELETE(BandLUMatrix)
   void CleanUp();                                // to clear store
};

/**************************** GenericMatrix class ************************/

class GenericMatrix : public BaseMatrix
{
   GeneralMatrix* gm;
   int search(const BaseMatrix* bm) const;
   friend class BaseMatrix;
public:
   GenericMatrix() : gm(0) {}
   GenericMatrix(const BaseMatrix& bm)
      { gm = ((BaseMatrix&)bm).Evaluate(); gm = gm->Image(); }
   GenericMatrix(const GenericMatrix& bm)
      { gm = bm.gm->Image(); }
   void operator=(const GenericMatrix&);
   void operator=(const BaseMatrix&);
   ~GenericMatrix() { delete gm; }
   void CleanUp() { delete gm; gm = 0; }
   void Release() { gm->Release(); } 
   GeneralMatrix* Evaluate(MatrixType) { return gm; };
   MatrixBandWidth BandWidth() const;
   NEW_DELETE(GenericMatrix)
};

/***************************** temporary classes *************************/

class MultipliedMatrix : public BaseMatrix
{
protected:
   // if these union statements cause problems, simply remove them
   // and declare the items individually
   union { const BaseMatrix* bm1; GeneralMatrix* gm1; };
						  // pointers to summands
   union { const BaseMatrix* bm2; GeneralMatrix* gm2; };
   MultipliedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
      : bm1(bm1x),bm2(bm2x) {}
   int search(const BaseMatrix*) const;
   friend class BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
   MatrixBandWidth BandWidth() const;
   NEW_DELETE(MultipliedMatrix)
};

class AddedMatrix : public MultipliedMatrix
{
protected:
   AddedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
      : MultipliedMatrix(bm1x,bm2x) {}

   friend class BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
   MatrixBandWidth BandWidth() const;
#ifdef __GNUG__
   void SelectVersion(MatrixType, int&, int&) const;
#else
   void SelectVersion(MatrixType, Boolean&, Boolean&) const;
#endif
   NEW_DELETE(AddedMatrix)
};

class SPMatrix : public AddedMatrix
{
protected:
   SPMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
      : AddedMatrix(bm1x,bm2x) {}

   friend class BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
   MatrixBandWidth BandWidth() const;
#ifdef __GNUG__
   void SelectVersion(MatrixType, int&, int&) const;
#else
   void SelectVersion(MatrixType, Boolean&, Boolean&) const;
#endif
   NEW_DELETE(AddedMatrix)

#ifndef TEMPS_DESTROYED_QUICKLY
   friend SPMatrix SP(const BaseMatrix&, const BaseMatrix&);
#else
   friend SPMatrix& SP(const BaseMatrix&, const BaseMatrix&);
#endif
};

class ConcatenatedMatrix : public MultipliedMatrix
{
protected:
   ConcatenatedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
      : MultipliedMatrix(bm1x,bm2x) {}

   friend class BaseMatrix;
public:
   MatrixBandWidth BandWidth() const;
   GeneralMatrix* Evaluate(MatrixType);
   NEW_DELETE(ConcatenatedMatrix)
};

class StackedMatrix : public ConcatenatedMatrix
{
protected:
   StackedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
      : ConcatenatedMatrix(bm1x,bm2x) {}

   friend class BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
   NEW_DELETE(StackedMatrix)
};

class SolvedMatrix : public MultipliedMatrix
{
   SolvedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
      : MultipliedMatrix(bm1x,bm2x) {}
   friend class BaseMatrix;
   friend class InvertedMatrix;                        // for operator*
public:
   GeneralMatrix* Evaluate(MatrixType);
   MatrixBandWidth BandWidth() const;
   NEW_DELETE(SolvedMatrix)
};

class SubtractedMatrix : public AddedMatrix
{
   SubtractedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
      : AddedMatrix(bm1x,bm2x) {}
   friend class BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
   NEW_DELETE(SubtractedMatrix)
};

class ShiftedMatrix : public BaseMatrix
{
protected:
   Real f;
   union { const BaseMatrix* bm; GeneralMatrix* gm; };
   ShiftedMatrix(const BaseMatrix* bmx, Real fx) : bm(bmx),f(fx) {}
   int search(const BaseMatrix*) const;
   friend class BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
   NEW_DELETE(ShiftedMatrix)
};

class ScaledMatrix : public ShiftedMatrix
{
   ScaledMatrix(const BaseMatrix* bmx, Real fx) : ShiftedMatrix(bmx,fx) {}
   friend class BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
   MatrixBandWidth BandWidth() const;
   NEW_DELETE(ScaledMatrix)