// Gestion de block de donnees avec partage de references
// malheureusement tres mal concu...  C.Magneville 04/99
// LAL (Orsay) / IN2P3-CNRS  DAPNIA/SPP (Saclay) / CEA
#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include <complex>
#include "pexceptions.h"
#include "ndatablock.h"
#include "objfio.h"

using namespace PlanckDPC;

// define DEBUG_NDATABLOCK

#ifdef DEBUG_NDATABLOCK
  static size_t NallocData = 0;
  static size_t NallocSRef = 0;
#endif

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

template <class T>
NDataBlock<T>::NDataBlock(size_t n)
// Createur d'une structure de "n" donnees
: mSz(0), mSRef(NULL), mIsTemp(false)
{
#ifdef DEBUG_NDATABLOCK
cout<<"?_NDataBlock::NDataBlock("<<this<<",n="<<n<<")"<<endl;
#endif

Alloc(n);
}

template <class T>
NDataBlock<T>::NDataBlock(size_t n, T* data, Bridge* br)
// Createur d'une structure de "n" donnees, avec donnees preallouees.
// Attention createur TRES DANGEREUX (Voir explications dans Alloc()).
: mSz(0), mSRef(NULL), mIsTemp(false)
{
#ifdef DEBUG_NDATABLOCK
cout<<"?_NDataBlock::NDataBlock("<<this
    <<",data="<<data<<",br="<<br<<")"<<endl;
#endif

Alloc(n,data,br);
}

template <class T>
NDataBlock<T>::NDataBlock()
// Createur par default
: mSz(0), mSRef(NULL), mIsTemp(false)
{
#ifdef DEBUG_NDATABLOCK
cout<<"?_NDataBlock::NDataBlock("<<this<<") default"<<endl;
#endif
}

template <class T>
NDataBlock<T>::NDataBlock(const NDataBlock<T>& a)
// Createur par copie: partage les donnees si "a" temporaire, clone sinon.
: mSz(0), mSRef(NULL), mIsTemp(false)
{
#ifdef DEBUG_NDATABLOCK
cout<<"?_NDataBlock::NDataBlock("<<this<<",&a="<<&a<<")"<<endl;
#endif

CloneOrShare(a);
}

template <class T>
NDataBlock<T>::NDataBlock(const NDataBlock<T>& a,bool share)
// Createur avec choix de partager ou non selon "share"
: mSz(0), mSRef(NULL), mIsTemp(false)
{
#ifdef DEBUG_NDATABLOCK
cout<<"?_NDataBlock::NDataBlock("<<this<<",&a="<<&a<<",sh=<<"<<share<<")"<<endl;
#endif

if(share) Share(a); else Clone(a);
}

template <class T>
NDataBlock<T>::~NDataBlock()
// Destructeur
{
#ifdef DEBUG_NDATABLOCK
cout<<"?_NDataBlock::~NDataBlock("<<this<<")"<<endl;
#endif

Delete();
}

////////////////////////////////////////////////////////////////
//************ Gestion de donnees

