newmat4.cpp

非常好用的用C编写的矩阵类,可在不同编译器下编译使用.

   if (nr < nrows_val)
   {
      REPORT
      SquareMatrix X = sym_submatrix(1,nr);
      swap(X);
   }
   else if (nr > nrows_val)
   {
      REPORT
      SquareMatrix X(nr); X = 0;
      X.sym_submatrix(1,nrows_val) = *this;
      swap(X);
   }
}

void SquareMatrix::resize_keep(int nr, int nc)
{
   Tracer tr("SquareMatrix::resize_keep 2");
   REPORT
   if (nr != nc) Throw(NotSquareException(*this));
   resize_keep(nr);
}
 

void SymmetricMatrix::resize_keep(int nr)
{
   Tracer tr("SymmetricMatrix::resize_keep");
   if (nr < nrows_val)
   {
      REPORT
      SymmetricMatrix X = sym_submatrix(1,nr);
      swap(X);
   }
   else if (nr > nrows_val)
   {
      REPORT
      SymmetricMatrix X(nr); X = 0;
      X.sym_submatrix(1,nrows_val) = *this;
      swap(X);
   }
} 

void UpperTriangularMatrix::resize_keep(int nr)
{
   Tracer tr("UpperTriangularMatrix::resize_keep");
   if (nr < nrows_val)
   {
      REPORT
      UpperTriangularMatrix X = sym_submatrix(1,nr);
      swap(X);
   }
   else if (nr > nrows_val)
   {
      REPORT
      UpperTriangularMatrix X(nr); X = 0;
      X.sym_submatrix(1,nrows_val) = *this;
      swap(X);
   }
} 

void LowerTriangularMatrix::resize_keep(int nr)
{
   Tracer tr("LowerTriangularMatrix::resize_keep");
   if (nr < nrows_val)
   {
      REPORT
      LowerTriangularMatrix X = sym_submatrix(1,nr);
      swap(X);
   }
   else if (nr > nrows_val)
   {
      REPORT
      LowerTriangularMatrix X(nr); X = 0;
      X.sym_submatrix(1,nrows_val) = *this;
      swap(X);
   }
} 

void DiagonalMatrix::resize_keep(int nr)
{
   Tracer tr("DiagonalMatrix::resize_keep");
   if (nr < nrows_val)
   {
      REPORT
      DiagonalMatrix X = sym_submatrix(1,nr);
      swap(X);
   }
   else if (nr > nrows_val)
   {
      REPORT
      DiagonalMatrix X(nr); X = 0;
      X.sym_submatrix(1,nrows_val) = *this;
      swap(X);
   }
} 

void RowVector::resize_keep(int nc)
{
   Tracer tr("RowVector::resize_keep");
   if (nc < ncols_val)
   {
      REPORT
      RowVector X = columns(1,nc);
      swap(X);
   }
   else if (nc > ncols_val)
   {
      REPORT
      RowVector X(nc); X = 0;
      X.columns(1,ncols_val) = *this;
      swap(X);
   }
}

void RowVector::resize_keep(int nr, int nc)
{
   Tracer tr("RowVector::resize_keep 2");
   REPORT
   if (nr != 1) Throw(VectorException(*this));
   resize_keep(nc);
}

void ColumnVector::resize_keep(int nr)
{
   Tracer tr("ColumnVector::resize_keep");
   if (nr < nrows_val)
   {
      REPORT
      ColumnVector X = rows(1,nr);
      swap(X);
   }
   else if (nr > nrows_val)
   {
      REPORT
      ColumnVector X(nr); X = 0;
      X.rows(1,nrows_val) = *this;
      swap(X);
   }
} 

void ColumnVector::resize_keep(int nr, int nc)
{
   Tracer tr("ColumnVector::resize_keep 2");
   REPORT
   if (nc != 1) Throw(VectorException(*this));
   resize_keep(nr);
}


