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


kni
Author(s): Neuronics AG (see AUTHORS.txt); ROS wrapper by Martin Günther
autogenerated on Mon Oct 6 2014 10:45:33