template <class T>
void NDataBlock<T>::Alloc(size_t n,T* data,Bridge* br)
// Allocation d'un NOUVEL espace de stoquage de "n" donnees
// Si data==NULL : allocation de l'espace memoire (remplit de zeros)
//    data!=NULL : partage des donnees avec l'adresse data
// Si br==NULL   : les donnees nous appartiennent
//    br!=NULL   : les donnees ne nous appartiennent pas (ex: Blitz)
//
// Exemple: on veut connecter a un tableau de T*
// > float *x = new float[5]; ... remplissage de x[] ...;
// 1- On veut que NDataBlock NE DESALLOUE PAS le tableau "x[]"
//    a- Premiere solution
//       > NDataBlock A(5,x,new Bridge);
//       ......
//       > delete [] x;
//          - Il faut deleter x[] explicitement.
//          - Le destructeur de "A" ne detruit pas x[].
//          ATTENTION: Une fois x[] detruit, "A" ne peut
//                     plus acceder les donnees!
//          - Bridge est detruit par le destructeur de "A"
//    b- Autre solution:
//       > NDataBlock A(5); A.FillFrom(5,x);
//       > delete [] x;
//       ......
//          - Il faut deleter x[] explicitement.
//          - "A" possede une copie en local de x[].
//          - Le destructeur de "A" ne detruit pas x[] mais la copie locale.
// 2- On veut que NDataBlock desalloue le tableau
//       > NDataBlock A(5,x);
//          - Ne Pas Faire "delete [] x;"
//          - "A" partage les donnees avec x[].
//          - Le destructeur de "A" detruit x[].
//
// --- REMARQUE SUR LE DANGER DE CERTAINES SITUATIONS (CMV):
// 1-/ x = new float[n1]; NDataBlock A(n2,x);
//     1er danger: si n2>n1 depassement de tableaux (core dump)
//     2sd danger: celui qui alloue x[] ne doit pas faire le "delete"
//                 en desaccord avec toutes les regles de bonne conduite.
// 2-/ float x[5]={1,2,3,4,5}; {NDataBlock A(n2,&x[0]);} cout<<x[2];
//     Ici, a la sortie du bloc {}, le destructeur de "A" va detruire
//     l'adresse de &x[0]: je n'ose imaginer que ca se fasse sans probleme
//     et de toute facon, cout<<x[2]; va surement faire des etincelles.
// 3-/ x = new float[n1]; NDataBlock A(n2,x,new Bridge);
//     1er danger: si n2>n1 depassement de tableaux (core dump)
//     2sd danger: si la methode bridgee (blitz?) detruit x[]
//                 "A" n'a plus de donnees connectees!
// --- CONCLUSION
// Cette classe est franchement merdique.
// - On peut accepter la prise de risque liee a NDataBlock(n2,x,new Bridge);
//   car je ne vois pas comment on pourrait faire autrement pour connecter
//   un tableau de type blitz par exemple.
// - Par contre le createur NDataBlock(n2,x); doit etre interdit
//   dans sa forme actelle car trop dangereux et il me semble inutile.
// - Dans cette nouvelle optique:
//   NDataBlock(n2,x,new Bridge) et NDataBlock(n2,x) disparaissent
//   On remplace par NDataBlock(n2,x) {Alloc(n2,x,new Bridge);}
//      qui force le Bridge dans tout les cas puisque NDataBlock
//      ne possede pas les donnees.
//   Mais puis-je encore le faire vu que NDataBlock est a la base
//   de TVector,TMatrix et qu'il faut donc reprendre tout le code DPC
// - Quoiqu'il arrive Alloc est une methode privee et peut donc rester
//   sous sa forme actuelle.
//            
{
#ifdef DEBUG_NDATABLOCK
cout<<"?_NDataBlock::Alloc("<<this<<","
    <<n<<","<<data<<","<<br<<") mSz="<<mSz
    <<" mSRef="<<mSRef<<" IsTemp="<<mIsTemp<<endl;
#endif

if(br && !data)
  throw(NullPtrError("NDataBlock::Alloc br!=NULL && data==NULL\n"));
if(n==0) throw(SzMismatchError("NDataBlock::Alloc n==0\n"));
if(mSRef) Delete();
mSz = n;
mSRef = new NDREF;
mSRef->nref = 1;
if(data) mSRef->data = data;
else {mSRef->data = new T[n]; memset(mSRef->data,0,n*sizeof(T));}
mSRef->bridge = br;

#ifdef DEBUG_NDATABLOCK
// Meme dans le cas data!=0 et br==0 (connexion d'un tableau
// avec destruction geree par ~NDataBlock (cas 2-) on compte
// comme si on avait fait une allocation du tableau (ce qui a ete
// fait au niveau du dessus!).
if(!br) NallocData++; NallocSRef++;
cout<<"...?_NDataBlock::Alloc mSz="<<mSz<<" mSRef="<<mSRef
    <<" mSRef->nref="<<mSRef->nref<<" mSRef->data="<< mSRef->data
    <<" mSRef->bridge="<<mSRef->bridge<<" IsTemp="<<mIsTemp
    <<" Total("<<NallocData<<","<<NallocSRef<<")"<<endl;
#endif
}

