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


kni
Author(s): Martin Günther
autogenerated on Mon Aug 14 2017 02:44:12