newmat4.cpp

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

/// \ingroup newmat
///@{

/// \file newmat4.cpp
/// Constructors, resize, basic utilities, SimpleIntArray.


// Copyright (C) 1991,2,3,4,8,9: 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__,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_val=0; ncols_val=0; tag_val=-1; }

GeneralMatrix::GeneralMatrix(ArrayLengthSpecifier s)
{
   REPORT
   storage=s.Value(); tag_val=-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_val=m; ncols_val=n; }

SquareMatrix::SquareMatrix(ArrayLengthSpecifier n)
   : Matrix(n.Value(),n.Value())
{ REPORT }

SymmetricMatrix::SymmetricMatrix(ArrayLengthSpecifier n)
   : GeneralMatrix(tristore(n.Value()))
{ REPORT nrows_val=n.Value(); ncols_val=n.Value(); }

UpperTriangularMatrix::UpperTriangularMatrix(ArrayLengthSpecifier n)
   : GeneralMatrix(tristore(n.Value()))
{ REPORT nrows_val=n.Value(); ncols_val=n.Value(); }

LowerTriangularMatrix::LowerTriangularMatrix(ArrayLengthSpecifier n)
   : GeneralMatrix(tristore(n.Value()))
{ REPORT nrows_val=n.Value(); ncols_val=n.Value(); }

DiagonalMatrix::DiagonalMatrix(ArrayLengthSpecifier m) : GeneralMatrix(m)
{ REPORT nrows_val=m.Value(); ncols_val=m.Value(); }

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

SquareMatrix::SquareMatrix(const BaseMatrix& M) : Matrix(M)
{
   REPORT
   if (ncols_val != nrows_val)
   {
      Tracer tr("SquareMatrix");
      Throw(NotSquareException(*this));
   }
}


SquareMatrix::SquareMatrix(const Matrix& gm)
{
   REPORT
   if (gm.ncols_val != gm.nrows_val)
   {
      Tracer tr("SquareMatrix(Matrix)");
      Throw(NotSquareException(gm));
   }
   GetMatrix(&gm);
}


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

ColumnVector::ColumnVector(const BaseMatrix& M) : Matrix(M)
{
   REPORT
   if (ncols_val!=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();
   if (gm->nrows_val!=gm->ncols_val)
      { gm->tDelete(); Throw(NotSquareException(*gm)); }
   if (gm->type() == MatrixType::Ct)
      { REPORT  ((CroutMatrix*)gm)->get_aux(*this); GetMatrix(gm); }
   else
   {
      REPORT
      GeneralMatrix* gm1 = gm->Evaluate(MatrixType::Rt);
      GetMatrix(gm1);
      d=true; sing=false;
      indx=new int [nrows_val]; MatrixErrorNoSpace(indx);
      MONITOR_INT_NEW("Index (CroutMat)",nrows_val,indx)
      ludcmp();
   }
}

// could we use SetParameters instead of this
void CroutMatrix::get_aux(CroutMatrix& X)
{
   X.d = d; X.sing = sing;
   if (tag_val == 0 || tag_val == 1) // reuse the array 
      { REPORT  X.indx = indx; indx = 0; d = true; sing = true; return; }
   else if (nrows_val == 0)
      { REPORT indx = 0; d = true; sing = true; return; }
   else                              // copy the array
   { 
      REPORT
      Tracer tr("CroutMatrix::get_aux");
      int *ix = new int [nrows_val]; MatrixErrorNoSpace(ix);
      MONITOR_INT_NEW("Index (CM::get_aux)", nrows_val, ix)
      int n = nrows_val; int* i = ix; int* j = indx;
      while(n--) *i++ = *j++;
      X.indx = ix;
   }
}

CroutMatrix::CroutMatrix(const CroutMatrix& gm) : GeneralMatrix()
{
   REPORT
   Tracer tr("CroutMatrix(const CroutMatrix&)");
   ((CroutMatrix&)gm).get_aux(*this);
   GetMatrix(&gm);
}

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

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


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



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

// ************ Constructors for use with NR in C++ interface ***********

#ifdef SETUP_C_SUBSCRIPTS

Matrix::Matrix(Real a, int m, int n) : GeneralMatrix(m * n)
   { REPORT nrows_val=m; ncols_val=n; operator=(a); }
   
Matrix::Matrix(const Real* a, int m, int n) : GeneralMatrix(m * n)
   { REPORT nrows_val=m; ncols_val=n; *this << a; }

#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_val=nr; ncols_val=nc; tag_val=-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 SquareMatrix::resize(int n)
{ REPORT GeneralMatrix::resize(n,n,n*n); }

void SquareMatrix::resize(int nr, int nc)
{
   REPORT
   Tracer tr("SquareMatrix::resize");
   if (nc != nr) Throw(NotSquareException(*this));
   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 SquareMatrix::resize(const GeneralMatrix& A)
{
   REPORT
   int n = A.Nrows();
   if (n != A.Ncols())
   {
      Tracer tr("SquareMatrix::resize(GM)");
      Throw(NotSquareException(*this));
   }
   resize(n);
}

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"));
}

//*********************** resize_keep *******************************

void Matrix::resize_keep(int nr, int nc)
{
   Tracer tr("Matrix::resize_keep");
   if (nr == nrows_val && nc == ncols_val) { REPORT return; }
   
   if (nr <= nrows_val && nc <= ncols_val)
   {
      REPORT
      Matrix X = submatrix(1,nr,1,nc);
      swap(X);
   }
   else if (nr >= nrows_val && nc >= ncols_val)
   {
      REPORT
      Matrix X(nr, nc); X = 0;
      X.submatrix(1,nrows_val,1,ncols_val) = *this;
      swap(X);
   }
   else
   {
      REPORT
      Matrix X(nr, nc); X = 0;
      if (nr > nrows_val) nr = nrows_val;
      if (nc > ncols_val) nc = ncols_val;
      X.submatrix(1,nr,1,nc) = submatrix(1,nr,1,nc);
      swap(X);
   }
} 

void SquareMatrix::resize_keep(int nr)
{
   Tracer tr("SquareMatrix::resize_keep");