newmat.h

矩阵计算库

   friend class GetSubMatrix;
   friend class ConstMatrix;
   friend class ReturnMatrixX;
   friend class LinearEquationSolver;
	friend class GenericMatrix;
   NEW_DELETE(BaseMatrix)
};


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

class GeneralMatrix : public BaseMatrix         // declarable matrix types
{
   virtual GeneralMatrix* Image() const;        // copy of matrix
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
public:
   GeneralMatrix* Evaluate(MatrixType);
   virtual 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 Sum() 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
   operator ReturnMatrix() const;               // for building a ReturnMatrix
   ReturnMatrix ForReturn() const;
   MatrixInput operator<<(Real);                // for loading a list
   void CleanUp();                              // to clear store

   friend class Matrix;
   friend class nricMatrix;
   friend class SymmetricMatrix;
   friend class UpperTriangularMatrix;
   friend class LowerTriangularMatrix;
   friend class DiagonalMatrix;
   friend class CroutMatrix;
   friend class RowVector;
   friend class ColumnVector;
   friend class BandMatrix;
   friend class LowerBandMatrix;
   friend class UpperBandMatrix;
   friend class SymmetricBandMatrix;
   friend class BaseMatrix;
   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;
   friend class GetSubMatrix;
   friend class ConstMatrix;
   friend class ReturnMatrixX;
   friend class LinearEquationSolver;
   friend class GenericMatrix;
   NEW_DELETE(GeneralMatrix)
};