template <class T>
void NDataBlock<T>::Clone(const NDataBlock<T>& a)
// Clone: copie de donnees a partir de "a"
{
#ifdef DEBUG_NDATABLOCK
cout<<"?_NDataBlock::Clone("<<this<<","<<&a<<") a.(mSz="
    <<a.mSz<<" mSRef="<<a.mSRef<<" IsTemp="<<a.IsTemp()
    <<"), mSz="<<mSz<<" mSRef="<<mSRef<<" IsTemp="<<mIsTemp<<endl;
#endif

if(&a==NULL) throw(NullPtrError("NDataBlock::Clone  &a==NULL\n"));
if(!a.mSRef || a.mSz==0) throw(SzMismatchError("NDataBlock::Clone a.mSz==0\n"));
Alloc(a.mSz);
memcpy(Data(),a.Data(),mSz*sizeof(T));
}

template <class T>
void NDataBlock<T>::CloneOrShare(const NDataBlock<T>& a)
// CloneOrShare: Share si "a" temporaire, Clone sinon.
{
#ifdef DEBUG_NDATABLOCK
cout<<"?_NDataBlock::CloneOrShare("<<this<<","<<&a<<")"<<endl;
#endif

if(&a==NULL) throw(NullPtrError("NDataBlock::CloneOrShare  &a==NULL\n"));
if(a.IsTemp()) Share(a); else Clone(a);
}

template <class T>
void NDataBlock<T>::Share(const NDataBlock<T>& a)
// Share: Partage les donnees avec "a"
{
#ifdef DEBUG_NDATABLOCK
cout<<"?_NDataBlock::Share("<<this<<","<<&a<<")";
if(&a!=NULL) cout<<" a.(mSz="<<a.mSz<<" mSRef="<<a.mSRef
                 <<" IsTemp="<<a.IsTemp()<<")";
cout<<", mSz="<<mSz<<" mSRef="<<mSRef<<" IsTemp="<<mIsTemp<<endl;
#endif

if(&a==NULL) throw(NullPtrError("NDataBlock::Share  &a==NULL\n"));
if(!a.mSRef || a.mSz==0) throw(NullPtrError("NDataBlock::Share a.mSz=0\n"));
if(mSRef) Delete();
mSz = a.mSz; mSRef = a.mSRef; mSRef->nref++;

#ifdef DEBUG_NDATABLOCK
cout<<"...?_NDataBlock::Share mSz="<<mSz<<" mSRef="<<mSRef
    <<" mSRef->nref="<<mSRef->nref<<" mSRef->data="<< mSRef->data
    <<" mSRef->bridge="<<mSRef->bridge<<" IsTemp="<<mIsTemp<<endl;
#endif
}

template <class T>
void NDataBlock<T>::Delete(void)
// Pour detruire les pointeurs en tenant compte des references
{
#ifdef DEBUG_NDATABLOCK
cout<<"?_NDataBlock::Delete("<<this<<") mSz="<<mSz
    <<" mSRef="<<mSRef<<" IsTemp="<<mIsTemp;
if(mSRef)
  cout<<" mSRef->nref="<<mSRef->nref<<" mSRef->data="
      <<mSRef->data<<" mSRef->bridge="<<mSRef->bridge;
cout<<endl;
#endif

if(mSRef==NULL) return;
mSRef->nref--;
if(mSRef->nref != 0) {

#ifdef DEBUG_NDATABLOCK
  cout<<"...?_NDataBlock::Delete() pas de desallocation il reste nref="
      <<mSRef->nref<<" Total("<<NallocData<<","<<NallocSRef<<")"<<endl;
#endif

  mSz = 0; mSRef=NULL;
  return;
}

#ifdef DEBUG_NDATABLOCK
if(!mSRef->bridge) NallocData--; NallocSRef--;
cout<<"...?_NDataBlock::Delete() desallocation complete il reste nref="
    <<mSRef->nref<<" Total("<<NallocData<<","<<NallocSRef<<")"<<endl;
#endif

// Si il y a un Bridge les donnees ne n'appartiennent pas, on detruit le Bridge
// sinon, les donnees ont ete allouees par nos soins, on libere l'espace
if(mSRef->bridge) delete mSRef->bridge; else delete [] mSRef->data;
mSRef->bridge=NULL; mSRef->data=NULL;
delete mSRef; mSRef=NULL; mSz = 0;
}

