Index: trunk/SophyaLib/TArray/array.h
===================================================================
--- trunk/SophyaLib/TArray/array.h	(revision 804)
+++ trunk/SophyaLib/TArray/array.h	(revision 804)
@@ -0,0 +1,15 @@
+// This may look like C code, but it is really -*- C++ -*-
+//  All the necessari include files for array manipulation
+//                     R. Ansari, C.Magneville   03/2000
+
+#ifndef AllTArray_SEEN
+#define AllTArray_SEEN
+
+#include "tarray.h"
+#include "tmatrix.h"
+#include "tvector.h"
+#include "matharr.h"
+#include "fioarr.h"
+#include "sopemtx.h"
+
+#endif
Index: trunk/SophyaLib/TArray/basarr.cc
===================================================================
--- trunk/SophyaLib/TArray/basarr.cc	(revision 803)
+++ trunk/SophyaLib/TArray/basarr.cc	(revision 804)
@@ -14,4 +14,5 @@
 				 "Stride(int )", "ElemCheckBound()", "operator()" };
 uint_4 BaseArray::max_nprt_ = 50;
+short  BaseArray::default_memory_mapping = CMemoryMapping;
 
 // Methodes statiques globales
@@ -22,5 +23,5 @@
 
 
-uint_8 BaseArray::ComputeTotalSize(uint_4 ndim, uint_4 * siz, uint_4 step, uint_8 offset)
+uint_8 BaseArray::ComputeTotalSize(uint_4 ndim, const uint_4 * siz, uint_4 step, uint_8 offset)
 {
   uint_8 rs = step;
@@ -28,4 +29,46 @@
   return(rs+offset);
 }
+
+short BaseArray::SetDefaultMemoryMapping(short mm)
+{
+  default_memory_mapping = ( (mm == CMemoryMapping) ? CMemoryMapping : FortranMemoryMapping) ; 
+  return default_memory_mapping;
+}
+
+
+short BaseArray::SelectMemoryMapping(short mm)
+// si mm == CMemoryMapping, elements d'une ligne suivant X (consecutif)
+// sinon (mm == FortranMemoryMapping)  elements d'une colonne suivant X 
+{
+  if ( (mm == CMemoryMapping) || (mm == FortranMemoryMapping) )  return (mm) ; 
+  else return (default_memory_mapping);
+}
+
+void BaseArray::UpdateMemoryMapping(short mm, uint_4 & nx, uint_4 & ny)
+{
+  if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) mm = default_memory_mapping;
+  if (mm == CMemoryMapping) {
+    marowi_ = veceli_ = 1;  macoli_ = 0;
+    uint_4 nr = nx;
+    nx = ny;   ny = nr;
+  }
+  else {
+    marowi_ = veceli_ = 0;  macoli_ = 1;
+  }  
+}
+
+void BaseArray::UpdateMemoryMapping(BaseArray const & a, short mm)
+{
+  if (mm == SameMemoryMapping) 
+    mm = ((a.marowi_ == 1) ? CMemoryMapping : FortranMemoryMapping);
+  if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) mm = default_memory_mapping;
+  if (mm == CMemoryMapping) {
+    marowi_ = veceli_ = 1;  macoli_ = 0;
+  }
+  else {
+    marowi_ = veceli_ = 0;  macoli_ = 1;
+  }  
+}
+
 
 // -------------------------------------------------------
@@ -44,7 +87,7 @@
   moystep_ = 0;
   offset_ = 0;
-  // Default for matrices : Lines = Y , Columns = X
+  // Default for matrices : Lines = Y , Columns = X, Default Vectore= ColumnVector
   macoli_ = 0;
-  marowi_ = 1;
+  veceli_ = marowi_ = 1;
 }
 
@@ -60,5 +103,6 @@
   for(int k=0; k<ndim_; k++) 
     if (size_[k] != a.size_[k])  return(false);
-  if ( (macoli_ != a.macoli_) || (marowi_ != a.marowi_) )  return(false);
+  if ( (macoli_ != a.macoli_) || (marowi_ != a.marowi_) ||
+       (veceli_ != a.veceli_) )  return(false);
   return(true);
 }
@@ -196,5 +240,5 @@
   // Default for matrices : Lines = Y , Columns = X
   macoli_ = 0;
-  marowi_ = 1;
+  marowi_ = veceli_ = 1;
   // Update OK 
   ndim_ = ndim;
@@ -238,5 +282,5 @@
   // Default for matrices : Lines = Y , Columns = X
   macoli_ = 0;
-  marowi_ = 1;
+  marowi_ = veceli_ = 1;
   // Update OK 
   ndim_ = ndim;
@@ -278,4 +322,5 @@
   macoli_ = a.macoli_;
   marowi_ = a.marowi_;
+  veceli_ = a.veceli_;
   // Update OK 
   ndim_ = a.ndim_;
