home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / World_Of_Computer_Software-02-385-Vol-1of3.iso / n / newmat06.zip / NEWMAT3.CXX < prev    next >
C/C++ Source or Header  |  1992-07-14  |  15KB  |  504 lines
  1. //$$ newmat3.cxx        Matrix get and restore rows and columns
  2.  
  3. // Copyright (C) 1991,2: 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.    if ( +(mrc.cw*IsACopy) && (!(mrc.cw*StoreHere)) )
  69.    {
  70.       REPORT                                             // not accessed
  71.       Real* s = mrc.store + mrc.skip;
  72.       MONITOR_REAL_DELETE("Free   (NextCol)",mrc.storage,s) 
  73. #ifdef Version21
  74.       delete [] s;
  75. #else
  76.       delete [mrc.storage] s;
  77. #endif
  78.    }
  79.    mrc.rowcol++;
  80.    if (mrc.rowcol<ncols) { REPORT this->GetCol(mrc); }
  81.    else { REPORT mrc.cw -= (StoreOnExit+IsACopy); }
  82. }
  83.  
  84.  
  85. // routines for matrix
  86.  
  87. void Matrix::GetRow(MatrixRowCol& mrc)
  88. {
  89.    REPORT
  90.    mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=ncols;
  91.    mrc.store=store+mrc.rowcol*ncols;
  92. }
  93.  
  94.  
  95. void Matrix::GetCol(MatrixRowCol& mrc)
  96. {
  97.    REPORT
  98.    mrc.skip=0; mrc.storage=nrows;
  99.    if ( ncols==1 && !(mrc.cw*StoreHere) )
  100.       { REPORT mrc.cw-=IsACopy; mrc.store=store; }           // not accessed
  101.    else
  102.    {
  103.       mrc.cw+=IsACopy; Real* ColCopy;
  104.       if ( !(mrc.cw*StoreHere) )
  105.       {
  106.          REPORT
  107.          ColCopy = new Real [nrows]; MatrixErrorNoSpace(ColCopy);
  108.          MONITOR_REAL_NEW("Make (MatGetCol)",nrows,ColCopy)
  109.          mrc.store = ColCopy;
  110.       }
  111.       else { REPORT ColCopy = mrc.store; }
  112.       if (+(mrc.cw*LoadOnEntry))
  113.       {
  114.          REPORT
  115.          Real* Mstore = store+mrc.rowcol; int i=nrows;
  116.          while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
  117.       }
  118.    }
  119. }
  120.  
  121. void Matrix::RestoreCol(MatrixRowCol& mrc)
  122. {
  123. //  if (mrc.cw*StoreOnExit)
  124.    REPORT                                   // 429
  125.    if (+(mrc.cw*IsACopy))
  126.    {
  127.       REPORT                                // 426
  128.       Real* Mstore = store+mrc.rowcol; int i=nrows; Real* Cstore = mrc.store;
  129.       while (i--) { *Mstore = *Cstore++; Mstore+=ncols; }
  130.    }
  131. }
  132.  
  133. void Matrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrMat(); }  // 1808
  134.  
  135. void Matrix::NextCol(MatrixRowCol& mrc)
  136. {
  137.    REPORT                                        // 632
  138.    if (+(mrc.cw*StoreOnExit)) { REPORT RestoreCol(mrc); }
  139.    mrc.rowcol++;
  140.    if (mrc.rowcol<ncols)
  141.    {
  142.       if (+(mrc.cw*LoadOnEntry))
  143.       {
  144.      REPORT
  145.          Real* ColCopy = mrc.store;
  146.          Real* Mstore = store+mrc.rowcol; int i=nrows;
  147.          while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
  148.       }
  149.    }
  150.    else { REPORT mrc.cw -= StoreOnExit; }
  151. }
  152.  
  153. // routines for diagonal matrix
  154.  
  155. void DiagonalMatrix::GetRow(MatrixRowCol& mrc)
  156. {
  157.    REPORT
  158.    mrc.skip=mrc.rowcol; mrc.cw-=IsACopy; mrc.storage=1; mrc.store=store;
  159. }
  160.  
  161. void DiagonalMatrix::GetCol(MatrixRowCol& mrc)
  162. {
  163.    REPORT 
  164.    mrc.skip=mrc.rowcol; mrc.storage=1;
  165.    if (+(mrc.cw*StoreHere))
  166.       { REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); mrc.cw+=IsACopy; }
  167.    else { REPORT mrc.store = store; mrc.cw-=IsACopy; }     // not accessed
  168. }
  169.  
  170. void DiagonalMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
  171.                               // 800
  172. void DiagonalMatrix::NextCol(MatrixRowCol& mrc)
  173. {
  174.    REPORT
  175.    if (+(mrc.cw*StoreHere))
  176.    {
  177.       if (+(mrc.cw*StoreOnExit))
  178.          { REPORT *(store+mrc.rowcol)=*(mrc.store+mrc.rowcol); }
  179.       mrc.IncrDiag();
  180.       if (+(mrc.cw*LoadOnEntry) && mrc.rowcol < ncols)
  181.          { REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); }
  182.    }
  183.    else { REPORT mrc.IncrDiag(); }                     // not accessed
  184. }
  185.  
  186. // routines for upper triangular matrix
  187.  
  188. void UpperTriangularMatrix::GetRow(MatrixRowCol& mrc)
  189. {
  190.    REPORT
  191.    int row = mrc.rowcol; mrc.skip=row; mrc.cw-=IsACopy;
  192.    mrc.storage=ncols-row; mrc.store=store+(row*(2*ncols-row-1))/2;
  193. }
  194.  
  195.  
  196. void UpperTriangularMatrix::GetCol(MatrixRowCol& mrc)
  197. {
  198.    REPORT
  199.    mrc.skip=0; mrc.cw+=IsACopy; int i=mrc.rowcol+1; mrc.storage=i;
  200.    Real* ColCopy;
  201.    if ( !(mrc.cw*StoreHere) )
  202.    {
  203.       REPORT                                              // not accessed
  204.       ColCopy = new Real [i]; MatrixErrorNoSpace(ColCopy);
  205.       MONITOR_REAL_NEW("Make (UT GetCol)",i,ColCopy) 
  206.       mrc.store = ColCopy;
  207.    }
  208.    else { REPORT ColCopy = mrc.store; }
  209.    if (+(mrc.cw*LoadOnEntry))
  210.    {
  211.       REPORT
  212.       Real* Mstore = store+mrc.rowcol; int j = ncols;
  213.       while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
  214.    }
  215. }
  216.  
  217. void UpperTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
  218. {
  219. //  if (mrc.cw*StoreOnExit)
  220.   {
  221.      REPORT
  222.      Real* Mstore = store+mrc.rowcol; int i=mrc.rowcol+1; int j = ncols;
  223.      Real* Cstore = mrc.store;
  224.      while (i--) { *Mstore = *Cstore++; Mstore += --j; }
  225.   }
  226. }
  227.  
  228. void UpperTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrUT(); }
  229.                               // 722
  230.  
  231. // routines for lower triangular matrix
  232.  
  233. void LowerTriangularMatrix::GetRow(MatrixRowCol& mrc)
  234. {
  235.    REPORT
  236.    int row=mrc.rowcol; mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=row+1;
  237.    mrc.store=store+(row*(row+1))/2;
  238. }
  239.  
  240. void LowerTriangularMatrix::GetCol(MatrixRowCol& mrc)
  241. {
  242.    REPORT
  243.    int col=mrc.rowcol; mrc.skip=col; mrc.cw+=IsACopy;
  244.    int i=nrows-col; mrc.storage=i; Real* ColCopy;
  245.    if ( !(mrc.cw*StoreHere) )
  246.    {
  247.       REPORT                                            // not accessed
  248.       ColCopy = new Real [i]; MatrixErrorNoSpace(ColCopy);
  249.       MONITOR_REAL_NEW("Make (LT GetCol)",i,ColCopy) 
  250.       mrc.store = ColCopy-col;
  251.    }
  252.    else { REPORT ColCopy = mrc.store+col; }
  253.    if (+(mrc.cw*LoadOnEntry))
  254.    {
  255.       REPORT
  256.       Real* Mstore = store+(col*(col+3))/2;
  257.       while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
  258.    }
  259. }
  260.  
  261. void LowerTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
  262. {
  263. //  if (mrc.cw*StoreOnExit)
  264.    {
  265.       REPORT
  266.       int col=mrc.rowcol; Real* Cstore = mrc.store+col;
  267.       Real* Mstore = store+(col*(col+3))/2; int i=nrows-col;
  268.       while (i--) { *Mstore = *Cstore++; Mstore += ++col; }
  269.    }
  270. }
  271.  
  272. void LowerTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrLT(); }
  273.                                      //712
  274. // routines for symmetric matrix
  275.  
  276. void SymmetricMatrix::GetRow(MatrixRowCol& mrc)
  277. {
  278.    REPORT                                                //571
  279.    mrc.skip=0; int row=mrc.rowcol;
  280.    if (+(mrc.cw*DirectPart))
  281.    {
  282.       REPORT
  283.       mrc.cw-=IsACopy; mrc.storage=row+1; mrc.store=store+(row*(row+1))/2;
  284.    }
  285.    else
  286.    {
  287.       mrc.cw+=IsACopy; mrc.storage=ncols;
  288.       Real* RowCopy = new Real [ncols]; MatrixErrorNoSpace(RowCopy);
  289.       MONITOR_REAL_NEW("Make (SymGetRow)",ncols,RowCopy) 
  290.       mrc.store = RowCopy;
  291.       if (+(mrc.cw*LoadOnEntry))
  292.       {
  293.      REPORT                                         // 544
  294.          Real* Mstore = store+(row*(row+1))/2; int i = row;
  295.          while (i--) *RowCopy++ = *Mstore++;
  296.          i = ncols-row;
  297.      while (i--) { *RowCopy++ = *Mstore; Mstore += ++row; }
  298.       }
  299.    }
  300. }
  301.  
  302. // need to check this out under StoreHere
  303.  
  304. void SymmetricMatrix::GetCol(MatrixRowCol& mrc)
  305. {
  306.    REPORT
  307.    mrc.skip=0; int col=mrc.rowcol;
  308.    if (+(mrc.cw*DirectPart))
  309.    {
  310.       REPORT                                         // not accessed
  311.       mrc.cw-=IsACopy; mrc.storage=col+1; mrc.store=store+(col*(col+1))/2;
  312.    }
  313.    else
  314.    {
  315.       mrc.cw+=IsACopy; mrc.storage=ncols; Real* ColCopy;
  316.       if ( !(mrc.cw*StoreHere) )
  317.       {
  318.          REPORT                                      // not accessed
  319.          ColCopy = new Real [ncols]; MatrixErrorNoSpace(ColCopy);
  320.          MONITOR_REAL_NEW("Make (SymGetCol)",ncols,ColCopy) 
  321.          mrc.store = ColCopy;
  322.       }
  323.       else { REPORT ColCopy = mrc.store; }
  324.       if (+(mrc.cw*LoadOnEntry))
  325.       {
  326.          REPORT
  327.          Real* Mstore = store+(col*(col+1))/2; int i = col;
  328.          while (i--) *ColCopy++ = *Mstore++;
  329.          i = ncols-col;
  330.      while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
  331.       }
  332.    }
  333. }
  334.  
  335. //void SymmetricMatrix::RestoreRow(int row, Real* Rstore)
  336. //{
  337. ////   if (cw*IsACopy && cw*StoreOnExit)
  338. //   {
  339. //      Real* Mstore = store+(row*(row+1))/2; int i = row+1;
  340. //      while (i--) *Mstore++ = *Rstore++;
  341. //   }
  342. //}
  343.  
  344. //void SymmetricMatrix::RestoreCol(int col, Real* Cstore)
  345. //{
  346. ////   if (cw*IsACopy && cw*StoreOnExit)
  347. //   {
  348. //      Real* Mstore = store+(col*(col+3))/2;
  349. //      int i = nrows-col; int j = col;
  350. //      while (i--) { *Mstore = *Cstore++; Mstore+= ++j; }
  351. //   }
  352. //}
  353.  
  354. // routines for row vector
  355.  
  356. void RowVector::GetCol(MatrixRowCol& mrc)
  357. {
  358.    REPORT 
  359.    mrc.skip=0; mrc.storage=1;
  360.    if (mrc.cw >= StoreHere)
  361.    {
  362.       if (mrc.cw >= LoadOnEntry) { REPORT *(mrc.store) = *(store+mrc.rowcol); }
  363.       mrc.cw+=IsACopy;
  364.    }
  365.    else  { REPORT mrc.store = store+mrc.rowcol; mrc.cw-=IsACopy; }
  366.                                                          // not accessed
  367. }
  368.  
  369. void RowVector::NextCol(MatrixRowCol& mrc) 
  370. {
  371.    REPORT
  372.    if (+(mrc.cw*StoreHere))
  373.    {
  374.       if (+(mrc.cw*StoreOnExit)) { REPORT *(store+mrc.rowcol)=*(mrc.store); }
  375.                              // not accessed
  376.       mrc.rowcol++;
  377.       if (mrc.rowcol < ncols)
  378.       {
  379.      if (+(mrc.cw*LoadOnEntry)) { REPORT *(mrc.store)=*(store+mrc.rowcol); }
  380.       }
  381.       else { REPORT mrc.cw -= StoreOnExit; }
  382.    }
  383.    else  { REPORT mrc.rowcol++; mrc.store++; }             // not accessed
  384. }
  385.  
  386. void RowVector::RestoreCol(MatrixRowCol& mrc)
  387. {
  388.    REPORT                                            // not accessed
  389.    if (mrc.cw>=IsACopy)  { REPORT *(store+mrc.rowcol)=*(mrc.store); }
  390. }
  391.  
  392.  
  393. // routines for band matrices
  394.  
  395. void BandMatrix::GetRow(MatrixRowCol& mrc)
  396. {
  397.    REPORT
  398.    mrc.cw -= IsACopy; int r = mrc.rowcol; int w = lower+1+upper;
  399.    int s = r-lower; mrc.store = store+(r*w-s); if (s<0) { w += s; s = 0; }
  400.    mrc.skip = s; s += w-ncols; if (s>0) w -= s; mrc.storage = w;
  401. }
  402.  
  403. // make special versions of this for upper and lower band matrices
  404.  
  405. void BandMatrix::NextRow(MatrixRowCol& mrc)
  406. {
  407.    REPORT
  408.    int r = ++mrc.rowcol; mrc.store += lower+upper;
  409.    if (r<=lower) mrc.storage++; else mrc.skip++;
  410.    if (r>=ncols-upper) mrc.storage--;
  411. }
  412.  
  413. void BandMatrix::GetCol(MatrixRowCol& mrc)
  414. {
  415.    REPORT
  416.    mrc.cw += IsACopy; int c = mrc.rowcol; int n = lower+upper; int w = n+1;
  417.    int b; int s = c-upper; Real* ColCopy;
  418.    if (s<=0) { w += s; s = 0; b = c+lower; } else b = s*w+n;
  419.    mrc.skip = s; s += w-nrows; if (s>0) w -= s; mrc.storage = w;
  420.    if ( !(mrc.cw*StoreHere) )
  421.    {
  422.       REPORT
  423.       ColCopy = new Real [w]; MatrixErrorNoSpace(ColCopy);
  424.       MONITOR_REAL_NEW("Make (BMGetCol)",w,ColCopy)
  425.       mrc.store = ColCopy-mrc.skip;
  426.    }
  427.    else { REPORT ColCopy = mrc.store+mrc.skip; }
  428.    if (+(mrc.cw*LoadOnEntry))
  429.    {
  430.       REPORT
  431.       Real* Mstore = store+b;
  432.       while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
  433.    }
  434. }
  435.  
  436. void BandMatrix::RestoreCol(MatrixRowCol& mrc)
  437. {
  438. //  if (mrc.cw*StoreOnExit)
  439.    REPORT
  440.    int c = mrc.rowcol; int n = lower+upper; int s = c-upper;
  441.    Real* Mstore = store + ((s<=0) ? c+lower : s*n+s+n);
  442.    Real* Cstore = mrc.store+mrc.skip; int w = mrc.storage;
  443.    while (w--) { *Mstore = *Cstore++; Mstore += n; }
  444. }
  445.  
  446. // routines for symmetric band matrix
  447.  
  448. void SymmetricBandMatrix::GetRow(MatrixRowCol& mrc)
  449. {
  450.    REPORT
  451.    int r=mrc.rowcol; int s = r-lower; int w1 = lower+1; int o = r*w1;
  452.    if (s<0) { w1 += s; o -= s; s = 0; }  mrc.skip = s;
  453.  
  454.    if (+(mrc.cw*DirectPart))
  455.       { REPORT  mrc.cw -= IsACopy; mrc.store = store+o-s; mrc.storage = w1; }
  456.    else
  457.    {
  458.       mrc.cw += IsACopy; int w = w1+lower; s += w-ncols;
  459.       if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
  460.       Real* RowCopy = new Real [w]; MatrixErrorNoSpace(RowCopy);
  461.       MONITOR_REAL_NEW("Make (SmBGetRow)",w,RowCopy) 
  462.       mrc.store = RowCopy-mrc.skip;
  463.       if (+(mrc.cw*LoadOnEntry))
  464.       {
  465.      REPORT
  466.          Real* Mstore = store+o;
  467.          while (w1--) *RowCopy++ = *Mstore++;   Mstore--;
  468.          while (w2--) { Mstore += lower; *RowCopy++ = *Mstore; }
  469.       }
  470.    }
  471. }
  472.  
  473. // need to check this out under StoreHere
  474.  
  475. void SymmetricBandMatrix::GetCol(MatrixRowCol& mrc)
  476. {
  477.    REPORT
  478.    int c=mrc.rowcol; int s = c-lower; int w1 = lower+1; int o = c*w1;
  479.    if (s<0) { w1 += s; o -= s; s = 0; }  mrc.skip = s;
  480.  
  481.    if (+(mrc.cw*DirectPart))
  482.       { REPORT  mrc.cw -= IsACopy; mrc.store = store+o-s; mrc.storage = w1; }
  483.    else
  484.    {
  485.       mrc.cw += IsACopy; int w = w1+lower; s += w-ncols;
  486.       if (s>0) w -= s; mrc.storage = w; int w2 = w-w1; Real* ColCopy;
  487.       if ( !(mrc.cw*StoreHere) )
  488.       {
  489.          ColCopy = new Real [w]; MatrixErrorNoSpace(ColCopy);
  490.          MONITOR_REAL_NEW("Make (SmBGetCol)",w,ColCopy) 
  491.          mrc.store = ColCopy-mrc.skip;
  492.       }
  493.       else { REPORT ColCopy = mrc.store+mrc.skip; }
  494.       if (+(mrc.cw*LoadOnEntry))
  495.       {
  496.      REPORT
  497.          Real* Mstore = store+o;
  498.          while (w1--) *ColCopy++ = *Mstore++;   Mstore--;
  499.          while (w2--) { Mstore += lower; *ColCopy++ = *Mstore; }
  500.       }
  501.    }
  502. }
  503.  
  504.