newmat7.cpp

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

//$$ newmat7.cpp     Invert, solve, binary operations

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

#include "include.h"

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

#ifdef use_namespace
namespace NEWMAT {
#endif


#ifdef DO_REPORT
#define REPORT { static ExeCounter ExeCount(__LINE__,7); ++ExeCount; }
#else
#define REPORT {}
#endif


/***************************** solve routines ******************************/

GeneralMatrix* GeneralMatrix::MakeSolver()
{
   REPORT
   GeneralMatrix* gm = new CroutMatrix(*this);
   MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
}

GeneralMatrix* Matrix::MakeSolver()
{
   REPORT
   GeneralMatrix* gm = new CroutMatrix(*this);
   MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
}

void CroutMatrix::Solver(MatrixColX& mcout, const MatrixColX& mcin)
{
   REPORT
   int i = mcin.skip; Real* el = mcin.data-i; Real* el1 = el;
   while (i--) *el++ = 0.0;
   el += mcin.storage; i = nrows - mcin.skip - mcin.storage;
   while (i--) *el++ = 0.0;
   lubksb(el1, mcout.skip);
}


// Do we need check for entirely zero output?

void UpperTriangularMatrix::Solver(MatrixColX& mcout,
   const MatrixColX& mcin)
{
   REPORT
   int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
   while (i-- > 0) *elx++ = 0.0;
   int nr = mcin.skip+mcin.storage;
   elx = mcin.data+mcin.storage; Real* el = elx;
   int j = mcout.skip+mcout.storage-nr; int nc = ncols-nr; i = nr-mcout.skip;
   while (j-- > 0) *elx++ = 0.0;
   Real* Ael = store + (nr*(2*ncols-nr+1))/2; j = 0;
   while (i-- > 0)
   {
      elx = el; Real sum = 0.0; int jx = j++; Ael -= nc;
      while (jx--) sum += *(--Ael) * *(--elx);
      elx--; *elx = (*elx - sum) / *(--Ael);
   }
}

void LowerTriangularMatrix::Solver(MatrixColX& mcout,
   const MatrixColX& mcin)
{
   REPORT
   int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
   while (i-- > 0) *elx++ = 0.0;
   int nc = mcin.skip; i = nc+mcin.storage; elx = mcin.data+mcin.storage;
   int nr = mcout.skip+mcout.storage; int j = nr-i; i = nr-nc;
   while (j-- > 0) *elx++ = 0.0;
   Real* el = mcin.data; Real* Ael = store + (nc*(nc+1))/2; j = 0;
   while (i-- > 0)
   {
      elx = el; Real sum = 0.0; int jx = j++; Ael += nc;
      while (jx--) sum += *Ael++ * *elx++;
      *elx = (*elx - sum) / *Ael++;
   }
}

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

static GeneralMatrix*
   GeneralAdd(GeneralMatrix*,GeneralMatrix*,AddedMatrix*,MatrixType);
static GeneralMatrix*
   GeneralSub(GeneralMatrix*,GeneralMatrix*,SubtractedMatrix*,MatrixType);
static GeneralMatrix*
   GeneralMult(GeneralMatrix*,GeneralMatrix*,MultipliedMatrix*,MatrixType);
static GeneralMatrix*
   GeneralSolv(GeneralMatrix*,GeneralMatrix*,BaseMatrix*,MatrixType);
static GeneralMatrix*
   GeneralSolvI(GeneralMatrix*,BaseMatrix*,MatrixType);
static GeneralMatrix*
   GeneralSP(GeneralMatrix*,GeneralMatrix*,SPMatrix*,MatrixType);

GeneralMatrix* AddedMatrix::Evaluate(MatrixType mt)
{
   REPORT
   gm1=((BaseMatrix*&)bm1)->Evaluate();
   gm2=((BaseMatrix*&)bm2)->Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
   GeneralMatrix* gmx;
   Try { gmx = GeneralAdd(gm1,gm2,this,mt); }
   CatchAll { delete this; ReThrow; }
   delete this; return gmx;
#else
   return GeneralAdd(gm1,gm2,this,mt);
#endif   
}

GeneralMatrix* SubtractedMatrix::Evaluate(MatrixType mt)
{
   REPORT
   gm1=((BaseMatrix*&)bm1)->Evaluate();
   gm2=((BaseMatrix*&)bm2)->Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
   GeneralMatrix* gmx;
   Try { gmx = GeneralSub(gm1,gm2,this,mt); }
   CatchAll { delete this; ReThrow; }
   delete this; return gmx;
#else
   return GeneralSub(gm1,gm2,this,mt);
#endif   
}

GeneralMatrix* MultipliedMatrix::Evaluate(MatrixType mt)
{
   REPORT
   gm2 = ((BaseMatrix*&)bm2)->Evaluate();
   gm2 = gm2->Evaluate(gm2->Type().MultRHS());     // no symmetric on RHS
   gm1=((BaseMatrix*&)bm1)->Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
   GeneralMatrix* gmx;
   Try { gmx = GeneralMult(gm1, gm2, this, mt); }
   CatchAll { delete this; ReThrow; }
   delete this; return gmx;
#else
   return GeneralMult(gm1, gm2, this, mt);
#endif   
}

GeneralMatrix* SolvedMatrix::Evaluate(MatrixType mt)
{
   REPORT
   gm1=((BaseMatrix*&)bm1)->Evaluate();
   gm2=((BaseMatrix*&)bm2)->Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
   GeneralMatrix* gmx;
   Try { gmx = GeneralSolv(gm1,gm2,this,mt); }
   CatchAll { delete this; ReThrow; }
   delete this; return gmx;
#else
   return GeneralSolv(gm1,gm2,this,mt);
#endif   
}

GeneralMatrix* SPMatrix::Evaluate(MatrixType mt)
{
   REPORT
   gm1=((BaseMatrix*&)bm1)->Evaluate();
   gm2=((BaseMatrix*&)bm2)->Evaluate();
#ifdef TEMPS_DESTROYED_QUICKLY
   GeneralMatrix* gmx;
   Try { gmx = GeneralSP(gm1,gm2,this,mt); }
   CatchAll { delete this; ReThrow; }
   delete this; return gmx;
#else
   return GeneralSP(gm1,gm2,this,mt);
#endif   
}

// routines for adding or subtracting matrices of identical storage structure

static void Add(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
{
   REPORT
   Real* s1=gm1->Store(); Real* s2=gm2->Store();
   Real* s=gm->Store(); int i=gm->Storage() >> 2;
   while (i--)
   {
       *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
       *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
   }
   i=gm->Storage() & 3; while (i--) *s++ = *s1++ + *s2++;
}
   
static void Add(GeneralMatrix* gm, GeneralMatrix* gm2)
{
   REPORT
   Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
   while (i--)
   { *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; }
   i=gm->Storage() & 3; while (i--) *s++ += *s2++;
}

static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
{
   REPORT
   Real* s1=gm1->Store(); Real* s2=gm2->Store();
   Real* s=gm->Store(); int i=gm->Storage() >> 2;
   while (i--)
   {
       *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
       *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
   }
   i=gm->Storage() & 3; while (i--) *s++ = *s1++ - *s2++;
}

static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm2)
{
   REPORT
   Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
   while (i--)
   { *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; }
   i=gm->Storage() & 3; while (i--) *s++ -= *s2++;
}

static void ReverseSubtract(GeneralMatrix* gm, GeneralMatrix* gm2)
{
   REPORT
   Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
   while (i--)
   {
      *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
      *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
   }
   i=gm->Storage() & 3; while (i--) { *s = *s2++ - *s; s++; }
}

static void SP(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
{
   REPORT
   Real* s1=gm1->Store(); Real* s2=gm2->Store();
   Real* s=gm->Store(); int i=gm->Storage() >> 2;
   while (i--)
   {
       *s++ = *s1++ * *s2++; *s++ = *s1++ * *s2++;
       *s++ = *s1++ * *s2++; *s++ = *s1++ * *s2++;
   }
   i=gm->Storage() & 3; while (i--) *s++ = *s1++ * *s2++;
}
   
static void SP(GeneralMatrix* gm, GeneralMatrix* gm2)
{
   REPORT
   Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
   while (i--)
   { *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; }
   i=gm->Storage() & 3; while (i--) *s++ *= *s2++;
}

// routines for adding or subtracting matrices of different storage structure

static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
{
   int nr = gm->Nrows();
   MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
   MatrixRow mr(gm, StoreOnExit+DirectPart);
   while (nr--) { mr.Add(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
}
   
static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm2)
// Add into first argument
{
   int nr = gm->Nrows();
   MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
   MatrixRow mr2(gm2, LoadOnEntry);
   while (nr--) { mr.Add(mr2); mr.Next(); mr2.Next(); }
}

static void SubtractDS
   (GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
{
   int nr = gm->Nrows();
   MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
   MatrixRow mr(gm, StoreOnExit+DirectPart);
   while (nr--) { mr.Sub(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
}

static void SubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
{
   int nr = gm->Nrows();
   MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
   MatrixRow mr2(gm2, LoadOnEntry);
   while (nr--) { mr.Sub(mr2); mr.Next(); mr2.Next(); }
}

static void ReverseSubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
{
   int nr = gm->Nrows();
   MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
   MatrixRow mr2(gm2, LoadOnEntry);
   while (nr--) { mr.RevSub(mr2); mr2.Next(); mr.Next(); }
}

static void SPDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
{
   int nr = gm->Nrows();
   MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
   MatrixRow mr(gm, StoreOnExit+DirectPart);
   while (nr--) { mr.Multiply(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
}
   
static void SPDS(GeneralMatrix* gm, GeneralMatrix* gm2)
// SP into first argument
{
   int nr = gm->Nrows();
   MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
   MatrixRow mr2(gm2, LoadOnEntry);
   while (nr--) { mr.Multiply(mr2); mr.Next(); mr2.Next(); }
}

#ifdef __GNUG__
void AddedMatrix::SelectVersion
   (MatrixType mtx, int& c1, int& c2) const
#else
void AddedMatrix::SelectVersion
   (MatrixType mtx, bool& c1, bool& c2) const
#endif
// for determining version of add and subtract 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 + bw2;
   c1 = (mtx == gm1->Type()) && (bwx == bw1);
   c2 = (mtx == gm2->Type()) && (bwx == bw2);
}

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


static GeneralMatrix* GeneralSub(GeneralMatrix* gm1, GeneralMatrix* gm2,
   SubtractedMatrix* sm, MatrixType mtx)
{
   Tracer tr("GeneralSub");
   Compare(gm1->Type() + gm2->Type(),mtx);
   int nr=gm1->Nrows(); int nc=gm1->Ncols();
   if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
      Throw(IncompatibleDimensionsException(*gm1, *gm2));
#ifdef __GNUG__
   int c1,c2; sm->SelectVersion(mtx,c1,c2);
#else
   bool 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(*gm1, *gm2));
   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.Data();                              // 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(*gm1, *gm2));
   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.Data();                            // 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(*gm1,*gm2));