Index: /trunk/SophyaLib/TArray/fioarr.cc
===================================================================
--- /trunk/SophyaLib/TArray/fioarr.cc	(revision 772)
+++ /trunk/SophyaLib/TArray/fioarr.cc	(revision 772)
@@ -0,0 +1,143 @@
+//     Persistence manager for template numerical arrays
+//               R. Ansari, C.Magneville   03/2000
+
+#include "pexceptions.h"
+#include "fioarr.h"
+
+// --------------------------------------------------------
+//   Les objets delegues pour la gestion de persistance 
+// --------------------------------------------------------
+///////////////////////////////////////////////////////////
+
+template <class T>
+FIO_TArray<T>::FIO_TArray()
+{
+  dobj=new TArray<T>;
+  ownobj=true;
+}
+
+template <class T>
+FIO_TArray<T>::FIO_TArray(string const & filename) 
+{
+  dobj=new TArray<T>;
+  ownobj=true; 
+  Read(filename);
+}
+
+template <class T>
+FIO_TArray<T>::FIO_TArray(const TArray<T> & obj) 
+{ 
+  dobj = new TArray<T>(obj);
+  ownobj=true; 
+}
+
+template <class T>
+FIO_TArray<T>::FIO_TArray(TArray<T> * obj) 
+{ 
+  dobj = obj;
+  ownobj=false; 
+}
+
+template <class T>
+FIO_TArray<T>::~FIO_TArray()
+{
+  if (ownobj && dobj) delete dobj;
+}
+
+template <class T>
+AnyDataObj* FIO_TArray<T>::DataObj()
+{
+  return(dobj);
+}
+
+template <class T>
+void FIO_TArray<T>::SetDataObj(AnyDataObj & o)
+{
+  TArray<T> * po = dynamic_cast< TArray<T> * >(&o);
+  if (po == NULL) return;
+  if (ownobj && dobj) delete dobj;
+  dobj = po; ownobj = false;
+} 
+
+template <class T>
+void FIO_TArray<T>::ReadSelf(PInPersist& is)
+{
+// On lit les 3 premiers uint_4
+//  0: Numero de version,   : NRows,  2 : NCol
+  uint_4 itab[5];
+  is.Get(itab,5);
+  if (dobj == NULL) dobj = new TArray<T>;
+// On lit les tailles, etc ...
+  is.Get(dobj->ndim_);
+  is.Get(dobj->size_, TARRAY_MAXNDIMS);
+  is.Get(dobj->totsize_);
+  is.Get(dobj->step_, TARRAY_MAXNDIMS);
+  is.Get(dobj->minstep_);
+  is.Get(dobj->offset_);
+  is.Get(dobj->marowi_);
+  is.Get(dobj->macoli_);
+// On lit le datablock
+  is >> dobj->DataBlock();
+// On ecrit le DVList info si necessaire
+  if (itab[2] != 0)  is >> dobj->Info();  
+}
+
+template <class T>
+void FIO_TArray<T>::WriteSelf(POutPersist& os) const
+{
+  if (dobj == NULL)   return;
+//  On ecrit 4 uint_4 .... 
+//  0: Numero de version,  1 : Type (Array, matrix, Vector, ...)  2 != 0 , has Info
+  uint_4 itab[5];
+  itab[0] = 1;  // Numero de version a 1
+  itab[1] = 0;
+  itab[2] = (dobj->mInfo != NULL) ? 1 : 0;
+  itab[3] = itab[4] = 0;
+  os.Put(itab,3);
+// On ecrit les tailles, etc ...
+  os.Put(dobj->ndim_);
+  os.Put(dobj->size_, TARRAY_MAXNDIMS);
+  os.Put(dobj->totsize_);
+  os.Put(dobj->step_, TARRAY_MAXNDIMS);
+  os.Put(dobj->minstep_);
+  os.Put(dobj->offset_);
+  os.Put(dobj->marowi_);
+  os.Put(dobj->macoli_);
+// On ecrit le datablock
+  os << dobj->DataBlock();
+// On ecrit le DVList info si necessaire
+  if (itab[2] != 0)  os << dobj->Info();
+}
+
+
+
+///////////////////////////////////////////////////////////////
+#ifdef __CXX_PRAGMA_TEMPLATES__
+// Instances des delegues FileIO (PPersist)
+#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_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<r_8>
+#pragma define_template FIO_TArray<r_4>
+#pragma define_template FIO_TArray< complex<r_4> >
+#pragma define_template FIO_TArray< complex<r_8> >
+#endif
+
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+// Instances des delegues FileIO (PPersist)
+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_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>;
+template class FIO_TArray< complex<r_4> >;
+template class FIO_TArray< complex<r_8> >;
+#endif
Index: /trunk/SophyaLib/TArray/fioarr.h
===================================================================
--- /trunk/SophyaLib/TArray/fioarr.h	(revision 772)
+++ /trunk/SophyaLib/TArray/fioarr.h	(revision 772)
@@ -0,0 +1,43 @@
+// This may look like C code, but it is really -*- C++ -*-
+//      template array class for numerical types
+//                     R. Ansari, C.Magneville   03/2000
+
+#ifndef FIOTArray_SEEN
+#define FIOTArray_SEEN
+
+#include "machdefs.h"
+#include <iostream.h>
+#include "tarray.h"
+#include "ppersist.h"
+
+namespace SOPHYA {
+/////////////////////////////////////////////////////////////////////////
+// Classe pour la gestion de persistance
+template <class T>
+class FIO_TArray : public  PPersist  {
+public:
+  FIO_TArray();
+  FIO_TArray(string const & filename); 
+  FIO_TArray(const TArray<T> & obj);
+  FIO_TArray(TArray<T> * obj);
+  virtual ~FIO_TArray();
+  virtual AnyDataObj* DataObj();
+  virtual void        SetDataObj(AnyDataObj & o);
+  inline operator TArray<T>() { return(*dobj); }
+protected :
+  virtual void ReadSelf(PInPersist&);           
+  virtual void WriteSelf(POutPersist&) const;  
+  TArray<T> * dobj;
+  bool ownobj;
+};
+
+template <class T>
+inline POutPersist& operator << (POutPersist& os, TArray<T> & obj)
+{ FIO_TArray<T> fio(&obj);  fio.Write(os);  return(os); }
+template <class T>
+inline PInPersist& operator >> (PInPersist& is, TArray<T> & obj)
+{ FIO_TArray<T> fio(&obj);  fio.Read(is);  return(is); }
+
+} // Fin du namespace
+
+#endif
Index: /trunk/SophyaLib/TArray/sopemtx.cc
===================================================================
--- /trunk/SophyaLib/TArray/sopemtx.cc	(revision 772)
+++ /trunk/SophyaLib/TArray/sopemtx.cc	(revision 772)
@@ -0,0 +1,367 @@
+//                         C.Magneville, R.Ansari        03/2000
+
+#include "machdefs.h"
+#include <stdio.h>
+#include <iostream.h>
+#include <complex>
+#include <math.h>
+#include "sopemtx.h"
+#include "smathconst.h"
+
+////////////////////////////////////////////////////////////////
+// -------------------------------------------------------------
+//   La classe de gestion des lignes et colonnes d'une matrice
+// -------------------------------------------------------------
+////////////////////////////////////////////////////////////////
+
+template <class T> TMatrixRC<T>::TMatrixRC()
+: matrix(NULL), data(NULL), index(0), step(0)
+{}
+
+template <class T> TMatrixRC<T>::TMatrixRC(TMatrix<T>& m,TRCKind rckind,uint_4 ind)
+: matrix(&m), data(Org(m,rckind,ind)),
+  index(ind), step(Step(m,rckind)), kind(rckind)
+{
+if (kind == TmatrixDiag && m.NCols() != m.NRows())
+  throw(SzMismatchError("TMatrixRC::TMatrixRC(...,TmatrixDiag,...) size mismatch\n"));
+}
+
+////////////////////////////////////////////////////////////////
+// Acces aux rangees et colonnes de matrices
+
+template <class T> 
+TMatrixRC<T> TMatrixRC<T>::Row(TMatrix<T> & m, uint_4 r) 
+{
+TMatrixRC<T> rc(m, TmatrixRow, r);
+return rc;
+}
+
+template <class T> 
+TMatrixRC<T> TMatrixRC<T>::Col(TMatrix<T> & m, uint_4 c) 
+{
+TMatrixRC<T> rc(m, TmatrixCol, c);
+return rc;
+}
+
+template <class T> 
+TMatrixRC<T> TMatrixRC<T>::Diag(TMatrix<T> & m) 
+{
+TMatrixRC<T> rc(m, TmatrixDiag);
+return rc;
+}
+
+
+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)
+{
+if(!matrix) return -1;
+if(c<0 || c>(int_4)matrix->NCols()) return -1;
+kind = TmatrixCol;
+index = c;
+step = Step(*matrix, TmatrixCol);
+data = Org(*matrix, TmatrixCol, c);
+return c;
+}
+
+template <class T> int_4 TMatrixRC<T>::SetRow(int_4 r)
+{
+if(!matrix) return -1;
+if(r<0 && r>(int_4)matrix->NRows()) return -1;
+kind = TmatrixRow;
+index = r;
+step = Step(*matrix, TmatrixRow);
+data = Org(*matrix, TmatrixRow, r);
+return r;
+}
+
+template <class T> int_4 TMatrixRC<T>::SetDiag()
+{
+if (!matrix) return -1;
+if (matrix->NCols() != matrix->NRows())
+  throw(SzMismatchError("TMatrixRC::SetDiag size mismatch\n"));
+kind = TmatrixDiag;
+index = 0;
+step = Step(*matrix, TmatrixDiag);
+data = Org(*matrix, TmatrixDiag);
+return 0;
+}
+
+
+template <class T> TMatrixRC<T>& TMatrixRC<T>::operator = (const TMatrixRC<T>& rc)
+{
+matrix = rc.matrix;
+data   = rc.data;
+index  = rc.index;
+step   = rc.step;
+kind   = rc.kind;
+return *this;
+}
+
+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;
+}
+
+
+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> 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>
+TMatrixRC<T>& TMatrixRC<T>::LinComb(T a, T b, const TMatrixRC<T>& rc, uint_4 first)
+{
+if ( NElts() != rc.NElts() )
+  throw(SzMismatchError("TMatrixRC::LinComb size mismatch\n"));
+if ( kind != rc.kind )
+  throw(SzMismatchError("TMatrixRC::LinComb type mismatch\n"));
+for(uint_4 i=first; i<NElts(); i++) (*this)(i) = (*this)(i)*a + rc(i)*b;
+return *this;
+}
+
+template <class T>
+TMatrixRC<T>& TMatrixRC<T>::LinComb(T b, const TMatrixRC<T>& rc, uint_4 first)
+{
+if ( NElts() != rc.NElts() )
+  throw(SzMismatchError("TMatrixRC::LinComb size mismatch\n"));
+if ( kind != rc.kind )
+  throw(SzMismatchError("TMatrixRC::LinComb type mismatch\n"));
+for(uint_4 i=first; i<NElts(); i++) (*this)(i) += rc(i)*b;
+return *this;
+}
+
+template <class T> uint_4 TMatrixRC<T>::IMaxAbs(uint_4 first)
+{
+if (first>NElts())
+  throw(SzMismatchError("TMatrixRC::IMaxAbs size mismatch\n"));
+uint_4 imax = first;
+double vmax = Abs_Value((*this)(first));
+for(uint_4 i=first+1; i<NElts(); i++) {
+  double v = Abs_Value((*this)(i));
+  if(v > vmax) {vmax = v; imax = i;}
+}
+return imax;
+}
+
+template <class T>
+void TMatrixRC<T>::Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2)
+{
+if(rc1.NElts() != rc2.NElts())
+  throw(SzMismatchError("TMatrixRC::Swap size mismatch\n"));
+if(rc1.kind != rc2.kind)
+  throw(SzMismatchError("TMatrixRC::Swap type mismatch\n"));
+if(rc1.data == rc2.data) return;
+for(uint_4 i=0; i<rc1.NElts(); i++)
+  {T tmp = rc1(i); rc1(i) = rc2(i); rc2(i) = tmp;}
+}
+
+
+
+
+////////////////////////////////////////////////////////////////
+//**** Pour inversion
+#ifndef M_LN2
+#define M_LN2 0.69314718055994530942
+#endif
+#ifndef LN_MINDOUBLE
+#define LN_MINDOUBLE  (M_LN2 * (DMINEXP - 1))
+#endif
+#ifndef LN_MAXDOUBLE
+#define LN_MAXDOUBLE  (M_LN2 * DMAXEXP)
+#endif
+
+template <class T>
+T SimpleMatrixOperation<T>::GausPiv(TMatrix<T>& a, TMatrix<T>& b)
+// Pivot de Gauss
+// * Attention: egcs impose que cette fonction soit mise dans le .cc
+//              avant ::Inverse() (car Inverse() l'utilise)
+// {TMatrix A(a); TMatrix B(b); return (T) TMatrix::GausPiv(A,B);}
+{
+uint_4 n = a.NRows();
+if(n!=b.NRows())
+  throw(SzMismatchError("TMatrix::GausPiv size mismatch\n"));
+// On fait une normalisation un peu brutale...
+double vmin=MAXDOUBLE;
+double vmax=0;
+for(uint_4 iii=0; iii<a.NRows(); iii++)
+  for(uint_4 jjj=0; jjj<a.NCols(); jjj++) {
+    double v = TMatrixRC<T>::Abs_Value(a(iii,jjj));
+    if(v>vmax) vmax = v;
+    if(v<vmin && v>0) vmin = v;
+}
+double nrm = sqrt(vmin*vmax);
+if(nrm > 1.e5 || nrm < 1.e-5) {
+  a /= nrm;
+  b /= nrm;
+  //cout << "normalisation matrice " << nrm << endl;
+} else nrm=1;
+
+double det = 1.0;
+if(nrm != 1) {
+  double ld = a.NRows() * log(nrm);
+  if (ld <= LN_MINDOUBLE || ld >= LN_MAXDOUBLE) {
+   // cerr << "TMatrix warning, overflow for det" << endl;
+  } else {
+    det = exp(ld);
+  }
+}
+
+TMatrixRC<T> pivRowa(a,TmatrixRow);
+TMatrixRC<T> pivRowb(b,TmatrixRow);
+
+for(uint_4 k=0; k<n-1; k++) {
+  uint_4 iPiv = TMatrixRC<T>::Col(a, k).IMaxAbs(k);
+  if(iPiv != k) {
+    TMatrixRC<T> aIPiv(TMatrixRC<T>::Row(a,iPiv));
+    TMatrixRC<T> aK(TMatrixRC<T>::Row(a, k));
+    TMatrixRC<T>::Swap(aIPiv,aK);
+    TMatrixRC<T> bIPiv(TMatrixRC<T>::Row(b, iPiv));
+    TMatrixRC<T> bK(TMatrixRC<T>::Row(b, k));
+    TMatrixRC<T>::Swap(bIPiv,bK);
+  }
+  double pivot = a(k,k);
+  if (fabs(pivot) < 1.e-50) return 0.0;
+  //det *= pivot;
+  pivRowa.SetRow(k); // to avoid constructors
+  pivRowb.SetRow(k);
+  for (uint_4 i=k+1; i<n; i++) {
+    double r = -a(i,k)/pivot;
+    TMatrixRC<T>::Row(a, i).LinComb(r, pivRowa); // + rapide que -= r * pivRowa
+    TMatrixRC<T>::Row(b, i).LinComb(r, pivRowb);
+  }
+}
+det *= a(n-1, n-1);
+
+// on remonte
+for(uint_4 kk=n-1; kk>0; kk--) {
+  double pivot = a(kk,kk);
+  if (fabs(pivot) <= 1.e-50) return 0.0;
+  pivRowa.SetRow(kk); // to avoid constructors
+  pivRowb.SetRow(kk);
+  for(uint_4 jj=0; jj<kk; jj++) {
+    double r = -a(jj,kk)/pivot;
+    TMatrixRC<T>::Row(a, jj).LinComb(r, pivRowa);
+    TMatrixRC<T>::Row(b, jj).LinComb(r, pivRowb);
+  }
+}
+
+for(uint_4 l=0; l<n; l++) {
+  if (fabs(a(l,l)) <= 1.e-50) return 0.0;
+  TMatrixRC<T>::Row(b, l) /= a(l,l);
+}
+
+return det;
+}
+
+template <class T>
+TMatrix<T> SimpleMatrixOperation<T>::Inverse(TMatrix<T> const & A) 
+// Inversion
+{
+TMatrix<T> a(A);
+TMatrix<T> b(a.NCols(),a.NRows());  b = 1.;
+if(fabs(GausPiv(a,b)) < 1.e-50)
+  throw(MathExc("TMatrix Inverse() Singular OMatrix"));
+return b;
+}
+
+
+
+///////////////////////////////////////////////////////////////
+#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<r_8>
+#endif
+
+#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<r_4>;
+template class SimpleMatrixOperation<r_8>;
+#endif
Index: /trunk/SophyaLib/TArray/sopemtx.h
===================================================================
--- /trunk/SophyaLib/TArray/sopemtx.h	(revision 772)
+++ /trunk/SophyaLib/TArray/sopemtx.h	(revision 772)
@@ -0,0 +1,154 @@
+// This may look like C code, but it is really -*- C++ -*-
+#ifndef SOpeMatrix_SEEN
+#define SOpeMatrix_SEEN
+
+#include "machdefs.h"
+#include "tmatrix.h"
+#include "tvector.h"
+
+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;
+
+////////////////////////////////////////////////////////////////
+template <class T>
+class SimpleMatrixOperation {
+public:
+  // Pivot de Gauss : diagonalise la matrice A, en effectuant les memes
+  // operations sur la matrice B
+  static TMatrix<T> Inverse(TMatrix<T> const & A);
+  static T GausPiv(TMatrix<T>& A, TMatrix<T>& B);
+};
+
+// Resolution du systeme A*C = B
+inline r_8 LinSolveInPlace(TMatrix<r_8>& a, TVector<r_8>& b)
+{
+if(a.NCols() != b.NRows() || a.NCols() != a.NRows())
+  throw(SzMismatchError("LinSolveInPlace(TMatrix<r_8>,TVector<r_8>) size mismatch"));
+return SimpleMatrixOperation<r_8>::GausPiv(a,b);
+}
+
+// Resolution du systeme A*C = B, avec C retourne dans B
+inline r_8 LinSolve(const TMatrix<r_8>& a, const TVector<r_8>& b, TVector<r_8>& c)
+{ 
+if(a.NCols() != b.NRows() || a.NCols() != a.NRows())
+  throw(SzMismatchError("LinSolve(TMatrix<r_8>,TVector<r_8>) size mismatch"));
+c = b;
+TMatrix<r_8> a1(a);
+return SimpleMatrixOperation<r_8>::GausPiv(a1,c);
+}
+
+
+} // Fin du namespace
+
+#endif
Index: /trunk/SophyaLib/TArray/tarray.cc
===================================================================
--- /trunk/SophyaLib/TArray/tarray.cc	(revision 772)
+++ /trunk/SophyaLib/TArray/tarray.cc	(revision 772)
@@ -0,0 +1,586 @@
+//      template array class for numerical types
+//                     R. Ansari, C.Magneville   03/2000
+
+#include "machdefs.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include "pexceptions.h"
+#include "tarray.h"
+
+
+// Classe utilitaires 
+Sequence::Sequence(double start, double step)
+{
+  start_ = start;
+  step_ = step;
+}
+
+// Variables statiques globales
+template <class T>
+char * TArray<T>::ck_op_msg_[6] = {"???", "Size(int )", "IsPacked(int )", 
+				 "Stride(int )", "ElemCheckBound()", "operator()" };
+template <class T>
+uint_4 TArray<T>::max_nprt_ = 50;
+
+// Methodes statiques globales
+template <class T>
+void TArray<T>::SetMaxPrint(uint_4 nprt)
+{
+  max_nprt_ = nprt;
+}
+
+
+template <class T>
+uint_8 TArray<T>::ComputeTotalSize(uint_4 ndim, uint_4 * siz, uint_4 step, uint_8 offset)
+{
+  uint_8 rs = step;
+  for(int k=0; k<ndim; k++) rs *= siz[k];
+  return(rs+offset);
+}
+
+// -------------------------------------------------------
+//                Methodes de la classe
+// -------------------------------------------------------
+
+// Les constructeurs 
+template <class T>
+TArray<T>::TArray()
+  : mNDBlock() , mInfo(NULL)
+// Default constructor
+{
+  ndim_ = 0;
+  for(int k=0; k<TARRAY_MAXNDIMS; k++) step_[k] = size_[k] = 0;
+  totsize_ = 0;
+  minstep_ = 0;
+  offset_ = 0;
+  // Default for matrices : Lines = Y , Columns = X
+  macoli_ = 0;
+  marowi_ = 1;
+}
+
+template <class T>
+TArray<T>::TArray(uint_4 ndim, uint_4 * siz, uint_4 step)
+  : mNDBlock(ComputeTotalSize(ndim, siz, step, 1)) , mInfo(NULL)
+{
+  string exmsg = "TArray<T>::TArray(uint_4, uint_4 *, uint_4)";
+  if (!UpdateSizes(ndim, siz, step, 0, exmsg))  throw( ParmError(exmsg) );
+}
+
+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)) , mInfo(NULL)
+{
+  uint_4 size[TARRAY_MAXNDIMS];
+  size[0] = nx;   size[1] = ny;   size[2] = nz;
+  int ndim = 1;
+  if ((size[1] > 0) && (size[2] > 0)) ndim = 3;
+  else if (size[1] > 0)  ndim = 2;
+  else ndim = 1;
+  string exmsg = "TArray<T>::TArray(uint_4, uint_4, uint_4)";
+  if (!UpdateSizes(ndim, size, 1, 0, exmsg))  throw( ParmError(exmsg) );
+}
+
+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) , mInfo(NULL)
+{
+  string exmsg = "TArray<T>::TArray(uint_4, uint_4 *,  NDataBlock<T> & ... )";
+  if (!UpdateSizes(ndim, siz, step, offset, exmsg))  throw( ParmError(exmsg) );
+  
+}
+
+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) , mInfo(NULL)
+{
+  string exmsg = "TArray<T>::TArray(uint_4, uint_4 *, T* ... )";
+  if (!UpdateSizes(ndim, siz, step, offset, exmsg))  throw( ParmError(exmsg) );
+}
+
+template <class T>
+TArray<T>::TArray(const TArray<T>& a)
+  : mNDBlock(a.mNDBlock) , mInfo(NULL)
+{
+  string exmsg = "TArray<T>::TArray(const TArray<T>&)";
+  if (!UpdateSizes(a, exmsg))  throw( ParmError(exmsg) );
+  if (a.mInfo) mInfo = new DVList(*(a.mInfo));
+}
+
+template <class T>
+TArray<T>::TArray(const TArray<T>& a, bool share)
+  : mNDBlock(a.mNDBlock, share) , mInfo(NULL)
+{
+  string exmsg = "TArray<T>::TArray(const TArray<T>&, bool)";
+  if (!UpdateSizes(a, exmsg))  throw( ParmError(exmsg) );
+  if (a.mInfo) mInfo = new DVList(*(a.mInfo));
+}
+
+// Destructeur 
+template <class T>
+TArray<T>::~TArray()
+{
+}
+
+template <class T>
+TArray<T> TArray<T>::SubArray(uint_4 ndim, uint_4 * siz, uint_4 * pos)
+{
+  if ( (ndim > ndim_) || (ndim < 1) ) throw(SzMismatchError("TArray<T>::SubArray( ... ) NDim Error") );
+  int k;
+  for(k=0; k<ndim; k++) 
+    if ( (siz[k]+pos[k]) > size_[k] )  
+      throw(SzMismatchError("TArray<T>::SubArray( ... ) Size/Pos Error (1)") );
+  for(k=ndim; k<ndim_; k++)   
+    if ( (pos[k]) >= size_[k] )  
+      throw(SzMismatchError("TArray<T>::SubArray( ... ) Size/Pos Error (2)") );
+  TArray<T> ra(*this);
+  ra.ndim_ = ndim;
+  for(k=ndim; k<TARRAY_MAXNDIMS; k++) {
+    ra.size_[k] = 1;
+    ra.step_[k] = 0;
+  }
+  ra.offset_ = offset_;
+  for(k=0; k<ndim_; k++) ra.offset_ += pos[k]*step_[k]; 
+  ra.SetTemp(true);
+  return(ra);   
+}
+
+template <class T>
+TArray<T>& TArray<T>::operator = (const TArray<T>& a)
+{
+  if (this != &a) { 
+    CloneOrShare(a);
+    if (mInfo) delete mInfo;
+    if (a.mInfo) mInfo = new DVList(*(a.mInfo));
+  }
+  return(*this);
+}
+
+template <class T>
+void TArray<T>::Clone(const TArray<T>& a)
+{
+  string exmsg = "TArray<T>::Clone()";
+  if (!UpdateSizes(a, exmsg))  throw( ParmError(exmsg) );
+  mNDBlock.Clone(a.mNDBlock);
+  if (mInfo) delete mInfo;
+  if (a.mInfo) mInfo = new DVList(*(a.mInfo));
+}
+
+template <class T>
+void TArray<T>::ReSize(uint_4 ndim, uint_4 * siz, uint_4 step)
+{
+  string exmsg = "TArray<T>::ReSize()";
+  if (!UpdateSizes(ndim, siz, step, 0, exmsg))  throw( ParmError(exmsg) );
+  mNDBlock.ReSize(totsize_);    
+}
+
+template <class T>
+void TArray<T>::Realloc(uint_4 ndim, uint_4 * siz, uint_4 step, bool force)
+{
+  string exmsg = "TArray<T>::Realloc()";
+  if (!UpdateSizes(ndim, siz, step, 0, exmsg))  throw( ParmError(exmsg) );
+  mNDBlock.ReSize(totsize_);    
+}
+
+template <class T>
+bool TArray<T>::CompareSizes(const TArray<T>& a)
+{
+  if (ndim_ != a.ndim_)  return(false);
+  for(int k=0; k<ndim_; k++) 
+    if (size_[k] != a.size_[k])  return(false);
+  if ( (macoli_ != a.macoli_) || (marowi_ != a.marowi_) )  return(false);
+  return(true);
+}
+
+//   ...... Operation de calcul sur les tableaux ......
+//              ------- Attention --------
+//  Boucles normales prenant en compte les steps ....
+//  Possibilite de // , vectorisation 
+template <class T>
+TArray<T>& TArray<T>::operator = (Sequence & seq)
+{
+  T * pe = Data();
+  uint_8 step = Step();
+  for(uint_8 k=0; k<totsize_; k++ )  pe[k*step] = seq(k);
+  return(*this);
+}
+
+//  >>>> Operations avec 2nd membre de type scalaire 
+
+template <class T>
+TArray<T>& TArray<T>::operator = (T x)
+{
+  T * pe = Data();
+  uint_8 step = Step();
+  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] = x;
+  return(*this);
+}
+
+template <class T>
+TArray<T>& TArray<T>::operator += (T x)
+{
+  T * pe = Data();
+  uint_8 step = Step();
+  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] += x;
+  return(*this);
+}
+
+template <class T>
+TArray<T>& TArray<T>::operator -= (T x)
+{
+  T * pe = Data();
+  uint_8 step = Step();
+  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] -= x;
+  return(*this);
+}
+
+template <class T>
+TArray<T>& TArray<T>::operator *= (T x)
+{
+  T * pe = Data();
+  uint_8 step = Step();
+  for(uint_8 k=0; k<totsize_; k += step)  pe[k] *= x;
+  return(*this);
+}
+
+template <class T>
+TArray<T>& TArray<T>::operator /= (T x)
+{
+  T * pe = Data();
+  uint_8 step = Step();
+  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] /= x;
+  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>&)")) ;
+  T * pe = Data();
+  const T * pea = a.Data();
+  uint_8 step = Step();
+  uint_8 stepa = a.Step();
+  uint_8 k, ka;
+  k = ka = 0;
+  for(k=0; k<totsize_; k += step)  {
+    pe[k] += pea[ka];
+    ka += stepa;
+  }
+  return(*this);
+}
+
+template <class T>
+TArray<T>& TArray<T>::operator -= (const TArray<T>& a)
+{
+  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::operator -= (const TArray<T>&)")) ;
+  T * pe = Data();
+  const T * pea = a.Data();
+  uint_8 step = Step();
+  uint_8 stepa = a.Step();
+  uint_8 k, ka;
+  k = ka = 0;
+  for(k=0; k<totsize_; k += step)  {
+    pe[k] -= pea[ka];
+    ka += stepa;
+  }
+  return(*this);
+}
+
+
+template <class T>
+TArray<T>& TArray<T>::MultElt(const TArray<T>& a)
+{
+  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::MultElt(const TArray<T>&)")) ;
+  T * pe = Data();
+  const T * pea = a.Data();
+  uint_8 step = Step();
+  uint_8 stepa = a.Step();
+  uint_8 k, ka;
+  k = ka = 0;
+  for(k=0; k<totsize_; k += step)  {
+    pe[k] *= pea[ka];
+    ka += stepa;
+  }
+  return(*this);
+}
+
+template <class T>
+TArray<T>& TArray<T>::DivElt(const TArray<T>& a)
+{
+  if (!CompareSizes(a)) throw(SzMismatchError("TArray<T>::DivElt(const TArray<T>&)")) ;
+  T * pe = Data();
+  const T * pea = a.Data();
+  uint_8 step = Step();
+  uint_8 stepa = a.Step();
+  uint_8 k, ka;
+  k = ka = 0;
+  for(k=0; k<totsize_; k += step)  {
+    pe[k] /= pea[ka];
+    ka += stepa;
+  }
+  return(*this);
+}
+
+
+// ----------------------------------------------------
+//          Application d'une fonction
+// ----------------------------------------------------
+template <class T>
+TArray<T>& TArray<T>::ApplyFunction(Arr_DoubleFunctionOfX f)
+{
+  T * pe = Data();
+  uint_8 step = Step();
+  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] = (T)(f((double)pe[k]));
+  return(*this);
+}
+
+template <class T>
+TArray<T>& TArray<T>::ApplyFunction(Arr_FloatFunctionOfX f)
+{
+  T * pe = Data();
+  uint_8 step = Step();
+  for(uint_8 k=0; k<totsize_*step; k += step)  pe[k] = (T)(f((float)pe[k]));
+  return(*this);
+}
+
+TArray< complex<r_4> >& TArray< complex<r_4> >::ApplyFunction(Arr_DoubleFunctionOfX f)
+{
+  throw( NotAvailableOperation("TArray< complex<r_4> >::ApplyFunction() - Unsupported operation ") );
+}
+
+TArray< complex<r_4> >& TArray< complex<r_4> >::ApplyFunction(Arr_FloatFunctionOfX f)
+{
+  throw(NotAvailableOperation( "TArray< complex<r_4> >::ApplyFunction() - Unsupported operation ") );
+}
+
+TArray< complex<r_8> >& TArray< complex<r_8> >::ApplyFunction(Arr_DoubleFunctionOfX f)
+{
+  throw(NotAvailableOperation("TArray< complex<r_8> >::ApplyFunction() - Unsupported operation ") );
+}
+
+TArray< complex<r_8> >& TArray< complex<r_8> >::ApplyFunction(Arr_FloatFunctionOfX f)
+{
+  throw(NotAvailableOperation("TArray< complex<r_8> >::ApplyFunction() - Unsupported operation ") );
+}
+
+// ----------------------------------------------------
+//       Impression, etc ...
+// ----------------------------------------------------
+
+template <class T>
+void TArray<T>::Show(ostream& os, bool si) const
+{
+  os << " TArray< " << typeid(T).name() << " NDim= " << ndim_
+     << "(" << typeid(*this).name() << " )" << endl;
+  os << "TotSize=" << totsize_ << " Size(X*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 << "\n " << endl;
+  if (si && (mInfo != NULL)) os << (*mInfo) << endl;
+ 
+}
+
+template <class T>
+void TArray<T>::Print(ostream& os, int_4 maxprt, bool si) const
+{
+  if (maxprt < 0)  maxprt = max_nprt_;
+  uint_4 npr = 0;
+  Show(os, si);
+  int k0,k1,k2,k3,k4;
+  for(k4=0; k4<size_[4]; k4++) {
+    if (size_[4] > 1) cout << "\n ----- Dimension 5 (U) K4= " << k4;
+    for(k3=0; k3<size_[3]; k3++) {
+      if (size_[3] > 1) cout << "\n ----- Dimension 4 (T) K3= " << k3;
+      for(k2=0; k2<size_[2]; k2++) {
+	if (size_[2] > 1) cout << "\n ----- Dimension 3 (Z) K2= " << k2;
+	for(k1=0; k1<size_[1]; k1++) {
+	  if ( (size_[1] > 1) && (size_[0] > 10) ) cout << "----- Dimension 2 (Y) K1= " << k1;
+	  for(k0=0; k0<size_[0]; k0++) {
+            os << Elem(k0, k1, k2, k3, k4) << ", ";    npr++;
+	    if (npr >= maxprt) {
+	      if (npr < totsize_)  os << "\n     .... " << endl; return;
+	    }
+	  }
+          os << endl;
+	}
+      }
+    }
+  }
+  os << "\n" << endl;
+}
+
+template <class T>
+DVList& TArray<T>::Info()
+{
+if (mInfo == NULL)  mInfo = new DVList;
+return(*mInfo);
+}
+
+
+template <class T>
+bool TArray<T>::UpdateSizes(uint_4 ndim, const uint_4 * siz, uint_4 step, uint_8 offset, string & exmsg)
+{
+  if (ndim >= TARRAY_MAXNDIMS) {
+    exmsg += " NDim Error";  return false;
+  }
+  if (step < 1) {
+    exmsg += " Step(=0) Error";  return false;
+  }
+
+  minstep_ = 1;
+
+  // Flagging bad updates ...
+  ndim_ = 0;
+
+  totsize_ = 1;
+  int k;
+  for(k=0; k<TARRAY_MAXNDIMS; k++) {
+    size_[k] = 1;
+    step_[k] = 0;
+  }
+  for(k=0; k<ndim; k++) {
+    size_[k] = siz[k] ;
+    step_[k] = totsize_*step;
+    totsize_ *= size_[k];
+  }
+  if (totsize_ < 1) {
+    exmsg += " Size Error";  return false;
+  }
+  offset_ = offset;
+  // Default for matrices : Lines = Y , Columns = X
+  macoli_ = 0;
+  marowi_ = 1;
+  // Update OK 
+  ndim_ = ndim;
+  return true;
+}
+
+template <class T>
+bool TArray<T>::UpdateSizes(uint_4 ndim, const uint_4 * siz, const uint_4 * step, uint_8 offset, string & exmsg)
+{
+  if (ndim >= TARRAY_MAXNDIMS) {
+    exmsg += " NDim Error";  return false;
+  }
+
+  // Flagging bad updates ...
+  ndim_ = 0;
+
+  totsize_ = 1;
+  int k;
+  for(k=0; k<TARRAY_MAXNDIMS; k++) {
+    size_[k] = 1;
+    step_[k] = 0;
+  }
+  minstep_ = step[0];
+  for(k=0; k<ndim; k++) {
+    size_[k] = siz[k] ;
+    step_[k] = step[k];
+    totsize_ *= size_[k];
+    if (step_[k] < minstep_)  minstep_ = step_[k];
+  }
+  if (minstep_ < 1) {
+    exmsg += " Step(=0) Error";  return false;
+  }
+  if (totsize_ < 1) {
+    exmsg += " Size Error";  return false;
+  }
+  offset_ = offset;
+  // Default for matrices : Lines = Y , Columns = X
+  macoli_ = 0;
+  marowi_ = 1;
+  // Update OK 
+  ndim_ = ndim;
+  return true;
+}
+template <class T>
+bool TArray<T>::UpdateSizes(const TArray<T>& a, string & exmsg)
+{
+  if (a.ndim_ >= TARRAY_MAXNDIMS) {
+    exmsg += " NDim Error";  return false;
+  }
+
+  // Flagging bad updates ...
+  ndim_ = 0;
+
+  totsize_ = 1;
+  int k;
+  for(k=0; k<TARRAY_MAXNDIMS; k++) {
+    size_[k] = 1;
+    step_[k] = 0;
+  }
+  minstep_ = a.step_[0];
+  for(k=0; k<a.ndim_; k++) {
+    size_[k] = a.size_[k] ;
+    step_[k] = a.step_[k];
+    totsize_ *= size_[k];
+    if (step_[k] < minstep_)  minstep_ = step_[k];
+  }
+  if (minstep_ < 1) {
+    exmsg += " Step(=0) Error";  return false;
+  }
+  if (totsize_ < 1) {
+    exmsg += " Size Error";  return false;
+  }
+
+  offset_ = a.offset_;
+  macoli_ = a.macoli_;
+  marowi_ = a.marowi_;
+  // Update OK 
+  ndim_ = a.ndim_;
+  return true;
+}
+
+template <class T>
+void TArray<T>::CloneOrShare(const TArray<T>& a)
+{
+  string exmsg = "TArray<T>::CloneOrShare()";
+  if (!UpdateSizes(a.ndim_, a.size_, a.step_, a.offset_, exmsg))  throw( ParmError(exmsg) );
+  mNDBlock.CloneOrShare(a.mNDBlock);
+}
+
+template <class T>
+void TArray<T>::Share(const TArray<T>& a)
+{
+  string exmsg = "TArray<T>::Share()";
+  if (!UpdateSizes(a.ndim_, a.size_, a.step_, a.offset_, exmsg))  throw( ParmError(exmsg) );
+  mNDBlock.Share(a.mNDBlock);
+}
+
+
+///////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template TArray<uint_1>
+#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>
+#pragma define_template TArray< complex<r_4> > 
+#pragma define_template TArray< complex<r_8> > 
+#endif
+
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class TArray<uint_1>;
+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>;
+template class TArray< complex<r_4> >;
+template class TArray< complex<r_8> >;
+#endif
+
+
Index: /trunk/SophyaLib/TArray/tarray.h
===================================================================
--- /trunk/SophyaLib/TArray/tarray.h	(revision 772)
+++ /trunk/SophyaLib/TArray/tarray.h	(revision 772)
@@ -0,0 +1,358 @@
+// This may look like C code, but it is really -*- C++ -*-
+//      template array class for numerical types
+//                     R. Ansari, C.Magneville   03/2000
+
+#ifndef TArray_SEEN
+#define TArray_SEEN
+
+#include "machdefs.h"
+#include <math.h>
+#include <iostream.h>
+#include "ndatablock.h"
+#include <complex>
+#include "anydataobj.h"
+#include "ndatablock.h"
+#include "dvlist.h"
+
+
+//  Maximum number of dimensions for array
+#define TARRAY_MAXNDIMS  5
+
+namespace SOPHYA {
+
+/* Quelques utilitaires */
+typedef double (* Arr_DoubleFunctionOfX) (double x);
+typedef float  (* Arr_FloatFunctionOfX)  (float x);
+
+class Sequence { 
+public:
+  Sequence (double start=0., double step=1.);
+  inline double & Start() { return start_; }
+  inline double & Step() { return step_; }
+  inline double operator () (uint_4 k) { return(start_+(double)k*step_); }
+protected:
+  double start_, step_;
+};
+
+// Forward declaration
+template <class T> class FIO_TArray;
+
+//   ------------ classe template Array ----------- 
+template <class T>
+class TArray : public AnyDataObj {
+public:
+  // 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(const TArray<T>& a);
+  TArray(const TArray<T>& a, bool share);
+
+  virtual ~TArray();
+
+  // A = B : partage les donnees si "a" est temporaire, clone sinon.
+  virtual TArray<T>& operator = (const TArray<T>& a);
+
+  // Extraction de sous-tableau
+  virtual TArray<T> SubArray(uint_4 ndim, uint_4 * siz, uint_4 * pos);
+
+  // 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);
+  // Returns true if ndim and sizes are equal
+  virtual bool CompareSizes(const TArray<T>& a);
+ 
+  // Array dimensions
+  inline uint_4 NbDimensions() const { return( ndim_ ); }
+
+  inline uint_8 Size() const { return(totsize_); }
+  inline uint_4 SizeX() const { return(size_[0]); }
+  inline uint_4 SizeY() const { return(size_[1]); }
+  inline uint_4 SizeZ() const { return(size_[2]); }
+  inline uint_4 Size(int ka) const { return(size_[CheckDI(ka,1)]); }
+
+  // memory organisation - packing information
+  inline bool   IsPacked() const { return(minstep_ == 1); }
+  inline bool   IsPackedX() const { return(step_[0] == 1); }
+  inline bool   IsPackedY() const { return(step_[1] == 1); }
+  inline bool   IsPackedZ() const { return(step_[2] == 1); }
+  inline bool   IsPacked(int ka) const { return(step_[CheckDI(ka,2)] == 1); }
+
+  inline uint_4 Step() const  { return(minstep_); }
+  inline uint_4 StepX() const { return(step_[0]); }
+  inline uint_4 StepZ() const { return(step_[1]); }
+  inline uint_4 StepY() const { return(step_[2]); }
+  inline uint_4 Step(int ka) const { return(step_[CheckDI(ka,3)]); }
+
+  inline uint_8 Offset() const { return(offset_); }
+
+
+  // Temporaire?
+  inline bool   IsTemp(void) const {return mNDBlock.IsTemp();}
+  inline void   SetTemp(bool temp=false) const {mNDBlock.SetTemp(temp);}
+
+  // Acces to data
+  inline T const& operator()(uint_4 ix, uint_4 iy, uint_4 iz) const ;
+  inline T&       operator()(uint_4 ix, uint_4 iy, uint_4 iz);
+  inline T const& operator()(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu=0) const ;
+  inline T&       operator()(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu=0);
+  inline T const& operator[](uint_4 ip) const ;
+  inline T&       operator[](uint_4 ip);
+
+  inline T const& Elem(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it=0, uint_4 iu=0) const ;
+  inline T&       Elem(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it=0, uint_4 iu=0);
+  inline T const& ElemCheckBound(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it=0, uint_4 iu=0) const ;
+  inline T&       ElemCheckBound(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it=0, uint_4 iu=0);
+
+  inline        T* Data()       {return mNDBlock.Begin()+offset_;}
+  inline const  T* Data() const {return mNDBlock.Begin()+offset_;}
+  inline        NDataBlock<T>& DataBlock()       {return mNDBlock;}
+  inline const  NDataBlock<T>& DataBlock() const {return mNDBlock;}
+
+// Operations diverses  = , +=, ...
+// Conversion en type T, if Size() == 1
+  inline operator T();
+// Met les elements a une suite de valeurs
+  virtual TArray<T>&  operator = (Sequence & seq);
+// A = x (tous les elements a x)
+  virtual TArray<T>&  operator = (T 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);
+// 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);
+// Multilication, division element par element les deux tableaux 
+  virtual TArray<T>&  MultElt(const TArray<T>& a);
+  virtual TArray<T>&  DivElt(const TArray<T>& a);
+
+// Application d'une fonction
+  virtual TArray<T>&  ApplyFunction(Arr_DoubleFunctionOfX f);
+  virtual TArray<T>&  ApplyFunction(Arr_FloatFunctionOfX f);
+
+// 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          Print(ostream& os, int_4 maxprt=-1, bool si=false) const ;
+
+//  Objet DVList info 
+  DVList&       Info();
+
+//  Pour la gestion de persistance
+  friend class  FIO_TArray<T>;
+
+protected:
+  // Verifie la compatibilite de l'index de dimension
+  inline int CheckDI(int ka, int msg) const ;
+  // Verifie la compatibilite des bornes d'index
+  inline void CheckBound(int ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu, int msg) const ;
+  // Changing Sizes/NDim ... return true if OK
+  bool UpdateSizes(uint_4 ndim, const uint_4 * siz, uint_4 step, uint_8 offset, string & exmsg);
+  bool UpdateSizes(uint_4 ndim, const uint_4 * siz, const uint_4 * step, uint_8 offset, string & exmsg);
+  bool UpdateSizes(const TArray<T>& a, string & exmsg);
+  static uint_8 ComputeTotalSize(uint_4 ndim, uint_4 * siz, uint_4 step, uint_8 offset) ;
+  // partage les donnees si "a" temporaire, clone sinon.
+  void CloneOrShare(const TArray<T>& a);
+  // Share: partage les donnees de "a"
+  void Share(const TArray<T>& a);
+
+  uint_4 ndim_;                   // nb of dimensions
+  uint_4 size_[TARRAY_MAXNDIMS];     // array size in each dimension
+  uint_8 totsize_;                   // Total number of elements
+  uint_4 step_[TARRAY_MAXNDIMS];     // two consecutive elements distance in a given dimension 
+  uint_4 minstep_;                   // minimal step (in any axes)
+  uint_8 offset_;              // global offset -> position of elem[0] in DataBlock
+  NDataBlock<T> mNDBlock;      // Le bloc des donnees
+  uint_4 marowi_, macoli_;     // For matrices, Row index and column index in dimensions
+
+  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
+};
+
+////////////////////////////////////////////////////////////////
+//   Surcharge d'operateur << 
+template <class T>
+inline ostream& operator << (ostream& os, const TArray<T>& a)
+                            { a.Print(os);    return(os);    }
+
+////////////////////////////////////////////////////////////////
+// Surcharge d'operateurs A (+,-,*,/) (T) x
+
+template <class T> inline TArray<T> operator + (const TArray<T>& a, T b)
+    {TArray<T> result(a); result.SetTemp(true); result += 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;}
+
+template <class T> inline TArray<T> operator - (const TArray<T>& a, T b)
+    {TArray<T> result(a); result.SetTemp(true); result -= 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;}
+
+template <class T> inline TArray<T> operator * (const TArray<T>& a, T b)
+    {TArray<T> result(a); result.SetTemp(true); result *= 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;}
+
+template <class T> inline TArray<T> operator / (const TArray<T>& a, T b)
+    {TArray<T> result(a); result.SetTemp(true); result /= b; return result;}
+
+////////////////////////////////////////////////////////////////
+// Surcharge d'operateurs C = A (+,-) B
+
+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;}
+
+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;}
+
+////////////////////////////////////////////////
+//  Calcul de fonction 
+template <class T>
+inline TArray<T> sin(const TArray<T>& a)
+    {TArray<T> result(a); result.SetTemp(true); result.ApplyFunction(sin); return result;}
+
+template <class T>
+inline TArray<T> cos(const TArray<T>& a)
+    {TArray<T> result(a); result.SetTemp(true); result.ApplyFunction(cos); return result;}
+
+// --------------------------------------------------
+//        Methodes inline de verification
+// --------------------------------------------------
+template <class T>
+inline int TArray<T>::CheckDI(int ka, int msg)  const
+{
+  if ( (ka < 0) || (ka >= ndim_) ) {
+    string txt = "TArray<T>::CheckDimensionIndex/Error ";   txt += ck_op_msg_[msg];
+    throw(ParmError(txt));
+  }
+  return(ka);
+}
+
+template <class T>
+inline void TArray<T>::CheckBound(int ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu, int msg)  const
+{
+  if ( (ix >= size_[0]) || (iy >= size_[1]) || (iz > size_[2]) ||
+       (it >= size_[3]) || (iu >= size_[4]) ) {
+    string txt = "TArray<T>::CheckArrayBound/Error ";   txt += ck_op_msg_[msg];
+    throw(ParmError(txt));
+  }
+  return;
+}
+
+// --------------------------------------------------
+//        inline element acces methods
+// --------------------------------------------------
+template <class T>
+inline T const& TArray<T>::Elem(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu) const
+{
+  return ( *( mNDBlock.Begin()+ offset_+ 
+	      ix*step_[0] + iy*step_[1] + iz*step_[2] +
+	      it*step_[3] + iu*step_[4]) );
+}
+
+template <class T>
+inline T & TArray<T>::Elem(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu)
+{
+  return ( *( mNDBlock.Begin()+ offset_+ 
+	      ix*step_[0] + iy*step_[1] + iz*step_[2] +
+	      it*step_[3] + iu*step_[4]) );
+}
+
+template <class T>
+inline T const& TArray<T>::ElemCheckBound(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu) const
+{
+  CheckBound(ix, iy, iz, it, iu, 4);
+  Elem(ix, iy, iz, it, iu);
+}
+
+template <class T>
+inline T & TArray<T>::ElemCheckBound(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu)
+{
+  CheckBound(ix, iy, iz, it, iu, 4);
+  Elem(ix, iy, iz, it, iu);
+}
+
+template <class T>
+inline T const& TArray<T>::operator()(uint_4 ix, uint_4 iy, uint_4 iz) const
+{
+#ifdef SO_BOUNDCHECKING
+  CheckBound(ix, iy, iz, 0, 0, 4);
+#endif
+  return ( *( mNDBlock.Begin()+ offset_+ 
+	      ix*step_[0] + iy*step_[1] + iz*step_[2]) );
+}
+
+template <class T>
+inline T & TArray<T>::operator()(uint_4 ix, uint_4 iy, uint_4 iz) 
+{
+#ifdef SO_BOUNDCHECKING
+  CheckBound(ix, iy, iz, 0, 0, 4);
+#endif
+  return ( *( mNDBlock.Begin()+ offset_+ 
+	      ix*step_[0] + iy*step_[1] + iz*step_[2]) );
+}
+
+template <class T>
+inline T const& TArray<T>::operator()(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu) const
+{
+#ifdef SO_BOUNDCHECKING
+  CheckBound(ix, iy, iz, it, iu, 4);
+#endif
+  return ( *( mNDBlock.Begin()+ offset_+ 
+	      ix*step_[0] + iy*step_[1] + iz*step_[2] +
+	      it*step_[3] + iu*step_[4]) );
+}
+
+template <class T>
+inline T & TArray<T>::operator()(uint_4 ix, uint_4 iy, uint_4 iz, uint_4 it, uint_4 iu)
+{
+#ifdef SO_BOUNDCHECKING
+  CheckBound(ix, iy, iz, it, iu, 4);
+#endif
+  return ( *( mNDBlock.Begin()+ offset_+ 
+	      ix*step_[0] + iy*step_[1] + iz*step_[2] +
+	      it*step_[3] + iu*step_[4]) );
+}
+
+template <class T>
+inline T const& TArray<T>::operator[](uint_4 ip) const
+{
+#ifdef SO_BOUNDCHECKING
+  if (ip >= totsize_)  throw( ParmError("TArray<T>::operator[] Out-of-bound Error") );
+#endif
+return *(mNDBlock.Begin()+offset_+ip*minstep_);
+}
+
+template <class T>
+inline T & TArray<T>::operator[](uint_4 ip) 
+{
+#ifdef SO_BOUNDCHECKING
+  if (ip >= totsize_)  throw( ParmError("TArray<T>::operator[] Out-of-bound Error") );
+#endif
+return *(mNDBlock.Begin()+offset_+ip*minstep_);
+}
+
+template <class T>
+inline TArray<T>::operator T()
+{ 
+  if (Size() != 1) throw(SzMismatchError("TArray<T>::operator T() Size() != 1")) ;
+  return ( (*this)[0] );
+}
+
+} // Fin du namespace
+
+#endif
Index: /trunk/SophyaLib/TArray/tarrinit.cc
===================================================================
--- /trunk/SophyaLib/TArray/tarrinit.cc	(revision 771)
+++ /trunk/SophyaLib/TArray/tarrinit.cc	(revision 772)
@@ -5,4 +5,5 @@
 #include "tmatrix.h"
 #include "tvector.h"
