newmat3.cpp

矩阵计算库

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

// Copyright (C) 1991,2,3,4: R B Davies

#include "include.h"

#include "newmat.h"
#include "newmatrc.h"

//#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.

// For symmetric matrices, treat columns as rows unless StoreHere is set;
// then stick to columns as this will give better performance for doing
// inverses


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

// Assume pointer arithmetic works for pointers out of range - not strict C++.


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_REAL_DELETE("Free   (NextRow)",mrc.storage,s)
#ifdef Version21
		delete [] s;
#else
		delete [mrc.storage] s;
#endif
	}
	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); }
	int t1 = +(mrc.cw*IsACopy); int t2 = !(mrc.cw*StoreHere);
	if ( t1 && t2 )
	{
		REPORT                                             // not accessed
		Real* s = mrc.store + mrc.skip;
		MONITOR_REAL_DELETE("Free   (NextCol)",mrc.storage,s)
#ifdef Version21
		delete [] s;
#else
		delete [mrc.storage] s;
#endif
	}
	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.length=ncols;
	mrc.store=store+mrc.rowcol*ncols;
}


void Matrix::GetCol(MatrixRowCol& mrc)
{
	REPORT
	mrc.skip=0; mrc.storage=nrows; mrc.length=nrows;
        int t1 = !(mrc.cw*StoreHere);
	if ( ncols==1 && t1 )
		{ REPORT mrc.cw-=IsACopy; mrc.store=store; }
	else
	{
		mrc.cw+=IsACopy; Real* ColCopy;
		if ( !(mrc.cw*StoreHere) )
		{
			REPORT
			ColCopy = new Real [nrows]; MatrixErrorNoSpace(ColCopy);
			MONITOR_REAL_NEW("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;
   mrc.length=ncols;
}

void DiagonalMatrix::GetCol(MatrixRowCol& mrc)
{
   REPORT 
   mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows;
   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();
      int t1 = +(mrc.cw*LoadOnEntry);
      if (t1 && 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.length=ncols;
   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;
   mrc.length=nrows;
   Real* ColCopy;
   if ( !(mrc.cw*StoreHere) )
   {
      REPORT                                              // not accessed
      ColCopy = new Real [i]; MatrixErrorNoSpace(ColCopy);
      MONITOR_REAL_NEW("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.length=ncols; mrc.store=store+(row*(row+1))/2;
}

void LowerTriangularMatrix::GetCol(MatrixRowCol& mrc)
{
   REPORT
   int col=mrc.rowcol; mrc.skip=col; mrc.cw+=IsACopy; mrc.length=nrows; 
   int i=nrows-col; mrc.storage=i; Real* ColCopy;
   if ( !(mrc.cw*StoreHere) )
   {
      REPORT                                            // not accessed
      ColCopy = new Real [i]; MatrixErrorNoSpace(ColCopy);
      MONITOR_REAL_NEW("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)
{