template <class T>
void NDataBlock<T>::FillFrom(size_t n,T* data)
// Remplissage par un tableau de donnees
// - Si classe vide : creation de l'espace memoire
// - Si classe connectee : on ecrit selon la longueur minimale
//                         (cad this->mSz ou "n")
{
if(data==NULL) throw(NullPtrError("NDataBlock::FillFrom  data==NULL\n"));
if(n==0) throw(ParmError("NDataBlock::FillFrom  n<=0\n"));
if(mSRef==NULL) Alloc(n);
if(mSz<n) n = mSz;
memcpy(Data(),data,n*sizeof(T));
}

template <class T>
void NDataBlock<T>::Realloc(size_t nnew,bool force)
// Re-allocation de "nnew" place memoire pour les donnees
// avec conservation des "nold" donnees precedentes si possible.
// "force" gere la re-allocation de la place memoire pour les donnees.
// Divers cas se presentent:
// a-/ *** nnew>nold force=quelconque ***
//     place re-allouee, donnees [0,nold[ copiees, surplus [nold,new[ mis a zero
// b-/ *** nnew<=nold force=true ***
//     place re-allouee, donnees [0,nnew[ copiees, pas de surplus
// c-/ *** nnew<=nold force=false ***
//     place non re-allouee, seule la valeur de la taille est diminuee
// - On tient compte du partage des donnees dans tous les cas.
// - Si il n'y a pas de donnees connectees a la classe, on re-alloue
//   dans tous les cas
{
if(nnew==0) throw(ParmError("NDataBlock::Realloc  n<=0\n"));

// Cas sans re-allocation memoire
if(mSRef && nnew<=mSz && ! force) { mSz=nnew; return;}

// Cas avec re-allocation memoire
size_t ncop;
if(!mSRef || mSz==0) ncop=0; else if(mSz<nnew) ncop=mSz; else ncop=nnew;
T* dataloc = new T[nnew];
if(ncop>0) memcpy(dataloc,mSRef->data,ncop*sizeof(T));
if(nnew>ncop) memset(dataloc+ncop,0,(nnew-ncop)*sizeof(T));
Alloc(nnew,dataloc,NULL); //Alloc gere partage de reference et bridge
}

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

template <class T>
void NDataBlock<T>::Print(ostream& os,size_t i1,size_t n) const
// Impression de n elements a partir de i1
{
size_t nr = 0;
T* p = NULL; Bridge* br = NULL;
if(mSRef) {nr = mSRef->nref; p = mSRef->data; br = mSRef->bridge;}
os<<"NDataBlock::Print("<<this<<",Sz="<<mSz<<",IsTemp="<<mIsTemp<<")\n"
  <<"            mSRef="<<mSRef<<"(nref="<<nr<<",data="<<p
  <<",bridge="<<br<<")"<<endl;
if(i1>=mSz || n<=0 || !p) return;
size_t i2 = i1+n; if(i2>mSz) i2=mSz;
size_t im = 1; bool enl=false;
while(i1<i2) {
  enl = false;
  os<<" "<<(*this)(i1);  i1++;
  if(im==8) {os<<"\n"; im=1; enl=true;} else im++;
}
if(!enl) os<<endl;
}

////////////////////////////////////////////////////////////////

template <class T>
T NDataBlock<T>::Sum(size_t i1,size_t n) const
// Somme des elements de i1 a i1+n-1
{
if(i1>=mSz) return 0;
if(n>mSz) n = mSz; if(n==0) n = mSz-i1;
T const *p=Begin()+i1, *pe=p+n;
T val = 0;
while (p<pe) val += *p++;
return val;
}

template <class T>
T NDataBlock<T>::Product(size_t i1,size_t n) const
// Produit des elements de i1 a i1+n-1
{
if(i1>=mSz) return 0;
if(n>mSz) n = mSz; if(n==0) n = mSz-i1;
T const *p=Begin()+i1, *pe=p+n;
T val = 0;
while (p<pe) val *= *p++;
return val;
}

////////////////////////////////////////////////////////////////
//**** Surcharge de = : NDataBlock=NDataBlock; NDataBlock=<T> b; 

template <class T>
NDataBlock<T>& NDataBlock<T>::operator = (const NDataBlock<T>& a)
// Affectation: partage des donnees si "a" temporaire, clone sinon.
{
#ifdef DEBUG_NDATABLOCK
cout<<"?_NDataBlock::operator=("<<this<<","<<&a<<") a.(mSz="
    <<a.mSz<<" mSRef="<<a.mSRef<<" IsTemp="<<a.IsTemp()
    <<"), mSz="<<mSz<<" mSRef="<<mSRef<<" IsTemp="<<mIsTemp<<endl;
#endif

if(this == &a) return *this;
if(a.mSz==0)
  throw(SzMismatchError("NDataBlock::operator=A null size\n"));
CloneOrShare(a);
return *this;
}

