// $Id: tmatrix.cc,v 1.1 1999-04-29 08:57:15 ansari Exp $
//                         C.Magneville          04/99
#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include <complex>
#include "pexceptions.h"
#include "tmatrix.h"
#include "objfio.h"

using namespace PlanckDPC;

////////////////////////////////////////////////////////////////
//************ Createur, Destructeur, gestion des donnees

template <class T>
TMatrix<T>::TMatrix()
// Constructeur par defaut.
: mNr(0), mNc(0), mNDBlock(NULL)
{
}

template <class T>
TMatrix<T>::TMatrix(uint_4 r,uint_4 c)
// Construit une matrice de r lignes et c colonnes.
: mNr(r), mNc(c), mNDBlock(new NDataBlock<T>(r*c))
{
}

template <class T>
TMatrix<T>::TMatrix(uint_4 r,uint_4 c,T* values,Bridge* br)
// Construit une matrice de r lignes et c colonnes. On fournit
// le tableau des valeurs et eventuellement un Bridge.
: mNr(r), mNc(c), mNDBlock(new NDataBlock<T>(r*c,values,br))
{
}

template <class T>
TMatrix<T>::TMatrix(const TMatrix<T>& a)
// Constructeur par copie.
: mNr(a.mNr), mNc(a.mNc)
, mNDBlock(new NDataBlock<T>(*(a.DataBlock())))
{
}

template <class T>
TMatrix<T>::TMatrix(const TMatrix<T>& a,bool share)
// Constructeur par copie.
: mNr(a.mNr), mNc(a.mNc)
, mNDBlock(new NDataBlock<T>(*(a.DataBlock()),share))
{
}

template <class T>
TMatrix<T>::~TMatrix()
// Destructeur
{
Delete();
}

template <class T>
void TMatrix<T>::Delete(void)
// Pour desallouer
{
mNr = mNc = 0;
if(mNDBlock) delete mNDBlock; mNDBlock = NULL;
}

////////////////////////////////////////////////////////////////
template <class T>
void TMatrix<T>::Clone(const TMatrix<T>& a)
// Clone (copie de donnee) a partir de "a"
{
if(a.mNr==0 || a.mNc==0 || a.DataBlock()==NULL) {Delete(); return;}
if(mNDBlock) Delete();
else NDataBlock<T>* mNDBlock = new NDataBlock<T>(1); // cas "a" vide
mNr = a.mNr; mNc = a.mNc;
mNDBlock->Clone(*(a.DataBlock()));
}

template <class T>
void TMatrix<T>::Reset(T v)
// Reset de la matrice a "v"
{
if(mNDBlock==NULL)
  throw(NullPtrError("TMatrix::Reset mNDBlock==NULL\n"));
if(mNr==0 || mNc==0)
  throw(SzMismatchError("TMatrix::Reset mNr==0 || mNc==0\n"));
mNDBlock->Reset(v);
}

template <class T>
void TMatrix<T>::Realloc(uint_4 r,uint_4 c,bool force_alloc)
// Reallocation de place
{
if(r==0 || c==0) throw(SzMismatchError("TMatrix::ReSize r ou c==0\n"));
if(mNDBlock==NULL) NDataBlock<T>* mNDBlock = new NDataBlock<T>(1);
if(!force_alloc && mNr==r && mNc==c)  return;
mNr = r; mNc = c;
mNDBlock->ReSize(r*c);
}

////////////////////////////////////////////////////////////////
//**** Surcharge de TMatrix=TMatrix; TMatrix=<T> b; 

template <class T>
TMatrix<T>& TMatrix<T>::operator = (T x)
// Operateur d'affectation depuis scalaire : identite * scalaire.
{
if(mNr!=mNc || mNr==0 || mNDBlock==NULL)
  throw(SzMismatchError("TMatrix::operator= mNc!=mNr ou ==0 ou mNDBlock==NULL\n"));
for(int r=0;r<mNr;r++) for(int c=0;c<mNc;c++) (*this)(r,c) = (r==c)? x: 0;
return *this;
}

template <class T>
TMatrix<T>& TMatrix<T>::operator = (const TMatrix<T>& a)
// Operateur d'affectation A=B (sans reaffectation)
{
if(this == &a) return *this;
if(mNr!=a.mNc || mNc!=a.mNc || mNr==0 || mNc==0)
  throw(SzMismatchError("TMatrix::operator= mNc!=a.mNc ou mNr!=a.mNr ou ==0\n"));
if(mNDBlock==NULL)
  throw(SzMismatchError("TMatrix::operator= mNDBlock==NULL\n"));
*(mNDBlock) = *(a.DataBlock());
return *this;
}

////////////////////////////////////////////////////////////////
//**** Impression

template <class T>
void TMatrix<T>::Print(ostream& os,int lp
                 ,uint_4 i0,uint_4 ni,uint_4 j0,uint_4 nj) const
// Impression de la sous-matrice (i:i+ni-1,i:j+nj-1)
{
  os<<"TMatrix::Print("<<mNr<<","<<mNc<<")"<<endl;
if(lp>0)
  {os<<"       this="<<this<<" mNDBlock="<<mNDBlock<<endl;
   if(mNDBlock) mNDBlock->Print(0,0);}
if(mNr==0 || mNc==0) return;
if(i0>=mNr || j0>=mNc || ni==0 || nj==0) return;
uint_4 i1 = i0+ni; if(i1>mNr) i1 = mNr;
uint_4 j1 = j0+nj; if(j1>mNc) j1 = mNc;
for(int i=i0;i<i1;i++) {
  for(int j=j0;j<j1;j++) cout<<" "<<(*this)(i,j)<<endl;
}
}

