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


rl_agent
Author(s): Todd Hester
autogenerated on Thu Jun 6 2019 22:00:13