Index: /trunk/SophyaLib/TArray/basarr.cc
===================================================================
--- /trunk/SophyaLib/TArray/basarr.cc	(revision 812)
+++ /trunk/SophyaLib/TArray/basarr.cc	(revision 813)
@@ -14,10 +14,21 @@
 				 "Stride(int )", "ElemCheckBound()", "operator()" };
 uint_4 BaseArray::max_nprt_ = 50;
+uint_4 BaseArray::prt_lev_ = 0;
 short  BaseArray::default_memory_mapping = CMemoryMapping;
-
-// Methodes statiques globales
-void BaseArray::SetMaxPrint(uint_4 nprt)
+short  BaseArray::default_vector_type = ColumnVector;
+uint_8 BaseArray::openmp_size_threshold = 200000;
+
+// ------ Methodes statiques globales --------
+
+//  1/ Gestion des impressions 
+void BaseArray::SetMaxPrint(uint_4 nprt, uint_4 lev)
 {
   max_nprt_ = nprt;
+  prt_lev_ = (lev < 3) ? lev : 3;
+}
+
+void BaseArray::SetOpenMPSizeThreshold(uint_8 thr)
+{
+  openmp_size_threshold = thr;
 }
 
@@ -32,9 +43,15 @@
 short BaseArray::SetDefaultMemoryMapping(short mm)
 {
-  default_memory_mapping = ( (mm == CMemoryMapping) ? CMemoryMapping : FortranMemoryMapping) ; 
+  default_memory_mapping = (mm != CMemoryMapping) ? FortranMemoryMapping : CMemoryMapping; 
   return default_memory_mapping;
 }
 
-
+short BaseArray::SetDefaultVectorType(short vt)
+{
+  default_vector_type = (vt != ColumnVector) ? RowVector : ColumnVector ; 
+  return default_vector_type;
+}
+
+// Memory organisation 
 short BaseArray::SelectMemoryMapping(short mm)
 // si mm == CMemoryMapping, elements d'une ligne suivant X (consecutif)
@@ -44,31 +61,88 @@
   else return (default_memory_mapping);
 }
-
-void BaseArray::UpdateMemoryMapping(short mm, uint_4 & nx, uint_4 & ny)
-{
+short BaseArray::SelectVectorType(short vt)
+{
+  if ((vt == ColumnVector) || (vt == RowVector))  return(vt);
+  else return(default_vector_type);
+}
+
+void BaseArray::UpdateMemoryMapping(short mm)
+{
+  short vt = default_vector_type;
   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;
-  }  
+    marowi_  = 1;  macoli_ = 0;
+  }
+  else {
+    marowi_ = 0;  macoli_ = 1;
+  }
+
+  if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
+    // Choix automatique Vecteur ligne ou colonne
+    if ( size_[macoli_] == 1)  veceli_ = marowi_;
+    else veceli_ = macoli_;
+  }
+  else veceli_ = (vt ==  ColumnVector ) ?  marowi_ : macoli_;
+  ck_memo_vt_ = true;   // Check MemMapping and VectorType for CompareSize  
 }
 
 void BaseArray::UpdateMemoryMapping(BaseArray const & a, short mm)
 {
-  if (mm == SameMemoryMapping) 
+  short vt = default_vector_type;
+  if (mm == SameMemoryMapping) {
     mm = ((a.marowi_ == 1) ? CMemoryMapping : FortranMemoryMapping);
+    vt = (a.marowi_ == a.veceli_) ? ColumnVector : RowVector;
+  }
+  else if (mm == AutoMemoryMapping)   mm = default_memory_mapping;
+
   if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) mm = default_memory_mapping;
   if (mm == CMemoryMapping) {
-    marowi_ = veceli_ = 1;  macoli_ = 0;
-  }
-  else {
-    marowi_ = veceli_ = 0;  macoli_ = 1;
-  }  
-}
-
+    marowi_  = 1;  macoli_ = 0;
+  }
+  else {
+    marowi_ = 0;  macoli_ = 1;
+  }
+  if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
+    // Choix automatique Vecteur ligne ou colonne
+    if ( size_[macoli_] == 1)  veceli_ = marowi_;
+    else veceli_ = marowi_;
+  }
+  else veceli_ = (vt ==  ColumnVector ) ?  marowi_ : macoli_;
+  ck_memo_vt_ = true;   // Check MemMapping and VectorType for CompareSize  
+}
+
+void BaseArray::SetMemoryMapping(short mm)
+{
+  if (mm == SameMemoryMapping) return;
+  if (mm == AutoMemoryMapping)  mm = default_memory_mapping;
+  if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) return;
+  short vt = (marowi_ == veceli_) ? ColumnVector : RowVector;
+  if (mm == CMemoryMapping) {
+    marowi_ =  1;  macoli_ = 0;
+  }
+  else {
+    marowi_ =  0;  macoli_ = 1;
+  }
+  if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
+    // Choix automatique Vecteur ligne ou colonne
+    if ( size_[macoli_] == 1)  veceli_ = marowi_;
+    else veceli_ = macoli_;
+  }
+  else   veceli_ = (vt ==  ColumnVector ) ?  marowi_ : macoli_;
+  ck_memo_vt_ = true;   // Check MemMapping and VectorType for CompareSize  
+}
+
+void BaseArray::SetVectorType(short vt)
+{
+  if (vt == SameVectorType) return;
+  if (vt == AutoVectorType) vt =  default_vector_type;
+  if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
+    // Choix automatique Vecteur ligne ou colonne
+    if ( size_[macoli_] == 1)  veceli_ = marowi_;
+    else veceli_ = macoli_;
+  }
+  else  veceli_ = (vt ==  ColumnVector ) ?  marowi_ : macoli_;
+  ck_memo_vt_ = true;   // Check MemMapping and VectorType for CompareSize  
+}
 
 // -------------------------------------------------------
@@ -87,7 +161,13 @@
   moystep_ = 0;
   offset_ = 0;
-  // Default for matrices : Lines = Y , Columns = X, Default Vectore= ColumnVector
-  macoli_ = 0;
-  veceli_ = marowi_ = 1;
+  // Default for matrices : Memory organisation and Vector type 
+  if (default_memory_mapping == CMemoryMapping) {
+    marowi_ =  1;  macoli_ = 0;
+  }
+  else {
+    marowi_ =  0;  macoli_ = 1;
+  }
+  veceli_ = (default_vector_type ==  ColumnVector ) ?  marowi_ : macoli_;
+  ck_memo_vt_ = false;   // Default : Don't Check MemMapping and VectorType for CompareSize  
 }
 
@@ -103,6 +183,8 @@
   for(int k=0; k<ndim_; k++) 
     if (size_[k] != a.size_[k])  return(false);
