newmat4.cxx

科学和工程计算中使用统计功能开发工具包.

int ConstMatrix::NrowsV() const { return cgm->Nrows(); }
int ReturnMatrixX::NrowsV() const { return gm->Nrows(); }

int GeneralMatrix::NcolsV() const { return ncols; }
int ColedMatrix::NcolsV() const { return 1; }
int MatedMatrix::NcolsV() const { return nc; }
int GetSubMatrix::NcolsV() const { return col_number; }
int MultipliedMatrix::NcolsV() const  { return bm2->NcolsV(); }
int ShiftedMatrix::NcolsV() const { return bm->NcolsV(); }
int TransposedMatrix::NcolsV() const { return bm->NrowsV(); }
int NegatedMatrix::NcolsV() const { return bm->NcolsV(); }
int RowedMatrix::NcolsV() const { return bm->NrowsV() * bm->NcolsV(); }
int DiagedMatrix::NcolsV() const { return bm->NrowsV() * bm->NcolsV(); }
int ConstMatrix::NcolsV() const { return cgm->Ncols(); }
int ReturnMatrixX::NcolsV() const { return gm->Ncols(); }
*/

MatrixBandWidth BaseMatrix::BandWidth() const { return -1; }
MatrixBandWidth DiagonalMatrix::BandWidth() const { return 0; }

MatrixBandWidth BandMatrix::BandWidth() const
   { return MatrixBandWidth(lower,upper); }

MatrixBandWidth AddedMatrix::BandWidth() const
{ return gm1->BandWidth() + gm2->BandWidth(); }

MatrixBandWidth MultipliedMatrix::BandWidth() const
{ return gm1->BandWidth() * gm2->BandWidth(); }

MatrixBandWidth SolvedMatrix::BandWidth() const { return -1; }
MatrixBandWidth ScaledMatrix::BandWidth() const { return gm->BandWidth(); }
MatrixBandWidth NegatedMatrix::BandWidth() const { return gm->BandWidth(); }

MatrixBandWidth TransposedMatrix::BandWidth() const
{ return gm->BandWidth().t(); }

MatrixBandWidth InvertedMatrix::BandWidth() const { return -1; }
MatrixBandWidth RowedMatrix::BandWidth() const { return -1; }
MatrixBandWidth ColedMatrix::BandWidth() const { return -1; }
MatrixBandWidth DiagedMatrix::BandWidth() const { return 0; }
MatrixBandWidth MatedMatrix::BandWidth() const { return -1; }
MatrixBandWidth ConstMatrix::BandWidth() const { return cgm->BandWidth(); }
MatrixBandWidth ReturnMatrixX::BandWidth() const { return gm->BandWidth(); }

MatrixBandWidth GetSubMatrix::BandWidth() const
{

   if (row_skip==col_skip && row_number==col_number) return gm->BandWidth();
   else return MatrixBandWidth(-1);
}



//  Rules regarding tDelete, reuse, GetStore
//    All matrices processed during expression evaluation must be subject
//    to exactly one of reuse(), tDelete(), GetStore() or BorrowStore().
//    If reuse returns TRUE the matrix must be reused.
//    GetMatrix(gm) always calls gm->GetStore()
//    gm->Evaluate obeys rules
//    bm->Evaluate obeys rules for matrices in bm structure

void GeneralMatrix::tDelete()
{
   if (tag<0)
   {
      if (tag<-1) { REPORT store=0; delete this; return; }  // borrowed
      else { REPORT return; }
   }
   if (tag==1)
   {
      REPORT  MONITOR_REAL_DELETE("Free   (tDelete)",storage,store)
#ifdef Version21
      if (store) delete [] store;
#else
      if (store) delete [storage] store;
#endif
      store=0; tag=-1; return;
   }
   if (tag==0) { REPORT delete this; return; }
   REPORT tag--; return;
}

static void BlockCopy(int n, Real* from, Real* to)
{
   REPORT
   int i = (n >> 3);
   while (i--)
   {
      *to++ = *from++; *to++ = *from++; *to++ = *from++; *to++ = *from++;
      *to++ = *from++; *to++ = *from++; *to++ = *from++; *to++ = *from++;
   }
   i = n & 7; while (i--) *to++ = *from++;
}

Boolean GeneralMatrix::reuse()
{
   if (tag<-1)
   {
      REPORT
      Real* s = new Real [storage]; MatrixErrorNoSpace(s);
      MONITOR_REAL_NEW("Make     (reuse)",storage,s)
      BlockCopy(storage, store, s); store=s; tag=0; return TRUE;
   }
   if (tag<0) { REPORT return FALSE; }
   if (tag<=1)  { REPORT return TRUE; }
   REPORT tag--; return FALSE;
}

Real* GeneralMatrix::GetStore()
{
   if (tag<0 || tag>1)
   {
      Real* s = new Real [storage]; MatrixErrorNoSpace(s);
      MONITOR_REAL_NEW("Make  (GetStore)",storage,s)
      BlockCopy(storage, store, s);
      if (tag>1) { REPORT tag--; }
      else if (tag < -1) { REPORT store=0; delete this; } // borrowed store
      else { REPORT }
      return s;
   }
   Real* s=store; store=0;
   if (tag==0) { REPORT delete this; }
   else { REPORT tag=-1; }
   return s;
}

