Index: trunk/SophyaLib/NTools/tmatrix.cc
===================================================================
--- trunk/SophyaLib/NTools/tmatrix.cc	(revision 285)
+++ trunk/SophyaLib/NTools/tmatrix.cc	(revision 286)
@@ -1,3 +1,3 @@
-// $Id: tmatrix.cc,v 1.1 1999-04-29 08:57:15 ansari Exp $
+// $Id: tmatrix.cc,v 1.2 1999-04-30 11:02:52 ansari Exp $
 //                         C.Magneville          04/99
 #include "machdefs.h"
@@ -18,5 +18,5 @@
 TMatrix<T>::TMatrix()
 // Constructeur par defaut.
-: mNr(0), mNc(0), mNDBlock(NULL)
+: mNr(0), mNc(0), mNDBlock()
 {
 }
@@ -25,5 +25,5 @@
 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))
+: mNr(r), mNc(c), mNDBlock(r*c)
 {
 }
@@ -33,5 +33,5 @@
 // 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))
+: mNr(r), mNc(c), mNDBlock(r*c,values,br)
 {
 }
@@ -39,7 +39,6 @@
 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())))
+// Constructeur par copie (partage des donnees).
+: mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock)
 {
 }
@@ -48,6 +47,5 @@
 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))
+: mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock,share)
 {
 }
@@ -57,13 +55,4 @@
 // Destructeur
 {
-Delete();
-}
-
-template <class T>
-void TMatrix<T>::Delete(void)
-// Pour desallouer
-{
-mNr = mNc = 0;
-if(mNDBlock) delete mNDBlock; mNDBlock = NULL;
 }
 
@@ -71,11 +60,8 @@
 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()));
+// 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;
 }
 
@@ -84,9 +70,5 @@
 // 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);
+mNDBlock.Reset(v);
 }
 
@@ -96,8 +78,6 @@
 {
 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);
+mNr = r; mNc = c; mNDBlock.ReSize(r*c,force_alloc);
 }
 
@@ -109,7 +89,8 @@
 // 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;
+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;
 }
@@ -117,12 +98,10 @@
 template <class T>
 TMatrix<T>& TMatrix<T>::operator = (const TMatrix<T>& a)
-// Operateur d'affectation A=B (sans reaffectation)
+// Operateur d'affectation A=a
+// (Re)allocation de la place memoire
+// sauf si "a" est cree temporairement (partage)
 {
 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());
+Clone(a);
 return *this;
 }
@@ -138,12 +117,12 @@
   os<<"TMatrix::Print("<<mNr<<","<<mNc<<")"<<endl;
 if(lp>0)
-  {os<<"       this="<<this<<" mNDBlock="<<mNDBlock<<endl;
-   if(mNDBlock) mNDBlock->Print(0,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(int i=i0;i<i1;i++) {
-  for(int j=j0;j<j1;j++) cout<<" "<<(*this)(i,j)<<endl;
+for(uint_4 i=i0;i<i1;i++) {
+  for(uint_4 j=j0;j<j1;j++) cout<<" "<<(*this)(i,j);
+  cout<<endl;
 }
 }
@@ -155,7 +134,5 @@
 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;
+mNDBlock += b;
 return *this;
 }
@@ -164,7 +141,5 @@
 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;
+mNDBlock -= b;
 return *this;
 }
@@ -173,7 +148,5 @@
 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;
+mNDBlock *= b;
 return *this;
 }
@@ -182,8 +155,5 @@
 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;
+mNDBlock /= b;
 return *this;
 }
@@ -195,7 +165,7 @@
 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);
+if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
+  throw(SzMismatchError("TMatrix::operator+=A size mismatch"));
+mNDBlock += a.mNDBlock;
 return *this;
 }
@@ -204,7 +174,7 @@
 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);
+if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
+  throw(SzMismatchError("TMatrix::operator-=A size mismatch"));
+mNDBlock -= a.mNDBlock;
 return *this;
 }
@@ -215,6 +185,6 @@
 {
 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"));
+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
@@ -238,7 +208,137 @@
 }
 
+////////////////////////////////////////////////////////////////
+//**** 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>
@@ -249,13 +349,20 @@
 #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> >
-*/
+#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)
@@ -272,3 +379,14 @@
 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
Index: trunk/SophyaLib/NTools/tmatrix.h
===================================================================
--- trunk/SophyaLib/NTools/tmatrix.h	(revision 285)
+++ trunk/SophyaLib/NTools/tmatrix.h	(revision 286)
@@ -24,4 +24,9 @@
   virtual ~TMatrix();
 
+  // Temporaire?
+  inline bool IsTemp(void) const {return mNDBlock.IsTemp();}
+  inline void SetTemp(bool temp=false) const
+                     {mNDBlock.SetTemp(temp);}
+
   // Gestion taille/Remplissage
   void Clone(const TMatrix<T>& a);
@@ -33,11 +38,11 @@
   inline int NCol()  const {return mNc;}
   T const& operator()(uint_4 r,uint_4 c) const