+#include "fioarr.h"
 
 
@@ -17,4 +18,15 @@
 //   Enregistrement des classes PPersist du modules TArray
   
+
+  PPRegister(FIO_TArray<uint_1>);
+  PPRegister(FIO_TArray<uint_2>);
+  PPRegister(FIO_TArray<int_2>);
+  PPRegister(FIO_TArray<int_4>);
+  PPRegister(FIO_TArray<int_8>);
+  PPRegister(FIO_TArray<uint_8>);
+  PPRegister(FIO_TArray<r_4>);
+  PPRegister(FIO_TArray<r_8>);
+  PPRegister(FIO_TArray< complex<r_4> >);
+  PPRegister(FIO_TArray< complex<r_8> >);
 
   PPRegister(FIO_TMatrix<uint_1>);
Index: /trunk/SophyaLib/TArray/tmatrix.cc
===================================================================
--- /trunk/SophyaLib/TArray/tmatrix.cc	(revision 771)
+++ /trunk/SophyaLib/TArray/tmatrix.cc	(revision 772)
@@ -1,3 +1,3 @@
-// $Id: tmatrix.cc,v 1.1.1.1 2000-03-02 16:48:24 ansari Exp $
+// $Id: tmatrix.cc,v 1.2 2000-03-10 15:13:22 ansari Exp $
 //                         C.Magneville          04/99
 #include "machdefs.h"