@@ -284,11 +329,12 @@
 
 
-void BaseArray::SubArray(BaseArray & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step)
-{
-  if ( (ndim > ndim_) || (ndim < 1) ) throw(SzMismatchError("BaseArray::SubArray( ... ) NDim Error") );
+void BaseArray::UpdateSubArraySizes(BaseArray & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step) const
+{
+  if ( (ndim > ndim_) || (ndim < 1) ) 
+    throw(SzMismatchError("BaseArray::UpdateSubArraySizes( ... ) NDim Error") );
   int k;
   for(k=0; k<ndim; k++) 
     if ( (siz[k]*step[k]+pos[k]) > size_[k] )  
-      throw(SzMismatchError("BaseArray::SubArray( ... ) Size/Pos Error") );
+      throw(SzMismatchError("BaseArray::UpdateSubArraySizes( ... ) Size/Pos Error") );
   uint_8 offset = offset_;
   for(k=0; k<ndim_; k++) { 
@@ -296,8 +342,7 @@
     step[k] *= step_[k];
   }
-  string exm = "BaseArray::SubArray() ";
+  string exm = "BaseArray::UpdateSubArraySizes() ";
   if (!ra.UpdateSizes(ndim, siz, step, offset,  exm))
      throw( ParmError(exm) );
-  //  (ra.mNDBlock).Share(mNDBlock);
   return;
 }
Index: trunk/SophyaLib/TArray/basarr.h
===================================================================
--- trunk/SophyaLib/TArray/basarr.h	(revision 803)
+++ trunk/SophyaLib/TArray/basarr.h	(revision 804)
@@ -21,4 +21,14 @@
 class BaseArray : public AnyDataObj {
 public:
+  //  To define Array / Matrix memory mapping
+  enum { AutoMemoryMapping = -1, SameMemoryMapping = 0, 
+	 CMemoryMapping = 1, FortranMemoryMapping = 2 };
+
+  //    Max Nb of printed elements
+  static void    SetMaxPrint(uint_4 nprt=50);   
+  //   Memory organisation (for matrices)
+  static short   SetDefaultMemoryMapping(short mm=CMemoryMapping);
+  static inline short GetDefaultMemoryMapping() { return default_memory_mapping; }
+  
   // Creation / destruction 
   BaseArray();
@@ -44,4 +54,11 @@
 
   // memory organisation - packing information
+  inline short  GetMemoryMapping() const 
+                {return ( (marowi_ == 1) ? CMemoryMapping : FortranMemoryMapping) ; }
+
+  inline uint_4 RowsKA() const {return marowi_; }    // Dimension correspondant aux lignes
+  inline uint_4 ColsKA() const {return macoli_; }    // Dimension correspondant aux colonnes
+  inline uint_4 VectKA() const {return veceli_; }    // Dimension corr. aux elts d'un vecteur
+
   inline bool   IsPacked() const { return(moystep_ == 1); }
   inline bool   IsPackedX() const { return(step_[0] == 1); }
@@ -66,7 +83,6 @@
 
 // Impression, I/O, ...
-  static void   SetMaxPrint(uint_4 nprt=50);   
-  void          Show(ostream& os, bool si=false) const;
-  inline void   Show() const { Show(cout); }
+  void           Show(ostream& os, bool si=false) const;
+  inline void    Show() const { Show(cout); }
   virtual string DataType() const = 0;
 
@@ -83,7 +99,12 @@
   bool UpdateSizes(uint_4 ndim, const uint_4 * siz, const uint_4 * step, uint_8 offset, string & exmsg);
   bool UpdateSizes(const BaseArray& a, string & exmsg);
-  static uint_8 ComputeTotalSize(uint_4 ndim, uint_4 * siz, uint_4 step, uint_8 offset) ;
-  // Pour Extraction de sous-tableau
-  virtual void SubArray(BaseArray & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step);
+  static uint_8 ComputeTotalSize(uint_4 ndim, const uint_4 * siz, uint_4 step, uint_8 offset) ;
+  //  Organisation memoire
+  static short SelectMemoryMapping(short mm);
+  void UpdateMemoryMapping(short mm, uint_4 & nx, uint_4 & ny);
+  void UpdateMemoryMapping(BaseArray const & a, short mm);
+
+  // Pour Extraction de sous-tableau 
+  virtual void UpdateSubArraySizes(BaseArray & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step) const;
 
   uint_4 ndim_;                   // nb of dimensions
@@ -95,4 +116,5 @@
   uint_8 offset_;              // global offset -> position of elem[0] in DataBlock
   uint_4 marowi_, macoli_;     // For matrices, Row index and column index in dimensions
+  uint_4 veceli_;              // For vectors,  dimension index = marowi_/macoli_ (Row/Col vectors)
 
   DVList* mInfo;               // Infos (variables) attachees au tableau
@@ -100,4 +122,5 @@
   static char * ck_op_msg_[6];  // Operation messages for CheckDI() CheckBound()
   static uint_4 max_nprt_;      // Nb maxi d'elements imprimes
+  static short default_memory_mapping;  // Default memory mapping
 };
 
Index: trunk/SophyaLib/TArray/fioarr.cc
===================================================================
--- trunk/SophyaLib/TArray/fioarr.cc	(revision 803)
+++ trunk/SophyaLib/TArray/fioarr.cc	(revision 804)
@@ -3,5 +3,8 @@
 
 #include "pexceptions.h"
+#include "fiondblock.h"
 #include "fioarr.h"
+#include "tmatrix.h"
+#include "tvector.h"
 
 // --------------------------------------------------------
@@ -13,6 +16,6 @@
 FIO_TArray<T>::FIO_TArray()
 {
-  dobj=new TArray<T>;
-  ownobj=true;
+  dobj=NULL;
+  ownobj=false;
 }
 
@@ -20,6 +23,6 @@
 FIO_TArray<T>::FIO_TArray(string const & filename) 
 {
-  dobj=new TArray<T>;
-  ownobj=true; 
+  dobj=NULL;
+  ownobj=false;
   Read(filename);
 }
@@ -28,5 +31,11 @@
 FIO_TArray<T>::FIO_TArray(const TArray<T> & obj) 
 { 
-  dobj = new TArray<T>(obj);
+  const TVector<T> * tv = dynamic_cast<const TVector<T> *>(&obj);
+  if (tv != NULL)  dobj = new TVector<T>(*tv, true);
+  else {
+  const TMatrix<T> * tm = dynamic_cast<const TMatrix<T> *>(&obj);
+  if (tm != NULL)  dobj = new TMatrix<T>(*tm, true);
+  else dobj = new TArray<T>(obj, true);
+  }
   ownobj=true; 
 }
@@ -63,9 +72,14 @@
 void FIO_TArray<T>::ReadSelf(PInPersist& is)
 {
-// On lit les 3 premiers uint_4
-//  0: Numero de version,   : NRows,  2 : NCol
+// On lit les 5 premiers uint_4
+//  0: Numero de version,  1 : Type (Array, matrix, Vector, ...)  2 != 0 , has Info
+//          1:Type = 0 TArray , 12=(4+8) TMatrix , 48=(16+32) TVector
   uint_4 itab[5];
   is.Get(itab,5);
-  if (dobj == NULL) dobj = new TArray<T>;
+  if (dobj == NULL) {
+    if (itab[1] == 12)   dobj = new TMatrix<T>;
+    else if (itab[1] == 48)   dobj = new TVector<T>;
+    else   dobj = new TArray<T>;
+  }
 // On lit les tailles, etc ...
   is.Get(dobj->ndim_);
@@ -78,4 +92,5 @@
   is.Get(dobj->marowi_);
   is.Get(dobj->macoli_);
+  is.Get(dobj->veceli_);
 // On lit le datablock
   is >> dobj->DataBlock();
@@ -88,12 +103,22 @@
 {
   if (dobj == NULL)   return;
-//  On ecrit 4 uint_4 .... 
+//  On ecrit 5 uint_4 .... 
 //  0: Numero de version,  1 : Type (Array, matrix, Vector, ...)  2 != 0 , has Info
+//          1:Type = 0 TArray , 12=(4+8) TMatrix , 48=(16+32) TVector
+  uint_4 typa = 0;
+  TVector<T> * tv = dynamic_cast<TVector<T> *>(dobj);
+  if (tv != NULL)  typa = 48;
+  else {
+  TMatrix<T> * tm = dynamic_cast<TMatrix<T> *>(dobj);
+  if (tm != NULL)  typa = 12;
+  else  typa = 0;
+  }
+
   uint_4 itab[5];
   itab[0] = 1;  // Numero de version a 1
-  itab[1] = 0;
+  itab[1] = typa;  // Real object type
   itab[2] = (dobj->mInfo != NULL) ? 1 : 0;
   itab[3] = itab[4] = 0;
-  os.Put(itab,3);
+  os.Put(itab,5);
 // On ecrit les tailles, etc ...
   os.Put(dobj->ndim_);
@@ -106,4 +131,5 @@
   os.Put(dobj->marowi_);
   os.Put(dobj->macoli_);
+  os.Put(dobj->veceli_);
 // On ecrit le datablock
   os << dobj->DataBlock();
@@ -117,11 +143,11 @@
 #ifdef __CXX_PRAGMA_TEMPLATES__
 // Instances des delegues FileIO (PPersist)
-#pragma define_template FIO_TArray<uint_1>
+// #pragma define_template FIO_TArray<uint_1>
 #pragma define_template FIO_TArray<uint_2>
-#pragma define_template FIO_TArray<int_2>
+// #pragma define_template FIO_TArray<int_2>
 #pragma define_template FIO_TArray<int_4>
 #pragma define_template FIO_TArray<int_8>
-#pragma define_template FIO_TArray<uint_4>
-#pragma define_template FIO_TArray<uint_8>
+// #pragma define_template FIO_TArray<uint_4>
+// #pragma define_template FIO_TArray<uint_8>
 #pragma define_template FIO_TArray<r_8>
 #pragma define_template FIO_TArray<r_4>
@@ -132,11 +158,11 @@
 #if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
 // Instances des delegues FileIO (PPersist)
-template class FIO_TArray<uint_1>;
+// template class FIO_TArray<uint_1>;
 template class FIO_TArray<uint_2>;
-template class FIO_TArray<int_2>;
-template class FIO_TArray<int_4>;
+// template class FIO_TArray<int_2>;
+// template class FIO_TArray<int_4>;
 template class FIO_TArray<int_8>;
-template class FIO_TArray<uint_4>;
-template class FIO_TArray<uint_8>;
+// template class FIO_TArray<uint_4>;
+// template class FIO_TArray<uint_8>;
 template class FIO_TArray<r_8>;
 template class FIO_TArray<r_4>;
Index: trunk/SophyaLib/TArray/fioarr.h
===================================================================
--- trunk/SophyaLib/TArray/fioarr.h	(revision 803)
+++ trunk/SophyaLib/TArray/fioarr.h	(revision 804)
@@ -9,4 +9,6 @@
 #include <iostream.h>
 #include "tarray.h"
+#include "tmatrix.h"
+#include "tvector.h"
 #include "ppersist.h"
 
@@ -39,4 +41,5 @@
 { FIO_TArray<T> fio(&obj);  fio.Read(is);  return(is); }
 
+
 } // Fin du namespace
 
Index: trunk/SophyaLib/TArray/matharr.cc
===================================================================
--- trunk/SophyaLib/TArray/matharr.cc	(revision 803)
+++ trunk/SophyaLib/TArray/matharr.cc	(revision 804)
@@ -13,4 +13,6 @@
 TArray<T>& MathArray<T>::ApplyFunctionInPlace(TArray<T> & a, Arr_DoubleFunctionOfX f)
 {
+  if (a.NbDimensions() < 1) 
+    throw RangeCheckError("MathArray<T>::ApplyFunctionInPlace(TArray<T> & a..) Not Allocated Array a !");
   T * pe;
   uint_8 j,k;
@@ -36,4 +38,6 @@
 TArray<T>& MathArray<T>::ApplyFunctionInPlace(TArray<T> & a, Arr_FloatFunctionOfX f)
 {
+  if (a.NbDimensions() < 1) 
+    throw RangeCheckError("MathArray<T>::ApplyFunctionInPlace(TArray<T> & a..) Not Allocated Array a !");
   T * pe;
   uint_8 j,k;
@@ -75,4 +79,42 @@
 }
 
+template <class T>
+double MathArray<T>::MeanSigma(TArray<T> const & a, double & mean, double & sig)
+{
+  if (a.NbDimensions() < 1) 
+    throw RangeCheckError("MathArray<T>::MeanSigma(TArray<T> const & a..) Not Allocated Array a !");
+  const T * pe;
+  uint_8 j,k;
+  mean=0.;
+  sig = 0.;
+  double valok;
+  if (a.AvgStep() > 0)   {  // regularly spaced elements
+    uint_8 step = a.AvgStep(); 
+    uint_8 maxx = a.Size()*step;
+    pe = a.Data();
+    for(k=0; k<maxx; k+=step )  { 
+      valok = (double) pe[k]; 
+      mean += valok;  sig += valok*valok;
+    }
+  }
+  else {    // Non regular data spacing ...
+    uint_4 ka = a.MaxSizeKA();
+    uint_8 step = a.Step(ka);
+    uint_8 gpas = a.Size(ka)*step;
+    for(j=0; j<a.Size(); j += a.Size(ka))  {
+      pe = a.DataBlock().Begin()+a.Offset(j);
+      for(k=0; k<gpas; k+=step) { 
+	valok = (double) pe[k]; 
+	mean += valok;  sig += valok*valok;
+      }
+    }
+  }
+  double dsz = (double)(a.Size());
+  mean /= dsz;
+  sig = sig/dsz - mean*mean;
+  if (sig >= 0.) sig = sqrt(sig);
+  return(mean);
+}
+
 ///////////////////////////////////////////////////////////////
 #ifdef __CXX_PRAGMA_TEMPLATES__
Index: trunk/SophyaLib/TArray/matharr.h
===================================================================
--- trunk/SophyaLib/TArray/matharr.h	(revision 803)
+++ trunk/SophyaLib/TArray/matharr.h	(revision 804)
@@ -25,4 +25,6 @@
   virtual TArray<T>  ApplyFunction(TArray<T> const & a, Arr_DoubleFunctionOfX f);
   virtual TArray<T>  ApplyFunction(TArray<T> const & a, Arr_FloatFunctionOfX f);
+  // Computing Mean-Sigma 
+  virtual double     MeanSigma(TArray<T> const & a, double & mean, double & sig);
 };
 
@@ -41,4 +43,37 @@
   { MathArray<T> ma;   return( ma.ApplyFunction(a, tan) ); }
 
+template <class T>
+inline TArray<T> asin(const TArray<T>& a)
+  { MathArray<T> ma;   return( ma.ApplyFunction(a, asin) ); }
+
+template <class T>
+inline TArray<T> acos(const TArray<T>& a)
+  { MathArray<T> ma;   return( ma.ApplyFunction(a, acos) ); }
+
+template <class T>
+inline TArray<T> atan(const TArray<T>& a)
+  { MathArray<T> ma;   return( ma.ApplyFunction(a, atan) ); }
+
+template <class T>
+inline TArray<T> exp(const TArray<T>& a)
+  { MathArray<T> ma;   return( ma.ApplyFunction(a, exp) ); }
+
+template <class T>
+inline TArray<T> log(const TArray<T>& a)
+  { MathArray<T> ma;   return( ma.ApplyFunction(a, log) ); }
+
+template <class T>
+inline TArray<T> log10(const TArray<T>& a)
+  { MathArray<T> ma;   return( ma.ApplyFunction(a, log10) ); }
+
+template <class T>
+inline double MeanSigma(const TArray<T>& a, double & mean, double & sig)
+  { MathArray<T> ma;   return( ma.MeanSigma(a, mean, sig) ); }
+
+template <class T>
+inline double Mean(const TArray<T>& a)
+  { MathArray<T> ma;  double mean, sig;  return( ma.MeanSigma(a, mean, sig) ); }
+
+
 } // Fin du namespace
 
Index: trunk/SophyaLib/TArray/sopemtx.cc
===================================================================
--- trunk/SophyaLib/TArray/sopemtx.cc	(revision 803)
+++ trunk/SophyaLib/TArray/sopemtx.cc	(revision 804)
@@ -14,4 +14,123 @@
 // -------------------------------------------------------------
 ////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////
+// Classe de lignes/colonnes de matrices
+enum TRCKind {TmatrixRow=0, TmatrixCol=1, TmatrixDiag=2};
+template <class T>
+class TMatrixRC {
+public:
+  TMatrixRC();
+  TMatrixRC(TMatrix<T>& m, TRCKind kind, uint_4 index=0);
+  virtual ~TMatrixRC() {}
+
+  // Acces aux rangees et colonnes de matrices
+  static TMatrixRC<T> Row(TMatrix<T> & m, uint_4 r);
+  static TMatrixRC<T> Col(TMatrix<T> & m, uint_4 c);
+  static TMatrixRC<T> Diag(TMatrix<T> & m);
+
+  // ---- A virer   $CHECK$ Reza 03/2000
+  //  int_4 Next();
+  //  int_4 Prev();
+  int_4 SetCol(int_4 c);
+  int_4 SetRow(int_4 r);
+  int_4 SetDiag();
+
+  static uint_4 Step(const TMatrix<T>& m, TRCKind rckind);
+  static T* Org(const TMatrix<T>&, TRCKind rckind, uint_4 ind=0);
+
+  TRCKind Kind() const { return kind; }
+  uint_4 NElts() const;
+  T& operator()(uint_4 i);
+  T  operator()(uint_4 i) const;
+
+
+   TMatrixRC<T>& operator = (const TMatrixRC<T>& rc);
+
+// ---- A virer   $CHECK$ Reza 03/2000
+//   TVector<T> GetVect() const;
+//   TMatrixRC<T>& operator += (const TMatrixRC<T>& rc);
+//   TMatrixRC<T>& operator -= (const TMatrixRC<T>& rc);
+
+  TMatrixRC<T>& operator *= (T x);
+  TMatrixRC<T>& operator /= (T x);
+
+//   TMatrixRC<T>& operator -= (T x);
+//   TMatrixRC<T>& operator += (T x);
+  
+
+  TMatrixRC<T>& LinComb(T a, T b, const TMatrixRC& rc, uint_4 first=0);
+  TMatrixRC<T>& LinComb(T b, const TMatrixRC<T>& rc, uint_4 first=0);
+
+  uint_4 IMaxAbs(uint_4 first=0);
+
+  static void Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2);
+
+  inline static double Abs_Value(uint_1 v) {return (double) v;}
+  inline static double Abs_Value(uint_2 v) {return (double) v;}
+  inline static double Abs_Value(int_2 v)  {return (v>0)? (double) v: (double) -v;}
+  inline static double Abs_Value(int_4 v)  {return (v>0)? (double) v: (double) -v;}
+  inline static double Abs_Value(int_8 v)  {return (v>0)? (double) v: (double) -v;}
+  inline static double Abs_Value(uint_4 v) {return (double) v;}
+  inline static double Abs_Value(uint_8 v) {return (double) v;}
+  inline static double Abs_Value(r_4 v)    {return (double) fabsf(v);}
+  inline static double Abs_Value(r_8 v)    {return fabs(v);}
+  inline static double Abs_Value(complex<float> v)
+                {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
+  inline static double Abs_Value(complex<double> v)
+                {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
+
+protected:
+  TMatrix<T>* matrix;
+  T*          data;
+  int_4       index;
+  uint_4      step;
+  TRCKind     kind;
+};
+
+
+
+template <class T>
+inline T operator * (const TMatrixRC<T>& a, const TMatrixRC<T>& b)
+  {
+  if ( a.NElts() != b.NElts() )
+    throw(SzMismatchError("TMatrixRC::operator * size mismatch\n"));
+  if ( a.Kind() != b.Kind() )
+    throw(SzMismatchError("TMatrixRC::operator * type mismatch\n"));
+  T sum = 0;
+  for(uint_4 i=0; i<a.NElts(); i++) sum += a(i)*b(i);
+  return sum;
+  }
+
+
+template <class T>
+inline uint_4 TMatrixRC<T>::Step(const TMatrix<T>& m, TRCKind rckind)
+  { switch (rckind) { case TmatrixRow  : return m.Step(m.RowsKA());
+                      case TmatrixCol  : return m.Step(m.ColsKA());
+                      case TmatrixDiag : return (m.Step(m.RowsKA())+m.Step(m.ColsKA())); }
+    return 0; }
+
+template <class T>
+inline T* TMatrixRC<T>::Org(const TMatrix<T>& m, TRCKind rckind, uint_4 index)
+  { switch (rckind) { case TmatrixRow  : return const_cast<T *>(&(m(index,0)));
+                      case TmatrixCol  : return const_cast<T *>(&(m(0,index)));
+                      case TmatrixDiag : return const_cast<T *>(m.Data()); }
+    return NULL; }
+
+template <class T> inline uint_4 TMatrixRC<T>::NElts() const
+  { if (!matrix) return 0;
+    switch (kind) { case TmatrixRow  : return matrix->NCols();
+                    case TmatrixCol  : return matrix->NRows();
+                    case TmatrixDiag : return matrix->NCols(); }
+    return 0; }
+
+template <class T>
+inline T& TMatrixRC<T>::operator()(uint_4 i) {return data[i*step];}
+template <class T>
+inline T  TMatrixRC<T>::operator()(uint_4 i) const {return data[i*step];}
+
+////////////////////////////////////////////////////////////////
+// Typedef pour simplifier et compatibilite Peida
+typedef TMatrixRC<r_8> MatrixRC;
+
 
 template <class T> TMatrixRC<T>::TMatrixRC()
@@ -52,28 +171,30 @@
 
 
-template <class T> int_4 TMatrixRC<T>::Next()
-{
-if (!matrix || kind==TmatrixDiag) return -1;
-index++;
-if(kind == TmatrixRow) {
-  if(index > (int_4)matrix->NRows()) {index = (int_4)matrix->NRows(); return -1;}
-  data += matrix->NCols();
-} else {
-  if (index > (int_4)matrix->NCols()) {index = (int_4)matrix->NCols(); return -1;}
-  data++;
-}
-return index;
-}
-
-template <class T> int_4 TMatrixRC<T>::Prev()
-{
-if (!matrix || kind == TmatrixDiag) return -1;
-index--;
-if(index < 0) {index = 0; return -1;}
-if(kind == TmatrixRow) data -= matrix->NCols();
-else data--;
-return index;
-}
-
+// ---- A virer   $CHECK$ Reza 03/2000
+// template <class T> int_4 TMatrixRC<T>::Next()
+// {
+// if (!matrix || kind==TmatrixDiag) return -1;
+// index++;
+// if(kind == TmatrixRow) {
+//   if(index > (int_4)matrix->NRows()) {index = (int_4)matrix->NRows(); return -1;}
+//   data += matrix->NCols();
+// } else {
+//   if (index > (int_4)matrix->NCols()) {index = (int_4)matrix->NCols(); return -1;}
+//   data++;
+// }
+// return index;
+// }
+
+// template <class T> int_4 TMatrixRC<T>::Prev()
+// {
+// if (!matrix || kind == TmatrixDiag) return -1;
+// index--;
+// if(index < 0) {index = 0; return -1;}
+// if(kind == TmatrixRow) data -= matrix->NCols();
+// else data--;
+// return index;
+// }
+
+ 
 template <class T> int_4 TMatrixRC<T>::SetCol(int_4 c)
 {
@@ -121,30 +242,31 @@
 }
 
-template <class T> TVector<T> TMatrixRC<T>::GetVect() const
-{
-TVector<T> v(NElts());
-for (uint_4 i=0; i<NElts(); i++) v(i) = (*this)(i);
-return v;
-}
-
-template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (const TMatrixRC<T>& rc)
-{
-if ( NElts() != rc.NElts() )
-  throw(SzMismatchError("TMatrixRC::operator+= size mismatch\n"));
-if ( kind != rc.kind )
-  throw(SzMismatchError("TMatrixRC::operator+= type mismatch\n"));
-for (uint_4 i=0; i<NElts(); i++) (*this)(i) += rc(i);
-return *this;
-}
-
-template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (const TMatrixRC<T>& rc)
-{
-if( NElts() != rc.NElts() )
-  throw(SzMismatchError("TMatrixRC::operator-= size mismatch\n"));
-if( kind != rc.kind )
-  throw(SzMismatchError("TMatrixRC::operator-= type mismatch\n"));
-for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= rc(i);
-return *this;
-}
+// ---- A virer   $CHECK$ Reza 03/2000
+// template <class T> TVector<T> TMatrixRC<T>::GetVect() const
+// {
+// TVector<T> v(NElts());
+// for (uint_4 i=0; i<NElts(); i++) v(i) = (*this)(i);
+// return v;
+// }
+
+// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (const TMatrixRC<T>& rc)
+// {
+// if ( NElts() != rc.NElts() )
+//   throw(SzMismatchError("TMatrixRC::operator+= size mismatch\n"));
+// if ( kind != rc.kind )
+//   throw(SzMismatchError("TMatrixRC::operator+= type mismatch\n"));
+// for (uint_4 i=0; i<NElts(); i++) (*this)(i) += rc(i);
+// return *this;
+// }
+
+// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (const TMatrixRC<T>& rc)
+// {
+// if( NElts() != rc.NElts() )
+//   throw(SzMismatchError("TMatrixRC::operator-= size mismatch\n"));
+// if( kind != rc.kind )
+//   throw(SzMismatchError("TMatrixRC::operator-= type mismatch\n"));
+// for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= rc(i);
+// return *this;
+// }
 
 
@@ -162,15 +284,16 @@
 
 
-template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (T x)
-{
-for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= x;
-return *this;
-}
-
-template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (T x)
-{
-for(uint_4 i=0; i<NElts(); i++) (*this)(i) += x;
-return *this;
-}
+// ---- A virer   $CHECK$ Reza 03/2000
+// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (T x)
+// {
+// for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= x;
+// return *this;
+// }
+
+// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (T x)
+// {
+// for(uint_4 i=0; i<NElts(); i++) (*this)(i) += x;
+// return *this;
+// }
 
 template <class T>
@@ -331,20 +454,138 @@
 
 
+LinFitter::LinFitter()
+{
+}
+
+LinFitter::~LinFitter()
+{
+}
+
+double LinFitter::LinFit(const Vector& x, const Vector& y, int nf, 
+                         double (*f)(int, double), Vector& c)
+{
+  int n = x.NElts();
+  if (n != y.NElts()) THROW(sizeMismatchErr);
+  
+  Matrix fx(nf, n);
+  for (int i=0; i<nf; i++)
+    for (int j=0; j<n; j++)
+      fx(i,j) = f(i,x(j));
+
+  return LinFit(fx,y,c);
+}
+
+
+
+double LinFitter::LinFit(const Matrix& fx, const Vector& y, Vector& c)
+{
+  int n = y.NElts();
+  if ( n != fx.NCol()) THROW(sizeMismatchErr);
+
+  int nf = fx.NRows();
+
+  Matrix a(nf,nf);
+
+  for (int j=0; j<nf; j++)
+    for (int k=j; k<nf; k++)
+      a(j,k) = a(k,j) = TMatrixRC<r_8>::Row(const_cast<Matrix &>(fx), j) 
+
+	* TMatrixRC<r_8>::Row(const_cast<Matrix &>(fx), k); /* $CHECK$ Reza 10/3/2000 */
+   
+  c = fx * y;
+
+  if (SimpleMatrixOperation<r_8>::GausPiv(a,c) == 0) THROW(singMatxErr); /* $CHECK$ Reza 10/3/2000 */
+
+  double xi2 = 0;
+
+  for (int k=0; k<n; k++) {
+    double x = 0;
+    for (int i=0; i<nf; i++)
+      x += c(i) * fx(i,k);
+    x -= y(k);
+    xi2 += x*x;
+  }
+  return xi2;
+}
+
+
+
+double LinFitter::LinFit(const Vector& x, const Vector& y, const Vector& errY2, int nf,
+	                 double (*f)(int, double), Vector& c, Vector& errC)
+{
+  int n = x.NElts();
+  if (n != y.NElts()) THROW(sizeMismatchErr);
+  
+  Matrix fx(nf, n);
+  for (int i=0; i<nf; i++)
+    for (int j=0; j<n; j++)
+      fx(i,j) = f(i,x(j));
+
+  return LinFit(fx,y,errY2,c,errC);
+}
+
+
+double LinFitter::LinFit(const Matrix& fx, const Vector& y, const Vector& errY2,
+	   Vector& c, Vector& errC)
+{
+  int n = y.NElts();
+  if ( n != errY2.NElts()) THROW(sizeMismatchErr);
+  if ( n != fx.NCol()) THROW(sizeMismatchErr);
+
+  int nf = fx.NRows();
+
+  Matrix a(nf,nf);
+
+  c.Realloc(nf);
+  errC.Realloc(nf);
+
+  for (int j=0; j<nf; j++)
+    for (int k=j; k<nf; k++) {
+      double x=0;
+      for (int l=0; l<n; l++)
+	x += fx(j,l) * fx(k,l) / errY2(l);        // Matrice a inverser
+      a(j,k) = a(k,j) = x;
+    }
+ 
+  Matrix d(nf,nf+1);
+  for (int k=0; k<nf; k++) {
+    double x=0;
+    for (int l=0; l<n; l++)
+      x += fx(k,l) * y(l) / errY2(l);             // Second membre 1ere colonne
+    d(k,0) = x;
+    for (int m=1; m<=nf; m++)
+      d(k,m) = (k==m) ? 1 : 0;                // Reste second membre = Id.
+  }
+
+  if (SimpleMatrixOperation<r_8>::GausPiv(a,d) == 0) THROW(singMatxErr); /* $CHECK$ Reza 10/3/2000 */
+
+  for (int l=0; l<nf; l++) {
+    c(l) = d(l,0);                            // Parametres = 1ere colonne
+    errC(l) = d(l,l+1);                        // Erreurs = diag inverse.
+  }
+
+  double xi2 = 0;
+
+  for (int jj=0; jj<n; jj++) {
+    double x = 0;
+    for (int ii=0; ii<nf; ii++)
+      x += c(ii) * fx(ii,jj);
+    x -= y(jj);
+    xi2 += x*x/errY2(jj);
+  }
+  return xi2;
+}
+
+
+
 
 ///////////////////////////////////////////////////////////////
 #ifdef __CXX_PRAGMA_TEMPLATES__
 // Instances gestion lignes/colonnes
-#pragma define_template TMatrixRC<uint_1>
-#pragma define_template TMatrixRC<uint_2>
-#pragma define_template TMatrixRC<int_2>
 #pragma define_template TMatrixRC<int_4>
-#pragma define_template TMatrixRC<int_8>
-#pragma define_template TMatrixRC<uint_4>
-#pragma define_template TMatrixRC<uint_8>
 #pragma define_template TMatrixRC<r_4>
 #pragma define_template TMatrixRC<r_8>
-#pragma define_template TMatrixRC< complex<r_4> > 
-#pragma define_template TMatrixRC< complex<r_8> >
-// #pragma define_template SimpleMatrixOperation<r_4>
+#pragma define_template SimpleMatrixOperation<int_4>
+#pragma define_template SimpleMatrixOperation<r_4>
 #pragma define_template SimpleMatrixOperation<r_8>
 #endif
@@ -352,15 +593,8 @@
 #if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
 // Instances gestion lignes/colonnes
-template class TMatrixRC<uint_1>;
-template class TMatrixRC<uint_2>;
-template class TMatrixRC<int_2>;
 template class TMatrixRC<int_4>;
-template class TMatrixRC<int_8>;
-template class TMatrixRC<uint_4>;
-template class TMatrixRC<uint_8>;
 template class TMatrixRC<r_4>;
 template class TMatrixRC<r_8>;
-template class TMatrixRC< complex<r_4> >;
-template class TMatrixRC< complex<r_8> >;
+template class SimpleMatrixOperation<int_4>;
 template class SimpleMatrixOperation<r_4>;
 template class SimpleMatrixOperation<r_8>;
Index: trunk/SophyaLib/TArray/sopemtx.h
===================================================================
--- trunk/SophyaLib/TArray/sopemtx.h	(revision 803)
+++ trunk/SophyaLib/TArray/sopemtx.h	(revision 804)
@@ -9,115 +9,5 @@
 namespace SOPHYA {
 
-/////////////////////////////////////////////////////////////////////////
-// Classe de lignes/colonnes de matrices
-enum TRCKind {TmatrixRow=0, TmatrixCol=1, TmatrixDiag=2};
-template <class T>
-class TMatrixRC {
-public:
-  TMatrixRC();
-  TMatrixRC(TMatrix<T>& m, TRCKind kind, uint_4 index=0);
-  virtual ~TMatrixRC() {}
 
-  // Acces aux rangees et colonnes de matrices
-  static TMatrixRC<T> Row(TMatrix<T> & m, uint_4 r);
-  static TMatrixRC<T> Col(TMatrix<T> & m, uint_4 c);
-  static TMatrixRC<T> Diag(TMatrix<T> & m);
-
-  int_4 Next();
-  int_4 Prev();
-  int_4 SetCol(int_4 c);
-  int_4 SetRow(int_4 r);
-  int_4 SetDiag();
-
-  static uint_4 Step(const TMatrix<T>& m, TRCKind rckind);
-  static T* Org(const TMatrix<T>&, TRCKind rckind, uint_4 ind=0);
-
-  TRCKind Kind() const { return kind; }
-  uint_4 NElts() const;
-  T& operator()(uint_4 i);
-  T  operator()(uint_4 i) const;
-
-  TMatrixRC<T>& operator = (const TMatrixRC<T>& rc);
-  TVector<T> GetVect() const;
-
-  TMatrixRC<T>& operator += (const TMatrixRC<T>& rc);
-  TMatrixRC<T>& operator -= (const TMatrixRC<T>& rc);
-
-  TMatrixRC<T>& operator *= (T x);
-  TMatrixRC<T>& operator /= (T x);
-  TMatrixRC<T>& operator -= (T x);
-  TMatrixRC<T>& operator += (T x);
-
-  TMatrixRC<T>& LinComb(T a, T b, const TMatrixRC& rc, uint_4 first=0);
-  TMatrixRC<T>& LinComb(T b, const TMatrixRC<T>& rc, uint_4 first=0);
-
-  uint_4 IMaxAbs(uint_4 first=0);
-
-  static void Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2);
-
-  inline static double Abs_Value(uint_1 v) {return (double) v;}
-  inline static double Abs_Value(uint_2 v) {return (double) v;}
-  inline static double Abs_Value(int_2 v)  {return (v>0)? (double) v: (double) -v;}
-  inline static double Abs_Value(int_4 v)  {return (v>0)? (double) v: (double) -v;}
-  inline static double Abs_Value(int_8 v)  {return (v>0)? (double) v: (double) -v;}
-  inline static double Abs_Value(uint_4 v) {return (double) v;}
-  inline static double Abs_Value(uint_8 v) {return (double) v;}
-  inline static double Abs_Value(r_4 v)    {return (double) fabsf(v);}
-  inline static double Abs_Value(r_8 v)    {return fabs(v);}
-  inline static double Abs_Value(complex<float> v)
-                {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
-  inline static double Abs_Value(complex<double> v)
-                {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
-
-protected:
-  TMatrix<T>* matrix;
-  T*          data;
-  int_4       index;
-  uint_4      step;
-  TRCKind     kind;
-};
-
-
-template <class T>
-inline T operator * (const TMatrixRC<T>& a, const TMatrixRC<T>& b)
-  {
-  if ( a.NElts() != b.NElts() )
-    throw(SzMismatchError("TMatrixRC::operator * size mismatch\n"));
-  if ( a.Kind() != b.Kind() )
-    throw(SzMismatchError("TMatrixRC::operator * type mismatch\n"));
-  T sum = 0;
-  for(uint_4 i=0; i<a.NElts(); i++) sum += a(i)*b(i);
-  return sum;
-  }
-
-template <class T>
-inline uint_4 TMatrixRC<T>::Step(const TMatrix<T>& m, TRCKind rckind)
-  { switch (rckind) { case TmatrixRow  : return 1;
-                      case TmatrixCol  : return m.NCols();
-                      case TmatrixDiag : return m.NCols()+1; }
-    return 0; }
-
-template <class T>
-inline T* TMatrixRC<T>::Org(const TMatrix<T>& m, TRCKind rckind, uint_4 index)
-  { switch (rckind) { case TmatrixRow  : return const_cast<T *>(m.Data()) + index * m.NCols();
-                      case TmatrixCol  : return const_cast<T *>(m.Data()) + index;
-                      case TmatrixDiag : return const_cast<T *>(m.Data()); }
-    return NULL; }
-
-template <class T> inline uint_4 TMatrixRC<T>::NElts() const
-  { if (!matrix) return 0;
-    switch (kind) { case TmatrixRow  : return matrix->NCols();
-                    case TmatrixCol  : return matrix->NRows();
-                    case TmatrixDiag : return matrix->NCols(); }
-    return 0; }
-
-template <class T>
-inline T& TMatrixRC<T>::operator()(uint_4 i) {return data[i*step];}
-template <class T>
-inline T  TMatrixRC<T>::operator()(uint_4 i) const {return data[i*step];}
-
-////////////////////////////////////////////////////////////////
-// Typedef pour simplifier et compatibilite Peida
-typedef TMatrixRC<r_8> MatrixRC;
 
 ////////////////////////////////////////////////////////////////
@@ -149,4 +39,37 @@
 }
 
+//--------------------------------------
+//        Linear fitting 
+//--------------------------------------
+
+class LinFitter {
+public :
+		LinFitter();
+  virtual	~LinFitter();
+
+  double	LinFit(const Vector& x, const Vector& y, int nf,
+		       double (*f)(int, double), Vector& c);
+// fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1;
+
+  double	LinFit(const Matrix& fx, const Vector& y, Vector& c);
+// fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1,
+// la matrice fx contient les valeurs des f:
+// fx(i,j) = f(i, x(j)).
+		      
+  double	LinFit(const Vector& x, const Vector& y, const Vector& errY2, int nf,
+		       double (*f)(int, double), Vector& c, Vector& errC);
+// fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1,
+// errY2 contient les carres des erreurs sur les Y.
+// au retour, errC contient les erreurs sur les coefs.
+
+  double	LinFit(const Matrix& fx, const Vector& y, const Vector& errY2,
+		       Vector& c, Vector& errC);
+// fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1,
+// la matrice fx contient les valeurs des f:
+// fx(i,j) = f(i, x(j)).
+// errY2 contient les carres des erreurs sur les Y.
+// au retour, errC contient les erreurs sur les coefs.
+};
+
 
 } // Fin du namespace
Index: trunk/SophyaLib/TArray/tarray.cc
===================================================================
--- trunk/SophyaLib/TArray/tarray.cc	(revision 803)
+++ trunk/SophyaLib/TArray/tarray.cc	(revision 804)
@@ -18,5 +18,5 @@
 template <class T>
 TArray<T>::TArray()
-  : mNDBlock() , BaseArray()
+  : BaseArray() , mNDBlock() 
 // Default constructor
 {
@@ -24,6 +24,6 @@
 
 template <class T>
-TArray<T>::TArray(uint_4 ndim, uint_4 * siz, uint_4 step)
-  : mNDBlock(ComputeTotalSize(ndim, siz, step, 1)) , BaseArray()
+TArray<T>::TArray(uint_4 ndim, const uint_4 * siz, uint_4 step)
+  : BaseArray() , mNDBlock(ComputeTotalSize(ndim, siz, step, 1))
 {
   string exmsg = "TArray<T>::TArray(uint_4, uint_4 *, uint_4)";
@@ -32,11 +32,14 @@
 
 template <class T>
-TArray<T>::TArray(uint_4 nx, uint_4 ny=1, uint_4 nz=1)
-  : mNDBlock(nx*((ny>0)?ny:1)*((nz>0)?nz:1)) , BaseArray()
+TArray<T>::TArray(uint_4 nx, uint_4 ny, uint_4 nz, uint_4 nt, uint_4 nu)
+  : BaseArray() , mNDBlock(nx*((ny>0)?ny:1)*((nz>0)?nz:1)*((nt>0)?nt:1)*((nu>0)?nu:1)) 
 {
   uint_4 size[BASEARRAY_MAXNDIMS];
   size[0] = nx;   size[1] = ny;   size[2] = nz;
+  size[3] = nt;   size[4] = nu;
   int ndim = 1;
-  if ((size[1] > 0) && (size[2] > 0)) ndim = 3;
+  if ((size[1] > 0) && (size[2] > 0) && (size[3] > 0) && (size[4] > 0) ) ndim = 5;
+  else if ((size[1] > 0) && (size[2] > 0) && (size[3] > 0) ) ndim = 4;
+  else if ((size[1] > 0) && (size[2] > 0)) ndim = 3;
   else if (size[1] > 0)  ndim = 2;
   else ndim = 1;
@@ -46,6 +49,6 @@
 
 template <class T>
-TArray<T>::TArray(uint_4 ndim, uint_4 * siz, NDataBlock<T> & db, bool share, uint_4 step, uint_8 offset)
-  : mNDBlock(db, share) , BaseArray()
+TArray<T>::TArray(uint_4 ndim, const uint_4 * siz, NDataBlock<T> & db, bool share, uint_4 step, uint_8 offset)
+  : BaseArray() , mNDBlock(db, share) 
 {
   string exmsg = "TArray<T>::TArray(uint_4, uint_4 *,  NDataBlock<T> & ... )";
@@ -55,6 +58,6 @@
 
 template <class T>
-TArray<T>::TArray(uint_4 ndim, uint_4 * siz, T* values, uint_4 step, uint_8 offset, Bridge* br)
-  : mNDBlock(ComputeTotalSize(ndim, siz, step, 1), values, br) , BaseArray()
+TArray<T>::TArray(uint_4 ndim, const uint_4 * siz, T* values, uint_4 step, uint_8 offset, Bridge* br)
+  : BaseArray() , mNDBlock(ComputeTotalSize(ndim, siz, step, 1), values, br) 
 {
   string exmsg = "TArray<T>::TArray(uint_4, uint_4 *, T* ... )";
@@ -64,5 +67,5 @@
 template <class T>
 TArray<T>::TArray(const TArray<T>& a)
-  : mNDBlock(a.mNDBlock) , BaseArray()
+  : BaseArray() , mNDBlock(a.mNDBlock) 
 {
   string exmsg = "TArray<T>::TArray(const TArray<T>&)";
@@ -73,5 +76,5 @@
 template <class T>
 TArray<T>::TArray(const TArray<T>& a, bool share)
-  : mNDBlock(a.mNDBlock, share) , BaseArray()
+  : BaseArray() , mNDBlock(a.mNDBlock, share) 
 {
   string exmsg = "TArray<T>::TArray(const TArray<T>&, bool)";
@@ -87,5 +90,5 @@
 
 template <class T>
-TArray<T>& TArray<T>::operator = (const TArray<T>& a)
+TArray<T>& TArray<T>::Set(const TArray<T>& a)
 {
   if (this != &a) { 
@@ -171,8 +174,29 @@
 
 
+template <class T>
+TArray<T> TArray<T>::PackElements(bool force) const
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::PackElements() - Not Allocated Array ! ");
+  if ( !force && (AvgStep() == 1) ) {
+    TArray<T> ra(*this, true);
+    ra.SetTemp(true);
+    return(ra);
+  }
+  else {
+    TArray<T> ra(ndim_, size_, 1);
+    ra.CopyElt(*this);
+    ra.SetTemp(true);
+    return(ra);
+  }
+}
+
 // SubArrays 
-template <class T>
-TArray<T> TArray<T>::operator () (Range rx, Range ry, Range rz, Range rt, Range ru)
-{
+// $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
+template <class T>
+TArray<T> TArray<T>::SubArray(Range rx, Range ry, Range rz, Range rt, Range ru) const
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::operator () (Range, ...) - Not Allocated Array ! ");
   uint_4 ndim = 0;
   uint_4 size[BASEARRAY_MAXNDIMS];
@@ -199,5 +223,5 @@
   ndim = ndim_;
   TArray<T> ra;
-  SubArray(ra, ndim, size, pos, step); 
+  UpdateSubArraySizes(ra, ndim, size, pos, step); 
   ra.DataBlock().Share(this->DataBlock());
   ra.SetTemp(true);
@@ -210,6 +234,8 @@
 //  Possibilite de // , vectorisation 
 template <class T>
-TArray<T>& TArray<T>::operator = (Sequence seq)
-{
+TArray<T>& TArray<T>::Set(Sequence seq)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::Set(Sequence ) - Not Allocated Array ! ");
   T * pe;
   uint_8 j,k;
@@ -234,6 +260,8 @@
 
 template <class T>
-TArray<T>& TArray<T>::operator = (T x)
-{
+TArray<T>& TArray<T>::Set(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::Set(T )  - Not Allocated Array ! ");
   T * pe;
   uint_8 j,k;
@@ -257,6 +285,8 @@
 
 template <class T>
-TArray<T>& TArray<T>::operator += (T x)
-{
+TArray<T>& TArray<T>::Add(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::Add(T )  - Not Allocated Array ! ");
   T * pe;
   uint_8 j,k;
@@ -280,6 +310,8 @@
 
 template <class T>
-TArray<T>& TArray<T>::operator -= (T x)
-{
+TArray<T>& TArray<T>::Sub(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::Sub(T )  - Not Allocated Array ! ");
   T * pe;
   uint_8 j,k;
@@ -303,6 +335,8 @@
 
 template <class T>
-TArray<T>& TArray<T>::operator *= (T x)
-{
+TArray<T>& TArray<T>::Mul(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::Mul(T )  - Not Allocated Array ! ");
   T * pe;
   uint_8 j,k;
@@ -326,6 +360,8 @@
 
 template <class T>
-TArray<T>& TArray<T>::operator /= (T x)
-{
+TArray<T>& TArray<T>::Div(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::Div(T )  - Not Allocated Array ! ");
   T * pe;
   uint_8 j,k;
@@ -349,9 +385,62 @@
 
 
+template <class T>
+TArray<T>& TArray<T>::SubInv(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::SubInv(T )  - Not Allocated Array ! ");
+  T * pe;
+  uint_8 j,k;
+  if (AvgStep() > 0)   {  // regularly spaced elements
+    uint_8 step = AvgStep();
+    uint_8 maxx = totsize_*step;
+    pe = Data();
+    for(k=0; k<maxx; k+=step )  pe[k] = x-pe[k];
+  }
+  else {    // Non regular data spacing ...
+    uint_4 ka = MaxSizeKA();
+    uint_8 step = Step(ka);
+    uint_8 gpas = Size(ka)*step;
+    for(j=0; j<totsize_; j += Size(ka))  {
+      pe = mNDBlock.Begin()+Offset(j);
+      for(k=0; k<gpas; k+=step)  pe[k] = x-pe[k];
+    }
+  }
+  return(*this);
+}
+
+template <class T>
+TArray<T>& TArray<T>::DivInv(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::DivInv(T )  - Not Allocated Array ! ");
+  T * pe;
+  uint_8 j,k;
+  if (AvgStep() > 0)   {  // regularly spaced elements
+    uint_8 step = AvgStep();
+    uint_8 maxx = totsize_*step;
+    pe = Data();
+    for(k=0; k<maxx; k+=step )  pe[k] = x/pe[k];
+  }
+  else {    // Non regular data spacing ...
+    uint_4 ka = MaxSizeKA();
+    uint_8 step = Step(ka);
+    uint_8 gpas = Size(ka)*step;
+    for(j=0; j<totsize_; j += Size(ka))  {
+      pe = mNDBlock.Begin()+Offset(j);
+      for(k=0; k<gpas; k+=step)  pe[k] = x/pe[k];
+    }
+  }
+  return(*this);
+}
+
+
 //  >>>> Operations avec 2nd membre de type tableau
 template <class T>
-TArray<T>& TArray<T>::operator += (const TArray<T>& a)
-{
-  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::operator += (const TArray<T>&)")) ;
+TArray<T>& TArray<T>::AddElt(const TArray<T>& a)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::AddElt(const TArray<T>& )  - Not Allocated Array ! ");
+  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::AddElt(const TArray<T>&)")) ;
 
   T * pe;
@@ -381,7 +470,9 @@
 
 template <class T>
-TArray<T>& TArray<T>::operator -= (const TArray<T>& a)
-{
-  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::operator -= (const TArray<T>&)")) ;
+TArray<T>& TArray<T>::SubElt(const TArray<T>& a)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::SubElt(const TArray<T>& )  - Not Allocated Array ! ");
+  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::SubElt(const TArray<T>&)")) ;
 
   T * pe;
@@ -412,7 +503,9 @@
 
 template <class T>
-TArray<T>& TArray<T>::MultElt(const TArray<T>& a)
-{
-  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::MultElt(const TArray<T>&)")) ;
+TArray<T>& TArray<T>::MulElt(const TArray<T>& a)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::MulElt(const TArray<T>& )  - Not Allocated Array ! ");
+  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::MulElt(const TArray<T>&)")) ;
 
   T * pe;
@@ -441,7 +534,10 @@
 }
 
+
 template <class T>
 TArray<T>& TArray<T>::DivElt(const TArray<T>& a)
 {
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::DivElt(const TArray<T>& )  - Not Allocated Array ! ");
   if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::DivElt(const TArray<T>&)")) ;
 
@@ -471,4 +567,91 @@
 }
 
+template <class T>
+TArray<T>& TArray<T>::CopyElt(const TArray<T>& a)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::CopyElt(const TArray<T>& )  - Not Allocated Array ! ");
+  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::MultElt(const TArray<T>&)")) ;
+
+  T * pe;
+  const T * pea;
+  uint_8 j,k,ka;
+  if ((AvgStep() > 0) && (a.AvgStep() > 0) )   {  // regularly spaced elements
+    uint_8 step = AvgStep();
+    uint_8 stepa = a.AvgStep();
+    uint_8 maxx = totsize_*step;
+    pe = Data();
+    pea = a.Data();
+    for(k=0, ka=0;  k<maxx;  k+=step, ka+=stepa )  pe[k] = pea[ka] ;
+  }
+  else {    // Non regular data spacing ...
+    uint_4 ax = MaxSizeKA();
+    uint_8 step = Step(ax);
+    uint_8 stepa = a.Step(ax);
+    uint_8 gpas = Size(ax)*step;
+    for(j=0; j<totsize_; j += Size(ax))  {
+      pe = mNDBlock.Begin()+Offset(j);
+      pea = a.DataBlock().Begin()+a.Offset(j);
+      for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] = pea[ka]; 
+    }
+  }
+  return(*this);
+}
+
+
+// Somme et produit des elements
+template <class T>
+T TArray<T>::Sum() const
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::Sum()  - Not Allocated Array ! ");
+  T ret=0;
+  const T * pe;
+  uint_8 j,k;
+  if (AvgStep() > 0)   {  // regularly spaced elements
+    uint_8 step = AvgStep();
+    uint_8 maxx = totsize_*step;
+    pe = Data();
+    for(k=0; k<maxx; k+=step )  ret += pe[k];
+  }
+  else {    // Non regular data spacing ...
+    uint_4 ka = MaxSizeKA();
+    uint_8 step = Step(ka);
+    uint_8 gpas = Size(ka)*step;
+    for(j=0; j<totsize_; j += Size(ka))  {
+      pe = mNDBlock.Begin()+Offset(j);
+      for(k=0; k<gpas; k+=step)  ret += pe[k] ;
+    }
+  }
+  return ret;
+}
+
+template <class T>
+T TArray<T>::Product() const
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::Product()  - Not Allocated Array ! ");
+  T ret=0;
+  const T * pe;
+  uint_8 j,k;
+  if (AvgStep() > 0)   {  // regularly spaced elements
+    uint_8 step = AvgStep();
+    uint_8 maxx = totsize_*step;
+    pe = Data();
+    for(k=0; k<maxx; k+=step )  ret *= pe[k];
+  }
+  else {    // Non regular data spacing ...
+    uint_4 ka = MaxSizeKA();
+    uint_8 step = Step(ka);
+    uint_8 gpas = Size(ka)*step;
+    for(j=0; j<totsize_; j += Size(ka))  {
+      pe = mNDBlock.Begin()+Offset(j);
+      for(k=0; k<gpas; k+=step)  ret *= pe[k] ;
+    }
+  }
+  return ret;
+}
+
+
 
 // ----------------------------------------------------
@@ -487,7 +670,7 @@
 {
   if (maxprt < 0)  maxprt = max_nprt_;
-  uint_4 npr = 0;
+  int_4 npr = 0;
   Show(os, si);
-  int k0,k1,k2,k3,k4;
+  uint_4 k0,k1,k2,k3,k4;
   for(k4=0; k4<size_[4]; k4++) {
     if (size_[4] > 1) cout << "\n ----- Dimension 5 (U) K4= " << k4 << endl;
@@ -535,11 +718,13 @@
 ///////////////////////////////////////////////////////////////
 #ifdef __CXX_PRAGMA_TEMPLATES__
+/*
 #pragma define_template TArray<uint_1>
+#pragma define_template TArray<int_2>
+#pragma define_template TArray<uint_4>
+#pragma define_template TArray<uint_8>
+*/
 #pragma define_template TArray<uint_2>
-#pragma define_template TArray<int_2>
 #pragma define_template TArray<int_4>
 #pragma define_template TArray<int_8>
-#pragma define_template TArray<uint_4>
-#pragma define_template TArray<uint_8>
 #pragma define_template TArray<r_4>
 #pragma define_template TArray<r_8>
@@ -549,11 +734,13 @@
 
 #if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
-template class TArray<uint_1>;
+/*
+template class TArray<uint_1>; 
+template class TArray<int_2>;  
+template class TArray<uint_4>;
+template class TArray<uint_8>;
+*/
 template class TArray<uint_2>;
-template class TArray<int_2>;
 template class TArray<int_4>;
 template class TArray<int_8>;
-template class TArray<uint_4>;
-template class TArray<uint_8>;
 template class TArray<r_4>;
 template class TArray<r_8>;
Index: trunk/SophyaLib/TArray/tarray.h
===================================================================
--- trunk/SophyaLib/TArray/tarray.h	(revision 803)
+++ trunk/SophyaLib/TArray/tarray.h	(revision 804)
@@ -28,8 +28,8 @@
   // Creation / destruction 
   TArray();
-  TArray(uint_4 ndim, uint_4 * siz, uint_4 step =1);
-  TArray(uint_4 nx, uint_4 ny=0, uint_4 nz=0);
-  TArray(uint_4 ndim, uint_4 * siz, NDataBlock<T> & db, bool share=false, uint_4 step=1, uint_8 offset=0);
-  TArray(uint_4 ndim, uint_4 * siz, T* values, uint_4 step=1, uint_8 offset=0, Bridge* br=NULL);
+  TArray(uint_4 ndim, const uint_4 * siz, uint_4 step =1);
+  TArray(uint_4 nx, uint_4 ny=0, uint_4 nz=0, uint_4 nt=0, uint_4 nu=0);
+  TArray(uint_4 ndim, const uint_4 * siz, NDataBlock<T> & db, bool share=false, uint_4 step=1, uint_8 offset=0);
+  TArray(uint_4 ndim, const uint_4 * siz, T* values, uint_4 step=1, uint_8 offset=0, Bridge* br=NULL);
   TArray(const TArray<T>& a);
   TArray(const TArray<T>& a, bool share);
@@ -38,10 +38,11 @@
 
   // A = B : partage les donnees si "a" est temporaire, clone sinon.
-  virtual TArray<T>& operator = (const TArray<T>& a);
+  inline  TArray<T>& operator = (const TArray<T>& a) { return Set(a); }
+  virtual TArray<T>& Set(const TArray<T>& a);
 
   // Gestion taille/Remplissage
   virtual void Clone(const TArray<T>& a); 
-  virtual void ReSize(uint_4 ndim, uint_4 * siz, uint_4 step=1);
-  virtual void Realloc(uint_4 ndim, uint_4 * siz, uint_4 step=1, bool force=false);
+  void ReSize(uint_4 ndim, uint_4 * siz, uint_4 step=1);
+  void Realloc(uint_4 ndim, uint_4 * siz, uint_4 step=1, bool force=false);
 
   // Compacts size=1 array dimensions
@@ -49,7 +50,11 @@
   virtual TArray<T>& CompactTrailingDimensions();
 
-  // SubArrays
-  virtual TArray<T> operator () (Range rx, Range ry, Range rz=0, Range rt=0, Range ru=0); 
-
+  // Packing array elements in memory
+  virtual TArray<T> PackElements(bool force=false) const ;
+
+  // SubArrays - $CHECK$ Reza 03/2000 je ne sais pas s'il faut declarer ca const ??
+  TArray<T> SubArray(Range rx, Range ry, Range rz, Range rt, Range ru) const ; 
+  inline TArray<T> operator () (Range rx, Range ry, Range rz, Range rt=0, Range ru=0) const
+                   { return  SubArray(rx, ry, rz, rt, ru); }
 
   // ---- Access to data
@@ -83,18 +88,35 @@
   inline T toScalar();
 // Met les elements a une suite de valeurs
-  virtual TArray<T>&  operator = (Sequence seq);
+  virtual TArray<T>&  Set(Sequence seq);
+  inline  TArray<T>&  operator = (Sequence seq)    { return Set(seq); }
 // A = x (tous les elements a x)
-  virtual TArray<T>&  operator = (T x);
+  virtual TArray<T>&  Set(T x);
+  inline  TArray<T>&  operator = (T x)             { return Set(x); }
 // A += -= *= /= x (ajoute, soustrait, ... x a tous les elements)
-  virtual TArray<T>&  operator += (T x);
-  virtual TArray<T>&  operator -= (T x);
-  virtual TArray<T>&  operator *= (T x);
-  virtual TArray<T>&  operator /= (T x);
+  virtual TArray<T>&  Add(T x);
+  inline  TArray<T>&  operator += (T x)            { return Add(x); }
+  virtual TArray<T>&  Sub(T x);
+  inline  TArray<T>&  operator -= (T x)            { return Sub(x); }
+  virtual TArray<T>&  Mul(T x);
+  inline  TArray<T>&  operator *= (T x)            { return Mul(x); }
+  virtual TArray<T>&  Div(T x);
+  inline  TArray<T>&  operator /= (T x)            { return Div(x); }
+  virtual TArray<T>&  SubInv(T x);  //   A ---> x-A
+  virtual TArray<T>&  DivInv(T x);  //   A ---> x/A
+
 // A += -=  (ajoute, soustrait element par element les deux tableaux )
-  virtual TArray<T>&  operator += (const TArray<T>& a);
-  virtual TArray<T>&  operator -= (const TArray<T>& a);
+  virtual TArray<T>&  AddElt(const TArray<T>& a);
+  inline  TArray<T>&  operator += (const TArray<T>& a)  { return AddElt(a); }
+  virtual TArray<T>&  SubElt(const TArray<T>& a);
+  inline  TArray<T>&  operator -= (const TArray<T>& a)  { return SubElt(a); }
 // Multiplication, division element par element les deux tableaux 
-  virtual TArray<T>&  MultElt(const TArray<T>& a);
+  virtual TArray<T>&  MulElt(const TArray<T>& a);
   virtual TArray<T>&  DivElt(const TArray<T>& a);
+// Recopie des valeurs, element par element
+  virtual TArray<T>&  CopyElt(const TArray<T>& a);
+
+// Somme et produit des elements
+  virtual T Sum() const ;
+  virtual T Product() const ;
 
 // Impression, I/O, ...
@@ -124,24 +146,23 @@
 
 template <class T> inline TArray<T> operator + (const TArray<T>& a, T b)
-    {TArray<T> result(a); result.SetTemp(true); result += b; return result;}
+    {TArray<T> result(a); result.SetTemp(true); result.Add(b); return result;}
 
 template <class T> inline TArray<T> operator + (T b,const TArray<T>& a)
-    {TArray<T> result(a); result.SetTemp(true); result += b; return result;}
+    {TArray<T> result(a); result.SetTemp(true); result.Add(b); return result;}
 
 template <class T> inline TArray<T> operator - (const TArray<T>& a, T b)
-    {TArray<T> result(a); result.SetTemp(true); result -= b; return result;}
+    {TArray<T> result(a); result.SetTemp(true); result.Sub(b); return result;}
 
 template <class T> inline TArray<T> operator - (T b,const TArray<T>& a)
-    {TArray<T> result(a); result.SetTemp(true);
-     result.DataBlock() = b-result.DataBlock(); return result;}
+    {TArray<T> result(a); result.SetTemp(true); result.SubInv(b); return result;}
 
 template <class T> inline TArray<T> operator * (const TArray<T>& a, T b)
-    {TArray<T> result(a); result.SetTemp(true); result *= b; return result;}
+    {TArray<T> result(a); result.SetTemp(true); result.Mul(b); return result;}
 
 template <class T> inline TArray<T> operator * (T b,const TArray<T>& a)
-    {TArray<T> result(a); result.SetTemp(true); result *= b; return result;}
+    {TArray<T> result(a); result.SetTemp(true); result.Mul(b); return result;}
 
 template <class T> inline TArray<T> operator / (const TArray<T>& a, T b)
-    {TArray<T> result(a); result.SetTemp(true); result /= b; return result;}
+    {TArray<T> result(a); result.SetTemp(true); result.DivInv(b); return result;}
 
 ////////////////////////////////////////////////////////////////
@@ -150,9 +171,9 @@
 template <class T>
 inline TArray<T> operator + (const TArray<T>& a,const TArray<T>& b)
-    {TArray<T> result(a); result.SetTemp(true); result += b; return result;}
+    {TArray<T> result(a); result.SetTemp(true); result.AddElt(b); return result;}
 
 template <class T>
 inline TArray<T> operator - (const TArray<T>& a,const TArray<T>& b)
-    {TArray<T> result(a); result.SetTemp(true); result += b; return result;}
+    {TArray<T> result(a); result.SetTemp(true); result.SubElt(b); return result;}
 
 
@@ -180,5 +201,5 @@
 {
   CheckBound(ix, iy, iz, it, iu, 4);
-  Elem(ix, iy, iz, it, iu);
+  return(Elem(ix, iy, iz, it, iu));
 }
 
@@ -187,5 +208,5 @@
 {
   CheckBound(ix, iy, iz, it, iu, 4);
-  Elem(ix, iy, iz, it, iu);
+  return(Elem(ix, iy, iz, it, iu));
 }
 
@@ -259,4 +280,7 @@
 }
 
+// Typedef pour simplifier
+typedef TArray<r_8> Array;
+
 } // Fin du namespace
 
Index: trunk/SophyaLib/TArray/tarrinit.cc
===================================================================
--- trunk/SophyaLib/TArray/tarrinit.cc	(revision 803)
+++ trunk/SophyaLib/TArray/tarrinit.cc	(revision 804)
@@ -3,6 +3,4 @@
 #include "tarrinit.h"
 
-#include "tmatrix.h"
-#include "tvector.h"
 #include "fioarr.h"
 
@@ -19,38 +17,26 @@
   
 
-  PPRegister(FIO_TArray<uint_1>);
+  //  PPRegister(FIO_TArray<uint_1>);
+  //  DObjRegister(FIO_TArray<uint_1>, TArray<uint_1>);
   PPRegister(FIO_TArray<uint_2>);
-  PPRegister(FIO_TArray<int_2>);
+  DObjRegister(FIO_TArray<uint_2>, TArray<uint_2>);
+  //  PPRegister(FIO_TArray<int_2>);
+  //  DObjRegister(FIO_TArray<int_2>, TArray<int_2>);
   PPRegister(FIO_TArray<int_4>);
+  DObjRegister(FIO_TArray<int_4>, TArray<int_4>);
   PPRegister(FIO_TArray<int_8>);
-  PPRegister(FIO_TArray<uint_8>);
+  DObjRegister(FIO_TArray<int_8>, TArray<int_8>);
+  //  PPRegister(FIO_TArray<uint_4>);
+  // DObjRegister(FIO_TArray<uint_4>, TArray<uint_4>);
+  //  PPRegister(FIO_TArray<uint_8>);
+  //  DObjRegister(FIO_TArray<uint_8>, TArray<uint_8>);
   PPRegister(FIO_TArray<r_4>);
+  DObjRegister(FIO_TArray<r_4>, TArray<r_4>);
   PPRegister(FIO_TArray<r_8>);
+  DObjRegister(FIO_TArray<r_8>, TArray<r_8>);
   PPRegister(FIO_TArray< complex<r_4> >);
+  DObjRegister(FIO_TArray< complex<r_4> >, TArray< complex<r_4> >);
   PPRegister(FIO_TArray< complex<r_8> >);
-
-  PPRegister(FIO_TMatrix<uint_1>);
-  PPRegister(FIO_TMatrix<uint_2>);
-  PPRegister(FIO_TMatrix<int_2>);
-  PPRegister(FIO_TMatrix<int_4>);
-  PPRegister(FIO_TMatrix<int_8>);
-  PPRegister(FIO_TMatrix<uint_8>);
-  PPRegister(FIO_TMatrix<r_4>);
-  PPRegister(FIO_TMatrix<r_8>);
-  PPRegister(FIO_TMatrix< complex<r_4> >);
-  PPRegister(FIO_TMatrix< complex<r_8> >);
-
-  PPRegister(FIO_TVector<uint_1>);
-  PPRegister(FIO_TVector<uint_2>);
-  PPRegister(FIO_TVector<int_2>);
-  PPRegister(FIO_TVector<int_4>);
-  PPRegister(FIO_TVector<int_8>);
-  PPRegister(FIO_TVector<uint_8>);
-  PPRegister(FIO_TVector<r_4>);
-  PPRegister(FIO_TVector<r_8>);
-  PPRegister(FIO_TVector< complex<r_4> >);
-  PPRegister(FIO_TVector< complex<r_8> >);
-
-
+  DObjRegister(FIO_TArray< complex<r_8> >, TArray< complex<r_8> >);
 
 }
Index: trunk/SophyaLib/TArray/tmatrix.cc
===================================================================
--- trunk/SophyaLib/TArray/tmatrix.cc	(revision 803)
+++ trunk/SophyaLib/TArray/tmatrix.cc	(revision 804)
@@ -1,37 +1,30 @@
-// $Id: tmatrix.cc,v 1.2 2000-03-10 15:13:22 ansari Exp $
+// $Id: tmatrix.cc,v 1.3 2000-04-03 17:35:58 ansari Exp $
 //                         C.Magneville          04/99
 #include "machdefs.h"
 #include <stdio.h>
 #include <stdlib.h>
-#include <iostream.h>
-#include <values.h>
-#include <complex>
 #include "pexceptions.h"
 #include "tmatrix.h"
-#include "objfio.h"
+
+
+
 
 ////////////////////////////////////////////////////////////////
 //**** 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)
+  : TArray<T>()
+{
+}
+
+template <class T>
+TMatrix<T>::TMatrix(uint_4 r,uint_4 c, short mm)
 // 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)
-{
+  :  TArray<T>() 
+{
+  if ( (r == 0) || (c == 0) )
+    throw ParmError("TMatrix<T>::TMatrix(uint_4 r,uint_4 c) NRows or NCols = 0");
+  ReSize(r, c, mm);
 }
 
@@ -39,13 +32,31 @@
 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)
+  : TArray<T>(a)
+{
+}
+
+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)
-{
+: TArray<T>(a, share)
+{
+}
+
+
+template <class T>
+TMatrix<T>::TMatrix(const TArray<T>& a)
+: TArray<T>(a)
+{
+  if (a.NbDimensions() != 2) 
+    throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a) a.NbDimensions() != 2 ");
+}
+
+template <class T>
+TMatrix<T>::TMatrix(const TArray<T>& a, bool share, short mm )
+: TArray<T>(a, share)
+{
+  if (a.NbDimensions() != 2) 
+    throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a, ...) a.NbDimensions() != 2 ");
+  UpdateMemoryMapping(a, mm);
 }
 
@@ -54,256 +65,193 @@
 // Destructeur
 {
-}
-
-////////////////////////////////////////////////////////////////
-// Operations matricielles
-template <class T>
-TMatrix<T> TMatrix<T>::Transpose(void) const
+
+}
+
+template <class T>
+TArray<T>& TMatrix<T>::Set(const TArray<T>& a)
+{
+  if (a.NbDimensions() != 2) 
+    throw SzMismatchError("TMatrix<T>::Set(const TArray<T>& a) a.NbDimensions() != 2 ");
+  return(TArray<T>::Set(a));
+}
+
+template <class T>
+void TMatrix<T>::ReSize(uint_4 r, uint_4 c, short mm)
+{
+  if(r==0||c==0) 
+    throw(SzMismatchError("TMatrix::ReSize r or c==0 "));
+  uint_4 size[BASEARRAY_MAXNDIMS];
+  for(int kk=0; kk<BASEARRAY_MAXNDIMS; kk++)  size[kk] = 0;
+  size[0] = r;  size[1] = c;
+  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
+  UpdateMemoryMapping(mm, size[0], size[1]);
+  TArray<T>::ReSize(2, size, 1);
+  UpdateMemoryMapping(mm, size[0], size[1]);
+}
+
+template <class T>
+void TMatrix<T>::Realloc(uint_4 r,uint_4 c, short mm, bool force)
+{
+  if(r==0||c==0) 
+    throw(SzMismatchError("TMatrix::Realloc r or c==0 "));
+  uint_4 size[BASEARRAY_MAXNDIMS];
+  for(int kk=0; kk<BASEARRAY_MAXNDIMS; kk++)  size[kk] = 0;
+  UpdateMemoryMapping(mm, size[0], size[1]);
+  TArray<T>::Realloc(2, size, 1, force);
+  UpdateMemoryMapping(mm, size[0], size[1]);
+}
+
+// $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
+template <class T>
+TMatrix<T> TMatrix<T>::operator () (Range rline, Range rcol) const
+{
+  uint_4 nx, ny;
+  nx = ny = 0;
+  if (GetMemoryMapping() == CMemoryMapping ) {
+    TMatrix sm(SubArray(rcol, rline, 0, 0, 0),true,CMemoryMapping);
+    sm.UpdateMemoryMapping(CMemoryMapping, nx, ny);
+    sm.SetTemp(true);
+    return(sm);
+  }
+  else {
+    TMatrix sm(SubArray(rline, rcol, 0, 0, 0),true,CMemoryMapping);
+    sm.UpdateMemoryMapping(FortranMemoryMapping, nx, ny);
+    sm.SetTemp(true);
+    return(sm);
+  }
+}
+
+////////////////////////////////////////////////////////////////
 // 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;
-}
+template <class T>
+TMatrix<T>& TMatrix<T>::Transpose() 
+{
+  uint_4 rci = macoli_;
+  macoli_ = marowi_;
+  marowi_ = rci;
+  return(*this);
+}
+
+
+template <class T>
+TMatrix<T> TMatrix<T>::Transpose(short mm)
+{
+  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
+  TMatrix<T> tm(NCols(), NRows(), mm);
+  for(uint_4 i=0; i<NRows(); i++)
+    for(uint_4 j=0; j<NCols(); j++) 
+      tm(j,i) = (*this)(i,j);
+  tm.SetTemp(true);
+  return tm;
+}
+
+// Rearrangement memoire
+template <class T>
+TMatrix<T> TMatrix<T>::Rearrange(short mm)
+{
+  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
+  TMatrix<T> tm(NRows(), NCols(), mm);
+  for(uint_4 i=0; i<NRows(); i++)
+    for(uint_4 j=0; j<NCols(); j++) 
+      tm(i,j) = (*this)(i,j);
+  tm.SetTemp(true);
+  return tm;
+}
+
+template <class T>
+TMatrix<T>& TMatrix<T>::SetIdentity(IdentityMatrix imx)
+{
+  if (ndim_ == 0) {
+    uint_4 sz = imx.Size();
+    if (sz < 1) sz = 1;
+    ReSize(sz, sz);
+  }
+  T diag = (T)imx.Diag();
+  if (NRows() != NCols()) 
+    throw SzMismatchError("TMatrix::operator= (IdentityMatrix) NRows() != NCols()") ;
+  for(uint_4 i=0; i<NRows(); i++) (*this)(i,i) = diag;
+
+  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 *= 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;
-}
-
-
-///////////////////////////////////////////////////////////
-// --------------------------------------------------------
-//   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>::SetDataObj(AnyDataObj & o)
-{
-TMatrix<T> * po = dynamic_cast< TMatrix<T> * >(&o);
-if (po == NULL) return;
-if (ownobj && dobj) delete dobj;
-dobj = po; ownobj = false;
-} 
-
-template <class T>
-void FIO_TMatrix<T>::ReadSelf(PInPersist& is)
-{
-// On lit les 3 premiers uint_4
-//  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);
-}
-
-////////////////////////////////////////////////////////////////
-// -------------------------------------------------------------
-//   La classe de gestion des lignes et colonnes d'une matrice
-// -------------------------------------------------------------
-////////////////////////////////////////////////////////////////
-#include "tvector.h"
+template <class T>
+void TMatrix<T>::Print(ostream& os, int_4 maxprt, bool si) const
+{
+  os << "... TMatrix<T>(NRows=" << NRows() << " x NCols=" << NCols() 
+     << ")::Print() (MemOrg=" << GetMemoryMapping() << ")" << endl;
+  os << " RowsKA= " << RowsKA() << " ColsKA= " << ColsKA() << " VectKA=" << VectKA() << endl;
+  if (maxprt < 0)  maxprt = max_nprt_;
+  int_4 npr = 0;
+  Show(os, si);
+  uint_4 kc,kr;  
+  for(kr=0; kr<size_[marowi_]; kr++) {
+    if ( (size_[marowi_] > 1) && (size_[macoli_] > 10) ) cout << "----- Ligne Line= " << kr << endl;
+    for(kc=0; kc<size_[macoli_]; kc++) {
+      if(kc > 0) os << ", ";  
+      os << (*this)(kr, kc);   npr++; 
+      if (npr >= maxprt) {
+	if (npr < totsize_)  os << "\n     .... " << endl; return;
+      }
+    }
+    os << endl;
+  }
+  os << "\n" << endl;
+}
+
+////////////////////////////////////////////////////////////////
+//****  Multiplication matricielle  *****
+////////////////////////////////////////////////////////////////
+
+template <class T>
+TMatrix<T> TMatrix<T>::Multiply(const TMatrix<T>& b, short mm) const
+{
+  if (NCols() != b.NRows()) 
+    throw(SzMismatchError("TMatrix<T>::Multiply(b) NCols() != b.NRows() ") );
+  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
+  TMatrix<T> rm(NRows(), b.NCols(), mm);
+
+  const T * pea;
+  const T * peb;
+  T sum;
+  uint_4 r,c,k;
+  uint_4 stepa = Step(ColsKA());
+  uint_4 stepb = b.Step(RowsKA());
+  // Calcul de C=rm = A*B   (A=*this)
+  for(r=0; r<rm.NRows(); r++)      // Boucle sur les lignes de A
+    for(c=0; c<rm.NCols(); c++) {     // Boucle sur les colonnes de B
+      sum = 0;
+      pea = &((*this)(r,0));       // 1er element de la ligne r de A
+      peb = &(b(0,c));                // 1er element de la colonne c de B
+      for(k=0; k<NCols(); k++)  sum += pea[k*stepa]*peb[k*stepb];
+      rm(r,c) = sum;
+    }
+
+  rm.SetTemp(true);
+  return rm;
+}
 
 
 ///////////////////////////////////////////////////////////////
 #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<r_8> 
 #pragma define_template TMatrix< complex<r_4> > 
 #pragma define_template TMatrix< complex<r_8> > 
-// 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<r_4> >
-#pragma define_template FIO_TMatrix< complex<r_8> >
 #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<r_4> >;
 template class TMatrix< complex<r_8> >;
-// 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<r_4> >;
-template class FIO_TMatrix< complex<r_8> >;
 #endif
Index: trunk/SophyaLib/TArray/tmatrix.h
===================================================================
--- trunk/SophyaLib/TArray/tmatrix.h	(revision 803)
+++ trunk/SophyaLib/TArray/tmatrix.h	(revision 804)
@@ -5,188 +5,103 @@
 
 #include "machdefs.h"
-#include <stdio.h>
-#include <iostream.h>
-#include <complex>
-#include "ppersist.h"
-#include "anydataobj.h"
-#include "ndatablock.h"
+#include "tarray.h"
 
 namespace SOPHYA {
 
 template <class T>
-class TMatrix : public AnyDataObj {
+class TMatrix : public TArray<T> {
 public:
-
   // Creation / destruction 
   TMatrix();
-  TMatrix(uint_4 r,uint_4 c);
-  TMatrix(uint_4 r,uint_4 c,T* values,Bridge* br=NULL);
+  TMatrix(uint_4 r,uint_4 c, short mm=AutoMemoryMapping);
   TMatrix(const TMatrix<T>& a);
-  TMatrix(const TMatrix<T>& a,bool share);
+  TMatrix(const TMatrix<T>& a, bool share);
+  TMatrix(const TArray<T>& a); 
+  TMatrix(const TArray<T>& a,  bool share, short mm=CMemoryMapping);
   virtual ~TMatrix();
 
-  // Temporaire?
-  inline bool IsTemp(void) const {return mNDBlock.IsTemp();}
-  inline void SetTemp(bool temp=false) const {mNDBlock.SetTemp(temp);}
+  //  Pour verifiez la compatibilite de dimensions lors de l'affectation
+  virtual TArray<T>& Set(const TArray<T>& a);
+  inline  TMatrix<T>& operator = (const TMatrix<T>& a)
+                       { Set(a);  return(*this); }
 
-  // Gestion taille/Remplissage
-  inline void Clone(const TMatrix<T>& a)  // Clone: copie des donnees de "a"
-              {mNDBlock.Clone(a.mNDBlock); mNr = a.mNr; mNc = a.mNc;}
-  inline void Reset(T v=0) {mNDBlock.Reset(v);}
-  inline void ReSize(uint_4 r,uint_4 c)  // Reallocation de place
-       {if(r==0||c==0) throw(SzMismatchError("TMatrix::ReSize r ou c==0\n"));
-        mNDBlock.ReSize(r*c); mNr = r; mNc = c;}
-  inline void Realloc(uint_4 r,uint_4 c,bool force=false)
-       {if(r==0||c==0) throw(SzMismatchError("TMatrix::Realloc r ou c==0\n"));
-        mNDBlock.Realloc(r*c,force); mNr = r; mNc = c;}
+  // Size - Changing the Size 
+  inline uint_4 NRows() const {return Size(marowi_); }
+  inline uint_4 NCols() const {return Size(macoli_); }
+  inline uint_4 NCol()  const {return Size(macoli_); } // back-compat Peida
 
-  // Informations pointeur/data
-  inline uint_4 NRows() const {return mNr;}
-  inline uint_4 NCols() const {return mNc;}
-  inline uint_4 NCol() const {return mNc;} // back-compat Peida
-  inline T const& operator()(uint_4 r,uint_4 c) const
-                            {return *(mNDBlock.Begin()+r*mNc+c);}
-  inline T&       operator()(uint_4 r,uint_4 c)
-                            {return *(mNDBlock.Begin()+r*mNc+c);}
-  inline T const& operator[](uint_4 ip) const
-                            {return *(mNDBlock.Begin()+ip);}
-  inline T&       operator[](uint_4 ip)
-                            {return *(mNDBlock.Begin()+ip);}
-  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;}
+  void ReSize(uint_4 r,uint_4 c, short mm=SameMemoryMapping);  // Reallocation de place
+  void Realloc(uint_4 r,uint_4 c, short mm=SameMemoryMapping, bool force=false);
+
+  // Sub-matrix extraction $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
+  TMatrix<T> operator () (Range rline, Range rcol) const ;
+
+  // Inline element acces methods 
+  inline T const& operator()(uint_4 r,uint_4 c) const;
+  inline T&       operator()(uint_4 r,uint_4 c);
 
   // Operations matricielles
-  TMatrix<T> Transpose(void) const;
+  TMatrix<T>& Transpose();
+  //mm = SameMemoryMapping or CMemoryMapping or FortranMemoryMapping 
+  TMatrix<T>  Transpose(short mm); 
+  // Rearranging Matrix Elements
+  TMatrix<T>  Rearrange(short mm);
 
   // Operateur d'affectation
-  // A = x (matrice diagonale x*Identite)
-  inline TMatrix<T>& operator = (T x)
-    {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:(T)0;
-     return *this;}
-  // A = B : partage les donnees si "a" est temporaire, clone sinon.
-  inline TMatrix<T>& operator = (const TMatrix<T>& a)
-                {if(this == &a) return *this; CloneOrShare(a); return *this;}
+  // A = x (matrice diagonale Identite)
+  virtual TMatrix<T>& SetIdentity(IdentityMatrix imx);
+  inline  TMatrix<T>& operator = (IdentityMatrix imx) { return SetIdentity(imx); }
 
-  // Impression
-  void Print(ostream& os,int lp=0,uint_4 i0=0,uint_4 ni=10,uint_4 j0=0,uint_4 nj=10) const;
-  inline void Print(int lp=0,uint_4 i0=0,uint_4 ni=10,uint_4 j0=0,uint_4 nj=10) const
-              {Print(cout,lp,i0,ni,j0,nj);}
+  // Operations diverses  avec une constante
+  inline  TMatrix<T>&  operator = (T x)             { Set(x); return(*this); }
+  inline  TMatrix<T>&  operator += (T x)            { Add(x); return(*this); }
+  inline  TMatrix<T>&  operator -= (T x)            { Sub(x); return(*this); }
+  inline  TMatrix<T>&  operator *= (T x)            { Mul(x); return(*this); }
+  inline  TMatrix<T>&  operator /= (T x)            { Div(x); return(*this); }
 
-  // Surcharge d'operateurs INPLACE: A (+=,-=,*=,/=) (T) x
-  inline TMatrix<T>& operator += (T b) {mNDBlock += b; return *this;}
-  inline TMatrix<T>& operator -= (T b) {mNDBlock -= b; return *this;}
-  inline TMatrix<T>& operator *= (T b) {mNDBlock *= b; return *this;}
-  inline TMatrix<T>& operator /= (T b) {mNDBlock /= b; return *this;}
+  //  operations avec matrices 
+  inline  TMatrix<T>&  operator += (const TMatrix<T>& a)  { return AddElt(a); }
+  inline  TMatrix<T>&  operator -= (const TMatrix<T>& a)  { return SubElt(a); }
 
-  // Surcharge d'operateurs INPLACE: A (+=,-=,*=,/=) B
-  inline 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;}
-  inline 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;}
-  TMatrix<T>& operator *= (const TMatrix<T>& a);
+  // Input-Output
+  virtual void  Print(ostream& os, int_4 maxprt=-1, bool si=false) const ;
 
-  // Pour surcharge d'operateurs C = A (+,-,*) B
-  TMatrix<T> Add(const TMatrix<T>& b) const;
-  TMatrix<T> Sub(const TMatrix<T>& b) const;
-  TMatrix<T> Mul(const TMatrix<T>& b) const;
+  // Produit matriciel 
+  TMatrix<T>  Multiply(const TMatrix<T>& b, short mm=SameMemoryMapping) const;
+  //  inline TMatrix<T>& operator *= (const TMatrix<T>& b) 
+
+protected:
+};
+
+//  ---- inline acces methods ------
+template <class T>
+inline T const& TMatrix<T>::operator()(uint_4 r, uint_4 c) const
+{
+#ifdef SO_BOUNDCHECKING
+  if (marowi_ == 0)   CheckBound(r, c, 0, 0, 0, 4);
+  else   CheckBound(c, r, 0, 0, 0, 4);
+#endif
+  return ( *( mNDBlock.Begin()+ offset_+ 
+	      r*step_[marowi_] + c*step_[macoli_] ) );
+}
+
+template <class T>
+inline T & TMatrix<T>::operator()(uint_4 r, uint_4 c) 
+{
+#ifdef SO_BOUNDCHECKING
+  if (marowi_ == 0)   CheckBound(r, c, 0, 0, 0, 4);
+  else   CheckBound(c, r, 0, 0, 0, 4);
+#endif
+  return ( *( mNDBlock.Begin()+ offset_+ 
+	      r*step_[marowi_] + c*step_[macoli_] ) );
+}
 
 
-protected:
-  // partage les donnees si "a" temporaire, clone sinon.
-  inline void CloneOrShare(const TMatrix<T>& a)
-              {mNDBlock.CloneOrShare(a.mNDBlock); mNr=a.mNr; mNc=a.mNc;}
-  // Share: partage les donnees de "a"
-  inline void Share(const TMatrix<T>& a)
-              {mNDBlock.Share(a.mNDBlock); mNr=a.mNr; mNc=a.mNc;}
+// Surcharge d'operateurs C = A * B
 
-  uint_4 mNr,mNc;
-  NDataBlock<T> mNDBlock;
-};
+template <class T> inline TMatrix<T> operator * (const TMatrix<T>& a, const TMatrix<T>& b)
+{ TMatrix<T> result(a); result.SetTemp(true); result.Multiply(b);  return result;}
 
-////////////////////////////////////////////////////////////////
-// Impression
-
-template <class T>
-inline ostream& operator << (ostream& os, const TMatrix<T>& a)
-                            {a.Print(os); return(os);}
-
-////////////////////////////////////////////////////////////////
-// Surcharge d'operateurs A (+,-,*,/) (T) x
-
-template <class T> inline TMatrix<T> operator + (const TMatrix<T>& a, T b)
-    {TMatrix<T> result(a); result.SetTemp(true); result += b; return result;}
-
-template <class T> inline TMatrix<T> operator + (T b,const TMatrix<T>& a)
-    {TMatrix<T> result(a); result.SetTemp(true); result += b; return result;}
-
-template <class T> inline TMatrix<T> operator - (const TMatrix<T>& a, T b)
-    {TMatrix<T> result(a); result.SetTemp(true); result -= b; return result;}
-
-template <class T> inline TMatrix<T> operator - (T b,const TMatrix<T>& a)
-    {TMatrix<T> result(a); result.SetTemp(true);
-     result.DataBlock() = b-result.DataBlock(); return result;}
-
-template <class T> inline TMatrix<T> operator * (const TMatrix<T>& a, T b)
-    {TMatrix<T> result(a); result.SetTemp(true); result *= b; return result;}
-
-template <class T> inline TMatrix<T> operator * (T b,const TMatrix<T>& a)
-    {TMatrix<T> result(a); result.SetTemp(true); result *= b; return result;}
-
-template <class T> inline TMatrix<T> operator / (const TMatrix<T>& a, T b)
-    {TMatrix<T> result(a); result.SetTemp(true); result /= b; return result;}
-
-////////////////////////////////////////////////////////////////
-// Surcharge d'operateurs C = A (+,-,*,/) B
-
-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);}
-
-////////////////////////////////////////////////////////////////
-// Typedef pour simplifier et compatibilite Peida
 typedef TMatrix<r_8> Matrix;
-
-/////////////////////////////////////////////////////////////////////////
-// 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();
-  virtual void        SetDataObj(AnyDataObj & o);
-  inline operator TMatrix<T>() { return(*dobj); }
-protected :
-  virtual void ReadSelf(PInPersist&);           
-  virtual void WriteSelf(POutPersist&) const;  
-  TMatrix<T> * dobj;
-  bool ownobj;
-};
-
-template <class T>
-inline POutPersist& operator << (POutPersist& os, TMatrix<T> & obj)
-{ FIO_TMatrix<T> fio(&obj);  fio.Write(os);  return(os); }
-template <class T>
-inline PInPersist& operator >> (PInPersist& is, TMatrix<T> & obj)
-{ FIO_TMatrix<T> fio(&obj);  fio.Read(is);  return(is); }
-
 
 } // Fin du namespace
Index: trunk/SophyaLib/TArray/tvector.cc
===================================================================
--- trunk/SophyaLib/TArray/tvector.cc	(revision 803)
+++ trunk/SophyaLib/TArray/tvector.cc	(revision 804)
@@ -1,188 +1,132 @@
-// $Id: tvector.cc,v 1.1.1.1 2000-03-02 16:48:24 ansari Exp $
+// $Id: tvector.cc,v 1.2 2000-04-03 17:35:59 ansari Exp $
 //                         C.Magneville          04/99
 #include "machdefs.h"
-#include <stdio.h>
 #include <stdlib.h>
-#include <iostream.h>
-#include <complex>
 #include "pexceptions.h"
 #include "tvector.h"
-#include "objfio.h"
+
+////////////////////////////////////////////////////////////////
+//**** Createur, Destructeur
 
 template <class T>
-TVector<T>::TVector(uint_4 n)
-// Constructeur d'un vecteur de "n" elements
-: TMatrix<T>(n,1)
+TVector<T>::TVector()
+  : TMatrix<T>()
 {
 }
 
 template <class T>
-TVector<T>::TVector(uint_4 n, T* values,Bridge* br)
-// Construit un vecteur de n elements. On fournit
-// le tableau des valeurs et eventuellement un Bridge.
-: TMatrix<T>(n,1,values,br)
+TVector<T>::TVector(uint_4 n, short lcv, short mm)
+// Constructeur 
+  : TMatrix<T>((lcv == ColumnVector) ? n : 1,(lcv == ColumnVector) ? 1 : 0, mm)
+{
+  if (lcv == ColumnVector) veceli_ = marowi_;
+  else veceli_ = macoli_;
+}
+
+
+template <class T>
+TVector<T>::TVector(const TVector<T>& a)
+// Constructeur par copie (partage si "a" temporaire).
+  : TMatrix<T>(a)
 {
 }
 
 template <class T>
-TVector<T>::TVector(const TVector<T>& v)
-// Constructeur par copie (partage si "v" temporaire).
-: TMatrix<T>(v)
+TVector<T>::TVector(const TVector<T>& a, bool share)
+// Constructeur par copie avec possibilite de forcer le partage ou non.
+: TMatrix<T>(a, share)
 {
 }
 
 template <class T>
-TVector<T>::TVector(const TVector<T>& v,bool share)
-// Constructeur par copie avec possibilite de forcer le partage ou non.
-: TMatrix<T>(v,share)
+TVector<T>::TVector(const TArray<T>& a)
+: TMatrix<T>(a)
 {
 }
 
+
 template <class T>
-TVector<T>::TVector(const TMatrix<T>& a)
-// Constructeur a partir d'une matrice "a" de dimension (n,1)
-: TMatrix<T>(a)
+TVector<T>::TVector(const TArray<T>& a, bool share, short lcv, short mm )
+: TMatrix<T>(a, share, mm)
 {
-if(a.NCols() != 1)
-  throw(SzMismatchError("TVector(const TMatrix<T>& a) size mismatch, ncol!=1"));
+  if (lcv == SameTypeVector) veceli_ = a.VectKA();
+  else {
+    if (lcv == ColumnVector) veceli_ = marowi_;
+    else veceli_ = macoli_;
+  }
+}
+
+template <class T>
+TVector<T>::~TVector()
+// Destructeur
+{
+
+}
+
+template <class T>
+void TVector<T>::ReSize(uint_4 n, short lcv)
+{
+  if( n == 0 ) 
+    throw(SzMismatchError("TVector::ReSize() n = 0 "));
+  uint_4 r,c;
+  if (lcv == SameTypeVector)  lcv = GetVectorType();
+  if (lcv == ColumnVector) { r = n;  c = 1; }
+  else { c = n; r = 1; }
+  TMatrix<T>::ReSize(r,c);
+}
+
+template <class T>
+void TVector<T>::Realloc(uint_4 n, short lcv, bool force)
+{
+  if( n == 0 ) 
+    throw(SzMismatchError("TVector::Realloc() n = 0 "));
+  uint_4 r,c;
+  if (lcv == SameTypeVector)  lcv = GetVectorType();
+  if (lcv == ColumnVector) { r = n;  c = 1; }
+  else { c = n; r = 1; }
+  TMatrix<T>::Realloc(r,c,SameMemoryMapping,force);
+}
+
+// $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
+template <class T>
+TVector<T> TVector<T>::operator () (Range relt) const
+{
+  Range rr, cr;
+  if (GetVectorType() == ColumnVector )  rr = relt; 
+  else cr = relt;
+  TMatrix<T> const & mtx = (*this);
+  TVector sv( mtx(rr, cr) , true, GetVectorType(), GetMemoryMapping() );
+  sv.SetTemp(true);
+  return(sv);  
+}
+
+template <class T>
+T TVector<T>::Norm2() const
+{
+  T ret = 0;
+  for(uint_8 k=0; k<Size(); k++)   ret += (*this)(k)*(*this)(k);
+  return ret;
 }
 
 
-///////////////////////////////////////////////////////////
-// --------------------------------------------------------
-//   Les objets delegues pour la gestion de persistance 
-// --------------------------------------------------------
-///////////////////////////////////////////////////////////
-
-template <class T>
-FIO_TVector<T>::FIO_TVector()
-{
-dobj=new TVector<T>;
-ownobj=true;
-}
-
-template <class T>
-FIO_TVector<T>::FIO_TVector(string const & filename) 
-{
-dobj=new TVector<T>;
-ownobj=true; 
-Read(filename);
-}
-
-template <class T>
-FIO_TVector<T>::FIO_TVector(const TVector<T> & obj) 
-{ 
-dobj = new TVector<T>(obj);
-ownobj=true; 
-}
-
-template <class T>
-FIO_TVector<T>::FIO_TVector(TVector<T> * obj) 
-{ 
-dobj = obj;
-ownobj=false; 
-}
-
-template <class T>
-FIO_TVector<T>::~FIO_TVector()
-{
-if (ownobj && dobj) delete dobj;
-}
-
-template <class T>
-AnyDataObj* FIO_TVector<T>::DataObj()
-{
-return(dobj);
-}
-
-template <class T>
-void FIO_TVector<T>::SetDataObj(AnyDataObj & o)
-{
-TVector<T> * po = dynamic_cast< TVector<T> * >(&o);
-if (po == NULL) return;
-if (ownobj && dobj) delete dobj;
-dobj = po; ownobj = false;
-} 
-
-template <class T>
-void FIO_TVector<T>::ReadSelf(PInPersist& is)
-{
-// On lit les 3 premiers uint_4
-//  0: Numero de version,  1 : NRows=NElts,  2 : NCol=1
-uint_4 itab[3];
-is.Get(itab,3);
-if (dobj == NULL) dobj = new TVector<T>(itab[1]);
-else dobj->ReSize(itab[1]);
-// On lit le NDataBlock
-FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
-fio_nd.Read(is);
-}
-
-template <class T>
-void FIO_TVector<T>::WriteSelf(POutPersist& os) const
-{
-if (dobj == NULL)   return;
-//  On ecrit 3 uint_4 .... 
-//  0: Numero de version,  1 : NRows=NElts,  2 : NCol=1
-uint_4 itab[3];
- itab[0] = 1;  // Numero de version a 1
-itab[1] = dobj->NElts();
-itab[2] = 1;
-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 TVector<uint_1>
 #pragma define_template TVector<uint_2>
-#pragma define_template TVector<int_2>
 #pragma define_template TVector<int_4>
 #pragma define_template TVector<int_8>
-#pragma define_template TVector<uint_4>
-#pragma define_template TVector<uint_8>
 #pragma define_template TVector<r_4>
 #pragma define_template TVector<r_8>
 #pragma define_template TVector< complex<r_4> > 
 #pragma define_template TVector< complex<r_8> > 
-// Instances des delegues FileIO (PPersist)
-#pragma define_template FIO_TVector<uint_1>
-#pragma define_template FIO_TVector<uint_2>
-#pragma define_template FIO_TVector<int_2>
-#pragma define_template FIO_TVector<int_4>
-#pragma define_template FIO_TVector<int_8>
-#pragma define_template FIO_TVector<uint_4>
-#pragma define_template FIO_TVector<uint_8>
-#pragma define_template FIO_TVector<r_8>
-#pragma define_template FIO_TVector<r_4>
-#pragma define_template FIO_TVector< complex<r_4> >
-#pragma define_template FIO_TVector< complex<r_8> >
 #endif
 
 #if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
-template class TVector<uint_1>;
 template class TVector<uint_2>;
-template class TVector<int_2>;
 template class TVector<int_4>;
 template class TVector<int_8>;
-template class TVector<uint_4>;
-template class TVector<uint_8>;
 template class TVector<r_4>;
 template class TVector<r_8>;
 template class TVector< complex<r_4> >;
 template class TVector< complex<r_8> >;
-// Instances des delegues FileIO (PPersist)
-template class FIO_TVector<uint_1>;
-template class FIO_TVector<uint_2>;
-template class FIO_TVector<int_2>;
-template class FIO_TVector<int_4>;
-template class FIO_TVector<int_8>;
-template class FIO_TVector<uint_4>;
-template class FIO_TVector<uint_8>;
-template class FIO_TVector<r_8>;
-template class FIO_TVector<r_4>;
-template class FIO_TVector< complex<r_4> >;
-template class FIO_TVector< complex<r_8> >;
 #endif
+
Index: trunk/SophyaLib/TArray/tvector.h
===================================================================
--- trunk/SophyaLib/TArray/tvector.h	(revision 803)
+++ trunk/SophyaLib/TArray/tvector.h	(revision 804)
@@ -11,79 +11,75 @@
 class TVector : public TMatrix<T> {
 public:
+  enum {SameTypeVector=0, ColumnVector=1, LineVector=2};
+  // Creation / destruction 
+  TVector();
+  TVector(uint_4 n, short lcv=ColumnVector, short mm=AutoMemoryMapping);
+  TVector(const TVector<T>& v);
+  TVector(const TVector<T>& v, bool share);
+  TVector(const TArray<T>& a);
+  TVector(const TArray<T>& a,  bool share, short lcv=ColumnVector, short mm=CMemoryMapping);
 
-  // Creation / destruction 
-  TVector(uint_4 n=1);
-  TVector(uint_4 n, T* values,Bridge* br=NULL);
-  TVector(const TVector<T>& v);
-  TVector(const TVector<T>& v,bool share);
-  TVector(const TMatrix<T>& a);
+  virtual ~TVector();
+
+  inline  TVector<T>& operator = (const TVector<T>& a)
+                       { Set(a);  return(*this); }
+
+  // Vector type   Line or Column vector
+  inline short GetVectorType() const 
+               { return((marowi_ == veceli_) ? ColumnVector : LineVector); }
 
   // Gestion taille/Remplissage
-  inline void ReSize(uint_4 n) {TMatrix<T>::ReSize(n,1);} // Reallocation de place
-  inline void Realloc(uint_4 n,bool force=false) {TMatrix<T>::Realloc(n,1,force);}
+  void ReSize(uint_4 n, short lcv=SameTypeVector ); 
+  void Realloc(uint_4 n, short lcv=SameTypeVector, bool force=false);
+
+  // Sub-Vector extraction $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
+  TVector<T> operator () (Range relt) const ;
 
   // Informations pointeur/data
-  inline uint_4 NElts() const {return NRows();}
+  inline uint_4 NElts() const {return Size(); }
   
-  // Acces aux elements
-  inline T& operator()(uint_4 n) {return (*this)[n];}
-  inline T const& operator()(uint_4 n) const {return (*this)[n];}
-  inline T& Element(uint_4 n) {return (*this)[n];}
-  inline T const& Element(uint_4 n) const {return (*this)[n];}
+  // Inline element acces methods 
+  inline T const& operator()(uint_4 n) const;
+  inline T&       operator()(uint_4 n);
 
-  // Operateur d'affectation
-  inline TVector& operator = (const TVector& v)
-                  {TMatrix<T>::operator =(v); return *this;}
-  inline TVector& operator = (T x)
-                  {for(uint_4 i=0;i<NRows();i++) (*this)(i)=x; return *this;}
+  // Operations diverses  avec une constante
+  inline  TVector<T>&  operator = (T x)             { Set(x); return(*this); }
+  inline  TVector<T>&  operator += (T x)            { Add(x); return(*this); }
+  inline  TVector<T>&  operator -= (T x)            { Sub(x); return(*this); }
+  inline  TVector<T>&  operator *= (T x)            { Mul(x); return(*this); }
+  inline  TVector<T>&  operator /= (T x)            { Div(x); return(*this); }
 
+  //  operations avec matrices 
+  inline  TVector<T>&  operator += (const TVector<T>& a)  { AddElt(a); return(*this); }
+  inline  TVector<T>&  operator -= (const TVector<T>& a)  { SubElt(a); return(*this); }
 
+  // Norme(^2) 
+  T Norm2() const ;
 };
 
-// produit scalaire, matrice*vecteur
+//  ---- inline acces methods ------
 template <class T>
-inline T operator* (const TVector<T>& v1, const TVector<T>& v2)
-  {if(v1.NRows() != v2.NRows())
-     throw(SzMismatchError("TVector::operator*(TVector& v1,TVector v2) size mismatch"));
-   T *p = const_cast<T *>(v1.Data()), *pEnd = p+v1.NElts(),
-     *q = const_cast<T *>(v2.Data()), r = 0;
-   while (p<pEnd) r += *p++ * *q++;
-   return r;}
+inline T const& TVector<T>::operator()(uint_4 n) const
+{
+#ifdef SO_BOUNDCHECKING
+  if (veceli__ == 0)   CheckBound(n, 0, 0, 0, 0, 4);
+  else   CheckBound(0, n, 0, 0, 0, 4);
+#endif
+  return ( *( mNDBlock.Begin()+ offset_ + n*step_[veceli_] ) );
+}
 
 template <class T>
-inline TVector<T> operator* (const TMatrix<T>& a, const TVector<T>& b)
-   {return TVector<T>(a * ((TMatrix<T> const&)(b)));}
-
+inline T & TVector<T>::operator()(uint_4 n) 
+{
+#ifdef SO_BOUNDCHECKING
+  if (veceli__ == 0)   CheckBound(n, 0, 0, 0, 0, 4);
+  else   CheckBound(0, n, 0, 0, 0, 4);
+#endif
+  return ( *( mNDBlock.Begin()+ offset_ + n*step_[veceli_] ) );
+}
 
 // Typedef pour simplifier et compatibilite Peida
 typedef TVector<r_8> Vector;
 
-/////////////////////////////////////////////////////////////////////////
-// Classe pour la gestion de persistance
-template <class T>
-class FIO_TVector : public  PPersist  {
-public:
-  FIO_TVector();
-  FIO_TVector(string const & filename); 
-  FIO_TVector(const TVector<T> & obj);
-  FIO_TVector(TVector<T> * obj);
-  virtual ~FIO_TVector();
-  virtual AnyDataObj* DataObj();
-  virtual void        SetDataObj(AnyDataObj & o);
-  inline operator TVector<T>() { return(*dobj); }
-protected :
-  virtual void ReadSelf(PInPersist&);           
-  virtual void WriteSelf(POutPersist&) const;  
-  TVector<T> * dobj;
-  bool ownobj;
-};
-
-template <class T>
-inline POutPersist& operator << (POutPersist& os, TVector<T> & obj)
-{ FIO_TVector<T> fio(&obj);  fio.Write(os);  return(os); }
-template <class T>
-inline PInPersist& operator >> (PInPersist& is, TVector<T> & obj)
-{ FIO_TVector<T> fio(&obj);  fio.Read(is);  return(is); }
-
 } // Fin du namespace
 
Index: trunk/SophyaLib/TArray/utilarr.cc
===================================================================
--- trunk/SophyaLib/TArray/utilarr.cc	(revision 803)
+++ trunk/SophyaLib/TArray/utilarr.cc	(revision 804)
@@ -27,4 +27,5 @@
 }
 
+/*
 Range & Range::operator = (uint_4 start)
 {
@@ -34,2 +35,10 @@
   return (*this);
 }
+*/
+
+
+IdentityMatrix::IdentityMatrix(double diag, uint_4 n)
+{
+  size_ = n;
+  diag_ = diag;
+}
Index: trunk/SophyaLib/TArray/utilarr.h
===================================================================
--- trunk/SophyaLib/TArray/utilarr.h	(revision 803)
+++ trunk/SophyaLib/TArray/utilarr.h	(revision 804)
@@ -18,5 +18,5 @@
 class Sequence {
 public:
-  Sequence (double start=0., double step=1., Arr_DoubleFunctionOfX f=NULL);
+  explicit Sequence (double start=0., double step=1., Arr_DoubleFunctionOfX f=NULL);
   inline double & Start() { return start_; }
   inline double & Step() { return step_; }
@@ -29,11 +29,20 @@
 class Range {
 public:
-  Range(uint_4 start=0, uint_4 size=1, uint_4 step=1);
+  explicit Range(uint_4 start=0, uint_4 size=1, uint_4 step=1);
   inline uint_4 & Start()  {  return start_; }
   inline uint_4 & Size()  {  return size_; }
   inline uint_4 & Step()  {  return step_; }
-  Range & operator = (uint_4 start);
 protected:
   uint_4 start_, size_, step_;
+};
+
+class IdentityMatrix {
+public:
+  explicit IdentityMatrix(double diag=1., uint_4 n=0);
+  inline uint_4 Size() { return size_; }
+  inline double Diag() { return diag_; }
+protected:
+  uint_4 size_;
+  double diag_;
 };
 
