// $Id: tmatrix.cc,v 1.2 1999-04-30 11:02:52 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()
{
}

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(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(r*c,values,br)
{
}

template <class T>
TMatrix<T>::TMatrix(const TMatrix<T>& a)
// Constructeur par copie (partage des donnees).
: mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock)
{
}

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

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

////////////////////////////////////////////////////////////////
template <class T>
void TMatrix<T>::Clone(const TMatrix<T>& a)
// Clone (copie de donnees) a partir de "a"
// Partage des donnees si "a" est temporaire.
{
mNDBlock.Clone(a.mNDBlock); mNr = a.mNr; mNc = a.mNc;
}

template <class T>
void TMatrix<T>::Reset(T v)
// Reset de la matrice a "v"
{
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(!force_alloc && mNr==r && mNc==c)  return;
mNr = r; mNc = c; mNDBlock.ReSize(r*c,force_alloc);
}

////////////////////////////////////////////////////////////////
//**** 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)
  throw(SzMismatchError("TMatrix::operator= mNc!=mNr ou ==0\n"));
for(uint_4 r=0;r<mNr;r++) for(uint_4 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=a
// (Re)allocation de la place memoire
// sauf si "a" est cree temporairement (partage)
{
if(this == &a) return *this;
Clone(a);
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<<endl; 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(uint_4 i=i0;i<i1;i++) {
  for(uint_4 j=j0;j<j1;j++) cout<<" "<<(*this)(i,j);
  cout<<endl;
}
}

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

template <class T>
TMatrix<T>& TMatrix<T>::operator += (T b)
{
mNDBlock += b;
return *this;
}

template <class T>
TMatrix<T>& TMatrix<T>::operator -= (T b)
{
mNDBlock -= b;
return *this;
}

template <class T>
TMatrix<T>& TMatrix<T>::operator *= (T b)
{
mNDBlock *= b;
return *this;
}

template <class T>
TMatrix<T>& TMatrix<T>::operator /= (T b)
{
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)
  throw(SzMismatchError("TMatrix::operator+=A size mismatch"));
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)
  throw(SzMismatchError("TMatrix::operator-=A size mismatch"));
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)
  throw(SzMismatchError("TMatrix::operator*=A size mismatch"));
// 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;
}

////////////////////////////////////////////////////////////////
//**** Surcharge de +,-,*: TMatrix = TMatrix + TMatrix;

template <class T>
TMatrix<T> TMatrix<T>::Add(const TMatrix<T>& b) const
// C = A(this)+B
{
if(mNr!=b.mNr || mNc!=b.mNc)
  throw(SzMismatchError("NDataBlock operator C=A+B size mismatch\n"));
TMatrix<T> result; result.SetTemp(true);
if(b.IsTemp()) {
  result.Clone(b);
  result += *this;
} else {
  result.Clone(*this);
  result += b;
}
return result;
}

template <class T>
TMatrix<T> TMatrix<T>::Sub(const TMatrix<T>& b) const
// C = A(this)-B
{
if(mNr!=b.mNr || mNc!=b.mNc)
  throw(SzMismatchError("NDataBlock operator C=A-B size mismatch\n"));
TMatrix<T> result; result.SetTemp(true);
if(b.IsTemp()) {
  result.Clone(b);
  result.mNDBlock = (*this).mNDBlock - result.mNDBlock;
} else {
  result.Clone(*this);
  result -= b;
}
return result;
}

