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__,4); ++ExeCount; }
00018 #else
00019 #define REPORT {}
00020 #endif
00021
00022
00023 #define DO_SEARCH // search for LHS of = in RHS
00024
00025
00026
00027 static int tristore(int n)
00028 { return (n*(n+1))/2; }
00029
00030
00031
00032
00033 GeneralMatrix::GeneralMatrix()
00034 { store=0; storage=0; nrows_value=0; ncols_value=0; tag=-1; }
00035
00036 GeneralMatrix::GeneralMatrix(ArrayLengthSpecifier s)
00037 {
00038 REPORT
00039 storage=s.Value(); tag=-1;
00040 if (storage)
00041 {
00042 store = new Real [storage]; MatrixErrorNoSpace(store);
00043 MONITOR_REAL_NEW("Make (GenMatrix)",storage,store)
00044 }
00045 else store = 0;
00046 }
00047
00048 Matrix::Matrix(int m, int n) : GeneralMatrix(m*n)
00049 { REPORT nrows_value=m; ncols_value=n; }
00050
00051 SquareMatrix::SquareMatrix(ArrayLengthSpecifier n)
00052 : Matrix(n.Value(),n.Value())
00053 { REPORT }
00054
00055 SymmetricMatrix::SymmetricMatrix(ArrayLengthSpecifier n)
00056 : GeneralMatrix(tristore(n.Value()))
00057 { REPORT nrows_value=n.Value(); ncols_value=n.Value(); }
00058
00059 UpperTriangularMatrix::UpperTriangularMatrix(ArrayLengthSpecifier n)
00060 : GeneralMatrix(tristore(n.Value()))
00061 { REPORT nrows_value=n.Value(); ncols_value=n.Value(); }
00062
00063 LowerTriangularMatrix::LowerTriangularMatrix(ArrayLengthSpecifier n)
00064 : GeneralMatrix(tristore(n.Value()))
00065 { REPORT nrows_value=n.Value(); ncols_value=n.Value(); }
00066
00067 DiagonalMatrix::DiagonalMatrix(ArrayLengthSpecifier m) : GeneralMatrix(m)
00068 { REPORT nrows_value=m.Value(); ncols_value=m.Value(); }
00069
00070 Matrix::Matrix(const BaseMatrix& M)
00071 {
00072 REPORT
00073
00074 GeneralMatrix* gmx=(const_cast<BaseMatrix&>(M)).Evaluate(MatrixType::Rt);
00075 GetMatrix(gmx);
00076 }
00077
00078 SquareMatrix::SquareMatrix(const BaseMatrix& M) : Matrix(M)
00079 {
00080 REPORT
00081 if (ncols_value != nrows_value)
00082 {
00083 Tracer tr("SquareMatrix");
00084 Throw(NotSquareException(*this));
00085 }
00086 }
00087
00088
00089 SquareMatrix::SquareMatrix(const Matrix& gm)
00090 {
00091 REPORT
00092 if (gm.ncols_value != gm.nrows_value)
00093 {
00094 Tracer tr("SquareMatrix(Matrix)");
00095 Throw(NotSquareException(gm));
00096 }
00097 GetMatrix(&gm);
00098 }
00099
00100
00101 RowVector::RowVector(const BaseMatrix& M) : Matrix(M)
00102 {
00103 REPORT
00104 if (nrows_value!=1)
00105 {
00106 Tracer tr("RowVector");
00107 Throw(VectorException(*this));
00108 }
00109 }
00110
00111 ColumnVector::ColumnVector(const BaseMatrix& M) : Matrix(M)
00112 {
00113 REPORT
00114 if (ncols_value!=1)
00115 {
00116 Tracer tr("ColumnVector");
00117 Throw(VectorException(*this));
00118 }
00119 }
00120
00121 SymmetricMatrix::SymmetricMatrix(const BaseMatrix& M)
00122 {
00123 REPORT
00124
00125 GeneralMatrix* gmx=(const_cast<BaseMatrix&>(M)).Evaluate(MatrixType::Sm);
00126 GetMatrix(gmx);
00127 }
00128
00129 UpperTriangularMatrix::UpperTriangularMatrix(const BaseMatrix& M)
00130 {
00131 REPORT
00132
00133 GeneralMatrix* gmx=(const_cast<BaseMatrix&>(M)).Evaluate(MatrixType::UT);
00134 GetMatrix(gmx);
00135 }
00136
00137 LowerTriangularMatrix::LowerTriangularMatrix(const BaseMatrix& M)
00138 {
00139 REPORT
00140
00141 GeneralMatrix* gmx=(const_cast<BaseMatrix&>(M)).Evaluate(MatrixType::LT);
00142 GetMatrix(gmx);
00143 }
00144
00145 DiagonalMatrix::DiagonalMatrix(const BaseMatrix& M)
00146 {
00147 REPORT
00148
00149 GeneralMatrix* gmx=(const_cast<BaseMatrix&>(M)).Evaluate(MatrixType::Dg);
00150 GetMatrix(gmx);
00151 }
00152
00153 IdentityMatrix::IdentityMatrix(const BaseMatrix& M)
00154 {
00155 REPORT
00156
00157 GeneralMatrix* gmx=(const_cast<BaseMatrix&>(M)).Evaluate(MatrixType::Id);
00158 GetMatrix(gmx);
00159 }
00160
00161 GeneralMatrix::~GeneralMatrix()
00162 {
00163 if (store)
00164 {
00165 MONITOR_REAL_DELETE("Free (GenMatrix)",storage,store)
00166 delete [] store;
00167 }
00168 }
00169
00170 CroutMatrix::CroutMatrix(const BaseMatrix& m)
00171 {
00172 REPORT
00173 Tracer tr("CroutMatrix");
00174 indx = 0;
00175 GeneralMatrix* gm = (const_cast<BaseMatrix&>(m)).Evaluate(MatrixType::Rt);
00176 GetMatrix(gm);
00177 if (nrows_value!=ncols_value) { CleanUp(); Throw(NotSquareException(*gm)); }
00178 d=true; sing=false;
00179 indx=new int [nrows_value]; MatrixErrorNoSpace(indx);
00180 MONITOR_INT_NEW("Index (CroutMat)",nrows_value,indx)
00181 ludcmp();
00182 }
00183
00184 CroutMatrix::~CroutMatrix()
00185 {
00186 MONITOR_INT_DELETE("Index (CroutMat)",nrows_value,indx)
00187 delete [] indx;
00188 }
00189
00190
00191
00192
00193
00194
00195
00196
00197 GeneralMatrix::operator ReturnMatrix() const
00198 {
00199 REPORT
00200 GeneralMatrix* gm = Image(); gm->ReleaseAndDelete();
00201 return ReturnMatrix(gm);
00202 }
00203
00204
00205
00206 ReturnMatrix GeneralMatrix::ForReturn() const
00207 {
00208 REPORT
00209 GeneralMatrix* gm = Image(); gm->ReleaseAndDelete();
00210 return ReturnMatrix(gm);
00211 }
00212
00213
00214
00215 #ifdef SETUP_C_SUBSCRIPTS
00216
00217 Matrix::Matrix(Real a, int m, int n) : GeneralMatrix(m * n)
00218 { REPORT nrows_value=m; ncols_value=n; operator=(a); }
00219
00220 Matrix::Matrix(const Real* a, int m, int n) : GeneralMatrix(m * n)
00221 { REPORT nrows_value=m; ncols_value=n; *this << a; }
00222
00223 #endif
00224
00225
00226
00227
00228
00229 void GeneralMatrix::ReSize(int nr, int nc, int s)
00230 {
00231 REPORT
00232 if (store)
00233 {
00234 MONITOR_REAL_DELETE("Free (ReDimensi)",storage,store)
00235 delete [] store;
00236 }
00237 storage=s; nrows_value=nr; ncols_value=nc; tag=-1;
00238 if (s)
00239 {
00240 store = new Real [storage]; MatrixErrorNoSpace(store);
00241 MONITOR_REAL_NEW("Make (ReDimensi)",storage,store)
00242 }
00243 else store = 0;
00244 }
00245
00246 void Matrix::ReSize(int nr, int nc)
00247 { REPORT GeneralMatrix::ReSize(nr,nc,nr*nc); }
00248
00249 void SquareMatrix::ReSize(int n)
00250 { REPORT GeneralMatrix::ReSize(n,n,n*n); }
00251
00252 void SquareMatrix::ReSize(int nr, int nc)
00253 {
00254 REPORT
00255 Tracer tr("SquareMatrix::ReSize");
00256 if (nc != nr) Throw(NotSquareException(*this));
00257 GeneralMatrix::ReSize(nr,nc,nr*nc);
00258 }
00259
00260 void SymmetricMatrix::ReSize(int nr)
00261 { REPORT GeneralMatrix::ReSize(nr,nr,tristore(nr)); }
00262
00263 void UpperTriangularMatrix::ReSize(int nr)
00264 { REPORT GeneralMatrix::ReSize(nr,nr,tristore(nr)); }
00265
00266 void LowerTriangularMatrix::ReSize(int nr)
00267 { REPORT GeneralMatrix::ReSize(nr,nr,tristore(nr)); }
00268
00269 void DiagonalMatrix::ReSize(int nr)
00270 { REPORT GeneralMatrix::ReSize(nr,nr,nr); }
00271
00272 void RowVector::ReSize(int nc)
00273 { REPORT GeneralMatrix::ReSize(1,nc,nc); }
00274
00275 void ColumnVector::ReSize(int nr)
00276 { REPORT GeneralMatrix::ReSize(nr,1,nr); }
00277
00278 void RowVector::ReSize(int nr, int nc)
00279 {
00280 Tracer tr("RowVector::ReSize");
00281 if (nr != 1) Throw(VectorException(*this));
00282 REPORT GeneralMatrix::ReSize(1,nc,nc);
00283 }
00284
00285 void ColumnVector::ReSize(int nr, int nc)
00286 {
00287 Tracer tr("ColumnVector::ReSize");
00288 if (nc != 1) Throw(VectorException(*this));
00289 REPORT GeneralMatrix::ReSize(nr,1,nr);
00290 }
00291
00292 void IdentityMatrix::ReSize(int nr)
00293 { REPORT GeneralMatrix::ReSize(nr,nr,1); *store = 1; }
00294
00295
00296 void Matrix::ReSize(const GeneralMatrix& A)
00297 { REPORT ReSize(A.Nrows(), A.Ncols()); }
00298
00299 void SquareMatrix::ReSize(const GeneralMatrix& A)
00300 {
00301 REPORT
00302 int n = A.Nrows();
00303 if (n != A.Ncols())
00304 {
00305 Tracer tr("SquareMatrix::ReSize(GM)");
00306 Throw(NotSquareException(*this));
00307 }
00308 ReSize(n);
00309 }
00310
00311 void nricMatrix::ReSize(const GeneralMatrix& A)
00312 { REPORT ReSize(A.Nrows(), A.Ncols()); }
00313
00314 void ColumnVector::ReSize(const GeneralMatrix& A)
00315 { REPORT ReSize(A.Nrows(), A.Ncols()); }
00316
00317 void RowVector::ReSize(const GeneralMatrix& A)
00318 { REPORT ReSize(A.Nrows(), A.Ncols()); }
00319
00320 void SymmetricMatrix::ReSize(const GeneralMatrix& A)
00321 {
00322 REPORT
00323 int n = A.Nrows();
00324 if (n != A.Ncols())
00325 {
00326 Tracer tr("SymmetricMatrix::ReSize(GM)");
00327 Throw(NotSquareException(*this));
00328 }
00329 ReSize(n);
00330 }
00331
00332 void DiagonalMatrix::ReSize(const GeneralMatrix& A)
00333 {
00334 REPORT
00335 int n = A.Nrows();
00336 if (n != A.Ncols())
00337 {
00338 Tracer tr("DiagonalMatrix::ReSize(GM)");
00339 Throw(NotSquareException(*this));
00340 }
00341 ReSize(n);
00342 }
00343
00344 void UpperTriangularMatrix::ReSize(const GeneralMatrix& A)
00345 {
00346 REPORT
00347 int n = A.Nrows();
00348 if (n != A.Ncols())
00349 {
00350 Tracer tr("UpperTriangularMatrix::ReSize(GM)");
00351 Throw(NotSquareException(*this));
00352 }
00353 ReSize(n);
00354 }
00355
00356 void LowerTriangularMatrix::ReSize(const GeneralMatrix& A)
00357 {
00358 REPORT
00359 int n = A.Nrows();
00360 if (n != A.Ncols())
00361 {
00362 Tracer tr("LowerTriangularMatrix::ReSize(GM)");
00363 Throw(NotSquareException(*this));
00364 }
00365 ReSize(n);
00366 }
00367
00368 void IdentityMatrix::ReSize(const GeneralMatrix& A)
00369 {
00370 REPORT
00371 int n = A.Nrows();
00372 if (n != A.Ncols())
00373 {
00374 Tracer tr("IdentityMatrix::ReSize(GM)");
00375 Throw(NotSquareException(*this));
00376 }
00377 ReSize(n);
00378 }
00379
00380 void GeneralMatrix::ReSize(const GeneralMatrix&)
00381 {
00382 REPORT
00383 Tracer tr("GeneralMatrix::ReSize(GM)");
00384 Throw(NotDefinedException("ReSize", "this type of matrix"));
00385 }
00386
00387 void GeneralMatrix::ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix&)
00388 { REPORT ReSize(A); }
00389
00390 void GeneralMatrix::ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix&)
00391 { REPORT ReSize(A); }
00392
00393
00394
00395
00396
00397
00398
00399 bool GeneralMatrix::SameStorageType(const GeneralMatrix& A) const
00400 {
00401 REPORT
00402 return Type() == A.Type();
00403 }
00404
00405
00406
00407
00408 int GeneralMatrix::search(const BaseMatrix* s) const
00409 { REPORT return (s==this) ? 1 : 0; }
00410
00411 int GenericMatrix::search(const BaseMatrix* s) const
00412 { REPORT return gm->search(s); }
00413
00414 int MultipliedMatrix::search(const BaseMatrix* s) const
00415 { REPORT return bm1->search(s) + bm2->search(s); }
00416
00417 int ShiftedMatrix::search(const BaseMatrix* s) const
00418 { REPORT return bm->search(s); }
00419
00420 int NegatedMatrix::search(const BaseMatrix* s) const
00421 { REPORT return bm->search(s); }
00422
00423 int ReturnMatrix::search(const BaseMatrix* s) const
00424 { REPORT return (s==gm) ? 1 : 0; }
00425
00426 MatrixType Matrix::Type() const { return MatrixType::Rt; }
00427 MatrixType SquareMatrix::Type() const { return MatrixType::Sq; }
00428 MatrixType SymmetricMatrix::Type() const { return MatrixType::Sm; }
00429 MatrixType UpperTriangularMatrix::Type() const { return MatrixType::UT; }
00430 MatrixType LowerTriangularMatrix::Type() const { return MatrixType::LT; }
00431 MatrixType DiagonalMatrix::Type() const { return MatrixType::Dg; }
00432 MatrixType RowVector::Type() const { return MatrixType::RV; }
00433 MatrixType ColumnVector::Type() const { return MatrixType::CV; }
00434 MatrixType CroutMatrix::Type() const { return MatrixType::Ct; }
00435 MatrixType BandMatrix::Type() const { return MatrixType::BM; }
00436 MatrixType UpperBandMatrix::Type() const { return MatrixType::UB; }
00437 MatrixType LowerBandMatrix::Type() const { return MatrixType::LB; }
00438 MatrixType SymmetricBandMatrix::Type() const { return MatrixType::SB; }
00439
00440 MatrixType IdentityMatrix::Type() const { return MatrixType::Id; }
00441
00442
00443
00444 MatrixBandWidth BaseMatrix::BandWidth() const { REPORT return -1; }
00445 MatrixBandWidth DiagonalMatrix::BandWidth() const { REPORT return 0; }
00446 MatrixBandWidth IdentityMatrix::BandWidth() const { REPORT return 0; }
00447
00448 MatrixBandWidth UpperTriangularMatrix::BandWidth() const
00449 { REPORT return MatrixBandWidth(0,-1); }
00450
00451 MatrixBandWidth LowerTriangularMatrix::BandWidth() const
00452 { REPORT return MatrixBandWidth(-1,0); }
00453
00454 MatrixBandWidth BandMatrix::BandWidth() const
00455 { REPORT return MatrixBandWidth(lower,upper); }
00456
00457 MatrixBandWidth GenericMatrix::BandWidth()const
00458 { REPORT return gm->BandWidth(); }
00459
00460 MatrixBandWidth AddedMatrix::BandWidth() const
00461 { REPORT return gm1->BandWidth() + gm2->BandWidth(); }
00462
00463 MatrixBandWidth SPMatrix::BandWidth() const
00464 { REPORT return gm1->BandWidth().minimum(gm2->BandWidth()); }
00465
00466 MatrixBandWidth KPMatrix::BandWidth() const
00467 {
00468 int lower, upper;
00469 MatrixBandWidth bw1 = gm1->BandWidth(), bw2 = gm2->BandWidth();
00470 if (bw1.Lower() < 0)
00471 {
00472 if (bw2.Lower() < 0) { REPORT lower = -1; }
00473 else { REPORT lower = bw2.Lower() + (gm1->Nrows() - 1) * gm2->Nrows(); }
00474 }
00475 else
00476 {
00477 if (bw2.Lower() < 0)
00478 { REPORT lower = (1 + bw1.Lower()) * gm2->Nrows() - 1; }
00479 else { REPORT lower = bw2.Lower() + bw1.Lower() * gm2->Nrows(); }
00480 }
00481 if (bw1.Upper() < 0)
00482 {
00483 if (bw2.Upper() < 0) { REPORT upper = -1; }
00484 else { REPORT upper = bw2.Upper() + (gm1->Nrows() - 1) * gm2->Nrows(); }
00485 }
00486 else
00487 {
00488 if (bw2.Upper() < 0)
00489 { REPORT upper = (1 + bw1.Upper()) * gm2->Nrows() - 1; }
00490 else { REPORT upper = bw2.Upper() + bw1.Upper() * gm2->Nrows(); }
00491 }
00492 return MatrixBandWidth(lower, upper);
00493 }
00494
00495 MatrixBandWidth MultipliedMatrix::BandWidth() const
00496 { REPORT return gm1->BandWidth() * gm2->BandWidth(); }
00497
00498 MatrixBandWidth ConcatenatedMatrix::BandWidth() const { REPORT return -1; }
00499
00500 MatrixBandWidth SolvedMatrix::BandWidth() const
00501 {
00502 if (+gm1->Type() & MatrixType::Diagonal)
00503 { REPORT return gm2->BandWidth(); }
00504 else { REPORT return -1; }
00505 }
00506
00507 MatrixBandWidth ScaledMatrix::BandWidth() const
00508 { REPORT return gm->BandWidth(); }
00509
00510 MatrixBandWidth NegatedMatrix::BandWidth() const
00511 { REPORT return gm->BandWidth(); }
00512
00513 MatrixBandWidth TransposedMatrix::BandWidth() const
00514 { REPORT return gm->BandWidth().t(); }
00515
00516 MatrixBandWidth InvertedMatrix::BandWidth() const
00517 {
00518 if (+gm->Type() & MatrixType::Diagonal)
00519 { REPORT return MatrixBandWidth(0,0); }
00520 else { REPORT return -1; }
00521 }
00522
00523 MatrixBandWidth RowedMatrix::BandWidth() const { REPORT return -1; }
00524 MatrixBandWidth ColedMatrix::BandWidth() const { REPORT return -1; }
00525 MatrixBandWidth DiagedMatrix::BandWidth() const { REPORT return 0; }
00526 MatrixBandWidth MatedMatrix::BandWidth() const { REPORT return -1; }
00527 MatrixBandWidth ReturnMatrix::BandWidth() const
00528 { REPORT return gm->BandWidth(); }
00529
00530 MatrixBandWidth GetSubMatrix::BandWidth() const
00531 {
00532
00533 if (row_skip==col_skip && row_number==col_number)
00534 { REPORT return gm->BandWidth(); }
00535 else { REPORT return MatrixBandWidth(-1); }
00536 }
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548 void GeneralMatrix::tDelete()
00549 {
00550 if (tag<0)
00551 {
00552 if (tag<-1) { REPORT store = 0; delete this; return; }
00553 else { REPORT return; }
00554 }
00555 if (tag==1)
00556 {
00557 if (store)
00558 {
00559 REPORT MONITOR_REAL_DELETE("Free (tDelete)",storage,store)
00560 delete [] store;
00561 }
00562 MiniCleanUp(); return;
00563 }
00564 if (tag==0) { REPORT delete this; return; }
00565
00566 REPORT tag--; return;
00567 }
00568
00569 static void BlockCopy(int n, Real* from, Real* to)
00570 {
00571 REPORT
00572 int i = (n >> 3);
00573 while (i--)
00574 {
00575 *to++ = *from++; *to++ = *from++; *to++ = *from++; *to++ = *from++;
00576 *to++ = *from++; *to++ = *from++; *to++ = *from++; *to++ = *from++;
00577 }
00578 i = n & 7; while (i--) *to++ = *from++;
00579 }
00580
00581 bool GeneralMatrix::reuse()
00582 {
00583 if (tag < -1)
00584 {
00585 if (storage)
00586 {
00587 REPORT
00588 Real* s = new Real [storage]; MatrixErrorNoSpace(s);
00589 MONITOR_REAL_NEW("Make (reuse)",storage,s)
00590 BlockCopy(storage, store, s); store = s;
00591 }
00592 else { REPORT MiniCleanUp(); }
00593 tag = 0; return true;
00594 }
00595 if (tag < 0 ) { REPORT return false; }
00596 if (tag <= 1 ) { REPORT return true; }
00597 REPORT tag--; return false;
00598 }
00599
00600 Real* GeneralMatrix::GetStore()
00601 {
00602 if (tag<0 || tag>1)
00603 {
00604 Real* s;
00605 if (storage)
00606 {
00607 s = new Real [storage]; MatrixErrorNoSpace(s);
00608 MONITOR_REAL_NEW("Make (GetStore)",storage,s)
00609 BlockCopy(storage, store, s);
00610 }
00611 else s = 0;
00612 if (tag > 1) { REPORT tag--; }
00613 else if (tag < -1) { REPORT store = 0; delete this; }
00614 else { REPORT }
00615 return s;
00616 }
00617 Real* s = store;
00618 if (tag==0) { REPORT store = 0; delete this; }
00619 else { REPORT MiniCleanUp(); }
00620 return s;
00621 }
00622
00623 void GeneralMatrix::GetMatrix(const GeneralMatrix* gmx)
00624 {
00625 REPORT tag=-1; nrows_value=gmx->Nrows(); ncols_value=gmx->Ncols();
00626 storage=gmx->storage; SetParameters(gmx);
00627 store=(const_cast<GeneralMatrix*>(gmx))->GetStore();
00628 }
00629
00630 GeneralMatrix* GeneralMatrix::BorrowStore(GeneralMatrix* gmx, MatrixType mt)
00631
00632
00633 {
00634 if (!mt)
00635 {
00636 if (tag == -1) { REPORT gmx->tag = -2; gmx->store = store; }
00637 else { REPORT gmx->tag = 0; gmx->store = GetStore(); }
00638 }
00639 else if (Compare(gmx->Type(),mt))
00640 { REPORT gmx->tag = 0; gmx->store = GetStore(); }
00641 else
00642 {
00643 REPORT gmx->tag = -2; gmx->store = store;
00644 gmx = gmx->Evaluate(mt); gmx->tag = 0; tDelete();
00645 }
00646
00647 return gmx;
00648 }
00649
00650 void GeneralMatrix::Eq(const BaseMatrix& X, MatrixType mt)
00651
00652
00653
00654
00655 {
00656 #ifdef DO_SEARCH
00657 int counter=X.search(this);
00658 if (counter==0)
00659 {
00660 REPORT
00661 if (store)
00662 {
00663 MONITOR_REAL_DELETE("Free (operator=)",storage,store)
00664 REPORT delete [] store; storage = 0; store = 0;
00665 }
00666 }
00667 else { REPORT Release(counter); }
00668 GeneralMatrix* gmx = (const_cast<BaseMatrix&>(X)).Evaluate(mt);
00669 if (gmx!=this) { REPORT GetMatrix(gmx); }
00670 else { REPORT }
00671 Protect();
00672 #else
00673 GeneralMatrix* gmx = (const_cast<BaseMatrix&>(X)).Evaluate(mt);
00674 if (gmx!=this)
00675 {
00676 REPORT
00677 if (store)
00678 {
00679 MONITOR_REAL_DELETE("Free (operator=)",storage,store)
00680 REPORT delete [] store; storage = 0; store = 0;
00681 }
00682 GetMatrix(gmx);
00683 }
00684 else { REPORT }
00685 Protect();
00686 #endif
00687 }
00688
00689
00690 void GeneralMatrix::Eq(const GeneralMatrix& X)
00691 {
00692 GeneralMatrix* gmx = const_cast<GeneralMatrix*>(&X);
00693 if (gmx!=this)
00694 {
00695 REPORT
00696 if (store)
00697 {
00698 MONITOR_REAL_DELETE("Free (operator=)",storage,store)
00699 REPORT delete [] store; storage = 0; store = 0;
00700 }
00701 GetMatrix(gmx);
00702 }
00703 else { REPORT }
00704 Protect();
00705 }
00706
00707
00708 void GeneralMatrix::Eq(const BaseMatrix& X, MatrixType mt, bool ldok)
00709 {
00710 REPORT
00711 if (ldok) mt.SetDataLossOK();
00712 Eq(X, mt);
00713 }
00714
00715 void GeneralMatrix::Eq2(const BaseMatrix& X, MatrixType mt)
00716
00717
00718
00719
00720 {
00721 GeneralMatrix* gmx = (const_cast<BaseMatrix&>(X)).Evaluate(mt);
00722 if (gmx!=this) { REPORT GetMatrix(gmx); }
00723 else { REPORT }
00724 Protect();
00725 }
00726
00727 void GeneralMatrix::Inject(const GeneralMatrix& X)
00728
00729 {
00730 REPORT
00731 Tracer tr("Inject");
00732 if (nrows_value != X.nrows_value || ncols_value != X.ncols_value)
00733 Throw(IncompatibleDimensionsException());
00734 MatrixRow mr(const_cast<GeneralMatrix*>(&X), LoadOnEntry);
00735 MatrixRow mrx(this, LoadOnEntry+StoreOnExit+DirectPart);
00736 int i=nrows_value;
00737 while (i--) { mrx.Inject(mr); mrx.Next(); mr.Next(); }
00738 }
00739
00740
00741
00742 bool Compare(const MatrixType& source, MatrixType& destination)
00743 {
00744 if (!destination) { destination=source; return true; }
00745 if (destination==source) return true;
00746 if (!destination.DataLossOK && !(destination>=source))
00747 Throw(ProgramException("Illegal Conversion", source, destination));
00748 return false;
00749 }
00750
00751
00752
00753 GeneralMatrix* Matrix::Image() const
00754 {
00755 REPORT
00756 GeneralMatrix* gm = new Matrix(*this); MatrixErrorNoSpace(gm);
00757 return gm;
00758 }
00759
00760 GeneralMatrix* SquareMatrix::Image() const
00761 {
00762 REPORT
00763 GeneralMatrix* gm = new SquareMatrix(*this); MatrixErrorNoSpace(gm);
00764 return gm;
00765 }
00766
00767 GeneralMatrix* SymmetricMatrix::Image() const
00768 {
00769 REPORT
00770 GeneralMatrix* gm = new SymmetricMatrix(*this); MatrixErrorNoSpace(gm);
00771 return gm;
00772 }
00773
00774 GeneralMatrix* UpperTriangularMatrix::Image() const
00775 {
00776 REPORT
00777 GeneralMatrix* gm = new UpperTriangularMatrix(*this);
00778 MatrixErrorNoSpace(gm); return gm;
00779 }
00780
00781 GeneralMatrix* LowerTriangularMatrix::Image() const
00782 {
00783 REPORT
00784 GeneralMatrix* gm = new LowerTriangularMatrix(*this);
00785 MatrixErrorNoSpace(gm); return gm;
00786 }
00787
00788 GeneralMatrix* DiagonalMatrix::Image() const
00789 {
00790 REPORT
00791 GeneralMatrix* gm = new DiagonalMatrix(*this); MatrixErrorNoSpace(gm);
00792 return gm;
00793 }
00794
00795 GeneralMatrix* RowVector::Image() const
00796 {
00797 REPORT
00798 GeneralMatrix* gm = new RowVector(*this); MatrixErrorNoSpace(gm);
00799 return gm;
00800 }
00801
00802 GeneralMatrix* ColumnVector::Image() const
00803 {
00804 REPORT
00805 GeneralMatrix* gm = new ColumnVector(*this); MatrixErrorNoSpace(gm);
00806 return gm;
00807 }
00808
00809 GeneralMatrix* BandMatrix::Image() const
00810 {
00811 REPORT
00812 GeneralMatrix* gm = new BandMatrix(*this); MatrixErrorNoSpace(gm);
00813 return gm;
00814 }
00815
00816 GeneralMatrix* UpperBandMatrix::Image() const
00817 {
00818 REPORT
00819 GeneralMatrix* gm = new UpperBandMatrix(*this); MatrixErrorNoSpace(gm);
00820 return gm;
00821 }
00822
00823 GeneralMatrix* LowerBandMatrix::Image() const
00824 {
00825 REPORT
00826 GeneralMatrix* gm = new LowerBandMatrix(*this); MatrixErrorNoSpace(gm);
00827 return gm;
00828 }
00829
00830 GeneralMatrix* SymmetricBandMatrix::Image() const
00831 {
00832 REPORT
00833 GeneralMatrix* gm = new SymmetricBandMatrix(*this); MatrixErrorNoSpace(gm);
00834 return gm;
00835 }
00836
00837 GeneralMatrix* nricMatrix::Image() const
00838 {
00839 REPORT
00840 GeneralMatrix* gm = new nricMatrix(*this); MatrixErrorNoSpace(gm);
00841 return gm;
00842 }
00843
00844 GeneralMatrix* IdentityMatrix::Image() const
00845 {
00846 REPORT
00847 GeneralMatrix* gm = new IdentityMatrix(*this); MatrixErrorNoSpace(gm);
00848 return gm;
00849 }
00850
00851 GeneralMatrix* GeneralMatrix::Image() const
00852 {
00853 bool dummy = true;
00854 if (dummy)
00855 Throw(InternalException("Cannot apply Image to this matrix type"));
00856 return 0;
00857 }
00858
00859
00860
00861
00862 void nricMatrix::MakeRowPointer()
00863 {
00864 REPORT
00865 if (nrows_value > 0)
00866 {
00867 row_pointer = new Real* [nrows_value]; MatrixErrorNoSpace(row_pointer);
00868 Real* s = Store() - 1; int i = nrows_value; Real** rp = row_pointer;
00869 if (i) for (;;) { *rp++ = s; if (!(--i)) break; s+=ncols_value; }
00870 }
00871 else row_pointer = 0;
00872 }
00873
00874 void nricMatrix::DeleteRowPointer()
00875 { REPORT if (nrows_value) delete [] row_pointer; }
00876
00877 void GeneralMatrix::CheckStore() const
00878 {
00879 REPORT
00880 if (!store)
00881 Throw(ProgramException("NRIC accessing matrix with unset dimensions"));
00882 }
00883
00884
00885
00886
00887 void GeneralMatrix::CleanUp()
00888 {
00889
00890 REPORT
00891 if (store && storage)
00892 {
00893 MONITOR_REAL_DELETE("Free (CleanUp) ",storage,store)
00894 REPORT delete [] store;
00895 }
00896 store=0; storage=0; nrows_value=0; ncols_value=0; tag = -1;
00897 }
00898
00899 void nricMatrix::CleanUp()
00900 { REPORT DeleteRowPointer(); GeneralMatrix::CleanUp(); }
00901
00902 void nricMatrix::MiniCleanUp()
00903 { REPORT DeleteRowPointer(); GeneralMatrix::MiniCleanUp(); }
00904
00905 void RowVector::CleanUp()
00906 { REPORT GeneralMatrix::CleanUp(); nrows_value=1; }
00907
00908 void ColumnVector::CleanUp()
00909 { REPORT GeneralMatrix::CleanUp(); ncols_value=1; }
00910
00911 void CroutMatrix::CleanUp()
00912 {
00913 REPORT
00914 if (nrows_value) delete [] indx;
00915 GeneralMatrix::CleanUp();
00916 }
00917
00918 void CroutMatrix::MiniCleanUp()
00919 {
00920 REPORT
00921 if (nrows_value) delete [] indx;
00922 GeneralMatrix::MiniCleanUp();
00923 }
00924
00925 void BandLUMatrix::CleanUp()
00926 {
00927 REPORT
00928 if (nrows_value) delete [] indx;
00929 if (storage2) delete [] store2;
00930 GeneralMatrix::CleanUp();
00931 }
00932
00933 void BandLUMatrix::MiniCleanUp()
00934 {
00935 REPORT
00936 if (nrows_value) delete [] indx;
00937 if (storage2) delete [] store2;
00938 GeneralMatrix::MiniCleanUp();
00939 }
00940
00941
00942
00943
00944
00945
00946 SimpleIntArray::SimpleIntArray(int xn) : n(xn)
00947 {
00948 if (n < 0) Throw(Logic_error("invalid array length"));
00949 else if (n == 0) { REPORT a = 0; }
00950 else { REPORT a = new int [n]; if (!a) Throw(Bad_alloc()); }
00951 }
00952
00953
00954
00955 SimpleIntArray::~SimpleIntArray() { REPORT if (a) delete [] a; }
00956
00957
00958
00959
00960
00961
00962 int& SimpleIntArray::operator[](int i)
00963 {
00964 REPORT
00965 if (i < 0 || i >= n) Throw(Logic_error("array index out of range"));
00966 return a[i];
00967 }
00968
00969
00970
00971
00972 int SimpleIntArray::operator[](int i) const
00973 {
00974 REPORT
00975 if (i < 0 || i >= n) Throw(Logic_error("array index out of range"));
00976 return a[i];
00977 }
00978
00979
00980
00981 void SimpleIntArray::operator=(int ai)
00982 { REPORT for (int i = 0; i < n; i++) a[i] = ai; }
00983
00984
00985
00986
00987 void SimpleIntArray::operator=(const SimpleIntArray& b)
00988 {
00989 REPORT
00990 if (b.n != n) ReSize(b.n);
00991 for (int i = 0; i < n; i++) a[i] = b.a[i];
00992 }
00993
00994
00995
00996
00997 SimpleIntArray::SimpleIntArray(const SimpleIntArray& b) : Janitor(), n(b.n)
00998 {
00999 if (n == 0) { REPORT a = 0; }
01000 else
01001 {
01002 REPORT
01003 a = new int [n]; if (!a) Throw(Bad_alloc());
01004 for (int i = 0; i < n; i++) a[i] = b.a[i];
01005 }
01006 }
01007
01008
01009
01010
01011 void SimpleIntArray::ReSize(int n1, bool keep)
01012 {
01013 if (n1 == n) { REPORT return; }
01014 else if (n1 == 0) { REPORT n = 0; delete [] a; a = 0; }
01015 else if (n == 0)
01016 { REPORT a = new int [n1]; if (!a) Throw(Bad_alloc()); n = n1; }
01017 else
01018 {
01019 int* a1 = a;
01020 if (keep)
01021 {
01022 REPORT
01023 a = new int [n1]; if (!a) Throw(Bad_alloc());
01024 if (n > n1) n = n1;
01025 for (int i = 0; i < n; i++) a[i] = a1[i];
01026 n = n1; delete [] a1;
01027 }
01028 else
01029 {
01030 REPORT n = n1; delete [] a1;
01031 a = new int [n]; if (!a) Throw(Bad_alloc());
01032 }
01033 }
01034 }
01035
01036
01037
01038
01039
01040 void GeneralMatrix::swap(GeneralMatrix& gm)
01041 {
01042 REPORT
01043 int t;
01044 t = tag; tag = gm.tag; gm.tag = t;
01045 t = nrows_value; nrows_value = gm.nrows_value; gm.nrows_value = t;
01046 t = ncols_value; ncols_value = gm.ncols_value; gm.ncols_value = t;
01047 t = storage; storage = gm.storage; gm.storage = t;
01048 Real* s = store; store = gm.store; gm.store = s;
01049 }
01050
01051 void nricMatrix::swap(nricMatrix& gm)
01052 {
01053 REPORT
01054 GeneralMatrix::swap((GeneralMatrix&)gm);
01055 Real** rp = row_pointer; row_pointer = gm.row_pointer; gm.row_pointer = rp;
01056 }
01057
01058 void CroutMatrix::swap(CroutMatrix& gm)
01059 {
01060 REPORT
01061 GeneralMatrix::swap((GeneralMatrix&)gm);
01062 int* i = indx; indx = gm.indx; gm.indx = i;
01063 bool b;
01064 b = d; d = gm.d; gm.d = b;
01065 b = sing; sing = gm.sing; gm.sing = b;
01066 }
01067
01068 void BandMatrix::swap(BandMatrix& gm)
01069 {
01070 REPORT
01071 GeneralMatrix::swap((GeneralMatrix&)gm);
01072 int i;
01073 i = lower; lower = gm.lower; gm.lower = i;
01074 i = upper; upper = gm.upper; gm.upper = i;
01075 }
01076
01077 void SymmetricBandMatrix::swap(SymmetricBandMatrix& gm)
01078 {
01079 REPORT
01080 GeneralMatrix::swap((GeneralMatrix&)gm);
01081 int i;
01082 i = lower; lower = gm.lower; gm.lower = i;
01083 }
01084
01085 void BandLUMatrix::swap(BandLUMatrix& gm)
01086 {
01087 REPORT
01088 GeneralMatrix::swap((GeneralMatrix&)gm);
01089 int* i = indx; indx = gm.indx; gm.indx = i;
01090 bool b;
01091 b = d; d = gm.d; gm.d = b;
01092 b = sing; sing = gm.sing; gm.sing = b;
01093 int m;
01094 m = storage2; storage2 = gm.storage2; gm.storage2 = m;
01095 m = m1; m1 = gm.m1; gm.m1 = m;
01096 m = m2; m2 = gm.m2; gm.m2 = m;
01097 Real* s = store2; store2 = gm.store2; gm.store2 = s;
01098 }
01099
01100 void GenericMatrix::swap(GenericMatrix& x)
01101 {
01102 REPORT
01103 GeneralMatrix* tgm = gm; gm = x.gm; x.gm = tgm;
01104 }
01105
01106
01107
01108 RealStarStar::RealStarStar(Matrix& A)
01109 {
01110 REPORT
01111 Tracer tr("RealStarStar");
01112 int n = A.ncols();
01113 int m = A.nrows();
01114 a = new Real*[m];
01115 MatrixErrorNoSpace(a);
01116 Real* d = A.data();
01117 for (int i = 0; i < m; ++i) a[i] = d + i * n;
01118 }
01119
01120 ConstRealStarStar::ConstRealStarStar(const Matrix& A)
01121 {
01122 REPORT
01123 Tracer tr("ConstRealStarStar");
01124 int n = A.ncols();
01125 int m = A.nrows();
01126 a = new const Real*[m];
01127 MatrixErrorNoSpace(a);
01128 const Real* d = A.data();
01129 for (int i = 0; i < m; ++i) a[i] = d + i * n;
01130 }
01131
01132
01133
01134 #ifdef use_namespace
01135 }
01136 #endif
01137