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

// Copyright (C) 1991: R B Davies and DSIR

#ifndef MATRIX_LIB
#define MATRIX_LIB 0

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

#include "boolean.hxx"


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

void MatrixError(char*);                        // error handler
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; }
   void operator*=(real);
   void ChangeSign() { sign = -sign; }
   real LogValue() { return log_value; }
   int Sign() { return sign; }
   real Value();
};

// 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 MatrixType
{
public:
   enum Type { UnSp,UT,LT,Rect,Sym,Diag,RowV,ColV,EqEl,Crout };
   static nTypes() { return 8; }               // number of different types
					       // exclude Crout, UnSp
private:
   Type type;
public:
   MatrixType operator+(const MatrixType&) const;
   MatrixType operator-(const MatrixType&) const;
   MatrixType operator*(const MatrixType&) const;
   BOOL operator>=(const MatrixType&) const;
   BOOL operator==(const MatrixType& t) const; // { return (type == t.type); }
   BOOL operator!=(const MatrixType& t) const; // { return !(*this == t); }
   BOOL operator!() const { return type == UnSp; }
   MatrixType operator-() const;               // type of negative
   MatrixType i() const;                       // type of inverse
   MatrixType t() const;                       // type of transpose
   MatrixType sub() const;                     // type of submatrix
   MatrixType ssub() const;                    // type of sym submatrix
   MatrixType (Type tx) : type(tx) {}          // (& doesn't work with AT&T)
   MatrixType () {}
   GeneralMatrix* New() const;                 // new matrix of given type
   GeneralMatrix* New(int,int) const;          // new matrix of given type
   operator int() const { return (int)type; }
   operator char*() const;                     // for printing type
};

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