@@ -161,128 +161,4 @@
 }
 
-////////////////////////////////////////////////////////////////
-//**** Pour gestion des lignes et des colonnes
-
-template <class T> TMatrixRC<T> TMatrix<T>::Row(uint_4 r) const
-{
-TMatrixRC<T> rc((TMatrix<T>&)*this, TmatrixRow, r);
-return rc;
-}
-
-template <class T> TMatrixRC<T> TMatrix<T>::Col(uint_4 c) const
-{
-TMatrixRC<T> rc((TMatrix<T>&)*this, TmatrixCol, c);
-return rc;
-}
-
-template <class T> TMatrixRC<T> TMatrix<T>::Diag() const
-{
-TMatrixRC<T> rc((TMatrix<T>&)*this, TmatrixDiag);
-return rc;
-}
-
-////////////////////////////////////////////////////////////////
-//**** Pour inversion
-#ifndef M_LN2
-#define M_LN2 0.69314718055994530942
-#endif
-#ifndef LN_MINDOUBLE
-#define LN_MINDOUBLE  (M_LN2 * (DMINEXP - 1))
-#endif
-#ifndef LN_MAXDOUBLE
-#define LN_MAXDOUBLE  (M_LN2 * DMAXEXP)
-#endif
-
-r_8 TMatrix<r_8>::GausPiv(TMatrix<r_8>& a, TMatrix<r_8>& b)
-// Pivot de Gauss
-// * Attention: egcs impose que cette fonction soit mise dans le .cc
-//              avant ::Inverse() (car Inverse() l'utilise)
-// {TMatrix A(a); TMatrix B(b); return (r_8) TMatrix::GausPiv(A,B);}
-{
-uint_4 n = a.NRows();
-if(n!=b.NRows())
-  throw(SzMismatchError("TMatrix::GausPiv size mismatch\n"));
-// On fait une normalisation un peu brutale...
-double vmin=MAXDOUBLE;
-double vmax=0;
-for(uint_4 iii=0; iii<a.NRows(); iii++)
-  for(uint_4 jjj=0; jjj<a.NCols(); jjj++) {
-    double v = TMatrixRC<r_8>::Abs_Value(a(iii,jjj));
-    if(v>vmax) vmax = v;
-    if(v<vmin && v>0) vmin = v;
-}
-double nrm = sqrt(vmin*vmax);
-if(nrm > 1.e5 || nrm < 1.e-5) {
-  a /= nrm;
-  b /= nrm;
-  //cout << "normalisation matrice " << nrm << endl;
-} else nrm=1;
-
-double det = 1.0;
-if(nrm != 1) {
-  double ld = a.NRows() * log(nrm);
-  if (ld <= LN_MINDOUBLE || ld >= LN_MAXDOUBLE) {
-   // cerr << "TMatrix warning, overflow for det" << endl;
-  } else {
-    det = exp(ld);
-  }
-}
-
-TMatrixRC<r_8> pivRowa(a,TmatrixRow);
-TMatrixRC<r_8> pivRowb(b,TmatrixRow);
-
-for(uint_4 k=0; k<n-1; k++) {
-  uint_4 iPiv = a.Col(k).IMaxAbs(k);
-  if(iPiv != k) {
-    TMatrixRC<r_8> aIPiv(a.Row(iPiv));
-    TMatrixRC<r_8> aK(a.Row(k));
-    TMatrixRC<r_8>::Swap(aIPiv,aK);
-    TMatrixRC<r_8> bIPiv(b.Row(iPiv));
-    TMatrixRC<r_8> bK(b.Row(k));
-    TMatrixRC<r_8>::Swap(bIPiv,bK);
-  }
-  double pivot = a(k,k);
-  if (fabs(pivot) < 1.e-50) return 0.0;
-  //det *= pivot;
-  pivRowa.SetRow(k); // to avoid constructors
-  pivRowb.SetRow(k);
-  for (uint_4 i=k+1; i<n; i++) {
-    double r = -a(i,k)/pivot;
-    a.Row(i).LinComb(r, pivRowa); // + rapide que -= r * pivRowa
-    b.Row(i).LinComb(r, pivRowb);
-  }
-}
-det *= a(n-1, n-1);
-
-// on remonte
-for(uint_4 kk=n-1; kk>0; kk--) {
-  double pivot = a(kk,kk);
-  if (fabs(pivot) <= 1.e-50) return 0.0;
-  pivRowa.SetRow(kk); // to avoid constructors
-  pivRowb.SetRow(kk);
-  for(uint_4 jj=0; jj<kk; jj++) {
-    double r = -a(jj,kk)/pivot;
-    a.Row(jj).LinComb(r, pivRowa);
-    b.Row(jj).LinComb(r, pivRowb);
-  }
-}
-
-for(uint_4 l=0; l<n; l++) {
-  if (fabs(a(l,l)) <= 1.e-50) return 0.0;
-  b.Row(l) /= a(l,l);
-}
-
-return det;
-}
-
-TMatrix<r_8> TMatrix<r_8>::Inverse() const
-// Inversion
-{
-TMatrix<r_8> b(mNc,mNr); TMatrix<r_8> a(*this); b = 1.;
-if(fabs(TMatrix<r_8>::GausPiv(a,b)) < 1.e-50)
-  throw(MathExc("TMatrix Inverse() Singular OMatrix"));
-return b;
-}
-
 
 ///////////////////////////////////////////////////////////
