newmat5.cpp

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

/// \ingroup newmat
///@{

/// \file newmat5.cpp
/// Transpose, evaluate, operations with scalar, matrix input.

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

//#define WANT_STREAM

#include "include.h"

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

#ifdef use_namespace
namespace NEWMAT {
#endif


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


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


GeneralMatrix* GeneralMatrix::Transpose(TransposedMatrix* tm, MatrixType mt)
{
   GeneralMatrix* gm1;

   if (Compare(Type().t(),mt))
   {
      REPORT
      gm1 = mt.New(ncols_val,nrows_val,tm);
      for (int i=0; i<ncols_val; i++)
      {
         MatrixRow mr(gm1, StoreOnExit+DirectPart, i);
         MatrixCol mc(this, mr.Data(), LoadOnEntry, i);
      }
   }
   else
   {
      REPORT
      gm1 = mt.New(ncols_val,nrows_val,tm);
      MatrixRow mr(this, LoadOnEntry);
      MatrixCol mc(gm1, StoreOnExit+DirectPart);
      int i = nrows_val;
      while (i--) { mc.Copy(mr); mr.Next(); mc.Next(); }
   }
   tDelete(); gm1->ReleaseAndDelete(); return gm1;
}

GeneralMatrix* SymmetricMatrix::Transpose(TransposedMatrix*, MatrixType mt)
{ REPORT  return Evaluate(mt); }


GeneralMatrix* DiagonalMatrix::Transpose(TransposedMatrix*, MatrixType mt)
{ REPORT return Evaluate(mt); }

GeneralMatrix* ColumnVector::Transpose(TransposedMatrix*, MatrixType mt)
{
   REPORT
   GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
   gmx->nrows_val = 1; gmx->ncols_val = gmx->storage = storage;
   return BorrowStore(gmx,mt);
}

GeneralMatrix* RowVector::Transpose(TransposedMatrix*, MatrixType mt)
{
   REPORT
   GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
   gmx->ncols_val = 1; gmx->nrows_val = gmx->storage = storage;
   return BorrowStore(gmx,mt);
}

GeneralMatrix* IdentityMatrix::Transpose(TransposedMatrix*, MatrixType mt)
{ REPORT return Evaluate(mt); }

GeneralMatrix* GeneralMatrix::Evaluate(MatrixType mt)
{
   if (Compare(this->Type(),mt)) { REPORT return this; }
   REPORT
   GeneralMatrix* gmx = mt.New(nrows_val,ncols_val,this);
   MatrixRow mr(this, LoadOnEntry);
   MatrixRow mrx(gmx, StoreOnExit+DirectPart);
   int i=nrows_val;
   while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
   tDelete(); gmx->ReleaseAndDelete(); return gmx;
}

GeneralMatrix* CroutMatrix::Evaluate(MatrixType mt)
{
   if (Compare(this->Type(),mt)) { REPORT return this; }
   REPORT
   Tracer et("CroutMatrix::Evaluate");
   bool dummy = true;
   if (dummy) Throw(ProgramException("Illegal use of CroutMatrix", *this));
   return this;
}

GeneralMatrix* GenericMatrix::Evaluate(MatrixType mt)
   { REPORT  return gm->Evaluate(mt); }

GeneralMatrix* ShiftedMatrix::Evaluate(MatrixType mt)
{
   gm=((BaseMatrix*&)bm)->Evaluate();
   int nr=gm->Nrows(); int nc=gm->Ncols();
   Compare(gm->Type().AddEqualEl(),mt);
   if (!(mt==gm->Type()))
   {
      REPORT
      GeneralMatrix* gmx = mt.New(nr,nc,this);
      MatrixRow mr(gm, LoadOnEntry);
      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
      while (nr--) { mrx.Add(mr,f); mrx.Next(); mr.Next(); }
      gmx->ReleaseAndDelete(); gm->tDelete();
      return gmx;
   }
   else if (gm->reuse())
   {
      REPORT gm->Add(f);
      return gm;
   }
   else
   {
      REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc,this);
      gmy->ReleaseAndDelete(); gmy->Add(gm,f);
      return gmy;
   }
}

