Index: trunk/SophyaLib/TArray/sopemtx.cc
===================================================================
--- trunk/SophyaLib/TArray/sopemtx.cc	(revision 813)
+++ trunk/SophyaLib/TArray/sopemtx.cc	(revision 926)
@@ -14,10 +14,18 @@
 // -------------------------------------------------------------
 ////////////////////////////////////////////////////////////////
-/////////////////////////////////////////////////////////////////////////
-// Classe de lignes/colonnes de matrices
-enum TRCKind {TmatrixRow=0, TmatrixCol=1, TmatrixDiag=2};
+
+//! Class of line, column or diagonal of a TMatrix
+/*!
+  A TMatrixRC represents a line, a column or the diagonal of a TMatrix
+ */
 template <class T>
 class TMatrixRC {
 public:
+  //! Define type of TMatrixRC
+  enum TRCKind {
+    TmatrixRow=0,  //!< TMatrixRC ligne
+    TmatrixCol=1,  //!< TMatrixRC column
+    TmatrixDiag=2  //!< TMatrixRC diagonal
+  };
   TMatrixRC();
   TMatrixRC(TMatrix<T>& m, TRCKind kind, uint_4 index=0);
@@ -39,4 +47,5 @@
   static T* Org(const TMatrix<T>&, TRCKind rckind, uint_4 ind=0);
 
+  //! Return the kind of TMatrix (line,column,diagonal)
   TRCKind Kind() const { return kind; }
   uint_4 NElts() const;
@@ -67,28 +76,39 @@
   static void Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2);
 
+  //! Define Absolute value for uint_1
   inline static double Abs_Value(uint_1 v) {return (double) v;}
+  //! Define Absolute value for uint_2
   inline static double Abs_Value(uint_2 v) {return (double) v;}
+  //! Define Absolute value for int_2
   inline static double Abs_Value(int_2 v)  {return (v>0)? (double) v: (double) -v;}
+  //! Define Absolute value for int_4
   inline static double Abs_Value(int_4 v)  {return (v>0)? (double) v: (double) -v;}
+  //! Define Absolute value for int_8
   inline static double Abs_Value(int_8 v)  {return (v>0)? (double) v: (double) -v;}
+  //! Define Absolute value for uint_4
   inline static double Abs_Value(uint_4 v) {return (double) v;}
+  //! Define Absolute value for uint_8
   inline static double Abs_Value(uint_8 v) {return (double) v;}
+  //! Define Absolute value for r_4
   inline static double Abs_Value(r_4 v)    {return (double) fabsf(v);}
+  //! Define Absolute value for r_8
   inline static double Abs_Value(r_8 v)    {return fabs(v);}
-  inline static double Abs_Value(complex<float> v)
+  //! Define Absolute value for complex r_4
+  inline static double Abs_Value(complex<r_4> v)
                 {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
-  inline static double Abs_Value(complex<double> v)
+  //! Define Absolute value for complex r_8
+  inline static double Abs_Value(complex<r_8> 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;
+  TMatrix<T>* matrix;  //!< pointer to the TMatrix
+  T*          data;    //!< pointer to the beginnig of interesting datas
+  int_4       index;   //!< index of the line/column
+  uint_4      step;    //!< step of the line/column
+  TRCKind     kind;    //!< type: line, column or diagonal
 };
 
 
-
+//! Multiply two TMatrixRC
 template <class T>
 inline T operator * (const TMatrixRC<T>& a, const TMatrixRC<T>& b)
@@ -103,5 +123,5 @@
   }
 
-
+//! Get the step in datas for a TMatrix for type rckind (line/col/diag)
 template <class T>
 inline uint_4 TMatrixRC<T>::Step(const TMatrix<T>& m, TRCKind rckind)
@@ -111,4 +131,6 @@
     return 0; }
 
+/*! Get the origin of datas for a TMatrix for type rckind and
+  number index (line/col/diag). ex: origine for line "index". */
 template <class T>
 inline T* TMatrixRC<T>::Org(const TMatrix<T>& m, TRCKind rckind, uint_4 index)
@@ -118,4 +140,5 @@
     return NULL; }
 
+//! return number of elements for a TMatrixRC
 template <class T> inline uint_4 TMatrixRC<T>::NElts() const
   { if (!matrix) return 0;
@@ -125,18 +148,27 @@
     return 0; }
 