template <class T>
NDataBlock<T>& NDataBlock<T>::operator = (T v)
// Affectation de tous les elements a une constante "v"
{
#ifdef DEBUG_NDATABLOCK
cout<<"?_NDataBlock::operator=("<<this<<","<<v<<")"
    <<" mSz="<<mSz<<" mSRef="<<mSRef<<" IsTemp="<<mIsTemp<<endl;
#endif

if(mSz==0) throw(SzMismatchError("NDataBlock::operator=v null size\n"));
T *p=Begin(), *pe=End();
while (p<pe) *p++ = v;
return *this;
}

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

template <class T>
NDataBlock<T>& NDataBlock<T>::operator += (T b)
{
if(mSz==0) throw(SzMismatchError("NDataBlock::operator+=v null size\n"));
T *p=Begin(), *pe=End();
while (p<pe) *p++ += b;
return *this;
}

template <class T>
NDataBlock<T>& NDataBlock<T>::operator -= (T b)
{
if(mSz==0) throw(SzMismatchError("NDataBlock::operator-=v null size\n"));
T *p=Begin(), *pe=End();
while (p<pe) *p++ -= b;
return *this;
}

template <class T>
NDataBlock<T>& NDataBlock<T>::operator *= (T b)
{
if(mSz==0) throw(SzMismatchError("NDataBlock::operator*=v null size\n"));
T *p=Begin(), *pe=End();
while (p<pe) *p++ *= b;
return *this;
}

template <class T>
NDataBlock<T>& NDataBlock<T>::operator /= (T b)
{
if(b==(T) 0) throw(ParmError("NDataBlock::operator/=v divide by zero\n"));
if(mSz==0) throw(SzMismatchError("NDataBlock::operator/=v null size\n"));
T *p=Begin(), *pe=End();
while (p<pe) *p++ /= b;
return *this;
}

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

template <class T>
NDataBlock<T>& NDataBlock<T>::operator += (const NDataBlock<T>& a)
{
if(mSz==0 || mSz!=a.mSz)
  throw(SzMismatchError("NDataBlock::operator+=A size mismatch/null"));
T *p=Begin(), *pe=End();
T const * pa=a.Begin();
while (p<pe) *p++ += *pa++;
return *this;
}

template <class T>
NDataBlock<T>& NDataBlock<T>::operator -= (const NDataBlock<T>& a)
{
if(mSz==0 || mSz!=a.mSz)
  throw(SzMismatchError("NDataBlock::operator-=A size mismatch/null"));
T *p=Begin(), *pe=End();
T const *pa=a.Begin();
while (p<pe) *p++ -= *pa++;
return *this;
}

template <class T>
NDataBlock<T>& NDataBlock<T>::operator *= (const NDataBlock<T>& a)
{
if(mSz==0 || mSz!=a.mSz)
  throw(SzMismatchError("NDataBlock::operator*=A size mismatch/null"));
T *p=Begin(), *pe=End();
T const *pa=a.Begin();
while (p<pe) *p++ *= *pa++;
return *this;
}

template <class T>
NDataBlock<T>& NDataBlock<T>::operator /= (const NDataBlock<T>& a)
// Attention, aucune protection si un element de "a" est nul.
{
if(mSz==0 || mSz!=a.mSz)
  throw(SzMismatchError("NDataBlock::operator/=A size mismatch/null"));
T *p=Begin(), *pe=End();
T const *pa=a.Begin();
while (p<pe) *p++ /= *pa++;
return *this;
}

////////////////////////////////////////////////////////////////
//**** Surcharge de +,-,*,/ : NDataBlock = NDataBlock+<T>b;
//                            NDataBlock = <T>b+NDataBlock;

template <class T>
NDataBlock<T> NDataBlock<T>::Add(T b) const
// Pour A+b
{
NDataBlock<T> result(*this); result.SetTemp(true);
result += b;
return result;
}

