newmat3.cpp

矩阵计算库

   REPORT                                                //571
   mrc.skip=0; int row=mrc.rowcol; mrc.length=ncols;
   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_REAL_NEW("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)
{
	int col=mrc.rowcol; mrc.length=nrows;
	if (mrc.cw >= StoreHere+DirectPart)
	{
		REPORT
		mrc.skip=col; mrc.cw+=IsACopy;
		int i=nrows-col; mrc.storage=i;
		Real* ColCopy = mrc.store+col;
		if (+(mrc.cw*LoadOnEntry))
		{
			REPORT                           // not accessed
			Real* Mstore = store+(col*(col+3))/2;
			while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
		}
	}
	else
	{
		REPORT
		mrc.skip=0;
		if (+(mrc.cw*DirectPart))
		{
			REPORT
			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) )
			{
				REPORT                                      // not accessed
				ColCopy = new Real [ncols]; MatrixErrorNoSpace(ColCopy);
				MONITOR_REAL_NEW("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(MatrixRowCol& mrc)
{
if (+(mrc.cw*StoreHere))
	{
		REPORT
		// Really do restore column
		// Restoring row instead would be possible but would make
		// solve less efficient
		// Same as LowerTriangular::RestoreCol

		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; }
	}
}

// routines for row vector

void RowVector::GetCol(MatrixRowCol& mrc)
{
	REPORT
	mrc.skip=0; mrc.storage=1; mrc.length=nrows;
	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); }
}


// routines for band matrices

void BandMatrix::GetRow(MatrixRowCol& mrc)
{
   REPORT
   mrc.cw -= IsACopy; int r = mrc.rowcol; int w = lower+1+upper;
   mrc.length=ncols;
   int s = r-lower; mrc.store = store+(r*w-s); if (s<0) { w += s; s = 0; }
   mrc.skip = s; s += w-ncols; if (s>0) w -= s; mrc.storage = w;
}

// make special versions of this for upper and lower band matrices

void BandMatrix::NextRow(MatrixRowCol& mrc)
{
   REPORT
   int r = ++mrc.rowcol; mrc.store += lower+upper;
   if (r<=lower) mrc.storage++; else mrc.skip++;
   if (r>=ncols-upper) mrc.storage--;
}

void BandMatrix::GetCol(MatrixRowCol& mrc)
{
   REPORT
   mrc.cw += IsACopy; int c = mrc.rowcol; int n = lower+upper; int w = n+1;
   mrc.length=nrows;
   int b; int s = c-upper; Real* ColCopy;
   if (s<=0) { w += s; s = 0; b = c+lower; } else b = s*w+n;
   mrc.skip = s; s += w-nrows; if (s>0) w -= s; mrc.storage = w;
   if ( !(mrc.cw*StoreHere) )
   {
      REPORT
      ColCopy = new Real [w]; MatrixErrorNoSpace(ColCopy);
      MONITOR_REAL_NEW("Make (BMGetCol)",w,ColCopy)
      mrc.store = ColCopy-mrc.skip;
   }
   else { REPORT ColCopy = mrc.store+mrc.skip; }
   if (+(mrc.cw*LoadOnEntry))
   {
      REPORT
      Real* Mstore = store+b;
      while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
   }
}

void BandMatrix::RestoreCol(MatrixRowCol& mrc)
{
//  if (mrc.cw*StoreOnExit)
	REPORT
	int c = mrc.rowcol; int n = lower+upper; int s = c-upper;
	Real* Mstore = store + ((s<=0) ? c+lower : s*n+s+n);
	Real* Cstore = mrc.store+mrc.skip; int w = mrc.storage;
	while (w--) { *Mstore = *Cstore++; Mstore += n; }
}

// routines for symmetric band matrix

void SymmetricBandMatrix::GetRow(MatrixRowCol& mrc)
{
   REPORT
   int r=mrc.rowcol; int s = r-lower; int w1 = lower+1; int o = r*w1;
   mrc.length = ncols;
   if (s<0) { w1 += s; o -= s; s = 0; }  mrc.skip = s;

   if (+(mrc.cw*DirectPart))
      { REPORT  mrc.cw -= IsACopy; mrc.store = store+o-s; mrc.storage = w1; }
   else
   {
      mrc.cw += IsACopy; int w = w1+lower; s += w-ncols;
      if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
      Real* RowCopy = new Real [w]; MatrixErrorNoSpace(RowCopy);
      MONITOR_REAL_NEW("Make (SmBGetRow)",w,RowCopy) 
      mrc.store = RowCopy-mrc.skip;
      if (+(mrc.cw*LoadOnEntry))
      {
	 REPORT
         Real* Mstore = store+o;
         while (w1--) *RowCopy++ = *Mstore++;   Mstore--;
         while (w2--) { Mstore += lower; *RowCopy++ = *Mstore; }
      }
   }
}

// need to check this out under StoreHere

void SymmetricBandMatrix::GetCol(MatrixRowCol& mrc)
{
   int c=mrc.rowcol; int w1 = lower+1; mrc.length=nrows;
   if (mrc.cw >= StoreHere+DirectPart)
   {
      REPORT
      mrc.cw += IsACopy;
      Real* ColCopy; int b = c*w1+lower;
      mrc.skip = c; c += w1-nrows; w1 -= c; mrc.storage = w1;
      REPORT ColCopy = mrc.store+mrc.skip;
      if (+(mrc.cw*LoadOnEntry))
      {
         REPORT
         Real* Mstore = store+b;
         while (w1--) { *ColCopy++ = *Mstore; Mstore += lower; }
      }
   }
   else
   {
      REPORT
      int s = c-lower; int o = c*w1;
      if (s<0) { w1 += s; o -= s; s = 0; }  mrc.skip = s;

      if (+(mrc.cw*DirectPart))
         { REPORT  mrc.cw -= IsACopy; mrc.store = store+o-s; mrc.storage = w1; }
      else
      {
         mrc.cw += IsACopy; int w = w1+lower; s += w-ncols;
         if (s>0) w -= s; mrc.storage = w; int w2 = w-w1; Real* ColCopy;
         if ( !(mrc.cw*StoreHere) )
         {
            REPORT ColCopy = new Real [w]; MatrixErrorNoSpace(ColCopy);
            MONITOR_REAL_NEW("Make (SmBGetCol)",w,ColCopy)
            mrc.store = ColCopy-mrc.skip;
         }
         else { REPORT ColCopy = mrc.store+mrc.skip; }
         if (+(mrc.cw*LoadOnEntry))
         {
            REPORT
            Real* Mstore = store+o;
            while (w1--) *ColCopy++ = *Mstore++;   Mstore--;
            while (w2--) { Mstore += lower; *ColCopy++ = *Mstore; }
         }
      }
   }
}

void SymmetricBandMatrix::RestoreCol(MatrixRowCol& mrc)
{
if (+(mrc.cw*StoreHere))
   {
      REPORT
      int c = mrc.rowcol;
      Real* Mstore = store + c*lower+c+lower;
      Real* Cstore = mrc.store+mrc.skip; int w = mrc.storage;
      while (w--) { *Mstore = *Cstore++; Mstore += lower; }
   }
}