/*
void GeneralMatrix::resizeForAdd(const GeneralMatrix& A, const GeneralMatrix&)
{ REPORT resize(A); }

void GeneralMatrix::resizeForSP(const GeneralMatrix& A, const GeneralMatrix&)
{ REPORT resize(A); }


// ************************* SameStorageType ******************************

// SameStorageType checks A and B have same storage type including bandwidth
// It does not check same dimensions since we assume this is already done

bool GeneralMatrix::SameStorageType(const GeneralMatrix& A) const
{
   REPORT
   return type() == A.type();
}
*/

// ******************* manipulate types, storage **************************/

int GeneralMatrix::search(const BaseMatrix* s) const
{ REPORT return (s==this) ? 1 : 0; }

int GenericMatrix::search(const BaseMatrix* s) const
{ REPORT return gm->search(s); }

int MultipliedMatrix::search(const BaseMatrix* s) const
{ REPORT return bm1->search(s) + bm2->search(s); }

int ShiftedMatrix::search(const BaseMatrix* s) const
{ REPORT return bm->search(s); }

int NegatedMatrix::search(const BaseMatrix* s) const
{ REPORT return bm->search(s); }

int ReturnMatrix::search(const BaseMatrix* s) const
{ REPORT return (s==gm) ? 1 : 0; }

MatrixType Matrix::type() const { return MatrixType::Rt; }
MatrixType SquareMatrix::type() const { return MatrixType::Sq; }
MatrixType SymmetricMatrix::type() const { return MatrixType::Sm; }
MatrixType UpperTriangularMatrix::type() const { return MatrixType::UT; }
MatrixType LowerTriangularMatrix::type() const { return MatrixType::LT; }
MatrixType DiagonalMatrix::type() const { return MatrixType::Dg; }
MatrixType RowVector::type() const { return MatrixType::RV; }
MatrixType ColumnVector::type() const { return MatrixType::CV; }
MatrixType CroutMatrix::type() const { return MatrixType::Ct; }
MatrixType BandMatrix::type() const { return MatrixType::BM; }
MatrixType UpperBandMatrix::type() const { return MatrixType::UB; }
MatrixType LowerBandMatrix::type() const { return MatrixType::LB; }
MatrixType SymmetricBandMatrix::type() const { return MatrixType::SB; }

MatrixType IdentityMatrix::type() const { return MatrixType::Id; }



MatrixBandWidth BaseMatrix::bandwidth() const { REPORT return -1; }
MatrixBandWidth DiagonalMatrix::bandwidth() const { REPORT return 0; }
MatrixBandWidth IdentityMatrix::bandwidth() const { REPORT return 0; }

MatrixBandWidth UpperTriangularMatrix::bandwidth() const
   { REPORT return MatrixBandWidth(0,-1); }

MatrixBandWidth LowerTriangularMatrix::bandwidth() const
   { REPORT return MatrixBandWidth(-1,0); }

MatrixBandWidth BandMatrix::bandwidth() const
   { REPORT return MatrixBandWidth(lower_val,upper_val); }

MatrixBandWidth BandLUMatrix::bandwidth() const
   { REPORT return MatrixBandWidth(m1,m2); }
   
MatrixBandWidth GenericMatrix::bandwidth()const
   { REPORT return gm->bandwidth(); }

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

MatrixBandWidth SPMatrix::bandwidth() const
   { REPORT return gm1->bandwidth().minimum(gm2->bandwidth()); }

MatrixBandWidth KPMatrix::bandwidth() const
{
   int lower, upper;
   MatrixBandWidth bw1 = gm1->bandwidth(), bw2 = gm2->bandwidth();
   if (bw1.Lower() < 0)
   {
      if (bw2.Lower() < 0) { REPORT lower = -1; }
      else { REPORT lower = bw2.Lower() + (gm1->Nrows() - 1) * gm2->Nrows(); }
   }
   else
   {
      if (bw2.Lower() < 0)
         { REPORT lower = (1 + bw1.Lower()) * gm2->Nrows() - 1; }
      else { REPORT lower = bw2.Lower() + bw1.Lower() * gm2->Nrows(); }
   }
   if (bw1.Upper() < 0)
   {
      if (bw2.Upper() < 0) { REPORT upper = -1; }
      else { REPORT upper = bw2.Upper() + (gm1->Nrows() - 1) * gm2->Nrows(); }
   }
   else
   {
      if (bw2.Upper() < 0)
         { REPORT upper = (1 + bw1.Upper()) * gm2->Nrows() - 1; }
      else { REPORT upper = bw2.Upper() + bw1.Upper() * gm2->Nrows(); }
   }
   return MatrixBandWidth(lower, upper);
}

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

