home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume26 / newmat04 / part03 / newmat3.cxx < prev    next >
Encoding:
C/C++ Source or Header  |  1991-11-30  |  10.4 KB  |  381 lines
  1. //$$ newmat3.cxx        Matrix get and restore rows and columns
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5.  
  6. #include "include.hxx"
  7.  
  8. #include "newmat.hxx"
  9. #include "newmatrc.hxx"
  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.  
  42. void GeneralMatrix::NextRow(MatrixRowCol& mrc)
  43. {
  44.    REPORT
  45.    if (mrc.cw*StoreOnExit) { REPORT this->RestoreRow(mrc); }
  46.    if (mrc.cw*IsACopy)
  47.    {
  48.       REPORT
  49.       real* s = mrc.store + mrc.skip;
  50.       MONITOR("Free   (NextRow)",mrc.storage,s) 
  51.       delete [mrc.storage] s;
  52.    }
  53.    mrc.rowcol++;
  54.    if (mrc.rowcol<nrows) { REPORT this->GetRow(mrc); }
  55.    else { REPORT mrc.cw -= (StoreOnExit+IsACopy); }
  56. }
  57.  
  58. void GeneralMatrix::NextCol(MatrixRowCol& mrc)
  59. {
  60.    REPORT                                                // 423
  61.    if (mrc.cw*StoreOnExit) { REPORT this->RestoreCol(mrc); }
  62.    if (mrc.cw*IsACopy && (mrc.cw*StoreHere==0))
  63.    {
  64.       REPORT                                             // not accessed
  65.       real* s = mrc.store + mrc.skip;
  66.       MONITOR("Free   (NextCol)",mrc.storage,s) 
  67.       delete [mrc.storage] s;
  68.    }
  69.    mrc.rowcol++;
  70.    if (mrc.rowcol<ncols) { REPORT this->GetCol(mrc); }
  71.    else { REPORT mrc.cw -= (StoreOnExit+IsACopy); }
  72. }
  73.  
  74.  
  75. // routines for matrix
  76.  
  77. void Matrix::GetRow(MatrixRowCol& mrc)
  78. {
  79.    REPORT
  80.    mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=ncols;
  81.    mrc.store=store+mrc.rowcol*ncols;
  82. }
  83.  
  84.  
  85. void Matrix::GetCol(MatrixRowCol& mrc)
  86. {
  87.    REPORT
  88.    mrc.skip=0; mrc.storage=nrows;
  89.    if (ncols==1 && (mrc.cw*StoreHere)==0)
  90.       { REPORT mrc.cw-=IsACopy; mrc.store=store; }           // not accessed
  91.    else
  92.    {
  93.       mrc.cw+=IsACopy; real* ColCopy;
  94.       if ((mrc.cw*StoreHere)==0)
  95.       {
  96.          REPORT
  97.          ColCopy = new real [nrows]; MatrixErrorNoSpace(ColCopy);
  98.          MONITOR("Make (MatGetCol)",nrows,ColCopy)
  99.          mrc.store = ColCopy;
  100.       }
  101.       else { REPORT ColCopy = mrc.store; }
  102.       if (mrc.cw*LoadOnEntry)
  103.       {
  104.          REPORT
  105.          real* Mstore = store+mrc.rowcol; int i=nrows;
  106.          while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
  107.       }
  108.    }
  109. }
  110.  
  111. void Matrix::RestoreCol(MatrixRowCol& mrc)
  112. {
  113. //  if (mrc.cw*StoreOnExit)
  114.    REPORT                                   // 429
  115.    if (mrc.cw*IsACopy)
  116.    {
  117.       REPORT                                // 426
  118.       real* Mstore = store+mrc.rowcol; int i=nrows; real* Cstore = mrc.store;
  119.       while (i--) { *Mstore = *Cstore++; Mstore+=ncols; }
  120.    }
  121. }
  122.  
  123. void Matrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrMat(); }  // 1808
  124.  
  125. void Matrix::NextCol(MatrixRowCol& mrc)
  126. {
  127.    REPORT                                        // 632
  128.    if (mrc.cw*StoreOnExit) { REPORT RestoreCol(mrc); }
  129.    mrc.rowcol++;
  130.    if (mrc.rowcol<ncols)
  131.    {
  132.       if (mrc.cw*LoadOnEntry)
  133.       {
  134.          REPORT
  135.          real* ColCopy = mrc.store;
  136.          real* Mstore = store+mrc.rowcol; int i=nrows;
  137.          while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
  138.       }
  139.    }
  140.    else { REPORT mrc.cw -= StoreOnExit; }
  141. }
  142.  
  143. // routines for diagonal matrix
  144.  
  145. void DiagonalMatrix::GetRow(MatrixRowCol& mrc)
  146. {
  147.    REPORT
  148.    mrc.skip=mrc.rowcol; mrc.cw-=IsACopy; mrc.storage=1; mrc.store=store;
  149. }
  150.  
  151. void DiagonalMatrix::GetCol(MatrixRowCol& mrc)
  152. {
  153.    REPORT 
  154.    mrc.skip=mrc.rowcol; mrc.storage=1;
  155.    if (mrc.cw*StoreHere) 
  156.       { REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); mrc.cw+=IsACopy; }
  157.    else { REPORT mrc.store = store; mrc.cw-=IsACopy; }     // not accessed
  158. }
  159.  
  160. void DiagonalMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
  161.                               // 800
  162. void DiagonalMatrix::NextCol(MatrixRowCol& mrc)
  163. {
  164.    REPORT
  165.    if (mrc.cw*StoreHere)
  166.    {
  167.       if (mrc.cw*StoreOnExit) 
  168.          { REPORT *(store+mrc.rowcol)=*(mrc.store+mrc.rowcol); }
  169.       mrc.IncrDiag();
  170.       if (mrc.cw*LoadOnEntry && mrc.rowcol < ncols)
  171.          { REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); }
  172.    }
  173.    else { REPORT mrc.IncrDiag(); }                     // not accessed
  174. }
  175.  
  176. // routines for upper triangular matrix
  177.  
  178. void UpperTriangularMatrix::GetRow(MatrixRowCol& mrc)
  179. {
  180.    REPORT
  181.    int row = mrc.rowcol; mrc.skip=row; mrc.cw-=IsACopy;
  182.    mrc.storage=ncols-row; mrc.store=store+(row*(2*ncols-row-1))/2;
  183. }
  184.  
  185.  
  186. void UpperTriangularMatrix::GetCol(MatrixRowCol& mrc)
  187. {
  188.    REPORT
  189.    mrc.skip=0; mrc.cw+=IsACopy; int i=mrc.rowcol+1; mrc.storage=i;
  190.    real* ColCopy;
  191.    if ((mrc.cw*StoreHere)==0)
  192.    {
  193.       REPORT                                              // not accessed
  194.       ColCopy = new real [i]; MatrixErrorNoSpace(ColCopy);
  195.       MONITOR("Make (UT GetCol)",i,ColCopy) 
  196.       mrc.store = ColCopy;
  197.    }
  198.    else { REPORT ColCopy = mrc.store; }
  199.    if (mrc.cw*LoadOnEntry)
  200.    {
  201.       REPORT
  202.       real* Mstore = store+mrc.rowcol; int j = ncols;
  203.       while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
  204.    }
  205. }
  206.  
  207. void UpperTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
  208. {
  209. //  if (mrc.cw*StoreOnExit)
  210.   {
  211.      REPORT
  212.      real* Mstore = store+mrc.rowcol; int i=mrc.rowcol+1; int j = ncols;
  213.      real* Cstore = mrc.store;
  214.      while (i--) { *Mstore = *Cstore++; Mstore += --j; }
  215.   }
  216. }
  217.  
  218. void UpperTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrUT(); }
  219.                               // 722
  220.  
  221. // routines for lower triangular matrix
  222.  
  223. void LowerTriangularMatrix::GetRow(MatrixRowCol& mrc)
  224. {
  225.    REPORT
  226.    int row=mrc.rowcol; mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=row+1;
  227.    mrc.store=store+(row*(row+1))/2;
  228. }
  229.  
  230. void LowerTriangularMatrix::GetCol(MatrixRowCol& mrc)
  231. {
  232.    REPORT
  233.    int col=mrc.rowcol; mrc.skip=col; mrc.cw+=IsACopy;
  234.    int i=nrows-col; mrc.storage=i; real* ColCopy;
  235.    if ((mrc.cw*StoreHere)==0)
  236.    {
  237.       REPORT                                            // not accessed
  238.       ColCopy = new real [i]; MatrixErrorNoSpace(ColCopy);
  239.       MONITOR("Make (LT GetCol)",i,ColCopy) 
  240.       mrc.store = ColCopy-col;
  241.    }
  242.    else { REPORT ColCopy = mrc.store+col; }
  243.    if (mrc.cw*LoadOnEntry)
  244.    {
  245.       REPORT
  246.       real* Mstore = store+(col*(col+3))/2;
  247.       while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
  248.    }
  249. }
  250.  
  251. void LowerTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
  252. {
  253. //  if (mrc.cw*StoreOnExit)
  254.    {
  255.       REPORT
  256.       int col=mrc.rowcol; real* Cstore = mrc.store+col;
  257.       real* Mstore = store+(col*(col+3))/2; int i=nrows-col;
  258.       while (i--) { *Mstore = *Cstore++; Mstore += ++col; }
  259.    }
  260. }
  261.  
  262. void LowerTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrLT(); }
  263.                                      //712
  264. // routines for symmetric matrix
  265.  
  266. void SymmetricMatrix::GetRow(MatrixRowCol& mrc)
  267. {
  268.    REPORT                                                //571
  269.    mrc.skip=0; int row=mrc.rowcol;
  270.    if (mrc.cw*DirectPart)
  271.    {
  272.       REPORT
  273.       mrc.cw-=IsACopy; mrc.storage=row+1; mrc.store=store+(row*(row+1))/2;
  274.    }
  275.    else
  276.    {
  277.       mrc.cw+=IsACopy; mrc.storage=ncols;
  278.       real* RowCopy = new real [ncols]; MatrixErrorNoSpace(RowCopy);
  279.       MONITOR("Make (SymGetRow)",ncols,RowCopy) 
  280.       mrc.store = RowCopy;
  281.       if (mrc.cw*LoadOnEntry)
  282.       {
  283.      REPORT                                         // 544
  284.          real* Mstore = store+(row*(row+1))/2; int i = row;
  285.          while (i--) *RowCopy++ = *Mstore++;
  286.          i = ncols-row;
  287.      while (i--) { *RowCopy++ = *Mstore; Mstore += ++row; }
  288.       }
  289.    }
  290. }
  291.  
  292. // need to check this out under StoreHere
  293.  
  294. void SymmetricMatrix::GetCol(MatrixRowCol& mrc)
  295. {
  296.    REPORT
  297.    mrc.skip=0; int col=mrc.rowcol;
  298.    if (mrc.cw*DirectPart)
  299.    {
  300.       REPORT                                         // not accessed
  301.       mrc.cw-=IsACopy; mrc.storage=col+1; mrc.store=store+(col*(col+1))/2;
  302.    }
  303.    else
  304.    {
  305.       mrc.cw+=IsACopy; mrc.storage=ncols; real* ColCopy;
  306.       if ((mrc.cw*StoreHere)==0)
  307.       {
  308.          REPORT                                      // not accessed
  309.          ColCopy = new real [ncols]; MatrixErrorNoSpace(ColCopy);
  310.          MONITOR("Make (SymGetCol)",ncols,ColCopy) 
  311.          mrc.store = ColCopy;
  312.       }
  313.       else { REPORT ColCopy = mrc.store; }
  314.       if (mrc.cw*LoadOnEntry)
  315.       {
  316.          REPORT
  317.          real* Mstore = store+(col*(col+1))/2; int i = col;
  318.          while (i--) *ColCopy++ = *Mstore++;
  319.          i = ncols-col;
  320.      while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
  321.       }
  322.    }
  323. }
  324.  
  325. //void SymmetricMatrix::RestoreRow(int row, real* Rstore)
  326. //{
  327. ////   if (cw*IsACopy && cw*StoreOnExit)
  328. //   {
  329. //      real* Mstore = store+(row*(row+1))/2; int i = row+1;
  330. //      while (i--) *Mstore++ = *Rstore++;
  331. //   }
  332. //}
  333.  
  334. //void SymmetricMatrix::RestoreCol(int col, real* Cstore)
  335. //{
  336. ////   if (cw*IsACopy && cw*StoreOnExit)
  337. //   {
  338. //      real* Mstore = store+(col*(col+3))/2;
  339. //      int i = nrows-col; int j = col;
  340. //      while (i--) { *Mstore = *Cstore++; Mstore+= ++j; }
  341. //   }
  342. //}
  343.  
  344. // routines for row vector
  345.  
  346. void RowVector::GetCol(MatrixRowCol& mrc)
  347. {
  348.    REPORT 
  349.    mrc.skip=0; mrc.storage=1;
  350.    if (mrc.cw >= StoreHere)
  351.    {
  352.       if (mrc.cw >= LoadOnEntry) { REPORT *(mrc.store) = *(store+mrc.rowcol); }
  353.       mrc.cw+=IsACopy;
  354.    }
  355.    else  { REPORT mrc.store = store+mrc.rowcol; mrc.cw-=IsACopy; }
  356.                                                          // not accessed
  357. }
  358.  
  359. void RowVector::NextCol(MatrixRowCol& mrc) 
  360. {
  361.    REPORT
  362.    if (mrc.cw*StoreHere)
  363.    {
  364.       if (mrc.cw*StoreOnExit) { REPORT *(store+mrc.rowcol)=*(mrc.store); }
  365.                                                          // not accessed
  366.       mrc.rowcol++;
  367.       if (mrc.rowcol < ncols)
  368.       {
  369.          if (mrc.cw*LoadOnEntry) { REPORT *(mrc.store)=*(store+mrc.rowcol); }
  370.       }
  371.       else { REPORT mrc.cw -= StoreOnExit; }
  372.    }
  373.    else  { REPORT mrc.rowcol++; mrc.store++; }             // not accessed
  374. }
  375.  
  376. void RowVector::RestoreCol(MatrixRowCol& mrc)
  377. {
  378.    REPORT                                            // not accessed
  379.    if (mrc.cw>=IsACopy)  { REPORT *(store+mrc.rowcol)=*(mrc.store); }
  380. }
  381.