-  if ( (macoli_ != a.macoli_) || (marowi_ != a.marowi_) ||
-       (veceli_ != a.veceli_) )  return(false);
+  //  $CHECK$   Reza doit-on verifier ca 
+  if (ck_memo_vt_ && a.ck_memo_vt_)
+    if ( (macoli_ != a.macoli_) || (marowi_ != a.marowi_) ||
+         (veceli_ != a.veceli_) )  return(false);
   return(true);
 }
@@ -166,15 +248,55 @@
 
 //  Acces lineaire aux elements ....  Calcul d'offset
+// --------------------------------------------------
+// Position de l'element 0 du vecteur i selon l'axe ka
+// --------------------------------------------------
+uint_8 BaseArray::Offset(uint_4 ka, uint_8 i) const
+{
+
+  if ( (ndim_ < 1) || (i == 0) )  return(offset_);
+  //#ifdef SO_BOUNDCHECKING
+  if (ka >= ndim_) 
+    throw RangeCheckError("BaseArray::Offset(uint_4 ka, uint_8 i) Axe KA Error");
+  if ( i*size_[ka] >= totsize_ )  
+    throw RangeCheckError("BaseArray::Offset(uint_4 ka, uint_8 i) Index Error");
+  //#endif
+  uint_4 idx[BASEARRAY_MAXNDIMS];
+  int k;
+  uint_8 rest = i;
+  idx[ka] = 0;
+  for(k=0; k<ndim_; k++) {
+    if (k == ka) continue;
+    idx[k] = rest%size_[k];   rest /= size_[k];
+  }
+  uint_8 off = offset_;
+  for(k=0; k<ndim_; k++)  off += idx[k]*step_[k];
+  return (off);
+}
 
 uint_8 BaseArray::Offset(uint_8 ip) const
 {
-if (ip == 0)  return(offset_);
-uint_4 idx[BASEARRAY_MAXNDIMS];
-int k;
-uint_8 rest = ip;
-for(k=0; k<ndim_; k++)   { idx[k] = rest%size_[k];   rest /= size_[k];  }   
-uint_8 off = offset_;
-for(k=0; k<ndim_; k++)  off += idx[k]*step_[k];
-return (off);
+  if ( (ndim_ < 1) || (ip == 0) )  return(offset_);
+  //#ifdef SO_BOUNDCHECKING
+  if (ip >= totsize_)
+    throw RangeCheckError("BaseArray::Offset(uint_8 ip) Out of range index ip");
+  //#endif
+
+  uint_4 idx[BASEARRAY_MAXNDIMS];
+  int k;
+  uint_8 rest = ip;
+  for(k=0; k<ndim_; k++) {
+    idx[k] = rest%size_[k];   rest /= size_[k];
+  }
+  //#ifdef SO_BOUNDCHECKING
+  if (rest != 0) 
+    throw PError("BaseArray::Offset(uint_8 ip) GUG !!! rest != 0");
+  //#endif
+//   if (rest != 0) cerr << " BUG ---- BaseArray::Offset( " << ip << " )" << rest << endl;
+//   cerr << " DBG-Offset( " << ip << ")" ;
+//   for(k=0; k<ndim_; k++) cerr << idx[k] << "," ;
+//   cerr << " ZZZZ " << endl;
+  uint_8 off = offset_;
+  for(k=0; k<ndim_; k++)  off += idx[k]*step_[k];
+  return (off);
 }
 
@@ -186,19 +308,32 @@
 void BaseArray::Show(ostream& os, bool si) const
 {
-  os << " TArray< " << DataType() << " >  NDim= " << ndim_
-     << "(" << typeid(*this).name() << " )" << endl;
-  os << "TotSize=" << totsize_ << " Size(X*Y*...)= " ;
+  os << "\n--- " << InfoString() ; 
+  os << " ND=" << ndim_ << " SizeX*Y*...= " ;
   for(int k=0; k<ndim_; k++) { 
     os << size_[k];
-    if (k<ndim_-1)  os << " x ";
-  }
-  os << endl;
-  os <<  " Step(X Y ...)= " ;
-  for(int k=0; k<ndim_; k++)     os << step_[k] << "  " ;
-  os << " Offset= " << offset_  << "\n " << endl;
-  if (si && (mInfo != NULL)) os << (*mInfo) << endl;
+    if (k<ndim_-1)  os << "x";
+  }
+  os << " ---" << endl;
+  if (prt_lev_ > 0) {
+    os <<  " TotSize= " << totsize_ << " Step(X Y ...)="  ;
+    for(int k=0; k<ndim_; k++)     os << step_[k] << "  " ;
+    os << " Offset= " << offset_  << endl;
+  }
+  if (prt_lev_ > 1) {
+    os << " MemoryMapping=" << GetMemoryMapping() << " VecType= " << GetVectorType()
+       << " RowsKA= " << RowsKA() << " ColsKA= " << ColsKA() 
+       << " VectKA=" << VectKA() << endl;
+  }
+  if (!si && (prt_lev_ < 2)) return;
+  if (mInfo != NULL) os << (*mInfo) << endl;
  
 }
 
+string BaseArray::InfoString() const
+{
+  string rs = "BaseArray Type= ";
+  rs +=  typeid(*this).name() ;
+  return rs;
+}
 
 DVList& BaseArray::Info()
@@ -238,7 +373,13 @@
   }
   offset_ = offset;
-  // Default for matrices : Lines = Y , Columns = X
-  macoli_ = 0;
-  marowi_ = veceli_ = 1;
+  // Default for matrices : Memory organisation and Vector type 
+  if (default_memory_mapping == CMemoryMapping) {
+    marowi_ =  1;  macoli_ = 0;
+  }
+  else {
+    marowi_ =  0;  macoli_ = 1;
+  }
+  veceli_ = (default_vector_type ==  ColumnVector ) ?  marowi_ : macoli_;
+  ck_memo_vt_ = false;   // Default : Don't Check MemMapping and VectorType for CompareSize  
   // Update OK 
   ndim_ = ndim;
@@ -280,7 +421,13 @@
   minstep_ = minstep;
   offset_ = offset;
-  // Default for matrices : Lines = Y , Columns = X
-  macoli_ = 0;
-  marowi_ = veceli_ = 1;
+  // Default for matrices : Memory organisation and Vector type 
+  if (default_memory_mapping == CMemoryMapping) {
+    marowi_ =  1;  macoli_ = 0;
+  }
+  else {
+    marowi_ =  0;  macoli_ = 1;
+  }
+  veceli_ = (default_vector_type ==  ColumnVector ) ?  marowi_ : macoli_;
+  ck_memo_vt_ = false;   // Default : Don't Check MemMapping and VectorType for CompareSize  
   // Update OK 
   ndim_ = ndim;
@@ -323,4 +470,5 @@
   marowi_ = a.marowi_;
   veceli_ = a.veceli_;
+  ck_memo_vt_ = a.ck_memo_vt_;
   // Update OK 
   ndim_ = a.ndim_;
