//$$ newmat3.cxx        Matrix get and restore rows and columns

// Copyright (C) 1991: R B Davies and DSIR


#include "include.hxx"

#include "newmat.hxx"
#include "newmatrc.hxx"

//#define REPORT { static ExeCounter ExeCount(__LINE__,3); ExeCount++; }

#define REPORT {}

//#define MONITOR(what,storage,store) \
//   { cout << what << " " << storage << " at " << (long)store << "\n"; }

#define MONITOR(what,store,storage) {}


// Control bits codes for GetRow, GetCol, RestoreRow, RestoreCol
//
// LoadOnEntry:
//    Load data into MatrixRow or Col dummy array under GetRow or GetCol
// StoreOnExit:
//    Restore data to original matrix under RestoreRow or RestoreCol
// IsACopy:
//    Set by GetRow/Col: MatrixRow or Col array is a copy
// DirectPart:
//    Load or restore only part directly stored; must be set with StoreOnExit
//    Still have decide  how to handle this with symmetric
// StoreHere:
//    used in columns only - store data at supplied storage address, adjusted
//    for skip; used for GetCol, NextCol & RestoreCol. No need to fill out
//    zeros.


// These will work as a default
// but need to override NextRow for efficiency


void GeneralMatrix::NextRow(MatrixRowCol& mrc)
{
   REPORT
   if (mrc.cw*StoreOnExit) { REPORT this->RestoreRow(mrc); }
   if (mrc.cw*IsACopy)
   {
      REPORT
      real* s = mrc.store + mrc.skip;
      MONITOR("Free   (NextRow)",mrc.storage,s) 
      delete [mrc.storage] s;
   }
   mrc.rowcol++;
   if (mrc.rowcol<nrows) { REPORT this->GetRow(mrc); }
   else { REPORT mrc.cw -= (StoreOnExit+IsACopy); }
}

void GeneralMatrix::NextCol(MatrixRowCol& mrc)
{
   REPORT                                                // 423
   if (mrc.cw*StoreOnExit) { REPORT this->RestoreCol(mrc); }
   if (mrc.cw*IsACopy && (mrc.cw*StoreHere==0))
   {
      REPORT                                             // not accessed
      real* s = mrc.store + mrc.skip;
      MONITOR("Free   (NextCol)",mrc.storage,s) 
      delete [mrc.storage] s;
   }
   mrc.rowcol++;
   if (mrc.rowcol<ncols) { REPORT this->GetCol(mrc); }
   else { REPORT mrc.cw -= (StoreOnExit+IsACopy); }
}


// routines for matrix

void Matrix::GetRow(MatrixRowCol& mrc)
{
   REPORT
   mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=ncols;
   mrc.store=store+mrc.rowcol*ncols;
}


void Matrix::GetCol(MatrixRowCol& mrc)
{
   REPORT
   mrc.skip=0; mrc.storage=nrows;
   if (ncols==1 && (mrc.cw*StoreHere)==0)
      { REPORT mrc.cw-=IsACopy; mrc.store=store; }           // not accessed
   else
   {
      mrc.cw+=IsACopy; real* ColCopy;
      if ((mrc.cw*StoreHere)==0)
      {
         REPORT
         ColCopy = new real [nrows]; MatrixErrorNoSpace(ColCopy);
         MONITOR("Make (MatGetCol)",nrows,ColCopy)
         mrc.store = ColCopy;
      }
      else { REPORT ColCopy = mrc.store; }
      if (mrc.cw*LoadOnEntry)
      {
         REPORT
         real* Mstore = store+mrc.rowcol; int i=nrows;
         while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
      }
   }
}

void Matrix::RestoreCol(MatrixRowCol& mrc)
{
//  if (mrc.cw*StoreOnExit)
   REPORT                                   // 429
   if (mrc.cw*IsACopy)
   {
      REPORT                                // 426
      real* Mstore = store+mrc.rowcol; int i=nrows; real* Cstore = mrc.store;
      while (i--) { *Mstore = *Cstore++; Mstore+=ncols; }
   }
}

void Matrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrMat(); }  // 1808

