newmat5.cxx

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

//$$ newmat5.cxx         Transpose, evaluate etc

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

#include "include.h"

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

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

#define REPORT {}


/************************ carry out operations ******************************/


GeneralMatrix* GeneralMatrix::Transpose(TransposedMatrix* tm, MatrixType mt)
{
   if (!mt) mt = Type().t();           // type of transposed matrix
   GeneralMatrix* gm1 = mt.New(ncols,nrows,tm); gm1->ReleaseAndDelete();

   if (mt == Type().t())
   {
      REPORT
      for (int i=0; i<ncols; i++)
      {
	 MatrixRow mr(gm1, StoreOnExit+DirectPart, i);
         MatrixCol mc(this, mr.Store(), LoadOnEntry, i);
      }
   }
   else
   {
      REPORT
      MatrixRow mr(this, LoadOnEntry);
      MatrixCol mc(gm1, StoreOnExit+DirectPart);
      int i = nrows;
      while (i--) { mc.Copy(mr); mr.Next(); mc.Next(); }
   }
   tDelete(); return gm1;
}

GeneralMatrix* SymmetricMatrix::Transpose(TransposedMatrix*, MatrixType mt)
{ REPORT  return Evaluate(mt); }


GeneralMatrix* DiagonalMatrix::Transpose(TransposedMatrix*, MatrixType mt)
{ REPORT return Evaluate(mt); }

Boolean GeneralMatrix::IsZero() const
{
   REPORT
   Real* s=store; int i=storage;
   while (i--) { if (*s++) return FALSE; }
   return TRUE;
}

GeneralMatrix* ColumnVector::Transpose(TransposedMatrix*, MatrixType mt)
{
   REPORT
   GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
   gmx->nrows = 1; gmx->ncols = gmx->storage = storage;
   return BorrowStore(gmx,mt);
}

GeneralMatrix* RowVector::Transpose(TransposedMatrix*, MatrixType mt)
{
   REPORT
   GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
   gmx->ncols = 1; gmx->nrows = gmx->storage = storage;
   return BorrowStore(gmx,mt);
}

GeneralMatrix* GeneralMatrix::Evaluate(MatrixType mt)
{
   if (!mt) { REPORT return this; }
   if (mt==this->Type()) { REPORT return this; }
   REPORT
   GeneralMatrix* gmx = mt.New(nrows,ncols,this);
   MatrixRow mr(this, LoadOnEntry);
   MatrixRow mrx(gmx, StoreOnExit+DirectPart);
   int i=nrows;
   while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
   tDelete(); gmx->ReleaseAndDelete(); return gmx;
}

GeneralMatrix* ConstMatrix::Evaluate(MatrixType mt)
{
   if (!mt || mt==cgm->Type())
   {
      REPORT
#ifdef TEMPS_DESTROYED_QUICKLY
      GeneralMatrix* gmx = (GeneralMatrix*)cgm; delete this; return gmx;
#else
      return (GeneralMatrix*)cgm;
#endif
   }
   REPORT
   GeneralMatrix* gmx = cgm->Type().New(cgm->Nrows(),cgm->Ncols(),this);
   MatrixRow mr((GeneralMatrix*)cgm, LoadOnEntry);//assume won't change this
   MatrixRow mrx(gmx, StoreOnExit+DirectPart);
   int i=cgm->Nrows();
   while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
   gmx->ReleaseAndDelete();
#ifdef TEMPS_DESTROYED_QUICKLY
   delete this;
#endif
   return gmx; // no tDelete
}

GeneralMatrix* ShiftedMatrix::Evaluate(MatrixType mt)
{
   gm=((BaseMatrix*&)bm)->Evaluate();
   if (!mt) mt = gm->Type().AddEqualEl();
   int nr=gm->Nrows(); int nc=gm->Ncols();
   if (!(mt==gm->Type()))
   {
      REPORT
      GeneralMatrix* gmx = mt.New(nr,nc,this);
      MatrixRow mr(gm, LoadOnEntry); 
      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
      while (nr--) { mrx.Add(mr,f); mrx.Next(); mr.Next(); }
      gmx->ReleaseAndDelete(); gm->tDelete();
#ifdef TEMPS_DESTROYED_QUICKLY
      delete this;
#endif
      return gmx;
   }
   else if (gm->reuse())
   {
      REPORT gm->Add(f);
#ifdef TEMPS_DESTROYED_QUICKLY
      GeneralMatrix* gmx = gm; delete this; return gmx;
#else
      return gm;
#endif
   }
   else
   {
      REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc,this);
      gmy->ReleaseAndDelete(); gmy->Add(gm,f);
#ifdef TEMPS_DESTROYED_QUICKLY
      delete this;
#endif
      return gmy;
   }
}

