NEWMAT5.CPP Source File Main Page Namespace List Class Hierarchy Alphabetical List Compound List File List Namespace Members Compound Members File Members
NEWMAT5.CPP

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