home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #1 / NN_1993_1.iso / spool / comp / sources / misc / 4251 / newmat.h < prev    next >
Encoding:
C/C++ Source or Header  |  1993-01-11  |  46.3 KB  |  1,343 lines
  1. //$$ newmat.h           definition file for new version of matrix package
  2.  
  3. // Copyright (C) 1991,2,3: R B Davies
  4.  
  5. #ifndef MATRIX_LIB
  6. #define MATRIX_LIB 0
  7.  
  8. #ifdef NO_LONG_NAMES
  9. #define UpperTriangularMatrix UTMatrix
  10. #define LowerTriangularMatrix LTMatrix
  11. #define SymmetricMatrix SMatrix
  12. #define DiagonalMatrix DMatrix
  13. #endif
  14.  
  15. #ifndef TEMPS_DESTROYED_QUICKLY
  16. #define ReturnMatrix ReturnMatrixX
  17. #else
  18. #define ReturnMatrix ReturnMatrixX&
  19. #endif
  20.  
  21. #include "boolean.h"
  22. #include "except.h"
  23.  
  24. /**************************** general utilities ****************************/
  25.  
  26. class GeneralMatrix;
  27. void MatrixErrorNoSpace(void*);                 // no space handler
  28.  
  29. class LogAndSign
  30. // Return from LogDeterminant function
  31. //    - value of the log plus the sign (+, - or 0)
  32. {
  33.    Real log_value;
  34.    int sign;
  35. public:
  36.    LogAndSign() { log_value=0.0; sign=1; }
  37.    LogAndSign(Real);
  38.    void operator*=(Real);
  39.    void ChangeSign() { sign = -sign; }
  40.    Real LogValue() const { return log_value; }
  41.    int Sign() const { return sign; }
  42.    Real Value() const;
  43.    FREE_CHECK(LogAndSign)
  44. };
  45.  
  46. // the following class is for counting the number of times a piece of code
  47. // is executed. It is used for locating any code not executed by test
  48. // routines. Use turbo GREP locate all places this code is called and
  49. // check which ones are not accessed.
  50. // Somewhat implementation dependent as it relies on "cout" still being
  51. // present when ExeCounter objects are destructed.
  52.  
  53. class ExeCounter
  54. {
  55.    int line;                                    // code line number
  56.    int fileid;                                  // file identifier
  57.    long nexe;                                   // number of executions
  58.    static int nreports;                         // number of reports
  59. public:
  60.    ExeCounter(int,int);
  61.    void operator++() { nexe++; }
  62.    ~ExeCounter();                               // prints out reports
  63. };
  64.  
  65.  
  66. /************** Class to show whether to check for loss of data ************/
  67.  
  68. class MatrixConversionCheck : public Janitor
  69. {
  70.    static Boolean DoCheck;
  71. public:
  72.    MatrixConversionCheck() { DoCheck=TRUE; }          // turn check on
  73.    ~MatrixConversionCheck() { DoCheck=FALSE; }        // turn check off
  74.    void CleanUp() { DoCheck=FALSE; }
  75.    static Boolean IsOn() { return DoCheck; }
  76.    static void DataLoss();
  77. };
  78.  
  79. /**************************** class MatrixType *****************************/
  80.  
  81. // Is used for finding the type of a matrix resulting from the binary operations
  82. // +, -, * and identifying what conversions are permissible.
  83. // This class must be updated when new matrix types are added.
  84.  
  85. class GeneralMatrix;                            // defined later
  86. class BaseMatrix;                               // defined later
  87. class MatrixInput;                              // defined later
  88.  
  89. class MatrixType
  90. {
  91. public:
  92.    enum Attribute {  Valid     = 1,
  93.                      Symmetric = 2,
  94.                      Band      = 4,
  95.                      Upper     = 8,
  96.                      Lower     = 16,
  97.                      LUDeco    = 32 };
  98.  
  99.    enum            { US = 0,
  100.                      UT = Valid + Upper,
  101.                      LT = Valid + Lower,
  102.                      Rt = Valid,
  103.                      Sm = Valid + Symmetric,
  104.                      Dg = Valid + Band + Lower + Upper + Symmetric,
  105.              RV = Valid,     //   don't separate out
  106.              CV = Valid,     //   vectors
  107.              BM = Valid + Band,
  108.              UB = Valid + Band + Upper,
  109.              LB = Valid + Band + Lower,
  110.              SB = Valid + Band + Symmetric,
  111.              Ct = Valid + LUDeco,
  112.              BC = Valid + Band + LUDeco,
  113.                    };
  114.  
  115.  
  116.    static nTypes() { return 9; }               // number of different types
  117.                            // exclude Ct, US, BC
  118. public:
  119.    int attribute;
  120. public:
  121.    MatrixType () {}
  122.    MatrixType (int i) : attribute(i) {}
  123.    int operator+() const { return attribute; }
  124.    MatrixType operator+(MatrixType mt) const
  125.       { return MatrixType(attribute & mt.attribute); }
  126.    MatrixType operator*(const MatrixType&) const;
  127.    Boolean operator>=(MatrixType mt) const
  128.       { return ( attribute & mt.attribute ) == attribute; }
  129.    Boolean operator==(MatrixType t) const
  130.       { return (attribute == t.attribute); }
  131.    Boolean operator!=(MatrixType t) const
  132.       { return (attribute != t.attribute); }
  133.    Boolean operator!() const { return (attribute & Valid) == 0; }
  134.    MatrixType i() const                        // type of inverse
  135.       { return MatrixType(attribute & ~(Band+LUDeco)); }
  136.    MatrixType t() const;                       // type of transpose
  137.    MatrixType AddEqualEl() const               // Add constant to matrix
  138.       { return MatrixType(attribute & (Valid + Symmetric)); }
  139.    MatrixType MultRHS() const;                 // type for rhs of multiply
  140.    MatrixType sub() const                      // type of submatrix
  141.       { return MatrixType(attribute & Valid); }
  142.    MatrixType ssub() const                     // type of sym submatrix
  143.       { return MatrixType(attribute); }        // not for selection matrix
  144.    GeneralMatrix* New() const;                 // new matrix of given type
  145.    GeneralMatrix* New(int,int,BaseMatrix*) const;
  146.                                                // new matrix of given type
  147.    operator char*() const;                     // to print type
  148.    friend Boolean Rectangular(MatrixType a, MatrixType b, MatrixType c)
  149.       { return ((a.attribute | b.attribute | c.attribute)
  150.           & ~MatrixType::Valid) == 0; }
  151.    friend Boolean Compare(const MatrixType&, MatrixType&);
  152.                                                // compare and check conv.
  153.    FREE_CHECK(MatrixType)
  154. };
  155.  
  156. void TestTypeAdd();                            // test +
  157. void TestTypeMult();                           // test *
  158. void TestTypeOrder();                          // test >=
  159.  
  160.  
  161. /************************* class MatrixBandWidth ***********************/
  162.  
  163. class MatrixBandWidth
  164. {
  165. public:
  166.    int lower;
  167.    int upper;
  168.    MatrixBandWidth(const int l, const int u) : lower(l), upper (u) {}
  169.    MatrixBandWidth(const int i) : lower(i), upper(i) {}
  170.    MatrixBandWidth operator+(const MatrixBandWidth&) const;
  171.    MatrixBandWidth operator*(const MatrixBandWidth&) const;
  172.    MatrixBandWidth t() const { return MatrixBandWidth(upper,lower); }
  173.    Boolean operator==(const MatrixBandWidth& bw) const
  174.       { return (lower == bw.lower) && (upper == bw.upper); }
  175.    FREE_CHECK(MatrixBandWidth)
  176. };
  177.  
  178.  
  179. /*********************** Array length specifier ************************/
  180.  
  181. // This class is introduced to avoid constructors such as
  182. //   ColumnVector(int)
  183. // being used for conversions
  184.  
  185. class ArrayLengthSpecifier
  186. {
  187.    int value;
  188. public:
  189.    int Value() const { return value; }
  190.    ArrayLengthSpecifier(int l) : value(l) {}
  191.    FREE_CHECK(ArrayLengthSpecifier)
  192. };
  193.  
  194. /*************************** Matrix routines ***************************/
  195.  
  196.  
  197. class MatrixRowCol;                             // defined later
  198. class MatrixRow;
  199. class MatrixCol;
  200.  
  201. class GeneralMatrix;                            // defined later
  202. class AddedMatrix;
  203. class MultipliedMatrix;
  204. class SubtractedMatrix;
  205. class SolvedMatrix;
  206. class ShiftedMatrix;
  207. class ScaledMatrix;
  208. class TransposedMatrix;
  209. class NegatedMatrix;
  210. class InvertedMatrix;
  211. class RowedMatrix;
  212. class ColedMatrix;
  213. class DiagedMatrix;
  214. class MatedMatrix;
  215. class GetSubMatrix;
  216. class ConstMatrix;
  217. class ReturnMatrixX;
  218. class Matrix;
  219. class nricMatrix;
  220. class RowVector;
  221. class ColumnVector;
  222. class SymmetricMatrix;
  223. class UpperTriangularMatrix;
  224. class LowerTriangularMatrix;
  225. class DiagonalMatrix;
  226. class CroutMatrix;
  227. class BandMatrix;
  228. class LowerBandMatrix;
  229. class UpperBandMatrix;
  230. class SymmetricBandMatrix;
  231. class LinearEquationSolver;
  232.  
  233.  
  234.  
  235. static MatrixType MatrixTypeUnSp(MatrixType::US);
  236.                         // AT&T needs this
  237.  
  238. class BaseMatrix : public Janitor               // base of all matrix classes
  239. {
  240. protected:
  241.    virtual int search(const BaseMatrix*) const = 0;
  242.                         // count number of times matrix
  243.                         // is referred to
  244. public:
  245. #ifndef __GNUG__
  246.    virtual GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp) = 0;
  247.                         // evaluate temporary
  248. #else
  249.    virtual GeneralMatrix* Evaluate(MatrixType mt) = 0;
  250.    GeneralMatrix* Evaluate() { return Evaluate(MatrixTypeUnSp); }
  251. #endif
  252. #ifndef TEMPS_DESTROYED_QUICKLY
  253.    AddedMatrix operator+(const BaseMatrix&) const;    // results of operations
  254.    MultipliedMatrix operator*(const BaseMatrix&) const;
  255.    SubtractedMatrix operator-(const BaseMatrix&) const;
  256.    ShiftedMatrix operator+(Real) const;
  257.    ScaledMatrix operator*(Real) const;
  258.    ScaledMatrix operator/(Real) const;
  259.    ShiftedMatrix operator-(Real) const;
  260.    TransposedMatrix t() const;
  261. //   TransposedMatrix t;
  262.    NegatedMatrix operator-() const;                   // change sign of elements
  263.    InvertedMatrix i() const;
  264. //   InvertedMatrix i;
  265.    RowedMatrix AsRow() const;
  266.    ColedMatrix AsColumn() const;
  267.    DiagedMatrix AsDiagonal() const;
  268.    MatedMatrix AsMatrix(int,int) const;
  269.    GetSubMatrix SubMatrix(int,int,int,int) const;
  270.    GetSubMatrix SymSubMatrix(int,int) const;
  271.    GetSubMatrix Row(int) const;
  272.    GetSubMatrix Rows(int,int) const;
  273.    GetSubMatrix Column(int) const;
  274.    GetSubMatrix Columns(int,int) const;
  275. #else
  276.    AddedMatrix& operator+(const BaseMatrix&) const;    // results of operations
  277.    MultipliedMatrix& operator*(const BaseMatrix&) const;
  278.    SubtractedMatrix& operator-(const BaseMatrix&) const;
  279.    ShiftedMatrix& operator+(Real) const;
  280.    ScaledMatrix& operator*(Real) const;
  281.    ScaledMatrix& operator/(Real) const;
  282.    ShiftedMatrix& operator-(Real) const;
  283.    TransposedMatrix& t() const;
  284. //   TransposedMatrix& t;
  285.    NegatedMatrix& operator-() const;                   // change sign of elements
  286.    InvertedMatrix& i() const;
  287. //   InvertedMatrix& i;
  288.    RowedMatrix& AsRow() const;
  289.    ColedMatrix& AsColumn() const;
  290.    DiagedMatrix& AsDiagonal() const;
  291.    MatedMatrix& AsMatrix(int,int) const;
  292.    GetSubMatrix& SubMatrix(int,int,int,int) const;
  293.    GetSubMatrix& SymSubMatrix(int,int) const;
  294.    GetSubMatrix& Row(int) const;
  295.    GetSubMatrix& Rows(int,int) const;
  296.    GetSubMatrix& Column(int) const;
  297.    GetSubMatrix& Columns(int,int) const;
  298. #endif
  299.    Real AsScalar() const;                      // conversion of 1 x 1 matrix
  300.    virtual LogAndSign LogDeterminant() const;
  301.    virtual Real SumSquare() const;
  302.    virtual Real SumAbsoluteValue() const;
  303.    virtual Real MaximumAbsoluteValue() const;
  304.    virtual Real Trace() const;
  305.    Real Norm1() const;
  306.    Real NormInfinity() const;
  307.    virtual MatrixBandWidth BandWidth() const;  // bandwidths of band matrix
  308.    virtual void CleanUp() {}                     // to clear store
  309. //protected:
  310. //   BaseMatrix() : t(this), i(this) {}
  311.  
  312.    friend class GeneralMatrix;
  313.    friend class Matrix;
  314.    friend class nricMatrix;
  315.    friend class RowVector;
  316.    friend class ColumnVector;
  317.    friend class SymmetricMatrix;
  318.    friend class UpperTriangularMatrix;
  319.    friend class LowerTriangularMatrix;
  320.    friend class DiagonalMatrix;
  321.    friend class CroutMatrix;
  322.    friend class BandMatrix;
  323.    friend class LowerBandMatrix;
  324.    friend class UpperBandMatrix;
  325.    friend class SymmetricBandMatrix;
  326.    friend class AddedMatrix;
  327.    friend class MultipliedMatrix;
  328.    friend class SubtractedMatrix;
  329.    friend class SolvedMatrix;
  330.    friend class ShiftedMatrix;
  331.    friend class ScaledMatrix;
  332.    friend class TransposedMatrix;
  333.    friend class NegatedMatrix;
  334.    friend class InvertedMatrix;
  335.    friend class RowedMatrix;
  336.    friend class ColedMatrix;
  337.    friend class DiagedMatrix;
  338.    friend class MatedMatrix;
  339.    friend class GetSubMatrix;
  340.    friend class ConstMatrix;
  341.    friend class ReturnMatrixX;
  342.    friend class LinearEquationSolver;
  343.    NEW_DELETE(BaseMatrix)
  344. };
  345.  
  346.  
  347. /******************************* working classes **************************/
  348.  
  349. class GeneralMatrix : public BaseMatrix         // declarable matrix types
  350. {
  351.    virtual GeneralMatrix* Image() const;        // copy of matrix
  352. protected:
  353.    int tag;                                     // shows whether can reuse
  354.    int nrows, ncols;                            // dimensions
  355.    int storage;                                 // total store required
  356.    Real* store;                                 // point to store (0=not set)
  357.    GeneralMatrix();                             // initialise with no store
  358.    GeneralMatrix(ArrayLengthSpecifier);         // constructor getting store
  359.    void Add(GeneralMatrix*, Real);              // sum of GM and Real
  360.    void Add(Real);                              // add Real to this
  361.    void Multiply(GeneralMatrix*, Real);         // product of GM and Real
  362.    void Multiply(Real);                         // multiply this by Real
  363.    void Negate(GeneralMatrix*);                 // change sign
  364.    void Negate();                               // change sign
  365.    void operator=(Real);                        // set matrix to constant
  366.    Real* GetStore();                            // get store or copy
  367.    GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
  368.                                                 // temporarily access store
  369.    void GetMatrix(const GeneralMatrix*);        // used by = and initialise
  370.    void Eq(const BaseMatrix&, MatrixType);      // used by =
  371.    int search(const BaseMatrix*) const;
  372.    virtual GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  373.    void CheckConversion(const BaseMatrix&);     // check conversion OK
  374.    void ReDimension(int, int, int);             // change dimensions
  375. public:
  376.    GeneralMatrix* Evaluate(MatrixType);
  377.    virtual MatrixType Type() const = 0;       // type of a matrix
  378.    int Nrows() const { return nrows; }          // get dimensions
  379.    int Ncols() const { return ncols; }
  380.    int Storage() const { return storage; }
  381.    Real* Store() const { return store; }
  382.    virtual ~GeneralMatrix();                    // delete store if set
  383.    void tDelete();                              // delete if tag permits
  384.    Boolean reuse();                                // TRUE if tag allows reuse
  385.    void Protect() { tag=-1; }                   // can't delete or reuse
  386.    int Tag() const { return tag; }
  387.    Boolean IsZero() const;                         // test matrix has all zeros
  388.    void Release() { tag=1; }                    // del store after next use
  389.    void Release(int t) { tag=t; }               // del store after t accesses
  390.    void ReleaseAndDelete() { tag=0; }           // delete matrix after use
  391.    void operator<<(const Real*);                // assignment from an array
  392.    void operator<<(const BaseMatrix& X) { Eq(X,this->Type()); }
  393.                                                 // = without checking type
  394.    void Inject(const GeneralMatrix&);           // copy stored els only
  395.    virtual GeneralMatrix* MakeSolver();         // for solving
  396.    virtual void Solver(MatrixRowCol&, const MatrixRowCol&) {}
  397.    virtual void GetRow(MatrixRowCol&) = 0;      // Get matrix row
  398.    virtual void RestoreRow(MatrixRowCol&) {}    // Restore matrix row
  399.    virtual void NextRow(MatrixRowCol&);         // Go to next row
  400.    virtual void GetCol(MatrixRowCol&) = 0;      // Get matrix col
  401.    virtual void RestoreCol(MatrixRowCol&) {}    // Restore matrix col
  402.    virtual void NextCol(MatrixRowCol&);         // Go to next col
  403.    Real SumSquare() const;
  404.    Real SumAbsoluteValue() const;
  405.    Real MaximumAbsoluteValue() const;
  406.    LogAndSign LogDeterminant() const;
  407. #ifndef TEMPS_DESTROYED_QUICKLY
  408.    ConstMatrix c() const;                       // to access constant matrices
  409. #else
  410.    ConstMatrix& c() const;                      // to access constant matrices
  411. #endif
  412.    void CheckStore() const;                     // check store is non-zero
  413.    virtual void SetParameters(const GeneralMatrix*) {}
  414.                                                 // set parameters in GetMatrix
  415. #ifndef TEMPS_DESTROYED_QUICKLY
  416.    operator ReturnMatrixX() const;              // for building a ReturnMatrix
  417. #else
  418.    operator ReturnMatrixX&() const;             // for building a ReturnMatrix
  419. #endif
  420.    MatrixInput operator<<(Real);                // for loading a list
  421.    void CleanUp();                              // to clear store
  422.  
  423.    friend class Matrix;
  424.    friend class nricMatrix;
  425.    friend class SymmetricMatrix;
  426.    friend class UpperTriangularMatrix;
  427.    friend class LowerTriangularMatrix;
  428.    friend class DiagonalMatrix;
  429.    friend class CroutMatrix;
  430.    friend class RowVector;
  431.    friend class ColumnVector;
  432.    friend class BandMatrix;
  433.    friend class LowerBandMatrix;
  434.    friend class UpperBandMatrix;
  435.    friend class SymmetricBandMatrix;
  436.    friend class BaseMatrix;
  437.    friend class AddedMatrix;
  438.    friend class MultipliedMatrix;
  439.    friend class SubtractedMatrix;
  440.    friend class SolvedMatrix;
  441.    friend class ShiftedMatrix;
  442.    friend class ScaledMatrix;
  443.    friend class TransposedMatrix;
  444.    friend class NegatedMatrix;
  445.    friend class InvertedMatrix;
  446.    friend class RowedMatrix;
  447.    friend class ColedMatrix;
  448.    friend class DiagedMatrix;
  449.    friend class MatedMatrix;
  450.    friend class GetSubMatrix;
  451.    friend class ConstMatrix;
  452.    friend class ReturnMatrixX;
  453.    friend class LinearEquationSolver;
  454.    NEW_DELETE(GeneralMatrix)
  455. };
  456.  
  457. class Matrix : public GeneralMatrix             // usual rectangular matrix
  458. {
  459.    GeneralMatrix* Image() const;                // copy of matrix
  460. public:
  461.    Matrix() {}
  462.    Matrix(int, int);                            // standard declaration
  463.    Matrix(const BaseMatrix&);                   // evaluate BaseMatrix
  464.    void operator=(const BaseMatrix&);
  465.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  466.    void operator=(const Matrix& m) { operator=((const BaseMatrix&)m); }
  467.    MatrixType Type() const;
  468.    Real& operator()(int, int);                  // access element
  469.    Real& element(int, int);                     // access element
  470. #ifdef SETUP_C_SUBSCRIPTS
  471.    Real* operator[](int m) { return store+m*ncols; }
  472.    const Real* operator[](int m) const { return store+m*ncols; }
  473. #endif
  474.    Matrix(const Matrix& gm) { GetMatrix(&gm); }
  475. #ifndef __ZTC__
  476.    Real operator()(int, int) const;             // access element
  477. #endif
  478.    GeneralMatrix* MakeSolver();
  479.    Real Trace() const;
  480.    void GetRow(MatrixRowCol&);
  481.    void GetCol(MatrixRowCol&);
  482.    void RestoreCol(MatrixRowCol&);
  483.    void NextRow(MatrixRowCol&);
  484.    void NextCol(MatrixRowCol&);
  485.    void ReDimension(int,int);                   // change dimensions
  486.    NEW_DELETE(Matrix)
  487. };
  488.  
  489. class nricMatrix : public Matrix                // for use with Numerical
  490.                                                 // Recipes in C
  491. {
  492.    GeneralMatrix* Image() const;                // copy of matrix
  493.    Real** row_pointer;                          // points to rows
  494.    void MakeRowPointer();                       // build rowpointer
  495.    void DeleteRowPointer();
  496. public:
  497.    nricMatrix() {}
  498.    nricMatrix(int m, int n)                     // standard declaration
  499.       :  Matrix(m,n) { MakeRowPointer(); }
  500.    nricMatrix(const BaseMatrix& bm)             // evaluate BaseMatrix
  501.       :  Matrix(bm) { MakeRowPointer(); }
  502.    void operator=(const BaseMatrix& bm)
  503.       { DeleteRowPointer(); Matrix::operator=(bm); MakeRowPointer(); }
  504.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  505.    void operator=(const nricMatrix& m) { operator=((const BaseMatrix&)m); }
  506.    void operator<<(const BaseMatrix& X)
  507.       { DeleteRowPointer(); Eq(X,this->Type()); MakeRowPointer(); }
  508.    nricMatrix(const nricMatrix& gm) { GetMatrix(&gm); MakeRowPointer(); }
  509.    void ReDimension(int m, int n)               // change dimensions
  510.       { DeleteRowPointer(); Matrix::ReDimension(m,n); MakeRowPointer(); }
  511.    ~nricMatrix() { DeleteRowPointer(); }
  512. #ifndef __ZTC__
  513.    Real** nric() const { CheckStore(); return row_pointer-1; }
  514. #endif
  515.    void CleanUp();                                // to clear store
  516.    NEW_DELETE(nricMatrix)
  517. };
  518.  
  519. class SymmetricMatrix : public GeneralMatrix
  520. {
  521.    GeneralMatrix* Image() const;                // copy of matrix
  522. public:
  523.    SymmetricMatrix() {}
  524.    SymmetricMatrix(ArrayLengthSpecifier);
  525.    SymmetricMatrix(const BaseMatrix&);
  526.    void operator=(const BaseMatrix&);
  527.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  528.    void operator=(const SymmetricMatrix& m) { operator=((const BaseMatrix&)m); }
  529.    Real& operator()(int, int);                  // access element
  530.    Real& element(int, int);                     // access element
  531. #ifdef SETUP_C_SUBSCRIPTS
  532.    Real* operator[](int m) { return store+(m*(m+1))/2; }
  533.    const Real* operator[](int m) const { return store+(m*(m+1))/2; }
  534. #endif
  535.    MatrixType Type() const;
  536.    SymmetricMatrix(const SymmetricMatrix& gm) { GetMatrix(&gm); }
  537. #ifndef __ZTC__
  538.    Real operator()(int, int) const;             // access element
  539. #endif
  540.    Real SumSquare() const;
  541.    Real SumAbsoluteValue() const;
  542.    Real Trace() const;
  543.    void GetRow(MatrixRowCol&);
  544.    void GetCol(MatrixRowCol&);
  545.    GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  546.    void ReDimension(int);                       // change dimensions
  547.    NEW_DELETE(SymmetricMatrix)
  548. };
  549.  
  550. class UpperTriangularMatrix : public GeneralMatrix
  551. {
  552.    GeneralMatrix* Image() const;                // copy of matrix
  553. public:
  554.    UpperTriangularMatrix() {}
  555.    UpperTriangularMatrix(ArrayLengthSpecifier);
  556.    void operator=(const BaseMatrix&);
  557.    void operator=(const UpperTriangularMatrix& m)
  558.       { operator=((const BaseMatrix&)m); }
  559.    UpperTriangularMatrix(const BaseMatrix&);
  560.    UpperTriangularMatrix(const UpperTriangularMatrix& gm) { GetMatrix(&gm); }
  561. #ifndef __ZTC__
  562.    Real operator()(int, int) const;             // access element
  563. #endif
  564.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  565.    Real& operator()(int, int);                  // access element
  566.    Real& element(int, int);                     // access element
  567. #ifdef SETUP_C_SUBSCRIPTS
  568.    Real* operator[](int m) { return store+m*ncols-(m*(m+1))/2; }
  569.    const Real* operator[](int m) const { return store+m*ncols-(m*(m+1))/2; }
  570. #endif
  571.    MatrixType Type() const;
  572.    GeneralMatrix* MakeSolver() { return this; } // for solving
  573.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  574.    LogAndSign LogDeterminant() const;
  575.    Real Trace() const;
  576.    void GetRow(MatrixRowCol&);
  577.    void GetCol(MatrixRowCol&);
  578.    void RestoreCol(MatrixRowCol&);
  579.    void NextRow(MatrixRowCol&);
  580.    void ReDimension(int);                       // change dimensions
  581.    NEW_DELETE(UpperTriangularMatrix)
  582. };
  583.  
  584. class LowerTriangularMatrix : public GeneralMatrix
  585. {
  586.    GeneralMatrix* Image() const;                // copy of matrix
  587. public:
  588.    LowerTriangularMatrix() {}
  589.    LowerTriangularMatrix(ArrayLengthSpecifier);
  590.    LowerTriangularMatrix(const LowerTriangularMatrix& gm) { GetMatrix(&gm); }
  591. #ifndef __ZTC__
  592.    Real operator()(int, int) const;             // access element
  593. #endif
  594.    LowerTriangularMatrix(const BaseMatrix& M);
  595.    void operator=(const BaseMatrix&);
  596.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  597.    void operator=(const LowerTriangularMatrix& m)
  598.       { operator=((const BaseMatrix&)m); }
  599.    Real& operator()(int, int);                  // access element
  600.    Real& element(int, int);                     // access element
  601. #ifdef SETUP_C_SUBSCRIPTS
  602.    Real* operator[](int m) { return store+(m*(m+1))/2; }
  603.    const Real* operator[](int m) const { return store+(m*(m+1))/2; }
  604. #endif
  605.    MatrixType Type() const;
  606.    GeneralMatrix* MakeSolver() { return this; } // for solving
  607.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  608.    LogAndSign LogDeterminant() const;
  609.    Real Trace() const;
  610.    void GetRow(MatrixRowCol&);
  611.    void GetCol(MatrixRowCol&);
  612.    void RestoreCol(MatrixRowCol&);
  613.    void NextRow(MatrixRowCol&);
  614.    void ReDimension(int);                       // change dimensions
  615.    NEW_DELETE(LowerTriangularMatrix)
  616. };
  617.  
  618. class DiagonalMatrix : public GeneralMatrix
  619. {
  620.    GeneralMatrix* Image() const;                // copy of matrix
  621. public:
  622.    DiagonalMatrix() {}
  623.    DiagonalMatrix(ArrayLengthSpecifier);
  624.    DiagonalMatrix(const BaseMatrix&);
  625.    DiagonalMatrix(const DiagonalMatrix& gm) { GetMatrix(&gm); }
  626. #ifndef __ZTC__
  627.    Real operator()(int, int) const;             // access element
  628.    Real operator()(int) const;
  629. #endif
  630.    void operator=(const BaseMatrix&);
  631.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  632.    void operator=(const DiagonalMatrix& m) { operator=((const BaseMatrix&)m); }
  633.    Real& operator()(int, int);                  // access element
  634.    Real& operator()(int);                       // access element
  635.    Real& element(int, int);                     // access element
  636.    Real& element(int);                          // access element
  637. #ifdef SETUP_C_SUBSCRIPTS
  638.    Real& operator[](int m) { return store[m]; }
  639.    const Real& operator[](int m) const { return store[m]; }
  640. #endif
  641.    MatrixType Type() const;
  642.  
  643.    LogAndSign LogDeterminant() const;
  644.    Real Trace() const;
  645.    void GetRow(MatrixRowCol&);
  646.    void GetCol(MatrixRowCol&);
  647.    void NextRow(MatrixRowCol&);
  648.    void NextCol(MatrixRowCol&);
  649.    GeneralMatrix* MakeSolver() { return this; } // for solving
  650.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  651.    GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  652.    void ReDimension(int);                       // change dimensions
  653. #ifndef __ZTC__
  654.    Real* nric() const
  655.       { CheckStore(); return store-1; }         // for use by NRIC
  656. #endif
  657.    MatrixBandWidth BandWidth() const;
  658.    NEW_DELETE(DiagonalMatrix)
  659. };
  660.  
  661. class RowVector : public Matrix
  662. {
  663.    GeneralMatrix* Image() const;                // copy of matrix
  664. public:
  665.    RowVector() {}
  666.    RowVector(ArrayLengthSpecifier n) : Matrix(1,n.Value()) {}
  667.    RowVector(const BaseMatrix&);
  668.    RowVector(const RowVector& gm) { GetMatrix(&gm); }
  669. #ifndef __ZTC__
  670.    Real operator()(int) const;                  // access element
  671. #endif
  672.    void operator=(const BaseMatrix&);
  673.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  674.    void operator=(const RowVector& m) { operator=((const BaseMatrix&)m); }
  675.    Real& operator()(int);                       // access element
  676.    Real& element(int);                          // access element
  677. #ifdef SETUP_C_SUBSCRIPTS
  678.    Real& operator[](int m) { return store[m]; }
  679.    const Real& operator[](int m) const { return store[m]; }
  680. #endif
  681.    MatrixType Type() const;
  682.    void GetCol(MatrixRowCol&);
  683.    void NextCol(MatrixRowCol&);
  684.    void RestoreCol(MatrixRowCol&);
  685.    GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  686.    void ReDimension(int);                       // change dimensions
  687. #ifndef __ZTC__
  688.    Real* nric() const
  689.       { CheckStore(); return store-1; }         // for use by NRIC
  690. #endif
  691.    void CleanUp();                                // to clear store
  692.    NEW_DELETE(RowVector)
  693. };
  694.  
  695. class ColumnVector : public Matrix
  696. {
  697.    GeneralMatrix* Image() const;                // copy of matrix
  698. public:
  699.    ColumnVector() {}
  700.    ColumnVector(ArrayLengthSpecifier n) : Matrix(n.Value(),1) {}
  701.    ColumnVector(const BaseMatrix&);
  702.    ColumnVector(const ColumnVector& gm) { GetMatrix(&gm); }
  703. #ifndef __ZTC__
  704.    Real operator()(int) const;                  // access element
  705. #endif
  706.    void operator=(const BaseMatrix&);
  707.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  708.    void operator=(const ColumnVector& m) { operator=((const BaseMatrix&)m); }
  709.    Real& operator()(int);                       // access element
  710.    Real& element(int);                          // access element
  711. #ifdef SETUP_C_SUBSCRIPTS
  712.    Real& operator[](int m) { return store[m]; }
  713.    const Real& operator[](int m) const { return store[m]; }
  714. #endif
  715.    MatrixType Type() const;
  716.    GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  717.    void ReDimension(int);                       // change dimensions
  718. #ifndef __ZTC__
  719.    Real* nric() const
  720.       { CheckStore(); return store-1; }         // for use by NRIC
  721. #endif
  722.    void CleanUp();                                // to clear store
  723.    NEW_DELETE(ColumnVector)
  724. };
  725.  
  726. class CroutMatrix : public GeneralMatrix        // for LU decomposition
  727. {
  728.    int* indx;
  729.    Boolean d;
  730.    Boolean sing;
  731.    void ludcmp();
  732. public:
  733.    CroutMatrix(const BaseMatrix&);
  734.    MatrixType Type() const;
  735.    void lubksb(Real*, int=0);
  736.    ~CroutMatrix();
  737.    GeneralMatrix* MakeSolver() { return this; } // for solving
  738.    LogAndSign LogDeterminant() const;
  739.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  740.    void GetRow(MatrixRowCol&);
  741.    void GetCol(MatrixRowCol&);
  742.    void operator=(const BaseMatrix&);
  743.    void operator=(const CroutMatrix& m) { operator=((const BaseMatrix&)m); }
  744.    void CleanUp();                                // to clear store
  745.    NEW_DELETE(CroutMatrix)
  746. };
  747.  
  748. /******************************* band matrices ***************************/
  749.  
  750. class BandMatrix : public GeneralMatrix         // band matrix
  751. {
  752.    GeneralMatrix* Image() const;                // copy of matrix
  753. protected:
  754.    void CornerClear() const;                    // set unused elements to zero
  755. public:
  756.    int lower, upper;                            // band widths
  757.    BandMatrix() { lower=0; upper=0; CornerClear(); }
  758.    BandMatrix(int n,int lb,int ub) { ReDimension(n,lb,ub); CornerClear(); }
  759.                                                 // standard declaration
  760.    BandMatrix(const BaseMatrix&);               // evaluate BaseMatrix
  761.    void operator=(const BaseMatrix&);
  762.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  763.    void operator=(const BandMatrix& m) { operator=((const BaseMatrix&)m); }
  764.    MatrixType Type() const;
  765.    Real& operator()(int, int);                  // access element
  766.    Real& element(int, int);                     // access element
  767. #ifdef SETUP_C_SUBSCRIPTS
  768.    Real* operator[](int m) { return store+(upper+lower)*m+lower; }
  769.    const Real* operator[](int m) const { return store+(upper+lower)*m+lower; }
  770. #endif
  771.    BandMatrix(const BandMatrix& gm) { GetMatrix(&gm); }
  772. #ifndef __ZTC__
  773.    Real operator()(int, int) const;             // access element
  774. #endif
  775.    LogAndSign LogDeterminant() const;
  776.    GeneralMatrix* MakeSolver();
  777.    Real Trace() const;
  778.    Real SumSquare() const { CornerClear(); return GeneralMatrix::SumSquare(); }
  779.    Real SumAbsoluteValue() const
  780.       { CornerClear(); return GeneralMatrix::SumAbsoluteValue(); }
  781.    Real MaximumAbsoluteValue() const
  782.       { CornerClear(); return GeneralMatrix::MaximumAbsoluteValue(); }
  783.    void GetRow(MatrixRowCol&);
  784.    void GetCol(MatrixRowCol&);
  785.    void RestoreCol(MatrixRowCol&);
  786.    void NextRow(MatrixRowCol&);
  787.    void ReDimension(int, int, int);             // change dimensions
  788.    MatrixBandWidth BandWidth() const;
  789.    void SetParameters(const GeneralMatrix*);
  790.    MatrixInput operator<<(Real);                // will give error
  791.    void operator<<(const Real* r);              // will give error
  792.       // the next is included because Zortech and Borland
  793.       // can't find the copy in GeneralMatrix
  794.    void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
  795.    NEW_DELETE(BandMatrix)
  796. };
  797.  
  798. class UpperBandMatrix : public BandMatrix       // upper band matrix
  799. {
  800.    GeneralMatrix* Image() const;                // copy of matrix
  801. public:
  802.    UpperBandMatrix() {}
  803.    UpperBandMatrix(int n, int ubw)              // standard declaration
  804.       : BandMatrix(n, 0, ubw) {}
  805.    UpperBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
  806.    void operator=(const BaseMatrix&);
  807.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  808.    void operator=(const UpperBandMatrix& m)
  809.       { operator=((const BaseMatrix&)m); }
  810.    MatrixType Type() const;
  811.    UpperBandMatrix(const UpperBandMatrix& gm) { GetMatrix(&gm); }
  812.    GeneralMatrix* MakeSolver() { return this; }
  813.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  814.    LogAndSign LogDeterminant() const;
  815.    void ReDimension(int n,int ubw)              // change dimensions
  816.       { BandMatrix::ReDimension(n,0,ubw); }
  817.    Real& operator()(int, int);
  818.    Real& element(int, int);
  819. #ifdef SETUP_C_SUBSCRIPTS
  820.    Real* operator[](int m) { return store+upper*m; }
  821.    const Real* operator[](int m) const { return store+upper*m; }
  822. #endif
  823.    NEW_DELETE(UpperBandMatrix)
  824. };
  825.  
  826. class LowerBandMatrix : public BandMatrix       // upper band matrix
  827. {
  828.    GeneralMatrix* Image() const;                // copy of matrix
  829. public:
  830.    LowerBandMatrix() {}
  831.    LowerBandMatrix(int n, int lbw)              // standard declaration
  832.       : BandMatrix(n, lbw, 0) {}
  833.    LowerBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
  834.    void operator=(const BaseMatrix&);
  835.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  836.    void operator=(const LowerBandMatrix& m)
  837.       { operator=((const BaseMatrix&)m); }
  838.    MatrixType Type() const;
  839.    LowerBandMatrix(const LowerBandMatrix& gm) { GetMatrix(&gm); }
  840.    GeneralMatrix* MakeSolver() { return this; }
  841.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  842.    LogAndSign LogDeterminant() const;
  843.    void ReDimension(int n,int lbw)             // change dimensions
  844.       { BandMatrix::ReDimension(n,lbw,0); }
  845.    Real& operator()(int, int);
  846.    Real& element(int, int);
  847. #ifdef SETUP_C_SUBSCRIPTS
  848.    Real* operator[](int m) { return store+lower*(m+1); }
  849.    const Real* operator[](int m) const { return store+lower*(m+1); }
  850. #endif
  851.    NEW_DELETE(LowerBandMatrix)
  852. };
  853.  
  854. class SymmetricBandMatrix : public GeneralMatrix
  855. {
  856.    GeneralMatrix* Image() const;                // copy of matrix
  857.    void CornerClear() const;                    // set unused elements to zero
  858. public:
  859.    int lower;                                   // lower band width
  860.    SymmetricBandMatrix() { lower=0; CornerClear(); }
  861.    SymmetricBandMatrix(int n, int lb) { ReDimension(n,lb); CornerClear(); }
  862.    SymmetricBandMatrix(const BaseMatrix&);
  863.    void operator=(const BaseMatrix&);
  864.    void operator=(Real f) { GeneralMatrix::operator=(f); }
  865.    void operator=(const SymmetricBandMatrix& m)
  866.       { operator=((const BaseMatrix&)m); }
  867.    Real& operator()(int, int);                  // access element
  868.    Real& element(int, int);                     // access element
  869. #ifdef SETUP_C_SUBSCRIPTS
  870.    Real* operator[](int m) { return store+lower*(m+1); }
  871.    const Real* operator[](int m) const { return store+lower*(m+1); }
  872. #endif
  873.    MatrixType Type() const;
  874.    SymmetricBandMatrix(const SymmetricBandMatrix& gm) { GetMatrix(&gm); }
  875. #ifndef __ZTC__
  876.    Real operator()(int, int) const;             // access element
  877. #endif
  878.    GeneralMatrix* MakeSolver();
  879.    Real SumSquare() const;
  880.    Real SumAbsoluteValue() const;
  881.    Real MaximumAbsoluteValue() const
  882.       { CornerClear(); return GeneralMatrix::MaximumAbsoluteValue(); }
  883.    Real Trace() const;
  884.    LogAndSign LogDeterminant() const;
  885.    void GetRow(MatrixRowCol&);
  886.    void GetCol(MatrixRowCol&);
  887.    GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  888.    void ReDimension(int,int);                       // change dimensions
  889.    MatrixBandWidth BandWidth() const;
  890.    void SetParameters(const GeneralMatrix*);
  891.    NEW_DELETE(SymmetricBandMatrix)
  892. };
  893.  
  894. class BandLUMatrix : public GeneralMatrix
  895. // for LU decomposition of band matrix
  896. {
  897.    int* indx;
  898.    Boolean d;
  899.    Boolean sing;                                   // TRUE if singular
  900.    Real* store2;
  901.    int storage2;
  902.    void ludcmp();
  903.    int m1,m2;                                   // lower and upper
  904. public:
  905.    BandLUMatrix(const BaseMatrix&);
  906.    MatrixType Type() const;
  907.    void lubksb(Real*, int=0);
  908.    ~BandLUMatrix();
  909.    GeneralMatrix* MakeSolver() { return this; } // for solving
  910.    LogAndSign LogDeterminant() const;
  911.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  912.    void GetRow(MatrixRowCol&);
  913.    void GetCol(MatrixRowCol&);
  914.    void operator=(const BaseMatrix&);
  915.    void operator=(const BandLUMatrix& m) { operator=((const BaseMatrix&)m); }
  916.    NEW_DELETE(BandLUMatrix)
  917.    void CleanUp();                                // to clear store
  918. };
  919.  
  920.  
  921. /***************************** temporary classes *************************/
  922.  
  923. class MultipliedMatrix : public BaseMatrix
  924. {
  925. protected:
  926.    union { const BaseMatrix* bm1; GeneralMatrix* gm1; };
  927.                           // pointers to summands
  928.    union { const BaseMatrix* bm2; GeneralMatrix* gm2; };
  929.    MultipliedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
  930.       : bm1(bm1x),bm2(bm2x) {}
  931.    int search(const BaseMatrix*) const;
  932.    friend class BaseMatrix;
  933. public:
  934.    GeneralMatrix* Evaluate(MatrixType);
  935.    MatrixBandWidth BandWidth() const;
  936.    NEW_DELETE(MultipliedMatrix)
  937. };
  938.  
  939. class AddedMatrix : public MultipliedMatrix
  940. {
  941. protected:
  942.    AddedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
  943.       : MultipliedMatrix(bm1x,bm2x) {}
  944.  
  945.    friend class BaseMatrix;
  946. public:
  947.    GeneralMatrix* Evaluate(MatrixType);
  948.    MatrixBandWidth BandWidth() const;
  949. #ifdef __GNUG__
  950.    void SelectVersion(MatrixType, int&, int&) const;
  951. #else
  952.    void SelectVersion(MatrixType, Boolean&, Boolean&) const;
  953. #endif
  954.    NEW_DELETE(AddedMatrix)
  955. };
  956.  
  957. class SolvedMatrix : public MultipliedMatrix
  958. {
  959.    SolvedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
  960.       : MultipliedMatrix(bm1x,bm2x) {}
  961.    friend class BaseMatrix;
  962.    friend class InvertedMatrix;                        // for operator*
  963. public:
  964.    GeneralMatrix* Evaluate(MatrixType);
  965.    MatrixBandWidth BandWidth() const;
  966.    NEW_DELETE(SolvedMatrix)
  967. };
  968.  
  969. class SubtractedMatrix : public AddedMatrix
  970. {
  971.    SubtractedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
  972.       : AddedMatrix(bm1x,bm2x) {}
  973.    friend class BaseMatrix;
  974. public:
  975.    GeneralMatrix* Evaluate(MatrixType);
  976.    NEW_DELETE(SubtractedMatrix)
  977. };
  978.  
  979. class ShiftedMatrix : public BaseMatrix
  980. {
  981. protected:
  982.    Real f;
  983.    union { const BaseMatrix* bm; GeneralMatrix* gm; };
  984.    ShiftedMatrix(const BaseMatrix* bmx, Real fx) : bm(bmx),f(fx) {}
  985.    int search(const BaseMatrix*) const;
  986.    friend class BaseMatrix;
  987. public:
  988.    GeneralMatrix* Evaluate(MatrixType);
  989.    NEW_DELETE(ShiftedMatrix)
  990. };
  991.  
  992. class ScaledMatrix : public ShiftedMatrix
  993. {
  994.    ScaledMatrix(const BaseMatrix* bmx, Real fx) : ShiftedMatrix(bmx,fx) {}
  995.    friend class BaseMatrix;
  996. public:
  997.    GeneralMatrix* Evaluate(MatrixType);
  998.    MatrixBandWidth BandWidth() const;
  999.    NEW_DELETE(ScaledMatrix)
  1000. };
  1001.  
  1002. class NegatedMatrix : public BaseMatrix
  1003. {
  1004. protected:
  1005.    union { const BaseMatrix* bm; GeneralMatrix* gm; };
  1006.    NegatedMatrix(const BaseMatrix* bmx) : bm(bmx) {}
  1007.    int search(const BaseMatrix*) const;
  1008. private:
  1009.    friend class BaseMatrix;
  1010. public:
  1011.    GeneralMatrix* Evaluate(MatrixType);
  1012.    MatrixBandWidth BandWidth() const;
  1013.    NEW_DELETE(NegatedMatrix)
  1014. };
  1015.  
  1016. class TransposedMatrix : public NegatedMatrix
  1017. {
  1018.    TransposedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  1019.    friend class BaseMatrix;
  1020. public:
  1021.    GeneralMatrix* Evaluate(MatrixType);
  1022.    MatrixBandWidth BandWidth() const;
  1023.    NEW_DELETE(TransposedMatrix)
  1024. };
  1025.  
  1026. class InvertedMatrix : public NegatedMatrix
  1027. {
  1028.    InvertedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  1029. public:
  1030. #ifndef TEMPS_DESTROYED_QUICKLY
  1031.    SolvedMatrix operator*(const BaseMatrix&) const;  // inverse(A) * B
  1032. #else
  1033.    SolvedMatrix& operator*(const BaseMatrix&);       // inverse(A) * B
  1034. #endif
  1035.    friend class BaseMatrix;
  1036.    GeneralMatrix* Evaluate(MatrixType);
  1037.    MatrixBandWidth BandWidth() const;
  1038.    NEW_DELETE(InvertedMatrix)
  1039. };
  1040.  
  1041. class RowedMatrix : public NegatedMatrix
  1042. {
  1043.    RowedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  1044.    friend class BaseMatrix;
  1045. public:
  1046.    GeneralMatrix* Evaluate(MatrixType);
  1047.    MatrixBandWidth BandWidth() const;
  1048.    NEW_DELETE(RowedMatrix)
  1049. };
  1050.  
  1051. class ColedMatrix : public NegatedMatrix
  1052. {
  1053.    ColedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  1054.    friend class BaseMatrix;
  1055. public:
  1056.    GeneralMatrix* Evaluate(MatrixType);
  1057.    MatrixBandWidth BandWidth() const;
  1058.    NEW_DELETE(ColedMatrix)
  1059. };
  1060.  
  1061. class DiagedMatrix : public NegatedMatrix
  1062. {
  1063.    DiagedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  1064.    friend class BaseMatrix;
  1065. public:
  1066.    GeneralMatrix* Evaluate(MatrixType);
  1067.    MatrixBandWidth BandWidth() const;
  1068.    NEW_DELETE(DiagedMatrix)
  1069. };
  1070.  
  1071. class MatedMatrix : public NegatedMatrix
  1072. {
  1073.    int nr, nc;
  1074.    MatedMatrix(const BaseMatrix* bmx, int nrx, int ncx)
  1075.       : NegatedMatrix(bmx), nr(nrx), nc(ncx) {}
  1076.    friend class BaseMatrix;
  1077. public:
  1078.    GeneralMatrix* Evaluate(MatrixType);
  1079.    MatrixBandWidth BandWidth() const;
  1080.    NEW_DELETE(MatedMatrix)
  1081. };
  1082.  
  1083. class ConstMatrix : public BaseMatrix
  1084. {
  1085.    const GeneralMatrix* cgm;
  1086.    int search(const BaseMatrix*) const;
  1087.    ConstMatrix(const GeneralMatrix* cgmx) : cgm(cgmx) {}
  1088.    friend class BaseMatrix;
  1089.    friend class GeneralMatrix;
  1090. public:
  1091.    GeneralMatrix* Evaluate(MatrixType);
  1092.    MatrixBandWidth BandWidth() const;
  1093.    NEW_DELETE(ConstMatrix)
  1094. };
  1095.  
  1096. class ReturnMatrixX : public BaseMatrix    // for matrix return
  1097. {
  1098.    GeneralMatrix* gm;
  1099.    int search(const BaseMatrix*) const;
  1100. public:
  1101.    GeneralMatrix* Evaluate(MatrixType);
  1102.    friend class BaseMatrix;
  1103. #ifdef TEMPS_DESTROYED_QUICKLY
  1104.    ReturnMatrixX(const ReturnMatrixX& tm);
  1105. #else
  1106.    ReturnMatrixX(const ReturnMatrixX& tm) : gm(tm.gm) {}
  1107. #endif
  1108.    ReturnMatrixX(const GeneralMatrix* gmx) : gm((GeneralMatrix*&)gmx) {}
  1109. //   ReturnMatrixX(GeneralMatrix&);
  1110.    MatrixBandWidth BandWidth() const;
  1111.    NEW_DELETE(ReturnMatrixX)
  1112. };
  1113.  
  1114.  
  1115. /**************************** submatrices ******************************/
  1116.  
  1117. class GetSubMatrix : public NegatedMatrix
  1118. {
  1119.    int row_skip;
  1120.    int row_number;
  1121.    int col_skip;
  1122.    int col_number;
  1123.    Boolean IsSym;
  1124.  
  1125.    GetSubMatrix
  1126.       (const BaseMatrix* bmx, int rs, int rn, int cs, int cn, Boolean is)
  1127.       : NegatedMatrix(bmx),
  1128.       row_skip(rs), row_number(rn), col_skip(cs), col_number(cn), IsSym(is) {}
  1129.    GetSubMatrix(const GetSubMatrix& g)
  1130.       : NegatedMatrix(g.bm), row_skip(g.row_skip), row_number(g.row_number),
  1131.       col_skip(g.col_skip), col_number(g.col_number), IsSym(g.IsSym) {}
  1132.    void SetUpLHS();
  1133.    friend class BaseMatrix;
  1134. public:
  1135.    GeneralMatrix* Evaluate(MatrixType);
  1136.    void operator=(const BaseMatrix&);
  1137.    void operator=(const GetSubMatrix& m) { operator=((const BaseMatrix&)m); }
  1138.    void operator<<(const BaseMatrix&);
  1139.    void operator<<(const Real*);                // copy from array
  1140.    void operator=(Real);                        // copy from constant
  1141.    void Inject(const GeneralMatrix&);           // copy stored els only
  1142.    MatrixBandWidth BandWidth() const;
  1143.    NEW_DELETE(GetSubMatrix)
  1144. };
  1145.  
  1146. /**************************** matrix input *******************************/
  1147.  
  1148. class MatrixInput          // for reading a list of values into a matrix
  1149.                            // the difficult part is detecting a mismatch
  1150.                            // in the number of elements
  1151. {
  1152.    static int n;           // number values still to be read
  1153.    static Real* r;         // pointer to last thing read
  1154.    static int depth;       // number of objects of this class in existence
  1155. public:
  1156.    MatrixInput() { depth++; }
  1157.    MatrixInput(const MatrixInput&) { depth++; }
  1158.    ~MatrixInput();
  1159.    MatrixInput operator<<(Real);
  1160.                            // could return a reference if we don't have
  1161.                            // TEMPS_DESTROYED_QUICKLY
  1162.    friend GeneralMatrix;                             
  1163. };
  1164.  
  1165. /***************************** exceptions ********************************/
  1166.  
  1167. class MatrixDetails
  1168. {
  1169.    MatrixType type;
  1170.    int nrows;
  1171.    int ncols;
  1172.    int ubw;
  1173.    int lbw;
  1174. public:
  1175.    MatrixDetails(const GeneralMatrix& A);
  1176.    void PrintOut();
  1177. };
  1178.    
  1179.  
  1180.  
  1181. class SpaceException : public Exception
  1182. {
  1183. public:
  1184.    static long st_type() { return 2; }
  1185.    long type() const { return 2; }
  1186.    static int action;
  1187.    SpaceException();
  1188.    static void SetAction(int a) { action=a; }
  1189. };
  1190.  
  1191.  
  1192. class MatrixException : public Exception
  1193. {
  1194. public:
  1195.    static long st_type() { return 3; }
  1196.    long type() const { return 3; }
  1197.    MatrixException(int);
  1198.    MatrixException(int, const GeneralMatrix&);
  1199.    MatrixException(int, const GeneralMatrix&, const GeneralMatrix&);
  1200. };
  1201.  
  1202. class DataException : public MatrixException
  1203. {
  1204. public:
  1205.    static long st_type() { return 3*53; }
  1206.    long type() const { return 3*53; }
  1207.    static int action;
  1208.    DataException(const GeneralMatrix& A);
  1209.    static void SetAction(int a) { action=a; }
  1210. };
  1211.  
  1212. class SingularException : public DataException
  1213. {
  1214. public:
  1215.    static long st_type() { return 3*53*109; }
  1216.    long type() const { return 3*53*109; }
  1217.    SingularException(const GeneralMatrix& A);
  1218. };
  1219.  
  1220. class NPDException : public DataException     // Not positive definite
  1221. {
  1222. public:
  1223.    static long st_type() { return 3*53*113; }
  1224.    long type() const { return 3*53*113; }
  1225.    NPDException(const GeneralMatrix&);
  1226. };
  1227.  
  1228. class ConvergenceException : public MatrixException
  1229. {
  1230. public:
  1231.    static long st_type() { return 3*59; }
  1232.    long type() const { return 3*59; }
  1233.    static int action;
  1234.    ConvergenceException(const GeneralMatrix& A);
  1235.    static void SetAction(int a) { action=a; }
  1236. };
  1237.  
  1238. class ProgramException : public MatrixException
  1239. {
  1240. public:
  1241.    static long st_type() { return 3*61; }
  1242.    long type() const { return 3*61; }
  1243.    static int action;
  1244.    ProgramException(char* c);
  1245.    ProgramException(char* c, const GeneralMatrix&);
  1246.    ProgramException(char* c, const GeneralMatrix&, const GeneralMatrix&);
  1247.    ProgramException();
  1248.    ProgramException(const GeneralMatrix&);
  1249.    static void SetAction(int a) { action=a; }
  1250. };
  1251.  
  1252. class IndexException : public ProgramException
  1253. {
  1254. public:
  1255.    static long st_type() { return 3*61*101; }
  1256.    long type() const { return 3*61*101; }
  1257.    IndexException(int i, const GeneralMatrix& A);
  1258.    IndexException(int i, int j, const GeneralMatrix& A);
  1259.    // next two are for access via element function
  1260.    IndexException(int i, const GeneralMatrix& A, Boolean);
  1261.    IndexException(int i, int j, const GeneralMatrix& A, Boolean);
  1262. };
  1263.  
  1264. class VectorException : public ProgramException    // can't convert to vector
  1265. {
  1266. public:
  1267.    static long st_type() { return 3*61*107; }
  1268.    long type() const { return 3*61*107; }
  1269.    VectorException();
  1270.    VectorException(const GeneralMatrix& A);
  1271. };
  1272.  
  1273. class NotSquareException : public ProgramException
  1274. {
  1275. public:
  1276.    static long st_type() { return 3*61*109; }
  1277.    long type() const { return 3*61*109; }
  1278.    NotSquareException(const GeneralMatrix& A);
  1279. };
  1280.  
  1281. class SubMatrixDimensionException : public ProgramException
  1282. {
  1283. public:
  1284.    static long st_type() { return 3*61*113; }
  1285.    long type() const { return 3*61*113; }
  1286.    SubMatrixDimensionException();
  1287. };
  1288.  
  1289. class IncompatibleDimensionsException : public ProgramException
  1290. {
  1291. public:
  1292.    static long st_type() { return 3*61*127; }
  1293.    long type() const { return 3*61*127; }
  1294.    IncompatibleDimensionsException();
  1295. };
  1296.  
  1297. class NotDefinedException : public ProgramException
  1298. {
  1299. public:
  1300.    static long st_type() { return 3*61*131; }
  1301.    long type() const { return 3*61*131; }
  1302.    NotDefinedException(char* op, char* matrix);
  1303. };
  1304.  
  1305. class CannotBuildException : public ProgramException
  1306. {
  1307. public:
  1308.    static long st_type() { return 3*61*137; }
  1309.    long type() const { return 3*61*137; }
  1310.    CannotBuildException(char* matrix);
  1311. };
  1312.  
  1313.  
  1314. class InternalException : public MatrixException
  1315. {
  1316. public:
  1317.    static long st_type() { return 3*67; }
  1318.    long type() const { return 3*67; }
  1319.    static int action;
  1320.    InternalException(char* c);
  1321.    static void SetAction(int a) { action=a; }
  1322. };
  1323.  
  1324.  
  1325. /***************************** functions ***********************************/
  1326.  
  1327.  
  1328. inline LogAndSign LogDeterminant(const BaseMatrix& B)
  1329.    { return B.LogDeterminant(); }
  1330. inline Real SumSquare(const BaseMatrix& B) { return B.SumSquare(); }
  1331. inline Real Trace(const BaseMatrix& B) { return B.Trace(); }
  1332. inline Real SumAbsoluteValue(const BaseMatrix& B)
  1333.    { return B.SumAbsoluteValue(); }
  1334. inline Real MaximumAbsoluteValue(const BaseMatrix& B)
  1335.    { return B.MaximumAbsoluteValue(); }
  1336. inline Real Norm1(const BaseMatrix& B) { return B.Norm1(); }
  1337. inline Real Norm1(RowVector& RV) { return RV.MaximumAbsoluteValue(); }
  1338. inline Real NormInfinity(const BaseMatrix& B) { return B.NormInfinity(); }
  1339. inline Real NormInfinity(ColumnVector& CV)
  1340.    { return CV.MaximumAbsoluteValue(); } 
  1341.  
  1342. #endif
  1343.