MatrixBandWidth ConcatenatedMatrix::bandwidth() const { REPORT return -1; }

MatrixBandWidth SolvedMatrix::bandwidth() const
{
   if (+gm1->type() & MatrixType::Diagonal)
      { REPORT return gm2->bandwidth(); }
   else { REPORT return -1; }
}

MatrixBandWidth ScaledMatrix::bandwidth() const
   { REPORT return gm->bandwidth(); }

MatrixBandWidth NegatedMatrix::bandwidth() const
   { REPORT return gm->bandwidth(); }

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

MatrixBandWidth InvertedMatrix::bandwidth() const
{
   if (+gm->type() & MatrixType::Diagonal)
      { REPORT return MatrixBandWidth(0,0); }
   else { REPORT return -1; }
}

MatrixBandWidth RowedMatrix::bandwidth() const { REPORT return -1; }
MatrixBandWidth ColedMatrix::bandwidth() const { REPORT return -1; }
MatrixBandWidth DiagedMatrix::bandwidth() const { REPORT return 0; }
MatrixBandWidth MatedMatrix::bandwidth() const { REPORT return -1; }
MatrixBandWidth ReturnMatrix::bandwidth() const
   { REPORT return gm->bandwidth(); }

MatrixBandWidth GetSubMatrix::bandwidth() const
{

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

// ********************** the memory managment tools **********************/

//  Rules regarding tDelete, reuse, GetStore, BorrowStore
//    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

//  Meaning of tag_val
//    tag_val = -1          memory cannot be reused (default situation)
//    tag_val = -2          memory has been borrowed from another matrix
//                               (don't change values)
//    tag_val = i > 0       delete or reuse memory after i operations
//    tag_val = 0           like value 1 but matrix was created by new
//                               so delete it

void GeneralMatrix::tDelete()
{
   if (tag_val<0)
   {
      if (tag_val<-1) { REPORT store = 0; delete this; return; }  // borrowed
      else { REPORT return; }   // not a temporary matrix - leave alone
   }
   if (tag_val==1)
   {
      if (store)
      {
         REPORT  MONITOR_REAL_DELETE("Free   (tDelete)",storage,store)
         delete [] store;
      }
      MiniCleanUp(); return;                           // CleanUp
   }
   if (tag_val==0) { REPORT delete this; return; }

   REPORT tag_val--; return;
}

void newmat_block_copy(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++;
}

bool GeneralMatrix::reuse()
{
   if (tag_val < -1)                 // borrowed storage
   {
      if (storage)
      {
         REPORT
         Real* s = new Real [storage]; MatrixErrorNoSpace(s);
         MONITOR_REAL_NEW("Make     (reuse)",storage,s)
         newmat_block_copy(storage, store, s); store = s;
      }
      else { REPORT MiniCleanUp(); }                // CleanUp
      tag_val = 0; return true;
   }
   if (tag_val < 0 ) { REPORT return false; }
   if (tag_val <= 1 )  { REPORT return true; }
   REPORT tag_val--; return false;
}

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

void GeneralMatrix::GetMatrix(const GeneralMatrix* gmx)
{
   REPORT  tag_val=-1; nrows_val=gmx->Nrows(); ncols_val=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 *gmx point to storage of *this if tag_val==-1.
{
   if (!mt)
   {
      if (tag_val == -1) { REPORT gmx->tag_val = -2; gmx->store = store; }
      else { REPORT gmx->tag_val = 0; gmx->store = GetStore(); }
   }
   else if (Compare(gmx->type(),mt))
   { REPORT  gmx->tag_val = 0; gmx->store = GetStore(); }
   else
   {
      REPORT gmx->tag_val = -2; gmx->store = store;
      gmx = gmx->Evaluate(mt); gmx->tag_val = 0; tDelete();
   }

   return gmx;
}

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