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


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