+//! access of element \b i
 template <class T>
 inline T& TMatrixRC<T>::operator()(uint_4 i) {return data[i*step];}
+//! access of element \b i
 template <class T>
 inline T  TMatrixRC<T>::operator()(uint_4 i) const {return data[i*step];}
 
 ////////////////////////////////////////////////////////////////
-// Typedef pour simplifier et compatibilite Peida
+//! Typedef to simplifier TMatrixRC<r_8>
 typedef TMatrixRC<r_8> MatrixRC;
 
 
+//! Default constructor
 template <class T> TMatrixRC<T>::TMatrixRC()
 : matrix(NULL), data(NULL), index(0), step(0)
 {}
 
+//! Constructor
+/*!
+  \param m : matrix
+  \param rckind : select line, column or diagonal
+  \param ind : number of the line or column
+*/
 template <class T> TMatrixRC<T>::TMatrixRC(TMatrix<T>& m,TRCKind rckind,uint_4 ind)
 : matrix(&m), data(Org(m,rckind,ind)),
@@ -150,4 +182,5 @@
 // Acces aux rangees et colonnes de matrices
 
+//! Return TMatrixRC for line \b r of matrix \b m
 template <class T> 
 TMatrixRC<T> TMatrixRC<T>::Row(TMatrix<T> & m, uint_4 r) 
@@ -157,4 +190,5 @@
 }
 
+//! Return TMatrixRC for column \b r of matrix \b m
 template <class T> 
 TMatrixRC<T> TMatrixRC<T>::Col(TMatrix<T> & m, uint_4 c) 
@@ -164,4 +198,5 @@
 }
 
+//! Return TMatrixRC for diagonal of matrix \b m
 template <class T> 
 TMatrixRC<T> TMatrixRC<T>::Diag(TMatrix<T> & m) 
@@ -197,5 +232,5 @@
 // }
 
- 
+//! Set column \b c for this TMatrixRC
 template <class T> int_4 TMatrixRC<T>::SetCol(int_4 c)
 {
@@ -209,4 +244,5 @@
 }
 
+//! Set line \b r for this TMatrixRC
 template <class T> int_4 TMatrixRC<T>::SetRow(int_4 r)
 {
@@ -220,4 +256,5 @@
 }
 
+//! Set line diaginal for this TMatrixRC
 template <class T> int_4 TMatrixRC<T>::SetDiag()
 {
@@ -232,5 +269,5 @@
 }
 
-
+//! Operator =
 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator = (const TMatrixRC<T>& rc)
 {
@@ -271,5 +308,5 @@
 // }
 
-
+//! Operator to multiply by constant \b x
 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator *= (T x)
 {
@@ -278,4 +315,5 @@
 }
 
+//! Operator to divide by constant \b x
 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator /= (T x)
 {
@@ -298,4 +336,9 @@
 // }
 
+//! Linear combination
+/*!
+  Do : \f$ MRC(i) = MRC(i)*a + rc(i)*b \f$
+  \return *this
+ */
 template <class T>
 TMatrixRC<T>& TMatrixRC<T>::LinComb(T a, T b, const TMatrixRC<T>& rc, uint_4 first)
@@ -309,4 +352,8 @@
 }
 
+//! Linear combination
+/*!
+  Do : \f$ MRC(i) = MRC(i) + rc(i)*b \f$
+ */
 template <class T>
 TMatrixRC<T>& TMatrixRC<T>::LinComb(T b, const TMatrixRC<T>& rc, uint_4 first)
@@ -320,4 +367,5 @@
 }
 
+//! Find maximum absolute value in TMatrixRC, search begin at \b first
 template <class T> uint_4 TMatrixRC<T>::IMaxAbs(uint_4 first)
 {
@@ -333,4 +381,5 @@
 }
 
+//! Print on stream \b os
 template <class T>
 void TMatrixRC<T>::Print(ostream & os) const
@@ -346,4 +395,5 @@
 }
 
+//! Swap two TMatrixRC of the same kind
 template <class T>
 void TMatrixRC<T>::Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2)
@@ -359,7 +409,10 @@
 
 
-
-
 ////////////////////////////////////////////////////////////////
+// -------------------------------------------------------------
+//   La classe de calcul simple sur les TMatrix
+// -------------------------------------------------------------
+////////////////////////////////////////////////////////////////
+
 //**** Pour inversion
 #ifndef M_LN2
@@ -373,4 +426,9 @@
 #endif
 
