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