newmat4.cpp

C++矩阵算法库

//$$ newmat4.cpp       Constructors, ReSize, basic utilities

// Copyright (C) 1991,2,3,4,8,9: 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__,4); ++ExeCount; }
#else
#define REPORT {}
#endif


#define DO_SEARCH                   // search for LHS of = in RHS

// ************************* general utilities *************************/

static int tristore(int n)                    // elements in triangular matrix
{ return (n*(n+1))/2; }


// **************************** constructors ***************************/

GeneralMatrix::GeneralMatrix()
{ store=0; storage=0; nrows=0; ncols=0; tag=-1; }

GeneralMatrix::GeneralMatrix(ArrayLengthSpecifier s)
{
   REPORT
   storage=s.Value(); tag=-1;
   if (storage)
   {
      store = new Real [storage]; MatrixErrorNoSpace(store);
      MONITOR_REAL_NEW("Make (GenMatrix)",storage,store)
   }
   else store = 0;
}

Matrix::Matrix(int m, int n) : GeneralMatrix(m*n)
{ REPORT nrows=m; ncols=n; }

SymmetricMatrix::SymmetricMatrix(ArrayLengthSpecifier n)
   : GeneralMatrix(tristore(n.Value()))
{ REPORT nrows=n.Value(); ncols=n.Value(); }

UpperTriangularMatrix::UpperTriangularMatrix(ArrayLengthSpecifier n)
   : GeneralMatrix(tristore(n.Value()))
{ REPORT nrows=n.Value(); ncols=n.Value(); }

LowerTriangularMatrix::LowerTriangularMatrix(ArrayLengthSpecifier n)
   : GeneralMatrix(tristore(n.Value()))
{ REPORT nrows=n.Value(); ncols=n.Value(); }

DiagonalMatrix::DiagonalMatrix(ArrayLengthSpecifier m) : GeneralMatrix(m)
{ REPORT nrows=m.Value(); ncols=m.Value(); }

Matrix::Matrix(const BaseMatrix& M)
{
   REPORT // CheckConversion(M);
   // MatrixConversionCheck mcc;
   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Rt);
   GetMatrix(gmx);
}

RowVector::RowVector(const BaseMatrix& M) : Matrix(M)
{
   if (nrows!=1)
   {
      Tracer tr("RowVector");
      Throw(VectorException(*this));
   }
}

ColumnVector::ColumnVector(const BaseMatrix& M) : Matrix(M)
{
   if (ncols!=1)
   {
      Tracer tr("ColumnVector");
      Throw(VectorException(*this));
   }
}

SymmetricMatrix::SymmetricMatrix(const BaseMatrix& M)
{
   REPORT  // CheckConversion(M);
   // MatrixConversionCheck mcc;
   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Sm);
   GetMatrix(gmx);
}

UpperTriangularMatrix::UpperTriangularMatrix(const BaseMatrix& M)
{
   REPORT // CheckConversion(M);
   // MatrixConversionCheck mcc;
   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::UT);
   GetMatrix(gmx);
}

LowerTriangularMatrix::LowerTriangularMatrix(const BaseMatrix& M)
{
   REPORT // CheckConversion(M);
   // MatrixConversionCheck mcc;
   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::LT);
   GetMatrix(gmx);
}

DiagonalMatrix::DiagonalMatrix(const BaseMatrix& M)
{
   REPORT //CheckConversion(M);
   // MatrixConversionCheck mcc;
   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Dg);
   GetMatrix(gmx);
}

IdentityMatrix::IdentityMatrix(const BaseMatrix& M)
{
   REPORT //CheckConversion(M);
   // MatrixConversionCheck mcc;
   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Id);
   GetMatrix(gmx);
}

GeneralMatrix::~GeneralMatrix()
{
   if (store)
   {
      MONITOR_REAL_DELETE("Free (GenMatrix)",storage,store)
      delete [] store;
   }
}

