newmat6.cpp

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

/// \ingroup newmat
///@{

/// \file newmat6.cpp
/// Operators, element access.

// 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__,6); ++ExeCount; }
#else
#define REPORT {}
#endif

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

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


/****************************** operators *******************************/

Real& Matrix::operator()(int m, int n)
{
   REPORT
   if (m<=0 || m>nrows_val || n<=0 || n>ncols_val)
      Throw(IndexException(m,n,*this));
   return store[(m-1)*ncols_val+n-1];
}

Real& SymmetricMatrix::operator()(int m, int n)
{
   REPORT
   if (m<=0 || n<=0 || m>nrows_val || n>ncols_val)
      Throw(IndexException(m,n,*this));
   if (m>=n) return store[tristore(m-1)+n-1];
   else return store[tristore(n-1)+m-1];
}

Real& UpperTriangularMatrix::operator()(int m, int n)
{
   REPORT
   if (m<=0 || n<m || n>ncols_val)
      Throw(IndexException(m,n,*this));
   return store[(m-1)*ncols_val+n-1-tristore(m-1)];
}

Real& LowerTriangularMatrix::operator()(int m, int n)
{
   REPORT
   if (n<=0 || m<n || m>nrows_val)
      Throw(IndexException(m,n,*this));
   return store[tristore(m-1)+n-1];
}

Real& DiagonalMatrix::operator()(int m, int n)
{
   REPORT
   if (n<=0 || m!=n || m>nrows_val || n>ncols_val)
      Throw(IndexException(m,n,*this));
   return store[n-1];
}

Real& DiagonalMatrix::operator()(int m)
{
   REPORT
   if (m<=0 || m>nrows_val) Throw(IndexException(m,*this));
   return store[m-1];
}

Real& ColumnVector::operator()(int m)
{
   REPORT
   if (m<=0 || m> nrows_val) Throw(IndexException(m,*this));
   return store[m-1];
}

Real& RowVector::operator()(int n)
{
   REPORT
   if (n<=0 || n> ncols_val) Throw(IndexException(n,*this));
   return store[n-1];
}

Real& BandMatrix::operator()(int m, int n)
{
   REPORT
   int w = upper_val+lower_val+1; int i = lower_val+n-m;
   if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w)
      Throw(IndexException(m,n,*this));
   return store[w*(m-1)+i];
}

Real& UpperBandMatrix::operator()(int m, int n)
{
   REPORT
   int w = upper_val+1; int i = n-m;
   if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w)
      Throw(IndexException(m,n,*this));
   return store[w*(m-1)+i];
}

Real& LowerBandMatrix::operator()(int m, int n)
{
   REPORT
   int w = lower_val+1; int i = lower_val+n-m;
   if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w)
      Throw(IndexException(m,n,*this));
   return store[w*(m-1)+i];
}

Real& SymmetricBandMatrix::operator()(int m, int n)
{
   REPORT
   int w = lower_val+1;
   if (m>=n)
   {
      REPORT
      int i = lower_val+n-m;
      if ( m>nrows_val || n<=0 || i<0 )
         Throw(IndexException(m,n,*this));
      return store[w*(m-1)+i];
   }
   else
   {
      REPORT
      int i = lower_val+m-n;
      if ( n>nrows_val || m<=0 || i<0 )
         Throw(IndexException(m,n,*this));
      return store[w*(n-1)+i];
   }
}


Real Matrix::operator()(int m, int n) const
{
   REPORT
   if (m<=0 || m>nrows_val || n<=0 || n>ncols_val)
      Throw(IndexException(m,n,*this));
   return store[(m-1)*ncols_val+n-1];
}

Real SymmetricMatrix::operator()(int m, int n) const
{
   REPORT
   if (m<=0 || n<=0 || m>nrows_val || n>ncols_val)
      Throw(IndexException(m,n,*this));
   if (m>=n) return store[tristore(m-1)+n-1];
   else return store[tristore(n-1)+m-1];
}

Real UpperTriangularMatrix::operator()(int m, int n) const
{
   REPORT
   if (m<=0 || n<m || n>ncols_val)
      Throw(IndexException(m,n,*this));
   return store[(m-1)*ncols_val+n-1-tristore(m-1)];
}

