newmat4.cxx

科学和工程计算中使用统计功能开发工具包.

//$$ newmat4.cxx       Constructors, ReDimension, basic utilities

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

#include "include.h"

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

//#define REPORT { static ExeCounter ExeCount(__LINE__,4); ExeCount++; }

#define REPORT {}

//#define REPORT1 { static ExeCounter ExeCount(__LINE__,4); ExeCount++; }

// REPORT1 constructors only - doesn't work in turbo and Borland C++

#define REPORT1 {}


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

static int tristore(int n)                      // els 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)
{
   REPORT1
   storage=s.Value(); tag=-1;
   store = new Real [storage]; MatrixErrorNoSpace(store);
   MONITOR_REAL_NEW("Make (GenMatrix)",storage,store)
}

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

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

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

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

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

Matrix::Matrix(const BaseMatrix& M)
{
   REPORT1 CheckConversion(M);
   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)
{
   REPORT1 CheckConversion(M);
   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Sm);
   GetMatrix(gmx);
}

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

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

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

GeneralMatrix::~GeneralMatrix()
{
   if (store)
   {
      MONITOR_REAL_DELETE("Free (GenMatrix)",storage,store)
#ifdef Version21
      delete [] store;
#else
      delete [storage] store;
#endif
   }
}

CroutMatrix::CroutMatrix(const BaseMatrix& m)
{
   REPORT1
   Tracer tr("CroutMatrix");
   GeneralMatrix* gm = ((BaseMatrix&)m).Evaluate(MatrixType::Rt);
   GetMatrix(gm);
   if (nrows!=ncols) Throw(NotSquareException(*this));
   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)
#ifdef Version21
   delete [] indx;
#else
   delete [nrows] indx;
#endif
}

//ReturnMatrixX::ReturnMatrixX(GeneralMatrix& gmx)
//{
//   REPORT1
//   gm = gmx.Type().New();  MatrixErrorNoSpace(gm);
//   gm->GetMatrix(&gmx);  gm->ReleaseAndDelete();
//}

#ifndef TEMPS_DESTROYED_QUICKLY

GeneralMatrix::operator ReturnMatrixX() const
{
   REPORT
   GeneralMatrix* gm = Type().New();  MatrixErrorNoSpace(gm);
   gm->GetMatrix(this);  gm->ReleaseAndDelete();
   return ReturnMatrixX(gm);
}

#else

GeneralMatrix::operator ReturnMatrixX&() const
{
   REPORT
   GeneralMatrix* gm = Type().New();  MatrixErrorNoSpace(gm);
   gm->GetMatrix(this);  gm->ReleaseAndDelete();
   ReturnMatrixX* x = new ReturnMatrixX(gm);
   MatrixErrorNoSpace(x); return *x;
}

#endif

/**************************** ReDimension matrices ***************************/

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

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

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

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

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

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

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

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


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

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

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 ConstMatrix::search(const BaseMatrix* s) const
{ REPORT return (s==cgm) ? 1 : 0; }

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 RowedMatrix::Type() const { return MatrixType::RV; }
MatrixType ColedMatrix::Type() const { return MatrixType::CV; }
MatrixType DiagedMatrix::Type() const { return MatrixType::Dg; }
MatrixType MatedMatrix::Type() const { return MatrixType::Rt; }
MatrixType GetSubMatrix::Type() const { return mt; }

MatrixType AddedMatrix::Type() const
{ REPORT return bm1->Type() + bm2->Type(); }

MatrixType MultipliedMatrix::Type() const
{ REPORT return bm1->Type() * bm2->Type(); }

MatrixType SolvedMatrix::Type() const
{ REPORT return bm1->Type().i() * bm2->Type(); }

MatrixType SubtractedMatrix::Type() const
{ REPORT return bm1->Type() + bm2->Type(); }

MatrixType ShiftedMatrix::Type() const
{ REPORT return bm->Type().AddEqualEl(); }

MatrixType ScaledMatrix::Type() const { REPORT return bm->Type(); }
MatrixType TransposedMatrix::Type() const { REPORT return bm->Type().t(); }
MatrixType NegatedMatrix::Type() const { REPORT return bm->Type(); }
MatrixType InvertedMatrix::Type() const { REPORT return bm->Type().i(); }
MatrixType ConstMatrix::Type() const { REPORT return cgm->Type(); }
MatrixType ReturnMatrixX::Type() const { REPORT return gm->Type(); }

/*
int GeneralMatrix::NrowsV() const { return nrows; }
int RowedMatrix::NrowsV() const { return 1; }
int MatedMatrix::NrowsV() const { return nr; }
int GetSubMatrix::NrowsV() const { return row_number; }
int MultipliedMatrix::NrowsV() const  { return bm1->NrowsV(); }
int ShiftedMatrix::NrowsV() const { return bm->NrowsV(); }
int TransposedMatrix::NrowsV() const { return bm->NcolsV(); }
int NegatedMatrix::NrowsV() const { return bm->NrowsV(); }
int ColedMatrix::NrowsV() const { return bm->NrowsV() * bm->NcolsV(); }
int DiagedMatrix::NrowsV() const { return bm->NrowsV() * bm->NcolsV(); }