Index: /trunk/SophyaLib/TArray/basarr.h
===================================================================
--- /trunk/SophyaLib/TArray/basarr.h	(revision 812)
+++ /trunk/SophyaLib/TArray/basarr.h	(revision 813)
@@ -22,12 +22,23 @@
 public:
   //  To define Array / Matrix memory mapping
-  enum { AutoMemoryMapping = -1, SameMemoryMapping = 0, 
+  enum MemoryMapping { AutoMemoryMapping = -1, SameMemoryMapping = 0, 
 	 CMemoryMapping = 1, FortranMemoryMapping = 2 };
+  //  Vector type
+  enum VectorType {AutoVectorType = -1, SameVectorType = 0, 
+		   ColumnVector = 1, RowVector = 2};
 
-  //    Max Nb of printed elements
-  static void    SetMaxPrint(uint_4 nprt=50);   
-  //   Memory organisation (for matrices)
+  //    Size threshold for parallel routine call
+  static void	 SetOpenMPSizeThreshold(uint_8 thr=200000);
+  static inline uint_8 GetOpenMPSizeThreshold() { return openmp_size_threshold; }
+  //    Max Nb of printed elements and print level
+  static void    SetMaxPrint(uint_4 nprt=50, uint_4 lev=0);   
+  static inline uint_4 GetMaxPrint() { return max_nprt_; }
+  static inline uint_4 GetPrintLevel() { return prt_lev_; }
+
+  //   Memory organisation (for matrices) and vector type
   static short   SetDefaultMemoryMapping(short mm=CMemoryMapping);
   static inline short GetDefaultMemoryMapping() { return default_memory_mapping; }
+  static short   SetDefaultVectorType(short vt=ColumnVector);
+  static inline short GetDefaultVectorType() { return default_vector_type; }
   
   // Creation / destruction 
@@ -53,12 +64,19 @@
   uint_4 MaxSizeKA() const ;
 
-  // memory organisation - packing information
+  // memory organisation 
   inline short  GetMemoryMapping() const 
-                {return ( (marowi_ == 1) ? CMemoryMapping : FortranMemoryMapping) ; }
+                { 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 uint_4 RowsKA() const {return marowi_; }    // Dimension des index de lignes
+  inline uint_4 ColsKA() const {return macoli_; }    // Dimension des index de colonnes
+  inline uint_4 VectKA() const {return veceli_; }    // Dimension des index des elts d'un vecteur
+  void          SetMemoryMapping(short mm=AutoMemoryMapping);
 
+  // Vector type   Line or Column vector
+  inline short  GetVectorType() const 
+                { return((marowi_ == veceli_) ? ColumnVector : RowVector); }
+  void          SetVectorType(short vt=AutoVectorType);
+
+  // memory organisation  - packing information
   inline bool   IsPacked() const { return(moystep_ == 1); }
   inline bool   IsPackedX() const { return(step_[0] == 1); }
@@ -76,5 +94,8 @@
   uint_4 MinStepKA() const ;
 
+  // Offset of element ip
   uint_8 Offset(uint_8 ip=0) const ;
+  // Offset of the i'th vector along axe ka 
+  uint_8  Offset(uint_4 ka, uint_8 i) const ;
   inline uint_8  Offset(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it=0, uint_4 iu=0) const;  
 
@@ -85,5 +106,5 @@
   void           Show(ostream& os, bool si=false) const;
   inline void    Show() const { Show(cout); }
-  virtual string DataType() const = 0;
+  virtual string InfoString() const;
 
 //  Objet DVList info 
@@ -102,5 +123,6 @@
   //  Organisation memoire
   static short SelectMemoryMapping(short mm);
-  void UpdateMemoryMapping(short mm, uint_4 & nx, uint_4 & ny);
+  static short SelectVectorType(short vt);
+  void UpdateMemoryMapping(short mm);
   void UpdateMemoryMapping(BaseArray const & a, short mm);
 
@@ -113,14 +135,17 @@
   uint_4 step_[BASEARRAY_MAXNDIMS];     // two consecutive elements distance in a given dimension 
   uint_4 minstep_;                   // minimal step (in any axes)
-  uint_4 moystep_;                   // mean step 0 non regular steps
+  uint_4 moystep_;                   // mean step if == 0 --> non regular steps
   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)
-
+  bool ck_memo_vt_;            // if true, check MemoryOrg./VectorType for CompareSize
   DVList* mInfo;               // Infos (variables) attachees au tableau
 
   static char * ck_op_msg_[6];  // Operation messages for CheckDI() CheckBound()
   static uint_4 max_nprt_;      // Nb maxi d'elements imprimes
+  static uint_4 prt_lev_;        // Niveau de print 0 ou 1
   static short default_memory_mapping;  // Default memory mapping
+  static short default_vector_type;     // Default vector type Row/Column
+  static uint_8 openmp_size_threshold;  // Size limit for parallel routine calls
 };
 
@@ -132,5 +157,5 @@
   if ( (ka < 0) || (ka >= ndim_) ) {
     string txt = "BaseArray::CheckDimensionIndex/Error ";   txt += ck_op_msg_[msg];
-    throw(ParmError(txt));
+    throw(RangeCheckError(txt));
   }
   return(ka);
@@ -142,8 +167,9 @@
        (it >= size_[3]) || (iu >= size_[4]) ) {
     string txt = "BaseArray::CheckArrayBound/Error ";   txt += ck_op_msg_[msg];
-    throw(ParmError(txt));
+    throw(RangeCheckError(txt));
   }
   return;
 }
+
 
 
Index: /trunk/SophyaLib/TArray/matharr.cc
===================================================================
--- /trunk/SophyaLib/TArray/matharr.cc	(revision 812)
+++ /trunk/SophyaLib/TArray/matharr.cc	(revision 813)
@@ -27,6 +27,7 @@
     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);
+    uint_8 naxa = a.Size()/a.Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = a.DataBlock().Begin()+a.Offset(ka,j);
       for(k=0; k<gpas; k+=step)  pe[k] = (T)(f((double)pe[k]));
     }
@@ -52,6 +53,7 @@
     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);
+    uint_8 naxa = a.Size()/a.Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = a.DataBlock().Begin()+a.Offset(ka,j);
       for(k=0; k<gpas; k+=step)  pe[k] = (T)(f((float)pe[k]));
     }
@@ -102,6 +104,7 @@
     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);
+    uint_8 naxa = a.Size()/a.Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = a.DataBlock().Begin()+a.Offset(ka,j);
       for(k=0; k<gpas; k+=step) { 
 	valok = (double) pe[k]; 
Index: /trunk/SophyaLib/TArray/sopemtx.cc
===================================================================
--- /trunk/SophyaLib/TArray/sopemtx.cc	(revision 812)
+++ /trunk/SophyaLib/TArray/sopemtx.cc	(revision 813)
@@ -63,4 +63,5 @@
 
   uint_4 IMaxAbs(uint_4 first=0);
+  void Print(ostream & os) const ;
 
   static void Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2);