void Matrix::NextCol(MatrixRowCol& mrc)
{
   REPORT                                        // 632
   if (mrc.cw*StoreOnExit) { REPORT RestoreCol(mrc); }
   mrc.rowcol++;
   if (mrc.rowcol<ncols)
   {
      if (mrc.cw*LoadOnEntry)
      {
         REPORT
         real* ColCopy = mrc.store;
         real* Mstore = store+mrc.rowcol; int i=nrows;
         while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
      }
   }
   else { REPORT mrc.cw -= StoreOnExit; }
}

// routines for diagonal matrix

void DiagonalMatrix::GetRow(MatrixRowCol& mrc)
{
   REPORT
   mrc.skip=mrc.rowcol; mrc.cw-=IsACopy; mrc.storage=1; mrc.store=store;
}

void DiagonalMatrix::GetCol(MatrixRowCol& mrc)
{
   REPORT 
   mrc.skip=mrc.rowcol; mrc.storage=1;
   if (mrc.cw*StoreHere) 
      { REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); mrc.cw+=IsACopy; }
   else { REPORT mrc.store = store; mrc.cw-=IsACopy; }     // not accessed
}

void DiagonalMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
						      // 800
void DiagonalMatrix::NextCol(MatrixRowCol& mrc)
{
   REPORT
   if (mrc.cw*StoreHere)
   {
      if (mrc.cw*StoreOnExit) 
         { REPORT *(store+mrc.rowcol)=*(mrc.store+mrc.rowcol); }
      mrc.IncrDiag();
      if (mrc.cw*LoadOnEntry && mrc.rowcol < ncols)
         { REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); }
   }
   else { REPORT mrc.IncrDiag(); }                     // not accessed
}

// routines for upper triangular matrix

void UpperTriangularMatrix::GetRow(MatrixRowCol& mrc)
{
   REPORT
   int row = mrc.rowcol; mrc.skip=row; mrc.cw-=IsACopy;
   mrc.storage=ncols-row; mrc.store=store+(row*(2*ncols-row-1))/2;
}


void UpperTriangularMatrix::GetCol(MatrixRowCol& mrc)
{
   REPORT
   mrc.skip=0; mrc.cw+=IsACopy; int i=mrc.rowcol+1; mrc.storage=i;
   real* ColCopy;
   if ((mrc.cw*StoreHere)==0)
   {
      REPORT                                              // not accessed
      ColCopy = new real [i]; MatrixErrorNoSpace(ColCopy);
      MONITOR("Make (UT GetCol)",i,ColCopy) 
      mrc.store = ColCopy;
   }
   else { REPORT ColCopy = mrc.store; }
   if (mrc.cw*LoadOnEntry)
   {
      REPORT
      real* Mstore = store+mrc.rowcol; int j = ncols;
      while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
   }
}

void UpperTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
{
//  if (mrc.cw*StoreOnExit)
  {
     REPORT
     real* Mstore = store+mrc.rowcol; int i=mrc.rowcol+1; int j = ncols;
     real* Cstore = mrc.store;
     while (i--) { *Mstore = *Cstore++; Mstore += --j; }
  }
}

void UpperTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrUT(); }
						      // 722

// routines for lower triangular matrix

void LowerTriangularMatrix::GetRow(MatrixRowCol& mrc)
{
   REPORT
   int row=mrc.rowcol; mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=row+1;
   mrc.store=store+(row*(row+1))/2;
}

void LowerTriangularMatrix::GetCol(MatrixRowCol& mrc)
{
   REPORT
   int col=mrc.rowcol; mrc.skip=col; mrc.cw+=IsACopy;
   int i=nrows-col; mrc.storage=i; real* ColCopy;
   if ((mrc.cw*StoreHere)==0)
   {
      REPORT                                            // not accessed
      ColCopy = new real [i]; MatrixErrorNoSpace(ColCopy);
      MONITOR("Make (LT GetCol)",i,ColCopy) 
      mrc.store = ColCopy-col;
   }
   else { REPORT ColCopy = mrc.store+col; }
   if (mrc.cw*LoadOnEntry)
   {
      REPORT
      real* Mstore = store+(col*(col+3))/2;
      while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
   }
}

void LowerTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
{
//  if (mrc.cw*StoreOnExit)
   {
      REPORT
      int col=mrc.rowcol; real* Cstore = mrc.store+col;
      real* Mstore = store+(col*(col+3))/2; int i=nrows-col;
      while (i--) { *Mstore = *Cstore++; Mstore += ++col; }
   }
}

void LowerTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrLT(); }
					                 //712
// routines for symmetric matrix

void SymmetricMatrix::GetRow(MatrixRowCol& mrc)
{
   REPORT                                                //571
   mrc.skip=0; int row=mrc.rowcol;
   if (mrc.cw*DirectPart)
   {
      REPORT
      mrc.cw-=IsACopy; mrc.storage=row+1; mrc.store=store+(row*(row+1))/2;
   }
   else
   {
      mrc.cw+=IsACopy; mrc.storage=ncols;
      real* RowCopy = new real [ncols]; MatrixErrorNoSpace(RowCopy);
      MONITOR("Make (SymGetRow)",ncols,RowCopy) 
      mrc.store = RowCopy;
      if (mrc.cw*LoadOnEntry)
      {
	 REPORT                                         // 544
         real* Mstore = store+(row*(row+1))/2; int i = row;
         while (i--) *RowCopy++ = *Mstore++;
         i = ncols-row;
	 while (i--) { *RowCopy++ = *Mstore; Mstore += ++row; }
      }
   }
}

// need to check this out under StoreHere

void SymmetricMatrix::GetCol(MatrixRowCol& mrc)
{
   REPORT
   mrc.skip=0; int col=mrc.rowcol;
   if (mrc.cw*DirectPart)
   {
      REPORT                                         // not accessed
      mrc.cw-=IsACopy; mrc.storage=col+1; mrc.store=store+(col*(col+1))/2;
   }
   else
   {
      mrc.cw+=IsACopy; mrc.storage=ncols; real* ColCopy;
      if ((mrc.cw*StoreHere)==0)
      {
         REPORT                                      // not accessed
         ColCopy = new real [ncols]; MatrixErrorNoSpace(ColCopy);
         MONITOR("Make (SymGetCol)",ncols,ColCopy) 
         mrc.store = ColCopy;
      }
      else { REPORT ColCopy = mrc.store; }
      if (mrc.cw*LoadOnEntry)
      {
         REPORT
         real* Mstore = store+(col*(col+1))/2; int i = col;
         while (i--) *ColCopy++ = *Mstore++;
         i = ncols-col;
	 while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
      }
   }
}

//void SymmetricMatrix::RestoreRow(int row, real* Rstore)
//{
////   if (cw*IsACopy && cw*StoreOnExit)
//   {
//      real* Mstore = store+(row*(row+1))/2; int i = row+1;
//      while (i--) *Mstore++ = *Rstore++;
//   }
//}

//void SymmetricMatrix::RestoreCol(int col, real* Cstore)
//{
////   if (cw*IsACopy && cw*StoreOnExit)
//   {
//      real* Mstore = store+(col*(col+3))/2;
//      int i = nrows-col; int j = col;
//      while (i--) { *Mstore = *Cstore++; Mstore+= ++j; }
//   }
//}

// routines for row vector

void RowVector::GetCol(MatrixRowCol& mrc)
{
   REPORT 
   mrc.skip=0; mrc.storage=1;
   if (mrc.cw >= StoreHere)
   {
      if (mrc.cw >= LoadOnEntry) { REPORT *(mrc.store) = *(store+mrc.rowcol); }
      mrc.cw+=IsACopy;
   }
   else  { REPORT mrc.store = store+mrc.rowcol; mrc.cw-=IsACopy; }
                                                         // not accessed
}

void RowVector::NextCol(MatrixRowCol& mrc) 
{
   REPORT
   if (mrc.cw*StoreHere)
   {
      if (mrc.cw*StoreOnExit) { REPORT *(store+mrc.rowcol)=*(mrc.store); }
                                                         // not accessed
      mrc.rowcol++;
      if (mrc.rowcol < ncols)
      {
         if (mrc.cw*LoadOnEntry) { REPORT *(mrc.store)=*(store+mrc.rowcol); }
      }
      else { REPORT mrc.cw -= StoreOnExit; }
   }
   else  { REPORT mrc.rowcol++; mrc.store++; }             // not accessed
}

void RowVector::RestoreCol(MatrixRowCol& mrc)
{
   REPORT                                            // not accessed
   if (mrc.cw>=IsACopy)  { REPORT *(store+mrc.rowcol)=*(mrc.store); }
}
