home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume26 / newmat04 / part02 / newmat.hxx < prev    next >
Encoding:
Text File  |  1991-11-30  |  25.8 KB  |  782 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. #ifndef GXX
  148.    virtual GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp) = 0;
  149.                         // evaluate temporary
  150. #else
  151.    virtual GeneralMatrix* Evaluate(MatrixType mt) = 0;
  152.    GeneralMatrix* Evaluate() { return Evaluate(MatrixTypeUnSp); }
  153. #endif
  154.    void MatrixParameters(int& nr, int& nc, MatrixType& mt)
  155.       { nr = NrowsV(); nc = NcolsV(); mt = Type(); }
  156.    AddedMatrix operator+(BaseMatrix&);          // results of operations
  157.    MultipliedMatrix operator*(BaseMatrix&);
  158.    SubtractedMatrix operator-(BaseMatrix&);
  159.    ShiftedMatrix operator+(real);
  160.    ScaledMatrix operator*(real);
  161.    ScaledMatrix operator/(real);
  162.    ShiftedMatrix operator-(real);
  163.    TransposedMatrix t();
  164.    NegatedMatrix operator-();                   // change sign of elements
  165.    InvertedMatrix i();
  166.    RowedMatrix CopyToRow();
  167.    ColedMatrix CopyToColumn();
  168.    DiagedMatrix CopyToDiagonal();
  169.    MatedMatrix CopyToMatrix(int,int);
  170.    GetSubMatrix SubMatrix(int,int,int,int);
  171.    GetSubMatrix SymSubMatrix(int,int);
  172.    GetSubMatrix Row(int);
  173.    GetSubMatrix Rows(int,int);
  174.    GetSubMatrix Column(int);
  175.    GetSubMatrix Columns(int,int);
  176.    operator real();                            // conversion of 1 x 1 matrix
  177.    virtual LogAndSign LogDeterminant();
  178.    virtual real SumSquare();
  179.    virtual real SumAbsoluteValue();
  180.    virtual real MaximumAbsoluteValue();
  181.    virtual real Trace();
  182.    real Norm1();
  183.    real NormInfinity();
  184.  
  185.    friend GeneralMatrix;
  186.    friend Matrix;
  187.    friend nricMatrix;
  188.    friend RowVector;
  189.    friend ColumnVector;
  190.    friend SymmetricMatrix;
  191.    friend UpperTriangularMatrix;
  192.    friend LowerTriangularMatrix;
  193.    friend DiagonalMatrix;
  194.    friend CroutMatrix;
  195.    friend AddedMatrix;
  196.    friend MultipliedMatrix;
  197.    friend SubtractedMatrix;
  198.    friend SolvedMatrix;
  199.    friend ShiftedMatrix;
  200.    friend ScaledMatrix;
  201.    friend TransposedMatrix;
  202.    friend NegatedMatrix;
  203.    friend InvertedMatrix;
  204.    friend RowedMatrix;
  205.    friend ColedMatrix;
  206.    friend DiagedMatrix;
  207.    friend MatedMatrix;
  208.    friend GetSubMatrix;
  209.    friend ConstMatrix;
  210.    friend ReturnMatrix;
  211. };
  212.  
  213.  
  214. /******************************* working classes **************************/
  215.  
  216. class GeneralMatrix : public BaseMatrix         // declarable matrix types
  217. {
  218. protected:
  219.    int tag;                                     // shows whether can reuse
  220.    int nrows, ncols;                            // dimensions
  221.    int storage;                                 // total store required
  222.    real* store;                                 // point to store (0=not set)
  223.    GeneralMatrix();                             // initialise with no store
  224.    GeneralMatrix(int);                          // constructor getting store
  225.    void Add(GeneralMatrix*, real);              // sum of GM and real
  226.    void Add(real);                              // add real to this
  227.    void Multiply(GeneralMatrix*, real);         // product of GM and real
  228.    void Multiply(real);                         // multiply this by real
  229.    void Negate(GeneralMatrix*);                 // change sign
  230.    void Negate();                               // change sign
  231.    void operator=(real);                        // set matrix to constant
  232.    real* GetStore();                            // get store or copy
  233.    GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
  234.                                                 // temporarily access store
  235.    void GetMatrix(GeneralMatrix*);              // used by = and initialise
  236. #ifndef __ZTC__
  237.    void GetMatrixC(const GeneralMatrix*);       // used by = and initialise
  238. #endif
  239.    void Eq(BaseMatrix&, MatrixType);            // used by =
  240.    int search(const GeneralMatrix*) const;
  241.    virtual GeneralMatrix* Transpose(MatrixType);
  242.    void CheckConversion(const BaseMatrix&);     // check conversion OK
  243.    void ReDimension(int, int, int);             // change dimensions
  244.    int NrowsV() const;                          // get dimensions
  245.    int NcolsV() const;                          // virtual version
  246. public:
  247.    GeneralMatrix* Evaluate(MatrixType);
  248.    MatrixType Type() const = 0;                 // type of a matrix
  249.    int Nrows() const { return nrows; }          // get dimensions
  250.    int Ncols() const { return ncols; }
  251.    int Storage() const { return storage; }
  252.    real* Store() const { return store; }
  253.    virtual ~GeneralMatrix();                    // delete store if set
  254.    void tDelete();                              // delete if tag permits
  255.    BOOL reuse();                                // TRUE if tag allows reuse
  256.    void Protect() { tag=-1; }                   // can't delete or reuse
  257.    int Tag() const { return tag; }
  258.    BOOL IsZero() const;                         // test matrix has all zeros
  259.    void Release() { tag=1; }                    // del store after next use
  260.    void Release(int t) { tag=t; }               // del store after t accesses
  261.    void ReleaseAndDelete() { tag=0; }           // delete matrix after use
  262.    void operator<<(const real*);                // assignment from an array
  263.    void operator<<(BaseMatrix& X) { Eq(X,this->Type()); }
  264.                                                 // = without checking type
  265.    void Inject(const GeneralMatrix&);           // copy stored els only
  266.    virtual GeneralMatrix* MakeSolver();         // for solving
  267.    virtual void Solver(MatrixRowCol&, const MatrixRowCol&) {}
  268.    virtual void GetRow(MatrixRowCol&) = 0;      // Get matrix row
  269.    virtual void RestoreRow(MatrixRowCol&) {}    // Restore matrix row
  270.    virtual void NextRow(MatrixRowCol&);         // Go to next row
  271.    virtual void GetCol(MatrixRowCol&) = 0;      // Get matrix col
  272.    virtual void RestoreCol(MatrixRowCol&) {}    // Restore matrix col
  273.    virtual void NextCol(MatrixRowCol&);         // Go to next col
  274.    real SumSquare();
  275.    real SumAbsoluteValue();
  276.    real MaximumAbsoluteValue();
  277.    LogAndSign LogDeterminant();
  278.    ConstMatrix c() const;                       // to access constant matrices
  279.    void CheckStore() const;                     // check store is non-zero
  280.  
  281.    friend Matrix;
  282.    friend nricMatrix;
  283.    friend SymmetricMatrix;
  284.    friend UpperTriangularMatrix;
  285.    friend LowerTriangularMatrix;
  286.    friend DiagonalMatrix;
  287.    friend CroutMatrix;
  288.    friend RowVector;
  289.    friend ColumnVector;
  290.    friend BaseMatrix;
  291.    friend AddedMatrix;
  292.    friend MultipliedMatrix;
  293.    friend SubtractedMatrix;
  294.    friend SolvedMatrix;
  295.    friend ShiftedMatrix;
  296.    friend ScaledMatrix;
  297.    friend TransposedMatrix;
  298.    friend NegatedMatrix;
  299.    friend InvertedMatrix;
  300.    friend RowedMatrix;
  301.    friend ColedMatrix;
  302.    friend DiagedMatrix;
  303.    friend MatedMatrix;
  304.    friend GetSubMatrix;
  305.    friend ConstMatrix;
  306.    friend ReturnMatrix;
  307. };
  308.  
  309. class Matrix : public GeneralMatrix             // usual rectangular matrix
  310. {
  311. public:
  312.    Matrix() {}
  313.    Matrix(int, int);                            // standard declaration
  314.    Matrix(BaseMatrix&);                         // evaluate BaseMatrix
  315.    void operator=(BaseMatrix&);
  316.    void operator=(real f) { GeneralMatrix::operator=(f); }
  317.    MatrixType Type() const;
  318.    real& operator()(int, int);                  // access element
  319.    real& element(int, int);                     // access element
  320.    Matrix(Matrix& gm) { GetMatrix(&gm); }
  321. #ifndef __ZTC__
  322.    real operator()(int, int) const;             // access element
  323.    Matrix(const Matrix& gm) { GetMatrixC(&gm); }
  324. #endif
  325.    GeneralMatrix* MakeSolver();
  326.    real Trace();
  327.    void GetRow(MatrixRowCol&);
  328.    void GetCol(MatrixRowCol&);
  329.    void RestoreCol(MatrixRowCol&);
  330.    void NextRow(MatrixRowCol&);
  331.    void NextCol(MatrixRowCol&);
  332.    void ReDimension(int,int);                   // change dimensions
  333. };
  334.  
  335. class nricMatrix : public Matrix                // for use with Numerical
  336.                                                 // Recipes in C
  337. {
  338.    real** row_pointer;                          // points to rows
  339.    void MakeRowPointer();                       // build rowpointer
  340.    void DeleteRowPointer();
  341. public:
  342.    nricMatrix() {}
  343.    nricMatrix(int m, int n)                     // standard declaration
  344.       :  Matrix(m,n) { MakeRowPointer(); }
  345.    nricMatrix(BaseMatrix& bm)                   // evaluate BaseMatrix
  346.       :  Matrix(bm) { MakeRowPointer(); }
  347.    void operator=(BaseMatrix& bm)
  348.       { DeleteRowPointer(); Matrix::operator=(bm); MakeRowPointer(); }
  349.    void operator=(real f) { GeneralMatrix::operator=(f); }
  350.    void operator<<(BaseMatrix& X)
  351.       { DeleteRowPointer(); Eq(X,this->Type()); MakeRowPointer(); }
  352.    nricMatrix(nricMatrix& gm) { GetMatrix(&gm); MakeRowPointer(); }
  353. #ifndef __ZTC__
  354.    nricMatrix(const nricMatrix& gm) { GetMatrixC(&gm); MakeRowPointer(); }
  355. #endif
  356.    void ReDimension(int m, int n)               // change dimensions
  357.       { DeleteRowPointer(); Matrix::ReDimension(m,n); MakeRowPointer(); }
  358.    ~nricMatrix() { DeleteRowPointer(); }
  359. #ifndef __ZTC__
  360.    operator real**() const { CheckStore(); return row_pointer-1; }
  361. #endif
  362. };
  363.  
  364. class SymmetricMatrix : public GeneralMatrix
  365. {
  366. public:
  367.    SymmetricMatrix() {}
  368.    SymmetricMatrix(int);
  369.    SymmetricMatrix(BaseMatrix&);
  370.    void operator=(BaseMatrix&);
  371.    void operator=(real f) { GeneralMatrix::operator=(f); }
  372.    real& operator()(int, int);                  // access element
  373.    real& element(int, int);                     // access element
  374.    MatrixType Type() const;
  375.    SymmetricMatrix(SymmetricMatrix& gm) { GetMatrix(&gm); }
  376. #ifndef __ZTC__
  377.    real operator()(int, int) const;             // access element
  378.    SymmetricMatrix(const SymmetricMatrix& gm) { GetMatrixC(&gm); }
  379. #endif
  380.    real SumSquare();
  381.    real SumAbsoluteValue();
  382.    real Trace();
  383.    void GetRow(MatrixRowCol&);
  384.    void GetCol(MatrixRowCol&);
  385.    GeneralMatrix* Transpose(MatrixType);
  386.    void ReDimension(int);                       // change dimensions
  387. };
  388.  
  389. class UpperTriangularMatrix : public GeneralMatrix
  390. {
  391. public:
  392.    UpperTriangularMatrix() {}
  393.    UpperTriangularMatrix(int);
  394.    void operator=(BaseMatrix&);
  395.    UpperTriangularMatrix(BaseMatrix&);
  396.    UpperTriangularMatrix(UpperTriangularMatrix& gm) { GetMatrix(&gm); }
  397. #ifndef __ZTC__
  398.    real operator()(int, int) const;             // access element
  399.    UpperTriangularMatrix(const UpperTriangularMatrix& gm) { GetMatrixC(&gm); }
  400. #endif
  401.    void operator=(real f) { GeneralMatrix::operator=(f); }
  402.    real& operator()(int, int);                  // access element
  403.    real& element(int, int);                     // access element
  404.    MatrixType Type() const;
  405.    GeneralMatrix* MakeSolver() { return this; } // for solving
  406.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  407.    LogAndSign LogDeterminant();
  408.    real Trace();
  409.    void GetRow(MatrixRowCol&);
  410.    void GetCol(MatrixRowCol&);
  411.    void RestoreCol(MatrixRowCol&);
  412.    void NextRow(MatrixRowCol&);
  413.    void ReDimension(int);                       // change dimensions
  414. };
  415.  
  416. class LowerTriangularMatrix : public GeneralMatrix
  417. {
  418. public:
  419.    LowerTriangularMatrix() {}
  420.    LowerTriangularMatrix(int);
  421.    LowerTriangularMatrix(LowerTriangularMatrix& gm) { GetMatrix(&gm); }
  422. #ifndef __ZTC__
  423.    real operator()(int, int) const;             // access element
  424.    LowerTriangularMatrix(const LowerTriangularMatrix& gm) { GetMatrixC(&gm); }
  425. #endif
  426.    LowerTriangularMatrix(BaseMatrix& M);
  427.    void operator=(BaseMatrix&);
  428.    void operator=(real f) { GeneralMatrix::operator=(f); }
  429.    real& operator()(int, int);                  // access element
  430.    real& element(int, int);                     // access element
  431.    MatrixType Type() const;
  432.    GeneralMatrix* MakeSolver() { return this; } // for solving
  433.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  434.    LogAndSign LogDeterminant();
  435.    real Trace();
  436.    void GetRow(MatrixRowCol&);
  437.    void GetCol(MatrixRowCol&);
  438.    void RestoreCol(MatrixRowCol&);
  439.    void NextRow(MatrixRowCol&);
  440.    void ReDimension(int);                       // change dimensions
  441. };
  442.  
  443. class DiagonalMatrix : public GeneralMatrix
  444. {
  445. public:
  446.    DiagonalMatrix() {}
  447.    DiagonalMatrix(int);
  448.    DiagonalMatrix(BaseMatrix&);
  449.    DiagonalMatrix(DiagonalMatrix& gm) { GetMatrix(&gm); }
  450. #ifndef __ZTC__
  451.    real operator()(int, int) const;             // access element
  452.    real operator()(int) const;
  453.    DiagonalMatrix(const DiagonalMatrix& gm) { GetMatrixC(&gm); }
  454. #endif
  455.    void operator=(BaseMatrix&);
  456.    void operator=(real f) { GeneralMatrix::operator=(f); }
  457.    real& operator()(int, int);                  // access element
  458.    real& operator()(int);                       // access element
  459.    real& element(int, int);                     // access element
  460.    real& element(int);                          // access element
  461.    MatrixType Type() const;
  462.  
  463.    LogAndSign LogDeterminant();
  464.    real Trace();
  465.    void GetRow(MatrixRowCol&);
  466.    void GetCol(MatrixRowCol&);
  467.    void NextRow(MatrixRowCol&);
  468.    void NextCol(MatrixRowCol&);
  469.    GeneralMatrix* MakeSolver() { return this; } // for solving
  470.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  471.    GeneralMatrix* Transpose(MatrixType);
  472.    void ReDimension(int);                       // change dimensions
  473. #ifndef __ZTC__
  474.    operator real*() const
  475.       { CheckStore(); return store-1; }         // for use by NRIC
  476. #endif
  477. };
  478.  
  479. class RowVector : public Matrix
  480. {
  481. public:
  482.    RowVector() {}
  483.    RowVector(int n) : Matrix(1,n) {}
  484.    RowVector(BaseMatrix&);
  485.    RowVector(RowVector& gm) { GetMatrix(&gm); }
  486. #ifndef __ZTC__
  487.    real operator()(int) const;                  // access element
  488.    RowVector(const RowVector& gm) {GetMatrixC(&gm); }
  489. #endif
  490.    void operator=(BaseMatrix&);
  491.    void operator=(real f) { GeneralMatrix::operator=(f); }
  492.    real& operator()(int);                       // access element
  493.    real& element(int);                          // access element
  494.    MatrixType Type() const;
  495.    void GetCol(MatrixRowCol&);
  496.    void NextCol(MatrixRowCol&);
  497.    void RestoreCol(MatrixRowCol&);
  498.    GeneralMatrix* Transpose(MatrixType);
  499.    void ReDimension(int);                       // change dimensions
  500. #ifndef __ZTC__
  501.    operator real*() const
  502.       { CheckStore(); return store-1; }         // for use by NRIC
  503. #endif
  504. };
  505.  
  506. class ColumnVector : public Matrix
  507. {
  508. public:
  509.    ColumnVector() {}
  510.    ColumnVector(int n) : Matrix(n,1) {}
  511.    ColumnVector(BaseMatrix&);
  512.    ColumnVector(ColumnVector& gm) { GetMatrix(&gm); }
  513. #ifndef __ZTC__
  514.    real operator()(int) const;                  // access element
  515.    ColumnVector(const ColumnVector& gm) { GetMatrixC(&gm); }
  516. #endif
  517.    void operator=(BaseMatrix&);
  518.    void operator=(real f) { GeneralMatrix::operator=(f); }
  519.    real& operator()(int);                       // access element
  520.    real& element(int);                          // access element
  521.    MatrixType Type() const;
  522.    GeneralMatrix* Transpose(MatrixType);
  523.    void ReDimension(int);                       // change dimensions
  524. #ifndef __ZTC__
  525.    operator real*() const
  526.       { CheckStore(); return store-1; }         // for use by NRIC
  527. #endif
  528. };
  529.  
  530. class CroutMatrix : public GeneralMatrix        // for LU decomposition
  531. {
  532.    int* indx;
  533.    BOOL d;
  534.    void ludcmp();
  535. public:
  536.    CroutMatrix(BaseMatrix&);
  537.    MatrixType Type() const;
  538.    void lubksb(real*, int=0);
  539.    ~CroutMatrix() { delete [nrows] indx; }
  540.    GeneralMatrix* MakeSolver() { return this; } // for solving
  541.    LogAndSign LogDeterminant();
  542.    void Solver(MatrixRowCol&, const MatrixRowCol&);
  543.    void GetRow(MatrixRowCol&) { MatrixError("GetRow not defined for Crout"); }
  544.    void GetCol(MatrixRowCol&) { MatrixError("GetCol not defined for Crout"); }
  545.    void operator=(BaseMatrix&)
  546.       { MatrixError("operator= not defined for Crout"); }
  547. };
  548.  
  549.  
  550. /***************************** temporary classes *************************/
  551.  
  552. class AddedMatrix : public BaseMatrix
  553. {
  554. protected:
  555.    BaseMatrix* bm1;                             // pointers to summands
  556.    BaseMatrix* bm2;
  557.    AddedMatrix(BaseMatrix* bm1x, BaseMatrix* bm2x) : bm1(bm1x),bm2(bm2x) {}
  558.    int search(const GeneralMatrix*) const;
  559.    int NrowsV() const;
  560.    int NcolsV() const;
  561. private:
  562.    MatrixType Type() const;
  563.    friend BaseMatrix;
  564. public:
  565.    GeneralMatrix* Evaluate(MatrixType);
  566. };
  567.  
  568. class MultipliedMatrix : public AddedMatrix
  569. {
  570.    MultipliedMatrix(BaseMatrix* bm1x, BaseMatrix* bm2x)
  571.       : AddedMatrix(bm1x,bm2x) {}
  572.    MatrixType Type() const;
  573.    friend BaseMatrix;
  574. public:
  575.    GeneralMatrix* Evaluate(MatrixType);
  576. };
  577.  
  578. class SolvedMatrix : public AddedMatrix
  579. {
  580.    SolvedMatrix(BaseMatrix* bm1x, BaseMatrix* bm2x)
  581.       : AddedMatrix(bm1x,bm2x) {}
  582.    MatrixType Type() const;
  583.    friend BaseMatrix;
  584.    friend InvertedMatrix;                        // for operator*
  585. public:
  586.    GeneralMatrix* Evaluate(MatrixType);
  587. };
  588.  
  589. class SubtractedMatrix : public AddedMatrix
  590. {
  591.    SubtractedMatrix(BaseMatrix* bm1x, BaseMatrix* bm2x)
  592.       : AddedMatrix(bm1x,bm2x) {}
  593.    MatrixType Type() const;
  594.    friend BaseMatrix;
  595. public:
  596.    GeneralMatrix* Evaluate(MatrixType);
  597. };
  598.  
  599. class ShiftedMatrix : public BaseMatrix
  600. {
  601. protected:
  602.    real f;
  603.    BaseMatrix* bm;
  604.    ShiftedMatrix(BaseMatrix* bmx, real fx) : bm(bmx),f(fx) {}
  605.    int search(const GeneralMatrix*) const;
  606.    int NrowsV() const;
  607.    int NcolsV() const;
  608. private:
  609.    MatrixType Type() const;
  610.    friend BaseMatrix;
  611. public:
  612.    GeneralMatrix* Evaluate(MatrixType);
  613. };
  614.  
  615. class ScaledMatrix : public ShiftedMatrix
  616. {
  617.    ScaledMatrix(BaseMatrix* bmx, real fx) : ShiftedMatrix(bmx,fx) {}
  618.    MatrixType Type() const;
  619.    friend BaseMatrix;
  620. public:
  621.    GeneralMatrix* Evaluate(MatrixType);
  622. };
  623.  
  624. class NegatedMatrix : public BaseMatrix
  625. {
  626. protected:
  627.    BaseMatrix* bm;
  628.    NegatedMatrix(BaseMatrix* bmx) : bm(bmx) {}
  629.    int search(const GeneralMatrix*) const;
  630.    int NrowsV() const;
  631.    int NcolsV() const;
  632. private:
  633.    MatrixType Type() const;
  634.    friend BaseMatrix;
  635. public:
  636.    GeneralMatrix* Evaluate(MatrixType);
  637. };
  638.  
  639. class TransposedMatrix : public NegatedMatrix
  640. {
  641.    TransposedMatrix(BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  642.    int NrowsV() const;
  643.    int NcolsV() const;
  644.    MatrixType Type() const;
  645.    friend BaseMatrix;
  646. public:
  647.    GeneralMatrix* Evaluate(MatrixType);
  648. };
  649.  
  650. class InvertedMatrix : public NegatedMatrix
  651. {
  652.    InvertedMatrix(BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  653.    MatrixType Type() const;
  654. public:
  655.    SolvedMatrix operator*(BaseMatrix&);         // inverse(A) * B
  656.    friend BaseMatrix;
  657.    GeneralMatrix* Evaluate(MatrixType);
  658. };
  659.  
  660. class RowedMatrix : public NegatedMatrix
  661. {
  662.    MatrixType Type() const;
  663.    int NrowsV() const;
  664.    int NcolsV() const;
  665.    RowedMatrix(BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  666.    friend BaseMatrix;
  667. public:
  668.    GeneralMatrix* Evaluate(MatrixType);
  669. };
  670.  
  671. class ColedMatrix : public NegatedMatrix
  672. {
  673.    MatrixType Type() const;
  674.    int NrowsV() const;
  675.    int NcolsV() const;
  676.    ColedMatrix(BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  677.    friend BaseMatrix;
  678. public:
  679.    GeneralMatrix* Evaluate(MatrixType);
  680. };
  681.  
  682. class DiagedMatrix : public NegatedMatrix
  683. {
  684.    MatrixType Type() const;
  685.    int NrowsV() const;
  686.    int NcolsV() const;
  687.    DiagedMatrix(BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  688.    friend BaseMatrix;
  689. public:
  690.    GeneralMatrix* Evaluate(MatrixType);
  691. };
  692.  
  693. class MatedMatrix : public NegatedMatrix
  694. {
  695.    int nr, nc;
  696.    MatrixType Type() const;
  697.    int NrowsV() const;
  698.    int NcolsV() const;
  699.    MatedMatrix(BaseMatrix* bmx, int nrx, int ncx)
  700.       : NegatedMatrix(bmx), nr(nrx), nc(ncx) {}
  701.    friend BaseMatrix;
  702. public:
  703.    GeneralMatrix* Evaluate(MatrixType);
  704. };
  705.  
  706. class ConstMatrix : public BaseMatrix
  707. {
  708.    const GeneralMatrix* cgm;
  709.    MatrixType Type() const;
  710.    int NrowsV() const;
  711.    int NcolsV() const;
  712.    int search(const GeneralMatrix*) const;
  713.    ConstMatrix(const GeneralMatrix* cgmx) : cgm(cgmx) {}
  714.    friend BaseMatrix;
  715.    friend GeneralMatrix;
  716. public:
  717.    GeneralMatrix* Evaluate(MatrixType);
  718. };
  719.  
  720. class ReturnMatrix : public BaseMatrix    // for matrix return
  721. {
  722.    GeneralMatrix* gm;
  723.    MatrixType Type() const;
  724.    int NrowsV() const;
  725.    int NcolsV() const;
  726.    int search(const GeneralMatrix*) const;
  727. public:
  728.    GeneralMatrix* Evaluate(MatrixType);
  729.    friend BaseMatrix;
  730.    ReturnMatrix(const ReturnMatrix& tm) : gm(tm.gm) {}
  731.    ReturnMatrix(GeneralMatrix&);
  732. };
  733.  
  734.  
  735. /**************************** submatrices ******************************/
  736.  
  737. class GetSubMatrix : public NegatedMatrix
  738. {
  739.    int row_skip;
  740.    int row_number;
  741.    int col_skip;
  742.    int col_number;
  743.    MatrixType mt;
  744.  
  745.    GetSubMatrix
  746.       (BaseMatrix* bmx, int rs, int rn, int cs, int cn, MatrixType mtx)
  747.       : NegatedMatrix(bmx),
  748.       row_skip(rs), row_number(rn), col_skip(cs), col_number(cn), mt(mtx) {}
  749.    GetSubMatrix(const GetSubMatrix& g)
  750.       : NegatedMatrix(g.bm), row_skip(g.row_skip), row_number(g.row_number),
  751.       col_skip(g.col_skip), col_number(g.col_number), mt(g.mt) {}
  752.    MatrixType Type() const;
  753.    int NrowsV() const;
  754.    int NcolsV() const;
  755.    friend BaseMatrix;
  756. public:
  757.    GeneralMatrix* Evaluate(MatrixType);
  758.    void operator=(BaseMatrix&)
  759.       { MatrixError("= not defined for submatrix; use <<"); }
  760.    void operator<<(BaseMatrix&);
  761.    void operator<<(const real*);                // copy from array
  762.    void operator<<(real);                       // copy from constant
  763.    void Inject(const GeneralMatrix&);           // copy stored els only
  764. };
  765.  
  766. /***************************** functions ***********************************/
  767.  
  768.  
  769. inline LogAndSign LogDeterminant(BaseMatrix& B) { return B.LogDeterminant(); }
  770. inline real SumSquare(BaseMatrix& B) { return B.SumSquare(); }
  771. inline real Trace(BaseMatrix& B) { return B.Trace(); }
  772. inline real SumAbsoluteValue(BaseMatrix& B) { return B.SumAbsoluteValue(); }
  773. inline real MaximumAbsoluteValue(BaseMatrix& B)
  774.    { return B.MaximumAbsoluteValue(); }
  775. inline real Norm1(BaseMatrix& B) { return B.Norm1(); }
  776. inline real Norm1(RowVector& RV) { return RV.MaximumAbsoluteValue(); }
  777. inline real NormInfinity(BaseMatrix& B) { return B.NormInfinity(); }
  778. inline real NormInfinity(ColumnVector& CV)
  779.    { return CV.MaximumAbsoluteValue(); } 
  780.  
  781. #endif
  782.