CroutMatrix::CroutMatrix(const BaseMatrix& m)
{
   REPORT
   Tracer tr("CroutMatrix");
   indx = 0;                     // in case of exception at next line
   GeneralMatrix* gm = ((BaseMatrix&)m).Evaluate(MatrixType::Rt);
   GetMatrix(gm);
   if (nrows!=ncols) { CleanUp(); Throw(NotSquareException(*gm)); }
   d=true; sing=false;
   indx=new int [nrows]; MatrixErrorNoSpace(indx);
   MONITOR_INT_NEW("Index (CroutMat)",nrows,indx)
   ludcmp();
}

CroutMatrix::~CroutMatrix()
{
   MONITOR_INT_DELETE("Index (CroutMat)",nrows,indx)
   delete [] indx;
}

//ReturnMatrixX::ReturnMatrixX(GeneralMatrix& gmx)
//{
//   REPORT
//   gm = gmx.Image(); gm->ReleaseAndDelete();
//}

#ifndef TEMPS_DESTROYED_QUICKLY_R

GeneralMatrix::operator ReturnMatrixX() const
{
   REPORT
   GeneralMatrix* gm = Image(); gm->ReleaseAndDelete();
   return ReturnMatrixX(gm);
}

#else

GeneralMatrix::operator ReturnMatrixX&() const
{
   REPORT
   GeneralMatrix* gm = Image(); gm->ReleaseAndDelete();
   ReturnMatrixX* x = new ReturnMatrixX(gm);
   MatrixErrorNoSpace(x); return *x;
}

#endif

#ifndef TEMPS_DESTROYED_QUICKLY_R

ReturnMatrixX GeneralMatrix::ForReturn() const
{
   REPORT
   GeneralMatrix* gm = Image(); gm->ReleaseAndDelete();
   return ReturnMatrixX(gm);
}

#else

ReturnMatrixX& GeneralMatrix::ForReturn() const
{
   REPORT
   GeneralMatrix* gm = Image(); gm->ReleaseAndDelete();
   ReturnMatrixX* x = new ReturnMatrixX(gm);
   MatrixErrorNoSpace(x); return *x;
}

#endif

// ************************** ReSize matrices ***************************/

void GeneralMatrix::ReSize(int nr, int nc, int s)
{
   REPORT
   if (store)
   {
      MONITOR_REAL_DELETE("Free (ReDimensi)",storage,store)
      delete [] store;
   }
   storage=s; nrows=nr; ncols=nc; tag=-1;
   if (s)
   {
      store = new Real [storage]; MatrixErrorNoSpace(store);
      MONITOR_REAL_NEW("Make (ReDimensi)",storage,store)
   }
   else store = 0;
}

void Matrix::ReSize(int nr, int nc)
{ REPORT GeneralMatrix::ReSize(nr,nc,nr*nc); }

void SymmetricMatrix::ReSize(int nr)
{ REPORT GeneralMatrix::ReSize(nr,nr,tristore(nr)); }

void UpperTriangularMatrix::ReSize(int nr)
{ REPORT GeneralMatrix::ReSize(nr,nr,tristore(nr)); }

void LowerTriangularMatrix::ReSize(int nr)
{ REPORT GeneralMatrix::ReSize(nr,nr,tristore(nr)); }

void DiagonalMatrix::ReSize(int nr)
{ REPORT GeneralMatrix::ReSize(nr,nr,nr); }

void RowVector::ReSize(int nc)
{ REPORT GeneralMatrix::ReSize(1,nc,nc); }

void ColumnVector::ReSize(int nr)
{ REPORT GeneralMatrix::ReSize(nr,1,nr); }

void RowVector::ReSize(int nr, int nc)
{
   Tracer tr("RowVector::ReSize");
   if (nr != 1) Throw(VectorException(*this));
   REPORT GeneralMatrix::ReSize(1,nc,nc);
}

void ColumnVector::ReSize(int nr, int nc)
{
   Tracer tr("ColumnVector::ReSize");
   if (nc != 1) Throw(VectorException(*this));
   REPORT GeneralMatrix::ReSize(nr,1,nr);
}

