newmat7.cpp

C++矩阵算法库

   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))
   {
      REPORT
      return mmMult(gm1, gm2);
   }
   else
   {
      REPORT
      Compare(gm1->Type() * gm2->Type(),mtx);
      int nr = gm2->Nrows(); int nc = gm2->Ncols();
      if (nc <= 5 && nr > nc) { REPORT return GeneralMult1(gm1, gm2, mm, mtx); }
      else { REPORT return GeneralMult2(gm1, gm2, mm, mtx); }
   }
}

static GeneralMatrix* GeneralKP(GeneralMatrix* gm1, GeneralMatrix* gm2,
   KPMatrix* kp, MatrixType mtx)
{
   REPORT
   Tracer tr("GeneralKP");
   int nr1 = gm1->Nrows(); int nc1 = gm1->Ncols();
   int nr2 = gm2->Nrows(); int nc2 = gm2->Ncols();
   Compare((gm1->Type()).KP(gm2->Type()),mtx);
   GeneralMatrix* gmx = mtx.New(nr1*nr2, nc1*nc2, kp);
   MatrixRow mrx(gmx, LoadOnEntry+StoreOnExit+DirectPart);
   MatrixRow mr1(gm1, LoadOnEntry);
   for (int i = 1; i <= nr1; ++i)
   {
      MatrixRow mr2(gm2, LoadOnEntry);
      for (int j = 1; j <= nr2; ++j)
         { mrx.KP(mr1,mr2); mr2.Next(); mrx.Next(); }
      mr1.Next();
   }
   gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); 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();
   if (nr != gm1->Ncols()) Throw(NotSquareException(*gm1));
   int nc = gm2->Ncols();
   if (gm1->Ncols() != gm2->Nrows())
      Throw(IncompatibleDimensionsException(*gm1, *gm2));
   GeneralMatrix* gmx = mtx.New(nr,nc,sm); MatrixErrorNoSpace(gmx);
   Real* r = new Real [nr]; MatrixErrorNoSpace(r);
   MONITOR_REAL_NEW("Make   (GenSolv)",nr,r)
   GeneralMatrix* gms = gm1->MakeSolver();
   Try
   {

      MatrixColX mcx(gmx, r, StoreOnExit+DirectPart);   // copy to and from r
         // this must be inside Try so mcx is destroyed before gmx
      MatrixColX mc2(gm2, r, LoadOnEntry);
      while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
   }
   CatchAll
   {
      if (gms) gms->tDelete();
      delete gmx;                   // <--------------------
      gm2->tDelete();
      MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
                          // ATandT version 2.1 gives an internal error
      delete [] r;
      ReThrow;
   }
   gms->tDelete(); gmx->ReleaseAndDelete(); gm2->tDelete();
   MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
                          // ATandT version 2.1 gives an internal error
   delete [] r;
   return gmx;
}

// version for inverses - gm2 is identity
static GeneralMatrix* GeneralSolvI(GeneralMatrix* gm1, BaseMatrix* sm,
   MatrixType mtx)
{
   REPORT
   Tracer tr("GeneralSolvI");
   Compare(gm1->Type().i(),mtx);
   int nr = gm1->Nrows();
   if (nr != gm1->Ncols()) Throw(NotSquareException(*gm1));
   int nc = nr;
   // DiagonalMatrix I(nr); I = 1;
   IdentityMatrix I(nr);
   GeneralMatrix* gmx = mtx.New(nr,nc,sm); MatrixErrorNoSpace(gmx);
   Real* r = new Real [nr]; MatrixErrorNoSpace(r);
   MONITOR_REAL_NEW("Make   (GenSolvI)",nr,r)
   GeneralMatrix* gms = gm1->MakeSolver();
   Try
   {

      MatrixColX mcx(gmx, r, StoreOnExit+DirectPart);   // copy to and from r
         // this must be inside Try so mcx is destroyed before gmx
      MatrixColX mc2(&I, r, LoadOnEntry);
      while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
   }
   CatchAll
   {
      if (gms) gms->tDelete();
      delete gmx;
      MONITOR_REAL_DELETE("Delete (GenSolvI)",nr,r)
                          // ATandT version 2.1 gives an internal error
      delete [] r;
      ReThrow;
   }
   gms->tDelete(); gmx->ReleaseAndDelete();
   MONITOR_REAL_DELETE("Delete (GenSolvI)",nr,r)
                          // ATandT version 2.1 gives an internal error
   delete [] r;
   return gmx;
}