/*************************** 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 ReturnMatrix;
class Matrix;
class nricMatrix;
class RowVector;
class ColumnVector;
class SymmetricMatrix;
class UpperTriangularMatrix;
class LowerTriangularMatrix;
class DiagonalMatrix;
class CroutMatrix;

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

class BaseMatrix                                // base of all matrix classes
{
protected:
//   BaseMatrix() {}
   virtual int search(const GeneralMatrix*) 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(); }
   AddedMatrix operator+(BaseMatrix&);          // results of operations
   MultipliedMatrix operator*(BaseMatrix&);
   SubtractedMatrix operator-(BaseMatrix&);
   ShiftedMatrix operator+(real);
   ScaledMatrix operator*(real);
   ScaledMatrix operator/(real);
   ShiftedMatrix operator-(real);
   TransposedMatrix t();
   NegatedMatrix operator-();                   // change sign of elements
   InvertedMatrix i();
   RowedMatrix CopyToRow();
   ColedMatrix CopyToColumn();
   DiagedMatrix CopyToDiagonal();
   MatedMatrix CopyToMatrix(int,int);
   GetSubMatrix SubMatrix(int,int,int,int);
   GetSubMatrix SymSubMatrix(int,int);
   GetSubMatrix Row(int);
   GetSubMatrix Rows(int,int);
   GetSubMatrix Column(int);
   GetSubMatrix Columns(int,int);
   operator real();                            // conversion of 1 x 1 matrix
   virtual LogAndSign LogDeterminant();
   virtual real SumSquare();
   virtual real SumAbsoluteValue();
   virtual real MaximumAbsoluteValue();
   virtual real Trace();
   real Norm1();
   real NormInfinity();

   friend GeneralMatrix;
   friend Matrix;
   friend nricMatrix;
   friend RowVector;
   friend ColumnVector;
   friend SymmetricMatrix;
   friend UpperTriangularMatrix;
   friend LowerTriangularMatrix;
   friend DiagonalMatrix;
   friend CroutMatrix;
   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 ReturnMatrix;
};


/******************************* 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(int);                          // 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(GeneralMatrix*);              // used by = and initialise
#ifndef __ZTC__
   void GetMatrixC(const GeneralMatrix*);       // used by = and initialise
#endif
   void Eq(BaseMatrix&, MatrixType);            // used by =
   int search(const GeneralMatrix*) const;
   virtual GeneralMatrix* Transpose(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
   BOOL reuse();                                // TRUE if tag allows reuse
   void Protect() { tag=-1; }                   // can't delete or reuse
   int Tag() const { return tag; }
   BOOL 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<<(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();
   real SumAbsoluteValue();
   real MaximumAbsoluteValue();
   LogAndSign LogDeterminant();
   ConstMatrix c() const;                       // to access constant matrices
   void CheckStore() const;                     // check store is non-zero

   friend Matrix;
   friend nricMatrix;
   friend SymmetricMatrix;
   friend UpperTriangularMatrix;
   friend LowerTriangularMatrix;
   friend DiagonalMatrix;
   friend CroutMatrix;
   friend RowVector;
   friend ColumnVector;
   friend BaseMatrix;
   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 ReturnMatrix;
};

class Matrix : public GeneralMatrix             // usual rectangular matrix
{
public:
   Matrix() {}
   Matrix(int, int);                            // standard declaration
   Matrix(BaseMatrix&);                         // evaluate BaseMatrix
   void operator=(BaseMatrix&);
   void operator=(real f) { GeneralMatrix::operator=(f); }
   MatrixType Type() const;
   real& operator()(int, int);                  // access element
   real& element(int, int);                     // access element
   Matrix(Matrix& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   real operator()(int, int) const;             // access element
   Matrix(const Matrix& gm) { GetMatrixC(&gm); }
#endif
   GeneralMatrix* MakeSolver();
   real Trace();
   void GetRow(MatrixRowCol&);
   void GetCol(MatrixRowCol&);
   void RestoreCol(MatrixRowCol&);
   void NextRow(MatrixRowCol&);
   void NextCol(MatrixRowCol&);
   void ReDimension(int,int);                   // change dimensions
};

class nricMatrix : public Matrix                // for use with Numerical
                                                // Recipes in C
{
   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(BaseMatrix& bm)                   // evaluate BaseMatrix
      :  Matrix(bm) { MakeRowPointer(); }
   void operator=(BaseMatrix& bm)
      { DeleteRowPointer(); Matrix::operator=(bm); MakeRowPointer(); }
   void operator=(real f) { GeneralMatrix::operator=(f); }
   void operator<<(BaseMatrix& X)
      { DeleteRowPointer(); Eq(X,this->Type()); MakeRowPointer(); }
   nricMatrix(nricMatrix& gm) { GetMatrix(&gm); MakeRowPointer(); }
#ifndef __ZTC__
   nricMatrix(const nricMatrix& gm) { GetMatrixC(&gm); MakeRowPointer(); }
#endif
   void ReDimension(int m, int n)               // change dimensions
      { DeleteRowPointer(); Matrix::ReDimension(m,n); MakeRowPointer(); }
   ~nricMatrix() { DeleteRowPointer(); }
#ifndef __ZTC__
   operator real**() const { CheckStore(); return row_pointer-1; }
#endif
};

class SymmetricMatrix : public GeneralMatrix
{
public:
   SymmetricMatrix() {}
   SymmetricMatrix(int);
   SymmetricMatrix(BaseMatrix&);
   void operator=(BaseMatrix&);
   void operator=(real f) { GeneralMatrix::operator=(f); }
   real& operator()(int, int);                  // access element
   real& element(int, int);                     // access element
   MatrixType Type() const;
   SymmetricMatrix(SymmetricMatrix& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   real operator()(int, int) const;             // access element
   SymmetricMatrix(const SymmetricMatrix& gm) { GetMatrixC(&gm); }
#endif
   real SumSquare();
   real SumAbsoluteValue();
   real Trace();
   void GetRow(MatrixRowCol&);
   void GetCol(MatrixRowCol&);
   GeneralMatrix* Transpose(MatrixType);
   void ReDimension(int);                       // change dimensions
};

class UpperTriangularMatrix : public GeneralMatrix
{
public:
   UpperTriangularMatrix() {}
   UpperTriangularMatrix(int);
   void operator=(BaseMatrix&);
   UpperTriangularMatrix(BaseMatrix&);
   UpperTriangularMatrix(UpperTriangularMatrix& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   real operator()(int, int) const;             // access element
   UpperTriangularMatrix(const UpperTriangularMatrix& gm) { GetMatrixC(&gm); }
#endif
   void operator=(real f) { GeneralMatrix::operator=(f); }
   real& operator()(int, int);                  // access element
   real& element(int, int);                     // access element
   MatrixType Type() const;
   GeneralMatrix* MakeSolver() { return this; } // for solving
   void Solver(MatrixRowCol&, const MatrixRowCol&);
   LogAndSign LogDeterminant();
   real Trace();
   void GetRow(MatrixRowCol&);
   void GetCol(MatrixRowCol&);
   void RestoreCol(MatrixRowCol&);
   void NextRow(MatrixRowCol&);
   void ReDimension(int);                       // change dimensions
};

class LowerTriangularMatrix : public GeneralMatrix
{
public:
   LowerTriangularMatrix() {}
   LowerTriangularMatrix(int);
   LowerTriangularMatrix(LowerTriangularMatrix& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   real operator()(int, int) const;             // access element
   LowerTriangularMatrix(const LowerTriangularMatrix& gm) { GetMatrixC(&gm); }
#endif
   LowerTriangularMatrix(BaseMatrix& M);
   void operator=(BaseMatrix&);
   void operator=(real f) { GeneralMatrix::operator=(f); }
   real& operator()(int, int);                  // access element
   real& element(int, int);                     // access element
   MatrixType Type() const;
   GeneralMatrix* MakeSolver() { return this; } // for solving
   void Solver(MatrixRowCol&, const MatrixRowCol&);
   LogAndSign LogDeterminant();
   real Trace();
   void GetRow(MatrixRowCol&);
   void GetCol(MatrixRowCol&);
   void RestoreCol(MatrixRowCol&);
   void NextRow(MatrixRowCol&);
   void ReDimension(int);                       // change dimensions
};

class DiagonalMatrix : public GeneralMatrix
{
public:
   DiagonalMatrix() {}
   DiagonalMatrix(int);
   DiagonalMatrix(BaseMatrix&);
   DiagonalMatrix(DiagonalMatrix& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   real operator()(int, int) const;             // access element
   real operator()(int) const;
   DiagonalMatrix(const DiagonalMatrix& gm) { GetMatrixC(&gm); }
#endif
   void operator=(BaseMatrix&);
   void operator=(real f) { GeneralMatrix::operator=(f); }
   real& operator()(int, int);                  // access element
   real& operator()(int);                       // access element
   real& element(int, int);                     // access element
   real& element(int);                          // access element
   MatrixType Type() const;

   LogAndSign LogDeterminant();
   real Trace();
   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(MatrixType);
   void ReDimension(int);                       // change dimensions
#ifndef __ZTC__
   operator real*() const
      { CheckStore(); return store-1; }         // for use by NRIC
#endif
};

class RowVector : public Matrix
{
public:
   RowVector() {}
   RowVector(int n) : Matrix(1,n) {}
   RowVector(BaseMatrix&);
   RowVector(RowVector& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   real operator()(int) const;                  // access element
   RowVector(const RowVector& gm) {GetMatrixC(&gm); }
#endif
   void operator=(BaseMatrix&);
   void operator=(real f) { GeneralMatrix::operator=(f); }
   real& operator()(int);                       // access element
   real& element(int);                          // access element
   MatrixType Type() const;
   void GetCol(MatrixRowCol&);
   void NextCol(MatrixRowCol&);
   void RestoreCol(MatrixRowCol&);
   GeneralMatrix* Transpose(MatrixType);
   void ReDimension(int);                       // change dimensions
#ifndef __ZTC__
   operator real*() const
      { CheckStore(); return store-1; }         // for use by NRIC
#endif
};

class ColumnVector : public Matrix
{
public:
   ColumnVector() {}
   ColumnVector(int n) : Matrix(n,1) {}
   ColumnVector(BaseMatrix&);
   ColumnVector(ColumnVector& gm) { GetMatrix(&gm); }
#ifndef __ZTC__
   real operator()(int) const;                  // access element
   ColumnVector(const ColumnVector& gm) { GetMatrixC(&gm); }
#endif
   void operator=(BaseMatrix&);
   void operator=(real f) { GeneralMatrix::operator=(f); }
   real& operator()(int);                       // access element
   real& element(int);                          // access element
   MatrixType Type() const;
   GeneralMatrix* Transpose(MatrixType);
   void ReDimension(int);                       // change dimensions
#ifndef __ZTC__
   operator real*() const
      { CheckStore(); return store-1; }         // for use by NRIC
#endif
};

class CroutMatrix : public GeneralMatrix        // for LU decomposition
{
   int* indx;
   BOOL d;
   void ludcmp();
public:
   CroutMatrix(BaseMatrix&);
   MatrixType Type() const;
   void lubksb(real*, int=0);
   ~CroutMatrix() { delete [nrows] indx; }
   GeneralMatrix* MakeSolver() { return this; } // for solving
   LogAndSign LogDeterminant();
   void Solver(MatrixRowCol&, const MatrixRowCol&);
   void GetRow(MatrixRowCol&) { MatrixError("GetRow not defined for Crout"); }
   void GetCol(MatrixRowCol&) { MatrixError("GetCol not defined for Crout"); }
   void operator=(BaseMatrix&)
      { MatrixError("operator= not defined for Crout"); }
};


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

class AddedMatrix : public BaseMatrix
{
protected:
   BaseMatrix* bm1;                             // pointers to summands
   BaseMatrix* bm2;
   AddedMatrix(BaseMatrix* bm1x, BaseMatrix* bm2x) : bm1(bm1x),bm2(bm2x) {}
   int search(const GeneralMatrix*) const;
   int NrowsV() const;
   int NcolsV() const;
private:
   MatrixType Type() const;
   friend BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
};

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

class SolvedMatrix : public AddedMatrix
{
   SolvedMatrix(BaseMatrix* bm1x, BaseMatrix* bm2x)
      : AddedMatrix(bm1x,bm2x) {}
   MatrixType Type() const;
   friend BaseMatrix;
   friend InvertedMatrix;                        // for operator*
public:
   GeneralMatrix* Evaluate(MatrixType);
};

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

class ShiftedMatrix : public BaseMatrix
{
protected:
   real f;
   BaseMatrix* bm;
   ShiftedMatrix(BaseMatrix* bmx, real fx) : bm(bmx),f(fx) {}
   int search(const GeneralMatrix*) const;
   int NrowsV() const;
   int NcolsV() const;
private:
   MatrixType Type() const;
   friend BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
};

class ScaledMatrix : public ShiftedMatrix
{
   ScaledMatrix(BaseMatrix* bmx, real fx) : ShiftedMatrix(bmx,fx) {}
   MatrixType Type() const;
   friend BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
};

class NegatedMatrix : public BaseMatrix
{
protected:
   BaseMatrix* bm;
   NegatedMatrix(BaseMatrix* bmx) : bm(bmx) {}
   int search(const GeneralMatrix*) const;
   int NrowsV() const;
   int NcolsV() const;
private:
   MatrixType Type() const;
   friend BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
};

class TransposedMatrix : public NegatedMatrix
{
   TransposedMatrix(BaseMatrix* bmx) : NegatedMatrix(bmx) {}
   int NrowsV() const;
   int NcolsV() const;
   MatrixType Type() const;
   friend BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
};

class InvertedMatrix : public NegatedMatrix
{
   InvertedMatrix(BaseMatrix* bmx) : NegatedMatrix(bmx) {}
   MatrixType Type() const;
public:
   SolvedMatrix operator*(BaseMatrix&);         // inverse(A) * B
   friend BaseMatrix;
   GeneralMatrix* Evaluate(MatrixType);
};

class RowedMatrix : public NegatedMatrix
{
   MatrixType Type() const;
   int NrowsV() const;
   int NcolsV() const;
   RowedMatrix(BaseMatrix* bmx) : NegatedMatrix(bmx) {}
   friend BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
};

class ColedMatrix : public NegatedMatrix
{
   MatrixType Type() const;
   int NrowsV() const;
   int NcolsV() const;
   ColedMatrix(BaseMatrix* bmx) : NegatedMatrix(bmx) {}
   friend BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
};

class DiagedMatrix : public NegatedMatrix
{
   MatrixType Type() const;
   int NrowsV() const;
   int NcolsV() const;
   DiagedMatrix(BaseMatrix* bmx) : NegatedMatrix(bmx) {}
   friend BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
};

class MatedMatrix : public NegatedMatrix
{
   int nr, nc;
   MatrixType Type() const;
   int NrowsV() const;
   int NcolsV() const;
   MatedMatrix(BaseMatrix* bmx, int nrx, int ncx)
      : NegatedMatrix(bmx), nr(nrx), nc(ncx) {}
   friend BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
};

class ConstMatrix : public BaseMatrix
{
   const GeneralMatrix* cgm;
   MatrixType Type() const;
   int NrowsV() const;
   int NcolsV() const;
   int search(const GeneralMatrix*) const;
   ConstMatrix(const GeneralMatrix* cgmx) : cgm(cgmx) {}
   friend BaseMatrix;
   friend GeneralMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
};

class ReturnMatrix : public BaseMatrix    // for matrix return
{
   GeneralMatrix* gm;
   MatrixType Type() const;
   int NrowsV() const;
   int NcolsV() const;
   int search(const GeneralMatrix*) const;
public:
   GeneralMatrix* Evaluate(MatrixType);
   friend BaseMatrix;
   ReturnMatrix(const ReturnMatrix& tm) : gm(tm.gm) {}
   ReturnMatrix(GeneralMatrix&);
};


/**************************** submatrices ******************************/

