Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

newmat4.cpp

Go to the documentation of this file.
00001 //$$ newmat4.cpp       Constructors, ReSize, basic utilities
00002 
00003 // Copyright (C) 1991,2,3,4,8,9: R B Davies
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 // ************************* general utilities *************************/
00026 
00027 static int tristore(int n)                    // elements in triangular matrix
00028 { return (n*(n+1))/2; }
00029 
00030 
00031 // **************************** constructors ***************************/
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 // CheckConversion(M);
00073    // MatrixConversionCheck mcc;
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  // CheckConversion(M);
00124    // MatrixConversionCheck mcc;
00125    GeneralMatrix* gmx=(const_cast<BaseMatrix&>(M)).Evaluate(MatrixType::Sm);
00126    GetMatrix(gmx);
00127 }
00128 
00129 UpperTriangularMatrix::UpperTriangularMatrix(const BaseMatrix& M)
00130 {
00131    REPORT // CheckConversion(M);
00132    // MatrixConversionCheck mcc;
00133    GeneralMatrix* gmx=(const_cast<BaseMatrix&>(M)).Evaluate(MatrixType::UT);
00134    GetMatrix(gmx);
00135 }
00136 
00137 LowerTriangularMatrix::LowerTriangularMatrix(const BaseMatrix& M)
00138 {
00139    REPORT // CheckConversion(M);
00140    // MatrixConversionCheck mcc;
00141    GeneralMatrix* gmx=(const_cast<BaseMatrix&>(M)).Evaluate(MatrixType::LT);
00142    GetMatrix(gmx);
00143 }
00144 
00145 DiagonalMatrix::DiagonalMatrix(const BaseMatrix& M)
00146 {
00147    REPORT //CheckConversion(M);
00148    // MatrixConversionCheck mcc;
00149    GeneralMatrix* gmx=(const_cast<BaseMatrix&>(M)).Evaluate(MatrixType::Dg);
00150    GetMatrix(gmx);
00151 }
00152 
00153 IdentityMatrix::IdentityMatrix(const BaseMatrix& M)
00154 {
00155    REPORT //CheckConversion(M);
00156    // MatrixConversionCheck mcc;
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;                     // in case of exception at next line
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 //ReturnMatrix::ReturnMatrix(GeneralMatrix& gmx)
00191 //{
00192 //   REPORT
00193 //   gm = gmx.Image(); gm->ReleaseAndDelete();
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 // ************ Constructors for use with NR in C++ interface ***********
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 // ************************** ReSize matrices ***************************/
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 // ************************* SameStorageType ******************************/
00395 
00396 // SameStorageType checks A and B have same storage type including bandwidth
00397 // It does not check same dimensions since we assume this is already done
00398 
00399 bool GeneralMatrix::SameStorageType(const GeneralMatrix& A) const
00400 {
00401    REPORT
00402    return Type() == A.Type();
00403 }
00404 
00405 
00406 // ******************* manipulate types, storage **************************/
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 // ********************** the memory managment tools **********************/
00539 
00540 //  Rules regarding tDelete, reuse, GetStore, BorrowStore
00541 //    All matrices processed during expression evaluation must be subject
00542 //    to exactly one of reuse(), tDelete(), GetStore() or BorrowStore().
00543 //    If reuse returns true the matrix must be reused.
00544 //    GetMatrix(gm) always calls gm->GetStore()
00545 //    gm->Evaluate obeys rules
00546 //    bm->Evaluate obeys rules for matrices in bm structure
00547 
00548 void GeneralMatrix::tDelete()
00549 {
00550    if (tag<0)
00551    {
00552       if (tag<-1) { REPORT store = 0; delete this; return; }  // borrowed
00553       else { REPORT return; }   // not a temporary matrix - leave alone
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;                           // CleanUp
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)                 // borrowed storage
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(); }                // CleanUp
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; } // borrowed store
00614       else { REPORT }
00615       return s;
00616    }
00617    Real* s = store;                             // CleanUp - done later
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 // Copy storage of *this to storage of *gmx. Then convert to type mt.
00632 // If mt == 0 just let *gmx point to storage of *this if tag==-1.
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 // Count number of references to this in X.
00652 // If zero delete storage in this;
00653 // otherwise tag this to show when storage can be deleted
00654 // evaluate X and copy to this
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 // version with no conversion
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 // version to work with operator<<
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 // a cut down version of Eq for use with += etc.
00717 // we know BaseMatrix points to two GeneralMatrix objects,
00718 // the first being this (may be the same).
00719 // we know tag has been set correctly in each.
00720 {
00721    GeneralMatrix* gmx = (const_cast<BaseMatrix&>(X)).Evaluate(mt);
00722    if (gmx!=this) { REPORT GetMatrix(gmx); }  // simplify GetMatrix ?
00723    else { REPORT }
00724    Protect();
00725 }
00726 
00727 void GeneralMatrix::Inject(const GeneralMatrix& X)
00728 // copy stored values of X; otherwise leave els of *this unchanged
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 // ************* checking for data loss during conversion *******************/
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 // ************* Make a copy of a matrix on the heap *********************/
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)                                   // get rid of warning message
00855       Throw(InternalException("Cannot apply Image to this matrix type"));
00856    return 0;
00857 }
00858 
00859 
00860 // *********************** nricMatrix routines *****************************/
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 // *************************** CleanUp routines *****************************/
00886 
00887 void GeneralMatrix::CleanUp()
00888 {
00889    // set matrix dimensions to zero, delete storage
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 // ************************ simple integer array class ***********************
00942 
00943 // construct a new array of length xn. Check that xn is non-negative and
00944 // that space is available
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 // destroy an array - return its space to memory
00954 
00955 SimpleIntArray::~SimpleIntArray() { REPORT  if (a) delete [] a; }
00956 
00957 // access an element of an array; return a "reference" so elements
00958 // can be modified.
00959 // check index is within range
00960 // in this array class the index runs from 0 to n-1
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 // same thing again but for arrays declared constant so we can't
00970 // modify its elements
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 // set all the elements equal to a given value
00980 
00981 void SimpleIntArray::operator=(int ai)
00982    { REPORT  for (int i = 0; i < n; i++) a[i] = ai; }
00983 
00984 // set the elements equal to those of another array.
00985 // now allow length to be changed
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 // construct a new array equal to an existing array
00995 // check that space is available
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 // change the size of an array; optionally copy data from old array to
01009 // new array
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 //************************** swap values ********************************
01037 
01038 // swap values
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 // ********************** C subscript classes ****************************
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 

newmat11b
Generated Mon May 9 04:54:18 2016 by Doxygen 1.6.3