Tekkotsu Homepage Demos Overview Downloads Dev. Resources Reference Credits

# newmat6.cpp

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

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