@@ -379,184 +255,4 @@
 #include "tvector.h"
 
-template <class T> TMatrixRC<T>::TMatrixRC()
-: matrix(NULL), data(NULL), index(0), step(0)
-{}
-
-template <class T> TMatrixRC<T>::TMatrixRC(TMatrix<T>& m,TRCKind rckind,uint_4 ind)
-: matrix(&m), data(Org(m,rckind,ind)),
-  index(ind), step(Step(m,rckind)), kind(rckind)
-{
-if (kind == TmatrixDiag && m.mNc != m.mNr)
-  throw(SzMismatchError("TMatrixRC::TMatrixRC(...,TmatrixDiag,...) size mismatch\n"));
-}
-
-template <class T> int_4 TMatrixRC<T>::Next()
-{
-if (!matrix || kind==TmatrixDiag) return -1;
-index++;
-if(kind == TmatrixRow) {
-  if(index > (int_4)matrix->mNr) {index = (int_4)matrix->mNr; return -1;}
-  data += matrix->mNc;
-} else {
-  if (index > (int_4)matrix->mNc) {index = (int_4)matrix->mNc; 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->mNc;
-else data--;
-return index;
-}
-
-template <class T> int_4 TMatrixRC<T>::SetCol(int_4 c)
-{
-if(!matrix) return -1;
-if(c<0 || c>(int_4)matrix->mNc) return -1;
-kind = TmatrixCol;
-index = c;
-step = Step(*matrix, TmatrixCol);
-data = Org(*matrix, TmatrixCol, c);
-return c;
-}
-
-template <class T> int_4 TMatrixRC<T>::SetRow(int_4 r)
-{
-if(!matrix) return -1;
-if(r<0 && r>(int_4)matrix->mNr) return -1;
-kind = TmatrixRow;
-index = r;
-step = Step(*matrix, TmatrixRow);
-data = Org(*matrix, TmatrixRow, r);
-return r;
-}
-
-template <class T> int_4 TMatrixRC<T>::SetDiag()
-{
-if (!matrix) return -1;
-if (matrix->mNc != matrix->mNr)
-  throw(SzMismatchError("TMatrixRC::SetDiag size mismatch\n"));
-kind = TmatrixDiag;
-index = 0;
-step = Step(*matrix, TmatrixDiag);
-data = Org(*matrix, TmatrixDiag);
-return 0;
-}
-
-
-template <class T> TMatrixRC<T>& TMatrixRC<T>::operator = (const TMatrixRC<T>& rc)
-{
-matrix = rc.matrix;
-data   = rc.data;
-index  = rc.index;
-step   = rc.step;
-kind   = rc.kind;
-return *this;
-}
-
-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;
-}
-
-
-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> 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>
-TMatrixRC<T>& TMatrixRC<T>::LinComb(T a, T b, const TMatrixRC<T>& rc, uint_4 first)
-{
-if ( NElts() != rc.NElts() )
-  throw(SzMismatchError("TMatrixRC::LinComb size mismatch\n"));
-if ( kind != rc.kind )
-  throw(SzMismatchError("TMatrixRC::LinComb type mismatch\n"));
-for(uint_4 i=first; i<NElts(); i++) (*this)(i) = (*this)(i)*a + rc(i)*b;
-return *this;
-}
-
-template <class T>
-TMatrixRC<T>& TMatrixRC<T>::LinComb(T b, const TMatrixRC<T>& rc, uint_4 first)
-{
-if ( NElts() != rc.NElts() )
-  throw(SzMismatchError("TMatrixRC::LinComb size mismatch\n"));
-if ( kind != rc.kind )
-  throw(SzMismatchError("TMatrixRC::LinComb type mismatch\n"));
-for(uint_4 i=first; i<NElts(); i++) (*this)(i) += rc(i)*b;
-return *this;
-}
-
-template <class T> uint_4 TMatrixRC<T>::IMaxAbs(uint_4 first)
-{
-if (first>NElts())
-  throw(SzMismatchError("TMatrixRC::IMaxAbs size mismatch\n"));
-uint_4 imax = first;
-double vmax = Abs_Value((*this)(first));
-for(uint_4 i=first+1; i<NElts(); i++) {
-  double v = Abs_Value((*this)(i));
-  if(v > vmax) {vmax = v; imax = i;}
-}
-return imax;
-}
-
-template <class T>
-void TMatrixRC<T>::Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2)
-{
-if(rc1.NElts() != rc2.NElts())
-  throw(SzMismatchError("TMatrixRC::Swap size mismatch\n"));
-if(rc1.kind != rc2.kind)
-  throw(SzMismatchError("TMatrixRC::Swap type mismatch\n"));
-if(rc1.data == rc2.data) return;
-for(uint_4 i=0; i<rc1.NElts(); i++)
-  {T tmp = rc1(i); rc1(i) = rc2(i); rc2(i) = tmp;}
-}
 
 ///////////////////////////////////////////////////////////////
