home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume26 / newmat04 / part03 / newmat2.cxx < prev    next >
Encoding:
C/C++ Source or Header  |  1991-11-30  |  10.2 KB  |  349 lines
  1. //$$ newmat2.cxx      Matrix row and column operations
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5. #define WANT_MATH
  6.  
  7. #include "include.hxx"
  8.  
  9. #include "newmat.hxx"
  10. #include "newmatrc.hxx"
  11.  
  12. //#define REPORT { static ExeCounter ExeCount(__LINE__,2); ExeCount++; }
  13.  
  14. #define REPORT {}
  15.  
  16. //#define MONITOR(what,storage,store) \
  17. //   { cout << what << " " << storage << " at " << (long)store << "\n"; }
  18.  
  19. #define MONITOR(what,store,storage) {}
  20.  
  21. /************************** Matrix Row/Col functions ************************/
  22.  
  23. void MatrixRowCol::Add(const MatrixRowCol& mrc)
  24. {
  25.    REPORT
  26.    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  27.    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  28.    if (l<=0) return;
  29.    real* elx=store+f; real* el=mrc.store+f;
  30.    while (l--) *elx++ += *el++;
  31. }
  32.  
  33. void MatrixRowCol::AddScaled(const MatrixRowCol& mrc, real x)
  34. {
  35.    REPORT
  36.    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  37.    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  38.    if (l<=0) return;
  39.    real* elx=store+f; real* el=mrc.store+f;
  40.    while (l--) *elx++ += *el++ * x;
  41. }
  42.  
  43. void MatrixRowCol::Sub(const MatrixRowCol& mrc)
  44. {
  45.    REPORT
  46.    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  47.    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  48.    if (l<=0) return;
  49.    real* elx=store+f; real* el=mrc.store+f;
  50.    while (l--) *elx++ -= *el++;
  51. }
  52.  
  53. void MatrixRowCol::Inject(const MatrixRowCol& mrc)
  54. // copy stored elements only
  55. {
  56.    REPORT
  57.    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  58.    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  59.    if (l<=0) return;
  60.    real* elx = store+f; real* ely = mrc.store+f;
  61.    while (l--) *elx++ = *ely++;
  62. }
  63.  
  64. real DotProd(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
  65. {
  66.    REPORT                                         // not accessed
  67.    int f = mrc1.skip; int f2 = mrc2.skip;
  68.    int l = f + mrc1.storage; int l2 = f2 + mrc2.storage;
  69.    if (f < f2) f = f2; if (l > l2) l = l2; l -= f;
  70.    if (l<=0) return 0.0;
  71.    
  72.    real* el1=mrc1.store+f; real* el2=mrc2.store+f;
  73.    real sum = 0.0;
  74.    while (l--) sum += *el1++ * *el2++;
  75.    return sum;
  76. }
  77.  
  78. void MatrixRowCol::Add(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
  79. {
  80.    int f = skip; int l = skip + storage;
  81.    int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
  82.    if (f1<f) f1=f; if (l1>l) l1=l;
  83.    int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
  84.    if (f2<f) f2=f; if (l2>l) l2=l;
  85.    real* el = store + f;
  86.    real* el1 = mrc1.store+f1; real* el2 = mrc2.store+f2;
  87.    if (f1<f2)
  88.    {
  89.       register int i = f1-f; while (i--) *el++ = 0.0;
  90.       if (l1<=f2)                              // disjoint
  91.       {
  92.      REPORT                                // not accessed
  93.          i = l1-f1;     while (i--) *el++ = *el1++;
  94.          i = f2-l1;     while (i--) *el++ = 0.0;
  95.          i = l2-f2;     while (i--) *el++ = *el2++;
  96.          i = l-l2;      while (i--) *el++ = 0.0;
  97.       }
  98.       else
  99.       {
  100.          i = f2-f1;    while (i--) *el++ = *el1++;
  101.          if (l1<=l2)
  102.          {
  103.             REPORT
  104.             i = l1-f2; while (i--) *el++ = *el1++ + *el2++;
  105.             i = l2-l1; while (i--) *el++ = *el2++;
  106.             i = l-l2;  while (i--) *el++ = 0.0;
  107.          }
  108.          else
  109.          {
  110.             REPORT
  111.             i = l2-f2; while (i--) *el++ = *el1++ + *el2++;
  112.             i = l1-l2; while (i--) *el++ = *el1++;
  113.             i = l-l1;  while (i--) *el++ = 0.0;
  114.          }
  115.       }
  116.    }
  117.    else
  118.    {
  119.       register int i = f2-f; while (i--) *el++ = 0.0;
  120.       if (l2<=f1)                              // disjoint
  121.       {
  122.      REPORT                                // not accessed
  123.          i = l2-f2;     while (i--) *el++ = *el2++;
  124.          i = f1-l2;     while (i--) *el++ = 0.0;
  125.      i = l1-f1;     while (i--) *el++ = *el1++;
  126.      i = l-l1;      while (i--) *el++ = 0.0;
  127.       }
  128.       else
  129.       {
  130.          i = f1-f2;    while (i--) *el++ = *el2++;
  131.          if (l2<=l1)
  132.          {
  133.         REPORT
  134.             i = l2-f1; while (i--) *el++ = *el1++ + *el2++;
  135.             i = l1-l2; while (i--) *el++ = *el1++;
  136.             i = l-l1;  while (i--) *el++ = 0.0;
  137.          }
  138.          else
  139.          {
  140.         REPORT
  141.             i = l1-f1; while (i--) *el++ = *el1++ + *el2++;
  142.             i = l2-l1; while (i--) *el++ = *el2++;
  143.             i = l-l2;  while (i--) *el++ = 0.0;
  144.          }
  145.       }
  146.    }
  147. }
  148.  
  149. void MatrixRowCol::Sub(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
  150. {
  151.    int f = skip; int l = skip + storage;
  152.    int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
  153.    if (f1<f) f1=f; if (l1>l) l1=l;
  154.    int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
  155.    if (f2<f) f2=f; if (l2>l) l2=l;
  156.    real* el = store + f;
  157.    real* el1 = mrc1.store+f1; real* el2 = mrc2.store+f2;
  158.    if (f1<f2)
  159.    {
  160.       register int i = f1-f; while (i--) *el++ = 0.0;
  161.       if (l1<=f2)                              // disjoint
  162.       {
  163.      REPORT                                // not accessed
  164.          i = l1-f1;     while (i--) *el++ = *el1++;
  165.          i = f2-l1;     while (i--) *el++ = 0.0;
  166.          i = l2-f2;     while (i--) *el++ = - *el2++;
  167.          i = l-l2;      while (i--) *el++ = 0.0;
  168.       }
  169.       else
  170.       {
  171.          i = f2-f1;    while (i--) *el++ = *el1++;
  172.          if (l1<=l2)
  173.          {
  174.         REPORT
  175.             i = l1-f2; while (i--) *el++ = *el1++ - *el2++;
  176.             i = l2-l1; while (i--) *el++ = - *el2++;
  177.             i = l-l2;  while (i--) *el++ = 0.0;
  178.          }
  179.          else
  180.          {
  181.         REPORT
  182.             i = l2-f2; while (i--) *el++ = *el1++ - *el2++;
  183.             i = l1-l2; while (i--) *el++ = *el1++;
  184.             i = l-l1;  while (i--) *el++ = 0.0;
  185.          }
  186.       }
  187.    }
  188.    else
  189.    {
  190.       register int i = f2-f; while (i--) *el++ = 0.0;
  191.       if (l2<=f1)                              // disjoint
  192.       {
  193.      REPORT                                // not accessed
  194.          i = l2-f2;     while (i--) *el++ = - *el2++;
  195.          i = f1-l2;     while (i--) *el++ = 0.0;
  196.          i = l1-f1;     while (i--) *el++ = *el1++;
  197.          i = l-l1;      while (i--) *el++ = 0.0;
  198.       }
  199.       else
  200.       {
  201.          i = f1-f2;    while (i--) *el++ = - *el2++;
  202.          if (l2<=l1)
  203.          {
  204.         REPORT
  205.             i = l2-f1; while (i--) *el++ = *el1++ - *el2++;
  206.             i = l1-l2; while (i--) *el++ = *el1++;
  207.             i = l-l1;  while (i--) *el++ = 0.0;
  208.          }
  209.          else
  210.          {
  211.             REPORT
  212.             i = l1-f1; while (i--) *el++ = *el1++ - *el2++;
  213.             i = l2-l1; while (i--) *el++ = - *el2++;
  214.             i = l-l2;  while (i--) *el++ = 0.0;
  215.          }
  216.       }
  217.    }
  218. }
  219.  
  220.  
  221. void MatrixRowCol::Add(const MatrixRowCol& mrc1, real x)
  222. {
  223.    REPORT
  224.    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  225.    if (f < skip) { f = skip; if (l < f) l = f; }
  226.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  227.  
  228.    real* elx = store+skip; real* ely = mrc1.store+f;
  229.  
  230.    int l1 = f-skip;  while (l1--) *elx++ = x;
  231.        l1 = l-f;     while (l1--) *elx++ = *ely++ + x;
  232.        lx -= l;      while (lx--) *elx++ = x;
  233. }
  234.  
  235. void MatrixRowCol::RevSub(const MatrixRowCol& mrc1)
  236. {
  237.    REPORT
  238.    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  239.    if (f < skip) { f = skip; if (l < f) l = f; }
  240.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  241.  
  242.    real* elx = store+skip; real* ely = mrc1.store+f;
  243.  
  244.    int l1 = f-skip;  while (l1--) { *elx = - *elx; elx++; }
  245.        l1 = l-f;     while (l1--) { *elx = *ely++ - *elx; elx++; }
  246.        lx -= l;      while (lx--) { *elx = - *elx; elx++; }
  247. }
  248.  
  249. void MatrixRowCol::Copy(const MatrixRowCol& mrc1)
  250. {
  251.    REPORT
  252.    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  253.    if (f < skip) { f = skip; if (l < f) l = f; }
  254.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  255.  
  256.    real* elx = store+skip; real* ely = mrc1.store+f;
  257.  
  258.    int l1 = f-skip;  while (l1--) *elx++ = 0.0;
  259.        l1 = l-f;     while (l1--) *elx++ = *ely++;
  260.        lx -= l;      while (lx--) *elx++ = 0.0;
  261. }
  262.  
  263. void MatrixRowCol::Negate(const MatrixRowCol& mrc1)
  264. {
  265.    REPORT
  266.    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  267.    if (f < skip) { f = skip; if (l < f) l = f; }
  268.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  269.  
  270.    real* elx = store+skip; real* ely = mrc1.store+f;
  271.  
  272.    int l1 = f-skip;  while (l1--) *elx++ = 0.0;
  273.        l1 = l-f;     while (l1--) *elx++ = - *ely++;
  274.        lx -= l;      while (lx--) *elx++ = 0.0;
  275. }
  276.  
  277. void MatrixRowCol::Multiply(const MatrixRowCol& mrc1, real s)
  278. {
  279.    REPORT
  280.    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  281.    if (f < skip) { f = skip; if (l < f) l = f; }
  282.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  283.  
  284.    real* elx = store+skip; real* ely = mrc1.store+f;
  285.  
  286.    int l1 = f-skip;  while (l1--) *elx++ = 0.0;
  287.        l1 = l-f;     while (l1--) *elx++ = *ely++ * s;
  288.        lx -= l;      while (lx--) *elx++ = 0.0;
  289. }
  290.  
  291. void DiagonalMatrix::Solver(MatrixRowCol& mrc, const MatrixRowCol& mrc1)
  292. {
  293.    REPORT
  294.    int f = mrc1.skip; int f0 = mrc.skip;
  295.    int l = f + mrc1.storage; int lx = f0 + mrc.storage;
  296.    if (f < f0) { f = f0; if (l < f) l = f; }
  297.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  298.  
  299.    real* elx = mrc.store+f0; real* eld = store+f;
  300.  
  301.    int l1 = f-f0;    while (l1--) *elx++ = 0.0;
  302.        l1 = l-f;     while (l1--) *elx++ /= *eld++;
  303.        lx -= l;      while (lx--) *elx++ = 0.0;
  304.    // Solver makes sure input and output point to same memory
  305. }
  306.  
  307. void MatrixRowCol::Copy(const real*& r)
  308. {
  309.    REPORT
  310.    real* elx = store+skip; const real* ely = r+skip; r += length;
  311.    int l = storage; while (l--) *elx++ = *ely++;
  312. }
  313.  
  314. void MatrixRowCol::Copy(real r)
  315. {
  316.    REPORT
  317.    real* elx = store+skip; int l = storage; while (l--) *elx++ = r;
  318. }
  319.  
  320. real MatrixRowCol::SumAbsoluteValue()
  321. {
  322.    REPORT
  323.    real sum = 0.0; real* elx = store+skip; int l = storage;
  324.    while (l--) sum += fabs(*elx++);
  325.    return sum;
  326. }
  327.  
  328. void MatrixRowCol::SubRowCol(MatrixRowCol& mrc, int skip1, int l1) const
  329. {
  330.    mrc.length = l1;  mrc.store = store + skip1;  int d = skip - skip1;
  331.    mrc.skip = (d < 0) ? 0 : d;  d = skip + storage - skip1;
  332.    d = ((l1 < d) ? l1 : d) - mrc.skip;  mrc.storage = (d < 0) ? 0 : d;
  333.    mrc.cw = 0;
  334. }
  335.  
  336. MatrixRowCol::~MatrixRowCol()
  337. {
  338.    if (cw*IsACopy && (cw*StoreHere==0))
  339.    {
  340.       MONITOR("Free    (RowCol)",storage,f) 
  341.       real* f = store+skip;  delete [storage] f;
  342.    }
  343. }
  344.  
  345. MatrixRow::~MatrixRow() { if (cw*StoreOnExit) gm->RestoreRow(*this); }
  346.  
  347. MatrixCol::~MatrixCol() { if (cw*StoreOnExit) gm->RestoreCol(*this); }
  348.  
  349.