Real LowerTriangularMatrix::operator()(int m, int n) const
{
   REPORT
   if (n<=0 || m<n || m>nrows_val)
      Throw(IndexException(m,n,*this));
   return store[tristore(m-1)+n-1];
}

Real DiagonalMatrix::operator()(int m, int n) const
{
   REPORT
   if (n<=0 || m!=n || m>nrows_val || n>ncols_val)
      Throw(IndexException(m,n,*this));
   return store[n-1];
}

Real DiagonalMatrix::operator()(int m) const
{
   REPORT
   if (m<=0 || m>nrows_val) Throw(IndexException(m,*this));
   return store[m-1];
}

Real ColumnVector::operator()(int m) const
{
   REPORT
   if (m<=0 || m> nrows_val) Throw(IndexException(m,*this));
   return store[m-1];
}

Real RowVector::operator()(int n) const
{
   REPORT
   if (n<=0 || n> ncols_val) Throw(IndexException(n,*this));
   return store[n-1];
}

Real BandMatrix::operator()(int m, int n) const
{
   REPORT
   int w = upper_val+lower_val+1; int i = lower_val+n-m;
   if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w)
      Throw(IndexException(m,n,*this));
   return store[w*(m-1)+i];
}

Real UpperBandMatrix::operator()(int m, int n) const
{
   REPORT
   int w = upper_val+1; int i = n-m;
   if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w)
      Throw(IndexException(m,n,*this));
   return store[w*(m-1)+i];
}

Real LowerBandMatrix::operator()(int m, int n) const
{
   REPORT
   int w = lower_val+1; int i = lower_val+n-m;
   if (m<=0 || m>nrows_val || n<=0 || n>ncols_val || i<0 || i>=w)
      Throw(IndexException(m,n,*this));
   return store[w*(m-1)+i];
}

Real SymmetricBandMatrix::operator()(int m, int n) const
{
   REPORT
   int w = lower_val+1;
   if (m>=n)
   {
      REPORT
      int i = lower_val+n-m;
      if ( m>nrows_val || n<=0 || i<0 )
         Throw(IndexException(m,n,*this));
      return store[w*(m-1)+i];
   }
   else
   {
      REPORT
      int i = lower_val+m-n;
      if ( n>nrows_val || m<=0 || i<0 )
         Throw(IndexException(m,n,*this));
      return store[w*(n-1)+i];
   }
}


Real BaseMatrix::as_scalar() const
{
   REPORT
   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();

   if (gm->nrows_val!=1 || gm->ncols_val!=1)
   {
      Tracer tr("as_scalar");
      Try
         { Throw(ProgramException("Cannot convert to scalar", *gm)); }
      CatchAll { gm->tDelete(); ReThrow; }
   }

   Real x = *(gm->store); gm->tDelete(); return x;
}


AddedMatrix BaseMatrix::operator+(const BaseMatrix& bm) const
{ REPORT return AddedMatrix(this, &bm); }

SPMatrix SP(const BaseMatrix& bm1,const BaseMatrix& bm2)
{ REPORT return SPMatrix(&bm1, &bm2); }

KPMatrix KP(const BaseMatrix& bm1,const BaseMatrix& bm2)
{ REPORT return KPMatrix(&bm1, &bm2); }

MultipliedMatrix BaseMatrix::operator*(const BaseMatrix& bm) const
{ REPORT return MultipliedMatrix(this, &bm); }

ConcatenatedMatrix BaseMatrix::operator|(const BaseMatrix& bm) const
{ REPORT return ConcatenatedMatrix(this, &bm); }

StackedMatrix BaseMatrix::operator&(const BaseMatrix& bm) const
{ REPORT return StackedMatrix(this, &bm); }

SolvedMatrix InvertedMatrix::operator*(const BaseMatrix& bmx) const
{ REPORT return SolvedMatrix(bm, &bmx); }

SubtractedMatrix BaseMatrix::operator-(const BaseMatrix& bm) const
{ REPORT return SubtractedMatrix(this, &bm); }

ShiftedMatrix BaseMatrix::operator+(Real f) const
{ REPORT return ShiftedMatrix(this, f); }