////////////////////////////////////////////////////////////////
//**** Surcharge de +=,-=,*=,/= (INPLACE): TMatrix += <T> b;

template <class T>
TMatrix<T>& TMatrix<T>::operator += (T b)
{
if(mNr==0 || mNc==0 || !mNDBlock)
  throw(SzMismatchError("TMatrix::operator+=v mNr==0 || mNc==0 || !mNDBlock\n"));
*mNDBlock += b;
return *this;
}

template <class T>
TMatrix<T>& TMatrix<T>::operator -= (T b)
{
if(mNr==0 || mNc==0 || !mNDBlock)
  throw(SzMismatchError("TMatrix::operator-=v mNr==0 || mNc==0 || !mNDBlock\n"));
*mNDBlock -= b;
return *this;
}

template <class T>
TMatrix<T>& TMatrix<T>::operator *= (T b)
{
if(mNr==0 || mNc==0 || !mNDBlock)
  throw(SzMismatchError("TMatrix::operator*=v mNr==0 || mNc==0 || !mNDBlock\n"));
*mNDBlock *= b;
return *this;
}

template <class T>
TMatrix<T>& TMatrix<T>::operator /= (T b)
{
if(b==(T) 0) throw(ParmError("TMatrix::operator/=v divide by zero\n"));
if(mNr==0 || mNc==0 || !mNDBlock)
  throw(SzMismatchError("TMatrix::operator/=v mNr==0 || mNc==0 || !mNDBlock\n"));
*mNDBlock /= b;
return *this;
}

////////////////////////////////////////////////////////////////
//**** Surcharge de +=,-=,*= (INPLACE): TMatrix += TMatrix;

template <class T>
TMatrix<T>& TMatrix<T>::operator += (const TMatrix<T>& a)
{
if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc || !mNDBlock || !a.mNDBlock)
  throw(SzMismatchError("TMatrix::operator+=A size mismatch/null"));
*(mNDBlock) += *(a.mNDBlock);
return *this;
}

template <class T>
TMatrix<T>& TMatrix<T>::operator -= (const TMatrix<T>& a)
{
if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc || !mNDBlock || !a.mNDBlock)
  throw(SzMismatchError("TMatrix::operator-=A size mismatch/null"));
*(mNDBlock) -= *(a.mNDBlock);
return *this;
}

template <class T>
TMatrix<T>& TMatrix<T>::operator *= (const TMatrix<T>& a)
// A = A*B  -> A(n,m) = A(n,m)*B(m,m)
{
uint_4 ndata = mNr*mNc;
if(ndata==0 || mNc!=a.mNr || a.mNr!=a.mNc || !mNDBlock || !a.mNDBlock)
  throw(SzMismatchError("TMatrix::operator*=A size mismatch/null"));
// A(i,j) = Sum(k) A(i,k)*B(k,j) ... il faut sauver la ligne "i" de A
// Vecteur oi   : vecteur ou est sauve la ligne i de la matrice *this
//                oi,oe = pointeur de debut et de fin du vecteur temporaire
//                oij   = pointeur parcourant le vecteur oi
// Matrice *this: i     = pointeur du debut de la ligne i
//                ij    = pointeur parcourant la ligne i
// Matrice a    : aj    = pointeur de debut de la colonne j
//                aji   = pointeur parcourant la colonne j
T* oi = new T[mNc]; T* oe = oi+mNc;
for(T *i=Data(); i<Data()+ndata; i+=mNc) {
  {for(T *oij=oi,*ij=i; oij<oe;) *oij++ = *ij++;}
  {for(T *ij=i,*aj=const_cast<T *>(a.Data()); aj<a.Data()+a.mNc; ij++,aj++) {
    T sum = 0;
    for(T *oij=oi,*aji=aj; oij<oe; oij++,aji+=a.mNc) sum += *oij * *aji;
    *ij = sum;
  }}
}
delete [] oi;
return *this;
}

///////////////////////////////////////////////////////////////
#ifdef __CXX_PRAGMA_TEMPLATES__
/*
#pragma define_template TMatrix<uint_1>
#pragma define_template TMatrix<uint_2>
#pragma define_template TMatrix<int_2>
#pragma define_template TMatrix<int_4>
#pragma define_template TMatrix<int_8>
#pragma define_template TMatrix<uint_4>
#pragma define_template TMatrix<uint_8>
#pragma define_template TMatrix<r_4>
*/
#pragma define_template TMatrix<r_8>
/*
#pragma define_template TMatrix< complex<float> >
#pragma define_template TMatrix< complex<double> >
*/
// Instances des delegues FileIO (PPersist)
#endif


#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class TMatrix<uint_1>;
template class TMatrix<uint_2>;
template class TMatrix<int_2>;
template class TMatrix<int_4>;
template class TMatrix<int_8>;
template class TMatrix<uint_4>;
template class TMatrix<uint_8>;
template class TMatrix<r_4>;
template class TMatrix<r_8>;
template class TMatrix< complex<float> >;
template class TMatrix< complex<double> >;
// Instances des delegues FileIO (PPersist)
#endif