home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #1 / NN_1993_1.iso / spool / comp / sources / misc / 4252 / newmat3.cxx < prev    next >
Encoding:
C/C++ Source or Header  |  1993-01-11  |  13.9 KB  |  506 lines
  1. //$$ newmat3.cxx        Matrix get and restore rows and columns
  2.  
  3. // Copyright (C) 1991,2,3: R B Davies
  4.  
  5.  
  6. #include "include.h"
  7.  
  8. #include "newmat.h"
  9. #include "newmatrc.h"
  10.  
  11. //#define REPORT { static ExeCounter ExeCount(__LINE__,3); ExeCount++; }
  12.  
  13. #define REPORT {}
  14.  
  15. //#define MONITOR(what,storage,store) \
  16. //   { cout << what << " " << storage << " at " << (long)store << "\n"; }
  17.  
  18. #define MONITOR(what,store,storage) {}
  19.  
  20.  
  21. // Control bits codes for GetRow, GetCol, RestoreRow, RestoreCol
  22. //
  23. // LoadOnEntry:
  24. //    Load data into MatrixRow or Col dummy array under GetRow or GetCol
  25. // StoreOnExit:
  26. //    Restore data to original matrix under RestoreRow or RestoreCol
  27. // IsACopy:
  28. //    Set by GetRow/Col: MatrixRow or Col array is a copy
  29. // DirectPart:
  30. //    Load or restore only part directly stored; must be set with StoreOnExit
  31. //    Still have decide  how to handle this with symmetric
  32. // StoreHere:
  33. //    used in columns only - store data at supplied storage address, adjusted
  34. //    for skip; used for GetCol, NextCol & RestoreCol. No need to fill out
  35. //    zeros.
  36.  
  37.  
  38. // These will work as a default
  39. // but need to override NextRow for efficiency
  40.  
  41. // Assume pointer arithmetic works for pointers out of range - not strict C++.
  42.  
  43.  
  44. void GeneralMatrix::NextRow(MatrixRowCol& mrc)
  45. {
  46.    REPORT
  47.    if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreRow(mrc); }
  48.    if (+(mrc.cw*IsACopy))
  49.    {
  50.       REPORT
  51.       Real* s = mrc.store + mrc.skip;
  52.       MONITOR_REAL_DELETE("Free   (NextRow)",mrc.storage,s)
  53. #ifdef Version21
  54.       delete [] s;
  55. #else
  56.       delete [mrc.storage] s;
  57. #endif
  58.    }
  59.    mrc.rowcol++;
  60.    if (mrc.rowcol<nrows) { REPORT this->GetRow(mrc); }
  61.    else { REPORT mrc.cw -= (StoreOnExit+IsACopy); }
  62. }
  63.  
  64. void GeneralMatrix::NextCol(MatrixRowCol& mrc)
  65. {
  66.    REPORT                                                // 423
  67.    if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreCol(mrc); }
  68.    int t1 = +(mrc.cw*IsACopy); int t2 = !(mrc.cw*StoreHere);
  69.    if ( t1 && t2 )
  70.    {
  71.       REPORT                                             // not accessed
  72.       Real* s = mrc.store + mrc.skip;
  73.       MONITOR_REAL_DELETE("Free   (NextCol)",mrc.storage,s) 
  74. #ifdef Version21
  75.       delete [] s;
  76. #else
  77.       delete [mrc.storage] s;
  78. #endif
  79.    }
  80.    mrc.rowcol++;
  81.    if (mrc.rowcol<ncols) { REPORT this->GetCol(mrc); }
  82.    else { REPORT mrc.cw -= (StoreOnExit+IsACopy); }
  83. }
  84.  
  85.  
  86. // routines for matrix
  87.  
  88. void Matrix::GetRow(MatrixRowCol& mrc)
  89. {
  90.    REPORT
  91.    mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=ncols;
  92.    mrc.store=store+mrc.rowcol*ncols;
  93. }
  94.  
  95.  
  96. void Matrix::GetCol(MatrixRowCol& mrc)
  97. {
  98.    REPORT
  99.    mrc.skip=0; mrc.storage=nrows; int t1 = !(mrc.cw*StoreHere);
  100.    if ( ncols==1 && t1 )
  101.       { REPORT mrc.cw-=IsACopy; mrc.store=store; }           // not accessed
  102.    else
  103.    {
  104.       mrc.cw+=IsACopy; Real* ColCopy;
  105.       if ( !(mrc.cw*StoreHere) )
  106.       {
  107.          REPORT
  108.          ColCopy = new Real [nrows]; MatrixErrorNoSpace(ColCopy);
  109.          MONITOR_REAL_NEW("Make (MatGetCol)",nrows,ColCopy)
  110.          mrc.store = ColCopy;
  111.       }
  112.       else { REPORT ColCopy = mrc.store; }
  113.       if (+(mrc.cw*LoadOnEntry))
  114.       {
  115.          REPORT
  116.          Real* Mstore = store+mrc.rowcol; int i=nrows;
  117.          while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
  118.       }
  119.    }
  120. }
  121.  
  122. void Matrix::RestoreCol(MatrixRowCol& mrc)
  123. {
  124. //  if (mrc.cw*StoreOnExit)
  125.    REPORT                                   // 429
  126.    if (+(mrc.cw*IsACopy))
  127.    {
  128.       REPORT                                // 426
  129.       Real* Mstore = store+mrc.rowcol; int i=nrows; Real* Cstore = mrc.store;
  130.       while (i--) { *Mstore = *Cstore++; Mstore+=ncols; }
  131.    }
  132. }
  133.  
  134. void Matrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrMat(); }  // 1808
  135.  
  136. void Matrix::NextCol(MatrixRowCol& mrc)
  137. {
  138.    REPORT                                        // 632
  139.    if (+(mrc.cw*StoreOnExit)) { REPORT RestoreCol(mrc); }
  140.    mrc.rowcol++;
  141.    if (mrc.rowcol<ncols)
  142.    {
  143.       if (+(mrc.cw*LoadOnEntry))
  144.       {
  145.      REPORT
  146.          Real* ColCopy = mrc.store;
  147.          Real* Mstore = store+mrc.rowcol; int i=nrows;
  148.          while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
  149.       }
  150.    }
  151.    else { REPORT mrc.cw -= StoreOnExit; }
  152. }
  153.  
  154. // routines for diagonal matrix
  155.  
  156. void DiagonalMatrix::GetRow(MatrixRowCol& mrc)
  157. {
  158.    REPORT
  159.    mrc.skip=mrc.rowcol; mrc.cw-=IsACopy; mrc.storage=1; mrc.store=store;
  160. }
  161.  
  162. void DiagonalMatrix::GetCol(MatrixRowCol& mrc)
  163. {
  164.    REPORT 
  165.    mrc.skip=mrc.rowcol; mrc.storage=1;
  166.    if (+(mrc.cw*StoreHere))
  167.       { REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); mrc.cw+=IsACopy; }
  168.    else { REPORT mrc.store = store; mrc.cw-=IsACopy; }     // not accessed
  169. }
  170.  
  171. void DiagonalMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
  172.                               // 800
  173. void DiagonalMatrix::NextCol(MatrixRowCol& mrc)
  174. {
  175.    REPORT
  176.    if (+(mrc.cw*StoreHere))
  177.    {
  178.       if (+(mrc.cw*StoreOnExit))
  179.          { REPORT *(store+mrc.rowcol)=*(mrc.store+mrc.rowcol); }
  180.       mrc.IncrDiag();
  181.       int t1 = +(mrc.cw*LoadOnEntry);
  182.       if (t1 && mrc.rowcol < ncols)
  183.          { REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); }
  184.    }
  185.    else { REPORT mrc.IncrDiag(); }                     // not accessed
  186. }
  187.  
  188. // routines for upper triangular matrix
  189.  
  190. void UpperTriangularMatrix::GetRow(MatrixRowCol& mrc)
  191. {
  192.    REPORT
  193.    int row = mrc.rowcol; mrc.skip=row; mrc.cw-=IsACopy;
  194.    mrc.storage=ncols-row; mrc.store=store+(row*(2*ncols-row-1))/2;
  195. }
  196.  
  197.  
  198. void UpperTriangularMatrix::GetCol(MatrixRowCol& mrc)
  199. {
  200.    REPORT
  201.    mrc.skip=0; mrc.cw+=IsACopy; int i=mrc.rowcol+1; mrc.storage=i;
  202.    Real* ColCopy;
  203.    if ( !(mrc.cw*StoreHere) )
  204.    {
  205.       REPORT                                              // not accessed
  206.       ColCopy = new Real [i]; MatrixErrorNoSpace(ColCopy);
  207.       MONITOR_REAL_NEW("Make (UT GetCol)",i,ColCopy) 
  208.       mrc.store = ColCopy;
  209.    }
  210.    else { REPORT ColCopy = mrc.store; }
  211.    if (+(mrc.cw*LoadOnEntry))
  212.    {
  213.       REPORT
  214.       Real* Mstore = store+mrc.rowcol; int j = ncols;
  215.       while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
  216.    }
  217. }
  218.  
  219. void UpperTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
  220. {
  221. //  if (mrc.cw*StoreOnExit)
  222.   {
  223.      REPORT
  224.      Real* Mstore = store+mrc.rowcol; int i=mrc.rowcol+1; int j = ncols;
  225.      Real* Cstore = mrc.store;
  226.      while (i--) { *Mstore = *Cstore++; Mstore += --j; }
  227.   }
  228. }
  229.  
  230. void UpperTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrUT(); }
  231.                               // 722
  232.  
  233. // routines for lower triangular matrix
  234.  
  235. void LowerTriangularMatrix::GetRow(MatrixRowCol& mrc)
  236. {
  237.    REPORT
  238.    int row=mrc.rowcol; mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=row+1;
  239.    mrc.store=store+(row*(row+1))/2;
  240. }
  241.  
  242. void LowerTriangularMatrix::GetCol(MatrixRowCol& mrc)
  243. {
  244.    REPORT
  245.    int col=mrc.rowcol; mrc.skip=col; mrc.cw+=IsACopy;
  246.    int i=nrows-col; mrc.storage=i; Real* ColCopy;
  247.    if ( !(mrc.cw*StoreHere) )
  248.    {
  249.       REPORT                                            // not accessed
  250.       ColCopy = new Real [i]; MatrixErrorNoSpace(ColCopy);
  251.       MONITOR_REAL_NEW("Make (LT GetCol)",i,ColCopy) 
  252.       mrc.store = ColCopy-col;
  253.    }
  254.    else { REPORT ColCopy = mrc.store+col; }
  255.    if (+(mrc.cw*LoadOnEntry))
  256.    {
  257.       REPORT
  258.       Real* Mstore = store+(col*(col+3))/2;
  259.       while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
  260.    }
  261. }
  262.  
  263. void LowerTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
  264. {
  265. //  if (mrc.cw*StoreOnExit)
  266.    {
  267.       REPORT
  268.       int col=mrc.rowcol; Real* Cstore = mrc.store+col;
  269.       Real* Mstore = store+(col*(col+3))/2; int i=nrows-col;
  270.       while (i--) { *Mstore = *Cstore++; Mstore += ++col; }
  271.    }
  272. }
  273.  
  274. void LowerTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrLT(); }
  275.                                      //712
  276. // routines for symmetric matrix
  277.  
  278. void SymmetricMatrix::GetRow(MatrixRowCol& mrc)
  279. {
  280.    REPORT                                                //571
  281.    mrc.skip=0; int row=mrc.rowcol;
  282.    if (+(mrc.cw*DirectPart))
  283.    {
  284.       REPORT
  285.       mrc.cw-=IsACopy; mrc.storage=row+1; mrc.store=store+(row*(row+1))/2;
  286.    }
  287.    else
  288.    {
  289.       mrc.cw+=IsACopy; mrc.storage=ncols;
  290.       Real* RowCopy = new Real [ncols]; MatrixErrorNoSpace(RowCopy);
  291.       MONITOR_REAL_NEW("Make (SymGetRow)",ncols,RowCopy) 
  292.       mrc.store = RowCopy;
  293.       if (+(mrc.cw*LoadOnEntry))
  294.       {
  295.      REPORT                                         // 544
  296.          Real* Mstore = store+(row*(row+1))/2; int i = row;
  297.          while (i--) *RowCopy++ = *Mstore++;
  298.          i = ncols-row;
  299.      while (i--) { *RowCopy++ = *Mstore; Mstore += ++row; }
  300.       }
  301.    }
  302. }
  303.  
  304. // need to check this out under StoreHere
  305.  
  306. void SymmetricMatrix::GetCol(MatrixRowCol& mrc)
  307. {
  308.    REPORT
  309.    mrc.skip=0; int col=mrc.rowcol;
  310.    if (+(mrc.cw*DirectPart))
  311.    {
  312.       REPORT                                         // not accessed
  313.       mrc.cw-=IsACopy; mrc.storage=col+1; mrc.store=store+(col*(col+1))/2;
  314.    }
  315.    else
  316.    {
  317.       mrc.cw+=IsACopy; mrc.storage=ncols; Real* ColCopy;
  318.       if ( !(mrc.cw*StoreHere) )
  319.       {
  320.          REPORT                                      // not accessed
  321.          ColCopy = new Real [ncols]; MatrixErrorNoSpace(ColCopy);
  322.          MONITOR_REAL_NEW("Make (SymGetCol)",ncols,ColCopy) 
  323.          mrc.store = ColCopy;
  324.       }
  325.       else { REPORT ColCopy = mrc.store; }
  326.       if (+(mrc.cw*LoadOnEntry))
  327.       {
  328.          REPORT
  329.          Real* Mstore = store+(col*(col+1))/2; int i = col;
  330.          while (i--) *ColCopy++ = *Mstore++;
  331.          i = ncols-col;
  332.      while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
  333.       }
  334.    }
  335. }
  336.  
  337. //void SymmetricMatrix::RestoreRow(int row, Real* Rstore)
  338. //{
  339. ////   if (cw*IsACopy && cw*StoreOnExit)
  340. //   {
  341. //      Real* Mstore = store+(row*(row+1))/2; int i = row+1;
  342. //      while (i--) *Mstore++ = *Rstore++;
  343. //   }
  344. //}
  345.  
  346. //void SymmetricMatrix::RestoreCol(int col, Real* Cstore)
  347. //{
  348. ////   if (cw*IsACopy && cw*StoreOnExit)
  349. //   {
  350. //      Real* Mstore = store+(col*(col+3))/2;
  351. //      int i = nrows-col; int j = col;
  352. //      while (i--) { *Mstore = *Cstore++; Mstore+= ++j; }
  353. //   }
  354. //}
  355.  
  356. // routines for row vector
  357.  
  358. void RowVector::GetCol(MatrixRowCol& mrc)
  359. {
  360.    REPORT 
  361.    mrc.skip=0; mrc.storage=1;
  362.    if (mrc.cw >= StoreHere)
  363.    {
  364.       if (mrc.cw >= LoadOnEntry) { REPORT *(mrc.store) = *(store+mrc.rowcol); }
  365.       mrc.cw+=IsACopy;
  366.    }
  367.    else  { REPORT mrc.store = store+mrc.rowcol; mrc.cw-=IsACopy; }
  368.                                                          // not accessed
  369. }
  370.  
  371. void RowVector::NextCol(MatrixRowCol& mrc) 
  372. {
  373.    REPORT
  374.    if (+(mrc.cw*StoreHere))
  375.    {
  376.       if (+(mrc.cw*StoreOnExit)) { REPORT *(store+mrc.rowcol)=*(mrc.store); }
  377.                              // not accessed
  378.       mrc.rowcol++;
  379.       if (mrc.rowcol < ncols)
  380.       {
  381.      if (+(mrc.cw*LoadOnEntry)) { REPORT *(mrc.store)=*(store+mrc.rowcol); }
  382.       }
  383.       else { REPORT mrc.cw -= StoreOnExit; }
  384.    }
  385.    else  { REPORT mrc.rowcol++; mrc.store++; }             // not accessed
  386. }
  387.  
  388. void RowVector::RestoreCol(MatrixRowCol& mrc)
  389. {
  390.    REPORT                                            // not accessed
  391.    if (mrc.cw>=IsACopy)  { REPORT *(store+mrc.rowcol)=*(mrc.store); }
  392. }
  393.  
  394.  
  395. // routines for band matrices
  396.  
  397. void BandMatrix::GetRow(MatrixRowCol& mrc)
  398. {
  399.    REPORT
  400.    mrc.cw -= IsACopy; int r = mrc.rowcol; int w = lower+1+upper;
  401.    int s = r-lower; mrc.store = store+(r*w-s); if (s<0) { w += s; s = 0; }
  402.    mrc.skip = s; s += w-ncols; if (s>0) w -= s; mrc.storage = w;
  403. }
  404.  
  405. // make special versions of this for upper and lower band matrices
  406.  
  407. void BandMatrix::NextRow(MatrixRowCol& mrc)
  408. {
  409.    REPORT
  410.    int r = ++mrc.rowcol; mrc.store += lower+upper;
  411.    if (r<=lower) mrc.storage++; else mrc.skip++;
  412.    if (r>=ncols-upper) mrc.storage--;
  413. }
  414.  
  415. void BandMatrix::GetCol(MatrixRowCol& mrc)
  416. {
  417.    REPORT
  418.    mrc.cw += IsACopy; int c = mrc.rowcol; int n = lower+upper; int w = n+1;
  419.    int b; int s = c-upper; Real* ColCopy;
  420.    if (s<=0) { w += s; s = 0; b = c+lower; } else b = s*w+n;
  421.    mrc.skip = s; s += w-nrows; if (s>0) w -= s; mrc.storage = w;
  422.    if ( !(mrc.cw*StoreHere) )
  423.    {
  424.       REPORT
  425.       ColCopy = new Real [w]; MatrixErrorNoSpace(ColCopy);
  426.       MONITOR_REAL_NEW("Make (BMGetCol)",w,ColCopy)
  427.       mrc.store = ColCopy-mrc.skip;
  428.    }
  429.    else { REPORT ColCopy = mrc.store+mrc.skip; }
  430.    if (+(mrc.cw*LoadOnEntry))
  431.    {
  432.       REPORT
  433.       Real* Mstore = store+b;
  434.       while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
  435.    }
  436. }
  437.  
  438. void BandMatrix::RestoreCol(MatrixRowCol& mrc)
  439. {
  440. //  if (mrc.cw*StoreOnExit)
  441.    REPORT
  442.    int c = mrc.rowcol; int n = lower+upper; int s = c-upper;
  443.    Real* Mstore = store + ((s<=0) ? c+lower : s*n+s+n);
  444.    Real* Cstore = mrc.store+mrc.skip; int w = mrc.storage;
  445.    while (w--) { *Mstore = *Cstore++; Mstore += n; }
  446. }
  447.  
  448. // routines for symmetric band matrix
  449.  
  450. void SymmetricBandMatrix::GetRow(MatrixRowCol& mrc)
  451. {
  452.    REPORT
  453.    int r=mrc.rowcol; int s = r-lower; int w1 = lower+1; int o = r*w1;
  454.    if (s<0) { w1 += s; o -= s; s = 0; }  mrc.skip = s;
  455.  
  456.    if (+(mrc.cw*DirectPart))
  457.       { REPORT  mrc.cw -= IsACopy; mrc.store = store+o-s; mrc.storage = w1; }
  458.    else
  459.    {
  460.       mrc.cw += IsACopy; int w = w1+lower; s += w-ncols;
  461.       if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
  462.       Real* RowCopy = new Real [w]; MatrixErrorNoSpace(RowCopy);
  463.       MONITOR_REAL_NEW("Make (SmBGetRow)",w,RowCopy) 
  464.       mrc.store = RowCopy-mrc.skip;
  465.       if (+(mrc.cw*LoadOnEntry))
  466.       {
  467.      REPORT
  468.          Real* Mstore = store+o;
  469.          while (w1--) *RowCopy++ = *Mstore++;   Mstore--;
  470.          while (w2--) { Mstore += lower; *RowCopy++ = *Mstore; }
  471.       }
  472.    }
  473. }
  474.  
  475. // need to check this out under StoreHere
  476.  
  477. void SymmetricBandMatrix::GetCol(MatrixRowCol& mrc)
  478. {
  479.    REPORT
  480.    int c=mrc.rowcol; int s = c-lower; int w1 = lower+1; int o = c*w1;
  481.    if (s<0) { w1 += s; o -= s; s = 0; }  mrc.skip = s;
  482.  
  483.    if (+(mrc.cw*DirectPart))
  484.       { REPORT  mrc.cw -= IsACopy; mrc.store = store+o-s; mrc.storage = w1; }
  485.    else
  486.    {
  487.       mrc.cw += IsACopy; int w = w1+lower; s += w-ncols;
  488.       if (s>0) w -= s; mrc.storage = w; int w2 = w-w1; Real* ColCopy;
  489.       if ( !(mrc.cw*StoreHere) )
  490.       {
  491.          ColCopy = new Real [w]; MatrixErrorNoSpace(ColCopy);
  492.          MONITOR_REAL_NEW("Make (SmBGetCol)",w,ColCopy) 
  493.          mrc.store = ColCopy-mrc.skip;
  494.       }
  495.       else { REPORT ColCopy = mrc.store+mrc.skip; }
  496.       if (+(mrc.cw*LoadOnEntry))
  497.       {
  498.      REPORT
  499.          Real* Mstore = store+o;
  500.          while (w1--) *ColCopy++ = *Mstore++;   Mstore--;
  501.          while (w2--) { Mstore += lower; *ColCopy++ = *Mstore; }
  502.       }
  503.    }
  504. }
  505.  
  506.