Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

newmat.h

Go to the documentation of this file.
00001 //$$ newmat.h           definition file for new version of matrix package
00002 
00003 // Copyright (C) 1991,2,3,4,7,2000,2,3: R B Davies
00004 
00005 #ifndef NEWMAT_LIB
00006 #define NEWMAT_LIB 0
00007 
00008 #include "include.h"
00009 
00010 #include "myexcept.h"
00011 
00012 
00013 #ifdef use_namespace
00014 namespace NEWMAT { using namespace RBD_COMMON; }
00015 namespace RBD_LIBRARIES { using namespace NEWMAT; }
00016 namespace NEWMAT {
00017 #endif
00018 
00019 //#define DO_REPORT                     // to activate REPORT
00020 
00021 #ifdef NO_LONG_NAMES
00022 #define UpperTriangularMatrix UTMatrix
00023 #define LowerTriangularMatrix LTMatrix
00024 #define SymmetricMatrix SMatrix
00025 #define DiagonalMatrix DMatrix
00026 #define BandMatrix BMatrix
00027 #define UpperBandMatrix UBMatrix
00028 #define LowerBandMatrix LBMatrix
00029 #define SymmetricBandMatrix SBMatrix
00030 #define BandLUMatrix BLUMatrix
00031 #endif
00032 
00033 // ************************** general utilities ****************************/
00034 
00035 class GeneralMatrix;
00036 
00037 void MatrixErrorNoSpace(const void*);                 // no space handler
00038 
00039 class LogAndSign
00040 // Return from LogDeterminant function
00041 //    - value of the log plus the sign (+, - or 0)
00042 {
00043    Real log_value;
00044    int sign;
00045 public:
00046    LogAndSign() : log_value(0), sign(1) {}
00047    LogAndSign(Real);
00048    void operator*=(Real);
00049    void PowEq(int k);  // raise to power of k
00050    void ChangeSign() { sign = -sign; }
00051    Real LogValue() const { return log_value; }
00052    int Sign() const { return sign; }
00053    Real Value() const;
00054    FREE_CHECK(LogAndSign)
00055 };
00056 
00057 // the following class is for counting the number of times a piece of code
00058 // is executed. It is used for locating any code not executed by test
00059 // routines. Use turbo GREP locate all places this code is called and
00060 // check which ones are not accessed.
00061 // Somewhat implementation dependent as it relies on "cout" still being
00062 // present when ExeCounter objects are destructed.
00063 
00064 #ifdef DO_REPORT
00065 
00066 class ExeCounter
00067 {
00068    int line;                                    // code line number
00069    int fileid;                                  // file identifier
00070    long nexe;                                   // number of executions
00071    static int nreports;                         // number of reports
00072 public:
00073    ExeCounter(int,int);
00074    void operator++() { nexe++; }
00075    ~ExeCounter();                               // prints out reports
00076 };
00077 
00078 #endif
00079 
00080 
00081 // ************************** class MatrixType *****************************/
00082 
00083 // Is used for finding the type of a matrix resulting from the binary operations
00084 // +, -, * and identifying what conversions are permissible.
00085 // This class must be updated when new matrix types are added.
00086 
00087 class GeneralMatrix;                            // defined later
00088 class BaseMatrix;                               // defined later
00089 class MatrixInput;                              // defined later
00090 
00091 class MatrixType
00092 {
00093 public:
00094    enum Attribute {  Valid     = 1,
00095                      Diagonal  = 2,             // order of these is important
00096                      Symmetric = 4,
00097                      Band      = 8,
00098                      Lower     = 16,
00099                      Upper     = 32,
00100                      Square    = 64,
00101                      Skew      = 128,
00102                      LUDeco    = 256,
00103                      Ones      = 512 };
00104 
00105    enum            { US = 0,
00106                      UT = Valid + Upper + Square,
00107                      LT = Valid + Lower + Square,
00108                      Rt = Valid,
00109                      Sq = Valid + Square,
00110                      Sm = Valid + Symmetric + Square,
00111                      Sk = Valid + Skew + Square,
00112                      Dg = Valid + Diagonal + Band + Lower + Upper + Symmetric
00113                         + Square,
00114                      Id = Valid + Diagonal + Band + Lower + Upper + Symmetric
00115                         + Square + Ones,
00116                      RV = Valid,     //   do not separate out
00117                      CV = Valid,     //   vectors
00118                      BM = Valid + Band + Square,
00119                      UB = Valid + Band + Upper + Square,
00120                      LB = Valid + Band + Lower + Square,
00121                      SB = Valid + Band + Symmetric + Square,
00122                      Ct = Valid + LUDeco + Square,
00123                      BC = Valid + Band + LUDeco + Square,
00124                      Mask = ~Square
00125                    };
00126 
00127 
00128    static int nTypes() { return 12; }          // number of different types
00129                  // exclude Ct, US, BC
00130 public:
00131    int attribute;
00132    bool DataLossOK;                            // true if data loss is OK when
00133                                                // this represents a destination
00134 public:
00135    MatrixType () : attribute(0), DataLossOK(false) {}
00136    MatrixType (int n) : attribute(n), DataLossOK(false) {}
00137    MatrixType (int n, bool dlok) : attribute(n), DataLossOK(dlok) {}
00138    MatrixType (const MatrixType& mt)
00139       : attribute(mt.attribute), DataLossOK(mt.DataLossOK) {}
00140    MatrixType& operator=(const MatrixType& mt)
00141      { attribute = mt.attribute; DataLossOK = mt.DataLossOK; return *this; }
00142    void SetDataLossOK() { DataLossOK = true; }
00143    int operator+() const { return attribute; }
00144    MatrixType operator+(MatrixType mt) const
00145       { return MatrixType(attribute & mt.attribute); }
00146    MatrixType operator*(const MatrixType&) const;
00147    MatrixType SP(const MatrixType&) const;
00148    MatrixType KP(const MatrixType&) const;
00149    MatrixType operator|(const MatrixType& mt) const
00150       { return MatrixType(attribute & mt.attribute & Valid); }
00151    MatrixType operator&(const MatrixType& mt) const
00152       { return MatrixType(attribute & mt.attribute & Valid); }
00153    bool operator>=(MatrixType mt) const
00154       { return ( attribute & ~mt.attribute & Mask ) == 0; }
00155    bool operator<(MatrixType mt) const         // for MS Visual C++ 4
00156       { return ( attribute & ~mt.attribute & Mask ) != 0; }
00157    bool operator==(MatrixType mt) const
00158       { return (attribute == mt.attribute); }
00159    bool operator!=(MatrixType mt) const
00160       { return (attribute != mt.attribute); }
00161    bool operator!() const { return (attribute & Valid) == 0; }
00162    MatrixType i() const;                       // type of inverse
00163    MatrixType t() const;                       // type of transpose
00164    MatrixType AddEqualEl() const               // Add constant to matrix
00165       { return MatrixType(attribute & (Valid + Symmetric + Square)); }
00166    MatrixType MultRHS() const;                 // type for rhs of multiply
00167    MatrixType sub() const                      // type of submatrix
00168       { return MatrixType(attribute & Valid); }
00169    MatrixType ssub() const                     // type of sym submatrix
00170       { return MatrixType(attribute); }        // not for selection matrix
00171    GeneralMatrix* New() const;                 // new matrix of given type
00172    GeneralMatrix* New(int,int,BaseMatrix*) const;
00173                                                // new matrix of given type
00174    const char* Value() const;                  // to print type
00175    friend bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
00176    friend bool Compare(const MatrixType&, MatrixType&);
00177                                                // compare and check conv.
00178    bool IsBand() const { return (attribute & Band) != 0; }
00179    bool IsDiagonal() const { return (attribute & Diagonal) != 0; }
00180    bool IsSymmetric() const { return (attribute & Symmetric) != 0; }
00181    bool CannotConvert() const { return (attribute & LUDeco) != 0; }
00182                                                // used by operator== 
00183    FREE_CHECK(MatrixType)
00184 };
00185 
00186 
00187 // *********************** class MatrixBandWidth ***********************/
00188 
00189 class MatrixBandWidth
00190 {
00191 public:
00192    int lower;
00193    int upper;
00194    MatrixBandWidth(const int l, const int u) : lower(l), upper (u) {}
00195    MatrixBandWidth(const int n) : lower(n), upper(n) {}
00196    MatrixBandWidth operator+(const MatrixBandWidth&) const;
00197    MatrixBandWidth operator*(const MatrixBandWidth&) const;
00198    MatrixBandWidth minimum(const MatrixBandWidth&) const;
00199    MatrixBandWidth t() const { return MatrixBandWidth(upper,lower); }
00200    bool operator==(const MatrixBandWidth& bw) const
00201       { return (lower == bw.lower) && (upper == bw.upper); }
00202    bool operator!=(const MatrixBandWidth& bw) const { return !operator==(bw); }
00203    int Upper() const { return upper; }
00204    int Lower() const { return lower; }
00205    FREE_CHECK(MatrixBandWidth)
00206 };
00207 
00208 
00209 // ********************* Array length specifier ************************/
00210 
00211 // This class is introduced to avoid constructors such as
00212 //   ColumnVector(int)
00213 // being used for conversions
00214 
00215 class ArrayLengthSpecifier
00216 {
00217    int value;
00218 public:
00219    int Value() const { return value; }
00220    ArrayLengthSpecifier(int l) : value(l) {}
00221 };
00222 
00223 // ************************* Matrix routines ***************************/
00224 
00225 
00226 class MatrixRowCol;                             // defined later
00227 class MatrixRow;
00228 class MatrixCol;
00229 class MatrixColX;
00230 
00231 class GeneralMatrix;                            // defined later
00232 class AddedMatrix;
00233 class MultipliedMatrix;
00234 class SubtractedMatrix;
00235 class SPMatrix;
00236 class KPMatrix;
00237 class ConcatenatedMatrix;
00238 class StackedMatrix;
00239 class SolvedMatrix;
00240 class ShiftedMatrix;
00241 class NegShiftedMatrix;
00242 class ScaledMatrix;
00243 class TransposedMatrix;
00244 class ReversedMatrix;
00245 class NegatedMatrix;
00246 class InvertedMatrix;
00247 class RowedMatrix;
00248 class ColedMatrix;
00249 class DiagedMatrix;
00250 class MatedMatrix;
00251 class GetSubMatrix;
00252 class ReturnMatrix;
00253 class Matrix;
00254 class SquareMatrix;
00255 class nricMatrix;
00256 class RowVector;
00257 class ColumnVector;
00258 class SymmetricMatrix;
00259 class UpperTriangularMatrix;
00260 class LowerTriangularMatrix;
00261 class DiagonalMatrix;
00262 class CroutMatrix;
00263 class BandMatrix;
00264 class LowerBandMatrix;
00265 class UpperBandMatrix;
00266 class SymmetricBandMatrix;
00267 class LinearEquationSolver;
00268 class GenericMatrix;
00269 
00270 
00271 #define MatrixTypeUnSp 0
00272 //static MatrixType MatrixTypeUnSp(MatrixType::US);
00273 //            // AT&T needs this
00274 
00275 class BaseMatrix : public Janitor               // base of all matrix classes
00276 {
00277 protected:
00278    virtual int search(const BaseMatrix*) const = 0;
00279             // count number of times matrix
00280               // is referred to
00281 
00282 public:
00283    virtual GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp) = 0;
00284             // evaluate temporary
00285    // for old version of G++
00286    //   virtual GeneralMatrix* Evaluate(MatrixType mt) = 0;
00287    //   GeneralMatrix* Evaluate() { return Evaluate(MatrixTypeUnSp); }
00288    AddedMatrix operator+(const BaseMatrix&) const;    // results of operations
00289    MultipliedMatrix operator*(const BaseMatrix&) const;
00290    SubtractedMatrix operator-(const BaseMatrix&) const;
00291    ConcatenatedMatrix operator|(const BaseMatrix&) const;
00292    StackedMatrix operator&(const BaseMatrix&) const;
00293    ShiftedMatrix operator+(Real) const;
00294    ScaledMatrix operator*(Real) const;
00295    ScaledMatrix operator/(Real) const;
00296    ShiftedMatrix operator-(Real) const;
00297    TransposedMatrix t() const;
00298 //   TransposedMatrix t;
00299    NegatedMatrix operator-() const;                   // change sign of elements
00300    ReversedMatrix Reverse() const;
00301    InvertedMatrix i() const;
00302 //   InvertedMatrix i;
00303    RowedMatrix AsRow() const;
00304    ColedMatrix AsColumn() const;
00305    DiagedMatrix AsDiagonal() const;
00306    MatedMatrix AsMatrix(int,int) const;
00307    GetSubMatrix SubMatrix(int,int,int,int) const;
00308    GetSubMatrix SymSubMatrix(int,int) const;
00309    GetSubMatrix Row(int) const;
00310    GetSubMatrix Rows(int,int) const;
00311    GetSubMatrix Column(int) const;
00312    GetSubMatrix Columns(int,int) const;
00313    Real AsScalar() const;                      // conversion of 1 x 1 matrix
00314    virtual LogAndSign LogDeterminant() const;
00315    Real Determinant() const;
00316    virtual Real SumSquare() const;
00317    Real NormFrobenius() const;
00318    virtual Real SumAbsoluteValue() const;
00319    virtual Real Sum() const;
00320    virtual Real MaximumAbsoluteValue() const;
00321    virtual Real MaximumAbsoluteValue1(int& i) const;
00322    virtual Real MaximumAbsoluteValue2(int& i, int& j) const;
00323    virtual Real MinimumAbsoluteValue() const;
00324    virtual Real MinimumAbsoluteValue1(int& i) const;
00325    virtual Real MinimumAbsoluteValue2(int& i, int& j) const;
00326    virtual Real Maximum() const;
00327    virtual Real Maximum1(int& i) const;
00328    virtual Real Maximum2(int& i, int& j) const;
00329    virtual Real Minimum() const;
00330    virtual Real Minimum1(int& i) const;
00331    virtual Real Minimum2(int& i, int& j) const;
00332    virtual Real Trace() const;
00333    Real Norm1() const;
00334    Real NormInfinity() const;
00335    virtual MatrixBandWidth BandWidth() const;  // bandwidths of band matrix
00336    virtual void CleanUp() {}                   // to clear store
00337    void IEQND() const;                         // called by ineq. ops
00338 //   virtual ReturnMatrix Reverse() const;       // reverse order of elements
00339 //protected:
00340 //   BaseMatrix() : t(this), i(this) {}
00341 
00342    friend class GeneralMatrix;
00343    friend class Matrix;
00344    friend class SquareMatrix;
00345    friend class nricMatrix;
00346    friend class RowVector;
00347    friend class ColumnVector;
00348    friend class SymmetricMatrix;
00349    friend class UpperTriangularMatrix;
00350    friend class LowerTriangularMatrix;
00351    friend class DiagonalMatrix;
00352    friend class CroutMatrix;
00353    friend class BandMatrix;
00354    friend class LowerBandMatrix;
00355    friend class UpperBandMatrix;
00356    friend class SymmetricBandMatrix;
00357    friend class AddedMatrix;
00358    friend class MultipliedMatrix;
00359    friend class SubtractedMatrix;
00360    friend class SPMatrix;
00361    friend class KPMatrix;
00362    friend class ConcatenatedMatrix;
00363    friend class StackedMatrix;
00364    friend class SolvedMatrix;
00365    friend class ShiftedMatrix;
00366    friend class NegShiftedMatrix;
00367    friend class ScaledMatrix;
00368    friend class TransposedMatrix;
00369    friend class ReversedMatrix;
00370    friend class NegatedMatrix;
00371    friend class InvertedMatrix;
00372    friend class RowedMatrix;
00373    friend class ColedMatrix;
00374    friend class DiagedMatrix;
00375    friend class MatedMatrix;
00376    friend class GetSubMatrix;
00377    friend class ReturnMatrix;
00378    friend class LinearEquationSolver;
00379    friend class GenericMatrix;
00380    NEW_DELETE(BaseMatrix)
00381 };
00382 
00383 
00384 // ***************************** working classes **************************/
00385 
00386 class GeneralMatrix : public BaseMatrix         // declarable matrix types
00387 {
00388    virtual GeneralMatrix* Image() const;        // copy of matrix
00389 protected:
00390    int tag;                                     // shows whether can reuse
00391    int nrows_value, ncols_value;                // dimensions
00392    int storage;                                 // total store required
00393    Real* store;                                 // point to store (0=not set)
00394    GeneralMatrix();                             // initialise with no store
00395    GeneralMatrix(ArrayLengthSpecifier);         // constructor getting store
00396    GeneralMatrix(const GeneralMatrix&);
00397    void Add(GeneralMatrix*, Real);              // sum of GM and Real
00398    void Add(Real);                              // add Real to this
00399    void NegAdd(GeneralMatrix*, Real);           // Real - GM
00400    void NegAdd(Real);                           // this = this - Real
00401    void Multiply(GeneralMatrix*, Real);         // product of GM and Real
00402    void Multiply(Real);                         // multiply this by Real
00403    void Negate(GeneralMatrix*);                 // change sign
00404    void Negate();                               // change sign
00405    void ReverseElements();                      // internal reverse of elements
00406    void ReverseElements(GeneralMatrix*);        // reverse order of elements
00407    void operator=(Real);                        // set matrix to constant
00408    void operator=(const GeneralMatrix&);
00409    Real* GetStore();                            // get store or copy
00410    GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
00411                                                 // temporarily access store
00412    void GetMatrix(const GeneralMatrix*);        // used by = and initialise
00413    void Eq(const BaseMatrix&, MatrixType);      // used by =
00414    void Eq(const GeneralMatrix&);               // version with no conversion
00415    void Eq(const BaseMatrix&, MatrixType, bool);// used by <<
00416    void Eq2(const BaseMatrix&, MatrixType);     // cut down version of Eq
00417    int search(const BaseMatrix*) const;
00418    virtual GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
00419    void CheckConversion(const BaseMatrix&);     // check conversion OK
00420    void ReSize(int, int, int);                  // change dimensions
00421    virtual short SimpleAddOK(const GeneralMatrix* /*gm*/) { return 0; }
00422              // see bandmat.cpp for explanation
00423    virtual void MiniCleanUp()
00424       { store = 0; storage = 0; nrows_value = 0; ncols_value = 0; tag = -1;}
00425              // CleanUp when the data array has already been deleted
00426    void PlusEqual(const GeneralMatrix& gm);
00427    void MinusEqual(const GeneralMatrix& gm);
00428    void PlusEqual(Real f);
00429    void MinusEqual(Real f);
00430    void swap(GeneralMatrix& gm);                // swap values
00431 public:
00432    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
00433    virtual MatrixType Type() const = 0;         // type of a matrix
00434    int Nrows() const { return nrows_value; }    // get dimensions
00435    int Ncols() const { return ncols_value; }
00436    int Storage() const { return storage; }
00437    Real* Store() const { return store; }
00438    // updated names
00439    int nrows() const { return nrows_value; }    // get dimensions
00440    int ncols() const { return ncols_value; }
00441    int size() const { return storage; }
00442    Real* data() { return store; }
00443    const Real* data() const { return store; }
00444    const Real* const_data() const { return store; }
00445    virtual ~GeneralMatrix();                    // delete store if set
00446    void tDelete();                              // delete if tag permits
00447    bool reuse();                                // true if tag allows reuse
00448    void Protect() { tag=-1; }                   // cannot delete or reuse
00449    int Tag() const { return tag; }
00450    bool IsZero() const;                         // test matrix has all zeros
00451    void Release() { tag=1; }                    // del store after next use
00452    void Release(int tg) { tag=tg; }               // del store after t accesses
00453    void ReleaseAndDelete() { tag=0; }           // delete matrix after use
00454    void operator<<(const Real*);                // assignment from an array
00455    void operator<<(const int*);                // assignment from an array
00456    void operator<<(const BaseMatrix& X)
00457       { Eq(X,this->Type(),true); }              // = without checking type
00458    void Inject(const GeneralMatrix&);           // copy stored els only
00459    void operator+=(const BaseMatrix&);
00460    void operator-=(const BaseMatrix&);
00461    void operator*=(const BaseMatrix&);
00462    void operator|=(const BaseMatrix&);
00463    void operator&=(const BaseMatrix&);
00464    void operator+=(Real);
00465    void operator-=(Real r) { operator+=(-r); }
00466    void operator*=(Real);
00467    void operator/=(Real r) { operator*=(1/r); }
00468    virtual GeneralMatrix* MakeSolver();         // for solving
00469    virtual void Solver(MatrixColX&, const MatrixColX&) {}
00470    virtual void GetRow(MatrixRowCol&) = 0;      // Get matrix row
00471    virtual void RestoreRow(MatrixRowCol&) {}    // Restore matrix row
00472    virtual void NextRow(MatrixRowCol&);         // Go to next row
00473    virtual void GetCol(MatrixRowCol&) = 0;      // Get matrix col
00474    virtual void GetCol(MatrixColX&) = 0;        // Get matrix col
00475    virtual void RestoreCol(MatrixRowCol&) {}    // Restore matrix col
00476    virtual void RestoreCol(MatrixColX&) {}      // Restore matrix col
00477    virtual void NextCol(MatrixRowCol&);         // Go to next col
00478    virtual void NextCol(MatrixColX&);           // Go to next col
00479    Real SumSquare() const;
00480    Real SumAbsoluteValue() const;
00481    Real Sum() const;
00482    Real MaximumAbsoluteValue1(int& i) const;
00483    Real MinimumAbsoluteValue1(int& i) const;
00484    Real Maximum1(int& i) const;
00485    Real Minimum1(int& i) const;
00486    Real MaximumAbsoluteValue() const;
00487    Real MaximumAbsoluteValue2(int& i, int& j) const;
00488    Real MinimumAbsoluteValue() const;
00489    Real MinimumAbsoluteValue2(int& i, int& j) const;
00490    Real Maximum() const;
00491    Real Maximum2(int& i, int& j) const;
00492    Real Minimum() const;
00493    Real Minimum2(int& i, int& j) const;
00494    LogAndSign LogDeterminant() const;
00495    virtual bool IsEqual(const GeneralMatrix&) const;
00496                                                 // same type, same values
00497    void CheckStore() const;                     // check store is non-zero
00498    virtual void SetParameters(const GeneralMatrix*) {}
00499                                                 // set parameters in GetMatrix
00500    operator ReturnMatrix() const;               // for building a ReturnMatrix
00501    ReturnMatrix ForReturn() const;
00502    virtual bool SameStorageType(const GeneralMatrix& A) const;
00503    virtual void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
00504    virtual void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
00505    virtual void ReSize(const GeneralMatrix& A);
00506    MatrixInput operator<<(Real);                // for loading a list
00507    MatrixInput operator<<(int f);
00508 //   ReturnMatrix Reverse() const;                // reverse order of elements
00509    void CleanUp();                              // to clear store
00510 
00511    friend class Matrix;
00512    friend class SquareMatrix;
00513    friend class nricMatrix;
00514    friend class SymmetricMatrix;
00515    friend class UpperTriangularMatrix;
00516    friend class LowerTriangularMatrix;
00517    friend class DiagonalMatrix;
00518    friend class CroutMatrix;
00519    friend class RowVector;
00520    friend class ColumnVector;
00521    friend class BandMatrix;
00522    friend class LowerBandMatrix;
00523    friend class UpperBandMatrix;
00524    friend class SymmetricBandMatrix;
00525    friend class BaseMatrix;
00526    friend class AddedMatrix;
00527    friend class MultipliedMatrix;
00528    friend class SubtractedMatrix;
00529    friend class SPMatrix;
00530    friend class KPMatrix;
00531    friend class ConcatenatedMatrix;
00532    friend class StackedMatrix;
00533    friend class SolvedMatrix;
00534    friend class ShiftedMatrix;
00535    friend class NegShiftedMatrix;
00536    friend class ScaledMatrix;
00537    friend class TransposedMatrix;
00538    friend class ReversedMatrix;
00539    friend class NegatedMatrix;
00540    friend class InvertedMatrix;
00541    friend class RowedMatrix;
00542    friend class ColedMatrix;
00543    friend class DiagedMatrix;
00544    friend class MatedMatrix;
00545    friend class GetSubMatrix;
00546    friend class ReturnMatrix;
00547    friend class LinearEquationSolver;
00548    friend class GenericMatrix;
00549    NEW_DELETE(GeneralMatrix)
00550 };
00551 
00552 
00553 
00554 class Matrix : public GeneralMatrix             // usual rectangular matrix
00555 {
00556    GeneralMatrix* Image() const;                // copy of matrix
00557 public:
00558    Matrix() {}
00559    ~Matrix() {}
00560    Matrix(int, int);                            // standard declaration
00561    Matrix(const BaseMatrix&);                   // evaluate BaseMatrix
00562    void operator=(const BaseMatrix&);
00563    Matrix& operator=(Real f) { GeneralMatrix::operator=(f); return *this; }
00564    Matrix& operator=(const Matrix& m) { Eq(m); return *this; }
00565    MatrixType Type() const;
00566    Real& operator()(int, int);                  // access element
00567    Real& element(int, int);                     // access element
00568    Real operator()(int, int) const;            // access element
00569    Real element(int, int) const;               // access element
00570 #ifdef SETUP_C_SUBSCRIPTS
00571    Real* operator[](int m) { return store+m*ncols_value; }
00572    const Real* operator[](int m) const { return store+m*ncols_value; }
00573    // following for Numerical Recipes in C++
00574    Matrix(Real, int, int);
00575    Matrix(const Real*, int, int);
00576 #endif
00577    Matrix(const Matrix& gm) : GeneralMatrix() { GetMatrix(&gm); }
00578    GeneralMatrix* MakeSolver();
00579    Real Trace() const;
00580    void GetRow(MatrixRowCol&);
00581    void GetCol(MatrixRowCol&);
00582    void GetCol(MatrixColX&);
00583    void RestoreCol(MatrixRowCol&);
00584    void RestoreCol(MatrixColX&);
00585    void NextRow(MatrixRowCol&);
00586    void NextCol(MatrixRowCol&);
00587    void NextCol(MatrixColX&);
00588    virtual void ReSize(int,int);           // change dimensions
00589       // virtual so we will catch it being used in a vector called as a matrix
00590    void ReSize(const GeneralMatrix& A);
00591    Real MaximumAbsoluteValue2(int& i, int& j) const;
00592    Real MinimumAbsoluteValue2(int& i, int& j) const;
00593    Real Maximum2(int& i, int& j) const;
00594    Real Minimum2(int& i, int& j) const;
00595    void operator+=(const Matrix& M) { PlusEqual(M); }
00596    void operator-=(const Matrix& M) { MinusEqual(M); }
00597    void operator+=(Real f) { GeneralMatrix::Add(f); }
00598    void operator-=(Real f) { GeneralMatrix::Add(-f); }
00599    void swap(Matrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
00600    friend Real DotProduct(const Matrix& A, const Matrix& B);
00601    NEW_DELETE(Matrix)
00602 };
00603 
00604 class SquareMatrix : public Matrix              // square matrix
00605 {
00606    GeneralMatrix* Image() const;                // copy of matrix
00607 public:
00608    SquareMatrix() {}
00609    ~SquareMatrix() {}
00610    SquareMatrix(ArrayLengthSpecifier);          // standard declaration
00611    SquareMatrix(const BaseMatrix&);             // evaluate BaseMatrix
00612    void operator=(const BaseMatrix&);
00613    SquareMatrix& operator=(Real f) { GeneralMatrix::operator=(f); return *this; }
00614    SquareMatrix& operator=(const SquareMatrix& m) { Eq(m); return *this; }
00615    void operator=(const Matrix& m);
00616    MatrixType Type() const;
00617    SquareMatrix(const SquareMatrix& gm) : Matrix() { GetMatrix(&gm); }
00618    SquareMatrix(const Matrix& gm);
00619    void ReSize(int);                            // change dimensions
00620    virtual void ReSize(int,int);                // change dimensions
00621       // virtual so we will catch it being used in a vector called as a matrix
00622    void ReSize(const GeneralMatrix& A);
00623    void operator+=(const Matrix& M) { PlusEqual(M); }
00624    void operator-=(const Matrix& M) { MinusEqual(M); }
00625    void operator+=(Real f) { GeneralMatrix::Add(f); }
00626    void operator-=(Real f) { GeneralMatrix::Add(-f); }
00627    void swap(SquareMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
00628    NEW_DELETE(SquareMatrix)
00629 };
00630 
00631 class nricMatrix : public Matrix                // for use with Numerical
00632                                                 // Recipes in C
00633 {
00634    GeneralMatrix* Image() const;                // copy of matrix
00635    Real** row_pointer;                          // points to rows
00636    void MakeRowPointer();                       // build rowpointer
00637    void DeleteRowPointer();
00638 public:
00639    nricMatrix() : row_pointer(NULL) {}
00640    nricMatrix(int m, int n)                     // standard declaration
00641      :  Matrix(m,n), row_pointer(NULL) { MakeRowPointer(); }
00642    nricMatrix(const BaseMatrix& bm)             // evaluate BaseMatrix
00643      :  Matrix(bm), row_pointer(NULL) { MakeRowPointer(); }
00644    nricMatrix& operator=(const BaseMatrix& bm)
00645      { DeleteRowPointer(); Matrix::operator=(bm); MakeRowPointer(); return *this; }
00646    nricMatrix& operator=(Real f) { GeneralMatrix::operator=(f); return *this; }
00647    nricMatrix& operator=(const nricMatrix& m)
00648      { DeleteRowPointer(); Eq(m); MakeRowPointer(); return *this; }
00649    void operator<<(const BaseMatrix& X)
00650       { DeleteRowPointer(); Eq(X,this->Type(),true); MakeRowPointer(); }
00651    nricMatrix(const nricMatrix& gm) : Matrix(), row_pointer(NULL) { GetMatrix(&gm); MakeRowPointer(); }
00652    void ReSize(int m, int n)               // change dimensions
00653       { DeleteRowPointer(); Matrix::ReSize(m,n); MakeRowPointer(); }
00654    void ReSize(const GeneralMatrix& A);
00655    ~nricMatrix() { DeleteRowPointer(); }
00656    Real** nric() const { CheckStore(); return row_pointer-1; }
00657    void CleanUp();                                // to clear store
00658    void MiniCleanUp();
00659    void operator+=(const Matrix& M) { PlusEqual(M); }
00660    void operator-=(const Matrix& M) { MinusEqual(M); }
00661    void operator+=(Real f) { GeneralMatrix::Add(f); }
00662    void operator-=(Real f) { GeneralMatrix::Add(-f); }
00663    void swap(nricMatrix& gm);
00664    NEW_DELETE(nricMatrix)
00665 };
00666 
00667 class SymmetricMatrix : public GeneralMatrix
00668 {
00669    GeneralMatrix* Image() const;                // copy of matrix
00670 public:
00671    SymmetricMatrix() {}
00672    ~SymmetricMatrix() {}
00673    SymmetricMatrix(ArrayLengthSpecifier);
00674    SymmetricMatrix(const BaseMatrix&);
00675    void operator=(const BaseMatrix&);
00676    SymmetricMatrix& operator=(Real f) { GeneralMatrix::operator=(f); return *this; }
00677    SymmetricMatrix& operator=(const SymmetricMatrix& m) { Eq(m); return *this; }
00678    Real& operator()(int, int);                  // access element
00679    Real& element(int, int);                     // access element
00680    Real operator()(int, int) const;             // access element
00681    Real element(int, int) const;                // access element
00682 #ifdef SETUP_C_SUBSCRIPTS
00683    Real* operator[](int m) { return store+(m*(m+1))/2; }
00684    const Real* operator[](int m) const { return store+(m*(m+1))/2; }
00685 #endif
00686    MatrixType Type() const;
00687    SymmetricMatrix(const SymmetricMatrix& gm) : GeneralMatrix() { GetMatrix(&gm); }
00688    Real SumSquare() const;
00689    Real SumAbsoluteValue() const;
00690    Real Sum() const;
00691    Real Trace() const;
00692    void GetRow(MatrixRowCol&);
00693    void GetCol(MatrixRowCol&);
00694    void GetCol(MatrixColX&);
00695    void RestoreCol(MatrixRowCol&) {}
00696    void RestoreCol(MatrixColX&);
00697    GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
00698    void ReSize(int);                       // change dimensions
00699    void ReSize(const GeneralMatrix& A);
00700    void operator+=(const SymmetricMatrix& M) { PlusEqual(M); }
00701    void operator-=(const SymmetricMatrix& M) { MinusEqual(M); }
00702    void operator+=(Real f) { GeneralMatrix::Add(f); }
00703    void operator-=(Real f) { GeneralMatrix::Add(-f); }
00704    void swap(SymmetricMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
00705    NEW_DELETE(SymmetricMatrix)
00706 };
00707 
00708 class UpperTriangularMatrix : public GeneralMatrix
00709 {
00710    GeneralMatrix* Image() const;                // copy of matrix
00711 public:
00712    UpperTriangularMatrix() {}
00713    ~UpperTriangularMatrix() {}
00714    UpperTriangularMatrix(ArrayLengthSpecifier);
00715    void operator=(const BaseMatrix&);
00716    UpperTriangularMatrix& operator=(const UpperTriangularMatrix& m) { Eq(m); return *this; }
00717    UpperTriangularMatrix(const BaseMatrix&);
00718    UpperTriangularMatrix(const UpperTriangularMatrix& gm) : GeneralMatrix() { GetMatrix(&gm); }
00719    UpperTriangularMatrix& operator=(Real f) { GeneralMatrix::operator=(f); return *this; }
00720    Real& operator()(int, int);                  // access element
00721    Real& element(int, int);                     // access element
00722    Real operator()(int, int) const;             // access element
00723    Real element(int, int) const;                // access element
00724 #ifdef SETUP_C_SUBSCRIPTS
00725    Real* operator[](int m) { return store+m*ncols_value-(m*(m+1))/2; }
00726    const Real* operator[](int m) const
00727       { return store+m*ncols_value-(m*(m+1))/2; }
00728 #endif
00729    MatrixType Type() const;
00730    GeneralMatrix* MakeSolver() { return this; } // for solving
00731    void Solver(MatrixColX&, const MatrixColX&);
00732    LogAndSign LogDeterminant() const;
00733    Real Trace() const;
00734    void GetRow(MatrixRowCol&);
00735    void GetCol(MatrixRowCol&);
00736    void GetCol(MatrixColX&);
00737    void RestoreCol(MatrixRowCol&);
00738    void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
00739    void NextRow(MatrixRowCol&);
00740    void ReSize(int);                       // change dimensions
00741    void ReSize(const GeneralMatrix& A);
00742    MatrixBandWidth BandWidth() const;
00743    void operator+=(const UpperTriangularMatrix& M) { PlusEqual(M); }
00744    void operator-=(const UpperTriangularMatrix& M) { MinusEqual(M); }
00745    void operator+=(Real f) { GeneralMatrix::operator+=(f); }
00746    void operator-=(Real f) { GeneralMatrix::operator-=(f); }
00747    void swap(UpperTriangularMatrix& gm)
00748       { GeneralMatrix::swap((GeneralMatrix&)gm); }
00749    NEW_DELETE(UpperTriangularMatrix)
00750 };
00751 
00752 class LowerTriangularMatrix : public GeneralMatrix
00753 {
00754    GeneralMatrix* Image() const;                // copy of matrix
00755 public:
00756    LowerTriangularMatrix() {}
00757    ~LowerTriangularMatrix() {}
00758    LowerTriangularMatrix(ArrayLengthSpecifier);
00759    LowerTriangularMatrix(const LowerTriangularMatrix& gm) : GeneralMatrix() { GetMatrix(&gm); }
00760    LowerTriangularMatrix(const BaseMatrix& M);
00761    void operator=(const BaseMatrix&);
00762    LowerTriangularMatrix&  operator=(Real f) { GeneralMatrix::operator=(f); return *this; }
00763    LowerTriangularMatrix& operator=(const LowerTriangularMatrix& m) { Eq(m); return *this; }
00764    Real& operator()(int, int);                  // access element
00765    Real& element(int, int);                     // access element
00766    Real operator()(int, int) const;             // access element
00767    Real element(int, int) const;                // access element
00768 #ifdef SETUP_C_SUBSCRIPTS
00769    Real* operator[](int m) { return store+(m*(m+1))/2; }
00770    const Real* operator[](int m) const { return store+(m*(m+1))/2; }
00771 #endif
00772    MatrixType Type() const;
00773    GeneralMatrix* MakeSolver() { return this; } // for solving
00774    void Solver(MatrixColX&, const MatrixColX&);
00775    LogAndSign LogDeterminant() const;
00776    Real Trace() const;
00777    void GetRow(MatrixRowCol&);
00778    void GetCol(MatrixRowCol&);
00779    void GetCol(MatrixColX&);
00780    void RestoreCol(MatrixRowCol&);
00781    void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
00782    void NextRow(MatrixRowCol&);
00783    void ReSize(int);                       // change dimensions
00784    void ReSize(const GeneralMatrix& A);
00785    MatrixBandWidth BandWidth() const;
00786    void operator+=(const LowerTriangularMatrix& M) { PlusEqual(M); }
00787    void operator-=(const LowerTriangularMatrix& M) { MinusEqual(M); }
00788    void operator+=(Real f) { GeneralMatrix::operator+=(f); }
00789    void operator-=(Real f) { GeneralMatrix::operator-=(f); }
00790    void swap(LowerTriangularMatrix& gm)
00791       { GeneralMatrix::swap((GeneralMatrix&)gm); }
00792    NEW_DELETE(LowerTriangularMatrix)
00793 };
00794 
00795 class DiagonalMatrix : public GeneralMatrix
00796 {
00797    GeneralMatrix* Image() const;                // copy of matrix
00798 public:
00799    DiagonalMatrix() {}
00800    ~DiagonalMatrix() {}
00801    DiagonalMatrix(ArrayLengthSpecifier);
00802    DiagonalMatrix(const BaseMatrix&);
00803    DiagonalMatrix(const DiagonalMatrix& gm) : GeneralMatrix() { GetMatrix(&gm); }
00804    void operator=(const BaseMatrix&);
00805    DiagonalMatrix& operator=(Real f) { GeneralMatrix::operator=(f); return *this; }
00806    DiagonalMatrix& operator=(const DiagonalMatrix& m) { Eq(m); return *this; }
00807    Real& operator()(int, int);                  // access element
00808    Real& operator()(int);                       // access element
00809    Real operator()(int, int) const;             // access element
00810    Real operator()(int) const;
00811    Real& element(int, int);                     // access element
00812    Real& element(int);                          // access element
00813    Real element(int, int) const;                // access element
00814    Real element(int) const;                     // access element
00815 #ifdef SETUP_C_SUBSCRIPTS
00816    Real& operator[](int m) { return store[m]; }
00817    const Real& operator[](int m) const { return store[m]; }
00818 #endif
00819    MatrixType Type() const;
00820 
00821    LogAndSign LogDeterminant() const;
00822    Real Trace() const;
00823    void GetRow(MatrixRowCol&);
00824    void GetCol(MatrixRowCol&);
00825    void GetCol(MatrixColX&);
00826    void NextRow(MatrixRowCol&);
00827    void NextCol(MatrixRowCol&);
00828    void NextCol(MatrixColX&);
00829    GeneralMatrix* MakeSolver() { return this; } // for solving
00830    void Solver(MatrixColX&, const MatrixColX&);
00831    GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
00832    void ReSize(int);                       // change dimensions
00833    void ReSize(const GeneralMatrix& A);
00834    Real* nric() const
00835       { CheckStore(); return store-1; }         // for use by NRIC
00836    MatrixBandWidth BandWidth() const;
00837 //   ReturnMatrix Reverse() const;                // reverse order of elements
00838    void operator+=(const DiagonalMatrix& M) { PlusEqual(M); }
00839    void operator-=(const DiagonalMatrix& M) { MinusEqual(M); }
00840    void operator+=(Real f) { GeneralMatrix::operator+=(f); }
00841    void operator-=(Real f) { GeneralMatrix::operator-=(f); }
00842    void swap(DiagonalMatrix& gm)
00843       { GeneralMatrix::swap((GeneralMatrix&)gm); }
00844    NEW_DELETE(DiagonalMatrix)
00845 };
00846 
00847 class RowVector : public Matrix
00848 {
00849    GeneralMatrix* Image() const;                // copy of matrix
00850 public:
00851    RowVector() { nrows_value = 1; }
00852    ~RowVector() {}
00853    RowVector(ArrayLengthSpecifier n) : Matrix(1,n.Value()) {}
00854    RowVector(const BaseMatrix&);
00855    RowVector(const RowVector& gm) : Matrix() { GetMatrix(&gm); }
00856    void operator=(const BaseMatrix&);
00857    RowVector& operator=(Real f) { GeneralMatrix::operator=(f); return *this; }
00858    RowVector& operator=(const RowVector& m) { Eq(m); return *this; }
00859    Real& operator()(int);                       // access element
00860    Real& element(int);                          // access element
00861    Real operator()(int) const;                  // access element
00862    Real element(int) const;                     // access element
00863 #ifdef SETUP_C_SUBSCRIPTS
00864    Real& operator[](int m) { return store[m]; }
00865    const Real& operator[](int m) const { return store[m]; }
00866    // following for Numerical Recipes in C++
00867    RowVector(Real a, int n) : Matrix(a, 1, n) {}
00868    RowVector(const Real* a, int n) : Matrix(a, 1, n) {}
00869 #endif
00870    MatrixType Type() const;
00871    void GetCol(MatrixRowCol&);
00872    void GetCol(MatrixColX&);
00873    void NextCol(MatrixRowCol&);
00874    void NextCol(MatrixColX&);
00875    void RestoreCol(MatrixRowCol&) {}
00876    void RestoreCol(MatrixColX& c);
00877    GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
00878    void ReSize(int);                       // change dimensions
00879    void ReSize(int,int);                   // in case access is matrix
00880    void ReSize(const GeneralMatrix& A);
00881    Real* nric() const
00882       { CheckStore(); return store-1; }         // for use by NRIC
00883    void CleanUp();                              // to clear store
00884    void MiniCleanUp()
00885       { store = 0; storage = 0; nrows_value = 1; ncols_value = 0; tag = -1; }
00886    // friend ReturnMatrix GetMatrixRow(Matrix& A, int row);
00887    void operator+=(const Matrix& M) { PlusEqual(M); }
00888    void operator-=(const Matrix& M) { MinusEqual(M); }
00889    void operator+=(Real f) { GeneralMatrix::Add(f); }
00890    void operator-=(Real f) { GeneralMatrix::Add(-f); }
00891    void swap(RowVector& gm)
00892       { GeneralMatrix::swap((GeneralMatrix&)gm); }
00893    NEW_DELETE(RowVector)
00894 };
00895 
00896 class ColumnVector : public Matrix
00897 {
00898    GeneralMatrix* Image() const;                // copy of matrix
00899 public:
00900    ColumnVector() { ncols_value = 1; }
00901    ~ColumnVector() {}
00902    ColumnVector(ArrayLengthSpecifier n) : Matrix(n.Value(),1) {}
00903    ColumnVector(const BaseMatrix&);
00904    ColumnVector(const ColumnVector& gm) : Matrix() { GetMatrix(&gm); }
00905    void operator=(const BaseMatrix&);
00906    ColumnVector& operator=(Real f) { GeneralMatrix::operator=(f); return *this; }
00907    ColumnVector& operator=(const ColumnVector& m) { Eq(m); return *this; }
00908    Real& operator()(int);                       // access element
00909    Real& element(int);                          // access element
00910    Real operator()(int) const;                  // access element
00911    Real element(int) const;                     // access element
00912 #ifdef SETUP_C_SUBSCRIPTS
00913    Real& operator[](int m) { return store[m]; }
00914    const Real& operator[](int m) const { return store[m]; }
00915    // following for Numerical Recipes in C++
00916    ColumnVector(Real a, int m) : Matrix(a, m, 1) {}
00917    ColumnVector(const Real* a, int m) : Matrix(a, m, 1) {}
00918 #endif
00919    MatrixType Type() const;
00920    GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
00921    void ReSize(int);                       // change dimensions
00922    void ReSize(int,int);                   // in case access is matrix
00923    void ReSize(const GeneralMatrix& A);
00924    Real* nric() const
00925       { CheckStore(); return store-1; }         // for use by NRIC
00926    void CleanUp();                              // to clear store
00927    void MiniCleanUp()
00928       { store = 0; storage = 0; nrows_value = 0; ncols_value = 1; tag = -1; }
00929 //   ReturnMatrix Reverse() const;                // reverse order of elements
00930    void operator+=(const Matrix& M) { PlusEqual(M); }
00931    void operator-=(const Matrix& M) { MinusEqual(M); }
00932    void operator+=(Real f) { GeneralMatrix::Add(f); }
00933    void operator-=(Real f) { GeneralMatrix::Add(-f); }
00934    void swap(ColumnVector& gm)
00935       { GeneralMatrix::swap((GeneralMatrix&)gm); }
00936    NEW_DELETE(ColumnVector)
00937 };
00938 
00939 class CroutMatrix : public GeneralMatrix        // for LU decomposition
00940 {
00941    int* indx;
00942    bool d;
00943    bool sing;
00944    void ludcmp();
00945    void operator=(const CroutMatrix& /*m*/);     // not allowed
00946 public:
00947    CroutMatrix(const BaseMatrix&);
00948    CroutMatrix(const CroutMatrix&);
00949    MatrixType Type() const;
00950    void lubksb(Real*, int=0);
00951    ~CroutMatrix();
00952    GeneralMatrix* MakeSolver() { return this; } // for solving
00953    LogAndSign LogDeterminant() const;
00954    void Solver(MatrixColX&, const MatrixColX&);
00955    void GetRow(MatrixRowCol&);
00956    void GetCol(MatrixRowCol&);
00957    void GetCol(MatrixColX& c) { GetCol((MatrixRowCol&)c); }
00958    void CleanUp();                                // to clear store
00959    void MiniCleanUp();
00960    bool IsEqual(const GeneralMatrix&) const;
00961    bool IsSingular() const { return sing; }
00962    void swap(CroutMatrix& gm);
00963    NEW_DELETE(CroutMatrix)
00964 };
00965 
00966 // ***************************** band matrices ***************************/
00967 
00968 class BandMatrix : public GeneralMatrix         // band matrix
00969 {
00970    GeneralMatrix* Image() const;                // copy of matrix
00971 protected:
00972    void CornerClear() const;                    // set unused elements to zero
00973    short SimpleAddOK(const GeneralMatrix* gm);
00974 public:
00975    int lower, upper;                            // band widths
00976    BandMatrix() : lower(0), upper(0) { CornerClear(); }
00977    ~BandMatrix() {}
00978    BandMatrix(int n,int lb,int ub) : lower(0), upper(0) { ReSize(n,lb,ub); CornerClear(); }
00979                                                 // standard declaration
00980    BandMatrix(const BaseMatrix&);               // evaluate BaseMatrix
00981    void operator=(const BaseMatrix&);
00982    BandMatrix& operator=(Real f) { GeneralMatrix::operator=(f); return *this; }
00983    BandMatrix& operator=(const BandMatrix& m) { Eq(m); return *this; }
00984    MatrixType Type() const;
00985    Real& operator()(int, int);                  // access element
00986    Real& element(int, int);                     // access element
00987    Real operator()(int, int) const;             // access element
00988    Real element(int, int) const;                // access element
00989 #ifdef SETUP_C_SUBSCRIPTS
00990    Real* operator[](int m) { return store+(upper+lower)*m+lower; }
00991    const Real* operator[](int m) const { return store+(upper+lower)*m+lower; }
00992 #endif
00993    BandMatrix(const BandMatrix& gm) : GeneralMatrix(), lower(0), upper(0) { GetMatrix(&gm); }
00994    LogAndSign LogDeterminant() const;
00995    GeneralMatrix* MakeSolver();
00996    Real Trace() const;
00997    Real SumSquare() const { CornerClear(); return GeneralMatrix::SumSquare(); }
00998    Real SumAbsoluteValue() const
00999       { CornerClear(); return GeneralMatrix::SumAbsoluteValue(); }
01000    Real Sum() const
01001       { CornerClear(); return GeneralMatrix::Sum(); }
01002    Real MaximumAbsoluteValue() const
01003       { CornerClear(); return GeneralMatrix::MaximumAbsoluteValue(); }
01004    Real MinimumAbsoluteValue() const
01005       { int j, k; return GeneralMatrix::MinimumAbsoluteValue2(j, k); }
01006    Real Maximum() const { int j, k; return GeneralMatrix::Maximum2(j, k); }
01007    Real Minimum() const { int j, k; return GeneralMatrix::Minimum2(j, k); }
01008    void GetRow(MatrixRowCol&);
01009    void GetCol(MatrixRowCol&);
01010    void GetCol(MatrixColX&);
01011    void RestoreCol(MatrixRowCol&);
01012    void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
01013    void NextRow(MatrixRowCol&);
01014    virtual void ReSize(int, int, int);             // change dimensions
01015    void ReSize(const GeneralMatrix& A);
01016    bool SameStorageType(const GeneralMatrix& A) const;
01017    void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
01018    void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
01019    MatrixBandWidth BandWidth() const;
01020    void SetParameters(const GeneralMatrix*);
01021    MatrixInput operator<<(Real);                // will give error
01022    MatrixInput operator<<(int f);
01023    void operator<<(const Real* r);              // will give error
01024    void operator<<(const int* r);               // will give error
01025       // the next is included because Zortech and Borland
01026       // cannot find the copy in GeneralMatrix
01027    void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
01028    void swap(BandMatrix& gm);
01029    NEW_DELETE(BandMatrix)
01030 };
01031 
01032 class UpperBandMatrix : public BandMatrix       // upper band matrix
01033 {
01034    GeneralMatrix* Image() const;                // copy of matrix
01035 public:
01036    UpperBandMatrix() {}
01037    ~UpperBandMatrix() {}
01038    UpperBandMatrix(int n, int ubw)              // standard declaration
01039       : BandMatrix(n, 0, ubw) {}
01040    UpperBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
01041    void operator=(const BaseMatrix&);
01042    UpperBandMatrix& operator=(Real f) { GeneralMatrix::operator=(f); return *this;}
01043    UpperBandMatrix& operator=(const UpperBandMatrix& m) { Eq(m); return *this;}
01044    MatrixType Type() const;
01045    UpperBandMatrix(const UpperBandMatrix& gm) : BandMatrix() { GetMatrix(&gm); }
01046    GeneralMatrix* MakeSolver() { return this; }
01047    void Solver(MatrixColX&, const MatrixColX&);
01048    LogAndSign LogDeterminant() const;
01049    void ReSize(int, int, int);             // change dimensions
01050    void ReSize(int n,int ubw)              // change dimensions
01051       { BandMatrix::ReSize(n,0,ubw); }
01052    void ReSize(const GeneralMatrix& A) { BandMatrix::ReSize(A); }
01053    Real& operator()(int, int);
01054    Real operator()(int, int) const;
01055    Real& element(int, int);
01056    Real element(int, int) const;
01057 #ifdef SETUP_C_SUBSCRIPTS
01058    Real* operator[](int m) { return store+upper*m; }
01059    const Real* operator[](int m) const { return store+upper*m; }
01060 #endif
01061    void swap(UpperBandMatrix& gm)
01062       { BandMatrix::swap((BandMatrix&)gm); }
01063    NEW_DELETE(UpperBandMatrix)
01064 };
01065 
01066 class LowerBandMatrix : public BandMatrix       // upper band matrix
01067 {
01068    GeneralMatrix* Image() const;                // copy of matrix
01069 public:
01070    LowerBandMatrix() {}
01071    ~LowerBandMatrix() {}
01072    LowerBandMatrix(int n, int lbw)              // standard declaration
01073       : BandMatrix(n, lbw, 0) {}
01074    LowerBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
01075    void operator=(const BaseMatrix&);
01076    LowerBandMatrix& operator=(Real f) { GeneralMatrix::operator=(f); return *this;}
01077    LowerBandMatrix& operator=(const LowerBandMatrix& m) { Eq(m); return *this;}
01078    MatrixType Type() const;
01079    LowerBandMatrix(const LowerBandMatrix& gm) : BandMatrix() { GetMatrix(&gm); }
01080    GeneralMatrix* MakeSolver() { return this; }
01081    void Solver(MatrixColX&, const MatrixColX&);
01082    LogAndSign LogDeterminant() const;
01083    void ReSize(int, int, int);             // change dimensions
01084    void ReSize(int n,int lbw)             // change dimensions
01085       { BandMatrix::ReSize(n,lbw,0); }
01086    void ReSize(const GeneralMatrix& A) { BandMatrix::ReSize(A); }
01087    Real& operator()(int, int);
01088    Real operator()(int, int) const;
01089    Real& element(int, int);
01090    Real element(int, int) const;
01091 #ifdef SETUP_C_SUBSCRIPTS
01092    Real* operator[](int m) { return store+lower*(m+1); }
01093    const Real* operator[](int m) const { return store+lower*(m+1); }
01094 #endif
01095    void swap(LowerBandMatrix& gm)
01096       { BandMatrix::swap((BandMatrix&)gm); }
01097    NEW_DELETE(LowerBandMatrix)
01098 };
01099 
01100 class SymmetricBandMatrix : public GeneralMatrix
01101 {
01102    GeneralMatrix* Image() const;                // copy of matrix
01103    void CornerClear() const;                    // set unused elements to zero
01104    short SimpleAddOK(const GeneralMatrix* gm);
01105 public:
01106    int lower;                                   // lower band width
01107    SymmetricBandMatrix() : lower(0) { CornerClear(); }
01108    ~SymmetricBandMatrix() {}
01109    SymmetricBandMatrix(int n, int lb) : lower(0) { ReSize(n,lb); CornerClear(); }
01110    SymmetricBandMatrix(const BaseMatrix&);
01111    void operator=(const BaseMatrix&);
01112    SymmetricBandMatrix& operator=(Real f) { GeneralMatrix::operator=(f); return *this; }
01113    SymmetricBandMatrix& operator=(const SymmetricBandMatrix& m) { Eq(m); return *this; }
01114    Real& operator()(int, int);                  // access element
01115    Real& element(int, int);                     // access element
01116    Real operator()(int, int) const;             // access element
01117    Real element(int, int) const;                // access element
01118 #ifdef SETUP_C_SUBSCRIPTS
01119    Real* operator[](int m) { return store+lower*(m+1); }
01120    const Real* operator[](int m) const { return store+lower*(m+1); }
01121 #endif
01122    MatrixType Type() const;
01123    SymmetricBandMatrix(const SymmetricBandMatrix& gm) : GeneralMatrix(), lower(0) { GetMatrix(&gm); }
01124    GeneralMatrix* MakeSolver();
01125    Real SumSquare() const;
01126    Real SumAbsoluteValue() const;
01127    Real Sum() const;
01128    Real MaximumAbsoluteValue() const
01129       { CornerClear(); return GeneralMatrix::MaximumAbsoluteValue(); }
01130    Real MinimumAbsoluteValue() const
01131       { int j, k; return GeneralMatrix::MinimumAbsoluteValue2(j, k); }
01132    Real Maximum() const { int j, k; return GeneralMatrix::Maximum2(j, k); }
01133    Real Minimum() const { int j, k; return GeneralMatrix::Minimum2(j, k); }
01134    Real Trace() const;
01135    LogAndSign LogDeterminant() const;
01136    void GetRow(MatrixRowCol&);
01137    void GetCol(MatrixRowCol&);
01138    void GetCol(MatrixColX&);
01139    void RestoreCol(MatrixRowCol&) {}
01140    void RestoreCol(MatrixColX&);
01141    GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
01142    void ReSize(int,int);                       // change dimensions
01143    void ReSize(const GeneralMatrix& A);
01144    bool SameStorageType(const GeneralMatrix& A) const;
01145    void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
01146    void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
01147    MatrixBandWidth BandWidth() const;
01148    void SetParameters(const GeneralMatrix*);
01149    void operator<<(const Real* r);              // will give error
01150    void operator<<(const int* r);               // will give error
01151    void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
01152    void swap(SymmetricBandMatrix& gm);
01153    NEW_DELETE(SymmetricBandMatrix)
01154 };
01155 
01156 class BandLUMatrix : public GeneralMatrix
01157 // for LU decomposition of band matrix
01158 {
01159    int* indx;
01160    bool d;
01161    bool sing;                                   // true if singular
01162    Real* store2;
01163    int storage2;
01164    void ludcmp();
01165    int m1,m2;                                   // lower and upper
01166    void operator=(const BandLUMatrix& /*m*/);   // no allowed
01167 public:
01168    BandLUMatrix(const BaseMatrix&);
01169    BandLUMatrix(const BandLUMatrix&);
01170    MatrixType Type() const;
01171    void lubksb(Real*, int=0);
01172    ~BandLUMatrix();
01173    GeneralMatrix* MakeSolver() { return this; } // for solving
01174    LogAndSign LogDeterminant() const;
01175    void Solver(MatrixColX&, const MatrixColX&);
01176    void GetRow(MatrixRowCol&);
01177    void GetCol(MatrixRowCol&);
01178    void GetCol(MatrixColX& c) { GetCol((MatrixRowCol&)c); }
01179    void CleanUp();                                // to clear store
01180    void MiniCleanUp();
01181    bool IsEqual(const GeneralMatrix&) const;
01182    bool IsSingular() const { return sing; }
01183    void swap(BandLUMatrix& gm);
01184    NEW_DELETE(BandLUMatrix)
01185 };
01186 
01187 // ************************** special matrices ****************************
01188 
01189 class IdentityMatrix : public GeneralMatrix
01190 {
01191    GeneralMatrix* Image() const;          // copy of matrix
01192 public:
01193    IdentityMatrix() {}
01194    ~IdentityMatrix() {}
01195    IdentityMatrix(ArrayLengthSpecifier n) : GeneralMatrix(1)
01196       { nrows_value = ncols_value = n.Value(); *store = 1; }
01197    IdentityMatrix(const IdentityMatrix& gm) : GeneralMatrix() { GetMatrix(&gm); }
01198    IdentityMatrix(const BaseMatrix&);
01199    void operator=(const BaseMatrix&);
01200    IdentityMatrix& operator=(const IdentityMatrix& m) { Eq(m); return *this; }
01201    IdentityMatrix& operator=(Real f) { GeneralMatrix::operator=(f); return *this; }
01202    MatrixType Type() const;
01203 
01204    LogAndSign LogDeterminant() const;
01205    Real Trace() const;
01206    Real SumSquare() const;
01207    Real SumAbsoluteValue() const;
01208    Real Sum() const { return Trace(); }
01209    void GetRow(MatrixRowCol&);
01210    void GetCol(MatrixRowCol&);
01211    void GetCol(MatrixColX&);
01212    void NextRow(MatrixRowCol&);
01213    void NextCol(MatrixRowCol&);
01214    void NextCol(MatrixColX&);
01215    GeneralMatrix* MakeSolver() { return this; } // for solving
01216    void Solver(MatrixColX&, const MatrixColX&);
01217    GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
01218    void ReSize(int n);
01219    void ReSize(const GeneralMatrix& A);
01220    MatrixBandWidth BandWidth() const;
01221 //   ReturnMatrix Reverse() const;                // reverse order of elements
01222    void swap(IdentityMatrix& gm)
01223       { GeneralMatrix::swap((GeneralMatrix&)gm); }
01224    NEW_DELETE(IdentityMatrix)
01225 };
01226 
01227 
01228 
01229 
01230 // ************************** GenericMatrix class ************************/
01231 
01232 class GenericMatrix : public BaseMatrix
01233 {
01234    GeneralMatrix* gm;
01235    int search(const BaseMatrix* bm) const;
01236    friend class BaseMatrix;
01237 public:
01238    GenericMatrix() : gm(0) {}
01239    GenericMatrix(const BaseMatrix& bm) : gm((const_cast<BaseMatrix&>(bm)).Evaluate()) { gm = gm->Image(); }
01240    GenericMatrix(const GenericMatrix& bm) : BaseMatrix(), gm(bm.gm->Image()) {}
01241    void operator=(const GenericMatrix&);
01242    void operator=(const BaseMatrix&);
01243    void operator+=(const BaseMatrix&);
01244    void operator-=(const BaseMatrix&);
01245    void operator*=(const BaseMatrix&);
01246    void operator|=(const BaseMatrix&);
01247    void operator&=(const BaseMatrix&);
01248    void operator+=(Real);
01249    GenericMatrix& operator-=(Real r) { operator+=(-r); return *this; }
01250    void operator*=(Real);
01251    GenericMatrix& operator/=(Real r) { operator*=(1/r); return *this; }
01252    ~GenericMatrix() { delete gm; }
01253    void CleanUp() { delete gm; gm = 0; }
01254    void Release() { gm->Release(); }
01255    GeneralMatrix* Evaluate(MatrixType = MatrixTypeUnSp);
01256    MatrixBandWidth BandWidth() const;
01257    void swap(GenericMatrix& x);
01258    NEW_DELETE(GenericMatrix)
01259 };
01260 
01261 // *************************** temporary classes *************************/
01262 
01263 class MultipliedMatrix : public BaseMatrix
01264 {
01265 protected:
01266    // if these union statements cause problems, simply remove them
01267    // and declare the items individually
01268    union { const BaseMatrix* bm1; GeneralMatrix* gm1; };
01269               // pointers to summands
01270    union { const BaseMatrix* bm2; GeneralMatrix* gm2; };
01271    MultipliedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
01272       : bm1(bm1x),bm2(bm2x) {}
01273    int search(const BaseMatrix*) const;
01274    friend class BaseMatrix;
01275    friend class GeneralMatrix;
01276    friend class GenericMatrix;
01277 public:
01278    MultipliedMatrix(const MultipliedMatrix& mm) // needed a copy constructor to silence gcc4 warning
01279       : BaseMatrix(mm), bm1(mm.bm1),bm2(mm.bm2) {}
01280    ~MultipliedMatrix() {}
01281    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01282    MatrixBandWidth BandWidth() const;
01283    NEW_DELETE(MultipliedMatrix)
01284 };
01285 
01286 class AddedMatrix : public MultipliedMatrix
01287 {
01288 protected:
01289    AddedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
01290       : MultipliedMatrix(bm1x,bm2x) {}
01291 
01292    friend class BaseMatrix;
01293    friend class GeneralMatrix;
01294    friend class GenericMatrix;
01295 public:
01296    ~AddedMatrix() {}
01297    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01298    MatrixBandWidth BandWidth() const;
01299    NEW_DELETE(AddedMatrix)
01300 };
01301 
01302 class SPMatrix : public AddedMatrix
01303 {
01304 protected:
01305    SPMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
01306       : AddedMatrix(bm1x,bm2x) {}
01307 
01308    friend class BaseMatrix;
01309    friend class GeneralMatrix;
01310    friend class GenericMatrix;
01311 public:
01312    ~SPMatrix() {}
01313    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01314    MatrixBandWidth BandWidth() const;
01315 
01316    friend SPMatrix SP(const BaseMatrix&, const BaseMatrix&);
01317 
01318    NEW_DELETE(SPMatrix)
01319 };
01320 
01321 class KPMatrix : public MultipliedMatrix
01322 {
01323 protected:
01324    KPMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
01325       : MultipliedMatrix(bm1x,bm2x) {}
01326 
01327    friend class BaseMatrix;
01328    friend class GeneralMatrix;
01329    friend class GenericMatrix;
01330 public:
01331    ~KPMatrix() {}
01332    MatrixBandWidth BandWidth() const;
01333    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01334    friend KPMatrix KP(const BaseMatrix&, const BaseMatrix&);
01335    NEW_DELETE(KPMatrix)
01336 };
01337 
01338 class ConcatenatedMatrix : public MultipliedMatrix
01339 {
01340 protected:
01341    ConcatenatedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
01342       : MultipliedMatrix(bm1x,bm2x) {}
01343 
01344    friend class BaseMatrix;
01345    friend class GeneralMatrix;
01346    friend class GenericMatrix;
01347 public:
01348    ~ConcatenatedMatrix() {}
01349    MatrixBandWidth BandWidth() const;
01350    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01351    NEW_DELETE(ConcatenatedMatrix)
01352 };
01353 
01354 class StackedMatrix : public ConcatenatedMatrix
01355 {
01356 protected:
01357    StackedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
01358       : ConcatenatedMatrix(bm1x,bm2x) {}
01359 
01360    friend class BaseMatrix;
01361    friend class GeneralMatrix;
01362    friend class GenericMatrix;
01363 public:
01364    ~StackedMatrix() {}
01365    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01366    NEW_DELETE(StackedMatrix)
01367 };
01368 
01369 class SolvedMatrix : public MultipliedMatrix
01370 {
01371    SolvedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
01372       : MultipliedMatrix(bm1x,bm2x) {}
01373    friend class BaseMatrix;
01374    friend class InvertedMatrix;                        // for operator*
01375 public:
01376    ~SolvedMatrix() {}
01377    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01378    MatrixBandWidth BandWidth() const;
01379    NEW_DELETE(SolvedMatrix)
01380 };
01381 
01382 class SubtractedMatrix : public AddedMatrix
01383 {
01384    SubtractedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
01385       : AddedMatrix(bm1x,bm2x) {}
01386    friend class BaseMatrix;
01387    friend class GeneralMatrix;
01388    friend class GenericMatrix;
01389 public:
01390    ~SubtractedMatrix() {}
01391    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01392    NEW_DELETE(SubtractedMatrix)
01393 };
01394 
01395 class ShiftedMatrix : public BaseMatrix
01396 {
01397 protected:
01398    union { const BaseMatrix* bm; GeneralMatrix* gm; };
01399    Real f;
01400    ShiftedMatrix(const BaseMatrix* bmx, Real fx) : bm(bmx),f(fx) {}
01401    int search(const BaseMatrix*) const;
01402    friend class BaseMatrix;
01403    friend class GeneralMatrix;
01404    friend class GenericMatrix;
01405 public:
01406    ShiftedMatrix(const ShiftedMatrix& im) : BaseMatrix(im), bm(im.bm), f(im.f) {}
01407    ~ShiftedMatrix() {}
01408    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01409    friend ShiftedMatrix operator+(Real f, const BaseMatrix& BM);
01410    NEW_DELETE(ShiftedMatrix)
01411 };
01412 
01413 class NegShiftedMatrix : public ShiftedMatrix
01414 {
01415 protected:
01416    NegShiftedMatrix(Real fx, const BaseMatrix* bmx) : ShiftedMatrix(bmx,fx) {}
01417    friend class BaseMatrix;
01418    friend class GeneralMatrix;
01419    friend class GenericMatrix;
01420 public:
01421    ~NegShiftedMatrix() {}
01422    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01423    friend NegShiftedMatrix operator-(Real, const BaseMatrix&);
01424    NEW_DELETE(NegShiftedMatrix)
01425 };
01426 
01427 class ScaledMatrix : public ShiftedMatrix
01428 {
01429    ScaledMatrix(const BaseMatrix* bmx, Real fx) : ShiftedMatrix(bmx,fx) {}
01430    friend class BaseMatrix;
01431    friend class GeneralMatrix;
01432    friend class GenericMatrix;
01433 public:
01434    ScaledMatrix(const ScaledMatrix& im) : ShiftedMatrix(im) {}
01435    ~ScaledMatrix() {}
01436    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01437    MatrixBandWidth BandWidth() const;
01438    friend ScaledMatrix operator*(Real f, const BaseMatrix& BM);
01439    NEW_DELETE(ScaledMatrix)
01440 };
01441 
01442 class NegatedMatrix : public BaseMatrix
01443 {
01444 protected:
01445    union { const BaseMatrix* bm; GeneralMatrix* gm; };
01446    NegatedMatrix(const BaseMatrix* bmx) : bm(bmx) {}
01447    int search(const BaseMatrix*) const;
01448 private:
01449    friend class BaseMatrix;
01450 public:
01451    NegatedMatrix(const NegatedMatrix& im) : BaseMatrix(im), bm(im.bm) {}
01452    ~NegatedMatrix() {}
01453    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01454    MatrixBandWidth BandWidth() const;
01455    NEW_DELETE(NegatedMatrix)
01456 };
01457 
01458 class TransposedMatrix : public NegatedMatrix
01459 {
01460    TransposedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
01461    friend class BaseMatrix;
01462 public:
01463    ~TransposedMatrix() {}
01464    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01465    MatrixBandWidth BandWidth() const;
01466    NEW_DELETE(TransposedMatrix)
01467 };
01468 
01469 class ReversedMatrix : public NegatedMatrix
01470 {
01471    ReversedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
01472    friend class BaseMatrix;
01473 public:
01474    ~ReversedMatrix() {}
01475    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01476    NEW_DELETE(ReversedMatrix)
01477 };
01478 
01479 class InvertedMatrix : public NegatedMatrix
01480 {
01481    InvertedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
01482 public:
01483    InvertedMatrix(const InvertedMatrix& im) : NegatedMatrix(im) {}
01484    ~InvertedMatrix() {}
01485    SolvedMatrix operator*(const BaseMatrix&) const;       // inverse(A) * B
01486    ScaledMatrix operator*(Real sc) const { return BaseMatrix::operator*(sc); }
01487    friend class BaseMatrix;
01488    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01489    MatrixBandWidth BandWidth() const;
01490    NEW_DELETE(InvertedMatrix)
01491 };
01492 
01493 class RowedMatrix : public NegatedMatrix
01494 {
01495    RowedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
01496    friend class BaseMatrix;
01497 public:
01498    ~RowedMatrix() {}
01499    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01500    MatrixBandWidth BandWidth() const;
01501    NEW_DELETE(RowedMatrix)
01502 };
01503 
01504 class ColedMatrix : public NegatedMatrix
01505 {
01506    ColedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
01507    friend class BaseMatrix;
01508 public:
01509    ~ColedMatrix() {}
01510    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01511    MatrixBandWidth BandWidth() const;
01512    NEW_DELETE(ColedMatrix)
01513 };
01514 
01515 class DiagedMatrix : public NegatedMatrix
01516 {
01517    DiagedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
01518    friend class BaseMatrix;
01519 public:
01520    ~DiagedMatrix() {}
01521    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01522    MatrixBandWidth BandWidth() const;
01523    NEW_DELETE(DiagedMatrix)
01524 };
01525 
01526 class MatedMatrix : public NegatedMatrix
01527 {
01528    int nr, nc;
01529    MatedMatrix(const BaseMatrix* bmx, int nrx, int ncx)
01530       : NegatedMatrix(bmx), nr(nrx), nc(ncx) {}
01531    friend class BaseMatrix;
01532 public:
01533    ~MatedMatrix() {}
01534    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01535    MatrixBandWidth BandWidth() const;
01536    NEW_DELETE(MatedMatrix)
01537 };
01538 
01539 class ReturnMatrix : public BaseMatrix    // for matrix return
01540 {
01541    GeneralMatrix* gm;
01542    int search(const BaseMatrix*) const;
01543 public:
01544    ~ReturnMatrix() {}
01545    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01546    friend class BaseMatrix;
01547    ReturnMatrix(const ReturnMatrix& tm) : BaseMatrix(), gm(tm.gm) {}
01548    ReturnMatrix(const GeneralMatrix* gmx) : gm(const_cast<GeneralMatrix*>(gmx)) {}
01549 //   ReturnMatrix(GeneralMatrix&);
01550    MatrixBandWidth BandWidth() const;
01551    void operator=(const ReturnMatrix&);
01552    NEW_DELETE(ReturnMatrix)
01553 };
01554 
01555 
01556 // ************************** submatrices ******************************/
01557 
01558 class GetSubMatrix : public NegatedMatrix
01559 {
01560    int row_skip;
01561    int row_number;
01562    int col_skip;
01563    int col_number;
01564    bool IsSym;
01565 
01566    GetSubMatrix
01567       (const BaseMatrix* bmx, int rs, int rn, int cs, int cn, bool is)
01568       : NegatedMatrix(bmx),
01569       row_skip(rs), row_number(rn), col_skip(cs), col_number(cn), IsSym(is) {}
01570    void SetUpLHS();
01571    friend class BaseMatrix;
01572 public:
01573    GetSubMatrix(const GetSubMatrix& g)
01574       : NegatedMatrix(g.bm), row_skip(g.row_skip), row_number(g.row_number),
01575       col_skip(g.col_skip), col_number(g.col_number), IsSym(g.IsSym) {}
01576    ~GetSubMatrix() {}
01577    GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp);
01578    void operator=(const BaseMatrix&);
01579    void operator+=(const BaseMatrix&);
01580    void operator-=(const BaseMatrix&);
01581    GetSubMatrix& operator=(const GetSubMatrix& m) { operator=((const BaseMatrix&)m); return *this; }
01582    void operator<<(const BaseMatrix&);
01583    void operator<<(const Real*);                // copy from array
01584    void operator<<(const int*);                 // copy from array
01585    MatrixInput operator<<(Real);                // for loading a list
01586    MatrixInput operator<<(int f);
01587    void operator=(Real);                        // copy from constant
01588    void operator+=(Real);                       // add constant
01589    void operator-=(Real r) { operator+=(-r); }  // subtract constant
01590    void operator*=(Real);                       // multiply by constant
01591    void operator/=(Real r) { operator*=(1/r); } // divide by constant
01592    void Inject(const GeneralMatrix&);           // copy stored els only
01593    MatrixBandWidth BandWidth() const;
01594    NEW_DELETE(GetSubMatrix)
01595 };
01596 
01597 // ******************** linear equation solving ****************************/
01598 
01599 class LinearEquationSolver : public BaseMatrix
01600 {
01601    GeneralMatrix* gm;
01602    int search(const BaseMatrix*) const { return 0; }
01603    friend class BaseMatrix;
01604 public:
01605    LinearEquationSolver(const BaseMatrix& bm);
01606    LinearEquationSolver(const LinearEquationSolver&);
01607    void operator=(const LinearEquationSolver&);
01608    ~LinearEquationSolver() { delete gm; }
01609    void CleanUp() { delete gm; } 
01610    GeneralMatrix* Evaluate(MatrixType) { return gm; }
01611    // probably should have an error message if MatrixType != UnSp
01612    NEW_DELETE(LinearEquationSolver)
01613 };
01614 
01615 // ************************** matrix input *******************************/
01616 
01617 class MatrixInput          // for reading a list of values into a matrix
01618                            // the difficult part is detecting a mismatch
01619                            // in the number of elements
01620 {
01621    int n;                  // number values still to be read
01622    Real* r;                // pointer to next location to be read to
01623 public:
01624    MatrixInput(const MatrixInput& mi) : n(mi.n), r(mi.r) {}
01625    MatrixInput(int nx, Real* rx) : n(nx), r(rx) {}
01626    ~MatrixInput();
01627    MatrixInput operator<<(Real);
01628    MatrixInput operator<<(int f);
01629    void operator=(const MatrixInput&);
01630    friend class GeneralMatrix;
01631 };
01632 
01633 
01634 
01635 // **************** a very simple integer array class ********************/
01636 
01637 // A minimal array class to imitate a C style array but giving dynamic storage
01638 // mostly intended for internal use by newmat
01639 
01640 class SimpleIntArray : public Janitor
01641 {
01642 protected:
01643    int* a;                    // pointer to the array
01644    int n;                     // length of the array
01645 public:
01646    SimpleIntArray(int xn);    // build an array length xn
01647    SimpleIntArray() : a(0), n(0) {}  // build an array length 0
01648    ~SimpleIntArray();         // return the space to memory
01649    int& operator[](int i);    // access element of the array - start at 0
01650    int operator[](int i) const;
01651             // access element of constant array
01652    void operator=(int ai);    // set the array equal to a constant
01653    void operator=(const SimpleIntArray& b);
01654             // copy the elements of an array
01655    SimpleIntArray(const SimpleIntArray& b);
01656             // make a new array equal to an existing one
01657    int Size() const { return n; }
01658             // return the size of the array
01659    int size() const { return n; }
01660             // return the size of the array
01661    int* Data() { return a; }  // pointer to the data
01662    const int* Data() const { return a; }  // pointer to the data
01663    int* data() { return a; }  // pointer to the data
01664    const int* data() const { return a; }  // pointer to the data
01665    const int* const_data() const { return a; }  // pointer to the data
01666    void ReSize(int i, bool keep = false);
01667                               // change length, keep data if keep = true
01668    void resize(int i, bool keep = false) { ReSize(i, keep); }
01669                               // change length, keep data if keep = true
01670    void CleanUp() { ReSize(0); }
01671    NEW_DELETE(SimpleIntArray)
01672 };
01673 
01674 // ********************** C subscript classes ****************************
01675 
01676 class RealStarStar
01677 {
01678    Real** a;
01679 public:
01680    RealStarStar(Matrix& A);
01681    RealStarStar(const RealStarStar&);
01682    ~RealStarStar() { delete [] a; }
01683    void operator=(const RealStarStar&);
01684    operator Real**() { return a; }
01685 };
01686 
01687 class ConstRealStarStar
01688 {
01689    const Real** a;
01690 public:
01691    ConstRealStarStar(const Matrix& A);
01692    ConstRealStarStar(const ConstRealStarStar&);
01693    ~ConstRealStarStar() { delete [] a; }
01694    void operator=(const ConstRealStarStar&);
01695    operator const Real**() { return a; }
01696 };
01697 
01698 // *************************** exceptions ********************************/
01699 
01700 class NPDException : public Runtime_error     // Not positive definite
01701 {
01702 public:
01703    static unsigned long Select;          // for identifying exception
01704    NPDException(const GeneralMatrix&);
01705 };
01706 
01707 class ConvergenceException : public Runtime_error
01708 {
01709 public:
01710    static unsigned long Select;          // for identifying exception
01711    ConvergenceException(const GeneralMatrix& A);
01712    ConvergenceException(const char* c);
01713 };
01714 
01715 class SingularException : public Runtime_error
01716 {
01717 public:
01718    static unsigned long Select;          // for identifying exception
01719    SingularException(const GeneralMatrix& A);
01720 };
01721 
01722 class OverflowException : public Runtime_error
01723 {
01724 public:
01725    static unsigned long Select;          // for identifying exception
01726    OverflowException(const char* c);
01727 };
01728 
01729 class ProgramException : public Logic_error
01730 {
01731 protected:
01732    ProgramException();
01733 public:
01734    static unsigned long Select;          // for identifying exception
01735    ProgramException(const char* c);
01736    ProgramException(const char* c, const GeneralMatrix&);
01737    ProgramException(const char* c, const GeneralMatrix&, const GeneralMatrix&);
01738    ProgramException(const char* c, MatrixType, MatrixType);
01739 };
01740 
01741 class IndexException : public Logic_error
01742 {
01743 public:
01744    static unsigned long Select;          // for identifying exception
01745    IndexException(int i, const GeneralMatrix& A);
01746    IndexException(int i, int j, const GeneralMatrix& A);
01747    // next two are for access via element function
01748    IndexException(int i, const GeneralMatrix& A, bool);
01749    IndexException(int i, int j, const GeneralMatrix& A, bool);
01750 };
01751 
01752 class VectorException : public Logic_error    // cannot convert to vector
01753 {
01754 public:
01755    static unsigned long Select;          // for identifying exception
01756    VectorException();
01757    VectorException(const GeneralMatrix& A);
01758 };
01759 
01760 class NotSquareException : public Logic_error
01761 {
01762 public:
01763    static unsigned long Select;          // for identifying exception
01764    NotSquareException(const GeneralMatrix& A);
01765    NotSquareException();
01766 };
01767 
01768 class SubMatrixDimensionException : public Logic_error
01769 {
01770 public:
01771    static unsigned long Select;          // for identifying exception
01772    SubMatrixDimensionException();
01773 };
01774 
01775 class IncompatibleDimensionsException : public Logic_error
01776 {
01777 public:
01778    static unsigned long Select;          // for identifying exception
01779    IncompatibleDimensionsException();
01780    IncompatibleDimensionsException(const GeneralMatrix&, const GeneralMatrix&);
01781 };
01782 
01783 class NotDefinedException : public Logic_error
01784 {
01785 public:
01786    static unsigned long Select;          // for identifying exception
01787    NotDefinedException(const char* op, const char* matrix);
01788 };
01789 
01790 class CannotBuildException : public Logic_error
01791 {
01792 public:
01793    static unsigned long Select;          // for identifying exception
01794    CannotBuildException(const char* matrix);
01795 };
01796 
01797 
01798 class InternalException : public Logic_error
01799 {
01800 public:
01801    static unsigned long Select;          // for identifying exception
01802    InternalException(const char* c);
01803 };
01804 
01805 // ************************ functions ************************************ //
01806 
01807 bool operator==(const GeneralMatrix& A, const GeneralMatrix& B);
01808 bool operator==(const BaseMatrix& A, const BaseMatrix& B);
01809 inline bool operator!=(const GeneralMatrix& A, const GeneralMatrix& B)
01810    { return ! (A==B); }
01811 inline bool operator!=(const BaseMatrix& A, const BaseMatrix& B)
01812    { return ! (A==B); }
01813 
01814    // inequality operators are dummies included for compatibility
01815    // with STL. They throw an exception if actually called.
01816 inline bool operator<=(const BaseMatrix& A, const BaseMatrix&)
01817    { A.IEQND(); return true; }
01818 inline bool operator>=(const BaseMatrix& A, const BaseMatrix&)
01819    { A.IEQND(); return true; }
01820 inline bool operator<(const BaseMatrix& A, const BaseMatrix&)
01821    { A.IEQND(); return true; }
01822 inline bool operator>(const BaseMatrix& A, const BaseMatrix&)
01823    { A.IEQND(); return true; }
01824 
01825 // these had been prototypes via 'friend' declaration, but
01826 // apparently newer versions of gcc now require an explicit prototype
01827 bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
01828 bool Compare(const MatrixType&, MatrixType&);
01829 Real DotProduct(const Matrix& A, const Matrix& B);
01830 SPMatrix SP(const BaseMatrix&, const BaseMatrix&);
01831 KPMatrix KP(const BaseMatrix&, const BaseMatrix&);
01832 ShiftedMatrix operator+(Real f, const BaseMatrix& BM);
01833 NegShiftedMatrix operator-(Real, const BaseMatrix&);
01834 ScaledMatrix operator*(Real f, const BaseMatrix& BM);
01835 
01836 bool IsZero(const BaseMatrix& A);
01837 
01838 //! mutliply two quaternions together, producing a new quaternion (elements in order w,x,y,z)
01839 ReturnMatrix MultiplyQuaternions(const ColumnVector& quatA, const ColumnVector& quatB);
01840 //! rotate columns of a matrix @a m by quaternion @a q (elements in order w,x,y,z)
01841 /*! Only rotates first 3 dimensions, copies "higher" dimensions unchanged so return matrix is same dimensions as @a m. */
01842 ReturnMatrix ApplyQuaternion(const ColumnVector& q, const Matrix& m);
01843 
01844 Matrix CrossProduct(const Matrix& A, const Matrix& B);
01845 ReturnMatrix CrossProductRows(const Matrix& A, const Matrix& B);
01846 ReturnMatrix CrossProductColumns(const Matrix& A, const Matrix& B);
01847 
01848 
01849 // ********************* inline functions ******************************** //
01850 
01851 
01852 inline LogAndSign LogDeterminant(const BaseMatrix& B)
01853    { return B.LogDeterminant(); }
01854 inline Real Determinant(const BaseMatrix& B)
01855    { return B.Determinant(); }
01856 inline Real SumSquare(const BaseMatrix& B) { return B.SumSquare(); }
01857 inline Real NormFrobenius(const BaseMatrix& B) { return B.NormFrobenius(); }
01858 inline Real Trace(const BaseMatrix& B) { return B.Trace(); }
01859 inline Real SumAbsoluteValue(const BaseMatrix& B)
01860    { return B.SumAbsoluteValue(); }
01861 inline Real Sum(const BaseMatrix& B)
01862    { return B.Sum(); }
01863 inline Real MaximumAbsoluteValue(const BaseMatrix& B)
01864    { return B.MaximumAbsoluteValue(); }
01865 inline Real MinimumAbsoluteValue(const BaseMatrix& B)
01866    { return B.MinimumAbsoluteValue(); }
01867 inline Real Maximum(const BaseMatrix& B) { return B.Maximum(); }
01868 inline Real Minimum(const BaseMatrix& B) { return B.Minimum(); }
01869 inline Real Norm1(const BaseMatrix& B) { return B.Norm1(); }
01870 inline Real Norm1(RowVector& RV) { return RV.MaximumAbsoluteValue(); }
01871 inline Real NormInfinity(const BaseMatrix& B) { return B.NormInfinity(); }
01872 inline Real NormInfinity(ColumnVector& CV)
01873    { return CV.MaximumAbsoluteValue(); }
01874 inline bool IsZero(const GeneralMatrix& A) { return A.IsZero(); }
01875 
01876 
01877 inline MatrixInput MatrixInput::operator<<(int f) { return *this << (Real)f; }
01878 inline MatrixInput GeneralMatrix::operator<<(int f) { return *this << (Real)f; }
01879 inline MatrixInput BandMatrix::operator<<(int f) { return *this << (Real)f; }
01880 inline MatrixInput GetSubMatrix::operator<<(int f) { return *this << (Real)f; }
01881 
01882 inline void swap(Matrix& A, Matrix& B) { A.swap(B); }
01883 inline void swap(SquareMatrix& A, SquareMatrix& B) { A.swap(B); }
01884 inline void swap(nricMatrix& A, nricMatrix& B) { A.swap(B); }
01885 inline void swap(UpperTriangularMatrix& A, UpperTriangularMatrix& B)
01886    { A.swap(B); }
01887 inline void swap(LowerTriangularMatrix& A, LowerTriangularMatrix& B)
01888    { A.swap(B); }
01889 inline void swap(SymmetricMatrix& A, SymmetricMatrix& B) { A.swap(B); }
01890 inline void swap(DiagonalMatrix& A, DiagonalMatrix& B) { A.swap(B); }
01891 inline void swap(RowVector& A, RowVector& B) { A.swap(B); }
01892 inline void swap(ColumnVector& A, ColumnVector& B) { A.swap(B); }
01893 inline void swap(CroutMatrix& A, CroutMatrix& B) { A.swap(B); }
01894 inline void swap(BandMatrix& A, BandMatrix& B) { A.swap(B); }
01895 inline void swap(UpperBandMatrix& A, UpperBandMatrix& B) { A.swap(B); }
01896 inline void swap(LowerBandMatrix& A, LowerBandMatrix& B) { A.swap(B); }
01897 inline void swap(SymmetricBandMatrix& A, SymmetricBandMatrix& B) { A.swap(B); }
01898 inline void swap(BandLUMatrix& A, BandLUMatrix& B) { A.swap(B); }
01899 inline void swap(IdentityMatrix& A, IdentityMatrix& B) { A.swap(B); }
01900 inline void swap(GenericMatrix& A, GenericMatrix& B) { A.swap(B); }
01901 
01902 
01903 #ifdef OPT_COMPATIBLE                    // for compatibility with opt++
01904 
01905 inline Real Norm2(const ColumnVector& CV) { return CV.NormFrobenius(); }
01906 inline Real Dot(ColumnVector& CV1, ColumnVector& CV2)
01907    { return DotProduct(CV1, CV2); }
01908 
01909 #endif
01910 
01911 
01912 #ifdef use_namespace
01913 }
01914 #endif
01915 
01916 
01917 #endif
01918 
01919 // body file: newmat1.cpp
01920 // body file: newmat2.cpp
01921 // body file: newmat3.cpp
01922 // body file: newmat4.cpp
01923 // body file: newmat5.cpp
01924 // body file: newmat6.cpp
01925 // body file: newmat7.cpp
01926 // body file: newmat8.cpp
01927 // body file: newmatex.cpp
01928 // body file: bandmat.cpp
01929 // body file: submat.cpp
01930 
01931 
01932 
01933 
01934 
01935 
01936 

newmat11b
Generated Mon May 9 04:54:18 2016 by Doxygen 1.6.3