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

// Copyright (C) 1991: R B Davies and DSIR

#include "include.hxx"

#include "newmat.hxx"
#include "newmatrc.hxx"

//#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 {}

//#define MONITOR(what,storage,store) \
//   { cout << what << " " << storage << " at " << (long)store << "\n"; }

#define MONITOR(what,store,storage) {}

/*************************** 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(int s)
{
   REPORT1
   storage=s; tag=-1;
   store = new real [storage]; MatrixErrorNoSpace(store);
   MONITOR("Make (GenMatrix)",storage,store)
}

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

SymmetricMatrix::SymmetricMatrix(int n) : GeneralMatrix(tristore(n))
{ REPORT1 nrows=n; ncols=n; }

UpperTriangularMatrix::UpperTriangularMatrix(int n)
   : GeneralMatrix(tristore(n))
{ REPORT1 nrows=n; ncols=n; }

LowerTriangularMatrix::LowerTriangularMatrix(int n)
   : GeneralMatrix(tristore(n))
{ REPORT1 nrows=n; ncols=n; }

DiagonalMatrix::DiagonalMatrix(int m) : GeneralMatrix(m)
{ REPORT1 nrows=m; ncols=m; }

Matrix::Matrix(BaseMatrix& M)
{
   REPORT1 CheckConversion(M);
   GeneralMatrix* gmx=M.Evaluate(MatrixType::Rect); GetMatrix(gmx);
}

RowVector::RowVector(BaseMatrix& M) : Matrix(M)
{
   if (nrows!=1) MatrixError("Illegal conversion to row vector"); }

ColumnVector::ColumnVector(BaseMatrix& M) : Matrix(M)
{
   if (ncols!=1) MatrixError("Illegal conversion to column vector"); }

SymmetricMatrix::SymmetricMatrix(BaseMatrix& M)
{
   REPORT1 CheckConversion(M);
   GeneralMatrix* gmx=M.Evaluate(MatrixType::Sym); GetMatrix(gmx);
}

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

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

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

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

CroutMatrix::CroutMatrix(BaseMatrix& m)
{
   REPORT1
   GeneralMatrix* gm = m.Evaluate(MatrixType::Rect); GetMatrix(gm);
   if (nrows!=ncols) MatrixError("Matrix not square");
   d=TRUE; indx=new int [nrows]; MatrixErrorNoSpace(indx);
   ludcmp();
}

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


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

void GeneralMatrix::ReDimension(int nr, int nc, int s)
{
   REPORT 
   if (store)
   {
      MONITOR("Free (ReDimensi)",storage,store)
      delete [storage] store;
   }
   storage=s; nrows=nr; ncols=nc; tag=-1;
   store = new real [storage]; MatrixErrorNoSpace(store);
   MONITOR("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 GeneralMatrix* s) const
{ REPORT return (s==this) ? 1 : 0; }

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

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

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

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

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

MatrixType Matrix::Type() const { return MatrixType::Rect; }

MatrixType SymmetricMatrix::Type() const { return MatrixType::Sym; }

MatrixType UpperTriangularMatrix::Type() const { return MatrixType::UT; }

MatrixType LowerTriangularMatrix::Type() const { return MatrixType::LT; }

MatrixType DiagonalMatrix::Type() const { return MatrixType::Diag; }

MatrixType RowVector::Type() const { return MatrixType::RowV; }

MatrixType ColumnVector::Type() const { return MatrixType::ColV; }

MatrixType CroutMatrix::Type() const { return MatrixType::Crout; }

MatrixType RowedMatrix::Type() const { return MatrixType::RowV; }

MatrixType ColedMatrix::Type() const { return MatrixType::ColV; }

MatrixType DiagedMatrix::Type() const { return MatrixType::Diag; }

MatrixType MatedMatrix::Type() const { return MatrixType::Rect; }

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 MatrixType mteqel(MatrixType::EqEl); return bm->Type() + mteqel; }

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 ReturnMatrix::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 AddedMatrix::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(); }
int ConstMatrix::NrowsV() const { return cgm->Nrows(); }
int ReturnMatrix::NrowsV() const { return gm->Nrows(); }

int GeneralMatrix::NcolsV() const { return ncols; }
int ColedMatrix::NcolsV() const { return 1; }
int MatedMatrix::NcolsV() const { return nc; }
int GetSubMatrix::NcolsV() const { return col_number; }
int AddedMatrix::NcolsV() const  { return bm2->NcolsV(); }
int ShiftedMatrix::NcolsV() const { return bm->NcolsV(); }
int TransposedMatrix::NcolsV() const { return bm->NrowsV(); }
int NegatedMatrix::NcolsV() const { return bm->NcolsV(); }
int RowedMatrix::NcolsV() const { return bm->NrowsV() * bm->NcolsV(); }
int DiagedMatrix::NcolsV() const { return bm->NrowsV() * bm->NcolsV(); }
int ConstMatrix::NcolsV() const { return cgm->Ncols(); }
int ReturnMatrix::NcolsV() const { return gm->Ncols(); }


//  Rules regarding tDelete, reuse, GetStore
//    All matrices processed during expression evaluation must be subject
//    to exactly one of reuse(), tDelete(), GetStore() or BorrowStore().
//    If reuse returns TRUE the matrix must be reused.
//    GetMatrix(gm) always calls gm->GetStore()
//    gm->Evaluate obeys rules
//    bm->Evaluate obeys rules for matrices in bm structure

void GeneralMatrix::tDelete()
{
   if (tag<0)
   {
      if (tag<-1) { REPORT store=0; delete this; return; }  // borrowed
      else { REPORT return; }
   }
   if (tag==1)
   {
      REPORT  MONITOR("Free   (tDelete)",storage,store)
      if (store) delete [storage] store; store=0; tag=-1; return;
   }
   if (tag==0) { REPORT delete this; return; }
   REPORT tag--; return;
}

BOOL GeneralMatrix::reuse()
{
   if (tag<-1)
   {
      REPORT
      real* s = new real [storage]; MatrixErrorNoSpace(s);
      MONITOR("Make     (reuse)",storage,s)
      real* s1=store; real* s2=s; int i=storage;
      while(i--) *s2++ = *s1++;
      store=s; tag=0; return TRUE;
   }
   if (tag<0) { REPORT return FALSE; }
   if (tag<=1)  { REPORT return TRUE; }
   REPORT tag--; return FALSE;
}

real* GeneralMatrix::GetStore()
{
   if (tag<0 || tag>1)
   {
      real* s = new real [storage]; MatrixErrorNoSpace(s);
      MONITOR("Make  (GetStore)",storage,s)
      real* s1=store; real* s2=s; int i=storage;
      while(i--) *s2++ = *s1++;
      if (tag>1) { REPORT tag--; }
      else if (tag < -1) { REPORT store=0; delete this; } // borrowed store
      else { REPORT }
      return s;
   }
   real* s=store; store=0;
   if (tag==0) { REPORT delete this; }
   else { REPORT tag=-1; }
   return s;
}

#ifndef __ZTC__
void GeneralMatrix::GetMatrixC(const GeneralMatrix* gmx)
{
   REPORT tag=-1;
   nrows=gmx->nrows; ncols=gmx->ncols; storage=gmx->storage;

   store = new real [storage]; MatrixErrorNoSpace(store);
   MONITOR("Make (GetMatrix)",storage,store)
   real* s1=gmx->store; real* s2=store; int i=storage;
   while(i--) *s2++ = *s1++;
}
#endif

void GeneralMatrix::GetMatrix(GeneralMatrix* gmx)
{
   REPORT  tag=-1; nrows=gmx->Nrows(); ncols=gmx->Ncols();
   storage=gmx->storage; store=gmx->GetStore(); 
}

GeneralMatrix* GeneralMatrix::BorrowStore(GeneralMatrix* gmx, MatrixType mt)
// Copy storage of *this to storage of *gmx. Then convert to type mt.
// If mt == 0 just let *gm point to storage of *this if tag<0.
{
   if ((int)mt==0)
   {
      if (tag == -1) { REPORT gmx->tag = -2; gmx->store = store; }
      else { REPORT gmx->tag = 0; gmx->store = GetStore(); }
   }
   else if (mt!=gmx->Type())
   {
      REPORT gmx->tag = -2; gmx->store = store;
      gmx = gmx->Evaluate(mt); gmx->tag = 0; tDelete();
   }
   else { REPORT  gmx->tag = 0; gmx->store = GetStore(); }

   return gmx;
}   

void GeneralMatrix::Eq(BaseMatrix& X, MatrixType mt)
// Count number of references to this in X.
// If zero delete storage in X;
// otherwise tag X to show when storage can be deleted
// evaluate X and copy to current object
{
   int counter=X.search(this);
   if (counter==0)
   {
      REPORT  MONITOR("Free (operator=)",storage,store)
      if (store) { REPORT delete [storage] store; storage=0; }
   }
   else { REPORT Release(counter); }
   GeneralMatrix* gmx = X.Evaluate(mt);
   if (gmx!=this) { REPORT GetMatrix(gmx); }
   else { REPORT }
   Protect();
}

void GeneralMatrix::Inject(const GeneralMatrix& X)
// copy stored values of X; otherwise leave els of *this unchanged
{
   REPORT
   CheckConversion(X);
   if (nrows != X.nrows || ncols != X.ncols)
      MatrixError("Inject: dimensions don't match");
   MatrixRow mr((GeneralMatrix*)&X, LoadOnEntry);
   MatrixRow mrx(this, LoadOnEntry+StoreOnExit+DirectPart);
   int i=nrows;
   while (i--) { mrx.Inject(mr); mrx.Next(); mr.Next(); }
}  

void GeneralMatrix::CheckConversion(const BaseMatrix& M)
{ if (!(this->Type() >= M.Type())) MatrixError("Illegal Conversion"); }

/************************* nricMatrix routines *****************************/

void nricMatrix::MakeRowPointer()
{
   row_pointer = new real* [nrows]; MatrixErrorNoSpace(row_pointer);
   real* s = Store() - 1; int i = nrows; real** rp = row_pointer;
   while (i--) { *rp++ = s; s+=ncols; }
}

void nricMatrix::DeleteRowPointer()
{ if (nrows) delete [nrows] row_pointer; }

void GeneralMatrix::CheckStore() const
{ if (!store) MatrixError("NRIC accessing matrix with dimensions not set"); }