@@ -105,6 +106,6 @@
 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());
+  { switch (rckind) { case TmatrixRow  : return m.Step(m.ColsKA());
+                      case TmatrixCol  : return m.Step(m.RowsKA());
                       case TmatrixDiag : return (m.Step(m.RowsKA())+m.Step(m.ColsKA())); }
     return 0; }
@@ -333,4 +334,17 @@
 
 template <class T>
+void TMatrixRC<T>::Print(ostream & os) const
+{
+  os << " TMatrixRC<T>::Print(ostream & os) " << NElts() << " Kind=" 
+     << kind << " Index=" << index << " Step= " << step << endl;
+  for(uint_4 i=0; i<NElts(); i++) {
+    os << (*this)(i);
+    if (kind == TmatrixCol)  os << endl;
+    else os << ", "; 
+  }
+  os << endl;
+}
+
+template <class T>
 void TMatrixRC<T>::Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2)
 {
@@ -435,5 +449,5 @@
 
 for(uint_4 l=0; l<n; l++) {
-  if (fabs(a(l,l)) <= 1.e-50) return 0.0;
+  if (fabs((double)a(l,l)) <= 1.e-50) return 0.0;
   TMatrixRC<T>::Row(b, l) /= a(l,l);
 }
@@ -447,6 +461,6 @@
 {
 TMatrix<T> a(A);
-TMatrix<T> b(a.NCols(),a.NRows());  b = 1.;
-if(fabs(GausPiv(a,b)) < 1.e-50)
+TMatrix<T> b(a.NCols(),a.NRows());  b = IdentityMatrix(1.);
+if(fabs((double)GausPiv(a,b)) < 1.e-50)
   throw(MathExc("TMatrix Inverse() Singular OMatrix"));
 return b;
@@ -578,5 +592,31 @@
 
 
-
+void _ZZ_TestTMatrixRC_YY_(TMatrix<r_8> & m)
+{
+  cout << " ::::: _ZZ_TestTMatrixRC_YY_ :::: M= \n" << m << endl;
+  TMatrixRC<r_8> l0 = TMatrixRC<r_8>::Row(m,0);
+  cout << "TMatrixRC<r_8>::Row(m,0) = \n" ;
+  l0.Print(cout);
+  TMatrixRC<r_8> l1 = TMatrixRC<r_8>::Row(m,1);
+  cout << "TMatrixRC<r_8>::Row(m,1) = \n" ;
+  l1.Print(cout);
+  l0.SetRow(2);
+  cout << "TMatrixRC<r_8>::l0.SetRow(2 = \n" ;
+  l0.Print(cout);
+
+  TMatrixRC<r_8> c0 = TMatrixRC<r_8>::Col(m,0);
+  cout << "TMatrixRC<r_8>::Col(m,0) = \n" ;
+  c0.Print(cout);
+  TMatrixRC<r_8> c1 = TMatrixRC<r_8>::Col(m,1);
+  cout << "TMatrixRC<r_8>::Col(m,1) = \n" ;
+  c1.Print(cout);  
+  c0.SetCol(2);
+  cout << "TMatrixRC<r_8>::c0.SetCol(2) = \n" ;
+  c0.Print(cout);
+  TMatrixRC<r_8> c00 = TMatrixRC<r_8>::Col(m,0);
+  TMatrixRC<r_8>::Swap(c0, c00);
+  cout << " :::::  M Apres Swap = \n" << m << endl;
+  
+}
 
 ///////////////////////////////////////////////////////////////
Index: /trunk/SophyaLib/TArray/tarray.cc
===================================================================
--- /trunk/SophyaLib/TArray/tarray.cc	(revision 812)
+++ /trunk/SophyaLib/TArray/tarray.cc	(revision 813)
@@ -234,8 +234,8 @@
 //  Possibilite de // , vectorisation 
 template <class T>
-TArray<T>& TArray<T>::Set(Sequence seq)
-{
-  if (NbDimensions() < 1) 
-    throw RangeCheckError("TArray<T>::Set(Sequence ) - Not Allocated Array ! ");
+TArray<T>& TArray<T>::SetSeq(Sequence seq)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::SetSeq(Sequence ) - Not Allocated Array ! ");
   T * pe;
   uint_8 j,k;
@@ -246,10 +246,12 @@
   }
   else {    // Non regular data spacing ...
-    uint_4 ka = MaxSizeKA();
+    //    uint_4 ka = MaxSizeKA();
+    uint_4 ka = 0;
     uint_8 step = Step(ka);
     uint_8 gpas = Size(ka);
-    for(j=0; j<totsize_; j += Size(ka))  {
-      pe = mNDBlock.Begin()+Offset(j);
-      for(k=0; k<gpas; k++)  pe[k*step] = seq(j+k);
+    uint_8 naxa = Size()/Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ka,j);
+      for(k=0; k<gpas; k++)  pe[k*step] = seq(j*gpas+k);
     }
   }
@@ -260,8 +262,8 @@
 
 template <class T>
-TArray<T>& TArray<T>::Set(T x)
-{
-  if (NbDimensions() < 1) 
-    throw RangeCheckError("TArray<T>::Set(T )  - Not Allocated Array ! ");
+TArray<T>& TArray<T>::SetT(T x)
+{
+  if (NbDimensions() < 1) 
+    throw RangeCheckError("TArray<T>::SetT(T )  - Not Allocated Array ! ");
   T * pe;
   uint_8 j,k;
@@ -276,6 +278,7 @@
     uint_8 step = Step(ka);
     uint_8 gpas = Size(ka)*step;
-    for(j=0; j<totsize_; j += Size(ka))  {
-      pe = mNDBlock.Begin()+Offset(j);
+    uint_8 naxa = Size()/Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ka,j);
       for(k=0; k<gpas; k+=step)  pe[k] = x;
     }
@@ -301,6 +304,7 @@
     uint_8 step = Step(ka);
     uint_8 gpas = Size(ka)*step;
-    for(j=0; j<totsize_; j += Size(ka))  {
-      pe = mNDBlock.Begin()+Offset(j);
+    uint_8 naxa = Size()/Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ka,j);
       for(k=0; k<gpas; k+=step)  pe[k] += x;
     }
@@ -326,6 +330,7 @@
     uint_8 step = Step(ka);
     uint_8 gpas = Size(ka)*step;
-    for(j=0; j<totsize_; j += Size(ka))  {
-      pe = mNDBlock.Begin()+Offset(j);
+    uint_8 naxa = Size()/Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ka,j);
       for(k=0; k<gpas; k+=step)  pe[k] -= x;
     }
@@ -351,6 +356,7 @@
     uint_8 step = Step(ka);
     uint_8 gpas = Size(ka)*step;
