newmat4.cpp
Go to the documentation of this file.
00001 
00002 
00003 
00006 
00007 
00008 // Copyright (C) 1991,2,3,4,8,9: R B Davies
00009 
00010 //#define WANT_STREAM
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 // ************************* general utilities *************************/
00033 
00034 static int tristore(int n)                    // elements in triangular matrix
00035 { return (n*(n+1))/2; }
00036 
00037 
00038 // **************************** constructors ***************************/
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 // CheckConversion(M);
00080    // MatrixConversionCheck mcc;
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  // CheckConversion(M);
00131    // MatrixConversionCheck mcc;
00132    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Sm);
00133    GetMatrix(gmx);
00134 }
00135 
00136 UpperTriangularMatrix::UpperTriangularMatrix(const BaseMatrix& M)
00137 {
00138    REPORT // CheckConversion(M);
00139    // MatrixConversionCheck mcc;
00140    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::UT);
00141    GetMatrix(gmx);
00142 }
00143 
00144 LowerTriangularMatrix::LowerTriangularMatrix(const BaseMatrix& M)
00145 {
00146    REPORT // CheckConversion(M);
00147    // MatrixConversionCheck mcc;
00148    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::LT);
00149    GetMatrix(gmx);
00150 }
00151 
00152 DiagonalMatrix::DiagonalMatrix(const BaseMatrix& M)
00153 {
00154    REPORT //CheckConversion(M);
00155    // MatrixConversionCheck mcc;
00156    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Dg);
00157    GetMatrix(gmx);
00158 }
00159 
00160 IdentityMatrix::IdentityMatrix(const BaseMatrix& M)
00161 {
00162    REPORT //CheckConversion(M);
00163    // MatrixConversionCheck mcc;
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;                     // in case of exception at next line
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 // could we use SetParameters instead of this
00200 void CroutMatrix::get_aux(CroutMatrix& X)
00201 {
00202    X.d = d; X.sing = sing;
00203    if (tag_val == 0 || tag_val == 1) // reuse the array 
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                              // copy the array
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 //ReturnMatrix::ReturnMatrix(GeneralMatrix& gmx)
00234 //{
00235 //   REPORT
00236 //   gm = gmx.Image(); gm->ReleaseAndDelete();
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 // ************ Constructors for use with NR in C++ interface ***********
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 // ************************** resize matrices ***************************/
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 //*********************** resize_keep *******************************
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 void GeneralMatrix::resizeForAdd(const GeneralMatrix& A, const GeneralMatrix&)
00615 { REPORT resize(A); }
00616 
00617 void GeneralMatrix::resizeForSP(const GeneralMatrix& A, const GeneralMatrix&)
00618 { REPORT resize(A); }
00619 
00620 
00621 // ************************* SameStorageType ******************************
00622 
00623 // SameStorageType checks A and B have same storage type including bandwidth
00624 // It does not check same dimensions since we assume this is already done
00625 
00626 bool GeneralMatrix::SameStorageType(const GeneralMatrix& A) const
00627 {
00628    REPORT
00629    return type() == A.type();
00630 }
00631 */
00632 
00633 // ******************* manipulate types, storage **************************/
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 // ********************** the memory managment tools **********************/
00769 
00770 //  Rules regarding tDelete, reuse, GetStore, BorrowStore
00771 //    All matrices processed during expression evaluation must be subject
00772 //    to exactly one of reuse(), tDelete(), GetStore() or BorrowStore().
00773 //    If reuse returns true the matrix must be reused.
00774 //    GetMatrix(gm) always calls gm->GetStore()
00775 //    gm->Evaluate obeys rules
00776 //    bm->Evaluate obeys rules for matrices in bm structure
00777 
00778 //  Meaning of tag_val
00779 //    tag_val = -1          memory cannot be reused (default situation)
00780 //    tag_val = -2          memory has been borrowed from another matrix
00781 //                               (don't change values)
00782 //    tag_val = i > 0       delete or reuse memory after i operations
00783 //    tag_val = 0           like value 1 but matrix was created by new
00784 //                               so delete it
00785 
00786 void GeneralMatrix::tDelete()
00787 {
00788    if (tag_val<0)
00789    {
00790       if (tag_val<-1) { REPORT store = 0; delete this; return; }  // borrowed
00791       else { REPORT return; }   // not a temporary matrix - leave alone
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;                           // CleanUp
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)                 // borrowed storage
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(); }                // CleanUp
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; } // borrowed store
00852       else { REPORT }
00853       return s;
00854    }
00855    Real* s = store;                             // cleanup - done later
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 // Copy storage of *this to storage of *gmx. Then convert to type mt.
00870 // If mt == 0 just let *gmx point to storage of *this if tag_val==-1.
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 // Count number of references to this in X.
00890 // If zero delete storage in this;
00891 // otherwise tag this to show when storage can be deleted
00892 // evaluate X and copy to this
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 // version with no conversion
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 // version to work with operator<<
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 // a cut down version of Eq for use with += etc.
00955 // we know BaseMatrix points to two GeneralMatrix objects,
00956 // the first being this (may be the same).
00957 // we know tag_val has been set correctly in each.
00958 {
00959    GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate(mt);
00960    if (gmx!=this) { REPORT GetMatrix(gmx); }  // simplify GetMatrix ?
00961    else { REPORT }
00962    Protect();
00963 }
00964 
00965 void GeneralMatrix::inject(const GeneralMatrix& X)
00966 // copy stored values of X; otherwise leave els of *this unchanged
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 // ************* checking for data loss during conversion *******************/
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 // ************* Make a copy of a matrix on the heap *********************/
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)                                   // get rid of warning message
01072       Throw(InternalException("Cannot apply Image to this matrix type"));
01073    return 0;
01074 }
01075 
01076 
01077 // *********************** nricMatrix routines *****************************/
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 // *************************** CleanUp routines *****************************/
01103 
01104 void GeneralMatrix::cleanup()
01105 {
01106    // set matrix dimensions to zero, delete storage
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 // ************************ simple integer array class ***********************
01159 
01160 // construct a new array of length xn. Check that xn is non-negative and
01161 // that space is available
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 // destroy an array - return its space to memory
01171 
01172 SimpleIntArray::~SimpleIntArray() { REPORT  if (a) delete [] a; }
01173 
01174 // access an element of an array; return a "reference" so elements
01175 // can be modified.
01176 // check index is within range
01177 // in this array class the index runs from 0 to n-1
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 // same thing again but for arrays declared constant so we can't
01187 // modify its elements
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 // set all the elements equal to a given value
01197 
01198 void SimpleIntArray::operator=(int ai)
01199    { REPORT  for (int i = 0; i < n; i++) a[i] = ai; }
01200 
01201 // set the elements equal to those of another array.
01202 // now allow length to be changed
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 // construct a new array equal to an existing array
01212 // check that space is available
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 // change the size of an array; optionally copy data from old array to
01226 // new array
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 //************************** swap values ********************************
01261 
01262 // swap values
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 // ********************** C subscript classes ****************************
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 


kni
Author(s): Martin Günther
autogenerated on Mon Aug 14 2017 02:44:12