@@ -585,16 +281,4 @@
 #pragma define_template FIO_TMatrix< complex<r_4> >
 #pragma define_template FIO_TMatrix< complex<r_8> >
-// 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> >
 #endif
 
@@ -623,15 +307,3 @@
 template class FIO_TMatrix< complex<r_4> >;
 template class FIO_TMatrix< complex<r_8> >;
-// 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> >;
 #endif
Index: /trunk/SophyaLib/TArray/tmatrix.h
===================================================================
--- /trunk/SophyaLib/TArray/tmatrix.h	(revision 771)
+++ /trunk/SophyaLib/TArray/tmatrix.h	(revision 772)
@@ -14,11 +14,6 @@
 namespace SOPHYA {
 
-template <class T> class TVector;
-template <class T> class TMatrixRC;
-
 template <class T>
 class TMatrix : public AnyDataObj {
-  friend class TMatrixRC<T>;
-  friend class TVector<T>;
 public:
 
@@ -103,13 +98,4 @@
   TMatrix<T> Mul(const TMatrix<T>& b) const;
 
-  // Pivot de Gauss : diagonalise la matrice A, en effectuant les memes
-  // operations sur la matrice B
-  TMatrix<T> Inverse() const;
-  static T GausPiv(TMatrix<T>& A, TMatrix<T>& B);
-
-  // Acces aux rangees et colonnes
-  TMatrixRC<T> Row(uint_4 r) const;
-  TMatrixRC<T> Col(uint_4 c) const;
-  TMatrixRC<T> Diag() const;
 
 protected:
@@ -203,113 +189,4 @@
 { FIO_TMatrix<T> fio(&obj);  fio.Read(is);  return(is); }
 
-/////////////////////////////////////////////////////////////////////////
-// Classe de lignes/colonnes de matrices
-enum TRCKind {TmatrixRow=0, TmatrixCol=1, TmatrixDiag=2};
-template <class T>
-class TMatrixRC {
-  friend class TVector<T>;
-  friend class TMatrix<T>;
-public:
-  TMatrixRC();
-
-  virtual ~TMatrixRC() {}
-
-  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);
-
-protected:
-  TMatrixRC(TMatrix<T>& m, TRCKind kind, uint_4 index=0);
-  TMatrix<T>* matrix;
-  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());}
-
-  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.mNc;
-                      case TmatrixDiag : return m.mNc+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.mNc;
-                      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->mNc;
-                    case TmatrixCol  : return matrix->mNr;
-                    case TmatrixDiag : return matrix->mNc; }
-    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;
 
 } // Fin du namespace
