newmat6.cpp

各种矩阵算法库。支持UpperTriangularMatrix,LowerTriangularMatrix, DiagonalMatrix, SymmetricMatrix, BandMatrix,U

   Eq(X,MatrixType::RV);
   if (nrows!=1)
      { Tracer tr("RowVector(=)"); Throw(VectorException(*this)); }
}

void ColumnVector::operator=(const BaseMatrix& X)
{
   REPORT //CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::CV);
   if (ncols!=1)
      { Tracer tr("ColumnVector(=)"); Throw(VectorException(*this)); }
}

void SymmetricMatrix::operator=(const BaseMatrix& X)
{
   REPORT // CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::Sm);
}

void UpperTriangularMatrix::operator=(const BaseMatrix& X)
{
   REPORT //CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::UT);
}

void LowerTriangularMatrix::operator=(const BaseMatrix& X)
{
   REPORT //CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::LT);
}

void DiagonalMatrix::operator=(const BaseMatrix& X)
{
   REPORT // CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::Dg);
}

void GeneralMatrix::operator<<(const Real* r)
{
   REPORT
   int i = storage; Real* s=store;
   while(i--) *s++ = *r++;
}


void GenericMatrix::operator=(const GenericMatrix& bmx)
{
   if (&bmx != this) { REPORT if (gm) delete gm; gm = bmx.gm->Image();}
   else { REPORT }
   gm->Protect();
}

void GenericMatrix::operator=(const BaseMatrix& bmx)
{
   if (gm)
   {
      int counter=bmx.search(gm);
      if (counter==0) { REPORT delete gm; gm=0; }
      else { REPORT gm->Release(counter); }
   }
   else { REPORT }
   GeneralMatrix* gmx = ((BaseMatrix&)bmx).Evaluate();
   if (gmx != gm) { REPORT if (gm) delete gm; gm = gmx->Image(); }
   else { REPORT }
   gm->Protect();
}


/*************************** += etc ***************************************/

// will also need versions for SubMatrix



// GeneralMatrix operators

void GeneralMatrix::operator+=(const BaseMatrix& X)
{
   REPORT
   Tracer tr("GeneralMatrix::operator+=");
   // MatrixConversionCheck mcc;
   Protect();                                   // so it cannot get deleted
						// during Evaluate
   GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
   AddedMatrix* am = new AddedMatrix(this,gm);
   MatrixErrorNoSpace(am);
   if (gm==this) Release(2); else Release();
   Eq2(*am,Type());
#else
   AddedMatrix am(this,gm);
   if (gm==this) Release(2); else Release();
   Eq2(am,Type());
#endif
}

void GeneralMatrix::operator-=(const BaseMatrix& X)
{
   REPORT
   Tracer tr("GeneralMatrix::operator-=");
   // MatrixConversionCheck mcc;
   Protect();                                   // so it cannot get deleted
						// during Evaluate
   GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
   SubtractedMatrix* am = new SubtractedMatrix(this,gm);
   MatrixErrorNoSpace(am);
   if (gm==this) Release(2); else Release();
   Eq2(*am,Type());
#else
   SubtractedMatrix am(this,gm);
   if (gm==this) Release(2); else Release();
   Eq2(am,Type());
#endif
}

void GeneralMatrix::operator*=(const BaseMatrix& X)
{
   REPORT
   Tracer tr("GeneralMatrix::operator*=");
   // MatrixConversionCheck mcc;
   Protect();                                   // so it cannot get deleted
						// during Evaluate
   GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
   MultipliedMatrix* am = new MultipliedMatrix(this,gm);
   MatrixErrorNoSpace(am);
   if (gm==this) Release(2); else Release();
   Eq2(*am,Type());
#else
   MultipliedMatrix am(this,gm);
   if (gm==this) Release(2); else Release();
   Eq2(am,Type());
#endif
}

void GeneralMatrix::operator|=(const BaseMatrix& X)
{
   REPORT
   Tracer tr("GeneralMatrix::operator|=");
   // MatrixConversionCheck mcc;
   Protect();                                   // so it cannot get deleted
						// during Evaluate
   GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
   ConcatenatedMatrix* am = new ConcatenatedMatrix(this,gm);
   MatrixErrorNoSpace(am);
   if (gm==this) Release(2); else Release();
   Eq2(*am,Type());
#else
   ConcatenatedMatrix am(this,gm);
   if (gm==this) Release(2); else Release();
   Eq2(am,Type());
#endif
}