ShiftedMatrix operator+(Real f, const BaseMatrix& BM)
{ REPORT return ShiftedMatrix(&BM, f); }

NegShiftedMatrix operator-(Real f, const BaseMatrix& bm)
{ REPORT return NegShiftedMatrix(f, &bm); }

ScaledMatrix BaseMatrix::operator*(Real f) const
{ REPORT return ScaledMatrix(this, f); }

ScaledMatrix BaseMatrix::operator/(Real f) const
{ REPORT return ScaledMatrix(this, 1.0/f); }

ScaledMatrix operator*(Real f, const BaseMatrix& BM)
{ REPORT return ScaledMatrix(&BM, f); }

ShiftedMatrix BaseMatrix::operator-(Real f) const
{ REPORT return ShiftedMatrix(this, -f); }

TransposedMatrix BaseMatrix::t() const
{ REPORT return TransposedMatrix(this); }

NegatedMatrix BaseMatrix::operator-() const
{ REPORT return NegatedMatrix(this); }

ReversedMatrix BaseMatrix::reverse() const
{ REPORT return ReversedMatrix(this); }

InvertedMatrix BaseMatrix::i() const
{ REPORT return InvertedMatrix(this); }


RowedMatrix BaseMatrix::as_row() const
{ REPORT return RowedMatrix(this); }

ColedMatrix BaseMatrix::as_column() const
{ REPORT return ColedMatrix(this); }

DiagedMatrix BaseMatrix::as_diagonal() const
{ REPORT return DiagedMatrix(this); }

MatedMatrix BaseMatrix::as_matrix(int nrx, int ncx) const
{ REPORT return MatedMatrix(this,nrx,ncx); }


void GeneralMatrix::operator=(Real f)
{ REPORT int i=storage; Real* s=store; while (i--) { *s++ = f; } }

void Matrix::operator=(const BaseMatrix& X)
{
   REPORT //CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::Rt);
} 

void SquareMatrix::operator=(const BaseMatrix& X)
{
   REPORT //CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::Rt);
   if (nrows_val != ncols_val)
      { Tracer tr("SquareMatrix(=)"); Throw(NotSquareException(*this)); }
}

void SquareMatrix::operator=(const Matrix& m)
{
   REPORT
   if (m.nrows_val != m.ncols_val)
      { Tracer tr("SquareMatrix(=Matrix)"); Throw(NotSquareException(*this)); }
   Eq(m);
}

void RowVector::operator=(const BaseMatrix& X)
{
   REPORT  // CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::RV);
   if (nrows_val!=1)
      { Tracer tr("RowVector(=)"); Throw(VectorException(*this)); }
}

void ColumnVector::operator=(const BaseMatrix& X)
{
   REPORT //CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::CV);
   if (ncols_val!=1)
      { Tracer tr("ColumnVector(=)"); Throw(VectorException(*this)); }
}

void SymmetricMatrix::operator=(const BaseMatrix& X)
{
   REPORT // CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::Sm);
}

void UpperTriangularMatrix::operator=(const BaseMatrix& X)
{
   REPORT //CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::UT);
}

void LowerTriangularMatrix::operator=(const BaseMatrix& X)
{
   REPORT //CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::LT);
}

void DiagonalMatrix::operator=(const BaseMatrix& X)
{
   REPORT // CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::Dg);
}

void IdentityMatrix::operator=(const BaseMatrix& X)
{
   REPORT // CheckConversion(X);
   // MatrixConversionCheck mcc;
   Eq(X,MatrixType::Id);
}


void CroutMatrix::operator=(const CroutMatrix& gm)
{
   if (&gm == this) { REPORT tag_val = -1; return; }
   REPORT
   if (indx > 0) { delete [] indx; indx = 0; }
   ((CroutMatrix&)gm).get_aux(*this);
   Eq(gm);
}
   




void GeneralMatrix::operator<<(const double* r)
{
   REPORT
   int i = storage; Real* s=store;
   while(i--) *s++ = (Real)*r++;
}


void GeneralMatrix::operator<<(const float* r)
{
   REPORT
   int i = storage; Real* s=store;
   while(i--) *s++ = (Real)*r++;
}