Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

newmat2.cpp

Go to the documentation of this file.
00001 //$$ newmat2.cpp      Matrix row and column operations
00002 
00003 // Copyright (C) 1991,2,3,4: R B Davies
00004 
00005 #ifndef WANT_MATH
00006 #define WANT_MATH
00007 #endif
00008 
00009 #include "include.h"
00010 
00011 #include "newmat.h"
00012 #include "newmatrc.h"
00013 
00014 #ifdef use_namespace
00015 namespace NEWMAT {
00016 #endif
00017 
00018 
00019 #ifdef DO_REPORT
00020 #define REPORT { static ExeCounter ExeCount(__LINE__,2); ++ExeCount; }
00021 #else
00022 #define REPORT {}
00023 #endif
00024 
00025 //#define MONITOR(what,storage,store) { cout << what << " " << storage << " at " << (long)store << "\n"; }
00026 
00027 #define MONITOR(what,store,storage) {}
00028 
00029 /************************** Matrix Row/Col functions ************************/
00030 
00031 void MatrixRowCol::Add(const MatrixRowCol& mrc)
00032 {
00033    // THIS += mrc
00034    REPORT
00035    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
00036    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
00037    if (l<=0) return;
00038    Real* elx=data+(f-skip); Real* el=mrc.data+(f-mrc.skip);
00039    while (l--) *elx++ += *el++;
00040 }
00041 
00042 void MatrixRowCol::AddScaled(const MatrixRowCol& mrc, Real x)
00043 {
00044    REPORT
00045    // THIS += (mrc * x)
00046    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
00047    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
00048    if (l<=0) return;
00049    Real* elx=data+(f-skip); Real* el=mrc.data+(f-mrc.skip);
00050    while (l--) *elx++ += *el++ * x;
00051 }
00052 
00053 void MatrixRowCol::Sub(const MatrixRowCol& mrc)
00054 {
00055    REPORT
00056    // THIS -= mrc
00057    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
00058    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
00059    if (l<=0) return;
00060    Real* elx=data+(f-skip); Real* el=mrc.data+(f-mrc.skip);
00061    while (l--) *elx++ -= *el++;
00062 }
00063 
00064 void MatrixRowCol::Inject(const MatrixRowCol& mrc)
00065 // copy stored elements only
00066 {
00067    REPORT
00068    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
00069    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
00070    if (l<=0) return;
00071    Real* elx=data+(f-skip); Real* ely=mrc.data+(f-mrc.skip);
00072    while (l--) *elx++ = *ely++;
00073 }
00074 
00075 Real DotProd(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
00076 {
00077    REPORT                                         // not accessed
00078    int f = mrc1.skip; int f2 = mrc2.skip;
00079    int l = f + mrc1.storage; int l2 = f2 + mrc2.storage;
00080    if (f < f2) f = f2; if (l > l2) l = l2; l -= f;
00081    if (l<=0) return 0.0;
00082 
00083    Real* el1=mrc1.data+(f-mrc1.skip); Real* el2=mrc2.data+(f-mrc2.skip);
00084    Real sum = 0.0;
00085    while (l--) sum += *el1++ * *el2++;
00086    return sum;
00087 }
00088 
00089 void MatrixRowCol::Add(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
00090 {
00091    // THIS = mrc1 + mrc2
00092    int f = skip; int l = skip + storage;
00093    int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
00094    if (f1<f) f1=f; if (l1>l) l1=l;
00095    int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
00096    if (f2<f) f2=f; if (l2>l) l2=l;
00097    Real* el = data + (f-skip);
00098    Real* el1 = mrc1.data+(f1-mrc1.skip); Real* el2 = mrc2.data+(f2-mrc2.skip);
00099    if (f1<f2)
00100    {
00101       int i = f1-f; while (i--) *el++ = 0.0;
00102       if (l1<=f2)                              // disjoint
00103       {
00104          REPORT                                // not accessed
00105          i = l1-f1;     while (i--) *el++ = *el1++;
00106          i = f2-l1;     while (i--) *el++ = 0.0;
00107          i = l2-f2;     while (i--) *el++ = *el2++;
00108          i = l-l2;      while (i--) *el++ = 0.0;
00109       }
00110       else
00111       {
00112          i = f2-f1;    while (i--) *el++ = *el1++;
00113          if (l1<=l2)
00114          {
00115             REPORT
00116             i = l1-f2; while (i--) *el++ = *el1++ + *el2++;
00117             i = l2-l1; while (i--) *el++ = *el2++;
00118             i = l-l2;  while (i--) *el++ = 0.0;
00119          }
00120          else
00121          {
00122             REPORT
00123             i = l2-f2; while (i--) *el++ = *el1++ + *el2++;
00124             i = l1-l2; while (i--) *el++ = *el1++;
00125             i = l-l1;  while (i--) *el++ = 0.0;
00126          }
00127       }
00128    }
00129    else
00130    {
00131       int i = f2-f; while (i--) *el++ = 0.0;
00132       if (l2<=f1)                              // disjoint
00133       {
00134          REPORT                                // not accessed
00135          i = l2-f2;     while (i--) *el++ = *el2++;
00136          i = f1-l2;     while (i--) *el++ = 0.0;
00137          i = l1-f1;     while (i--) *el++ = *el1++;
00138          i = l-l1;      while (i--) *el++ = 0.0;
00139       }
00140       else
00141       {
00142          i = f1-f2;    while (i--) *el++ = *el2++;
00143          if (l2<=l1)
00144          {
00145             REPORT
00146             i = l2-f1; while (i--) *el++ = *el1++ + *el2++;
00147             i = l1-l2; while (i--) *el++ = *el1++;
00148             i = l-l1;  while (i--) *el++ = 0.0;
00149          }
00150          else
00151          {
00152             REPORT
00153             i = l1-f1; while (i--) *el++ = *el1++ + *el2++;
00154             i = l2-l1; while (i--) *el++ = *el2++;
00155             i = l-l2;  while (i--) *el++ = 0.0;
00156          }
00157       }
00158    }
00159 }
00160 
00161 void MatrixRowCol::Sub(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
00162 {
00163    // THIS = mrc1 - mrc2
00164    int f = skip; int l = skip + storage;
00165    int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
00166    if (f1<f) f1=f; if (l1>l) l1=l;
00167    int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
00168    if (f2<f) f2=f; if (l2>l) l2=l;
00169    Real* el = data + (f-skip);
00170    Real* el1 = mrc1.data+(f1-mrc1.skip); Real* el2 = mrc2.data+(f2-mrc2.skip);
00171    if (f1<f2)
00172    {
00173       int i = f1-f; while (i--) *el++ = 0.0;
00174       if (l1<=f2)                              // disjoint
00175       {
00176          REPORT                                // not accessed
00177          i = l1-f1;     while (i--) *el++ = *el1++;
00178          i = f2-l1;     while (i--) *el++ = 0.0;
00179          i = l2-f2;     while (i--) *el++ = - *el2++;
00180          i = l-l2;      while (i--) *el++ = 0.0;
00181       }
00182       else
00183       {
00184          i = f2-f1;    while (i--) *el++ = *el1++;
00185          if (l1<=l2)
00186          {
00187             REPORT
00188             i = l1-f2; while (i--) *el++ = *el1++ - *el2++;
00189             i = l2-l1; while (i--) *el++ = - *el2++;
00190             i = l-l2;  while (i--) *el++ = 0.0;
00191          }
00192          else
00193          {
00194             REPORT
00195             i = l2-f2; while (i--) *el++ = *el1++ - *el2++;
00196             i = l1-l2; while (i--) *el++ = *el1++;
00197             i = l-l1;  while (i--) *el++ = 0.0;
00198          }
00199       }
00200    }
00201    else
00202    {
00203       int i = f2-f; while (i--) *el++ = 0.0;
00204       if (l2<=f1)                              // disjoint
00205       {
00206          REPORT                                // not accessed
00207          i = l2-f2;     while (i--) *el++ = - *el2++;
00208          i = f1-l2;     while (i--) *el++ = 0.0;
00209          i = l1-f1;     while (i--) *el++ = *el1++;
00210          i = l-l1;      while (i--) *el++ = 0.0;
00211       }
00212       else
00213       {
00214          i = f1-f2;    while (i--) *el++ = - *el2++;
00215          if (l2<=l1)
00216          {
00217             REPORT
00218             i = l2-f1; while (i--) *el++ = *el1++ - *el2++;
00219             i = l1-l2; while (i--) *el++ = *el1++;
00220             i = l-l1;  while (i--) *el++ = 0.0;
00221          }
00222          else
00223          {
00224             REPORT
00225             i = l1-f1; while (i--) *el++ = *el1++ - *el2++;
00226             i = l2-l1; while (i--) *el++ = - *el2++;
00227             i = l-l2;  while (i--) *el++ = 0.0;
00228          }
00229       }
00230    }
00231 }
00232 
00233 
00234 void MatrixRowCol::Add(const MatrixRowCol& mrc1, Real x)
00235 {
00236    // THIS = mrc1 + x
00237    REPORT
00238    if (!storage) return;
00239    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00240    if (f < skip) { f = skip; if (l < f) l = f; }
00241    if (l > lx) { l = lx; if (f > lx) f = lx; }
00242 
00243    Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00244 
00245    int l1 = f-skip;  while (l1--) *elx++ = x;
00246        l1 = l-f;     while (l1--) *elx++ = *ely++ + x;
00247        lx -= l;      while (lx--) *elx++ = x;
00248 }
00249 
00250 void MatrixRowCol::NegAdd(const MatrixRowCol& mrc1, Real x)
00251 {
00252    // THIS = x - mrc1
00253    REPORT
00254    if (!storage) return;
00255    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00256    if (f < skip) { f = skip; if (l < f) l = f; }
00257    if (l > lx) { l = lx; if (f > lx) f = lx; }
00258 
00259    Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00260 
00261    int l1 = f-skip;  while (l1--) *elx++ = x;
00262        l1 = l-f;     while (l1--) *elx++ = x - *ely++;
00263        lx -= l;      while (lx--) *elx++ = x;
00264 }
00265 
00266 void MatrixRowCol::RevSub(const MatrixRowCol& mrc1)
00267 {
00268    // THIS = mrc - THIS
00269    REPORT
00270    if (!storage) return;
00271    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00272    if (f < skip) { f = skip; if (l < f) l = f; }
00273    if (l > lx) { l = lx; if (f > lx) f = lx; }
00274 
00275    Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00276 
00277    int l1 = f-skip;  while (l1--) { *elx = - *elx; elx++; }
00278        l1 = l-f;     while (l1--) { *elx = *ely++ - *elx; elx++; }
00279        lx -= l;      while (lx--) { *elx = - *elx; elx++; }
00280 }
00281 
00282 void MatrixRowCol::ConCat(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
00283 {
00284    // THIS = mrc1 | mrc2
00285    REPORT
00286    int f1 = mrc1.skip; int l1 = f1 + mrc1.storage; int lx = skip + storage;
00287    if (f1 < skip) { f1 = skip; if (l1 < f1) l1 = f1; }
00288    if (l1 > lx) { l1 = lx; if (f1 > lx) f1 = lx; }
00289 
00290    Real* elx = data;
00291 
00292    int i = f1-skip;  while (i--) *elx++ =0.0;
00293    i = l1-f1;
00294    if (i)                       // in case f1 would take ely out of range
00295       { Real* ely = mrc1.data+(f1-mrc1.skip);  while (i--) *elx++ = *ely++; }
00296 
00297    int f2 = mrc2.skip; int l2 = f2 + mrc2.storage; i = mrc1.length;
00298    int skipx = l1 - i; lx -= i; // addresses rel to second seg, maybe -ve
00299    if (f2 < skipx) { f2 = skipx; if (l2 < f2) l2 = f2; }
00300    if (l2 > lx) { l2 = lx; if (f2 > lx) f2 = lx; }
00301 
00302    i = f2-skipx; while (i--) *elx++ = 0.0;
00303    i = l2-f2;
00304    if (i)                       // in case f2 would take ely out of range
00305       { Real* ely = mrc2.data+(f2-mrc2.skip); while (i--) *elx++ = *ely++; }
00306    lx -= l2;     while (lx--) *elx++ = 0.0;
00307 }
00308 
00309 void MatrixRowCol::Multiply(const MatrixRowCol& mrc1)
00310 // element by element multiply into
00311 {
00312    REPORT
00313    if (!storage) return;
00314    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00315    if (f < skip) { f = skip; if (l < f) l = f; }
00316    if (l > lx) { l = lx; if (f > lx) f = lx; }
00317 
00318    Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00319 
00320    int l1 = f-skip;  while (l1--) *elx++ = 0;
00321        l1 = l-f;     while (l1--) *elx++ *= *ely++;
00322        lx -= l;      while (lx--) *elx++ = 0;
00323 }
00324 
00325 void MatrixRowCol::Multiply(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
00326 // element by element multiply
00327 {
00328    int f = skip; int l = skip + storage;
00329    int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
00330    if (f1<f) f1=f; if (l1>l) l1=l;
00331    int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
00332    if (f2<f) f2=f; if (l2>l) l2=l;
00333    Real* el = data + (f-skip); int i;
00334    if (f1<f2) f1 = f2; if (l1>l2) l1 = l2;
00335    if (l1<=f1) { REPORT i = l-f; while (i--) *el++ = 0.0; }  // disjoint
00336    else
00337    {
00338       REPORT
00339       Real* el1 = mrc1.data+(f1-mrc1.skip);
00340       Real* el2 = mrc2.data+(f1-mrc2.skip);
00341       i = f1-f ;    while (i--) *el++ = 0.0;
00342       i = l1-f1;    while (i--) *el++ = *el1++ * *el2++;
00343       i = l-l1;     while (i--) *el++ = 0.0;
00344    }
00345 }
00346 
00347 void MatrixRowCol::KP(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
00348 // row for Kronecker product
00349 {
00350    int f = skip; int s = storage; Real* el = data; int i;
00351 
00352    i = mrc1.skip * mrc2.length;
00353    if (i > f)
00354    {
00355       i -= f; f = 0; if (i > s) { i = s; s = 0; }  else s -= i;
00356       while (i--) *el++ = 0.0;
00357       if (s == 0) return;
00358    }
00359    else f -= i;
00360 
00361    i = mrc1.storage; Real* el1 = mrc1.data;
00362    int mrc2_skip = mrc2.skip; int mrc2_storage = mrc2.storage;
00363    int mrc2_length = mrc2.length;
00364    int mrc2_remain = mrc2_length - mrc2_skip - mrc2_storage;
00365    while (i--)
00366    {
00367       int j; Real* el2 = mrc2.data; Real vel1 = *el1;
00368       if (f == 0 && mrc2_length <= s)
00369       {
00370          j = mrc2_skip; s -= j;    while (j--) *el++ = 0.0;
00371          j = mrc2_storage; s -= j; while (j--) *el++ = vel1 * *el2++;
00372          j = mrc2_remain; s -= j;  while (j--) *el++ = 0.0;
00373       }
00374       else if (f >= mrc2_length) f -= mrc2_length;
00375       else
00376       {
00377          j = mrc2_skip;
00378          if (j > f)
00379          {
00380             j -= f; f = 0; if (j > s) { j = s; s = 0; } else s -= j;
00381             while (j--) *el++ = 0.0;
00382          }
00383          else f -= j;
00384 
00385          j = mrc2_storage;
00386          if (j > f)
00387          {
00388             j -= f; el2 += f; f = 0; if (j > s) { j = s; s = 0; } else s -= j;
00389             while (j--) *el++ = vel1 * *el2++;
00390          }
00391          else f -= j;
00392 
00393          j = mrc2_remain;
00394          if (j > f)
00395          {
00396             j -= f; f = 0; if (j > s) { j = s; s = 0; } else s -= j;
00397             while (j--) *el++ = 0.0;
00398          }
00399          else f -= j;
00400       }
00401       if (s == 0) return;
00402       ++el1;
00403    }
00404 
00405    i = (mrc1.length - mrc1.skip - mrc1.storage) * mrc2.length;
00406    if (i > f)
00407    {
00408       i -= f; if (i > s) i = s;
00409       while (i--) *el++ = 0.0;
00410    }
00411 }
00412 
00413 
00414 void MatrixRowCol::Copy(const MatrixRowCol& mrc1)
00415 {
00416    // THIS = mrc1
00417    REPORT
00418    if (!storage) return;
00419    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00420    if (f < skip) { f = skip; if (l < f) l = f; }
00421    if (l > lx) { l = lx; if (f > lx) f = lx; }
00422 
00423    Real* elx = data; Real* ely = 0;
00424 
00425    if (l-f) ely = mrc1.data+(f-mrc1.skip);
00426 
00427    int l1 = f-skip;  while (l1--) *elx++ = 0.0;
00428        l1 = l-f;     while (l1--) *elx++ = *ely++;
00429        lx -= l;      while (lx--) *elx++ = 0.0;
00430 }
00431 
00432 void MatrixRowCol::CopyCheck(const MatrixRowCol& mrc1)
00433 // Throw an exception if this would lead to a loss of data
00434 {
00435    REPORT
00436    if (!storage) return;
00437    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00438    if (f < skip || l > lx) Throw(ProgramException("Illegal Conversion"));
00439 
00440    Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00441 
00442    int l1 = f-skip;  while (l1--) *elx++ = 0.0;
00443        l1 = l-f;     while (l1--) *elx++ = *ely++;
00444        lx -= l;      while (lx--) *elx++ = 0.0;
00445 }
00446 
00447 void MatrixRowCol::Check(const MatrixRowCol& mrc1)
00448 // Throw an exception if +=, -=, copy etc would lead to a loss of data
00449 {
00450    REPORT
00451    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00452    if (f < skip || l > lx) Throw(ProgramException("Illegal Conversion"));
00453 }
00454 
00455 void MatrixRowCol::Check()
00456 // Throw an exception if +=, -= of constant would lead to a loss of data
00457 // that is: check full row is present
00458 // may not be appropriate for symmetric matrices
00459 {
00460    REPORT
00461    if (skip!=0 || storage!=length)
00462       Throw(ProgramException("Illegal Conversion"));
00463 }
00464 
00465 void MatrixRowCol::Negate(const MatrixRowCol& mrc1)
00466 {
00467    // THIS = -mrc1
00468    REPORT
00469    if (!storage) return;
00470    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00471    if (f < skip) { f = skip; if (l < f) l = f; }
00472    if (l > lx) { l = lx; if (f > lx) f = lx; }
00473 
00474    Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00475 
00476    int l1 = f-skip;  while (l1--) *elx++ = 0.0;
00477        l1 = l-f;     while (l1--) *elx++ = - *ely++;
00478        lx -= l;      while (lx--) *elx++ = 0.0;
00479 }
00480 
00481 void MatrixRowCol::Multiply(const MatrixRowCol& mrc1, Real s)
00482 {
00483    // THIS = mrc1 * s
00484    REPORT
00485    if (!storage) return;
00486    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
00487    if (f < skip) { f = skip; if (l < f) l = f; }
00488    if (l > lx) { l = lx; if (f > lx) f = lx; }
00489 
00490    Real* elx = data; Real* ely = mrc1.data+(f-mrc1.skip);
00491 
00492    int l1 = f-skip;  while (l1--) *elx++ = 0.0;
00493        l1 = l-f;     while (l1--) *elx++ = *ely++ * s;
00494        lx -= l;      while (lx--) *elx++ = 0.0;
00495 }
00496 
00497 void DiagonalMatrix::Solver(MatrixColX& mrc, const MatrixColX& mrc1)
00498 {
00499    // mrc = mrc / mrc1   (elementwise)
00500    REPORT
00501    int f = mrc1.skip; int f0 = mrc.skip;
00502    int l = f + mrc1.storage; int lx = f0 + mrc.storage;
00503    if (f < f0) { f = f0; if (l < f) l = f; }
00504    if (l > lx) { l = lx; if (f > lx) f = lx; }
00505 
00506    Real* elx = mrc.data; Real* eld = store+f;
00507 
00508    int l1 = f-f0;    while (l1--) *elx++ = 0.0;
00509        l1 = l-f;     while (l1--) *elx++ /= *eld++;
00510        lx -= l;      while (lx--) *elx++ = 0.0;
00511    // Solver makes sure input and output point to same memory
00512 }
00513 
00514 void IdentityMatrix::Solver(MatrixColX& mrc, const MatrixColX& mrc1)
00515 {
00516    // mrc = mrc / mrc1   (elementwise)
00517    REPORT
00518    int f = mrc1.skip; int f0 = mrc.skip;
00519    int l = f + mrc1.storage; int lx = f0 + mrc.storage;
00520    if (f < f0) { f = f0; if (l < f) l = f; }
00521    if (l > lx) { l = lx; if (f > lx) f = lx; }
00522 
00523    Real* elx = mrc.data; Real eldv = *store;
00524 
00525    int l1 = f-f0;    while (l1--) *elx++ = 0.0;
00526        l1 = l-f;     while (l1--) *elx++ /= eldv;
00527        lx -= l;      while (lx--) *elx++ = 0.0;
00528    // Solver makes sure input and output point to same memory
00529 }
00530 
00531 void MatrixRowCol::Copy(const Real*& r)
00532 {
00533    // THIS = *r
00534    REPORT
00535    Real* elx = data; const Real* ely = r+skip; r += length;
00536    int l = storage; while (l--) *elx++ = *ely++;
00537 }
00538 
00539 void MatrixRowCol::Copy(const int*& r)
00540 {
00541    // THIS = *r
00542    REPORT
00543    Real* elx = data; const int* ely = r+skip; r += length;
00544    int l = storage; while (l--) *elx++ = *ely++;
00545 }
00546 
00547 void MatrixRowCol::Copy(Real r)
00548 {
00549    // THIS = r
00550    REPORT  Real* elx = data; int l = storage; while (l--) *elx++ = r;
00551 }
00552 
00553 void MatrixRowCol::Zero()
00554 {
00555    // THIS = 0
00556    REPORT  Real* elx = data; int l = storage; while (l--) *elx++ = 0;
00557 }
00558 
00559 void MatrixRowCol::Multiply(Real r)
00560 {
00561    // THIS *= r
00562    REPORT  Real* elx = data; int l = storage; while (l--) *elx++ *= r;
00563 }
00564 
00565 void MatrixRowCol::Add(Real r)
00566 {
00567    // THIS += r
00568    REPORT
00569    Real* elx = data; int l = storage; while (l--) *elx++ += r;
00570 }
00571 
00572 Real MatrixRowCol::SumAbsoluteValue()
00573 {
00574    REPORT
00575    Real sum = 0.0; Real* elx = data; int l = storage;
00576    while (l--) sum += fabs(*elx++);
00577    return sum;
00578 }
00579 
00580 // max absolute value of r and elements of row/col
00581 // we use <= or >= in all of these so we are sure of getting
00582 // r reset at least once.
00583 Real MatrixRowCol::MaximumAbsoluteValue1(Real r, int& i)
00584 {
00585    REPORT
00586    Real* elx = data; int l = storage; int li = -1;
00587    while (l--) { Real f = fabs(*elx++); if (r <= f) { r = f; li = l; } }
00588    i = (li >= 0) ? storage - li + skip : 0;
00589    return r;
00590 }
00591 
00592 // min absolute value of r and elements of row/col
00593 Real MatrixRowCol::MinimumAbsoluteValue1(Real r, int& i)
00594 {
00595    REPORT
00596    Real* elx = data; int l = storage; int li = -1;
00597    while (l--) { Real f = fabs(*elx++); if (r >= f) { r = f; li = l; } }
00598    i = (li >= 0) ? storage - li + skip : 0;
00599    return r;
00600 }
00601 
00602 // max value of r and elements of row/col
00603 Real MatrixRowCol::Maximum1(Real r, int& i)
00604 {
00605    REPORT
00606    Real* elx = data; int l = storage; int li = -1;
00607    while (l--) { Real f = *elx++; if (r <= f) { r = f; li = l; } }
00608    i = (li >= 0) ? storage - li + skip : 0;
00609    return r;
00610 }
00611 
00612 // min value of r and elements of row/col
00613 Real MatrixRowCol::Minimum1(Real r, int& i)
00614 {
00615    REPORT
00616    Real* elx = data; int l = storage; int li = -1;
00617    while (l--) { Real f = *elx++; if (r >= f) { r = f; li = l; } }
00618    i = (li >= 0) ? storage - li + skip : 0;
00619    return r;
00620 }
00621 
00622 Real MatrixRowCol::Sum()
00623 {
00624    REPORT
00625    Real sum = 0.0; Real* elx = data; int l = storage;
00626    while (l--) sum += *elx++;
00627    return sum;
00628 }
00629 
00630 void MatrixRowCol::SubRowCol(MatrixRowCol& mrc, int skip1, int l1) const
00631 {
00632    mrc.length = l1;  int d = skip - skip1;
00633    if (d<0) { mrc.skip = 0; mrc.data = data - d; }
00634    else  { mrc.skip = d; mrc.data = data; }
00635    d = skip + storage - skip1;
00636    d = ((l1 < d) ? l1 : d) - mrc.skip;  mrc.storage = (d < 0) ? 0 : d;
00637    mrc.cw = 0;
00638 }
00639 
00640 #ifdef use_namespace
00641 }
00642 #endif
00643 

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