newmat5.cpp
Go to the documentation of this file.
00001 
00002 
00003 
00006 
00007 // Copyright (C) 1991,2,3,4: R B Davies
00008 
00009 //#define WANT_STREAM
00010 
00011 #include <iostream>
00012 
00013 #include "include.h"
00014 
00015 #include "newmat.h"
00016 #include "newmatrc.h"
00017 
00018 #ifdef use_namespace
00019 namespace NEWMAT {
00020 #endif
00021 
00022 
00023 #ifdef DO_REPORT
00024 #define REPORT { static ExeCounter ExeCount(__LINE__,5); ++ExeCount; }
00025 #else
00026 #define REPORT {}
00027 #endif
00028 
00029 
00030 /************************ carry out operations ******************************/
00031 
00032 
00033 GeneralMatrix* GeneralMatrix::Transpose(TransposedMatrix* tm, MatrixType mt)
00034 {
00035    GeneralMatrix* gm1;
00036 
00037    if (Compare(Type().t(),mt))
00038    {
00039       REPORT
00040       gm1 = mt.New(ncols_val,nrows_val,tm);
00041       for (int i=0; i<ncols_val; i++)
00042       {
00043          MatrixRow mr(gm1, StoreOnExit+DirectPart, i);
00044          MatrixCol mc(this, mr.Data(), LoadOnEntry, i);
00045       }
00046    }
00047    else
00048    {
00049       REPORT
00050       gm1 = mt.New(ncols_val,nrows_val,tm);
00051       MatrixRow mr(this, LoadOnEntry);
00052       MatrixCol mc(gm1, StoreOnExit+DirectPart);
00053       int i = nrows_val;
00054       while (i--) { mc.Copy(mr); mr.Next(); mc.Next(); }
00055    }
00056    tDelete(); gm1->ReleaseAndDelete(); return gm1;
00057 }
00058 
00059 GeneralMatrix* SymmetricMatrix::Transpose(TransposedMatrix*, MatrixType mt)
00060 { REPORT  return Evaluate(mt); }
00061 
00062 
00063 GeneralMatrix* DiagonalMatrix::Transpose(TransposedMatrix*, MatrixType mt)
00064 { REPORT return Evaluate(mt); }
00065 
00066 GeneralMatrix* ColumnVector::Transpose(TransposedMatrix*, MatrixType mt)
00067 {
00068    REPORT
00069    GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
00070    gmx->nrows_val = 1; gmx->ncols_val = gmx->storage = storage;
00071    return BorrowStore(gmx,mt);
00072 }
00073 
00074 GeneralMatrix* RowVector::Transpose(TransposedMatrix*, MatrixType mt)
00075 {
00076    REPORT
00077    GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
00078    gmx->ncols_val = 1; gmx->nrows_val = gmx->storage = storage;
00079    return BorrowStore(gmx,mt);
00080 }
00081 
00082 GeneralMatrix* IdentityMatrix::Transpose(TransposedMatrix*, MatrixType mt)
00083 { REPORT return Evaluate(mt); }
00084 
00085 GeneralMatrix* GeneralMatrix::Evaluate(MatrixType mt)
00086 {
00087    if (Compare(this->Type(),mt)) { REPORT return this; }
00088    REPORT
00089    GeneralMatrix* gmx = mt.New(nrows_val,ncols_val,this);
00090    MatrixRow mr(this, LoadOnEntry);
00091    MatrixRow mrx(gmx, StoreOnExit+DirectPart);
00092    int i=nrows_val;
00093    while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
00094    tDelete(); gmx->ReleaseAndDelete(); return gmx;
00095 }
00096 
00097 GeneralMatrix* CroutMatrix::Evaluate(MatrixType mt)
00098 {
00099    if (Compare(this->Type(),mt)) { REPORT return this; }
00100    REPORT
00101    Tracer et("CroutMatrix::Evaluate");
00102    bool dummy = true;
00103    if (dummy) Throw(ProgramException("Illegal use of CroutMatrix", *this));
00104    return this;
00105 }
00106 
00107 GeneralMatrix* GenericMatrix::Evaluate(MatrixType mt)
00108    { REPORT  return gm->Evaluate(mt); }
00109 
00110 GeneralMatrix* ShiftedMatrix::Evaluate(MatrixType mt)
00111 {
00112    gm=((BaseMatrix*&)bm)->Evaluate();
00113    int nr=gm->Nrows(); int nc=gm->Ncols();
00114    Compare(gm->Type().AddEqualEl(),mt);
00115    if (!(mt==gm->Type()))
00116    {
00117       REPORT
00118       GeneralMatrix* gmx = mt.New(nr,nc,this);
00119       MatrixRow mr(gm, LoadOnEntry);
00120       MatrixRow mrx(gmx, StoreOnExit+DirectPart);
00121       while (nr--) { mrx.Add(mr,f); mrx.Next(); mr.Next(); }
00122       gmx->ReleaseAndDelete(); gm->tDelete();
00123       return gmx;
00124    }
00125    else if (gm->reuse())
00126    {
00127       REPORT gm->Add(f);
00128       return gm;
00129    }
00130    else
00131    {
00132       REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc,this);
00133       gmy->ReleaseAndDelete(); gmy->Add(gm,f);
00134       return gmy;
00135    }
00136 }
00137 
00138 GeneralMatrix* NegShiftedMatrix::Evaluate(MatrixType mt)
00139 {
00140    gm=((BaseMatrix*&)bm)->Evaluate();
00141    int nr=gm->Nrows(); int nc=gm->Ncols();
00142    Compare(gm->Type().AddEqualEl(),mt);
00143    if (!(mt==gm->Type()))
00144    {
00145       REPORT
00146       GeneralMatrix* gmx = mt.New(nr,nc,this);
00147       MatrixRow mr(gm, LoadOnEntry);
00148       MatrixRow mrx(gmx, StoreOnExit+DirectPart);
00149       while (nr--) { mrx.NegAdd(mr,f); mrx.Next(); mr.Next(); }
00150       gmx->ReleaseAndDelete(); gm->tDelete();
00151       return gmx;
00152    }
00153    else if (gm->reuse())
00154    {
00155       REPORT gm->NegAdd(f);
00156       return gm;
00157    }
00158    else
00159    {
00160       REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc,this);
00161       gmy->ReleaseAndDelete(); gmy->NegAdd(gm,f);
00162       return gmy;
00163    }
00164 }
00165 
00166 GeneralMatrix* ScaledMatrix::Evaluate(MatrixType mt)
00167 {
00168    gm=((BaseMatrix*&)bm)->Evaluate();
00169    int nr=gm->Nrows(); int nc=gm->Ncols();
00170    if (Compare(gm->Type(),mt))
00171    {
00172       if (gm->reuse())
00173       {
00174          REPORT gm->Multiply(f);
00175          return gm;
00176       }
00177       else
00178       {
00179          REPORT GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
00180          gmx->ReleaseAndDelete(); gmx->Multiply(gm,f);
00181          return gmx;
00182       }
00183    }
00184    else
00185    {
00186       REPORT
00187       GeneralMatrix* gmx = mt.New(nr,nc,this);
00188       MatrixRow mr(gm, LoadOnEntry);
00189       MatrixRow mrx(gmx, StoreOnExit+DirectPart);
00190       while (nr--) { mrx.Multiply(mr,f); mrx.Next(); mr.Next(); }
00191       gmx->ReleaseAndDelete(); gm->tDelete();
00192       return gmx;
00193    }
00194 }
00195 
00196 GeneralMatrix* NegatedMatrix::Evaluate(MatrixType mt)
00197 {
00198    gm=((BaseMatrix*&)bm)->Evaluate();
00199    int nr=gm->Nrows(); int nc=gm->Ncols();
00200    if (Compare(gm->Type(),mt))
00201    {
00202       if (gm->reuse())
00203       {
00204          REPORT gm->Negate();
00205          return gm;
00206       }
00207       else
00208       {
00209          REPORT
00210          GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
00211          gmx->ReleaseAndDelete(); gmx->Negate(gm);
00212          return gmx;
00213       }
00214    }
00215    else
00216    {
00217       REPORT
00218       GeneralMatrix* gmx = mt.New(nr,nc,this);
00219       MatrixRow mr(gm, LoadOnEntry);
00220       MatrixRow mrx(gmx, StoreOnExit+DirectPart);
00221       while (nr--) { mrx.Negate(mr); mrx.Next(); mr.Next(); }
00222       gmx->ReleaseAndDelete(); gm->tDelete();
00223       return gmx;
00224    }
00225 }
00226 
00227 GeneralMatrix* ReversedMatrix::Evaluate(MatrixType mt)
00228 {
00229    gm=((BaseMatrix*&)bm)->Evaluate(); GeneralMatrix* gmx;
00230 
00231    if ((gm->Type()).is_band() && ! (gm->Type()).is_diagonal())
00232    {
00233       gm->tDelete();
00234       Throw(NotDefinedException("Reverse", "band matrices"));
00235    }
00236 
00237    if (gm->reuse()) { REPORT gm->ReverseElements(); gmx = gm; }
00238    else
00239    {
00240       REPORT
00241       gmx = gm->Type().New(gm->Nrows(), gm->Ncols(), this);
00242       gmx->ReverseElements(gm); gmx->ReleaseAndDelete();
00243    }
00244    return gmx->Evaluate(mt); // target matrix is different type?
00245 
00246 }
00247 
00248 GeneralMatrix* TransposedMatrix::Evaluate(MatrixType mt)
00249 {
00250    REPORT
00251    gm=((BaseMatrix*&)bm)->Evaluate();
00252    Compare(gm->Type().t(),mt);
00253    GeneralMatrix* gmx=gm->Transpose(this, mt);
00254    return gmx;
00255 }
00256 
00257 GeneralMatrix* RowedMatrix::Evaluate(MatrixType mt)
00258 {
00259    gm = ((BaseMatrix*&)bm)->Evaluate();
00260    GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
00261    gmx->nrows_val = 1; gmx->ncols_val = gmx->storage = gm->storage;
00262    return gm->BorrowStore(gmx,mt);
00263 }
00264 
00265 GeneralMatrix* ColedMatrix::Evaluate(MatrixType mt)
00266 {
00267    gm = ((BaseMatrix*&)bm)->Evaluate();
00268    GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
00269    gmx->ncols_val = 1; gmx->nrows_val = gmx->storage = gm->storage;
00270    return gm->BorrowStore(gmx,mt);
00271 }
00272 
00273 GeneralMatrix* DiagedMatrix::Evaluate(MatrixType mt)
00274 {
00275    gm = ((BaseMatrix*&)bm)->Evaluate();
00276    GeneralMatrix* gmx = new DiagonalMatrix; MatrixErrorNoSpace(gmx);
00277    gmx->nrows_val = gmx->ncols_val = gmx->storage = gm->storage;
00278    return gm->BorrowStore(gmx,mt);
00279 }
00280 
00281 GeneralMatrix* MatedMatrix::Evaluate(MatrixType mt)
00282 {
00283    Tracer tr("MatedMatrix::Evaluate");
00284    gm = ((BaseMatrix*&)bm)->Evaluate();
00285    GeneralMatrix* gmx = new Matrix; MatrixErrorNoSpace(gmx);
00286    gmx->nrows_val = nr; gmx->ncols_val = nc; gmx->storage = gm->storage;
00287    if (nr*nc != gmx->storage)
00288       Throw(IncompatibleDimensionsException());
00289    return gm->BorrowStore(gmx,mt);
00290 }
00291 
00292 GeneralMatrix* GetSubMatrix::Evaluate(MatrixType mt)
00293 {
00294    REPORT
00295    Tracer tr("SubMatrix(evaluate)");
00296    gm = ((BaseMatrix*&)bm)->Evaluate();
00297    if (row_number < 0) row_number = gm->Nrows();
00298    if (col_number < 0) col_number = gm->Ncols();
00299    if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
00300    {
00301       gm->tDelete();
00302       Throw(SubMatrixDimensionException());
00303    }
00304    if (IsSym) Compare(gm->Type().ssub(), mt);
00305    else Compare(gm->Type().sub(), mt);
00306    GeneralMatrix* gmx = mt.New(row_number, col_number, this);
00307    int i = row_number;
00308    MatrixRow mr(gm, LoadOnEntry, row_skip); 
00309    MatrixRow mrx(gmx, StoreOnExit+DirectPart);
00310    MatrixRowCol sub;
00311    while (i--)
00312    {
00313       mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
00314       mrx.Copy(sub); mrx.Next(); mr.Next();
00315    }
00316    gmx->ReleaseAndDelete(); gm->tDelete();
00317    return gmx;
00318 }
00319 
00320 
00321 GeneralMatrix* ReturnMatrix::Evaluate(MatrixType mt)
00322 {
00323    return gm->Evaluate(mt);
00324 }
00325 
00326 
00327 void GeneralMatrix::Add(GeneralMatrix* gm1, Real f)
00328 {
00329    REPORT
00330    Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
00331    while (i--)
00332    { *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; }
00333    i = storage & 3; while (i--) *s++ = *s1++ + f;
00334 }
00335    
00336 void GeneralMatrix::Add(Real f)
00337 {
00338    REPORT
00339    Real* s=store; int i=(storage >> 2);
00340    while (i--) { *s++ += f; *s++ += f; *s++ += f; *s++ += f; }
00341    i = storage & 3; while (i--) *s++ += f;
00342 }
00343    
00344 void GeneralMatrix::NegAdd(GeneralMatrix* gm1, Real f)
00345 {
00346    REPORT
00347    Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
00348    while (i--)
00349    { *s++ = f - *s1++; *s++ = f - *s1++; *s++ = f - *s1++; *s++ = f - *s1++; }
00350    i = storage & 3; while (i--) *s++ = f - *s1++;
00351 }
00352    
00353 void GeneralMatrix::NegAdd(Real f)
00354 {
00355    REPORT
00356    Real* s=store; int i=(storage >> 2);
00357    while (i--)
00358    {
00359       *s = f - *s; s++; *s = f - *s; s++;
00360       *s = f - *s; s++; *s = f - *s; s++;
00361    }
00362    i = storage & 3; while (i--)  { *s = f - *s; s++; }
00363 }
00364    
00365 void GeneralMatrix::Negate(GeneralMatrix* gm1)
00366 {
00367    // change sign of elements
00368    REPORT
00369    Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
00370    while (i--)
00371    { *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); }
00372    i = storage & 3; while(i--) *s++ = -(*s1++);
00373 }
00374    
00375 void GeneralMatrix::Negate()
00376 {
00377    REPORT
00378    Real* s=store; int i=(storage >> 2);
00379    while (i--)
00380    { *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; }
00381    i = storage & 3; while(i--) { *s = -(*s); s++; }
00382 }
00383    
00384 void GeneralMatrix::Multiply(GeneralMatrix* gm1, Real f)
00385 {
00386    REPORT
00387    Real* s1=gm1->store; Real* s=store;  int i=(storage >> 2);
00388    while (i--)
00389    { *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; }
00390    i = storage & 3; while (i--) *s++ = *s1++ * f;
00391 }
00392    
00393 void GeneralMatrix::Multiply(Real f)
00394 {
00395    REPORT
00396    Real* s=store; int i=(storage >> 2);
00397    while (i--) { *s++ *= f; *s++ *= f; *s++ *= f; *s++ *= f; }
00398    i = storage & 3; while (i--) *s++ *= f;
00399 }
00400    
00401 
00402 /************************ MatrixInput routines ****************************/
00403 
00404 // int MatrixInput::n;          // number values still to be read
00405 // Real* MatrixInput::r;        // pointer to next location to be read to
00406 
00407 MatrixInput MatrixInput::operator<<(double f)
00408 {
00409    REPORT
00410    Tracer et("MatrixInput");
00411    if (n<=0) Throw(ProgramException("List of values too long"));
00412    *r = (Real)f; int n1 = n-1; n=0;   // n=0 so we won't trigger exception
00413    return MatrixInput(n1, r+1);
00414 }
00415 
00416 
00417 MatrixInput GeneralMatrix::operator<<(double f)
00418 {
00419    REPORT
00420    Tracer et("MatrixInput");
00421    int n = Storage();
00422    if (n<=0) Throw(ProgramException("Loading data to zero length matrix"));
00423    Real* r; r = Store(); *r = (Real)f; n--;
00424    return MatrixInput(n, r+1);
00425 }
00426 
00427 MatrixInput GetSubMatrix::operator<<(double f)
00428 {
00429    REPORT
00430    Tracer et("MatrixInput (GetSubMatrix)");
00431    SetUpLHS();
00432    if (row_number != 1 || col_skip != 0 || col_number != gm->Ncols())
00433    {
00434       Throw(ProgramException("MatrixInput requires complete rows"));
00435    }
00436    MatrixRow mr(gm, DirectPart, row_skip);  // to pick up location and length
00437    int n = mr.Storage();
00438    if (n<=0)
00439    {
00440       Throw(ProgramException("Loading data to zero length row"));
00441    }
00442    Real* r; r = mr.Data(); *r = (Real)f; n--;
00443    if (+(mr.cw*HaveStore))
00444    {
00445       Throw(ProgramException("Fails with this matrix type"));
00446    }
00447    return MatrixInput(n, r+1);
00448 }
00449 
00450 MatrixInput MatrixInput::operator<<(float f)
00451 {
00452    REPORT
00453    Tracer et("MatrixInput");
00454    if (n<=0) Throw(ProgramException("List of values too long"));
00455    *r = (Real)f; int n1 = n-1; n=0;   // n=0 so we won't trigger exception
00456    return MatrixInput(n1, r+1);
00457 }
00458 
00459 
00460 MatrixInput GeneralMatrix::operator<<(float f)
00461 {
00462    REPORT
00463    Tracer et("MatrixInput");
00464    int n = Storage();
00465    if (n<=0) Throw(ProgramException("Loading data to zero length matrix"));
00466    Real* r; r = Store(); *r = (Real)f; n--;
00467    return MatrixInput(n, r+1);
00468 }
00469 
00470 MatrixInput GetSubMatrix::operator<<(float f)
00471 {
00472    REPORT
00473    Tracer et("MatrixInput (GetSubMatrix)");
00474    SetUpLHS();
00475    if (row_number != 1 || col_skip != 0 || col_number != gm->Ncols())
00476    {
00477       Throw(ProgramException("MatrixInput requires complete rows"));
00478    }
00479    MatrixRow mr(gm, DirectPart, row_skip);  // to pick up location and length
00480    int n = mr.Storage();
00481    if (n<=0)
00482    {
00483       Throw(ProgramException("Loading data to zero length row"));
00484    }
00485    Real* r; r = mr.Data(); *r = (Real)f; n--;
00486    if (+(mr.cw*HaveStore))
00487    {
00488       Throw(ProgramException("Fails with this matrix type"));
00489    }
00490    return MatrixInput(n, r+1);
00491 }
00492 MatrixInput::~MatrixInput()
00493 {
00494    REPORT
00495    Tracer et("MatrixInput");
00496    if (n!=0)
00497      std::cerr << "Error in destructor: A list of values was too short" << std::endl;
00498 }
00499 
00500 MatrixInput BandMatrix::operator<<(double)
00501 {
00502    Tracer et("MatrixInput");
00503    bool dummy = true;
00504    if (dummy)                                   // get rid of warning message
00505       Throw(ProgramException("Cannot use list read with a BandMatrix"));
00506    return MatrixInput(0, 0);
00507 }
00508 
00509 MatrixInput BandMatrix::operator<<(float)
00510 {
00511    Tracer et("MatrixInput");
00512    bool dummy = true;
00513    if (dummy)                                   // get rid of warning message
00514       Throw(ProgramException("Cannot use list read with a BandMatrix"));
00515    return MatrixInput(0, 0);
00516 }
00517 
00518 void BandMatrix::operator<<(const double*)
00519 { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
00520 
00521 void BandMatrix::operator<<(const float*)
00522 { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
00523 
00524 void BandMatrix::operator<<(const int*)
00525 { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
00526 
00527 void SymmetricBandMatrix::operator<<(const double*)
00528 { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
00529 
00530 void SymmetricBandMatrix::operator<<(const float*)
00531 { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
00532 
00533 void SymmetricBandMatrix::operator<<(const int*)
00534 { Throw(ProgramException("Cannot use array read with a BandMatrix")); }
00535 
00536 // ************************* Reverse order of elements ***********************
00537 
00538 void GeneralMatrix::ReverseElements(GeneralMatrix* gm)
00539 {
00540    // reversing into a new matrix
00541    REPORT
00542    int n = Storage(); Real* rx = Store() + n; Real* x = gm->Store();
00543    while (n--) *(--rx) = *(x++);
00544 }
00545 
00546 void GeneralMatrix::ReverseElements()
00547 {
00548    // reversing in place
00549    REPORT
00550    int n = Storage(); Real* x = Store(); Real* rx = x + n;
00551    n /= 2;
00552    while (n--) { Real t = *(--rx); *rx = *x; *(x++) = t; }
00553 }
00554 
00555 
00556 #ifdef use_namespace
00557 }
00558 #endif
00559 


kni
Author(s): Martin Günther
autogenerated on Thu Jun 6 2019 21:42:34