GeneralMatrix* InvertedMatrix::Evaluate(MatrixType mtx)
{
   // Matrix Inversion - use solve routines
   Tracer tr("InvertedMatrix::Evaluate");
   REPORT
   gm=((BaseMatrix*&)bm)->Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
	GeneralMatrix* gmx;
	Try { gmx = GeneralSolvI(gm,this,mtx); }
   CatchAll { delete this; ReThrow; }
   delete this; return gmx;
#else
   return GeneralSolvI(gm,this,mtx);
#endif
}

//*************************** New versions ************************

GeneralMatrix* AddedMatrix::Evaluate(MatrixType mtd)
{
   REPORT
   Tracer tr("AddedMatrix::Evaluate");
   gm1=((BaseMatrix*&)bm1)->Evaluate(); gm2=((BaseMatrix*&)bm2)->Evaluate();
   int nr=gm1->Nrows(); int nc=gm1->Ncols();
   if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
   {
      Try { Throw(IncompatibleDimensionsException(*gm1, *gm2)); }
      CatchAll
      {
         gm1->tDelete(); gm2->tDelete();
#ifdef TEMPS_DESTROYED_QUICKLY
         delete this;
#endif
         ReThrow;
      }
   }
   MatrixType mt1 = gm1->Type(), mt2 = gm2->Type(); MatrixType mts = mt1 + mt2;
   if (!mtd) { REPORT mtd = mts; }
   else if (!(mtd.DataLossOK || mtd >= mts))
   {
      REPORT
      gm1->tDelete(); gm2->tDelete();
#ifdef TEMPS_DESTROYED_QUICKLY
      delete this;
#endif
      Throw(ProgramException("Illegal Conversion", mts, mtd));
   }
   GeneralMatrix* gmx;
   bool c1 = (mtd == mt1), c2 = (mtd == mt2);
   if ( c1 && c2 && (gm1->SimpleAddOK(gm2) == 0) )
   {
      if (gm1->reuse()) { REPORT Add(gm1,gm2); gm2->tDelete(); gmx = gm1; }
      else if (gm2->reuse()) { REPORT Add(gm2,gm1); gmx = gm2; }
      else
      {
         REPORT
         // what if new throws an exception
         Try { gmx = mt1.New(nr,nc,this); }
         CatchAll
         {
#ifdef TEMPS_DESTROYED_QUICKLY
            delete this;
#endif
            ReThrow;
         }
         gmx->ReleaseAndDelete(); Add(gmx,gm1,gm2);
      }
   }
   else
   {
      if (c1 && c2)
      {
         short SAO = gm1->SimpleAddOK(gm2);
         if (SAO & 1) { REPORT c1 = false; }
         if (SAO & 2) { REPORT c2 = false; }
      }
      if (c1 && gm1->reuse() )               // must have type test first
         { REPORT AddDS(gm1,gm2); gm2->tDelete(); gmx = gm1; }
      else if (c2 && gm2->reuse() )
         { REPORT AddDS(gm2,gm1); if (!c1) gm1->tDelete(); gmx = gm2; }
      else
      {
         REPORT
         Try { gmx = mtd.New(nr,nc,this); }
         CatchAll
         {
            if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
#ifdef TEMPS_DESTROYED_QUICKLY
            delete this;
#endif
            ReThrow;
         }
         AddDS(gmx,gm1,gm2);
         if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
         gmx->ReleaseAndDelete();
      }
   }
#ifdef TEMPS_DESTROYED_QUICKLY
   delete this;
#endif
   return gmx;
}