-    for(j=0; j<totsize_; j += Size(ka))  {
-      pe = mNDBlock.Begin()+Offset(j);
+    uint_8 naxa = Size()/Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ka,j);
       for(k=0; k<gpas; k+=step)  pe[k] *= x;
     }
@@ -376,6 +382,7 @@
     uint_8 step = Step(ka);
     uint_8 gpas = Size(ka)*step;
-    for(j=0; j<totsize_; j += Size(ka))  {
-      pe = mNDBlock.Begin()+Offset(j);
+    uint_8 naxa = Size()/Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ka,j);
       for(k=0; k<gpas; k+=step)  pe[k] /= x;
     }
@@ -402,6 +409,7 @@
     uint_8 step = Step(ka);
     uint_8 gpas = Size(ka)*step;
-    for(j=0; j<totsize_; j += Size(ka))  {
-      pe = mNDBlock.Begin()+Offset(j);
+    uint_8 naxa = Size()/Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ka,j);
       for(k=0; k<gpas; k+=step)  pe[k] = x-pe[k];
     }
@@ -427,6 +435,7 @@
     uint_8 step = Step(ka);
     uint_8 gpas = Size(ka)*step;
-    for(j=0; j<totsize_; j += Size(ka))  {
-      pe = mNDBlock.Begin()+Offset(j);
+    uint_8 naxa = Size()/Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ka,j);
       for(k=0; k<gpas; k+=step)  pe[k] = x/pe[k];
     }
@@ -442,5 +451,6 @@
   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>&)")) ;
+  if (!CompareSizes(a)) 
+    throw(SzMismatchError("TArray<T>::AddElt(const TArray<T>&) SizeMismatch")) ;
 
   T * pe;
@@ -460,7 +470,8 @@
     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);
+    uint_8 naxa = Size()/Size(ax);
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ax,j);
+      pea = a.DataBlock().Begin()+a.Offset(ax,j);
       for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] += pea[ka]; 
     }
@@ -474,5 +485,6 @@
   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>&)")) ;
+  if (!CompareSizes(a)) 
+    throw(SzMismatchError("TArray<T>::SubElt(const TArray<T>&) SizeMismatch")) ;
 
   T * pe;
@@ -492,7 +504,8 @@
     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);
+    uint_8 naxa = Size()/Size(ax);
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ax,j);
+      pea = a.DataBlock().Begin()+a.Offset(ax,j);
       for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] -= pea[ka]; 
     }
@@ -507,5 +520,6 @@
   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>&)")) ;
+  if (!CompareSizes(a)) 
+    throw(SzMismatchError("TArray<T>::MulElt(const TArray<T>&) SizeMismatch")) ;
 
   T * pe;
@@ -525,7 +539,8 @@
     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);
+    uint_8 naxa = Size()/Size(ax);
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ax,j);
+      pea = a.DataBlock().Begin()+a.Offset(ax,j);
       for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] *= pea[ka]; 
     }
@@ -540,5 +555,6 @@
   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>&)")) ;
+  if (!CompareSizes(a)) 
+    throw(SzMismatchError("TArray<T>::DivElt(const TArray<T>&) SizeMismatch")) ;
 
   T * pe;
@@ -558,7 +574,8 @@
     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);
+    uint_8 naxa = Size()/Size(ax);
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ax,j);
+      pea = a.DataBlock().Begin()+a.Offset(ax,j);
       for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] /= pea[ka]; 
     }
@@ -572,5 +589,6 @@
   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>&)")) ;
+  if (!CompareSizes(a)) 
+    throw(SzMismatchError("TArray<T>::MultElt(const TArray<T>&) SizeMismatch")) ;
 
   T * pe;
@@ -590,7 +608,8 @@
     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);
+    uint_8 naxa = Size()/Size(ax);
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ax,j);
+      pea = a.DataBlock().Begin()+a.Offset(ax,j);
       for(k=0, ka=0;  k<gpas;  k+=step, ka+=stepa)  pe[k] = pea[ka]; 
     }
@@ -619,6 +638,7 @@
     uint_8 step = Step(ka);
     uint_8 gpas = Size(ka)*step;
-    for(j=0; j<totsize_; j += Size(ka))  {
-      pe = mNDBlock.Begin()+Offset(j);
+    uint_8 naxa = Size()/Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ka,j);
       for(k=0; k<gpas; k+=step)  ret += pe[k] ;
     }
@@ -645,6 +665,7 @@
     uint_8 step = Step(ka);
     uint_8 gpas = Size(ka)*step;
-    for(j=0; j<totsize_; j += Size(ka))  {
-      pe = mNDBlock.Begin()+Offset(j);
+    uint_8 naxa = Size()/Size(ka);
+    for(j=0; j<naxa; j++)  {
+      pe = mNDBlock.Begin()+Offset(ka,j);
       for(k=0; k<gpas; k+=step)  ret *= pe[k] ;
     }
@@ -660,7 +681,9 @@
 
 template <class T>
-string TArray<T>::DataType() const
-{
-  string rs = typeid(T).name();
+string TArray<T>::InfoString() const
+{
+  string rs = "TArray<" ;
+  rs += typeid(T).name();
+  rs += "> ";
   return(rs);
 }
@@ -693,5 +716,5 @@
     }
   }
-  os << "\n" << endl;
+  os <<  endl;
 }
 
Index: /trunk/SophyaLib/TArray/tarray.h
===================================================================
--- /trunk/SophyaLib/TArray/tarray.h	(revision 812)
+++ /trunk/SophyaLib/TArray/tarray.h	(revision 813)
@@ -88,9 +88,9 @@
   inline T toScalar();
 // Met les elements a une suite de valeurs
-  virtual TArray<T>&  Set(Sequence seq);
-  inline  TArray<T>&  operator = (Sequence seq)    { return Set(seq); }
+  virtual TArray<T>&  SetSeq(Sequence seq);
+  inline  TArray<T>&  operator = (Sequence seq)    { return SetSeq(seq); }
 // A = x (tous les elements a x)
-  virtual TArray<T>&  Set(T x);
-  inline  TArray<T>&  operator = (T x)             { return Set(x); }
+  virtual TArray<T>&  SetT(T x);
+  inline  TArray<T>&  operator = (T x)             { return SetT(x); }
 // A += -= *= /= x (ajoute, soustrait, ... x a tous les elements)
   virtual TArray<T>&  Add(T x);
@@ -121,6 +121,6 @@
 
 // Impression, I/O, ...
-  virtual string DataType() const;   // definition of abstract method inherited from BaseArray
-  virtual void  Print(ostream& os, int_4 maxprt=-1, bool si=false) const ;
+  virtual string InfoString() const;   
+  virtual void   Print(ostream& os, int_4 maxprt=-1, bool si=false) const ;
 
 //  Pour la gestion de persistance
