newmat7.cpp

矩阵计算库

      Throw(IncompatibleDimensionsException());
#ifdef __GNUG__
   int c1,c2; sm->SelectVersion(mtx,c1,c2);
#else
   Boolean c1,c2; sm->SelectVersion(mtx,c1,c2); // causes problems for g++
#endif
   if (c1 && c2)
   {
      if (gm1->reuse())
      { REPORT Subtract(gm1,gm2); gm2->tDelete(); return gm1; }
      else if (gm2->reuse()) { REPORT ReverseSubtract(gm2,gm1); return gm2; }
      else
      {
         REPORT
	 GeneralMatrix* gmx = gm1->Type().New(nr,nc,sm);
         gmx->ReleaseAndDelete(); Subtract(gmx,gm1,gm2); return gmx;
      }
   }
   else
   {
      if ( c1 && gm1->reuse() )
      { REPORT  SubtractDS(gm1,gm2); gm2->tDelete(); return gm1; }
      else if ( c2 && gm2->reuse() )
      {
         REPORT
         ReverseSubtractDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2;
      }
      else
      {
         REPORT
	 GeneralMatrix* gmx = mtx.New(nr,nc,sm); SubtractDS(gmx,gm1,gm2);
	 if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
	 gmx->ReleaseAndDelete(); return gmx;
      }
   }
}

static GeneralMatrix* GeneralMult1(GeneralMatrix* gm1, GeneralMatrix* gm2,
   MultipliedMatrix* mm, MatrixType mtx)
{
   REPORT
   Tracer tr("GeneralMult1");
   int nr=gm1->Nrows(); int nc=gm2->Ncols();
   if (gm1->Ncols() !=gm2->Nrows())
      Throw(IncompatibleDimensionsException());
   GeneralMatrix* gmx = mtx.New(nr,nc,mm);

   MatrixCol mcx(gmx, StoreOnExit+DirectPart);
   MatrixCol mc2(gm2, LoadOnEntry);
   while (nc--)
   {
      MatrixRow mr1(gm1, LoadOnEntry, mcx.Skip());
      Real* el = mcx();                              // pointer to an element
      int n = mcx.Storage();
      while (n--) { *(el++) = DotProd(mr1,mc2); mr1.Next(); }
      mc2.Next(); mcx.Next();
   }
   gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
}

static GeneralMatrix* GeneralMult2(GeneralMatrix* gm1, GeneralMatrix* gm2,
   MultipliedMatrix* mm, MatrixType mtx)
{
   // version that accesses by row only - not good for thin matrices
   // or column vectors in right hand term. Needs fixing
   REPORT
   Tracer tr("GeneralMult2");
   int nr=gm1->Nrows(); int nc=gm2->Ncols();
   if (gm1->Ncols() !=gm2->Nrows())
      Throw(IncompatibleDimensionsException());
   GeneralMatrix* gmx = mtx.New(nr,nc,mm);

   Real* el = gmx->Store(); int n = gmx->Storage();
   while (n--) *el++ = 0.0;
   MatrixRow mrx(gmx, LoadOnEntry+StoreOnExit+DirectPart);
   MatrixRow mr1(gm1, LoadOnEntry);
   while (nr--)
   {
      MatrixRow mr2(gm2, LoadOnEntry, mr1.Skip());
      el = mr1();                              // pointer to an element
      n = mr1.Storage();
      while (n--) { mrx.AddScaled(mr2, *el++); mr2.Next(); }
      mr1.Next(); mrx.Next();
   }
   gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
}

static GeneralMatrix* mmMult(GeneralMatrix* gm1, GeneralMatrix* gm2)
{
   // matrix multiplication for type Matrix only
   REPORT
   Tracer tr("MatrixMult");

   int nr=gm1->Nrows(); int ncr=gm1->Ncols(); int nc=gm2->Ncols();
   if (ncr != gm2->Nrows()) Throw(IncompatibleDimensionsException());

   Matrix* gm = new Matrix(nr,nc); MatrixErrorNoSpace(gm);

   Real* s1=gm1->Store(); Real* s2=gm2->Store(); Real* s=gm->Store();
   
   if (ncr)
   {
      while (nr--)
      {
         Real* s2x = s2; int j = ncr;
         Real* sx = s; Real f = *s1++; int k = nc;
         while (k--) *sx++ = f * *s2x++;
         while (--j)
            { sx = s; f = *s1++; k = nc; while (k--) *sx++ += f * *s2x++; }
         s = sx;
      }
   }
   else *gm = 0.0;

   gm->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gm;
}

static GeneralMatrix* GeneralMult(GeneralMatrix* gm1, GeneralMatrix* gm2,
   MultipliedMatrix* mm, MatrixType mtx)
{
   if ( Rectangular(gm1->Type(), gm2->Type(), mtx)) return mmMult(gm1, gm2);
   else
   {
      Compare(gm1->Type() * gm2->Type(),mtx);
      int nr = gm2->Nrows(); int nc = gm2->Ncols();
      if (nc <= 5 && nr > nc) return GeneralMult1(gm1, gm2, mm, mtx);
      else return GeneralMult2(gm1, gm2, mm, mtx);
   }
}