Index: /trunk/SophyaLib/TArray/tvector.h
===================================================================
--- /trunk/SophyaLib/TArray/tvector.h	(revision 771)
+++ /trunk/SophyaLib/TArray/tvector.h	(revision 772)
@@ -55,21 +55,4 @@
    {return TVector<T>(a * ((TMatrix<T> const&)(b)));}
 
-// Resolution du systeme A*C = B
-inline r_8 LinSolveInPlace(TMatrix<r_8>& a, TVector<r_8>& b)
-{
-if(a.NCols() != b.NRows() || a.NCols() != a.NRows())
-  throw(SzMismatchError("LinSolveInPlace(TMatrix<r_8>,TVector<r_8>) size mismatch"));
-return TMatrix<r_8>::GausPiv(a,b);
-}
-
-// Resolution du systeme A*C = B, avec C retourne dans B
-inline r_8 LinSolve(const TMatrix<r_8>& a, const TVector<r_8>& b, TVector<r_8>& c)
-{ 
-if(a.NCols() != b.NRows() || a.NCols() != a.NRows())
-  throw(SzMismatchError("LinSolve(TMatrix<r_8>,TVector<r_8>) size mismatch"));
-c = b;
-TMatrix<r_8> a1(a);
-return TMatrix<r_8>::GausPiv(a1,c);
-}
 
 // Typedef pour simplifier et compatibilite Peida
