home *** CD-ROM | disk | FTP | other *** search
/ Boston 2 / boston-2.iso / DOS / PROGRAM / C / NEWMAT / NEWMAT.HXX < prev    next >
Text File  |  1993-12-01  |  27KB  |  777 lines
  1. //$$ newmat.hxx         definition file for new version of matrix package
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  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. #include "boolean.hxx"
  16.  
  17.  
  18. /**************************** general utilities ****************************/
  19.  
  20. void MatrixError(char*);                        // error handler
  21. void MatrixErrorNoSpace(void*);                 // no space handler
  22.  
  23. class LogAndSign
  24. // Return from LogDeterminant function
  25. //    - value of the log plus the sign (+, - or 0)
  26. {
  27.    real log_value;
  28.    int sign;
  29. public:
  30.    LogAndSign() { log_value=0.0; sign=1; }
  31.    void operator*=(real);
  32.    void ChangeSign() { sign = -sign; }
  33.    real LogValue() { return log_value; }
  34.    int Sign() { return sign; }
  35.    real Value();
  36. };
  37.  
  38. // the following class is for counting the number of times a piece of code
  39. // is executed. It is used for locating any code not executed by test
  40. // routines. Use turbo GREP locate all places this code is called and
  41. // check which ones are not accessed.
  42. // Somewhat implementation dependent as it relies on "cout" still being
  43. // present when ExeCounter objects are destructed.
  44.  
  45. class ExeCounter
  46. {
  47.    int line;                                    // code line number
  48.    int fileid;                                  // file identifier
  49.    long nexe;                                   // number of executions
  50.    static int nreports;                         // number of reports
  51. public:
  52.    ExeCounter(int,int);
  53.    void operator++() { nexe++; }
  54.    ~ExeCounter();                               // prints out reports
  55. };
  56.  
  57.  
  58. /**************************** class MatrixType *****************************/
  59.  
  60. // Is used for finding the type of a matrix resulting from the binary operations
  61. // +, -, * and identifying what conversions are permissible.
  62. // This class must be updated when new matrix types are added.
  63.  
  64. class GeneralMatrix;                            // defined later
  65.  
  66. class MatrixType
  67. {
  68. public:
  69.    enum Type { UnSp,UT,LT,Rect,Sym,Diag,RowV,ColV,EqEl,Crout };
  70.    static nTypes() { return 8; }               // number of different types
  71.                            // exclude Crout, UnSp
  72. private:
  73.    Type type;
  74. public:
  75.    MatrixType operator+(const MatrixType&) const;
  76.    MatrixType operator-(const MatrixType&) const;
  77.    MatrixType operator*(const MatrixType&) const;
  78.    BOOL operator>=(const MatrixType&) const;
  79.    BOOL operator==(const MatrixType& t) const; // { return (type == t.type); }
  80.    BOOL operator!=(const MatrixType& t) const; // { return !(*this == t); }
  81.    BOOL operator!() const { return type == UnSp; }
  82.    MatrixType operator-() const;               // type of negative
  83.    MatrixType i() const;                       // type of inverse
  84.    MatrixType t() const;                       // type of transpose
  85.    MatrixType sub() const;                     // type of submatrix
  86.    MatrixType ssub() const;                    // type of sym submatrix
  87.    MatrixType (Type tx) : type(tx) {}          // (& doesn't work with AT&T)
  88.    MatrixType () {}
  89.    GeneralMatrix* New() const;                 // new matrix of given type
  90.    GeneralMatrix* New(int,int) const;          // new matrix of given type
  91.    operator int() const { return (int)type; }
  92.    operator char*() const;                     // for printing type
  93. };
  94.  
  95. void TestTypeAdd();                            // test +
  96. void TestTypeMult();                           // test *
  97. void TestTypeOrder();                          // test >=
  98.  
  99. /*************************** Matrix routines ***************************/
  100.  
  101.  
  102. class MatrixRowCol;                             // defined later
  103. class MatrixRow;
  104. class MatrixCol;
  105.  
  106. class GeneralMatrix;                            // defined later
  107. class AddedMatrix;
  108. class MultipliedMatrix;
  109. class SubtractedMatrix;
  110. class SolvedMatrix;
  111. class ShiftedMatrix;
  112. class ScaledMatrix;
  113. class TransposedMatrix;
  114. class NegatedMatrix;
  115. class InvertedMatrix;
  116. class RowedMatrix;
  117. class ColedMatrix;
  118. class DiagedMatrix;
  119. class MatedMatrix;
  120. class GetSubMatrix;
  121. class ConstMatrix;
  122. class ReturnMatrix;
  123. class Matrix;
  124. class nricMatrix;
  125. class RowVector;
  126. class ColumnVector;
  127. class SymmetricMatrix;
  128. class UpperTriangularMatrix;
  129. class LowerTriangularMatrix;
  130. class DiagonalMatrix;
  131. class CroutMatrix;
  132.  
  133. static MatrixType MatrixTypeUnSp(MatrixType::UnSp);
  134.                         // AT&T needs this
  135.  
  136. class BaseMatrix                                // base of all matrix classes
  137. {
  138. protected:
  139. //   BaseMatrix() {}
  140.    virtual int search(const GeneralMatrix*) const = 0;
  141.                         // count number of times matrix
  142.                         // is referred to
  143.    virtual MatrixType Type() const = 0;         // type of a matrix
  144.    virtual int NrowsV()  const = 0;
  145.    virtual int NcolsV()  const = 0;
  146. public:
  147.    virtual GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp) = 0;
  148.                         // evaluate temporary
  149.    void MatrixParameters(int& nr, int& nc, MatrixType& mt)
  150.       { nr = NrowsV(); nc = NcolsV(); mt = Type(); }
  151.    AddedMatrix operator+(BaseMatrix&);          // results of operations
  152.    MultipliedMatrix operator*(BaseMatrix&);
  153.    SubtractedMatrix operator-(BaseMatrix&);
  154.    ShiftedMatrix operator+(real);
  155.    ScaledMatrix operator*(real);
  156.    ScaledMatrix operator/(real);
  157.    ShiftedMatrix operator-(real);
  158.    TransposedMatrix t();
  159.    NegatedMatrix operator-();                   // change sign of elements
  160.    InvertedMatrix i();
  161.    RowedMatrix CopyToRow();
  162.    ColedMatrix CopyToColumn();
  163.    DiagedMatrix CopyToDiagonal();
  164.    MatedMatrix CopyToMatrix(int,int);
  165.    GetSubMatrix SubMatrix(int,int,int,int);
  166.    GetSubMatrix SymSubMatrix(int,int);
  167.    GetSubMatrix Row(int);
  168.    GetSubMatrix Rows(int,int);
  169.    GetSubMatrix Column(int);
  170.    GetSubMatrix Columns(int,int);
  171.    operator real();                            // conversion of 1 x 1 matrix
  172.    virtual LogAndSign LogDeterminant();
  173.    virtual real SumSquare();
  174.    virtual real SumAbsoluteValue();
  175.    virtual real MaximumAbsoluteValue();
  176.    virtual real Trace();
  177.    real Norm1();
  178.    real NormInfinity();
  179.  
  180.    friend GeneralMatrix;
  181.    friend Matrix;
  182.    friend nricMatrix;
  183.    friend RowVector;
  184.    friend ColumnVector;
  185.    friend SymmetricMatrix;
  186.    friend UpperTriangularMatrix;
  187.    friend LowerTriangularMatrix;
  188.    friend DiagonalMatrix;
  189.    friend CroutMatrix;
  190.    friend AddedMatrix;
  191.    friend MultipliedMatrix;
  192.    friend SubtractedMatrix;
  193.    friend SolvedMatrix;
  194.    friend ShiftedMatrix;
  195.    friend ScaledMatrix;
  196.    friend TransposedMatrix;
  197.    friend NegatedMatrix;
  198.    friend InvertedMatrix;
  199.    friend RowedMatrix;
  200.    friend ColedMatrix;
  201.    friend DiagedMatrix;
  202.    friend MatedMatrix;
  203.    friend GetSubMatrix;
  204.    friend ConstMatrix;
  205.    friend ReturnMatrix;
  206. };
  207.  
  208.  
  209. /******************************* working classes **************************/
  210.  
  211. class GeneralMatrix : public BaseMatrix         // declarable matrix types
  212. {
  213. protected:
  214.    int tag;                                     // shows whether can reuse
  215.    int nrows, ncols;                            // dimensions
  216.    int storage;                                 // total store required
  217.    real* store;                                 // point to store (0=not set)
  218.    GeneralMatrix();                             // initialise with no store
  219.    GeneralMatrix(int);                          // constructor getting store
  220.    void Add(GeneralMatrix*, real);              // sum of GM and real
  221.    void Add(real);                              // add real to this
  222.    void Multiply(GeneralMatrix*, real);         // product of GM and real
  223.    void Multiply(real);                         // multiply this by real
  224.    void Negate(GeneralMatrix*);                 // change sign
  225.    void Negate();                               // change sign
  226.    void operator=(real);                        // set matrix to constant
  227.    real* GetStore();                            // get store or copy
  228.    GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
  229.                                                 // temporarily access store
  230.    void GetMatrix(GeneralMatrix*);              // used by = and initialise
  231. #ifndef __ZTC__
  232.    void GetMatrixC(const GeneralMatrix*);       // used by = and initialise
  233. #endif
  234.    void Eq(BaseMatrix&, MatrixType);            // used by =
  235.    int search(const GeneralMatrix*) const;
  236.    virtual GeneralMatrix* Transpose(MatrixType);
  237.    void CheckConversion(const BaseMatrix&);     // check conversion OK
  238.    void ReDimension(int, int, int);             // change dimensions
  239.    int NrowsV() const;                          // get dimensions
  240.    int NcolsV() const;                          // virtual version
  241. public:
  242.    GeneralMatrix* Evaluate(MatrixType);
  243.    MatrixType Type() const = 0;                 // type of a matrix
  244.    int Nrows() const { return nrows; }          // get dimensions
  245.    int Ncols() const { return ncols; }
  246.    int Storage() const { return storage; }
  247.    real* Store() const { return store; }
  248.    virtual ~GeneralMatrix();                    // delete store if set
  249.    void tDelete();                              // delete if tag permits
  250.    BOOL reuse();                                // TRUE if tag allows reuse
  251.    void Protect() { tag=-1; }                   // can't delete or reuse
  252.    int Tag() const { return tag; }
  253.    BOOL IsZero() const;                         // test matrix has all zeros
  254.    void Release() { tag=1; }                    // del store after next use
  255.    void Release(int t) { tag=t; }               // del store after t accesses
  256.    void ReleaseAndDelete() { tag=0; }           // delete matrix after use
  257.    void operator<<(const real*);                // assignment from an array
  258.    void operator<<(BaseMatrix& X) { Eq(X,this->Type()); }
  259.                                                 // = without checking type
  260.    void Inject(const GeneralMatrix&);           // copy stored els only
  261.    virtual GeneralMatrix* MakeSolver();         // for solving
  262.    virtual void Solver(MatrixRowCol&, const MatrixRowCol&) {}
  263.    virtual void GetRow(MatrixRowCol&) = 0;      // Get matrix row
  264.    virtual void RestoreRow(MatrixRowCol&) {}    // Restore matrix row
  265.    virtual void NextRow(MatrixRowCol&);         // Go to next row
  266.    virtual void GetCol(MatrixRowCol&) = 0;      // Get matrix col
  267.    virtual void RestoreCol(MatrixRowCol&) {}    // Restore matrix col
  268.    virtual void NextCol(MatrixRowCol&);         // Go to next col
  269.    real SumSquare();
  270.    real SumAbsoluteValue();
  271.    real MaximumAbsoluteValue();
  272.    LogAndSign LogDeterminant();
  273.    ConstMatrix c() const;                       // to access constant matrices
  274.    void CheckStore() const;                     // check store is non-zero
  275.  
  276.    friend Matrix;
  277.    friend nricMatrix;
  278.    friend SymmetricMatrix;
  279.    friend UpperTriangularMatrix;
  280.    friend LowerTriangularMatrix;
  281.    friend DiagonalMatrix;
  282.    friend CroutMatrix;
  283.    friend RowVector;
  284.    friend ColumnVector;
  285.    friend BaseMatrix;
  286.    friend AddedMatrix;
  287.    friend MultipliedMatrix;
  288.    friend SubtractedMatrix;
  289.    friend SolvedMatrix;
  290.    friend ShiftedMatrix;
  291.    friend ScaledMatrix;
  292.    friend TransposedMatrix;
  293.    friend NegatedMatrix;
  294.    friend InvertedMatrix;
  295.    friend RowedMatrix;
  296.    friend ColedMatrix;
  297.    friend DiagedMatrix;
  298.    friend MatedMatrix;
  299.    friend GetSubMatrix;
  300.    friend ConstMatrix;
  301.    friend ReturnMatrix;
  302. };
  303.  
  304. class Matrix : public GeneralMatrix             // usual rectangular matrix
  305. {
  306. public:
  307.    Matrix() {}
  308.    Matrix(int, int);                            // standard declaration
  309.    Matrix(BaseMatrix&);                         // evaluate BaseMatrix
  310.    void operator=(BaseMatrix&);
  311.    void operator=(real f) { GeneralMatrix::operator=(f); }
  312.    MatrixType Type() const;
  313.    real& operator()(int, int);                  // access element
  314.    real& element(int, int);                     // access element
  315.    Matrix(Matrix& gm) { GetMatrix(&gm); }
  316. #ifndef __ZTC__
  317.    real operator()(int, int) const;             // access element
  318.    Matrix(const Matrix& gm) { GetMatrixC(&gm); }
  319. #endif
  320.    GeneralMatrix* MakeSolver();
  321.    real Trace();
  322.    void GetRow(MatrixRowCol&);
  323.    void GetCol(MatrixRowCol&);
  324.    void RestoreCol(MatrixRowCol&);
  325.    void NextRow(MatrixRowCol&);
  326.    void NextCol(MatrixRowCol&);
  327.    void ReDimension(int,int);                   // change dimensions
  328. };
  329.  
  330. class nricMatrix : public Matrix                // for use with Numerical
  331.                                                 // Recipes in C
  332. {
  333.    real** row_pointer;                          // points to rows
  334.    void MakeRowPointer();                       // build rowpointer
  335.    void DeleteRowPointer();
  336. public:
  337.    nricMatrix() {}
  338.    nricMatrix(int m, int n)                     // standard declaration
  339.       :  Matrix(m,n) { MakeRowPointer(); }
  340.    nricMatrix(BaseMatrix& bm)                   // evaluate BaseMatrix
  341.       :  Matrix(bm) { MakeRowPointer(); }
  342.    void operator=(BaseMatrix& bm)
  343.       { DeleteRowPointer(); Matrix::operator=(bm); MakeRowPointer(); }
  344.    void operator=(real f) { GeneralMatrix::operator=(f); }
  345.    void operator<<(BaseMatrix& X)
  346.       { DeleteRowPointer(); Eq(X,this->Type()); MakeRowPointer(); }
  347.    nricMatrix(nricMatrix& gm) { GetMatrix(&gm); MakeRowPointer(); }
  348. #ifndef __ZTC__
  349.    nricMatrix(const nricMatrix& gm) { GetMatrixC(&gm); MakeRowPointer(); }
  350. #endif
  351.    void ReDimension(int m, int n)               // change dimensions
  352.       { DeleteRowPointer(); Matrix::ReDimension(m,n); MakeRowPointer(); }
  353.    ~nricMatrix() { DeleteRowPointer(); }
  354. #ifndef __ZTC__
  355.    operator real**() const { CheckStore(); return row_pointer-1; }
  356. #endif
  357. };
  358.  
  359. class SymmetricMatrix : public GeneralMatrix
  360. {
  361. public:
  362.    SymmetricMatrix() {}
  363.    SymmetricMatrix(int);
  364.    SymmetricMatrix(BaseMatrix&);
  365.    void operator=(BaseMatrix&);
  366.    void operator=(real f) { GeneralMatrix::operator=(f); }
  367.    real& operator()(int, int);                  // access element
  368.    real& element(int, int);                     // access element
  369.    MatrixType Type() const;
  370.    SymmetricMatrix(SymmetricMatrix& gm) { GetMatrix(&gm); }
  371. #ifndef __ZTC__
  372.    real operator()(int, int) const;             // access element
  373.    SymmetricMatrix(const SymmetricMatrix& gm) { GetMatrixC(&gm); }
  374. #endif
  375.    real SumSquare();
  376.    real SumAbsoluteValue();
  377.    real Trace();
  378.    void GetRow(MatrixRowCol&);
  379.    void GetCol(MatrixRowCol&);
  380.    GeneralMatrix* Transpose(MatrixType);
  381.    void ReDimension(int);                       // change dimensions
  382. };
  383.  
  384. class UpperTriangularMatrix : public GeneralMatrix
  385. {
  386. public:
  387.    UpperTriangularMatrix() {}
  388.    UpperTriangularMatrix(int);
  389.    void operator=(BaseMatrix&);
  390.    UpperTriangularMatrix(BaseMatrix&);
  391.    UpperTriangularMatrix(UpperTriangularMatrix& gm) { GetMatrix(&gm); }
  392. #ifndef __ZTC__
  393.    real operator()(int, int) const;             // access element
  394.    UpperTriangularMatrix(const UpperTriangularMatrix& gm) { GetMatrixC(&gm); }
  395. #endif
  396.    void operator=(real f) { GeneralMatrix::operator=(f); }
  397.    real& operator()(int, int);                  // access element
  398.    real& element(int, int);                     // access element
  399.    MatrixType Type() const;
  400.    GeneralMatrix* MakeSolver() { return this; } // for solving
  401.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  402.    LogAndSign LogDeterminant();
  403.    real Trace();
  404.    void GetRow(MatrixRowCol&);
  405.    void GetCol(MatrixRowCol&);
  406.    void RestoreCol(MatrixRowCol&);
  407.    void NextRow(MatrixRowCol&);
  408.    void ReDimension(int);                       // change dimensions
  409. };
  410.  
  411. class LowerTriangularMatrix : public GeneralMatrix
  412. {
  413. public:
  414.    LowerTriangularMatrix() {}
  415.    LowerTriangularMatrix(int);
  416.    LowerTriangularMatrix(LowerTriangularMatrix& gm) { GetMatrix(&gm); }
  417. #ifndef __ZTC__
  418.    real operator()(int, int) const;             // access element
  419.    LowerTriangularMatrix(const LowerTriangularMatrix& gm) { GetMatrixC(&gm); }
  420. #endif
  421.    LowerTriangularMatrix(BaseMatrix& M);
  422.    void operator=(BaseMatrix&);
  423.    void operator=(real f) { GeneralMatrix::operator=(f); }
  424.    real& operator()(int, int);                  // access element
  425.    real& element(int, int);                     // access element
  426.    MatrixType Type() const;
  427.    GeneralMatrix* MakeSolver() { return this; } // for solving
  428.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  429.    LogAndSign LogDeterminant();
  430.    real Trace();
  431.    void GetRow(MatrixRowCol&);
  432.    void GetCol(MatrixRowCol&);
  433.    void RestoreCol(MatrixRowCol&);
  434.    void NextRow(MatrixRowCol&);
  435.    void ReDimension(int);                       // change dimensions
  436. };
  437.  
  438. class DiagonalMatrix : public GeneralMatrix
  439. {
  440. public:
  441.    DiagonalMatrix() {}
  442.    DiagonalMatrix(int);
  443.    DiagonalMatrix(BaseMatrix&);
  444.    DiagonalMatrix(DiagonalMatrix& gm) { GetMatrix(&gm); }
  445. #ifndef __ZTC__
  446.    real operator()(int, int) const;             // access element
  447.    real operator()(int) const;
  448.    DiagonalMatrix(const DiagonalMatrix& gm) { GetMatrixC(&gm); }
  449. #endif
  450.    void operator=(BaseMatrix&);
  451.    void operator=(real f) { GeneralMatrix::operator=(f); }
  452.    real& operator()(int, int);                  // access element
  453.    real& operator()(int);                       // access element
  454.    real& element(int, int);                     // access element
  455.    real& element(int);                          // access element
  456.    MatrixType Type() const;
  457.  
  458.    LogAndSign LogDeterminant();
  459.    real Trace();
  460.    void GetRow(MatrixRowCol&);
  461.    void GetCol(MatrixRowCol&);
  462.    void NextRow(MatrixRowCol&);
  463.    void NextCol(MatrixRowCol&);
  464.    GeneralMatrix* MakeSolver() { return this; } // for solving
  465.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  466.    GeneralMatrix* Transpose(MatrixType);
  467.    void ReDimension(int);                       // change dimensions
  468. #ifndef __ZTC__
  469.    operator real*() const
  470.       { CheckStore(); return store-1; }         // for use by NRIC
  471. #endif
  472. };
  473.  
  474. class RowVector : public Matrix
  475. {
  476. public:
  477.    RowVector() {}
  478.    RowVector(int n) : Matrix(1,n) {}
  479.    RowVector(BaseMatrix&);
  480.    RowVector(RowVector& gm) { GetMatrix(&gm); }
  481. #ifndef __ZTC__
  482.    real operator()(int) const;                  // access element
  483.    RowVector(const RowVector& gm) {GetMatrixC(&gm); }
  484. #endif
  485.    void operator=(BaseMatrix&);
  486.    void operator=(real f) { GeneralMatrix::operator=(f); }
  487.    real& operator()(int);                       // access element
  488.    real& element(int);                          // access element
  489.    MatrixType Type() const;
  490.    void GetCol(MatrixRowCol&);
  491.    void NextCol(MatrixRowCol&);
  492.    void RestoreCol(MatrixRowCol&);
  493.    GeneralMatrix* Transpose(MatrixType);
  494.    void ReDimension(int);                       // change dimensions
  495. #ifndef __ZTC__
  496.    operator real*() const
  497.       { CheckStore(); return store-1; }         // for use by NRIC
  498. #endif
  499. };
  500.  
  501. class ColumnVector : public Matrix
  502. {
  503. public:
  504.    ColumnVector() {}
  505.    ColumnVector(int n) : Matrix(n,1) {}
  506.    ColumnVector(BaseMatrix&);
  507.    ColumnVector(ColumnVector& gm) { GetMatrix(&gm); }
  508. #ifndef __ZTC__
  509.    real operator()(int) const;                  // access element
  510.    ColumnVector(const ColumnVector& gm) { GetMatrixC(&gm); }
  511. #endif
  512.    void operator=(BaseMatrix&);
  513.    void operator=(real f) { GeneralMatrix::operator=(f); }
  514.    real& operator()(int);                       // access element
  515.    real& element(int);                          // access element
  516.    MatrixType Type() const;
  517.    GeneralMatrix* Transpose(MatrixType);
  518.    void ReDimension(int);                       // change dimensions
  519. #ifndef __ZTC__
  520.    operator real*() const
  521.       { CheckStore(); return store-1; }         // for use by NRIC
  522. #endif
  523. };
  524.  
  525. class CroutMatrix : public GeneralMatrix        // for LU decomposition
  526. {
  527.    int* indx;
  528.    BOOL d;
  529.    void ludcmp();
  530. public:
  531.    CroutMatrix(BaseMatrix&);
  532.    MatrixType Type() const;
  533.    void lubksb(real*, int=0);
  534.    ~CroutMatrix() { delete [nrows] indx; }
  535.    GeneralMatrix* MakeSolver() { return this; } // for solving
  536.    LogAndSign LogDeterminant();
  537.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  538.    void GetRow(MatrixRowCol&) { MatrixError("GetRow not defined for Crout"); }
  539.    void GetCol(MatrixRowCol&) { MatrixError("GetCol not defined for Crout"); }
  540.    void operator=(BaseMatrix&)
  541.       { MatrixError("operator= not defined for Crout"); }
  542. };
  543.  
  544.  
  545. /***************************** temporary classes *************************/
  546.  
  547. class AddedMatrix : public BaseMatrix
  548. {
  549. protected:
  550.    BaseMatrix* bm1;                             // pointers to summands
  551.    BaseMatrix* bm2;
  552.    AddedMatrix(BaseMatrix* bm1x, BaseMatrix* bm2x) : bm1(bm1x),bm2(bm2x) {}
  553.    int search(const GeneralMatrix*) const;
  554.    int NrowsV() const;
  555.    int NcolsV() const;
  556. private:
  557.    MatrixType Type() const;
  558.    friend BaseMatrix;
  559. public:
  560.    GeneralMatrix* Evaluate(MatrixType);
  561. };
  562.  
  563. class MultipliedMatrix : public AddedMatrix
  564. {
  565.    MultipliedMatrix(BaseMatrix* bm1x, BaseMatrix* bm2x)
  566.       : AddedMatrix(bm1x,bm2x) {}
  567.    MatrixType Type() const;
  568.    friend BaseMatrix;
  569. public:
  570.    GeneralMatrix* Evaluate(MatrixType);
  571. };
  572.  
  573. class SolvedMatrix : public AddedMatrix
  574. {
  575.    SolvedMatrix(BaseMatrix* bm1x, BaseMatrix* bm2x)
  576.       : AddedMatrix(bm1x,bm2x) {}
  577.    MatrixType Type() const;
  578.    friend BaseMatrix;
  579.    friend InvertedMatrix;                        // for operator*
  580. public:
  581.    GeneralMatrix* Evaluate(MatrixType);
  582. };
  583.  
  584. class SubtractedMatrix : public AddedMatrix
  585. {
  586.    SubtractedMatrix(BaseMatrix* bm1x, BaseMatrix* bm2x)
  587.       : AddedMatrix(bm1x,bm2x) {}
  588.    MatrixType Type() const;
  589.    friend BaseMatrix;
  590. public:
  591.    GeneralMatrix* Evaluate(MatrixType);
  592. };
  593.  
  594. class ShiftedMatrix : public BaseMatrix
  595. {
  596. protected:
  597.    real f;
  598.    BaseMatrix* bm;
  599.    ShiftedMatrix(BaseMatrix* bmx, real fx) : bm(bmx),f(fx) {}
  600.    int search(const GeneralMatrix*) const;
  601.    int NrowsV() const;
  602.    int NcolsV() const;
  603. private:
  604.    MatrixType Type() const;
  605.    friend BaseMatrix;
  606. public:
  607.    GeneralMatrix* Evaluate(MatrixType);
  608. };
  609.  
  610. class ScaledMatrix : public ShiftedMatrix
  611. {
  612.    ScaledMatrix(BaseMatrix* bmx, real fx) : ShiftedMatrix(bmx,fx) {}
  613.    MatrixType Type() const;
  614.    friend BaseMatrix;
  615. public:
  616.    GeneralMatrix* Evaluate(MatrixType);
  617. };
  618.  
  619. class NegatedMatrix : public BaseMatrix
  620. {
  621. protected:
  622.    BaseMatrix* bm;
  623.    NegatedMatrix(BaseMatrix* bmx) : bm(bmx) {}
  624.    int search(const GeneralMatrix*) const;
  625.    int NrowsV() const;
  626.    int NcolsV() const;
  627. private:
  628.    MatrixType Type() const;
  629.    friend BaseMatrix;
  630. public:
  631.    GeneralMatrix* Evaluate(MatrixType);
  632. };
  633.  
  634. class TransposedMatrix : public NegatedMatrix
  635. {
  636.    TransposedMatrix(BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  637.    int NrowsV() const;
  638.    int NcolsV() const;
  639.    MatrixType Type() const;
  640.    friend BaseMatrix;
  641. public:
  642.    GeneralMatrix* Evaluate(MatrixType);
  643. };
  644.  
  645. class InvertedMatrix : public NegatedMatrix
  646. {
  647.    InvertedMatrix(BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  648.    MatrixType Type() const;
  649. public:
  650.    SolvedMatrix operator*(BaseMatrix&);         // inverse(A) * B
  651.    friend BaseMatrix;
  652.    GeneralMatrix* Evaluate(MatrixType);
  653. };
  654.  
  655. class RowedMatrix : public NegatedMatrix
  656. {
  657.    MatrixType Type() const;
  658.    int NrowsV() const;
  659.    int NcolsV() const;
  660.    RowedMatrix(BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  661.    friend BaseMatrix;
  662. public:
  663.    GeneralMatrix* Evaluate(MatrixType);
  664. };
  665.  
  666. class ColedMatrix : public NegatedMatrix
  667. {
  668.    MatrixType Type() const;
  669.    int NrowsV() const;
  670.    int NcolsV() const;
  671.    ColedMatrix(BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  672.    friend BaseMatrix;
  673. public:
  674.    GeneralMatrix* Evaluate(MatrixType);
  675. };
  676.  
  677. class DiagedMatrix : public NegatedMatrix
  678. {
  679.    MatrixType Type() const;
  680.    int NrowsV() const;
  681.    int NcolsV() const;
  682.    DiagedMatrix(BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  683.    friend BaseMatrix;
  684. public:
  685.    GeneralMatrix* Evaluate(MatrixType);
  686. };
  687.  
  688. class MatedMatrix : public NegatedMatrix
  689. {
  690.    int nr, nc;
  691.    MatrixType Type() const;
  692.    int NrowsV() const;
  693.    int NcolsV() const;
  694.    MatedMatrix(BaseMatrix* bmx, int nrx, int ncx)
  695.       : NegatedMatrix(bmx), nr(nrx), nc(ncx) {}
  696.    friend BaseMatrix;
  697. public:
  698.    GeneralMatrix* Evaluate(MatrixType);
  699. };
  700.  
  701. class ConstMatrix : public BaseMatrix
  702. {
  703.    const GeneralMatrix* cgm;
  704.    MatrixType Type() const;
  705.    int NrowsV() const;
  706.    int NcolsV() const;
  707.    int search(const GeneralMatrix*) const;
  708.    ConstMatrix(const GeneralMatrix* cgmx) : cgm(cgmx) {}
  709.    friend BaseMatrix;
  710.    friend GeneralMatrix;
  711. public:
  712.    GeneralMatrix* Evaluate(MatrixType);
  713. };
  714.  
  715. class ReturnMatrix : public BaseMatrix    // for matrix return
  716. {
  717.    GeneralMatrix* gm;
  718.    MatrixType Type() const;
  719.    int NrowsV() const;
  720.    int NcolsV() const;
  721.    int search(const GeneralMatrix*) const;
  722. public:
  723.    GeneralMatrix* Evaluate(MatrixType);
  724.    friend BaseMatrix;
  725.    ReturnMatrix(const ReturnMatrix& tm) : gm(tm.gm) {}
  726.    ReturnMatrix(GeneralMatrix&);
  727. };
  728.  
  729.  
  730. /**************************** submatrices ******************************/
  731.  
  732. class GetSubMatrix : public NegatedMatrix
  733. {
  734.    int row_skip;
  735.    int row_number;
  736.    int col_skip;
  737.    int col_number;
  738.    MatrixType mt;
  739.  
  740.    GetSubMatrix
  741.       (BaseMatrix* bmx, int rs, int rn, int cs, int cn, MatrixType mtx)
  742.       : NegatedMatrix(bmx),
  743.       row_skip(rs), row_number(rn), col_skip(cs), col_number(cn), mt(mtx) {}
  744.    GetSubMatrix(const GetSubMatrix& g)
  745.       : NegatedMatrix(g.bm), row_skip(g.row_skip), row_number(g.row_number),
  746.       col_skip(g.col_skip), col_number(g.col_number), mt(g.mt) {}
  747.    MatrixType Type() const;
  748.    int NrowsV() const;
  749.    int NcolsV() const;
  750.    friend BaseMatrix;
  751. public:
  752.    GeneralMatrix* Evaluate(MatrixType);
  753.    void operator=(BaseMatrix&)
  754.       { MatrixError("= not defined for submatrix; use <<"); }
  755.    void operator<<(BaseMatrix&);
  756.    void operator<<(const real*);                // copy from array
  757.    void operator<<(real);                       // copy from constant
  758.    void Inject(const GeneralMatrix&);           // copy stored els only
  759. };
  760.  
  761. /***************************** functions ***********************************/
  762.  
  763.  
  764. inline LogAndSign LogDeterminant(BaseMatrix& B) { return B.LogDeterminant(); }
  765. inline real SumSquare(BaseMatrix& B) { return B.SumSquare(); }
  766. inline real Trace(BaseMatrix& B) { return B.Trace(); }
  767. inline real SumAbsoluteValue(BaseMatrix& B) { return B.SumAbsoluteValue(); }
  768. inline real MaximumAbsoluteValue(BaseMatrix& B)
  769.    { return B.MaximumAbsoluteValue(); }
  770. inline real Norm1(BaseMatrix& B) { return B.Norm1(); }
  771. inline real Norm1(RowVector& RV) { return RV.MaximumAbsoluteValue(); }
  772. inline real NormInfinity(BaseMatrix& B) { return B.NormInfinity(); }
  773. inline real NormInfinity(ColumnVector& CV)
  774.    { return CV.MaximumAbsoluteValue(); } 
  775.  
  776. #endif
  777.