/*
#ifndef __ZTC__
void GeneralMatrix::GetMatrixC(const GeneralMatrix* gmx)
{
   REPORT tag=-1;
   nrows=gmx->nrows; ncols=gmx->ncols; storage=gmx->storage;
   SetParameters(gmx); 
   store = new Real [storage]; MatrixErrorNoSpace(store);
   MONITOR_REAL_NEW("Make (GetMatrix)",storage,store)
   BlockCopy(storage, gmx->store, store);
}
#endif
*/

void GeneralMatrix::GetMatrix(const GeneralMatrix* gmx)
{
   REPORT  tag=-1; nrows=gmx->Nrows(); ncols=gmx->Ncols();
   storage=gmx->storage; SetParameters(gmx);
   store=((GeneralMatrix*)gmx)->GetStore();
}

GeneralMatrix* GeneralMatrix::BorrowStore(GeneralMatrix* gmx, MatrixType mt)
// Copy storage of *this to storage of *gmx. Then convert to type mt.
// If mt == 0 just let *gm point to storage of *this if tag<0.
{
   if (!mt)
   {
      if (tag == -1) { REPORT gmx->tag = -2; gmx->store = store; }
      else { REPORT gmx->tag = 0; gmx->store = GetStore(); }
   }
   else if (mt!=gmx->Type())
   {
      REPORT gmx->tag = -2; gmx->store = store;
      gmx = gmx->Evaluate(mt); gmx->tag = 0; tDelete();
   }
   else { REPORT  gmx->tag = 0; gmx->store = GetStore(); }

   return gmx;
}

void GeneralMatrix::Eq(const BaseMatrix& X, MatrixType mt)
// Count number of references to this in X.
// If zero delete storage in X;
// otherwise tag X to show when storage can be deleted
// evaluate X and copy to current object
{
   int counter=X.search(this);
   if (counter==0)
   {
      REPORT
      if (store)
      {
         MONITOR_REAL_DELETE("Free (operator=)",storage,store)
#ifdef Version21
         REPORT delete [] store; storage=0;
#else
         REPORT delete [storage] store; storage=0;
#endif
      }
   }
   else { REPORT Release(counter); }
   GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate(mt);
   if (gmx!=this) { REPORT GetMatrix(gmx); }
   else { REPORT }
   Protect();
}

void GeneralMatrix::Inject(const GeneralMatrix& X)
// copy stored values of X; otherwise leave els of *this unchanged
{
   REPORT
   Tracer tr("Inject");
   if (nrows != X.nrows || ncols != X.ncols)
      Throw(IncompatibleDimensionsException());
   MatrixRow mr((GeneralMatrix*)&X, LoadOnEntry);
   MatrixRow mrx(this, LoadOnEntry+StoreOnExit+DirectPart);
   int i=nrows;
   while (i--) { mrx.Inject(mr); mrx.Next(); mr.Next(); }
}  

void GeneralMatrix::CheckConversion(const BaseMatrix& M)
{
   if (!(this->Type() >= M.Type()))
      Throw(ProgramException("Illegal Conversion"));
}


/************************* nricMatrix routines *****************************/

void nricMatrix::MakeRowPointer()
{
   row_pointer = new Real* [nrows]; MatrixErrorNoSpace(row_pointer);
   Real* s = Store() - 1; int i = nrows; Real** rp = row_pointer;
   while (i--) { *rp++ = s; s+=ncols; }
}

void nricMatrix::DeleteRowPointer()
#ifdef Version21
{ if (nrows) delete [] row_pointer; }
#else
{ if (nrows) delete [nrows] row_pointer; }
#endif

void GeneralMatrix::CheckStore() const
{
   if (!store) 
      Throw(ProgramException("NRIC accessing matrix with unset dimensions"));
}


/***************************** CleanUp routines *****************************/

void GeneralMatrix::CleanUp()
{
   // set matrix dimensions to zero, delete storage
   REPORT
   if (store && storage)
   {
      MONITOR_REAL_DELETE("Free (CleanUp)    ",storage,store)
#ifdef Version21
         REPORT delete [] store;
#else
         REPORT delete [storage] store;
#endif
   }
   store=0; storage=0; nrows=0; ncols=0;
}

void nricMatrix::CleanUp()
{ DeleteRowPointer(); GeneralMatrix::CleanUp(); }

void RowVector::CleanUp()
{ GeneralMatrix::CleanUp(); nrows=1; }

void ColumnVector::CleanUp()
{ GeneralMatrix::CleanUp(); ncols=1; }

void CroutMatrix::CleanUp()
{ 
#ifdef Version21
   if (nrows) delete [] indx;
#else
   if (nrows) delete [nrows] indx;
#endif
   GeneralMatrix::CleanUp();
}

void BandLUMatrix::CleanUp()
{ 
#ifdef Version21
   if (nrows) delete [] indx;
   if (storage2) delete [] store2;
#else
   if (nrows) delete [nrows] indx;
   if (storage2) delete [storage2] store2;
#endif
   GeneralMatrix::CleanUp();
}