Index: /trunk/SophyaLib/TArray/tarrinit.cc
===================================================================
--- /trunk/SophyaLib/TArray/tarrinit.cc	(revision 812)
+++ /trunk/SophyaLib/TArray/tarrinit.cc	(revision 813)
@@ -21,10 +21,20 @@
   PPRegister(FIO_TArray<uint_2>);
   DObjRegister(FIO_TArray<uint_2>, TArray<uint_2>);
+  DObjRegister(FIO_TArray<uint_2>, TMatrix<uint_2>);
+  DObjRegister(FIO_TArray<uint_2>, TVector<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>);
+  DObjRegister(FIO_TArray<int_4>, TMatrix<int_4>);
+  DObjRegister(FIO_TArray<int_4>, TVector<int_4>);
+
   PPRegister(FIO_TArray<int_8>);
   DObjRegister(FIO_TArray<int_8>, TArray<int_8>);
+  DObjRegister(FIO_TArray<int_8>, TMatrix<int_8>);
+  DObjRegister(FIO_TArray<int_8>, TVector<int_8>);
+
   //  PPRegister(FIO_TArray<uint_4>);
   // DObjRegister(FIO_TArray<uint_4>, TArray<uint_4>);
@@ -33,10 +43,21 @@
   PPRegister(FIO_TArray<r_4>);
   DObjRegister(FIO_TArray<r_4>, TArray<r_4>);
+  DObjRegister(FIO_TArray<r_4>, TMatrix<r_4>);
+  DObjRegister(FIO_TArray<r_4>, TVector<r_4>);
+
   PPRegister(FIO_TArray<r_8>);
   DObjRegister(FIO_TArray<r_8>, TArray<r_8>);
+  DObjRegister(FIO_TArray<r_8>, TMatrix<r_8>);
+  DObjRegister(FIO_TArray<r_8>, TVector<r_8>);
+
   PPRegister(FIO_TArray< complex<r_4> >);
   DObjRegister(FIO_TArray< complex<r_4> >, TArray< complex<r_4> >);
+  DObjRegister(FIO_TArray< complex<r_4> >, TMatrix< complex<r_4> >);
+  DObjRegister(FIO_TArray< complex<r_4> >, TVector< complex<r_4> >);
+
   PPRegister(FIO_TArray< complex<r_8> >);
   DObjRegister(FIO_TArray< complex<r_8> >, TArray< complex<r_8> >);
+  DObjRegister(FIO_TArray< complex<r_8> >, TMatrix< complex<r_8> >);
+  DObjRegister(FIO_TArray< complex<r_8> >, TVector< complex<r_8> >);
 
 }
Index: /trunk/SophyaLib/TArray/tmatrix.cc
===================================================================
--- /trunk/SophyaLib/TArray/tmatrix.cc	(revision 812)
+++ /trunk/SophyaLib/TArray/tmatrix.cc	(revision 813)
@@ -1,3 +1,3 @@
-// $Id: tmatrix.cc,v 1.3 2000-04-03 17:35:58 ansari Exp $
+// $Id: tmatrix.cc,v 1.4 2000-04-05 15:44:15 ansari Exp $
 //                         C.Magneville          04/99
 #include "machdefs.h"
@@ -17,4 +17,5 @@
   : TArray<T>()
 {
+  ck_memo_vt_ = true;
 }
 
@@ -48,6 +49,12 @@
 : TArray<T>(a)
 {
-  if (a.NbDimensions() != 2) 
-    throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a) a.NbDimensions() != 2 ");
+  if (a.NbDimensions() > 2) 
+    throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a) a.NbDimensions()>2 ");
+  if (a.NbDimensions() == 1) {
+    size_[1] = 1;
+    step_[1] = size_[0]*step_[0];
+    ndim_ = 2;
+  }
+  UpdateMemoryMapping(a, SameMemoryMapping);
 }
 
@@ -56,6 +63,11 @@
 : TArray<T>(a, share)
 {
-  if (a.NbDimensions() != 2) 
-    throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a, ...) a.NbDimensions() != 2 ");
+  if (a.NbDimensions() > 2) 
+    throw SzMismatchError("TMatrix<T>::TMatrix(const TArray<T>& a, ...) a.NbDimensions()>2");
+  if (a.NbDimensions() == 1) {
+    size_[1] = 1;
+    step_[1] = size_[0]*step_[0];
+    ndim_ = 2;
+  }
   UpdateMemoryMapping(a, mm);
 }
@@ -71,7 +83,14 @@
 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));
+  if (a.NbDimensions() > 2) 
+    throw SzMismatchError("TMatrix<T>::Set(const TArray<T>& a) a.NbDimensions() > 2");
+  TArray<T>::Set(a);
+  if (a.NbDimensions() == 1) {
+    size_[1] = 1;
+    step_[1] = size_[0]*step_[0];
+    ndim_ = 2;
+  }
+  UpdateMemoryMapping(a, SameMemoryMapping);
+  return(*this);
 }
 
@@ -83,9 +102,15 @@
   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]);
+  else if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) 
+    mm = GetDefaultMemoryMapping();
+  if (mm == CMemoryMapping) {
+    size[0] = c;  size[1] = r;
+  }
+  else {
+    size[0] = r;  size[1] = c;
+  }
   TArray<T>::ReSize(2, size, 1);
-  UpdateMemoryMapping(mm, size[0], size[1]);
+  UpdateMemoryMapping(mm);
 }
 
@@ -97,27 +122,29 @@
   uint_4 size[BASEARRAY_MAXNDIMS];
   for(int kk=0; kk<BASEARRAY_MAXNDIMS; kk++)  size[kk] = 0;
-  UpdateMemoryMapping(mm, size[0], size[1]);
+  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
+  else if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) 
+    mm = GetDefaultMemoryMapping();
+  if (mm == CMemoryMapping) {
+    size[0] = c;  size[1] = r;
+  }
+  else {
+    size[0] = r;  size[1] = c;
+  }
   TArray<T>::Realloc(2, size, 1, force);
-  UpdateMemoryMapping(mm, size[0], size[1]);
+  UpdateMemoryMapping(mm);
 }
 
 // $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);
-  }
+TMatrix<T> TMatrix<T>::SubMatrix(Range rline, Range rcol) const
+{
+  short mm = GetMemoryMapping();
+  Range rx, ry;
+  if (mm == CMemoryMapping)  { rx = rcol;  ry = rline; }
+  else { ry = rcol;  rx = rline; }
+  TMatrix sm(SubArray(rx, ry, Range(0), Range(0), Range(0)),true, mm);
+  sm.UpdateMemoryMapping(mm);
+  sm.SetTemp(true);
+  return(sm);
 }
 
@@ -127,7 +154,9 @@
 TMatrix<T>& TMatrix<T>::Transpose() 
 {
+  short vt = (marowi_ == veceli_) ? ColumnVector : RowVector;
   uint_4 rci = macoli_;
   macoli_ = marowi_;
   marowi_ = rci;
+  veceli_ = (vt ==  ColumnVector ) ?  marowi_ : macoli_;
   return(*this);
 }