class GetSubMatrix : public NegatedMatrix
{
   int row_skip;
   int row_number;
   int col_skip;
   int col_number;
   MatrixType mt;

   GetSubMatrix
      (BaseMatrix* bmx, int rs, int rn, int cs, int cn, MatrixType mtx)
      : NegatedMatrix(bmx),
      row_skip(rs), row_number(rn), col_skip(cs), col_number(cn), mt(mtx) {}
   GetSubMatrix(const GetSubMatrix& g)
      : NegatedMatrix(g.bm), row_skip(g.row_skip), row_number(g.row_number),
      col_skip(g.col_skip), col_number(g.col_number), mt(g.mt) {}
   MatrixType Type() const;
   int NrowsV() const;
   int NcolsV() const;
   friend BaseMatrix;
public:
   GeneralMatrix* Evaluate(MatrixType);
   void operator=(BaseMatrix&)
      { MatrixError("= not defined for submatrix; use <<"); }
   void operator<<(BaseMatrix&);
   void operator<<(const real*);                // copy from array
   void operator<<(real);                       // copy from constant
   void Inject(const GeneralMatrix&);           // copy stored els only
};

/***************************** functions ***********************************/


inline LogAndSign LogDeterminant(BaseMatrix& B) { return B.LogDeterminant(); }
inline real SumSquare(BaseMatrix& B) { return B.SumSquare(); }
inline real Trace(BaseMatrix& B) { return B.Trace(); }
inline real SumAbsoluteValue(BaseMatrix& B) { return B.SumAbsoluteValue(); }
inline real MaximumAbsoluteValue(BaseMatrix& B)
   { return B.MaximumAbsoluteValue(); }
inline real Norm1(BaseMatrix& B) { return B.Norm1(); }
inline real Norm1(RowVector& RV) { return RV.MaximumAbsoluteValue(); }
inline real NormInfinity(BaseMatrix& B) { return B.NormInfinity(); }
inline real NormInfinity(ColumnVector& CV)
   { return CV.MaximumAbsoluteValue(); } 

#endif
