Index: /trunk/SophyaLib/TArray/Makefile
===================================================================
--- /trunk/SophyaLib/TArray/Makefile	(revision 762)
+++ /trunk/SophyaLib/TArray/Makefile	(revision 762)
@@ -0,0 +1,37 @@
+include ../Mgr/Makefile.h
+ 
+all: $(LIB)libTArray.a
+clean:
+	rm -f $(OBJ)tarray.o $(OBJ)tarrinit.o $(OBJ)tmatrix.o $(OBJ)tvector.o
+	rm -f $(LIB)libTArray.a
+$(LIB)libTArray.a : $(OBJ)tarray.o $(OBJ)tarrinit.o $(OBJ)tmatrix.o $(OBJ)tvector.o
+	$(ARCXX) $(ARCXXFLAGS) $@ $?
+ifeq ($(CXX),cxx)
+	$(CXX) $? $(CPPFLAGS) $(CXXFLAGS) -L$(LIB) -lTArray -L$(EXTLIBPATH) -o $(OBJ)xx.x -Hf
+endif 
+ 
+$(OBJ)tarray.o: tarray.cc
+$(OBJ)tarrinit.o: tarrinit.cc $(INC)/machdefs.h \
+ $(INC)/pclassids.h tarrinit.h \
+ $(INC)/sophyainit.h tmatrix.h \
+ $(INC)/ppersist.h \
+ $(INC)/pexceptions.h \
+ $(INC)/gnumd5.h \
+ $(INC)/anydataobj.h \
+ $(INC)/ndatablock.h tvector.h
+$(OBJ)tmatrix.o: tmatrix.cc $(INC)/machdefs.h \
+ $(INC)/pclassids.h \
+ $(INC)/pexceptions.h tmatrix.h \
+ $(INC)/ppersist.h \
+ $(INC)/gnumd5.h \
+ $(INC)/anydataobj.h \
+ $(INC)/ndatablock.h \
+ $(INC)/objfio.h tvector.h
+$(OBJ)tvector.o: tvector.cc $(INC)/machdefs.h \
+ $(INC)/pclassids.h \
+ $(INC)/pexceptions.h tvector.h tmatrix.h \
+ $(INC)/ppersist.h \
+ $(INC)/gnumd5.h \
+ $(INC)/anydataobj.h \
+ $(INC)/ndatablock.h \
+ $(INC)/objfio.h
Index: /trunk/SophyaLib/TArray/objlist.list
===================================================================
--- /trunk/SophyaLib/TArray/objlist.list	(revision 762)
+++ /trunk/SophyaLib/TArray/objlist.list	(revision 762)
@@ -0,0 +1,3 @@
+tarrinit.o
+tmatrix.o
+tvector.o
Index: /trunk/SophyaLib/TArray/tarrinit.cc
===================================================================
--- /trunk/SophyaLib/TArray/tarrinit.cc	(revision 762)
+++ /trunk/SophyaLib/TArray/tarrinit.cc	(revision 762)
@@ -0,0 +1,55 @@
+// Classe d'initialisation du module TArray (template Arrays)
+#include "machdefs.h"
+#include "tarrinit.h"
+
+#include "tmatrix.h"
+#include "tvector.h"
+
+
+int TArrayInitiator::FgInit = 0;
+
+TArrayInitiator::TArrayInitiator()
+  : SophyaInitiator()
+{
+  FgInit++;
+  if (FgInit > 1)  return;
+
+//   Enregistrement des classes PPersist du modules TArray
+  
+
+  PPRegister(FIO_TMatrix<uint_1>);
+  PPRegister(FIO_TMatrix<uint_2>);
+  PPRegister(FIO_TMatrix<int_2>);
+  PPRegister(FIO_TMatrix<int_4>);
+  PPRegister(FIO_TMatrix<int_8>);
+  PPRegister(FIO_TMatrix<uint_8>);
+  PPRegister(FIO_TMatrix<r_4>);
+  PPRegister(FIO_TMatrix<r_8>);
+  PPRegister(FIO_TMatrix< complex<r_4> >);
+  PPRegister(FIO_TMatrix< complex<r_8> >);
+
+  PPRegister(FIO_TVector<uint_1>);
+  PPRegister(FIO_TVector<uint_2>);
+  PPRegister(FIO_TVector<int_2>);
+  PPRegister(FIO_TVector<int_4>);
+  PPRegister(FIO_TVector<int_8>);
+  PPRegister(FIO_TVector<uint_8>);
+  PPRegister(FIO_TVector<r_4>);
+  PPRegister(FIO_TVector<r_8>);
+  PPRegister(FIO_TVector< complex<r_4> >);
+  PPRegister(FIO_TVector< complex<r_8> >);
+
+
+
+}
+
+TArrayInitiator::~TArrayInitiator()
+{
+  FgInit--;
+}
+
+
+// On met un objet initiator en statique, pour les loaders qui savent 
+// appeler le constructeur des objets statiques   Reza 08/98
+static TArrayInitiator tarrayinit;
+
Index: /trunk/SophyaLib/TArray/tarrinit.h
===================================================================
--- /trunk/SophyaLib/TArray/tarrinit.h	(revision 762)
+++ /trunk/SophyaLib/TArray/tarrinit.h	(revision 762)
@@ -0,0 +1,27 @@
+// This may look like C code, but it is really -*- C++ -*-
+// Classe d'initialiseur du module TArray
+
+#ifndef TARRAYINIT_H_SEEN
+#define TARRAYINIT_H_SEEN
+
+
+// Classe d''initialisation pour le module de NTools
+#include "sophyainit.h"
+
+namespace SOPHYA {
+
+class TArrayInitiator : public SophyaInitiator {
+private:
+  static int FgInit;
+public:
+  			TArrayInitiator();
+  virtual		~TArrayInitiator();
+};
+
+#undef SophyaInit
+#define SophyaInit() TArrayInitiator tarrayinitiator
+
+} // Fin du namespace
+
+#endif
+
Index: /trunk/SophyaLib/TArray/tmatrix.cc
===================================================================
--- /trunk/SophyaLib/TArray/tmatrix.cc	(revision 762)
+++ /trunk/SophyaLib/TArray/tmatrix.cc	(revision 762)
@@ -0,0 +1,637 @@
+// $Id: tmatrix.cc,v 1.1.1.1 2000-03-02 16:48:24 ansari Exp $
+//                         C.Magneville          04/99
+#include "machdefs.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream.h>
+#include <values.h>
+#include <complex>
+#include "pexceptions.h"
+#include "tmatrix.h"
+#include "objfio.h"
+
+////////////////////////////////////////////////////////////////
+//**** Createur, Destructeur
+
+template <class T>
+TMatrix<T>::TMatrix()
+// Constructeur par defaut.
+: mNr(0), mNc(0), mNDBlock()
+{
+}
+
+template <class T>
+TMatrix<T>::TMatrix(uint_4 r,uint_4 c)
+// Construit une matrice de r lignes et c colonnes.
+: mNr(r), mNc(c), mNDBlock(r*c)
+{
+}
+
+template <class T>
+TMatrix<T>::TMatrix(uint_4 r,uint_4 c,T* values,Bridge* br)
+// Construit une matrice de r lignes et c colonnes. On fournit
+// le tableau des valeurs et eventuellement un Bridge.
+: mNr(r), mNc(c), mNDBlock(r*c,values,br)
+{
+}
+
+template <class T>
+TMatrix<T>::TMatrix(const TMatrix<T>& a)
+// Constructeur par copie (partage si "a" temporaire).
+: mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock)
+{
+}
+
+template <class T>
+TMatrix<T>::TMatrix(const TMatrix<T>& a,bool share)
+// Constructeur par copie avec possibilite de forcer le partage ou non.
+: mNr(a.mNr), mNc(a.mNc), mNDBlock(a.mNDBlock,share)
+{
+}
+
+template <class T>
+TMatrix<T>::~TMatrix()
+// Destructeur
+{
+}
+
+////////////////////////////////////////////////////////////////
+// Operations matricielles
+template <class T>
+TMatrix<T> TMatrix<T>::Transpose(void) const
+// Transposition
+{
+TMatrix<T> a; a.Clone(*this); a.SetTemp(true);
+a.mNr = mNc; a.mNc = mNr;
+{for(uint_4 i=0; i<mNr; i++)
+  for(uint_4 j=0; j<mNc; j++) {
+    a(j,i) = (*this)(i,j);
+}}
+return a;
+}
+
+////////////////////////////////////////////////////////////////
+//**** Impression
+
+template <class T>
+void TMatrix<T>::Print(ostream& os,int lp
+                 ,uint_4 i0,uint_4 ni,uint_4 j0,uint_4 nj) const
+// Impression de la sous-matrice (i:i+ni-1,i:j+nj-1)
+{
+os<<"TMatrix::Print("<<mNr<<","<<mNc<<")"<<endl;
+if(lp>0)
+  {os<<"       this="<<this<<endl; mNDBlock.Print(0,0);}
+if(mNr==0 || mNc==0) return;
+if(i0>=mNr || j0>=mNc || ni==0 || nj==0) return;
+uint_4 i1 = i0+ni; if(i1>mNr) i1 = mNr;
+uint_4 j1 = j0+nj; if(j1>mNc) j1 = mNc;
+for(uint_4 i=i0;i<i1;i++) {
+  for(uint_4 j=j0;j<j1;j++) cout<<" "<<(*this)(i,j);
+  cout<<endl;
+}
+}
+
+////////////////////////////////////////////////////////////////
+//**** Surcharge de *= (INPLACE): TMatrix *= TMatrix;
+
+template <class T>
+TMatrix<T>& TMatrix<T>::operator *= (const TMatrix<T>& a)
+// A = A*B  -> A(n,m) = A(n,m)*B(m,m)
+{
+uint_4 ndata = mNr*mNc;
+if(ndata==0 || mNc!=a.mNr || a.mNr!=a.mNc)
+  throw(SzMismatchError("TMatrix::operator*=A size mismatch"));
+// A(i,j) = Sum(k) A(i,k)*B(k,j) ... il faut sauver la ligne "i" de A
+// Vecteur oi   : vecteur ou est sauve la ligne i de la matrice *this
+//                oi,oe = pointeur de debut et de fin du vecteur temporaire
+//                oij   = pointeur parcourant le vecteur oi
+// Matrice *this: i     = pointeur du debut de la ligne i
+//                ij    = pointeur parcourant la ligne i
+// Matrice a    : aj    = pointeur de debut de la colonne j
+//                aji   = pointeur parcourant la colonne j
+T* oi = new T[mNc]; T* oe = oi+mNc;
+for(T *i=Data(); i<Data()+ndata; i+=mNc) {
+  {for(T *oij=oi,*ij=i; oij<oe;) *oij++ = *ij++;}
+  {for(T *ij=i,*aj=const_cast<T *>(a.Data()); aj<a.Data()+a.mNc; ij++,aj++) {
+    T sum = 0;
+    for(T *oij=oi,*aji=aj; oij<oe; oij++,aji+=a.mNc) sum += *oij * *aji;
+    *ij = sum;
+  }}
+}
+delete [] oi;
+return *this;
+}
+
+////////////////////////////////////////////////////////////////
+//**** Pour surcharge d'operateurs C = A (+,-,*,/) B
+
+template <class T> TMatrix<T> TMatrix<T>::Add(const TMatrix<T>& b) const
+{
+if(mNr!=b.mNr || mNc!=b.mNc)
+    throw(SzMismatchError("TMatrix operator C=A+B size mismatch\n"));
+TMatrix<T> result; result.SetTemp(true); result.mNr=mNr; result.mNc=mNc;
+result.mNDBlock = mNDBlock+b.mNDBlock;
+return result;
+}
+
+template <class T> TMatrix<T> TMatrix<T>::Sub(const TMatrix<T>& b) const
+{
+if(mNr!=b.mNr || mNc!=b.mNc)
+    throw(SzMismatchError("TMatrix operator C=A-B size mismatch\n"));
+TMatrix<T> result; result.SetTemp(true); result.mNr=mNr; result.mNc=mNc;
+result.mNDBlock = mNDBlock-b.mNDBlock;
+return result;
+}
+
+template <class T> TMatrix<T> TMatrix<T>::Mul(const TMatrix<T>& b) const
+// C = A(this)*B : Cij = Aik Bkj (allocation forcee dans tous les cas)
+{
+if(mNr==0 || mNc==0 || b.mNc==0 || mNc!=b.mNr)
+  throw(SzMismatchError("TMatrix operator C=A*B size mismatch\n"));
+TMatrix<T> r; r.SetTemp(true); r.ReSize(mNr,b.mNc);
+T *ai,*aik,*bj,*bkj,*ri,*rij;
+for(ri=const_cast<T *>(r.Data()),ai=const_cast<T *>(Data());
+    ri<r.Data()+r.mNr*r.mNc;ri+=r.mNc,ai+=mNc) {
+  for(rij=ri,bj=const_cast<T *>(b.Data());rij<ri+r.mNc;rij++,bj++) {
+    *rij = 0;
+    for(aik=ai,bkj=bj;aik<ai+mNc;aik++,bkj+=b.mNc) *rij += *aik * *bkj;
+  }
+}
+return r;
+}
+
+////////////////////////////////////////////////////////////////
+//**** 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;
+}
+
+
+///////////////////////////////////////////////////////////
+// --------------------------------------------------------
+//   Les objets delegues pour la gestion de persistance 
+// --------------------------------------------------------
+///////////////////////////////////////////////////////////
+
+template <class T>
+FIO_TMatrix<T>::FIO_TMatrix()
+{
+dobj=new TMatrix<T>;
+ownobj=true;
+}
+
+template <class T>
+FIO_TMatrix<T>::FIO_TMatrix(string const & filename) 
+{
+dobj=new TMatrix<T>;
+ownobj=true; 
+Read(filename);
+}
+
+template <class T>
+FIO_TMatrix<T>::FIO_TMatrix(const TMatrix<T> & obj) 
+{ 
+dobj = new TMatrix<T>(obj);
+ownobj=true; 
+}
+
+template <class T>
+FIO_TMatrix<T>::FIO_TMatrix(TMatrix<T> * obj) 
+{ 
+dobj = obj;
+ownobj=false; 
+}
+
+template <class T>
+FIO_TMatrix<T>::~FIO_TMatrix()
+{
+if (ownobj && dobj) delete dobj;
+}
+
+template <class T>
+AnyDataObj* FIO_TMatrix<T>::DataObj()
+{
+return(dobj);
+}
+
+template <class T>
+void FIO_TMatrix<T>::SetDataObj(AnyDataObj & o)
+{
+TMatrix<T> * po = dynamic_cast< TMatrix<T> * >(&o);
+if (po == NULL) return;
+if (ownobj && dobj) delete dobj;
+dobj = po; ownobj = false;
+} 
+
+template <class T>
+void FIO_TMatrix<T>::ReadSelf(PInPersist& is)
+{
+// On lit les 3 premiers uint_4
+//  0: Numero de version,  1 : NRows,  2 : NCol
+uint_4 itab[3];
+is.Get(itab,3);
+if (dobj == NULL) dobj = new TMatrix<T>(itab[1],itab[2]);
+else dobj->ReSize(itab[1],itab[2]);
+// On lit le NDataBlock
+FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
+fio_nd.Read(is);
+}
+
+template <class T>
+void FIO_TMatrix<T>::WriteSelf(POutPersist& os) const
+{
+if (dobj == NULL)   return;
+//  On ecrit 3 uint_4 .... 
+//  0: Numero de version,  1 : NRows,  2 : NCol
+uint_4 itab[3];
+ itab[0] = 1;  // Numero de version a 1
+itab[1] = dobj->NRows();
+itab[2] = dobj->NCols();
+os.Put(itab,3);
+// On ecrit le NDataBlock
+FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
+fio_nd.Write(os);
+}
+
+////////////////////////////////////////////////////////////////
+// -------------------------------------------------------------
+//   La classe de gestion des lignes et colonnes d'une matrice
+// -------------------------------------------------------------
+////////////////////////////////////////////////////////////////
+#include "tvector.h"
+
+template <class T> 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;}
+}
+
+///////////////////////////////////////////////////////////////
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template TMatrix<uint_1>
+#pragma define_template TMatrix<uint_2>
+#pragma define_template TMatrix<int_2>
+#pragma define_template TMatrix<int_4>
+#pragma define_template TMatrix<int_8>
+#pragma define_template TMatrix<uint_4>
+#pragma define_template TMatrix<uint_8>
+#pragma define_template TMatrix<r_4>
+#pragma define_template TMatrix<r_8>
+#pragma define_template TMatrix< complex<r_4> > 
+#pragma define_template TMatrix< complex<r_8> > 
+// Instances des delegues FileIO (PPersist)
+#pragma define_template FIO_TMatrix<uint_1>
+#pragma define_template FIO_TMatrix<uint_2>
+#pragma define_template FIO_TMatrix<int_2>
+#pragma define_template FIO_TMatrix<int_4>
+#pragma define_template FIO_TMatrix<int_8>
+#pragma define_template FIO_TMatrix<uint_4>
+#pragma define_template FIO_TMatrix<uint_8>
+#pragma define_template FIO_TMatrix<r_8>
+#pragma define_template FIO_TMatrix<r_4>
+#pragma define_template FIO_TMatrix< complex<r_4> >
+#pragma define_template FIO_TMatrix< complex<r_8> >
+// 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
+
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class TMatrix<uint_1>;
+template class TMatrix<uint_2>;
+template class TMatrix<int_2>;
+template class TMatrix<int_4>;
+template class TMatrix<int_8>;
+template class TMatrix<uint_4>;
+template class TMatrix<uint_8>;
+template class TMatrix<r_4>;
+template class TMatrix<r_8>;
+template class TMatrix< complex<r_4> >;
+template class TMatrix< complex<r_8> >;
+// Instances des delegues FileIO (PPersist)
+template class FIO_TMatrix<uint_1>;
+template class FIO_TMatrix<uint_2>;
+template class FIO_TMatrix<int_2>;
+template class FIO_TMatrix<int_4>;
+template class FIO_TMatrix<int_8>;
+template class FIO_TMatrix<uint_4>;
+template class FIO_TMatrix<uint_8>;
+template class FIO_TMatrix<r_8>;
+template class FIO_TMatrix<r_4>;
+template class FIO_TMatrix< complex<r_4> >;
+template class FIO_TMatrix< complex<r_8> >;
+// 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 762)
+++ /trunk/SophyaLib/TArray/tmatrix.h	(revision 762)
@@ -0,0 +1,317 @@
+// This may look like C code, but it is really -*- C++ -*-
+//                         C.Magneville          04/99
+#ifndef TMatrix_SEEN
+#define TMatrix_SEEN
+
+#include "machdefs.h"
+#include <stdio.h>
+#include <iostream.h>
+#include <complex>
+#include "ppersist.h"
+#include "anydataobj.h"
+#include "ndatablock.h"
+
+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:
+
+  // Creation / destruction 
+  TMatrix();
+  TMatrix(uint_4 r,uint_4 c);
+  TMatrix(uint_4 r,uint_4 c,T* values,Bridge* br=NULL);
+  TMatrix(const TMatrix<T>& a);
+  TMatrix(const TMatrix<T>& a,bool share);
+  virtual ~TMatrix();
+
+  // Temporaire?
+  inline bool IsTemp(void) const {return mNDBlock.IsTemp();}
+  inline void SetTemp(bool temp=false) const {mNDBlock.SetTemp(temp);}
+
+  // Gestion taille/Remplissage
+  inline void Clone(const TMatrix<T>& a)  // Clone: copie des donnees de "a"
+              {mNDBlock.Clone(a.mNDBlock); mNr = a.mNr; mNc = a.mNc;}
+  inline void Reset(T v=0) {mNDBlock.Reset(v);}
+  inline void ReSize(uint_4 r,uint_4 c)  // Reallocation de place
+       {if(r==0||c==0) throw(SzMismatchError("TMatrix::ReSize r ou c==0\n"));
+        mNDBlock.ReSize(r*c); mNr = r; mNc = c;}
+  inline void Realloc(uint_4 r,uint_4 c,bool force=false)
+       {if(r==0||c==0) throw(SzMismatchError("TMatrix::Realloc r ou c==0\n"));
+        mNDBlock.Realloc(r*c,force); mNr = r; mNc = c;}
+
+  // Informations pointeur/data
+  inline uint_4 NRows() const {return mNr;}
+  inline uint_4 NCols() const {return mNc;}
+  inline uint_4 NCol() const {return mNc;} // back-compat Peida
+  inline T const& operator()(uint_4 r,uint_4 c) const
+                            {return *(mNDBlock.Begin()+r*mNc+c);}
+  inline T&       operator()(uint_4 r,uint_4 c)
+                            {return *(mNDBlock.Begin()+r*mNc+c);}
+  inline T const& operator[](uint_4 ip) const
+                            {return *(mNDBlock.Begin()+ip);}
+  inline T&       operator[](uint_4 ip)
+                            {return *(mNDBlock.Begin()+ip);}
+  inline       T* Data()       {return mNDBlock.Begin();}
+  inline const T* Data() const {return mNDBlock.Begin();}
+  inline       NDataBlock<T>& DataBlock()       {return mNDBlock;}
+  inline const NDataBlock<T>& DataBlock() const {return mNDBlock;}
+
+  // Operations matricielles
+  TMatrix<T> Transpose(void) const;
+
+  // Operateur d'affectation
+  // A = x (matrice diagonale x*Identite)
+  inline TMatrix<T>& operator = (T x)
+    {if(mNr!=mNc || mNr==0) throw(SzMismatchError("TMatrix::operator= mNc!=mNr ou ==0\n"));
+     for(uint_4 r=0;r<mNr;r++) for(uint_4 c=0;c<mNc;c++) (*this)(r,c)=(r==c)?x:(T)0;
+     return *this;}
+  // A = B : partage les donnees si "a" est temporaire, clone sinon.
+  inline TMatrix<T>& operator = (const TMatrix<T>& a)
+                {if(this == &a) return *this; CloneOrShare(a); return *this;}
+
+  // Impression
+  void Print(ostream& os,int lp=0,uint_4 i0=0,uint_4 ni=10,uint_4 j0=0,uint_4 nj=10) const;
+  inline void Print(int lp=0,uint_4 i0=0,uint_4 ni=10,uint_4 j0=0,uint_4 nj=10) const
+              {Print(cout,lp,i0,ni,j0,nj);}
+
+  // Surcharge d'operateurs INPLACE: A (+=,-=,*=,/=) (T) x
+  inline TMatrix<T>& operator += (T b) {mNDBlock += b; return *this;}
+  inline TMatrix<T>& operator -= (T b) {mNDBlock -= b; return *this;}
+  inline TMatrix<T>& operator *= (T b) {mNDBlock *= b; return *this;}
+  inline TMatrix<T>& operator /= (T b) {mNDBlock /= b; return *this;}
+
+  // Surcharge d'operateurs INPLACE: A (+=,-=,*=,/=) B
+  inline TMatrix<T>& operator += (const TMatrix<T>& a)
+              {if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
+                  throw(SzMismatchError("TMatrix::operator+=A size mismatch"));
+               mNDBlock += a.mNDBlock; return *this;}
+  inline TMatrix<T>& operator -= (const TMatrix<T>& a)
+              {if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
+                  throw(SzMismatchError("TMatrix::operator-=A size mismatch"));
+               mNDBlock -= a.mNDBlock; return *this;}
+  TMatrix<T>& operator *= (const TMatrix<T>& a);
+
+  // Pour surcharge d'operateurs C = A (+,-,*) B
+  TMatrix<T> Add(const TMatrix<T>& b) const;
+  TMatrix<T> Sub(const TMatrix<T>& b) const;
+  TMatrix<T> Mul(const TMatrix<T>& b) const;
+
+  // 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:
+  // partage les donnees si "a" temporaire, clone sinon.
+  inline void CloneOrShare(const TMatrix<T>& a)
+              {mNDBlock.CloneOrShare(a.mNDBlock); mNr=a.mNr; mNc=a.mNc;}
+  // Share: partage les donnees de "a"
+  inline void Share(const TMatrix<T>& a)
+              {mNDBlock.Share(a.mNDBlock); mNr=a.mNr; mNc=a.mNc;}
+
+  uint_4 mNr,mNc;
+  NDataBlock<T> mNDBlock;
+};
+
+////////////////////////////////////////////////////////////////
+// Impression
+
+template <class T>
+inline ostream& operator << (ostream& os, const TMatrix<T>& a)
+                            {a.Print(os); return(os);}
+
+////////////////////////////////////////////////////////////////
+// Surcharge d'operateurs A (+,-,*,/) (T) x
+
+template <class T> inline TMatrix<T> operator + (const TMatrix<T>& a, T b)
+    {TMatrix<T> result(a); result.SetTemp(true); result += b; return result;}
+
+template <class T> inline TMatrix<T> operator + (T b,const TMatrix<T>& a)
+    {TMatrix<T> result(a); result.SetTemp(true); result += b; return result;}
+
+template <class T> inline TMatrix<T> operator - (const TMatrix<T>& a, T b)
+    {TMatrix<T> result(a); result.SetTemp(true); result -= b; return result;}
+
+template <class T> inline TMatrix<T> operator - (T b,const TMatrix<T>& a)
+    {TMatrix<T> result(a); result.SetTemp(true);
+     result.DataBlock() = b-result.DataBlock(); return result;}
+
+template <class T> inline TMatrix<T> operator * (const TMatrix<T>& a, T b)
+    {TMatrix<T> result(a); result.SetTemp(true); result *= b; return result;}
+
+template <class T> inline TMatrix<T> operator * (T b,const TMatrix<T>& a)
+    {TMatrix<T> result(a); result.SetTemp(true); result *= b; return result;}
+
+template <class T> inline TMatrix<T> operator / (const TMatrix<T>& a, T b)
+    {TMatrix<T> result(a); result.SetTemp(true); result /= b; return result;}
+
+////////////////////////////////////////////////////////////////
+// Surcharge d'operateurs C = A (+,-,*,/) B
+
+template <class T>
+inline TMatrix<T> operator + (const TMatrix<T>& a,const TMatrix<T>& b)
+                  {return a.Add(b);}
+
+template <class T>
+inline TMatrix<T> operator - (const TMatrix<T>& a,const TMatrix<T>& b)
+                  {return a.Sub(b);}
+
+template <class T>
+inline TMatrix<T> operator * (const TMatrix<T>& a,const TMatrix<T>& b)
+                  {return a.Mul(b);}
+
+////////////////////////////////////////////////////////////////
+// Typedef pour simplifier et compatibilite Peida
+typedef TMatrix<r_8> Matrix;
+
+/////////////////////////////////////////////////////////////////////////
+// Classe pour la gestion de persistance
+template <class T>
+class FIO_TMatrix : public  PPersist  {
+public:
+  FIO_TMatrix();
+  FIO_TMatrix(string const & filename); 
+  FIO_TMatrix(const TMatrix<T> & obj);
+  FIO_TMatrix(TMatrix<T> * obj);
+  virtual ~FIO_TMatrix();
+  virtual AnyDataObj* DataObj();
+  virtual void        SetDataObj(AnyDataObj & o);
+  inline operator TMatrix<T>() { return(*dobj); }
+protected :
+  virtual void ReadSelf(PInPersist&);           
+  virtual void WriteSelf(POutPersist&) const;  
+  TMatrix<T> * dobj;
+  bool ownobj;
+};
+
+template <class T>
+inline POutPersist& operator << (POutPersist& os, TMatrix<T> & obj)
+{ FIO_TMatrix<T> fio(&obj);  fio.Write(os);  return(os); }
+template <class T>
+inline PInPersist& operator >> (PInPersist& is, TMatrix<T> & obj)
+{ FIO_TMatrix<T> fio(&obj);  fio.Read(is);  return(is); }
+
+/////////////////////////////////////////////////////////////////////////
+// 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
+
+#endif
Index: /trunk/SophyaLib/TArray/triangmtx.h
===================================================================
--- /trunk/SophyaLib/TArray/triangmtx.h	(revision 762)
+++ /trunk/SophyaLib/TArray/triangmtx.h	(revision 762)
@@ -0,0 +1,65 @@
+// This may look like C code, but it is really -*- C++ -*-
+/*! Class for inferior triangular matrix (base class for the class Alm) */
+
+#ifndef TRIANGMTX_H_SEEN
+#define TRIANGMTX_H_SEEN
+
+#include "ndatablock.h"
+#include "pexceptions.h"
+
+namespace SOPHYA {
+
+template <class T>
+class TriangularMatrix 
+  {
+    
+  public :
+
+TriangularMatrix()   {};
+  /* instanciate a triangular matrix from the number of rows */
+TriangularMatrix(int rowSize)  : long_diag_((uint_4)rowSize) {elem_.ReSize((uint_4) (rowSize*(rowSize+1)/2) ); };
+TriangularMatrix(const TriangularMatrix<T>& a,  bool share=false)  : elem_(a.elem_, share),  long_diag_(a.long_diag_) {;}
+/*! resize the matrix with a new number of rows */
+inline void ReSizeRow(int rowSize) 
+  {
+    long_diag_=(uint_4)rowSize;
+    elem_.ReSize(long_diag_*(long_diag_+1)/2);
+  }
+ inline void SetTemp(bool temp=false) const {elem_.SetTemp(temp);}
+
+inline TriangularMatrix<T>& operator = (const TriangularMatrix<T>& a)
+{
+  elem_=a.elem_;
+  long_diag_ = a.long_diag_;
+  return *this;
+}
+inline T& operator()(int l, int m) 
+  {
+      return  elem_(adr_ij(l,m));
+  }
+inline T const& operator()(int l, int m) const 
+  {
+      return *(elem_.Begin()+ adr_ij(l,m));
+  }
+ inline  int_4  rowNumber() const {return (int_4)long_diag_;}
+  private: 
+ /*! compute the address of an element in the single array representing the matrix */
+inline uint_4 adr_ij(int i,int j) const 
+{
+  int adr= i*(i+1)/2+j;
+  if ( adr >= elem_.Size() || adr <0 )
+{
+ cout << " attention depassement dans triangularMatrix " << endl;
+ cout << " l= " << i << " m= " << j << " tableau reserve longueur " << elem_.Size() << endl; 
+}
+  return(i*(i+1)/2+j);
+}
+
+ uint_4 long_diag_;
+ NDataBlock<T> elem_;
+
+  };
+
+}   // namespace SOPHYA
+
+#endif
Index: /trunk/SophyaLib/TArray/tvector.cc
===================================================================
--- /trunk/SophyaLib/TArray/tvector.cc	(revision 762)
+++ /trunk/SophyaLib/TArray/tvector.cc	(revision 762)
@@ -0,0 +1,188 @@
+// $Id: tvector.cc,v 1.1.1.1 2000-03-02 16:48:24 ansari Exp $
+//                         C.Magneville          04/99
+#include "machdefs.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream.h>
+#include <complex>
+#include "pexceptions.h"
+#include "tvector.h"
+#include "objfio.h"
+
+template <class T>
+TVector<T>::TVector(uint_4 n)
+// Constructeur d'un vecteur de "n" elements
+: TMatrix<T>(n,1)
+{
+}
+
+template <class T>
+TVector<T>::TVector(uint_4 n, T* values,Bridge* br)
+// Construit un vecteur de n elements. On fournit
+// le tableau des valeurs et eventuellement un Bridge.
+: TMatrix<T>(n,1,values,br)
+{
+}
+
+template <class T>
+TVector<T>::TVector(const TVector<T>& v)
+// Constructeur par copie (partage si "v" temporaire).
+: TMatrix<T>(v)
+{
+}
+
+template <class T>
+TVector<T>::TVector(const TVector<T>& v,bool share)
+// Constructeur par copie avec possibilite de forcer le partage ou non.
+: TMatrix<T>(v,share)
+{
+}
+
+template <class T>
+TVector<T>::TVector(const TMatrix<T>& a)
+// Constructeur a partir d'une matrice "a" de dimension (n,1)
+: TMatrix<T>(a)
+{
+if(a.NCols() != 1)
+  throw(SzMismatchError("TVector(const TMatrix<T>& a) size mismatch, ncol!=1"));
+}
+
+
+///////////////////////////////////////////////////////////
+// --------------------------------------------------------
+//   Les objets delegues pour la gestion de persistance 
+// --------------------------------------------------------
+///////////////////////////////////////////////////////////
+
+template <class T>
+FIO_TVector<T>::FIO_TVector()
+{
+dobj=new TVector<T>;
+ownobj=true;
+}
+
+template <class T>
+FIO_TVector<T>::FIO_TVector(string const & filename) 
+{
+dobj=new TVector<T>;
+ownobj=true; 
+Read(filename);
+}
+
+template <class T>
+FIO_TVector<T>::FIO_TVector(const TVector<T> & obj) 
+{ 
+dobj = new TVector<T>(obj);
+ownobj=true; 
+}
+
+template <class T>
+FIO_TVector<T>::FIO_TVector(TVector<T> * obj) 
+{ 
+dobj = obj;
+ownobj=false; 
+}
+
+template <class T>
+FIO_TVector<T>::~FIO_TVector()
+{
+if (ownobj && dobj) delete dobj;
+}
+
+template <class T>
+AnyDataObj* FIO_TVector<T>::DataObj()
+{
+return(dobj);
+}
+
+template <class T>
+void FIO_TVector<T>::SetDataObj(AnyDataObj & o)
+{
+TVector<T> * po = dynamic_cast< TVector<T> * >(&o);
+if (po == NULL) return;
+if (ownobj && dobj) delete dobj;
+dobj = po; ownobj = false;
+} 
+
+template <class T>
+void FIO_TVector<T>::ReadSelf(PInPersist& is)
+{
+// On lit les 3 premiers uint_4
+//  0: Numero de version,  1 : NRows=NElts,  2 : NCol=1
+uint_4 itab[3];
+is.Get(itab,3);
+if (dobj == NULL) dobj = new TVector<T>(itab[1]);
+else dobj->ReSize(itab[1]);
+// On lit le NDataBlock
+FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
+fio_nd.Read(is);
+}
+
+template <class T>
+void FIO_TVector<T>::WriteSelf(POutPersist& os) const
+{
+if (dobj == NULL)   return;
+//  On ecrit 3 uint_4 .... 
+//  0: Numero de version,  1 : NRows=NElts,  2 : NCol=1
+uint_4 itab[3];
+ itab[0] = 1;  // Numero de version a 1
+itab[1] = dobj->NElts();
+itab[2] = 1;
+os.Put(itab,3);
+// On ecrit le NDataBlock
+FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
+fio_nd.Write(os);
+}
+
+///////////////////////////////////////////////////////////////
+#ifdef __CXX_PRAGMA_TEMPLATES__
+#pragma define_template TVector<uint_1>
+#pragma define_template TVector<uint_2>
+#pragma define_template TVector<int_2>
+#pragma define_template TVector<int_4>
+#pragma define_template TVector<int_8>
+#pragma define_template TVector<uint_4>
+#pragma define_template TVector<uint_8>
+#pragma define_template TVector<r_4>
+#pragma define_template TVector<r_8>
+#pragma define_template TVector< complex<r_4> > 
+#pragma define_template TVector< complex<r_8> > 
+// Instances des delegues FileIO (PPersist)
+#pragma define_template FIO_TVector<uint_1>
+#pragma define_template FIO_TVector<uint_2>
+#pragma define_template FIO_TVector<int_2>
+#pragma define_template FIO_TVector<int_4>
+#pragma define_template FIO_TVector<int_8>
+#pragma define_template FIO_TVector<uint_4>
+#pragma define_template FIO_TVector<uint_8>
+#pragma define_template FIO_TVector<r_8>
+#pragma define_template FIO_TVector<r_4>
+#pragma define_template FIO_TVector< complex<r_4> >
+#pragma define_template FIO_TVector< complex<r_8> >
+#endif
+
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class TVector<uint_1>;
+template class TVector<uint_2>;
+template class TVector<int_2>;
+template class TVector<int_4>;
+template class TVector<int_8>;
+template class TVector<uint_4>;
+template class TVector<uint_8>;
+template class TVector<r_4>;
+template class TVector<r_8>;
+template class TVector< complex<r_4> >;
+template class TVector< complex<r_8> >;
+// Instances des delegues FileIO (PPersist)
+template class FIO_TVector<uint_1>;
+template class FIO_TVector<uint_2>;
+template class FIO_TVector<int_2>;
+template class FIO_TVector<int_4>;
+template class FIO_TVector<int_8>;
+template class FIO_TVector<uint_4>;
+template class FIO_TVector<uint_8>;
+template class FIO_TVector<r_8>;
+template class FIO_TVector<r_4>;
+template class FIO_TVector< complex<r_4> >;
+template class FIO_TVector< complex<r_8> >;
+#endif
Index: /trunk/SophyaLib/TArray/tvector.h
===================================================================
--- /trunk/SophyaLib/TArray/tvector.h	(revision 762)
+++ /trunk/SophyaLib/TArray/tvector.h	(revision 762)
@@ -0,0 +1,107 @@
+// This may look like C code, but it is really -*- C++ -*-
+//                         C.Magneville          05/99
+#ifndef TVector_SEEN
+#define TVector_SEEN
+
+#include "tmatrix.h"
+
+namespace SOPHYA {
+
+template <class T>
+class TVector : public TMatrix<T> {
+public:
+
+  // Creation / destruction 
+  TVector(uint_4 n=1);
+  TVector(uint_4 n, T* values,Bridge* br=NULL);
+  TVector(const TVector<T>& v);
+  TVector(const TVector<T>& v,bool share);
+  TVector(const TMatrix<T>& a);
+
+  // Gestion taille/Remplissage
+  inline void ReSize(uint_4 n) {TMatrix<T>::ReSize(n,1);} // Reallocation de place
+  inline void Realloc(uint_4 n,bool force=false) {TMatrix<T>::Realloc(n,1,force);}
+
+  // Informations pointeur/data
+  inline uint_4 NElts() const {return NRows();}
+  
+  // Acces aux elements
+  inline T& operator()(uint_4 n) {return (*this)[n];}
+  inline T const& operator()(uint_4 n) const {return (*this)[n];}
+  inline T& Element(uint_4 n) {return (*this)[n];}
+  inline T const& Element(uint_4 n) const {return (*this)[n];}
+
+  // Operateur d'affectation
+  inline TVector& operator = (const TVector& v)
+                  {TMatrix<T>::operator =(v); return *this;}
+  inline TVector& operator = (T x)
+                  {for(uint_4 i=0;i<NRows();i++) (*this)(i)=x; return *this;}
+
+
+};
+
+// produit scalaire, matrice*vecteur
+template <class T>
+inline T operator* (const TVector<T>& v1, const TVector<T>& v2)
+  {if(v1.NRows() != v2.NRows())
+     throw(SzMismatchError("TVector::operator*(TVector& v1,TVector v2) size mismatch"));
+   T *p = const_cast<T *>(v1.Data()), *pEnd = p+v1.NElts(),
+     *q = const_cast<T *>(v2.Data()), r = 0;
+   while (p<pEnd) r += *p++ * *q++;
+   return r;}
+
+template <class T>
+inline TVector<T> operator* (const TMatrix<T>& a, const TVector<T>& b)
+   {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
+typedef TVector<r_8> Vector;
+
+/////////////////////////////////////////////////////////////////////////
+// Classe pour la gestion de persistance
+template <class T>
+class FIO_TVector : public  PPersist  {
+public:
+  FIO_TVector();
+  FIO_TVector(string const & filename); 
+  FIO_TVector(const TVector<T> & obj);
+  FIO_TVector(TVector<T> * obj);
+  virtual ~FIO_TVector();
+  virtual AnyDataObj* DataObj();
+  virtual void        SetDataObj(AnyDataObj & o);
+  inline operator TVector<T>() { return(*dobj); }
+protected :
+  virtual void ReadSelf(PInPersist&);           
+  virtual void WriteSelf(POutPersist&) const;  
+  TVector<T> * dobj;
+  bool ownobj;
+};
+
+template <class T>
+inline POutPersist& operator << (POutPersist& os, TVector<T> & obj)
+{ FIO_TVector<T> fio(&obj);  fio.Write(os);  return(os); }
+template <class T>
+inline PInPersist& operator >> (PInPersist& is, TVector<T> & obj)
+{ FIO_TVector<T> fio(&obj);  fio.Read(is);  return(is); }
+
+} // Fin du namespace
+
+#endif