GeneralMatrix* NegShiftedMatrix::Evaluate(MatrixType mt)
{
   gm=((BaseMatrix*&)bm)->Evaluate();
   int nr=gm->Nrows(); int nc=gm->Ncols();
   Compare(gm->Type().AddEqualEl(),mt);
   if (!(mt==gm->Type()))
   {
      REPORT
      GeneralMatrix* gmx = mt.New(nr,nc,this);
      MatrixRow mr(gm, LoadOnEntry);
      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
      while (nr--) { mrx.NegAdd(mr,f); mrx.Next(); mr.Next(); }
      gmx->ReleaseAndDelete(); gm->tDelete();
      return gmx;
   }
   else if (gm->reuse())
   {
      REPORT gm->NegAdd(f);
      return gm;
   }
   else
   {
      REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc,this);
      gmy->ReleaseAndDelete(); gmy->NegAdd(gm,f);
      return gmy;
   }
}

GeneralMatrix* ScaledMatrix::Evaluate(MatrixType mt)
{
   gm=((BaseMatrix*&)bm)->Evaluate();
   int nr=gm->Nrows(); int nc=gm->Ncols();
   if (Compare(gm->Type(),mt))
   {
      if (gm->reuse())
      {
         REPORT gm->Multiply(f);
         return gm;
      }
      else
      {
         REPORT GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
         gmx->ReleaseAndDelete(); gmx->Multiply(gm,f);
         return gmx;
      }
   }
   else
   {
      REPORT
      GeneralMatrix* gmx = mt.New(nr,nc,this);
      MatrixRow mr(gm, LoadOnEntry);
      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
      while (nr--) { mrx.Multiply(mr,f); mrx.Next(); mr.Next(); }
      gmx->ReleaseAndDelete(); gm->tDelete();
      return gmx;
   }
}

GeneralMatrix* NegatedMatrix::Evaluate(MatrixType mt)
{
   gm=((BaseMatrix*&)bm)->Evaluate();
   int nr=gm->Nrows(); int nc=gm->Ncols();
   if (Compare(gm->Type(),mt))
   {
      if (gm->reuse())
      {
         REPORT gm->Negate();
         return gm;
      }
      else
      {
         REPORT
         GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
         gmx->ReleaseAndDelete(); gmx->Negate(gm);
         return gmx;
      }
   }
   else
   {
      REPORT
      GeneralMatrix* gmx = mt.New(nr,nc,this);
      MatrixRow mr(gm, LoadOnEntry);
      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
      while (nr--) { mrx.Negate(mr); mrx.Next(); mr.Next(); }
      gmx->ReleaseAndDelete(); gm->tDelete();
      return gmx;
   }
}

GeneralMatrix* ReversedMatrix::Evaluate(MatrixType mt)
{
   gm=((BaseMatrix*&)bm)->Evaluate(); GeneralMatrix* gmx;

   if ((gm->Type()).is_band() && ! (gm->Type()).is_diagonal())
   {
      gm->tDelete();
      Throw(NotDefinedException("Reverse", "band matrices"));
   }

   if (gm->reuse()) { REPORT gm->ReverseElements(); gmx = gm; }
   else
   {
      REPORT
      gmx = gm->Type().New(gm->Nrows(), gm->Ncols(), this);
      gmx->ReverseElements(gm); gmx->ReleaseAndDelete();
   }
   return gmx->Evaluate(mt); // target matrix is different type?

}

GeneralMatrix* TransposedMatrix::Evaluate(MatrixType mt)
{
   REPORT
   gm=((BaseMatrix*&)bm)->Evaluate();
   Compare(gm->Type().t(),mt);
   GeneralMatrix* gmx=gm->Transpose(this, mt);
   return gmx;
}

GeneralMatrix* RowedMatrix::Evaluate(MatrixType mt)
{
   gm = ((BaseMatrix*&)bm)->Evaluate();
   GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
   gmx->nrows_val = 1; gmx->ncols_val = gmx->storage = gm->storage;
   return gm->BorrowStore(gmx,mt);
}

GeneralMatrix* ColedMatrix::Evaluate(MatrixType mt)
{
   gm = ((BaseMatrix*&)bm)->Evaluate();
   GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
   gmx->ncols_val = 1; gmx->nrows_val = gmx->storage = gm->storage;
   return gm->BorrowStore(gmx,mt);
}

GeneralMatrix* DiagedMatrix::Evaluate(MatrixType mt)
{
   gm = ((BaseMatrix*&)bm)->Evaluate();
   GeneralMatrix* gmx = new DiagonalMatrix; MatrixErrorNoSpace(gmx);
   gmx->nrows_val = gmx->ncols_val = gmx->storage = gm->storage;
   return gm->BorrowStore(gmx,mt);
}

GeneralMatrix* MatedMatrix::Evaluate(MatrixType mt)