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


kni
Author(s): Martin Günther
autogenerated on Mon Aug 14 2017 02:44:12