// $Id: tmatrix.cc,v 1.6 1999-05-18 12:23:13 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"
#include "generalfit.h"

using namespace PlanckDPC;

////////////////////////////////////////////////////////////////
//**** Createur, Destructeur

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 si "a" temporaire).
: mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock)
{
}

template <class T>
TMatrix<T>::TMatrix(const TMatrix<T>& a,bool share)
// Constructeur par copie avec possibilite de forcer le partage ou non.
: mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock,share)
{
}

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

////////////////////////////////////////////////////////////////
// Operations matricielles
template <class T>
TMatrix<T> TMatrix<T>::Transpose(void) const
// Transposition
{
TMatrix<T> a; a.Clone(*this); a.SetTemp(true);
a.mNr = mNc; a.mNc = mNr;
{for(uint_4 i=0; i<mNr; i++)
  for(uint_4 j=0; j<mNc; j++) {
    a(j,i) = (*this)(i,j);
}}
return a;
}

////////////////////////////////////////////////////////////////
//**** 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 *= TMatrix;

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

////////////////////////////////////////////////////////////////
//**** Pour surcharge d'operateurs C = A (+,-,*,/) B

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

template <class T> TMatrix<T> TMatrix<T>::Sub(const TMatrix<T>& b) const
{
if(mNr!=b.mNr || mNc!=b.mNc)
    throw(SzMismatchError("TMatrix operator C=A-B size mismatch\n"));
TMatrix<T> result; result.SetTemp(true); result.mNr=mNr; result.mNc=mNc;
result.mNDBlock = mNDBlock-b.mNDBlock;
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("TMatrix operator C=A*B size mismatch\n"));
TMatrix<T> r; r.SetTemp(true); r.ReSize(mNr,b.mNc);
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;
}

#include "matrix.h"
////////////////////////////////////////////////////////////////
//**** Pour inversion
TMatrix<r_8> TMatrix<r_8>::Inverse() const
// Inversion
{
TMatrix<r_8> b(mNc,mNr);
TMatrix<r_8> a(*this);
b = 1.;
if (fabs(TMatrix<r_8>::GausPiv(a,b)) < 1.e-50)
  throw(MathExc("TMatrix Inverse() Singular Matrix"));
return b;
}

double TMatrix<r_8>::GausPiv(TMatrix<r_8>& a, TMatrix<r_8>& b)
// Pivot de Gauss
{
Matrix A(a);
Matrix B(b);
return Matrix::GausPiv(A,B);
}

//////////////////////////////////////////////////////////
//**** Residus des fits
TMatrix<r_8> TMatrix<r_8>::FitResidus(GeneralFit& gfit
            ,double xorg,double yorg,double dx,double dy)
// Retourne une classe contenant les residus du fit ``gfit''.
// On suppose que x=j (colonnes) et y=i (lignes) pour m(i,j).
// Les coordonnees de l'element (i,j) sont :
//    (i,j) -> x = xorg+j*dx , y = yorg+i*dy
{
if(NCols()<=0||NRows()<=0)
  throw(SzMismatchError("TMatrix::FitResidus size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("TMatrix::FitResidus GeneraFit==NULL\n"));
int npar =  gfit.GetNPar();
if(npar==0)
  throw(SzMismatchError("TMatrix::FitResidus GeneraFit 0 parametre\n"));
double* par = new double[npar];
{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
TMatrix<r_8> m(*this);
for(uint_4 i=0;i<NRows();i++) for(uint_4 j=0;j<NCols();j++) {
  double x[2] = {xorg+j*dx,yorg+i*dy};
  m(i,j) -= f->Value(x,par);
}
delete [] par;
return m;
}

TMatrix<r_8> TMatrix<r_8>::FitFunction(GeneralFit& gfit
            ,double xorg,double yorg,double dx,double dy)
// Retourne une classe contenant la fonction du fit ``gfit''.
// On suppose que x=j (colonnes) et y=i (lignes) pour m(i,j).
// Les coordonnees de l'element (i,j) sont :
//    (i,j) -> x = xorg + j*dx , y = yorg + i*dy
{
if(NCols()<=0||NRows()<=0)
  throw(SzMismatchError("TMatrix::FitFunction size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("TMatrix::FitFunction GeneraFit==NULL\n"));
int npar =  gfit.GetNPar();
if(npar==0)
  throw(SzMismatchError("TMatrix::FitFunction GeneraFit 0 parametre\n"));
double* par = new double[npar];
{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
TMatrix<r_8> m(*this);
for(uint_4 i=0;i<NRows();i++) for(uint_4 j=0;j<NCols();j++) {
  double x[2] = {xorg+j*dx,yorg+i*dy};
  m(i,j) = f->Value(x,par);
}
delete [] par;
return m;
}

////////////////////////////////////////////////////////////////
// -------------------------------------------------------------------------
//   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
//  0: Numero de version,  1 : NRows,  2 : NCol
uint_4 itab[3];
is.Get(itab,3);
if (dobj == NULL) dobj = new TMatrix<T>(itab[1],itab[2]);
else dobj->ReSize(itab[1],itab[2]);
// On lit le NDataBlock
FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
fio_nd.Read(is);
}

template <class T>
void FIO_TMatrix<T>::WriteSelf(POutPersist& os) const
{
if (dobj == NULL)   return;
//  On ecrit 3 uint_4 .... 
//  0: Numero de version,  1 : NRows,  2 : NCol
uint_4 itab[3];
 itab[0] = 1;  // Numero de version a 1
itab[1] = dobj->NRows();
itab[2] = dobj->NCols();
os.Put(itab,3);
// On ecrit le NDataBlock
FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
fio_nd.Write(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