GeneralMatrix* SubtractedMatrix::Evaluate(MatrixType mtd)
{
   REPORT
   Tracer tr("SubtractedMatrix::Evaluate");
   gm1=((BaseMatrix*&)bm1)->Evaluate(); gm2=((BaseMatrix*&)bm2)->Evaluate();
   int nr=gm1->Nrows(); int nc=gm1->Ncols();
   if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
   {
      Try { Throw(IncompatibleDimensionsException(*gm1, *gm2)); }
      CatchAll
      {
         gm1->tDelete(); gm2->tDelete();
#ifdef TEMPS_DESTROYED_QUICKLY
         delete this;
#endif
         ReThrow;
      }
   }
   MatrixType mt1 = gm1->Type(), mt2 = gm2->Type(); MatrixType mts = mt1 + mt2;
   if (!mtd) { REPORT mtd = mts; }
   else if (!(mtd.DataLossOK || mtd >= mts))
   {
      gm1->tDelete(); gm2->tDelete();
#ifdef TEMPS_DESTROYED_QUICKLY
      delete this;
#endif
      Throw(ProgramException("Illegal Conversion", mts, mtd));
   }
   GeneralMatrix* gmx;
   bool c1 = (mtd == mt1), c2 = (mtd == mt2);
   if ( c1 && c2 && (gm1->SimpleAddOK(gm2) == 0) )
   {
      if (gm1->reuse()) { REPORT Subtract(gm1,gm2); gm2->tDelete(); gmx = gm1; }
      else if (gm2->reuse()) { REPORT ReverseSubtract(gm2,gm1); gmx = gm2; }
      else
      {
         REPORT
         Try { gmx = mt1.New(nr,nc,this); }
         CatchAll
         {
#ifdef TEMPS_DESTROYED_QUICKLY
            delete this;
#endif
            ReThrow;
         }
         gmx->ReleaseAndDelete(); Subtract(gmx,gm1,gm2);
      }
   }
   else
   {
      if (c1 && c2)
      {
         short SAO = gm1->SimpleAddOK(gm2);
         if (SAO & 1) { REPORT c1 = false; }
         if (SAO & 2) { REPORT c2 = false; }
      }
      if (c1 && gm1->reuse() )               // must have type test first
         { REPORT SubtractDS(gm1,gm2); gm2->tDelete(); gmx = gm1; }
      else if (c2 && gm2->reuse() )
      {
         REPORT ReverseSubtractDS(gm2,gm1);
         if (!c1) gm1->tDelete(); gmx = gm2;
      }
      else
      {
         REPORT
         // what if New throws and exception
         Try { gmx = mtd.New(nr,nc,this); }
         CatchAll
         {
            if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
#ifdef TEMPS_DESTROYED_QUICKLY
            delete this;
#endif
            ReThrow;
         }
         SubtractDS(gmx,gm1,gm2);
         if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
         gmx->ReleaseAndDelete();
      }
   }
#ifdef TEMPS_DESTROYED_QUICKLY
   delete this;
#endif
   return gmx;
}

GeneralMatrix* SPMatrix::Evaluate(MatrixType mtd)
{
   REPORT
   Tracer tr("SPMatrix::Evaluate");
   gm1=((BaseMatrix*&)bm1)->Evaluate(); gm2=((BaseMatrix*&)bm2)->Evaluate();
   int nr=gm1->Nrows(); int nc=gm1->Ncols();
   if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
   {
      Try { Throw(IncompatibleDimensionsException(*gm1, *gm2)); }
      CatchAll
      {
         gm1->tDelete(); gm2->tDelete();
#ifdef TEMPS_DESTROYED_QUICKLY
         delete this;
#endif
         ReThrow;
      }
   }
   MatrixType mt1 = gm1->Type(), mt2 = gm2->Type();
   MatrixType mts = mt1.SP(mt2);
   if (!mtd) { REPORT mtd = mts; }
   else if (!(mtd.DataLossOK || mtd >= mts))