NEWMAT6.CPP Source File Main Page Namespace List Class Hierarchy Alphabetical List Compound List File List Namespace Members Compound Members File Members
NEWMAT6.CPP

Go to the documentation of this file.
00001 //Owner: Fred
00002 //$$ newmat6.cpp            Operators, element access, submatrices
00003 
00004 // Copyright (C) 1991,2,3,4: R B Davies
00005 
00006 #include "include.h"
00007 
00008 #include "newmat.h"
00009 #include "newmatrc.h"
00010 
00011 #ifdef use_namespace
00012 namespace NEWMAT {
00013 #endif
00014 
00015 #ifdef DO_REPORT
00016 #define REPORT { static ExeCounter ExeCount(__LINE__,6); ++ExeCount; }
00017 #else
00018 #define REPORT {}
00019 #endif
00020 
00021     /*************************** general utilities *************************/
00022 
00023     static int tristore(int n)                      // els in triangular matrix
00024     { return (n*(n+1))/2; }
00025 
00026     /****************************** operators *******************************/
00027 
00028     Real& Matrix::operator()(int m, int n) {
00029         if (m<=0 || m>nrows || n<=0 || n>ncols)
00030             Throw(IndexException(m,n,*this));
00031         return store[(m-1)*ncols+n-1];
00032     }
00033 
00034     Real& SymmetricMatrix::operator()(int m, int n) {
00035         if (m<=0 || n<=0 || m>nrows || n>ncols)
00036             Throw(IndexException(m,n,*this));
00037         if (m>=n) return store[tristore(m-1)+n-1];
00038         else return store[tristore(n-1)+m-1];
00039     }
00040 
00041     Real& UpperTriangularMatrix::operator()(int m, int n) {
00042         if (m<=0 || n<m || n>ncols)
00043             Throw(IndexException(m,n,*this));
00044         return store[(m-1)*ncols+n-1-tristore(m-1)];
00045     }
00046 
00047     Real& LowerTriangularMatrix::operator()(int m, int n) {
00048         if (n<=0 || m<n || m>nrows)
00049             Throw(IndexException(m,n,*this));
00050         return store[tristore(m-1)+n-1];
00051     }
00052 
00053     Real& DiagonalMatrix::operator()(int m, int n) {
00054         if (n<=0 || m!=n || m>nrows || n>ncols)
00055             Throw(IndexException(m,n,*this));
00056         return store[n-1];
00057     }
00058 
00059     Real& DiagonalMatrix::operator()(int m) {
00060         if (m<=0 || m>nrows) Throw(IndexException(m,*this));
00061         return store[m-1];
00062     }
00063 
00064     Real& ColumnVector::operator()(int m) {
00065         if (m<=0 || m> nrows) Throw(IndexException(m,*this));
00066         return store[m-1];
00067     }
00068 
00069     Real& RowVector::operator()(int n) {
00070         if (n<=0 || n> ncols) Throw(IndexException(n,*this));
00071         return store[n-1];
00072     }
00073 
00074     Real& BandMatrix::operator()(int m, int n) {
00075         int w = upper+lower+1; int i = lower+n-m;
00076         if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
00077             Throw(IndexException(m,n,*this));
00078         return store[w*(m-1)+i];
00079     }
00080 
00081     Real& UpperBandMatrix::operator()(int m, int n) {
00082         int w = upper+1; int i = n-m;
00083         if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
00084             Throw(IndexException(m,n,*this));
00085         return store[w*(m-1)+i];
00086     }
00087 
00088     Real& LowerBandMatrix::operator()(int m, int n) {
00089         int w = lower+1; int i = lower+n-m;
00090         if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
00091             Throw(IndexException(m,n,*this));
00092         return store[w*(m-1)+i];
00093     }
00094 
00095     Real& SymmetricBandMatrix::operator()(int m, int n) {
00096         int w = lower+1;
00097         if (m>=n) {
00098             int i = lower+n-m;
00099             if ( m>nrows || n<=0 || i<0 )
00100                 Throw(IndexException(m,n,*this));
00101             return store[w*(m-1)+i];
00102         }
00103         else {
00104             int i = lower+m-n;
00105             if ( n>nrows || m<=0 || i<0 )
00106                 Throw(IndexException(m,n,*this));
00107             return store[w*(n-1)+i];
00108         }
00109     }
00110 
00111     Real Matrix::operator()(int m, int n) const {
00112         if (m<=0 || m>nrows || n<=0 || n>ncols)
00113             Throw(IndexException(m,n,*this));
00114         return store[(m-1)*ncols+n-1];
00115     }
00116 
00117     Real SymmetricMatrix::operator()(int m, int n) const {
00118         if (m<=0 || n<=0 || m>nrows || n>ncols)
00119             Throw(IndexException(m,n,*this));
00120         if (m>=n) return store[tristore(m-1)+n-1];
00121         else return store[tristore(n-1)+m-1];
00122     }
00123 
00124     Real UpperTriangularMatrix::operator()(int m, int n) const {
00125         if (m<=0 || n<m || n>ncols)
00126             Throw(IndexException(m,n,*this));
00127         return store[(m-1)*ncols+n-1-tristore(m-1)];
00128     }
00129 
00130     Real LowerTriangularMatrix::operator()(int m, int n) const {
00131         if (n<=0 || m<n || m>nrows)
00132             Throw(IndexException(m,n,*this));
00133         return store[tristore(m-1)+n-1];
00134     }
00135 
00136     Real DiagonalMatrix::operator()(int m, int n) const {
00137         if (n<=0 || m!=n || m>nrows || n>ncols)
00138             Throw(IndexException(m,n,*this));
00139         return store[n-1];
00140     }
00141 
00142     Real DiagonalMatrix::operator()(int m) const {
00143         if (m<=0 || m>nrows) Throw(IndexException(m,*this));
00144         return store[m-1];
00145     }
00146 
00147     Real ColumnVector::operator()(int m) const {
00148         if (m<=0 || m> nrows) Throw(IndexException(m,*this));
00149         return store[m-1];
00150     }
00151 
00152     Real RowVector::operator()(int n) const {
00153         if (n<=0 || n> ncols) Throw(IndexException(n,*this));
00154         return store[n-1];
00155     }
00156 
00157     Real BandMatrix::operator()(int m, int n) const {
00158         int w = upper+lower+1; int i = lower+n-m;
00159         if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
00160             Throw(IndexException(m,n,*this));
00161         return store[w*(m-1)+i];
00162     }
00163 
00164     Real SymmetricBandMatrix::operator()(int m, int n) const {
00165         int w = lower+1;
00166         if (m>=n) {
00167             int i = lower+n-m;
00168             if ( m>nrows || n<=0 || i<0 )
00169                 Throw(IndexException(m,n,*this));
00170             return store[w*(m-1)+i];
00171         }
00172         else {
00173             int i = lower+m-n;
00174             if ( n>nrows || m<=0 || i<0 )
00175                 Throw(IndexException(m,n,*this));
00176             return store[w*(n-1)+i];
00177         }
00178     }
00179 
00180     Real BaseMatrix::AsScalar() const {
00181         REPORT
00182             GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
00183 
00184         if (gm->nrows!=1 || gm->ncols!=1) {
00185             Tracer tr("AsScalar");
00186             Try {
00187                 Throw(ProgramException("Cannot convert to scalar", *gm));
00188             }
00189             CatchAll { gm->tDelete(); ReThrow; }
00190         }
00191 
00192         Real x = *(gm->store); gm->tDelete(); return x;
00193     }
00194 
00195 #ifdef TEMPS_DESTROYED_QUICKLY
00196 
00197     AddedMatrix& BaseMatrix::operator+(const BaseMatrix& bm) const {
00198         REPORT
00199             AddedMatrix* x = new AddedMatrix(this, &bm);
00200         MatrixErrorNoSpace(x);
00201         return *x;
00202     }
00203 
00204     SPMatrix& SP(const BaseMatrix& bm1,const BaseMatrix& bm2) {
00205         REPORT
00206             SPMatrix* x = new SPMatrix(&bm1, &bm2);
00207         MatrixErrorNoSpace(x);
00208         return *x;
00209     }
00210 
00211     ElmDivideMatrix& ElmDivide(const BaseMatrix& bm1,const BaseMatrix& bm2) {
00212         REPORT
00213             ElmDivideMatrix* x = new ElmDivideMatrix(&bm1, &bm2);
00214         MatrixErrorNoSpace(x);
00215         return *x;
00216     }
00217 
00218     MultipliedMatrix& BaseMatrix::operator*(const BaseMatrix& bm) const {
00219         REPORT
00220             MultipliedMatrix* x = new MultipliedMatrix(this, &bm);
00221         MatrixErrorNoSpace(x);
00222         return *x;
00223     }
00224 
00225     ConcatenatedMatrix& BaseMatrix::operator|(const BaseMatrix& bm) const {
00226         REPORT
00227             ConcatenatedMatrix* x = new ConcatenatedMatrix(this, &bm);
00228         MatrixErrorNoSpace(x);
00229         return *x;
00230     }
00231 
00232     StackedMatrix& BaseMatrix::operator&(const BaseMatrix& bm) const {
00233         REPORT
00234             StackedMatrix* x = new StackedMatrix(this, &bm);
00235         MatrixErrorNoSpace(x);
00236         return *x;
00237     }
00238 
00239     //SolvedMatrix& InvertedMatrix::operator*(const BaseMatrix& bmx) const
00240     SolvedMatrix& InvertedMatrix::operator*(const BaseMatrix& bmx) {
00241         REPORT
00242             SolvedMatrix* x;
00243         Try { x = new SolvedMatrix(bm, &bmx); MatrixErrorNoSpace(x); }
00244         CatchAll { delete this; ReThrow; }
00245         delete this;                                  // since we are using bm rather than this
00246         return *x;
00247     }
00248 
00249     SubtractedMatrix& BaseMatrix::operator-(const BaseMatrix& bm) const {
00250         REPORT
00251             SubtractedMatrix* x = new SubtractedMatrix(this, &bm);
00252         MatrixErrorNoSpace(x);
00253         return *x;
00254     }
00255 
00256     ShiftedMatrix& BaseMatrix::operator+(Real f) const {
00257         REPORT
00258             ShiftedMatrix* x = new ShiftedMatrix(this, f);
00259         MatrixErrorNoSpace(x);
00260         return *x;
00261     }
00262 
00263     NegShiftedMatrix& operator-(Real f,const BaseMatrix& bm1) {
00264         REPORT
00265             NegShiftedMatrix* x = new NegShiftedMatrix(f, &bm1);
00266         MatrixErrorNoSpace(x);
00267         return *x;
00268     }
00269 
00270     ScaledMatrix& BaseMatrix::operator*(Real f) const {
00271         REPORT
00272             ScaledMatrix* x = new ScaledMatrix(this, f);
00273         MatrixErrorNoSpace(x);
00274         return *x;
00275     }
00276 
00277     ScaledMatrix& BaseMatrix::operator/(Real f) const {
00278         REPORT
00279             ScaledMatrix* x = new ScaledMatrix(this, 1.0/f);
00280         MatrixErrorNoSpace(x);
00281         return *x;
00282     }
00283 
00284     ShiftedMatrix& BaseMatrix::operator-(Real f) const {
00285         REPORT
00286             ShiftedMatrix* x = new ShiftedMatrix(this, -f);
00287         MatrixErrorNoSpace(x);
00288         return *x;
00289     }
00290 
00291     TransposedMatrix& BaseMatrix::t() const {
00292         REPORT
00293             TransposedMatrix* x = new TransposedMatrix(this);
00294         MatrixErrorNoSpace(x);
00295         return *x;
00296     }
00297 
00298     NegatedMatrix& BaseMatrix::operator-() const {
00299         REPORT
00300             NegatedMatrix* x = new NegatedMatrix(this);
00301         MatrixErrorNoSpace(x);
00302         return *x;
00303     }
00304 
00305     ReversedMatrix& BaseMatrix::Reverse() const {
00306         REPORT
00307             ReversedMatrix* x = new ReversedMatrix(this);
00308         MatrixErrorNoSpace(x);
00309         return *x;
00310     }
00311 
00312     InvertedMatrix& BaseMatrix::i() const {
00313         REPORT
00314             InvertedMatrix* x = new InvertedMatrix(this);
00315         MatrixErrorNoSpace(x);
00316         return *x;
00317     }
00318 
00319     RowedMatrix& BaseMatrix::AsRow() const {
00320         REPORT
00321             RowedMatrix* x = new RowedMatrix(this);
00322         MatrixErrorNoSpace(x);
00323         return *x;
00324     }
00325 
00326     ColedMatrix& BaseMatrix::AsColumn() const {
00327         REPORT
00328             ColedMatrix* x = new ColedMatrix(this);
00329         MatrixErrorNoSpace(x);
00330         return *x;
00331     }
00332 
00333     DiagedMatrix& BaseMatrix::AsDiagonal() const {
00334         REPORT
00335             DiagedMatrix* x = new DiagedMatrix(this);
00336         MatrixErrorNoSpace(x);
00337         return *x;
00338     }
00339 
00340     MatedMatrix& BaseMatrix::AsMatrix(int nrx, int ncx) const {
00341         REPORT
00342             MatedMatrix* x = new MatedMatrix(this,nrx,ncx);
00343         MatrixErrorNoSpace(x);
00344         return *x;
00345     }
00346 
00347 #else
00348 
00349     AddedMatrix BaseMatrix::operator+(const BaseMatrix& bm) const
00350     { REPORT return AddedMatrix(this, &bm); }
00351 
00352     SPMatrix SP(const BaseMatrix& bm1,const BaseMatrix& bm2)
00353     { REPORT return SPMatrix(&bm1, &bm2); }
00354 
00355     ElmDivideMatrix ElmDivide(const BaseMatrix& bm1,const BaseMatrix& bm2)
00356     { REPORT return ElmDivideMatrix(&bm1, &bm2); }
00357 
00358     MultipliedMatrix BaseMatrix::operator*(const BaseMatrix& bm) const
00359     { REPORT return MultipliedMatrix(this, &bm); }
00360 
00361     ConcatenatedMatrix BaseMatrix::operator|(const BaseMatrix& bm) const
00362     { REPORT return ConcatenatedMatrix(this, &bm); }
00363 
00364     StackedMatrix BaseMatrix::operator&(const BaseMatrix& bm) const
00365     { REPORT return StackedMatrix(this, &bm); }
00366 
00367     SolvedMatrix InvertedMatrix::operator*(const BaseMatrix& bmx) const
00368     { REPORT return SolvedMatrix(bm, &bmx); }
00369 
00370     SubtractedMatrix BaseMatrix::operator-(const BaseMatrix& bm) const
00371     { REPORT return SubtractedMatrix(this, &bm); }
00372 
00373     ShiftedMatrix BaseMatrix::operator+(Real f) const
00374     { REPORT return ShiftedMatrix(this, f); }
00375 
00376     NegShiftedMatrix operator-(Real f, const BaseMatrix& bm)
00377     { REPORT return NegShiftedMatrix(f, &bm); }
00378 
00379     ScaledMatrix BaseMatrix::operator*(Real f) const
00380     { REPORT return ScaledMatrix(this, f); }
00381 
00382     ScaledMatrix BaseMatrix::operator/(Real f) const
00383     { REPORT return ScaledMatrix(this, 1.0/f); }
00384 
00385     ShiftedMatrix BaseMatrix::operator-(Real f) const
00386     { REPORT return ShiftedMatrix(this, -f); }
00387 
00388     TransposedMatrix BaseMatrix::t() const
00389     { REPORT return TransposedMatrix(this); }
00390 
00391     NegatedMatrix BaseMatrix::operator-() const
00392     { REPORT return NegatedMatrix(this); }
00393 
00394     ReversedMatrix BaseMatrix::Reverse() const
00395     { REPORT return ReversedMatrix(this); }
00396 
00397     InvertedMatrix BaseMatrix::i() const
00398     { REPORT return InvertedMatrix(this); }
00399 
00400     RowedMatrix BaseMatrix::AsRow() const
00401     { REPORT return RowedMatrix(this); }
00402 
00403     ColedMatrix BaseMatrix::AsColumn() const
00404     { REPORT return ColedMatrix(this); }
00405 
00406     DiagedMatrix BaseMatrix::AsDiagonal() const
00407     { REPORT return DiagedMatrix(this); }
00408 
00409     MatedMatrix BaseMatrix::AsMatrix(int nrx, int ncx) const
00410     { REPORT return MatedMatrix(this,nrx,ncx); }
00411 #endif
00412 
00413     void GeneralMatrix::operator=(Real f) {
00414         REPORT int i=storage; Real* s=store; while (i--) {
00415             *s++ = f;
00416         }
00417     }
00418 
00419     void Matrix::operator=(const BaseMatrix& X) {
00420         REPORT                                        //CheckConversion(X);
00421             MatrixConversionCheck mcc;
00422         Eq(X,MatrixType::Rt);
00423     }
00424 
00425     void RowVector::operator=(const BaseMatrix& X) {
00426         REPORT                                        // CheckConversion(X);
00427             MatrixConversionCheck mcc;
00428         Eq(X,MatrixType::RV);
00429         if (nrows!=1) {
00430             Tracer tr("RowVector(=)");
00431             Throw(VectorException(*this));
00432         }
00433     }
00434 
00435     void ColumnVector::operator=(const BaseMatrix& X) {
00436         REPORT                                        //CheckConversion(X);
00437             MatrixConversionCheck mcc;
00438         Eq(X,MatrixType::CV);
00439         if (ncols!=1) {
00440             Tracer tr("ColumnVector(=)");
00441             Throw(VectorException(*this));
00442         }
00443     }
00444 
00445     void SymmetricMatrix::operator=(const BaseMatrix& X) {
00446         REPORT                                        // CheckConversion(X);
00447             MatrixConversionCheck mcc;
00448         Eq(X,MatrixType::Sm);
00449     }
00450 
00451     void UpperTriangularMatrix::operator=(const BaseMatrix& X) {
00452         REPORT                                        //CheckConversion(X);
00453             MatrixConversionCheck mcc;
00454         Eq(X,MatrixType::UT);
00455     }
00456 
00457     void LowerTriangularMatrix::operator=(const BaseMatrix& X) {
00458         REPORT                                        //CheckConversion(X);
00459             MatrixConversionCheck mcc;
00460         Eq(X,MatrixType::LT);
00461     }
00462 
00463     void DiagonalMatrix::operator=(const BaseMatrix& X) {
00464         REPORT                                        // CheckConversion(X);
00465             MatrixConversionCheck mcc;
00466         Eq(X,MatrixType::Dg);
00467     }
00468 
00469     void GeneralMatrix::operator<<(const Real* r) {
00470         REPORT
00471             int i = storage; Real* s=store;
00472         while(i--) *s++ = *r++;
00473     }
00474 
00475     void GenericMatrix::operator=(const GenericMatrix& bmx) {
00476         if (&bmx != this) { REPORT if (gm) delete gm; gm = bmx.gm->Image();}
00477         else { REPORT }
00478         gm->Protect();
00479     }
00480 
00481     void GenericMatrix::operator=(const BaseMatrix& bmx) {
00482         if (gm) {
00483             int counter=bmx.search(gm);
00484             if (counter==0) { REPORT delete gm; gm=0; }
00485             else { REPORT gm->Release(counter); }
00486         }
00487         else { REPORT }
00488         GeneralMatrix* gmx = ((BaseMatrix&)bmx).Evaluate();
00489         if (gmx != gm) { REPORT if (gm) delete gm; gm = gmx->Image(); }
00490         else { REPORT }
00491         gm->Protect();
00492     }
00493 
00494     /*************************** += etc ***************************************/
00495 
00496     // will also need versions for SubMatrix
00497 
00498     // GeneralMatrix operators
00499 
00500     void GeneralMatrix::operator+=(const BaseMatrix& X) {
00501         REPORT
00502             Tracer tr("GeneralMatrix::operator+=");
00503         MatrixConversionCheck mcc;
00504         Protect();                                    // so it can't get deleted
00505         // during Evaluate
00506         GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
00507 #ifdef TEMPS_DESTROYED_QUICKLY
00508         AddedMatrix* am = new AddedMatrix(this,gm);
00509         MatrixErrorNoSpace(am);
00510         if (gm==this) Release(2); else Release();
00511         Eq2(*am,Type());
00512 #else
00513         AddedMatrix am(this,gm);
00514         if (gm==this) Release(2); else Release();
00515         Eq2(am,Type());
00516 #endif
00517     }
00518 
00519     void GeneralMatrix::operator-=(const BaseMatrix& X) {
00520         REPORT
00521             Tracer tr("GeneralMatrix::operator-=");
00522         MatrixConversionCheck mcc;
00523         Protect();                                    // so it can't get deleted
00524         // during Evaluate
00525         GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
00526 #ifdef TEMPS_DESTROYED_QUICKLY
00527         SubtractedMatrix* am = new SubtractedMatrix(this,gm);
00528         MatrixErrorNoSpace(am);
00529         if (gm==this) Release(2); else Release();
00530         Eq2(*am,Type());
00531 #else
00532         SubtractedMatrix am(this,gm);
00533         if (gm==this) Release(2); else Release();
00534         Eq2(am,Type());
00535 #endif
00536     }
00537 
00538     void GeneralMatrix::operator*=(const BaseMatrix& X) {
00539         REPORT
00540             Tracer tr("GeneralMatrix::operator*=");
00541         MatrixConversionCheck mcc;
00542         Protect();                                    // so it can't get deleted
00543         // during Evaluate
00544         GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
00545 #ifdef TEMPS_DESTROYED_QUICKLY
00546         MultipliedMatrix* am = new MultipliedMatrix(this,gm);
00547         MatrixErrorNoSpace(am);
00548         if (gm==this) Release(2); else Release();
00549         Eq2(*am,Type());
00550 #else
00551         MultipliedMatrix am(this,gm);
00552         if (gm==this) Release(2); else Release();
00553         Eq2(am,Type());
00554 #endif
00555     }
00556 
00557     void GeneralMatrix::operator|=(const BaseMatrix& X) {
00558         REPORT
00559             Tracer tr("GeneralMatrix::operator|=");
00560         MatrixConversionCheck mcc;
00561         Protect();                                    // so it can't get deleted
00562         // during Evaluate
00563         GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
00564 #ifdef TEMPS_DESTROYED_QUICKLY
00565         ConcatenatedMatrix* am = new ConcatenatedMatrix(this,gm);
00566         MatrixErrorNoSpace(am);
00567         if (gm==this) Release(2); else Release();
00568         Eq2(*am,Type());
00569 #else
00570         ConcatenatedMatrix am(this,gm);
00571         if (gm==this) Release(2); else Release();
00572         Eq2(am,Type());
00573 #endif
00574     }
00575 
00576     void GeneralMatrix::operator&=(const BaseMatrix& X) {
00577         REPORT
00578             Tracer tr("GeneralMatrix::operator&=");
00579         MatrixConversionCheck mcc;
00580         Protect();                                    // so it can't get deleted
00581         // during Evaluate
00582         GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
00583 #ifdef TEMPS_DESTROYED_QUICKLY
00584         StackedMatrix* am = new StackedMatrix(this,gm);
00585         MatrixErrorNoSpace(am);
00586         if (gm==this) Release(2); else Release();
00587         Eq2(*am,Type());
00588 #else
00589         StackedMatrix am(this,gm);
00590         if (gm==this) Release(2); else Release();
00591         Eq2(am,Type());
00592 #endif
00593     }
00594 
00595     void GeneralMatrix::operator+=(Real r) {
00596         REPORT
00597             Tracer tr("GeneralMatrix::operator+=(Real)");
00598         MatrixConversionCheck mcc;
00599 #ifdef TEMPS_DESTROYED_QUICKLY
00600         ShiftedMatrix* am = new ShiftedMatrix(this,r);
00601         MatrixErrorNoSpace(am);
00602         Release(); Eq2(*am,Type());
00603 #else
00604         ShiftedMatrix am(this,r);
00605         Release(); Eq2(am,Type());
00606 #endif
00607     }
00608 
00609     void GeneralMatrix::operator*=(Real r) {
00610         REPORT
00611             Tracer tr("GeneralMatrix::operator*=(Real)");
00612         MatrixConversionCheck mcc;
00613 #ifdef TEMPS_DESTROYED_QUICKLY
00614         ScaledMatrix* am = new ScaledMatrix(this,r);
00615         MatrixErrorNoSpace(am);
00616         Release(); Eq2(*am,Type());
00617 #else
00618         ScaledMatrix am(this,r);
00619         Release(); Eq2(am,Type());
00620 #endif
00621     }
00622 
00623     // Generic matrix operators
00624 
00625     void GenericMatrix::operator+=(const BaseMatrix& X) {
00626         REPORT
00627             Tracer tr("GenericMatrix::operator+=");
00628         if (!gm) Throw(ProgramException("GenericMatrix is null"));
00629         gm->Protect();                                // so it can't get deleted during Evaluate
00630         GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
00631 #ifdef TEMPS_DESTROYED_QUICKLY
00632         AddedMatrix* am = new AddedMatrix(gm,gmx);
00633         MatrixErrorNoSpace(am);
00634         if (gmx==gm) gm->Release(2); else gm->Release();
00635         GeneralMatrix* gmy = am->Evaluate();
00636 #else
00637         AddedMatrix am(gm,gmx);
00638         if (gmx==gm) gm->Release(2); else gm->Release();
00639         GeneralMatrix* gmy = am.Evaluate();
00640 #endif
00641         if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
00642         else { REPORT }
00643         gm->Protect();
00644     }
00645 
00646     void GenericMatrix::operator-=(const BaseMatrix& X) {
00647         REPORT
00648             Tracer tr("GenericMatrix::operator-=");
00649         if (!gm) Throw(ProgramException("GenericMatrix is null"));
00650         gm->Protect();                                // so it can't get deleted during Evaluate
00651         GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
00652 #ifdef TEMPS_DESTROYED_QUICKLY
00653         SubtractedMatrix* am = new SubtractedMatrix(gm,gmx);
00654         MatrixErrorNoSpace(am);
00655         if (gmx==gm) gm->Release(2); else gm->Release();
00656         GeneralMatrix* gmy = am->Evaluate();
00657 #else
00658         SubtractedMatrix am(gm,gmx);
00659         if (gmx==gm) gm->Release(2); else gm->Release();
00660         GeneralMatrix* gmy = am.Evaluate();
00661 #endif
00662         if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
00663         else { REPORT }
00664         gm->Protect();
00665     }
00666 
00667     void GenericMatrix::operator*=(const BaseMatrix& X) {
00668         REPORT
00669             Tracer tr("GenericMatrix::operator*=");
00670         if (!gm) Throw(ProgramException("GenericMatrix is null"));
00671         gm->Protect();                                // so it can't get deleted during Evaluate
00672         GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
00673 #ifdef TEMPS_DESTROYED_QUICKLY
00674         MultipliedMatrix* am = new MultipliedMatrix(gm,gmx);
00675         MatrixErrorNoSpace(am);
00676         if (gmx==gm) gm->Release(2); else gm->Release();
00677         GeneralMatrix* gmy = am->Evaluate();
00678 #else
00679         MultipliedMatrix am(gm,gmx);
00680         if (gmx==gm) gm->Release(2); else gm->Release();
00681         GeneralMatrix* gmy = am.Evaluate();
00682 #endif
00683         if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
00684         else { REPORT }
00685         gm->Protect();
00686     }
00687 
00688     void GenericMatrix::operator|=(const BaseMatrix& X) {
00689         REPORT
00690             Tracer tr("GenericMatrix::operator|=");
00691         if (!gm) Throw(ProgramException("GenericMatrix is null"));
00692         gm->Protect();                                // so it can't get deleted during Evaluate
00693         GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
00694 #ifdef TEMPS_DESTROYED_QUICKLY
00695         ConcatenatedMatrix* am = new ConcatenatedMatrix(gm,gmx);
00696         MatrixErrorNoSpace(am);
00697         if (gmx==gm) gm->Release(2); else gm->Release();
00698         GeneralMatrix* gmy = am->Evaluate();
00699 #else
00700         ConcatenatedMatrix am(gm,gmx);
00701         if (gmx==gm) gm->Release(2); else gm->Release();
00702         GeneralMatrix* gmy = am.Evaluate();
00703 #endif
00704         if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
00705         else { REPORT }
00706         gm->Protect();
00707     }
00708 
00709     void GenericMatrix::operator&=(const BaseMatrix& X) {
00710         REPORT
00711             Tracer tr("GenericMatrix::operator&=");
00712         if (!gm) Throw(ProgramException("GenericMatrix is null"));
00713         gm->Protect();                                // so it can't get deleted during Evaluate
00714         GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
00715 #ifdef TEMPS_DESTROYED_QUICKLY
00716         StackedMatrix* am = new StackedMatrix(gm,gmx);
00717         MatrixErrorNoSpace(am);
00718         if (gmx==gm) gm->Release(2); else gm->Release();
00719         GeneralMatrix* gmy = am->Evaluate();
00720 #else
00721         StackedMatrix am(gm,gmx);
00722         if (gmx==gm) gm->Release(2); else gm->Release();
00723         GeneralMatrix* gmy = am.Evaluate();
00724 #endif
00725         if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
00726         else { REPORT }
00727         gm->Protect();
00728     }
00729 
00730     void GenericMatrix::operator+=(Real r) {
00731         REPORT
00732             Tracer tr("GenericMatrix::operator+= (Real)");
00733         if (!gm) Throw(ProgramException("GenericMatrix is null"));
00734 #ifdef TEMPS_DESTROYED_QUICKLY
00735         ShiftedMatrix* am = new ShiftedMatrix(gm,r);
00736         MatrixErrorNoSpace(am);
00737         gm->Release();
00738         GeneralMatrix* gmy = am->Evaluate();
00739 #else
00740         ShiftedMatrix am(gm,r);
00741         gm->Release();
00742         GeneralMatrix* gmy = am.Evaluate();
00743 #endif
00744         if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
00745         else { REPORT }
00746         gm->Protect();
00747     }
00748 
00749     void GenericMatrix::operator*=(Real r) {
00750         REPORT
00751             Tracer tr("GenericMatrix::operator*= (Real)");
00752         if (!gm) Throw(ProgramException("GenericMatrix is null"));
00753 #ifdef TEMPS_DESTROYED_QUICKLY
00754         ScaledMatrix* am = new ScaledMatrix(gm,r);
00755         MatrixErrorNoSpace(am);
00756         gm->Release();
00757         GeneralMatrix* gmy = am->Evaluate();
00758 #else
00759         ScaledMatrix am(gm,r);
00760         gm->Release();
00761         GeneralMatrix* gmy = am.Evaluate();
00762 #endif
00763         if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
00764         else { REPORT }
00765         gm->Protect();
00766     }
00767 
00768     /************************* element access *********************************/
00769 
00770     Real& Matrix::element(int m, int n) {
00771         if (m<0 || m>= nrows || n<0 || n>= ncols)
00772             Throw(IndexException(m,n,*this,true));
00773         return store[m*ncols+n];
00774     }
00775 
00776     Real& SymmetricMatrix::element(int m, int n) {
00777         if (m<0 || n<0 || m >= nrows || n>=ncols)
00778             Throw(IndexException(m,n,*this,true));
00779         if (m>=n) return store[tristore(m)+n];
00780         else return store[tristore(n)+m];
00781     }
00782 
00783     Real& UpperTriangularMatrix::element(int m, int n) {
00784         if (m<0 || n<m || n>=ncols)
00785             Throw(IndexException(m,n,*this,true));
00786         return store[m*ncols+n-tristore(m)];
00787     }
00788 
00789     Real& LowerTriangularMatrix::element(int m, int n) {
00790         if (n<0 || m<n || m>=nrows)
00791             Throw(IndexException(m,n,*this,true));
00792         return store[tristore(m)+n];
00793     }
00794 
00795     Real& DiagonalMatrix::element(int m, int n) {
00796         if (n<0 || m!=n || m>=nrows || n>=ncols)
00797             Throw(IndexException(m,n,*this,true));
00798         return store[n];
00799     }
00800 
00801     Real& DiagonalMatrix::element(int m) {
00802         if (m<0 || m>=nrows) Throw(IndexException(m,*this,true));
00803         return store[m];
00804     }
00805 
00806     Real& ColumnVector::element(int m) {
00807         if (m<0 || m>= nrows) Throw(IndexException(m,*this,true));
00808         return store[m];
00809     }
00810 
00811     Real& RowVector::element(int n) {
00812         if (n<0 || n>= ncols)  Throw(IndexException(n,*this,true));
00813         return store[n];
00814     }
00815 
00816     Real& BandMatrix::element(int m, int n) {
00817         int w = upper+lower+1; int i = lower+n-m;
00818         if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
00819             Throw(IndexException(m,n,*this,true));
00820         return store[w*m+i];
00821     }
00822 
00823     Real& UpperBandMatrix::element(int m, int n) {
00824         int w = upper+1; int i = n-m;
00825         if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
00826             Throw(IndexException(m,n,*this,true));
00827         return store[w*m+i];
00828     }
00829 
00830     Real& LowerBandMatrix::element(int m, int n) {
00831         int w = lower+1; int i = lower+n-m;
00832         if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
00833             Throw(IndexException(m,n,*this,true));
00834         return store[w*m+i];
00835     }
00836 
00837     Real& SymmetricBandMatrix::element(int m, int n) {
00838         int w = lower+1;
00839         if (m>=n) {
00840             int i = lower+n-m;
00841             if ( m>=nrows || n<0 || i<0 )
00842                 Throw(IndexException(m,n,*this,true));
00843             return store[w*m+i];
00844         }
00845         else {
00846             int i = lower+m-n;
00847             if ( n>=nrows || m<0 || i<0 )
00848                 Throw(IndexException(m,n,*this,true));
00849             return store[w*n+i];
00850         }
00851     }
00852 
00853 #ifdef use_namespace
00854 }
00855 #endif
Generated on Fri Aug 23 01:37:09 2002 for VirtualU by
1.2.17