GeneralMatrix* ScaledMatrix::Evaluate(MatrixType mt)
{
   gm=((BaseMatrix*&)bm)->Evaluate();
   int nr=gm->Nrows(); int nc=gm->Ncols();
   if (!mt) mt = gm->Type();
   if (mt==gm->Type())
   {
      if (gm->reuse())
      {
         REPORT gm->Multiply(f);
#ifdef TEMPS_DESTROYED_QUICKLY
         GeneralMatrix* gmx = gm; delete this; return gmx;
#else
         return gm;
#endif
      }
      else
      {
         REPORT GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
         gmx->ReleaseAndDelete(); gmx->Multiply(gm,f);
#ifdef TEMPS_DESTROYED_QUICKLY
         delete this;
#endif
         return gmx;
      }
   }
   else
   {
      REPORT
      GeneralMatrix* gmx = mt.New(nr,nc,this);
      MatrixRow mr(gm, LoadOnEntry); 
      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
      while (nr--) { mrx.Multiply(mr,f); mrx.Next(); mr.Next(); }
      gmx->ReleaseAndDelete(); gm->tDelete();
#ifdef TEMPS_DESTROYED_QUICKLY
      delete this;
#endif
      return gmx;
   }
}

GeneralMatrix* NegatedMatrix::Evaluate(MatrixType mt)
{
   gm=((BaseMatrix*&)bm)->Evaluate();
   if (!mt) mt = gm->Type();
   int nr=gm->Nrows(); int nc=gm->Ncols();
   if (mt==gm->Type())
   {
      if (gm->reuse())
      {
         REPORT gm->Negate();
#ifdef TEMPS_DESTROYED_QUICKLY
         GeneralMatrix* gmx = gm; delete this; return gmx;
#else
         return gm;
#endif
      }
      else
      {
         REPORT
         GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
         gmx->ReleaseAndDelete(); gmx->Negate(gm);
#ifdef TEMPS_DESTROYED_QUICKLY
         delete this;
#endif
         return gmx;
      }
   }
   else
   {
      REPORT
      GeneralMatrix* gmx = mt.New(nr,nc,this);
      MatrixRow mr(gm, LoadOnEntry); 
      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
      while (nr--) { mrx.Negate(mr); mrx.Next(); mr.Next(); }
      gmx->ReleaseAndDelete(); gm->tDelete();
#ifdef TEMPS_DESTROYED_QUICKLY
      delete this;
#endif
      return gmx;
   }
}   

GeneralMatrix* TransposedMatrix::Evaluate(MatrixType mt)
{
   REPORT
   gm=((BaseMatrix*&)bm)->Evaluate();
   if (!mt) mt = gm->Type().t();
   GeneralMatrix* gmx=gm->Transpose(this, mt);
#ifdef TEMPS_DESTROYED_QUICKLY
   delete this;
#endif
   return gmx;
}
   
GeneralMatrix* RowedMatrix::Evaluate(MatrixType mt)
{
   gm = ((BaseMatrix*&)bm)->Evaluate();
   GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
   gmx->nrows = 1; gmx->ncols = gmx->storage = gm->storage;
#ifdef TEMPS_DESTROYED_QUICKLY
   GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
#else
   return gm->BorrowStore(gmx,mt);
#endif
}

GeneralMatrix* ColedMatrix::Evaluate(MatrixType mt)
{
   gm = ((BaseMatrix*&)bm)->Evaluate();
   GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
   gmx->ncols = 1; gmx->nrows = gmx->storage = gm->storage;
#ifdef TEMPS_DESTROYED_QUICKLY
   GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
#else
   return gm->BorrowStore(gmx,mt);
#endif
}