+//! Gaussian pivoting
+/*!
+  Diagonalize matrix \b a, doing the same opreations on matrix \b b
+  \return determinat of \b a
+ */  
 template <class T>
 T SimpleMatrixOperation<T>::GausPiv(TMatrix<T>& a, TMatrix<T>& b)
@@ -394,10 +452,10 @@
 double nrm = sqrt(vmin*vmax);
 if(nrm > 1.e5 || nrm < 1.e-5) {
-  a /= nrm;
-  b /= nrm;
+  a /= (T) nrm;
+  b /= (T) nrm;
   //cout << "normalisation matrice " << nrm << endl;
 } else nrm=1;
 
-double det = 1.0;
+T det = 1;
 if(nrm != 1) {
   double ld = a.NRows() * log(nrm);
@@ -405,10 +463,10 @@
    // cerr << "TMatrix warning, overflow for det" << endl;
   } else {
-    det = exp(ld);
-  }
-}
-
-TMatrixRC<T> pivRowa(a,TmatrixRow);
-TMatrixRC<T> pivRowb(b,TmatrixRow);
+    det = (T) exp(ld);
+  }
+}
+
+TMatrixRC<T> pivRowa(a,TMatrixRC<T>::TmatrixRow);
+TMatrixRC<T> pivRowb(b,TMatrixRC<T>::TmatrixRow);
 
 for(uint_4 k=0; k<n-1; k++) {
@@ -422,11 +480,11 @@
     TMatrixRC<T>::Swap(bIPiv,bK);
   }
-  double pivot = a(k,k);
-  if (fabs(pivot) < 1.e-50) return 0.0;
+  T pivot = a(k,k);
+  if( TMatrixRC<T>::Abs_Value(pivot) < 1.e-50 ) return (T) 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;
+    T 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);
@@ -437,10 +495,10 @@
 // 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;
+  T pivot = a(kk,kk);
+  if( TMatrixRC<T>::Abs_Value(pivot) <= 1.e-50 ) return (T) 0;
   pivRowa.SetRow(kk); // to avoid constructors
   pivRowb.SetRow(kk);
   for(uint_4 jj=0; jj<kk; jj++) {
-    double r = -a(jj,kk)/pivot;
+    T r = -a(jj,kk)/pivot;
     TMatrixRC<T>::Row(a, jj).LinComb(r, pivRowa);
     TMatrixRC<T>::Row(b, jj).LinComb(r, pivRowb);
@@ -449,5 +507,5 @@
 
 for(uint_4 l=0; l<n; l++) {
-  if (fabs((double)a(l,l)) <= 1.e-50) return 0.0;
+  if( TMatrixRC<T>::Abs_Value(a(l,l)) <= 1.e-50 ) return (T) 0;
   TMatrixRC<T>::Row(b, l) /= a(l,l);
 }
@@ -456,15 +514,21 @@
 }
 
+//! Return the inverse matrix of \b A
 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 = IdentityMatrix(1.);
-if(fabs((double)GausPiv(a,b)) < 1.e-50)
-  throw(MathExc("TMatrix Inverse() Singular OMatrix"));
+if( TMatrixRC<T>::Abs_Value(GausPiv(a,b)) < 1.e-50)
+  throw(MathExc("TMatrix Inverse() Singular Matrix"));
 return b;
 }
 
+
+////////////////////////////////////////////////////////////////
+// -------------------------------------------------------------
+//   La classe fit lineaire
+// -------------------------------------------------------------
+////////////////////////////////////////////////////////////////
 
 LinFitter::LinFitter()
@@ -626,7 +690,11 @@
 #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<int_4>
 #pragma define_template SimpleMatrixOperation<r_4>
 #pragma define_template SimpleMatrixOperation<r_8>
+#pragma define_template SimpleMatrixOperation< complex<r_4> >
+#pragma define_template SimpleMatrixOperation< complex<r_8> >
 #endif
 
@@ -636,6 +704,10 @@
 template class TMatrixRC<r_4>;
 template class TMatrixRC<r_8>;
+template class TMatrixRC< complex<r_4> >;
+template class TMatrixRC< complex<r_8> >;
 template class SimpleMatrixOperation<int_4>;
 template class SimpleMatrixOperation<r_4>;
 template class SimpleMatrixOperation<r_8>;
+template class SimpleMatrixOperation< complex<r_4> >;
+template class SimpleMatrixOperation< complex<r_8> >;
 #endif