-                     {return *(mNDBlock->Begin()+r*mNc+c);}
+                     {return *(mNDBlock.Begin()+r*mNc+c);}
   T&       operator()(uint_4 r,uint_4 c)
-                     {return *(mNDBlock->Begin()+r*mNc+c);}
-  inline       T* Data()       {return mNDBlock->Begin();}
-  inline const T* Data() const {return mNDBlock->Begin();}
-  inline       NDataBlock<T>* DataBlock()       {return mNDBlock;}
-  inline const NDataBlock<T>* DataBlock() const {return mNDBlock;}
+                     {return *(mNDBlock.Begin()+r*mNc+c);}
+  inline       T* Data()       {return mNDBlock.Begin();}
+  inline const T* Data() const {return mNDBlock.Begin();}
+  inline       NDataBlock<T>& DataBlock()       {return mNDBlock;}
+  inline const NDataBlock<T>& DataBlock() const {return mNDBlock;}
 
   // Operateur d'affectation
@@ -50,5 +55,5 @@
               {Print(cout,lp,i0,ni,j0,nj);}
 
-  // Surcharge d'operateurs
+  // Surcharge d'operateurs INPLACE
   TMatrix<T>& operator += (T b);
   TMatrix<T>& operator -= (T b);
@@ -60,15 +65,104 @@
   TMatrix<T>& operator *= (const TMatrix<T>& a);
 
+  // Pour surcharge d'operateurs
+  TMatrix<T> Add(const TMatrix<T>& b) const;
+  TMatrix<T> Sub(const TMatrix<T>& b) const;
+  TMatrix<T> Mul(const TMatrix<T>& b) const;
+
 protected:
-  void Delete(void);
-
   uint_4 mNr,mNc;
-  NDataBlock<T>* mNDBlock;
+  NDataBlock<T> mNDBlock;
 };
 
+/////////////////////////////////////////////////////////////////////////
 template<class T>
 inline ostream& operator << (ostream& os, const TMatrix<T>& a)
                       {a.Print(os); return(os);}
 
+/////////////////////////////////////////////////////////////////////////
+template<class T>
+inline TMatrix<T> operator + (const TMatrix<T>& a, double b) {
+  TMatrix<T> result(a,false); result.DataBlock().SetTemp(true);
+  result += b;
+  return result;
+}
+
+template<class T>
+inline TMatrix<T> operator + (double b,const TMatrix<T>& a) {
+  TMatrix<T> result(a,false); result.DataBlock().SetTemp(true);
+  result += b;
+  return result;
+}
+
+template<class T>
+inline TMatrix<T> operator - (const TMatrix<T>& a, double b) {
+  TMatrix<T> result(a,false); result.DataBlock().SetTemp(true);
+  result -= b;
+  return result;
+}
+
+template<class T>
+inline TMatrix<T> operator - (double b,const TMatrix<T>& a) {
+  TMatrix<T> result(a,false); result.DataBlock().SetTemp(true);
+  result.DataBlock() = b - result.DataBlock();
+  return result;
+}
+
+template<class T>
+inline TMatrix<T> operator * (const TMatrix<T>& a, double b) {
+  TMatrix<T> result(a,false); result.DataBlock().SetTemp(true);
+  result *= b;
+  return result;
+}
+
+template<class T>
+inline TMatrix<T> operator * (double b,const TMatrix<T>& a) {
+  TMatrix<T> result(a,false); result.DataBlock().SetTemp(true);
+  result *= b;
+  return result;
+}
+
+template<class T>
+inline TMatrix<T> operator / (const TMatrix<T>& a, double b) {
+  TMatrix<T> result(a,false); result.DataBlock().SetTemp(true);
+  result /= b;
+  return result;
+}
+
+/////////////////////////////////////////////////////////////////////////
+template<class T>
+inline TMatrix<T> operator + (const TMatrix<T>& a,const TMatrix<T>& b)
+                             {return a.Add(b);}
+
+template<class T>
+inline TMatrix<T> operator - (const TMatrix<T>& a,const TMatrix<T>& b)
+                             {return a.Sub(b);}
+
+template<class T>
+inline TMatrix<T> operator * (const TMatrix<T>& a,const TMatrix<T>& b)
+                             {return a.Mul(b);}
+
+/////////////////////////////////////////////////////////////////////////
+// Classe pour la gestion de persistance
+template <class T>
+class FIO_TMatrix : public  PPersist  {
+
+public:
+	    FIO_TMatrix();
+	    FIO_TMatrix(string const & filename); 
+	    FIO_TMatrix(const TMatrix<T> & obj);
+	    FIO_TMatrix(TMatrix<T> * obj);
+  virtual   ~FIO_TMatrix();
+
+  virtual   AnyDataObj* DataObj();
+  inline operator TMatrix<T>() { return(*dobj); }
+
+protected :
+  virtual void       ReadSelf(PInPersist&);           
+  virtual void       WriteSelf(POutPersist&) const;  
+  TMatrix<T> * dobj;
+  bool ownobj;
+};
+
 } // Fin du namespace
 