GeneralMatrix* DiagedMatrix::Evaluate(MatrixType mt)
{
   gm = ((BaseMatrix*&)bm)->Evaluate();
   GeneralMatrix* gmx = new DiagonalMatrix; MatrixErrorNoSpace(gmx);
   gmx->nrows = gmx->ncols = gmx->storage = gm->storage;
#ifdef TEMPS_DESTROYED_QUICKLY
   GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
#else
   return gm->BorrowStore(gmx,mt);
#endif
}

GeneralMatrix* MatedMatrix::Evaluate(MatrixType mt)
{
   Tracer tr("MatedMatrix::Evaluate");
   gm = ((BaseMatrix*&)bm)->Evaluate();
   GeneralMatrix* gmx = new Matrix; MatrixErrorNoSpace(gmx);
   gmx->nrows = nr; gmx->ncols = nc; gmx->storage = gm->storage;
   if (nr*nc != gmx->storage)
      Throw(IncompatibleDimensionsException());
#ifdef TEMPS_DESTROYED_QUICKLY
   GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
#else
   return gm->BorrowStore(gmx,mt);
#endif
}

GeneralMatrix* GetSubMatrix::Evaluate(MatrixType mt)
{
   REPORT
   Tracer tr("SubMatrix(evaluate)");
   gm = ((BaseMatrix*&)bm)->Evaluate();
   if (row_number < 0) row_number = gm->Nrows();
   if (col_number < 0) col_number = gm->Ncols();
   if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
      Throw(SubMatrixDimensionException());
   if (!mt) mt = MatrixType::Rt;
   GeneralMatrix* gmx = mt.New(row_number, col_number, this);
   int i = row_number;
   MatrixRow mr(gm, LoadOnEntry, row_skip); 
   MatrixRow mrx(gmx, StoreOnExit+DirectPart);
   MatrixRowCol sub;
   while (i--)
   {
      mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
      mrx.Copy(sub); mrx.Next(); mr.Next();
   }
   gmx->ReleaseAndDelete(); gm->tDelete();
#ifdef TEMPS_DESTROYED_QUICKLY
   delete this;
#endif
   return gmx;
}   


GeneralMatrix* ReturnMatrixX::Evaluate(MatrixType mt)
{
#ifdef TEMPS_DESTROYED_QUICKLY
   GeneralMatrix* gmx = gm; delete this; return gmx->Evaluate(mt);
#else
   return gm->Evaluate(mt);
#endif
}


void GeneralMatrix::Add(GeneralMatrix* gm1, Real f)
{
   REPORT
   Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
   while (i--)
   { *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; }
   i = storage & 3; while (i--) *s++ = *s1++ + f;
}
   
void GeneralMatrix::Add(Real f)
{
   REPORT
   Real* s=store; int i=(storage >> 2);
   while (i--) { *s++ += f; *s++ += f; *s++ += f; *s++ += f; }
   i = storage & 3; while (i--) *s++ += f;
}
   
void GeneralMatrix::Negate(GeneralMatrix* gm1)
{
   // change sign of elements
   REPORT
   Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
   while (i--)
   { *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); }
   i = storage & 3; while(i--) *s++ = -(*s1++);
}
   
void GeneralMatrix::Negate()
{
   REPORT
   Real* s=store; int i=(storage >> 2);
   while (i--)
   { *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; }
   i = storage & 3; while(i--) { *s = -(*s); s++; }
}
   
void GeneralMatrix::Multiply(GeneralMatrix* gm1, Real f)
{
   REPORT
   Real* s1=gm1->store; Real* s=store;  int i=(storage >> 2);
   while (i--)
   { *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; }
   i = storage & 3; while (i--) *s++ = *s1++ * f;
}
   
void GeneralMatrix::Multiply(Real f)
{
   REPORT
   Real* s=store; int i=(storage >> 2);
   while (i--) { *s++ *= f; *s++ *= f; *s++ *= f; *s++ *= f; }
   i = storage & 3; while (i--) *s++ *= f;
}