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

Go to the documentation of this file.
00001 //Owner: Fred
00002 //$$ newmat2.cpp      Matrix row and column operations
00003 
00004 // Copyright (C) 1991,2,3,4: R B Davies
00005 
00006 #define WANT_MATH
00007 
00008 #include "include.h"
00009 
00010 #include "newmat.h"
00011 #include "newmatrc.h"
00012 
00013 #ifdef use_namespace
00014 namespace NEWMAT {
00015 #endif
00016 
00017 #ifdef DO_REPORT
00018 #define REPORT { static ExeCounter ExeCount(__LINE__,2); ++ExeCount; }
00019 #else
00020 #define REPORT {}
00021 #endif
00022 
00023     //#define MONITOR(what,storage,store) { cout << what << " " << storage << " at " << (long)store << "\n"; }
00024 
00025 #define MONITOR(what,store,storage) {}
00026 
00027     /************************** Matrix Row/Col functions ************************/
00028 
00029     void MatrixRowCol::Add(const MatrixRowCol& mrc) {
00030         // THIS += mrc
00031         REPORT
00032             int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
00033             if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
00034             if (l<=0) return;
00035             Real* elx=data+(f-skip); Real* el=mrc.data+(f-mrc.skip);
00036             while (l--) *elx++ += *el++;
00037     }
00038 
00039     void MatrixRowCol::AddScaled(const MatrixRowCol& mrc, Real x) {
00040         REPORT
00041             // THIS += (mrc * x)
00042             int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
00043             if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
00044             if (l<=0) return;
00045             Real* elx=data+(f-skip); Real* el=mrc.data+(f-mrc.skip);
00046             while (l--) *elx++ += *el++ * x;
00047     }
00048 
00049     void MatrixRowCol::Sub(const MatrixRowCol& mrc) {
00050         REPORT
00051             // THIS -= mrc
00052             int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
00053             if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
00054             if (l<=0) return;
00055             Real* elx=data+(f-skip); Real* el=mrc.data+(f-mrc.skip);
00056             while (l--) *elx++ -= *el++;
00057     }
00058 
00059     void MatrixRowCol::Inject(const MatrixRowCol& mrc) {
00060         // copy stored elements only
00061         REPORT
00062             int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
00063             if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
00064             if (l<=0) return;
00065             Real* elx=data+(f-skip); Real* ely=mrc.data+(f-mrc.skip);
00066             while (l--) *elx++ = *ely++;
00067     }
00068 
00069     Real DotProd(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2) {
00070         REPORT                                        // not accessed
00071             int f = mrc1.skip; int f2 = mrc2.skip;
00072         int l = f + mrc1.storage; int l2 = f2 + mrc2.storage;
00073         if (f < f2) f = f2; if (l > l2) l = l2; l -= f;
00074         if (l<=0) return 0.0;
00075 
00076         Real* el1=mrc1.data+(f-mrc1.skip); Real* el2=mrc2.data+(f-mrc2.skip);
00077         Real sum = 0.0;
00078         while (l--) sum += *el1++ * *el2++;
00079         return sum;
00080     }
00081 
00082     void MatrixRowCol::Add(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2) {
00083         // THIS = mrc1 + mrc2
00084         int f = skip; int l = skip + storage;
00085         int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
00086         if (f1<f) f1=f; if (l1>l) l1=l;
00087         int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
00088         if (f2<f) f2=f; if (l2>l) l2=l;
00089         Real* el = data + (f-skip);
00090         Real* el1 = mrc1.data+(f1-mrc1.skip); Real* el2 = mrc2.data+(f2-mrc2.skip);
00091         if (f1<f2) {
00092             int i = f1-f; while (i--) *el++ = 0.0;
00093             if (l1<=f2) {                               // disjoint
00094                 REPORT                                    // not accessed
00095                     i = l1-f1;     while (i--) *el++ = *el1++;
00096                 i = f2-l1;     while (i--) *el++ = 0.0;
00097                 i = l2-f2;     while (i--) *el++ = *el2++;
00098                 i = l-l2;      while (i--) *el++ = 0.0;
00099             }
00100             else {
00101                 i = f2-f1;    while (i--) *el++ = *el1++;
00102                 if (l1<=l2) {
00103                     REPORT
00104                         i = l1-f2; while (i--) *el++ = *el1++ + *el2++;
00105                     i = l2-l1; while (i--) *el++ = *el2++;
00106                     i = l-l2;  while (i--) *el++ = 0.0;
00107                 }
00108                 else {
00109                     REPORT
00110                         i = l2-f2; while (i--) *el++ = *el1++ + *el2++;
00111                     i = l1-l2; while (i--) *el++ = *el1++;
00112                     i = l-l1;  while (i--) *el++ = 0.0;
00113                 }
00114             }
00115         }
00116         else {
00117             int i = f2-f; while (i--) *el++ = 0.0;
00118             if (l2<=f1) {                               // disjoint
00119                 REPORT                                    // not accessed
00120                     i = l2-f2;     while (i--) *el++ = *el2++;
00121                 i = f1-l2;     while (i--) *el++ = 0.0;
00122                 i = l1-f1;     while (i--) *el++ = *el1++;
00123                 i = l-l1;      while (i--) *el++ = 0.0;
00124             }
00125             else {
00126                 i = f1-f2;    while (i--) *el++ = *el2++;
00127                 if (l2<=l1) {
00128                     REPORT
00129                         i = l2-f1; while (i--) *el++ = *el1++ + *el2++;
00130                     i = l1-l2; while (i--) *el++ = *el1++;
00131                     i = l-l1;  while (i--) *el++ = 0.0;
00132                 }
00133                 else {
00134                     REPORT
00135                         i = l1-f1; while (i--) *el++ = *el1++ + *el2++;
00136                     i = l2-l1; while (i--) *el++ = *el2++;
00137                     i = l-l2;  while (i--) *el++ = 0.0;
00138                 }
00139             }
00140         }
00141     }
00142 
00143     void MatrixRowCol::Sub(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2) {
00144         // THIS = mrc1 - mrc2
00145         int f = skip; int l = skip + storage;
00146         int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
00147         if (f1<f) f1=f; if (l1>l) l1=l;
00148         int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
00149         if (f2<f) f2=f; if (l2>l) l2=l;
00150         Real* el = data + (f-skip);
00151         Real* el1 = mrc1.data+(f1-mrc1.skip); Real* el2 = mrc2.data+(f2-mrc2.skip);
00152         if (f1<f2) {
00153             int i = f1-f; while (i--) *el++ = 0.0;
00154             if (l1<=f2) {                               // disjoint
00155                 REPORT                                    // not accessed
00156                     i = l1-f1;     while (i--) *el++ = *el1++;
00157                 i = f2-l1;     while (i--) *el++ = 0.0;
00158                 i = l2-f2;     while (i--) *el++ = - *el2++;
00159                 i = l-l2;      while (i--) *el++ = 0.0;
00160             }
00161             else {
00162                 i = f2-f1;    while (i--) *el++ = *el1++;
00163                 if (l1<=l2) {
00164                     REPORT
00165                         i = l1-f2; while (i--) *el++ = *el1++ - *el2++;
00166                     i = l2-l1; while (i--) *el++ = - *el2++;
00167                     i = l-l2;  while (i--) *el++ = 0.0;
00168                 }
00169                 else {
00170                     REPORT
00171                         i = l2-f2; while (i--) *el++ = *el1++ - *el2++;
00172                     i = l1-l2; while (i--) *el++ = *el1++;
00173                     i = l-l1;  while (i--) *el++ = 0.0;
00174                 }
00175             }
00176         }
00177         else {
00178             int i = f2-f; while (i--) *el++ = 0.0;
00179             if (l2<=f1) {                               // disjoint
00180                 REPORT                                    // not accessed
00181                     i = l2-f2;     while (i--) *el++ = - *el2++;
00182                 i = f1-l2;     while (i--) *el++ = 0.0;
00183                 i = l1-f1;     while (i--) *el++ = *el1++;
00184                 i = l-l1;      while (i--) *el++ = 0.0;
00185             }
00186             else {
00187                 i = f1-f2;    while (i--) *el++ = - *el2++;
00188                 if (l2<=l1) {
00189                     REPORT
00190                         i = l2-f1; while (i--) *el++ = *el1++ - *el2++;
00191                     i = l1-l2; while (i--) *el++ = *el1++;
00192                     i = l-l1;  while (i--) *el++ = 0.0;
00193                 }
00194                 else {
00195                     REPORT
00196                         i = l1-f1; while (i--) *el++ = *el1++ - *el2++;
00197                     i = l2-l1; while (i--) *el++ = - *el2++;
00198                     i = l-l2;  while (i--) *el++ = 0.0;
00199                 }
00200             }
00201         }
00202     }
00203 
00204     void MatrixRowCol::Add(const MatrixRowCol& mrc1, Real x) {
00205         // THIS = mrc1 + x
00206         REPORT
00207             int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00208             if (f < skip) { f = skip; if (l < f) l = f; }
00209             if (l > lx) { l = lx; if (f > lx) f = lx; }
00210 
00211             Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00212 
00213             int l1 = f-skip;  while (l1--) *elx++ = x;
00214             l1 = l-f;     while (l1--) *elx++ = *ely++ + x;
00215             lx -= l;      while (lx--) *elx++ = x;
00216     }
00217 
00218     void MatrixRowCol::NegAdd(const MatrixRowCol& mrc1, Real x) {
00219         // THIS = x - mrc1
00220         REPORT
00221             int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00222             if (f < skip) { f = skip; if (l < f) l = f; }
00223             if (l > lx) { l = lx; if (f > lx) f = lx; }
00224 
00225             Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00226 
00227             int l1 = f-skip;  while (l1--) *elx++ = x;
00228             l1 = l-f;     while (l1--) *elx++ = x - *ely++;
00229             lx -= l;      while (lx--) *elx++ = x;
00230     }
00231 
00232     void MatrixRowCol::RevSub(const MatrixRowCol& mrc1) {
00233         // THIS = mrc - THIS
00234         REPORT
00235             int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00236             if (f < skip) { f = skip; if (l < f) l = f; }
00237             if (l > lx) { l = lx; if (f > lx) f = lx; }
00238 
00239             Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00240 
00241             int l1 = f-skip;  while (l1--) { *elx = - *elx; elx++; }
00242             l1 = l-f;     while (l1--) { *elx = *ely++ - *elx; elx++; }
00243             lx -= l;      while (lx--) { *elx = - *elx; elx++; }
00244     }
00245 
00246     void MatrixRowCol::ConCat(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2) {
00247         // THIS = mrc1 | mrc2
00248         REPORT
00249             int f1 = mrc1.skip; int l1 = f1 + mrc1.storage; int lx = skip + storage;
00250             if (f1 < skip) { f1 = skip; if (l1 < f1) l1 = f1; }
00251             if (l1 > lx) { l1 = lx; if (f1 > lx) f1 = lx; }
00252 
00253             Real* elx = data;
00254 
00255             int i = f1-skip;  while (i--) *elx++ =0.0;
00256             i = l1-f1;
00257             if (i)                                        // in case f1 would take ely out of range
00258             { Real* ely = mrc1.data+(f1-mrc1.skip);  while (i--) *elx++ = *ely++; }
00259 
00260             int f2 = mrc2.skip; int l2 = f2 + mrc2.storage; i = mrc1.length;
00261             int skipx = l1 - i; lx -= i;                  // addresses rel to second seg, maybe -ve
00262             if (f2 < skipx) { f2 = skipx; if (l2 < f2) l2 = f2; }
00263             if (l2 > lx) { l2 = lx; if (f2 > lx) f2 = lx; }
00264 
00265             i = f2-skipx; while (i--) *elx++ = 0.0;
00266             i = l2-f2;
00267             if (i)                                        // in case f2 would take ely out of range
00268             { Real* ely = mrc2.data+(f2-mrc2.skip); while (i--) *elx++ = *ely++; }
00269             lx -= l2;     while (lx--) *elx++ = 0.0;
00270     }
00271 
00272     void MatrixRowCol::Multiply(const MatrixRowCol& mrc1) {
00273         // element by element multiply into
00274         REPORT
00275             int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00276             if (f < skip) { f = skip; if (l < f) l = f; }
00277             if (l > lx) { l = lx; if (f > lx) f = lx; }
00278 
00279             Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00280 
00281             int l1 = f-skip;  while (l1--) *elx++ = 0;
00282             l1 = l-f;     while (l1--) *elx++ *= *ely++;
00283             lx -= l;      while (lx--) *elx++ = 0;
00284     }
00285 
00286     void MatrixRowCol::Multiply(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2) {
00287         // element by element multiply
00288         int f = skip; int l = skip + storage;
00289         int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
00290         if (f1<f) f1=f; if (l1>l) l1=l;
00291         int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
00292         if (f2<f) f2=f; if (l2>l) l2=l;
00293         Real* el = data + (f-skip); int i;
00294         if (f1<f2) f1 = f2; if (l1>l2) l1 = l2;
00295         if (l1<=f1) {                                 // disjoint
00296             REPORT i = l-f; while (i--) *el++ = 0.0;
00297         }
00298         else {
00299             REPORT
00300                 Real* el1 = mrc1.data+(f1-mrc1.skip);
00301             Real* el2 = mrc2.data+(f1-mrc2.skip);
00302             i = f1-f ;    while (i--) *el++ = 0.0;
00303             i = l1-f1;    while (i--) *el++ = *el1++ * *el2++;
00304             i = l-l1;     while (i--) *el++ = 0.0;
00305         }
00306     }
00307 
00308     void MatrixRowCol::Copy(const MatrixRowCol& mrc1) {
00309         // THIS = mrc1
00310         REPORT
00311             int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00312             if (f < skip) { f = skip; if (l < f) l = f; }
00313             if (l > lx) { l = lx; if (f > lx) f = lx; }
00314 
00315             Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00316 
00317             int l1 = f-skip;  while (l1--) *elx++ = 0.0;
00318             l1 = l-f;     while (l1--) *elx++ = *ely++;
00319             lx -= l;      while (lx--) *elx++ = 0.0;
00320     }
00321 
00322     void MatrixRowCol::CopyCheck(const MatrixRowCol& mrc1) {
00323         // Throw an exception if this would lead to a loss of data
00324         REPORT
00325             int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00326             if (f < skip || l > lx) Throw(ProgramException("Illegal Conversion"));
00327 
00328             Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00329 
00330             int l1 = f-skip;  while (l1--) *elx++ = 0.0;
00331             l1 = l-f;     while (l1--) *elx++ = *ely++;
00332             lx -= l;      while (lx--) *elx++ = 0.0;
00333     }
00334 
00335     void MatrixRowCol::Check(const MatrixRowCol& mrc1) {
00336         // Throw an exception if +=, -=, copy etc would lead to a loss of data
00337         REPORT
00338             int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00339             if (f < skip || l > lx) Throw(ProgramException("Illegal Conversion"));
00340     }
00341 
00342     void MatrixRowCol::Check() {
00343         // Throw an exception if +=, -= of constant would lead to a loss of data
00344         // that is: check full row is present
00345         // may not be appropriate for symmetric matrices
00346         REPORT
00347             if (skip!=0 || storage!=length)
00348                 Throw(ProgramException("Illegal Conversion"));
00349     }
00350 
00351     void MatrixRowCol::Negate(const MatrixRowCol& mrc1) {
00352         // THIS = -mrc1
00353         REPORT
00354             int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00355             if (f < skip) { f = skip; if (l < f) l = f; }
00356             if (l > lx) { l = lx; if (f > lx) f = lx; }
00357 
00358             Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00359 
00360             int l1 = f-skip;  while (l1--) *elx++ = 0.0;
00361             l1 = l-f;     while (l1--) *elx++ = - *ely++;
00362             lx -= l;      while (lx--) *elx++ = 0.0;
00363     }
00364 
00365     void MatrixRowCol::Multiply(const MatrixRowCol& mrc1, Real s) {
00366         // THIS = mrc1 * s
00367         REPORT
00368             int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00369             if (f < skip) { f = skip; if (l < f) l = f; }
00370             if (l > lx) { l = lx; if (f > lx) f = lx; }
00371 
00372             Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00373 
00374             int l1 = f-skip;  while (l1--) *elx++ = 0.0;
00375             l1 = l-f;     while (l1--) *elx++ = *ely++ * s;
00376             lx -= l;      while (lx--) *elx++ = 0.0;
00377     }
00378 
00379     void DiagonalMatrix::Solver(MatrixColX& mrc, const MatrixColX& mrc1) {
00380         // mrc = mrc / mrc1   (elementwise)
00381         REPORT
00382             int f = mrc1.skip; int f0 = mrc.skip;
00383         int l = f + mrc1.storage; int lx = f0 + mrc.storage;
00384         if (f < f0) { f = f0; if (l < f) l = f; }
00385         if (l > lx) { l = lx; if (f > lx) f = lx; }
00386 
00387         Real* elx = mrc.data; Real* eld = store+f;
00388 
00389         int l1 = f-f0;    while (l1--) *elx++ = 0.0;
00390         l1 = l-f;     while (l1--) *elx++ /= *eld++;
00391         lx -= l;      while (lx--) *elx++ = 0.0;
00392         // Solver makes sure input and output point to same memory
00393     }
00394 
00395     void MatrixRowCol::Copy(const Real*& r) {
00396         // THIS = *r
00397         REPORT
00398             Real* elx = data; const Real* ely = r+skip; r += length;
00399             int l = storage; while (l--) *elx++ = *ely++;
00400     }
00401 
00402     void MatrixRowCol::Copy(Real r) {
00403         // THIS = r
00404         REPORT
00405             Real* elx = data; int l = storage; while (l--) *elx++ = r;
00406     }
00407 
00408     void MatrixRowCol::Multiply(Real r) {
00409         // THIS *= r
00410         REPORT
00411             Real* elx = data; int l = storage; while (l--) *elx++ *= r;
00412     }
00413 
00414     void MatrixRowCol::Add(Real r) {
00415         // THIS += r
00416         REPORT
00417             Real* elx = data; int l = storage; while (l--) *elx++ += r;
00418     }
00419 
00420     Real MatrixRowCol::SumAbsoluteValue() {
00421         REPORT
00422             Real sum = 0.0; Real* elx = data; int l = storage;
00423             while (l--) sum += fabs(*elx++);
00424             return sum;
00425     }
00426 
00427     Real MatrixRowCol::Sum() {
00428         REPORT
00429             Real sum = 0.0; Real* elx = data; int l = storage;
00430             while (l--) sum += *elx++;
00431             return sum;
00432     }
00433 
00434     void MatrixRowCol::SubRowCol(MatrixRowCol& mrc, int skip1, int l1) const {
00435         mrc.length = l1;  int d = skip - skip1;
00436         if (d<0) { mrc.skip = 0; mrc.data = data - d; }
00437         else  { mrc.skip = d; mrc.data = data; }
00438         d = skip + storage - skip1;
00439         d = ((l1 < d) ? l1 : d) - mrc.skip;  mrc.storage = (d < 0) ? 0 : d;
00440         mrc.cw = 0;
00441     }
00442 
00443 #ifdef use_namespace
00444 }
00445 #endif
Generated on Fri Aug 23 01:37:05 2002 for VirtualU by
1.2.17