@@ -147,8 +176,15 @@
 
 // Rearrangement memoire
+// Si identique, renvoie une matrice qui partage les donnees 
 template <class T>
 TMatrix<T> TMatrix<T>::Rearrange(short mm)
 {
-  if (mm == SameMemoryMapping) mm = GetMemoryMapping();  
+  if ( mm == SameMemoryMapping)  mm = GetMemoryMapping();
+  else if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) 
+    mm = GetDefaultMemoryMapping();
+  
+  if  (mm == GetMemoryMapping())
+    return (TMatrix<T>(*this, true));
+  
   TMatrix<T> tm(NRows(), NCols(), mm);
   for(uint_4 i=0; i<NRows(); i++)
@@ -180,9 +216,17 @@
 //**** Impression
 template <class T>
+string TMatrix<T>::InfoString() const
+{
+  string rs = "TMatrix<";
+  rs += typeid(T).name();
+  char buff[64];
+  sprintf(buff, ">(NRows=%ld, NCols=%ld)", (long)NRows(), (long)NCols());
+  rs += buff;
+  return(rs);  
+}
+
+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;
@@ -200,5 +244,5 @@
     os << endl;
   }
-  os << "\n" << endl;
+  os << endl;
 }
 
@@ -235,5 +279,4 @@
 }
 
-
 ///////////////////////////////////////////////////////////////
 #ifdef __CXX_PRAGMA_TEMPLATES__
Index: /trunk/SophyaLib/TArray/tmatrix.h
===================================================================
--- /trunk/SophyaLib/TArray/tmatrix.h	(revision 812)
+++ /trunk/SophyaLib/TArray/tmatrix.h	(revision 813)
@@ -34,6 +34,13 @@
   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 ;
+  // Sub-matrix extraction $CHECK$ Reza 03/2000  Doit-on declarer ces methode const ?
+  TMatrix<T> SubMatrix(Range rline, Range rcol) const ;
+  inline TMatrix<T> operator () (Range rline, Range rcol) const 
+                    { return SubMatrix(rline, rcol); }
+  // Lignes et colonnes de la matrice
+  inline TMatrix<T> Row(uint_4 ir) const
+                    { return SubMatrix(Range(ir,ir), Range(0,NCols()-1)); }
+  inline TMatrix<T> Column(uint_4 ic) const
+                    { return SubMatrix(Range(0,NRows()-1), Range(ic,ic)); }
 
   // Inline element acces methods 
@@ -53,6 +60,8 @@
   inline  TMatrix<T>& operator = (IdentityMatrix imx) { return SetIdentity(imx); }
 
+  inline  TMatrix<T>&  operator = (Sequence seq)    { SetSeq(seq); return(*this); }
+
   // Operations diverses  avec une constante
-  inline  TMatrix<T>&  operator = (T x)             { Set(x); return(*this); }
+  inline  TMatrix<T>&  operator = (T x)             { SetT(x); return(*this); }
   inline  TMatrix<T>&  operator += (T x)            { Add(x); return(*this); }
   inline  TMatrix<T>&  operator -= (T x)            { Sub(x); return(*this); }
@@ -61,13 +70,14 @@
 
   //  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); }
+  inline  TMatrix<T>&  operator += (const TMatrix<T>& a)  { AddElt(a); return(*this); }
+  inline  TMatrix<T>&  operator -= (const TMatrix<T>& a)  { SubElt(a); return(*this); }
+  // Produit matriciel Multiply : C = (*this)*B
+  TMatrix<T>  Multiply(const TMatrix<T>& b, short mm=SameMemoryMapping) const;
+  inline  TMatrix<T>&  operator *= (const TMatrix<T>& b)
+          { this->Set(Multiply(b));  return(*this); }
 
-  // Input-Output
+  // I/O print, ...
+  virtual string InfoString() const;
   virtual void  Print(ostream& os, int_4 maxprt=-1, bool si=false) const ;
-
-  // Produit matriciel 
-  TMatrix<T>  Multiply(const TMatrix<T>& b, short mm=SameMemoryMapping) const;
-  //  inline TMatrix<T>& operator *= (const TMatrix<T>& b) 
 
 protected:
@@ -98,8 +108,21 @@
 
 
+// Surcharge d'operateurs C = A (+,-) B
+// $CHECK$ Reza 3/4/2000 Pas necessaire  de redefinir les operateurs
+// Defini au niveau de TArray<T> - Pour ameliorer l'efficacite
+// Doit-on le faire aussi pour les constantes ? - Fin de $CHECK$ Reza 3/4/2000
+
+template <class T>
+inline TMatrix<T> operator + (const TMatrix<T>& a,const TMatrix<T>& b)
+    {TMatrix<T> result(a); result.SetTemp(true); result.AddElt(b); return result;}
+
+template <class T>
+inline TMatrix<T> operator - (const TMatrix<T>& a,const TMatrix<T>& b)
+    {TMatrix<T> result(a); result.SetTemp(true); result.SubElt(b); return result;}
+
 // Surcharge d'operateurs C = A * B
 
 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;}
+{ TMatrix<T> result(a); result.SetTemp(true); return(result.Multiply(b)); }
 
 typedef TMatrix<r_8> Matrix;
Index: /trunk/SophyaLib/TArray/tvector.cc
===================================================================
--- /trunk/SophyaLib/TArray/tvector.cc	(revision 812)
+++ /trunk/SophyaLib/TArray/tvector.cc	(revision 813)
@@ -1,3 +1,3 @@
-// $Id: tvector.cc,v 1.2 2000-04-03 17:35:59 ansari Exp $
+// $Id: tvector.cc,v 1.3 2000-04-05 15:44:17 ansari Exp $
 //                         C.Magneville          04/99
 #include "machdefs.h"
@@ -18,8 +18,8 @@
 TVector<T>::TVector(uint_4 n, short lcv, short mm)
 // Constructeur 
-  : TMatrix<T>((lcv == ColumnVector) ? n : 1,(lcv == ColumnVector) ? 1 : 0, mm)
+  : TMatrix<T>(1,1,mm)
 {
-  if (lcv == ColumnVector) veceli_ = marowi_;
-  else veceli_ = macoli_;
+  lcv = SelectVectorType(lcv);
+  ReSize(n,lcv);
 }
 
@@ -43,15 +43,19 @@
 : TMatrix<T>(a)
 {
+  if ( (size_[0] != 1) && (size_[1] != 1) )
+    throw SzMismatchError("TVector<T>::TVector(const TArray<T>& a) NRows()!=1 && NCols()!=1 ");
 }
 
 
 template <class T>