void IdentityMatrix::ReSize(int nr)
{ REPORT GeneralMatrix::ReSize(nr,nr,1); *store = 1; }


void Matrix::ReSize(const GeneralMatrix& A)
{ REPORT  ReSize(A.Nrows(), A.Ncols()); }

void nricMatrix::ReSize(const GeneralMatrix& A)
{ REPORT  ReSize(A.Nrows(), A.Ncols()); }

void ColumnVector::ReSize(const GeneralMatrix& A)
{ REPORT  ReSize(A.Nrows(), A.Ncols()); }

void RowVector::ReSize(const GeneralMatrix& A)
{ REPORT  ReSize(A.Nrows(), A.Ncols()); }

void SymmetricMatrix::ReSize(const GeneralMatrix& A)
{
   REPORT
   int n = A.Nrows();
   if (n != A.Ncols())
   {
      Tracer tr("SymmetricMatrix::ReSize(GM)");
      Throw(NotSquareException(*this));
   }
   ReSize(n);
}

void DiagonalMatrix::ReSize(const GeneralMatrix& A)
{
   REPORT
   int n = A.Nrows();
   if (n != A.Ncols())
   {
      Tracer tr("DiagonalMatrix::ReSize(GM)");
      Throw(NotSquareException(*this));
   }
   ReSize(n);
}

void UpperTriangularMatrix::ReSize(const GeneralMatrix& A)
{
   REPORT
   int n = A.Nrows();
   if (n != A.Ncols())
   {
      Tracer tr("UpperTriangularMatrix::ReSize(GM)");
      Throw(NotSquareException(*this));
   }
   ReSize(n);
}

void LowerTriangularMatrix::ReSize(const GeneralMatrix& A)
{
   REPORT
   int n = A.Nrows();
   if (n != A.Ncols())
   {
      Tracer tr("LowerTriangularMatrix::ReSize(GM)");
      Throw(NotSquareException(*this));
   }
   ReSize(n);
}

void IdentityMatrix::ReSize(const GeneralMatrix& A)
{
   REPORT
   int n = A.Nrows();
   if (n != A.Ncols())
   {
      Tracer tr("IdentityMatrix::ReSize(GM)");
      Throw(NotSquareException(*this));
   }
   ReSize(n);
}

void GeneralMatrix::ReSize(const GeneralMatrix&)
{
   REPORT
   Tracer tr("GeneralMatrix::ReSize(GM)");
   Throw(NotDefinedException("ReSize", "this type of matrix"));
}

void GeneralMatrix::ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix&)
{ REPORT ReSize(A); }

void GeneralMatrix::ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix&)
{ REPORT ReSize(A); }


// ************************* SameStorageType ******************************/

// SameStorageType checks A and B have same storage type including bandwidth
// It does not check same dimensions since we assume this is already done

bool GeneralMatrix::SameStorageType(const GeneralMatrix& A) const
{
   REPORT
   return Type() == A.Type();
}


// ******************* manipulate types, storage **************************/

int GeneralMatrix::search(const BaseMatrix* s) const
{ REPORT return (s==this) ? 1 : 0; }

int GenericMatrix::search(const BaseMatrix* s) const
{ REPORT return gm->search(s); }

int MultipliedMatrix::search(const BaseMatrix* s) const
{ REPORT return bm1->search(s) + bm2->search(s); }

int ShiftedMatrix::search(const BaseMatrix* s) const
{ REPORT return bm->search(s); }

int NegatedMatrix::search(const BaseMatrix* s) const
{ REPORT return bm->search(s); }

int ReturnMatrixX::search(const BaseMatrix* s) const
{ REPORT return (s==gm) ? 1 : 0; }

