00001
00002
00003
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
00023
00024 static int tristore(int n)
00025 { return (n*(n+1))/2; }
00026
00027
00028
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
00348
00349 Eq(X,MatrixType::Rt);
00350 }
00351
00352 void SquareMatrix::operator=(const BaseMatrix& X)
00353 {
00354 REPORT
00355
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
00372
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
00381
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
00390
00391 Eq(X,MatrixType::Sm);
00392 }
00393
00394 void UpperTriangularMatrix::operator=(const BaseMatrix& X)
00395 {
00396 REPORT
00397
00398 Eq(X,MatrixType::UT);
00399 }
00400
00401 void LowerTriangularMatrix::operator=(const BaseMatrix& X)
00402 {
00403 REPORT
00404
00405 Eq(X,MatrixType::LT);
00406 }
00407
00408 void DiagonalMatrix::operator=(const BaseMatrix& X)
00409 {
00410 REPORT
00411
00412 Eq(X,MatrixType::Dg);
00413 }
00414
00415 void IdentityMatrix::operator=(const BaseMatrix& X)
00416 {
00417 REPORT
00418
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
00462
00463
00464
00465
00466
00467
00468
00469 void GeneralMatrix::operator+=(const BaseMatrix& X)
00470 {
00471 REPORT
00472 Tracer tr("GeneralMatrix::operator+=");
00473
00474 Protect();
00475
00476 GeneralMatrix* gm = (const_cast<BaseMatrix&>(X)).Evaluate();
00477 AddedMatrix am(this,gm);
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
00487 Protect();
00488
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
00500 Protect();
00501
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
00513 Protect();
00514
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
00526 Protect();
00527
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
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
00548 ScaledMatrix am(this,r);
00549 Release(); Eq2(am,Type());
00550 }
00551
00552
00553
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();
00561 GeneralMatrix* gmx = (const_cast<BaseMatrix&>(X)).Evaluate();
00562 AddedMatrix am(gm,gmx);
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();
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();
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();
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();
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
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