-TVector<T>::TVector(const TArray<T>& a, bool share, short lcv, short mm )
+TVector<T>::TVector(const TArray<T>& a, bool share, short mm, short lcv )
 : TMatrix<T>(a, share, mm)
 {
-  if (lcv == SameTypeVector) veceli_ = a.VectKA();
-  else {
-    if (lcv == ColumnVector) veceli_ = marowi_;
-    else veceli_ = macoli_;
+  if ( (size_[0] != 1) && (size_[1] != 1) )
+    throw SzMismatchError("TVector<T>::TVector(const TArray<T>& a) NRows()!=1 && NCols()!=1 ");
+  if ( (size_[0] == 1) && (size_[1] == 1) ) {
+    if (lcv == SameVectorType) lcv = a.GetVectorType();
+    if ( (lcv != ColumnVector) && (lcv != RowVector) ) lcv = GetDefaultVectorType();
+    veceli_ = (lcv ==  ColumnVector ) ?  marowi_ : macoli_;
   }
 }
@@ -70,8 +74,10 @@
     throw(SzMismatchError("TVector::ReSize() n = 0 "));
   uint_4 r,c;
-  if (lcv == SameTypeVector)  lcv = GetVectorType();
+  if (lcv == SameVectorType)  lcv = GetVectorType();
+  else if ( (lcv != ColumnVector) && (lcv != RowVector) ) lcv = GetDefaultVectorType();
   if (lcv == ColumnVector) { r = n;  c = 1; }
   else { c = n; r = 1; }
   TMatrix<T>::ReSize(r,c);
+  veceli_ = (lcv ==  ColumnVector ) ?  marowi_ : macoli_;
 }
 
@@ -82,13 +88,15 @@
     throw(SzMismatchError("TVector::Realloc() n = 0 "));
   uint_4 r,c;
-  if (lcv == SameTypeVector)  lcv = GetVectorType();
+  if (lcv == SameVectorType)  lcv = GetVectorType();
+  else if ( (lcv != ColumnVector) && (lcv != RowVector) ) lcv = GetDefaultVectorType();
   if (lcv == ColumnVector) { r = n;  c = 1; }
   else { c = n; r = 1; }
   TMatrix<T>::Realloc(r,c,SameMemoryMapping,force);
+  veceli_ = (lcv ==  ColumnVector ) ?  marowi_ : macoli_;
 }
 
 // $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
 template <class T>
-TVector<T> TVector<T>::operator () (Range relt) const
+TVector<T> TVector<T>::SubVector(Range relt) const
 {
   Range rr, cr;
@@ -109,4 +117,15 @@
 }
 
+template <class T>
+string TVector<T>::InfoString() const
+{
+  string rs = "TVector<";
+  rs += typeid(T).name();
+  char buff[64];
+  sprintf(buff, ">(%ld) (nr=%ld, nc=%ld)", (long)NElts(), (long)NRows(), (long)NCols());
+  rs += buff;
+  return(rs);  
+
+}
 
 ///////////////////////////////////////////////////////////////
Index: /trunk/SophyaLib/TArray/tvector.h
===================================================================
--- /trunk/SophyaLib/TArray/tvector.h	(revision 812)
+++ /trunk/SophyaLib/TArray/tvector.h	(revision 813)
@@ -11,5 +11,4 @@
 class TVector : public TMatrix<T> {
 public:
-  enum {SameTypeVector=0, ColumnVector=1, LineVector=2};
   // Creation / destruction 
   TVector();
@@ -18,5 +17,5 @@
   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);
+  TVector(const TArray<T>& a,  bool share, short mm=CMemoryMapping, short lcv=ColumnVector);
 
   virtual ~TVector();
@@ -25,14 +24,13 @@
                        { Set(a);  return(*this); }
 
-  // Vector type   Line or Column vector
-  inline short GetVectorType() const 
-               { return((marowi_ == veceli_) ? ColumnVector : LineVector); }
 
   // Gestion taille/Remplissage
-  void ReSize(uint_4 n, short lcv=SameTypeVector ); 
-  void Realloc(uint_4 n, short lcv=SameTypeVector, bool force=false);
+  void ReSize(uint_4 n, short lcv=SameVectorType ); 
+  void Realloc(uint_4 n, short lcv=SameVectorType, bool force=false);
 
   // Sub-Vector extraction $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
-  TVector<T> operator () (Range relt) const ;
+  TVector<T> SubVector(Range relt) const ;
+  inline TVector<T> operator () (Range relt) const 
+                    { return SubVector(relt); }
 
   // Informations pointeur/data
@@ -43,6 +41,9 @@
   inline T&       operator()(uint_4 n);
 
+  // Operateur d'affectation
+  inline  TMatrix<T>&  operator = (Sequence seq)    { SetSeq(seq); return(*this); }
+
   // Operations diverses  avec une constante
-  inline  TVector<T>&  operator = (T x)             { Set(x); return(*this); }
+  inline  TVector<T>&  operator = (T x)             { SetT(x); return(*this); }
   inline  TVector<T>&  operator += (T x)            { Add(x); return(*this); }
   inline  TVector<T>&  operator -= (T x)            { Sub(x); return(*this); }
@@ -56,4 +57,7 @@
   // Norme(^2) 
   T Norm2() const ;
+
+  virtual string InfoString() const;
+
 };
 
Index: /trunk/SophyaLib/TArray/utilarr.cc
===================================================================
--- /trunk/SophyaLib/TArray/utilarr.cc	(revision 812)
+++ /trunk/SophyaLib/TArray/utilarr.cc	(revision 813)
@@ -20,9 +20,17 @@
 }
 
-Range::Range(uint_4 start, uint_4 size, uint_4 step)
+Range::Range(uint_4 start, uint_4 end, uint_4 size, uint_4 step)
 {
   start_ = start;
-  size_ = (size > 0) ? size : 1;
   step_ = (step > 0) ? step : 1;
+  if (end > start) {  // Taille calcule automatiquement 
+    end_ = end;
+    if (step_ > ((end_-start_)+1))  size_ = 1;
+    else size_ = ((end-start)+1)/step_;
+  }
+  else {     // Taille fixee
+    size_ = size;
+    end_ = start_+size_*step_;
+  }
 }
 
Index: /trunk/SophyaLib/TArray/utilarr.h
===================================================================
--- /trunk/SophyaLib/TArray/utilarr.h	(revision 812)
+++ /trunk/SophyaLib/TArray/utilarr.h	(revision 813)
@@ -29,10 +29,11 @@
 class Range {
 public:
-  explicit Range(uint_4 start=0, uint_4 size=1, uint_4 step=1);
+  explicit Range(uint_4 start=0, uint_4 end=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_; }
+  inline uint_4 & End()    {  return end_; }
+  inline uint_4 & Size()   {  return size_; }
+  inline uint_4 & Step()   {  return step_; }
 protected:
-  uint_4 start_, size_, step_;
+  uint_4 start_, end_, size_, step_ ;
 };
 