void GeneralMatrix::operator&=(const BaseMatrix& X)
{
   REPORT
   Tracer tr("GeneralMatrix::operator&=");
   // MatrixConversionCheck mcc;
   Protect();                                   // so it cannot get deleted
						// during Evaluate
   GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
   StackedMatrix* am = new StackedMatrix(this,gm);
   MatrixErrorNoSpace(am);
   if (gm==this) Release(2); else Release();
   Eq2(*am,Type());
#else
   StackedMatrix am(this,gm);
   if (gm==this) Release(2); else Release();
   Eq2(am,Type());
#endif
}

void GeneralMatrix::operator+=(Real r)
{
   REPORT
   Tracer tr("GeneralMatrix::operator+=(Real)");
   // MatrixConversionCheck mcc;
#ifdef TEMPS_DESTROYED_QUICKLY
   ShiftedMatrix* am = new ShiftedMatrix(this,r);
   MatrixErrorNoSpace(am);
   Release(); Eq2(*am,Type());
#else
   ShiftedMatrix am(this,r);
   Release(); Eq2(am,Type());
#endif
}

void GeneralMatrix::operator*=(Real r)
{
   REPORT
   Tracer tr("GeneralMatrix::operator*=(Real)");
   // MatrixConversionCheck mcc;
#ifdef TEMPS_DESTROYED_QUICKLY
   ScaledMatrix* am = new ScaledMatrix(this,r);
   MatrixErrorNoSpace(am);
   Release(); Eq2(*am,Type());
#else
   ScaledMatrix am(this,r);
   Release(); Eq2(am,Type());
#endif
}


// Generic matrix operators

void GenericMatrix::operator+=(const BaseMatrix& X)
{
   REPORT
   Tracer tr("GenericMatrix::operator+=");
   if (!gm) Throw(ProgramException("GenericMatrix is null"));
   gm->Protect();            // so it cannot get deleted during Evaluate
   GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
   AddedMatrix* am = new AddedMatrix(gm,gmx);
   MatrixErrorNoSpace(am);
   if (gmx==gm) gm->Release(2); else gm->Release();
   GeneralMatrix* gmy = am->Evaluate();
#else
   AddedMatrix am(gm,gmx);
   if (gmx==gm) gm->Release(2); else gm->Release();
   GeneralMatrix* gmy = am.Evaluate();
#endif
   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
   else { REPORT }
   gm->Protect();
}

void GenericMatrix::operator-=(const BaseMatrix& X)
{
   REPORT
   Tracer tr("GenericMatrix::operator-=");
   if (!gm) Throw(ProgramException("GenericMatrix is null"));
   gm->Protect();            // so it cannot get deleted during Evaluate
   GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
   SubtractedMatrix* am = new SubtractedMatrix(gm,gmx);
   MatrixErrorNoSpace(am);
   if (gmx==gm) gm->Release(2); else gm->Release();
   GeneralMatrix* gmy = am->Evaluate();
#else
   SubtractedMatrix am(gm,gmx);
   if (gmx==gm) gm->Release(2); else gm->Release();
   GeneralMatrix* gmy = am.Evaluate();
#endif
   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
   else { REPORT }
   gm->Protect();
}

void GenericMatrix::operator*=(const BaseMatrix& X)
{
   REPORT
   Tracer tr("GenericMatrix::operator*=");
   if (!gm) Throw(ProgramException("GenericMatrix is null"));
   gm->Protect();            // so it cannot get deleted during Evaluate
   GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
   MultipliedMatrix* am = new MultipliedMatrix(gm,gmx);
   MatrixErrorNoSpace(am);
   if (gmx==gm) gm->Release(2); else gm->Release();
   GeneralMatrix* gmy = am->Evaluate();
#else
   MultipliedMatrix am(gm,gmx);
   if (gmx==gm) gm->Release(2); else gm->Release();
   GeneralMatrix* gmy = am.Evaluate();
#endif
   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
   else { REPORT }
   gm->Protect();
}

void GenericMatrix::operator|=(const BaseMatrix& X)
{
   REPORT
   Tracer tr("GenericMatrix::operator|=");
   if (!gm) Throw(ProgramException("GenericMatrix is null"));
   gm->Protect();            // so it cannot get deleted during Evaluate
   GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
   ConcatenatedMatrix* am = new ConcatenatedMatrix(gm,gmx);
   MatrixErrorNoSpace(am);
   if (gmx==gm) gm->Release(2); else gm->Release();
   GeneralMatrix* gmy = am->Evaluate();
#else
   ConcatenatedMatrix am(gm,gmx);
   if (gmx==gm) gm->Release(2); else gm->Release();
   GeneralMatrix* gmy = am.Evaluate();
#endif
   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
   else { REPORT }
   gm->Protect();
}