MatrixType Matrix::Type() const { return MatrixType::Rt; }
MatrixType SymmetricMatrix::Type() const { return MatrixType::Sm; }
MatrixType UpperTriangularMatrix::Type() const { return MatrixType::UT; }
MatrixType LowerTriangularMatrix::Type() const { return MatrixType::LT; }
MatrixType DiagonalMatrix::Type() const { return MatrixType::Dg; }
MatrixType RowVector::Type() const { return MatrixType::RV; }
MatrixType ColumnVector::Type() const { return MatrixType::CV; }
MatrixType CroutMatrix::Type() const { return MatrixType::Ct; }
MatrixType BandMatrix::Type() const { return MatrixType::BM; }
MatrixType UpperBandMatrix::Type() const { return MatrixType::UB; }
MatrixType LowerBandMatrix::Type() const { return MatrixType::LB; }
MatrixType SymmetricBandMatrix::Type() const { return MatrixType::SB; }

MatrixType IdentityMatrix::Type() const { return MatrixType::Id; }



MatrixBandWidth BaseMatrix::BandWidth() const { REPORT return -1; }
MatrixBandWidth DiagonalMatrix::BandWidth() const { REPORT return 0; }
MatrixBandWidth IdentityMatrix::BandWidth() const { REPORT return 0; }

MatrixBandWidth UpperTriangularMatrix::BandWidth() const
   { REPORT return MatrixBandWidth(0,-1); }

MatrixBandWidth LowerTriangularMatrix::BandWidth() const
   { REPORT return MatrixBandWidth(-1,0); }

MatrixBandWidth BandMatrix::BandWidth() const
   { REPORT return MatrixBandWidth(lower,upper); }

MatrixBandWidth GenericMatrix::BandWidth()const
   { REPORT return gm->BandWidth(); }

MatrixBandWidth AddedMatrix::BandWidth() const
   { REPORT return gm1->BandWidth() + gm2->BandWidth(); }

MatrixBandWidth SPMatrix::BandWidth() const
   { REPORT return gm1->BandWidth().minimum(gm2->BandWidth()); }

MatrixBandWidth KPMatrix::BandWidth() const
{
   int lower, upper;
   MatrixBandWidth bw1 = gm1->BandWidth(), bw2 = gm2->BandWidth();
   if (bw1.Lower() < 0)
   {
      if (bw2.Lower() < 0) { REPORT lower = -1; }
      else { REPORT lower = bw2.Lower() + (gm1->Nrows() - 1) * gm2->Nrows(); }
   }
   else
   {
      if (bw2.Lower() < 0)
         { REPORT lower = (1 + bw1.Lower()) * gm2->Nrows() - 1; }
      else { REPORT lower = bw2.Lower() + bw1.Lower() * gm2->Nrows(); }
   }
   if (bw1.Upper() < 0)
   {
      if (bw2.Upper() < 0) { REPORT upper = -1; }
      else { REPORT upper = bw2.Upper() + (gm1->Nrows() - 1) * gm2->Nrows(); }
   }
   else
   {
      if (bw2.Upper() < 0)
         { REPORT upper = (1 + bw1.Upper()) * gm2->Nrows() - 1; }
      else { REPORT upper = bw2.Upper() + bw1.Upper() * gm2->Nrows(); }
   }
   return MatrixBandWidth(lower, upper);
}

MatrixBandWidth MultipliedMatrix::BandWidth() const
{ REPORT return gm1->BandWidth() * gm2->BandWidth(); }

MatrixBandWidth ConcatenatedMatrix::BandWidth() const { REPORT return -1; }

MatrixBandWidth SolvedMatrix::BandWidth() const
{
   if (+gm1->Type() & MatrixType::Diagonal)
      { REPORT return gm2->BandWidth(); }
   else { REPORT return -1; }
}

MatrixBandWidth ScaledMatrix::BandWidth() const
   { REPORT return gm->BandWidth(); }

MatrixBandWidth NegatedMatrix::BandWidth() const
   { REPORT return gm->BandWidth(); }

MatrixBandWidth TransposedMatrix::BandWidth() const
   { REPORT return gm->BandWidth().t(); }

MatrixBandWidth InvertedMatrix::BandWidth() const