class Matrix : public GeneralMatrix             // usual rectangular matrix
{
   GeneralMatrix* Image() const;                // copy of matrix
public:
   Matrix() {}
   Matrix(int, int);                            // standard declaration
   Matrix(const BaseMatrix&);                   // evaluate BaseMatrix
   void operator=(const BaseMatrix&);
   void operator=(Real f) { GeneralMatrix::operator=(f); }
   void operator=(const Matrix& 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+m*ncols; }
   const Real* operator[](int m) const { return store+m*ncols; }
#endif
   Matrix(const Matrix& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   Real operator()(int, int) const;             // access element
#endif
   GeneralMatrix* MakeSolver();
   Real Trace() const;
   void GetRow(MatrixRowCol&);
   void GetCol(MatrixRowCol&);
   void RestoreCol(MatrixRowCol&);
   void NextRow(MatrixRowCol&);
   void NextCol(MatrixRowCol&);
   virtual void ReDimension(int,int);           // change dimensions
      // virtual so we will catch it being used in a vector called as a matrix
   NEW_DELETE(Matrix)
};

class nricMatrix : public Matrix                // for use with Numerical
                                                // Recipes in C
{
   GeneralMatrix* Image() const;                // copy of matrix
   Real** row_pointer;                          // points to rows
   void MakeRowPointer();                       // build rowpointer
   void DeleteRowPointer();
public:
   nricMatrix() {}
   nricMatrix(int m, int n)                     // standard declaration
      :  Matrix(m,n) { MakeRowPointer(); }
   nricMatrix(const BaseMatrix& bm)             // evaluate BaseMatrix
      :  Matrix(bm) { MakeRowPointer(); }
   void operator=(const BaseMatrix& bm)
      { DeleteRowPointer(); Matrix::operator=(bm); MakeRowPointer(); }
   void operator=(Real f) { GeneralMatrix::operator=(f); }
   void operator=(const nricMatrix& m) { operator=((const BaseMatrix&)m); }
   void operator<<(const BaseMatrix& X)
      { DeleteRowPointer(); Eq(X,this->Type()); MakeRowPointer(); }
   nricMatrix(const nricMatrix& gm) { GetMatrix(&gm); MakeRowPointer(); }
   void ReDimension(int m, int n)               // change dimensions
      { DeleteRowPointer(); Matrix::ReDimension(m,n); MakeRowPointer(); }
   ~nricMatrix() { DeleteRowPointer(); }
#ifndef __ZTC__
   Real** nric() const { CheckStore(); return row_pointer-1; }
#endif
   void CleanUp();                                // to clear store
   NEW_DELETE(nricMatrix)
};

class SymmetricMatrix : public GeneralMatrix
{
   GeneralMatrix* Image() const;                // copy of matrix
public:
   SymmetricMatrix() {}
   SymmetricMatrix(ArrayLengthSpecifier);
   SymmetricMatrix(const BaseMatrix&);
   void operator=(const BaseMatrix&);
   void operator=(Real f) { GeneralMatrix::operator=(f); }
   void operator=(const SymmetricMatrix& 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+(m*(m+1))/2; }
   const Real* operator[](int m) const { return store+(m*(m+1))/2; }
#endif
   MatrixType Type() const;
   SymmetricMatrix(const SymmetricMatrix& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   Real operator()(int, int) const;             // access element
#endif
   Real SumSquare() const;
   Real SumAbsoluteValue() const;
   Real Sum() const;
   Real Trace() const;
   void GetRow(MatrixRowCol&);
   void GetCol(MatrixRowCol&);
   void RestoreCol(MatrixRowCol&);
	GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
   void ReDimension(int);                       // change dimensions
   NEW_DELETE(SymmetricMatrix)
};

class UpperTriangularMatrix : public GeneralMatrix
{
   GeneralMatrix* Image() const;                // copy of matrix
public:
   UpperTriangularMatrix() {}
   UpperTriangularMatrix(ArrayLengthSpecifier);
   void operator=(const BaseMatrix&);
   void operator=(const UpperTriangularMatrix& m)
      { operator=((const BaseMatrix&)m); }
   UpperTriangularMatrix(const BaseMatrix&);
   UpperTriangularMatrix(const UpperTriangularMatrix& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   Real operator()(int, int) const;             // access element
#endif
   void operator=(Real f) { GeneralMatrix::operator=(f); }
   Real& operator()(int, int);                  // access element
   Real& element(int, int);                     // access element
#ifdef SETUP_C_SUBSCRIPTS
   Real* operator[](int m) { return store+m*ncols-(m*(m+1))/2; }
   const Real* operator[](int m) const { return store+m*ncols-(m*(m+1))/2; }
#endif
   MatrixType Type() const;
   GeneralMatrix* MakeSolver() { return this; } // for solving
   void Solver(MatrixRowCol&, const MatrixRowCol&);
   LogAndSign LogDeterminant() const;
   Real Trace() const;
   void GetRow(MatrixRowCol&);
   void GetCol(MatrixRowCol&);
   void RestoreCol(MatrixRowCol&);
   void NextRow(MatrixRowCol&);
   void ReDimension(int);                       // change dimensions
   NEW_DELETE(UpperTriangularMatrix)
};

class LowerTriangularMatrix : public GeneralMatrix
{
   GeneralMatrix* Image() const;                // copy of matrix
public:
   LowerTriangularMatrix() {}
   LowerTriangularMatrix(ArrayLengthSpecifier);
   LowerTriangularMatrix(const LowerTriangularMatrix& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   Real operator()(int, int) const;             // access element
#endif
   LowerTriangularMatrix(const BaseMatrix& M);
   void operator=(const BaseMatrix&);
   void operator=(Real f) { GeneralMatrix::operator=(f); }
   void operator=(const LowerTriangularMatrix& 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+(m*(m+1))/2; }
   const Real* operator[](int m) const { return store+(m*(m+1))/2; }
#endif
   MatrixType Type() const;
   GeneralMatrix* MakeSolver() { return this; } // for solving
   void Solver(MatrixRowCol&, const MatrixRowCol&);
   LogAndSign LogDeterminant() const;
   Real Trace() const;
   void GetRow(MatrixRowCol&);
   void GetCol(MatrixRowCol&);
   void RestoreCol(MatrixRowCol&);
   void NextRow(MatrixRowCol&);
   void ReDimension(int);                       // change dimensions
   NEW_DELETE(LowerTriangularMatrix)
};

class DiagonalMatrix : public GeneralMatrix
{
   GeneralMatrix* Image() const;                // copy of matrix
public:
   DiagonalMatrix() {}
   DiagonalMatrix(ArrayLengthSpecifier);
   DiagonalMatrix(const BaseMatrix&);
   DiagonalMatrix(const DiagonalMatrix& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   Real operator()(int, int) const;             // access element
   Real operator()(int) const;
#endif
   void operator=(const BaseMatrix&);
   void operator=(Real f) { GeneralMatrix::operator=(f); }
   void operator=(const DiagonalMatrix& m) { operator=((const BaseMatrix&)m); }
   Real& operator()(int, int);                  // access element
   Real& operator()(int);                       // access element
   Real& element(int, 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;

   LogAndSign LogDeterminant() const;
   Real Trace() const;
   void GetRow(MatrixRowCol&);
   void GetCol(MatrixRowCol&);
   void NextRow(MatrixRowCol&);
   void NextCol(MatrixRowCol&);
   GeneralMatrix* MakeSolver() { return this; } // for solving
   void Solver(MatrixRowCol&, const MatrixRowCol&);
   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
   void ReDimension(int);                       // change dimensions
#ifndef __ZTC__
   Real* nric() const
      { CheckStore(); return store-1; }         // for use by NRIC
#endif
   MatrixBandWidth BandWidth() const;
   NEW_DELETE(DiagonalMatrix)
};

class RowVector : public Matrix
{
   GeneralMatrix* Image() const;                // copy of matrix
public:
   RowVector() {}
   RowVector(ArrayLengthSpecifier n) : Matrix(1,n.Value()) {}
   RowVector(const BaseMatrix&);
   RowVector(const RowVector& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   Real operator()(int) const;                  // access element
#endif
   void operator=(const BaseMatrix&);
   void operator=(Real f) { GeneralMatrix::operator=(f); }
   void operator=(const RowVector& 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;
   void GetCol(MatrixRowCol&);
   void NextCol(MatrixRowCol&);
   void RestoreCol(MatrixRowCol&);
   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(RowVector)
};

class ColumnVector : public Matrix
{
   GeneralMatrix* Image() const;                // copy of matrix
public:
   ColumnVector() {}
   ColumnVector(ArrayLengthSpecifier n) : Matrix(n.Value(),1) {}
   ColumnVector(const BaseMatrix&);
   ColumnVector(const ColumnVector& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   Real operator()(int) const;                  // access element