Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

newmat5.cpp

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

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