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