template <class T>
TMatrix<T> TMatrix<T>::Mul(const TMatrix<T>& b) const
// C = A(this)*B : Cij = Aik Bkj (allocation forcee dans tous les cas)
{
if(mNr==0 || mNc==0 || b.mNc==0 || mNc!=b.mNr)
  throw(SzMismatchError("NDataBlock operator C=A*B size mismatch\n"));
TMatrix<T> r; r.SetTemp(true); r.Realloc(mNr,b.mNc,true);
T *ai,*aik,*bj,*bkj,*ri,*rij;
for(ri=const_cast<T *>(r.Data()),ai=const_cast<T *>(Data());
    ri<r.Data()+r.mNr*r.mNc;ri+=r.mNc,ai+=mNc) {
  for(rij=ri,bj=const_cast<T *>(b.Data());rij<ri+r.mNc;rij++,bj++) {
    *rij = 0;
    for(aik=ai,bkj=bj;aik<ai+mNc;aik++,bkj+=b.mNc) *rij += *aik * *bkj;
  }
}
return r;
}

////////////////////////////////////////////////////////////////
// -------------------------------------------------------------------------
//   Les objets delegues pour la gestion de persistance 
// -------------------------------------------------------------------------

template <class T>
FIO_TMatrix<T>::FIO_TMatrix()
{
dobj=new TMatrix<T>;
ownobj=true;
}

template <class T>
FIO_TMatrix<T>::FIO_TMatrix(string const & filename) 
{
dobj=new TMatrix<T>;
ownobj=true; 
Read(filename);
}

template <class T>
FIO_TMatrix<T>::FIO_TMatrix(const TMatrix<T> & obj) 
{ 
dobj = new TMatrix<T>(obj);
ownobj=true; 
}

template <class T>
FIO_TMatrix<T>::FIO_TMatrix(TMatrix<T> * obj) 
{ 
dobj = obj;
ownobj=false; 
}

template <class T>
FIO_TMatrix<T>::~FIO_TMatrix()
{
if (ownobj && dobj) delete dobj;
}

template <class T>
AnyDataObj* FIO_TMatrix<T>::DataObj()
{
return(dobj);
}


template <class T>
void FIO_TMatrix<T>::ReadSelf(PInPersist& is)
{
// On lit les 3 premiers uint_8
uint_4 itab[2];
is.Get(itab,2);
if (dobj == NULL) dobj = new TMatrix<T>(itab[0],itab[1]);
else dobj->Realloc(itab[0],itab[1],false);
// On lit le NDataBlock
//cmv   encore des problemes avec les consteries
//FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
//fio_nd.ReadSelf(is);
}

template <class T>
void FIO_TMatrix<T>::WriteSelf(POutPersist& os) const
{
if (dobj == NULL)   return;
//  On ecrit 2 uint_4 .... 0 : NRows,  1 : NCol
uint_4 itab[2];
itab[0] = dobj->NRows();
itab[1] = dobj->NCol();
os.Put(itab,2);
// On ecrit le NDataBlock
//cmv   encore des problemes avec les consteries
//FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
//fio_nd.WriteSelf(os);
}

///////////////////////////////////////////////////////////////
#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)
#pragma define_template FIO_TMatrix<uint_1>
#pragma define_template FIO_TMatrix<uint_2>
#pragma define_template FIO_TMatrix<int_2>
#pragma define_template FIO_TMatrix<int_4>
#pragma define_template FIO_TMatrix<int_8>
#pragma define_template FIO_TMatrix<uint_4>
#pragma define_template FIO_TMatrix<uint_8>
#pragma define_template FIO_TMatrix<r_8>
#pragma define_template FIO_TMatrix<r_4>
#pragma define_template FIO_TMatrix< complex<float> >
#pragma define_template FIO_TMatrix< complex<double> >
#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)
template class FIO_TMatrix<uint_1>;
template class FIO_TMatrix<uint_2>;
template class FIO_TMatrix<int_2>;
template class FIO_TMatrix<int_4>;
template class FIO_TMatrix<int_8>;
template class FIO_TMatrix<uint_4>;
template class FIO_TMatrix<uint_8>;
template class FIO_TMatrix<r_8>;
template class FIO_TMatrix<r_4>;
template class FIO_TMatrix< complex<float> >;
template class FIO_TMatrix< complex<double> >;
#endif