#ifdef __GNUG__
void SPMatrix::SelectVersion
   (MatrixType mtx, int& c1, int& c2) const
#else
void SPMatrix::SelectVersion
   (MatrixType mtx, Boolean& c1, Boolean& c2) const
#endif
// for determining version of SP routines
// will need to modify if further matrix structures are introduced
{
   MatrixBandWidth bw1 = gm1->BandWidth();
   MatrixBandWidth bw2 = gm2->BandWidth();
   MatrixBandWidth bwx(-1);  if (mtx.IsBand()) bwx = bw1.minimum(bw2);
   c1 = (mtx == gm1->Type()) && (bwx == bw1);
   c2 = (mtx == gm2->Type()) && (bwx == bw2);
}

static GeneralMatrix* GeneralSP(GeneralMatrix* gm1, GeneralMatrix* gm2,
   SPMatrix* am, MatrixType mtx)
{
   Tracer tr("GeneralSP");
   int nr=gm1->Nrows(); int nc=gm1->Ncols();
   if (nr!=gm2->Nrows() || nc!=gm2->Ncols()) 
      Throw(IncompatibleDimensionsException());
   Compare(gm1->Type().SP(gm2->Type()),mtx);
#ifdef __GNUG__
   int c1,c2; am->SelectVersion(mtx,c1,c2);
#else
   Boolean c1,c2; am->SelectVersion(mtx,c1,c2); // causes problems for g++
#endif
   if (c1 && c2)
   {
      if (gm1->reuse()) { REPORT SP(gm1,gm2); gm2->tDelete(); return gm1; }
      else if (gm2->reuse()) { REPORT SP(gm2,gm1); return gm2; }
      else
      {
         REPORT GeneralMatrix* gmx = gm1->Type().New(nr,nc,am);
         gmx->ReleaseAndDelete(); SP(gmx,gm1,gm2); return gmx;
      }
   }
   else
   {
      if (c1 && gm1->reuse() )               // must have type test first
      { REPORT SPDS(gm1,gm2); gm2->tDelete(); return gm1; }
      else if (c2 && gm2->reuse() )
      { REPORT SPDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2; }
      else
      {
         REPORT
	 GeneralMatrix* gmx = mtx.New(nr,nc,am); SPDS(gmx,gm1,gm2);
	 if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
         gmx->ReleaseAndDelete(); return gmx;
      }
   }
}


static GeneralMatrix* GeneralSolv(GeneralMatrix* gm1, GeneralMatrix* gm2,
   BaseMatrix* sm, MatrixType mtx)
{
   REPORT
   Tracer tr("GeneralSolv");
   Compare(gm1->Type().i() * gm2->Type(),mtx);
   int nr=gm1->Nrows(); int nc=gm2->Ncols();
   if (gm1->Ncols() !=gm2->Nrows())
      Throw(IncompatibleDimensionsException());
   GeneralMatrix* gmx = mtx.New(nr,nc,sm); MatrixErrorNoSpace(gmx);
   Real* r = new Real [nr]; MatrixErrorNoSpace(r);
#ifndef ATandT
   MONITOR_REAL_NEW("Make   (GenSolv)",nr,r)
                           // deleted for ATandT, to balance deletion below
#endif
   GeneralMatrix* gms = gm1->MakeSolver();
   Try
   {
      MatrixCol mcx(gmx, r, StoreOnExit+DirectPart);   // copy to and from r
         // this must be inside Try so mcx is destroyed before gmx
      MatrixCol mc2(gm2, r, LoadOnEntry);
      while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
   }
   CatchAll
   {
      gms->tDelete();
      delete gmx;                   // <--------------------
      gm2->tDelete();
#ifndef ATandT
      MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
                          // ATandT version 2.1 gives an internal error
#endif
#ifdef Version21
      delete [] r;
#else
      delete [nr] r;
#endif
      Bounce;
   }
   gms->tDelete(); gmx->ReleaseAndDelete(); gm2->tDelete();
#ifndef ATandT
   MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
                          // ATandT version 2.1 gives an internal error
#endif
#ifdef Version21
   delete [] r;
#else
   delete [nr] r;
#endif
   return gmx;
}

GeneralMatrix* InvertedMatrix::Evaluate(MatrixType mtx)
{
   // Matrix Inversion - use solve routines
   Tracer tr("InvertedMatrix::Evaluate");
   REPORT
   gm=((BaseMatrix*&)bm)->Evaluate();
   int n = gm->Nrows(); DiagonalMatrix I(n); I=1.0;
#ifdef TEMPS_DESTROYED_QUICKLY
	GeneralMatrix* gmx;
	Try
	{
		Compare(gm->Type().i(),mtx);
		SkipConversionCheck SCC;              // otherwise inverting
														  // symmetric causes a problem
		gmx = GeneralSolv(gm,&I,this,mtx);
	}
   CatchAll { delete this; Bounce; }
   delete this; return gmx;
#else
	Compare(gm->Type().i(),mtx);
	SkipConversionCheck SCC;              // otherwise inverting
													  // symmetric causes a problem
   return GeneralSolv(gm,&I,this,mtx);
#endif   
}