template <class T>
NDataBlock<T> NDataBlock<T>::Sub(T b) const
// Pour A-b
{
NDataBlock<T> result(*this); result.SetTemp(true);
return result -= b;
}

template <class T>
NDataBlock<T> NDataBlock<T>::SubInv(T b) const
// Pour b-A
{
NDataBlock<T> result(*this); result.SetTemp(true);
T *p=result.Begin(), *pe=result.End();
T const *pa=this->Begin();
while(p<pe) {*p++ = b - *pa++;}
return result;
}

template <class T>
NDataBlock<T> NDataBlock<T>::Mul(T b) const
// Pour A*b
{
NDataBlock<T> result(*this); result.SetTemp(true);
return result *= b;
}

template <class T>
NDataBlock<T> NDataBlock<T>::Div(T b) const
// Pour A/b
{
NDataBlock<T> result(*this); result.SetTemp(true);
return result /= b;
}

template <class T>
NDataBlock<T> NDataBlock<T>::DivInv(T b) const
// Pour b/A
{
NDataBlock<T> result(*this); result.SetTemp(true);
T *p=result.Begin(), *pe=result.End();
T const *pa = this->Begin();
while(p<pe) {*p++ = b / *pa++;}
return result;
}

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

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

template <class T>
NDataBlock<T> NDataBlock<T>::Mul(const NDataBlock<T>& b) const
// Pour A*B
{
if(mSz!=b.mSz)
  throw(SzMismatchError("NDataBlock operator C=A*B size mismatch/null\n"));
NDataBlock<T> result; result.SetTemp(true);
if(b.IsTemp()) {result.Share(b);            result *= *this;}
  else         {result.CloneOrShare(*this); result *= b;}
return result;
}

template <class T>
NDataBlock<T> NDataBlock<T>::Sub(const NDataBlock<T>& b) const
// Pour A-B
{
if(mSz!=b.mSz)
  throw(SzMismatchError("NDataBlock operator C=A-B size mismatch/null\n"));
NDataBlock<T> result; result.SetTemp(true);
if(b.IsTemp()) {
  result.Share(b);
  T *p=result.Begin(), *pe=result.End(); T const *pa=Begin();
  while(p<pe) {*p = *pa++  - *p; p++;}
} else {result.CloneOrShare(*this); result -= b;}
return result;
}

template <class T>
NDataBlock<T> NDataBlock<T>::Div(const NDataBlock<T>& b) const
// Pour A/B
{
if(mSz!=b.mSz)
  throw(SzMismatchError("NDataBlock operator C=A/B size mismatch/null\n"));
NDataBlock<T> result; result.SetTemp(true);
if(b.IsTemp()) {
  result.Share(b);
  T *p=result.Begin(), *pe=result.End(); T const *pa=Begin();
  while(p<pe) {*p = *pa++  / *p; p++;}
} else {result.CloneOrShare(*this); result /= b;}
return result;
}

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

/*
template <class T>
void ObjFileIO< NDataBlock<T> >::ReadSelf(PInPersist& is)
template <class T>
void ObjFileIO< NDataBlock<T> >::WriteSelf(POutPersist& os)
*/

// Pour pouvoir ecrire des tableaux de complex, en attendant 
// PIn/POutPersist::Get/Put(complex<>)
#include "piocmplx.h"

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

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

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

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

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

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


template <class T>
void FIO_NDataBlock<T>::ReadSelf(PInPersist& is)
{
// On lit les 3 premiers uint_8
uint_8 itab[3];
is.Get(itab, 3);
if (dobj == NULL) dobj = new NDataBlock<T>(itab[1]);
else if (itab[1] != dobj->Size()) dobj->ReSize(itab[1]);
// On lit le tableau de nombres
PIOSReadArray(is, dobj->Data(), dobj->Size()); 
}


template <class T>
void FIO_NDataBlock<T>::WriteSelf(POutPersist& os) const
{
if (dobj == NULL)   return;  // Attention - $CHECK$ Reza 26/04/99
//  On ecrit 3 uint_8 
//  0 : Numero de version,  1 : Taille,  2  reserve a l
uint_8 itab[3];
itab[0] = 1;
itab[1] = dobj->Size();
itab[2] = 0;
os.Put(itab, 3);
//  On ecrit le tableau de nombres 
PIOSWriteArray(os, dobj->Data(), dobj->Size()); 
}

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

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