void GenericMatrix::operator&=(const BaseMatrix& X)
{
   REPORT
   Tracer tr("GenericMatrix::operator&=");
   if (!gm) Throw(ProgramException("GenericMatrix is null"));
   gm->Protect();            // so it cannot get deleted during Evaluate
   GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
   StackedMatrix* am = new StackedMatrix(gm,gmx);
   MatrixErrorNoSpace(am);
   if (gmx==gm) gm->Release(2); else gm->Release();
   GeneralMatrix* gmy = am->Evaluate();
#else
   StackedMatrix am(gm,gmx);
   if (gmx==gm) gm->Release(2); else gm->Release();
   GeneralMatrix* gmy = am.Evaluate();
#endif
   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
   else { REPORT }
   gm->Protect();
}

void GenericMatrix::operator+=(Real r)
{
   REPORT
   Tracer tr("GenericMatrix::operator+= (Real)");
   if (!gm) Throw(ProgramException("GenericMatrix is null"));
#ifdef TEMPS_DESTROYED_QUICKLY
   ShiftedMatrix* am = new ShiftedMatrix(gm,r);
   MatrixErrorNoSpace(am);
   gm->Release();
   GeneralMatrix* gmy = am->Evaluate();
#else
   ShiftedMatrix am(gm,r);
   gm->Release();
   GeneralMatrix* gmy = am.Evaluate();
#endif
   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
   else { REPORT }
   gm->Protect();
}

void GenericMatrix::operator*=(Real r)
{
   REPORT
   Tracer tr("GenericMatrix::operator*= (Real)");
   if (!gm) Throw(ProgramException("GenericMatrix is null"));
#ifdef TEMPS_DESTROYED_QUICKLY
   ScaledMatrix* am = new ScaledMatrix(gm,r);
   MatrixErrorNoSpace(am);
   gm->Release();
   GeneralMatrix* gmy = am->Evaluate();
#else
   ScaledMatrix am(gm,r);
   gm->Release();
   GeneralMatrix* gmy = am.Evaluate();
#endif
   if (gmy != gm) { REPORT delete gm; gm = gmy->Image(); }
   else { REPORT }
   gm->Protect();
}


/************************* element access *********************************/

Real& Matrix::element(int m, int n)
{
   if (m<0 || m>= nrows || n<0 || n>= ncols)
      Throw(IndexException(m,n,*this,true));
   return store[m*ncols+n];
}

Real& SymmetricMatrix::element(int m, int n)
{
   if (m<0 || n<0 || m >= nrows || n>=ncols)
      Throw(IndexException(m,n,*this,true));
   if (m>=n) return store[tristore(m)+n];
   else return store[tristore(n)+m];
}

Real& UpperTriangularMatrix::element(int m, int n)
{
   if (m<0 || n<m || n>=ncols)
      Throw(IndexException(m,n,*this,true));
   return store[m*ncols+n-tristore(m)];
}

Real& LowerTriangularMatrix::element(int m, int n)
{
   if (n<0 || m<n || m>=nrows)
      Throw(IndexException(m,n,*this,true));
   return store[tristore(m)+n];
}

Real& DiagonalMatrix::element(int m, int n)
{
   if (n<0 || m!=n || m>=nrows || n>=ncols)
      Throw(IndexException(m,n,*this,true));
   return store[n];
}

Real& DiagonalMatrix::element(int m)
{
   if (m<0 || m>=nrows) Throw(IndexException(m,*this,true));
   return store[m];
}

Real& ColumnVector::element(int m)
{
   if (m<0 || m>= nrows) Throw(IndexException(m,*this,true));
   return store[m];
}

Real& RowVector::element(int n)
{
   if (n<0 || n>= ncols)  Throw(IndexException(n,*this,true));
   return store[n];
}

Real& BandMatrix::element(int m, int n)
{
   int w = upper+lower+1; int i = lower+n-m;
   if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
      Throw(IndexException(m,n,*this,true));
   return store[w*m+i];
}

Real& UpperBandMatrix::element(int m, int n)
{
   int w = upper+1; int i = n-m;
   if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
      Throw(IndexException(m,n,*this,true));
   return store[w*m+i];
}

Real& LowerBandMatrix::element(int m, int n)
{
   int w = lower+1; int i = lower+n-m;
   if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
      Throw(IndexException(m,n,*this,true));
   return store[w*m+i];
}

Real& SymmetricBandMatrix::element(int m, int n)
{
   int w = lower+1;
   if (m>=n)
   {
      int i = lower+n-m;
      if ( m>=nrows || n<0 || i<0 )
         Throw(IndexException(m,n,*this,true));
      return store[w*m+i];
   }
   else
   {
      int i = lower+m-n;
      if ( n>=nrows || m<0 || i<0 )
         Throw(IndexException(m,n,*this,true));
      return store[w*n+i];
   }
}

#ifdef use_namespace
}
#endif