/*************************** norm functions ****************************/

Real BaseMatrix::Norm1() const
{
   // maximum of sum of absolute values of a column
   REPORT
   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   int nc = gm->Ncols(); Real value = 0.0;
   MatrixCol mc(gm, LoadOnEntry);
   while (nc--)
      { Real v = mc.SumAbsoluteValue(); if (value < v) value = v; mc.Next(); }
   gm->tDelete(); return value;
}

Real BaseMatrix::NormInfinity() const
{
   // maximum of sum of absolute values of a row
   REPORT
   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
   int nr = gm->Nrows(); Real value = 0.0;
   MatrixRow mr(gm, LoadOnEntry);
   while (nr--)
      { Real v = mr.SumAbsoluteValue(); if (value < v) value = v; mr.Next(); }
   gm->tDelete(); return value;
}

/********************** Concatenation and stacking *************************/

GeneralMatrix* ConcatenatedMatrix::Evaluate(MatrixType mtx)
{
   REPORT
   Tracer tr("Concatenate");
#ifdef TEMPS_DESTROYED_QUICKLY
      Try
      {
         gm2 = ((BaseMatrix*&)bm2)->Evaluate();
         gm1 = ((BaseMatrix*&)bm1)->Evaluate();
         Compare(gm1->Type() | gm2->Type(),mtx);
         int nr=gm1->Nrows(); int nc = gm1->Ncols() + gm2->Ncols();
         if (nr != gm2->Nrows()) Throw(IncompatibleDimensionsException());
         GeneralMatrix* gmx = mtx.New(nr,nc,this);
         MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
         MatrixRow mr(gmx, StoreOnExit+DirectPart);
         while (nr--) { mr.ConCat(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
         gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); delete this;
         return gmx;
      }
      CatchAll  { delete this; Bounce; }
#ifndef UseExceptions
      return 0;
#endif 
#else
      gm2 = ((BaseMatrix*&)bm2)->Evaluate();
      gm1 = ((BaseMatrix*&)bm1)->Evaluate();
      Compare(gm1->Type() | gm2->Type(),mtx);
      int nr=gm1->Nrows(); int nc = gm1->Ncols() + gm2->Ncols();
      if (nr != gm2->Nrows()) Throw(IncompatibleDimensionsException());
      GeneralMatrix* gmx = mtx.New(nr,nc,this);
      MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
      MatrixRow mr(gmx, StoreOnExit+DirectPart);
      while (nr--) { mr.ConCat(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
      gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
#endif
}

GeneralMatrix* StackedMatrix::Evaluate(MatrixType mtx)
{
   REPORT
   Tracer tr("Stack");
#ifdef TEMPS_DESTROYED_QUICKLY
      Try
      {
         gm2 = ((BaseMatrix*&)bm2)->Evaluate();
         gm1 = ((BaseMatrix*&)bm1)->Evaluate();
         Compare(gm1->Type() & gm2->Type(),mtx);
         int nc=gm1->Ncols();
         int nr1 = gm1->Nrows(); int nr2 = gm2->Nrows();
         if (nc != gm2->Ncols()) Throw(IncompatibleDimensionsException());
         GeneralMatrix* gmx = mtx.New(nr1+nr2,nc,this);
         MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
         MatrixRow mr(gmx, StoreOnExit+DirectPart);
         while (nr1--) { mr.Copy(mr1); mr1.Next(); mr.Next(); }
         while (nr2--) { mr.Copy(mr2); mr2.Next(); mr.Next(); }
         gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); delete this;
         return gmx;
      }
      CatchAll  { delete this; Bounce; } 
#ifndef UseExceptions
      return 0;
#endif 
#else
      gm2 = ((BaseMatrix*&)bm2)->Evaluate();
      gm1 = ((BaseMatrix*&)bm1)->Evaluate();
      Compare(gm1->Type() & gm2->Type(),mtx);
      int nc=gm1->Ncols();
      int nr1 = gm1->Nrows(); int nr2 = gm2->Nrows();
      if (nc != gm2->Ncols()) Throw(IncompatibleDimensionsException());
      GeneralMatrix* gmx = mtx.New(nr1+nr2,nc,this);
      MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
      MatrixRow mr(gmx, StoreOnExit+DirectPart);
      while (nr1--) { mr.Copy(mr1); mr1.Next(); mr.Next(); }
      while (nr2--) { mr.Copy(mr2); mr2.Next(); mr.Next(); }
      gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
#endif
}