00001
00002
00003
00006
00007
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
00027
00028 static int tristore(int n)
00029 { return (n*(n+1))/2; }
00030
00031
00032
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
00352
00353 Eq(X,MatrixType::Rt);
00354 }
00355
00356 void SquareMatrix::operator=(const BaseMatrix& X)
00357 {
00358 REPORT
00359
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
00376
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
00385
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
00394
00395 Eq(X,MatrixType::Sm);
00396 }
00397
00398 void UpperTriangularMatrix::operator=(const BaseMatrix& X)
00399 {
00400 REPORT
00401
00402 Eq(X,MatrixType::UT);
00403 }
00404
00405 void LowerTriangularMatrix::operator=(const BaseMatrix& X)
00406 {
00407 REPORT
00408
00409 Eq(X,MatrixType::LT);
00410 }
00411
00412 void DiagonalMatrix::operator=(const BaseMatrix& X)
00413 {
00414 REPORT
00415
00416 Eq(X,MatrixType::Dg);
00417 }
00418
00419 void IdentityMatrix::operator=(const BaseMatrix& X)
00420 {
00421 REPORT
00422
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
00488
00489
00490
00491
00492 void GeneralMatrix::operator+=(const BaseMatrix& X)
00493 {
00494 REPORT
00495 Tracer tr("GeneralMatrix::operator+=");
00496
00497 Protect();
00498
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
00510 Protect();
00511
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
00523 Protect();
00524
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
00536 Protect();
00537
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
00549 Protect();
00550
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
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
00571 ScaledMatrix am(this,r);
00572 Release(); Eq2(am,type());
00573 }
00574
00575
00576
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();
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();
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();
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();
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();
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
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