Index: trunk/SophyaLib/TArray/sopemtx.cc
===================================================================
--- trunk/SophyaLib/TArray/sopemtx.cc	(revision 935)
+++ trunk/SophyaLib/TArray/sopemtx.cc	(revision 939)
@@ -10,8 +10,15 @@
 
 ////////////////////////////////////////////////////////////////
-// -------------------------------------------------------------
-//   La classe de gestion des lignes et colonnes d'une matrice
-// -------------------------------------------------------------
+// ---------------------------------------------------------- //
+//  La classe de gestion des lignes et colonnes d'une matrice //
+// ---------------------------------------------------------- //
 ////////////////////////////////////////////////////////////////
+
+/*!
+  \class SOPHYA::TMatrixRC
+  \ingroup TArray
+  Class for representing rows, columns and diagonal of a matrix.
+  \sa TMatrix TArray
+  */
 
 //! Class of line, column or diagonal of a TMatrix
@@ -109,5 +116,4 @@
 };
 
-
 //! Scalar product of two TMatrixRC
 /*!
@@ -423,7 +429,7 @@
 
 ////////////////////////////////////////////////////////////////
-// -------------------------------------------------------------
-//   La classe de calcul simple sur les TMatrix
-// -------------------------------------------------------------
+// ---------------------------------------------------------- //
+//      La classe de calcul simple sur les TMatrix            //
+// ---------------------------------------------------------- //
 ////////////////////////////////////////////////////////////////
 
@@ -542,133 +548,192 @@
 
 ////////////////////////////////////////////////////////////////
-// -------------------------------------------------------------
-//   La classe fit lineaire
-// -------------------------------------------------------------
+// ---------------------------------------------------------- //
+//               La classe fit lineaire                       //
+// ---------------------------------------------------------- //
 ////////////////////////////////////////////////////////////////
 
-LinFitter::LinFitter()
-{
-}
-
-LinFitter::~LinFitter()
-{
-}
-
-double LinFitter::LinFit(const Vector& x, const Vector& y, int nf, 
-                         double (*f)(int, double), Vector& c)
-{
-  int n = x.NElts();
-  if (n != y.NElts()) THROW(sizeMismatchErr);
-  
-  Matrix fx(nf, n);
-  for (int i=0; i<nf; i++)
-    for (int j=0; j<n; j++)
-      fx(i,j) = f(i,x(j));
-
-  return LinFit(fx,y,c);
-}
-
-
-
-double LinFitter::LinFit(const Matrix& fx, const Vector& y, Vector& c)
-{
-  int n = y.NElts();
-  if ( n != fx.NCol()) THROW(sizeMismatchErr);
-
-  int nf = fx.NRows();
-
-  Matrix a(nf,nf);
-
-  for (int j=0; j<nf; j++)
-    for (int k=j; k<nf; k++)
-      a(j,k) = a(k,j) = TMatrixRC<r_8>::Row(const_cast<Matrix &>(fx), j) 
-
-	* TMatrixRC<r_8>::Row(const_cast<Matrix &>(fx), k); /* $CHECK$ Reza 10/3/2000 */
-   
-  c = fx * y;
-
-  if (SimpleMatrixOperation<r_8>::GausPiv(a,c) == 0) THROW(singMatxErr); /* $CHECK$ Reza 10/3/2000 */
-
-  double xi2 = 0;
-
-  for (int k=0; k<n; k++) {
-    double x = 0;
-    for (int i=0; i<nf; i++)
-      x += c(i) * fx(i,k);
-    x -= y(k);
-    xi2 += x*x;
+//! Creator
+template <class T>
+LinFitter<T>::LinFitter()
+{
+}
+
+//! Destructor
+template <class T>
+LinFitter<T>::~LinFitter()
+{
+}
+
+// fit lineaire des y(k) en tant que somme de c(i)f(i,x(k)), i=0..nf-1;
+//! Linear fitting
+/*!
+  Linear fit of y(k) as the sum of \f$ c(i)f(i,x(k)), i=0..nf-1 \f$
+  \param x : vector of X values
+  \param y : vector of datas
+  \param nf: number of functions
+  \param f : f(i,x(k)), i=0..nf-1
+  \return c : vector of solutions
+  \return return chisquare
+ */
+template <class T>
+r_8 LinFitter<T>::LinFit(const TVector<T>& x, const TVector<T>& y,
+                       uint_4 nf, T (*f)(uint_4,T), TVector<T>& c)
+{
+uint_4 n = x.NElts();
+if (n != y.NElts()) THROW(sizeMismatchErr);
+ 
+TMatrix<T> fx(nf,n);
+
+for(uint_4 i=0; i<nf; i++)
+  for(uint_4 j=0; j<n; j++) fx(i,j) = f(i,x(j));
+
+return LinFit(fx,y,c);
+}
+
+// fit lineaire des y(k) en tant que somme de c(i)f(i,x(k)), i=0..nf-1,
+// la matrice fx contient les valeurs des f: fx(i,j) = f(i, x(j)).
+//! Linear fitting
+/*!
+  Linear fit of y(k) as the sum of \f$ c(i)f(i,x(k)), i=0..nf-1 \f$.
+  \param fx : matrix which contains \f$ fx(i,j) = f(i, x(j)) \f$.
+  \param y : vector of datas
+  \return c : vector of solutions
+  \return return chisquare
+ */
+template <class T>
+r_8 LinFitter<T>::LinFit(const TMatrix<T>& fx, const TVector<T>& y, TVector<T>& c)
+{
+uint_4 n = y.NElts();
+if (n != fx.NCol()) THROW(sizeMismatchErr);
+
+uint_4 nf = fx.NRows();
+
+TMatrix<T> a(nf,nf);
+
+for(uint_4 j=0; j<nf; j++)
+  for(uint_4 k=j; k<nf; k++)
+     a(j,k) = a(k,j) = TMatrixRC<T>::Row(const_cast<TMatrix<T> &>(fx), j) 
+	             * TMatrixRC<T>::Row(const_cast<TMatrix<T> &>(fx), k);
+                                /* $CHECK$ Reza 10/3/2000 */
+
+c = fx * y;
+
+if(SimpleMatrixOperation<T>::GausPiv(a,c)==(T) 0) THROW(singMatxErr);
+                   /* $CHECK$ Reza 10/3/2000 */
+
+r_8 xi2 = 0., ax;
+T x;
+for(uint_4 k=0; k<n; k++) {
+  x = (T) 0;
+  for(uint_4 i=0; i<nf; i++) x += c(i) * fx(i,k);
+  x -= y(k);
+  ax = TMatrixRC<T>::Abs_Value(x);
+  xi2 += ax*ax;
+}
+return xi2;
+}
+
+// fit lineaire des y(k) en tant que somme de c(i)f(i,x(k)), i=0..nf-1,
+// errY2 contient les carres des erreurs sur les Y.
+// au retour, errC contient les erreurs sur les coefs.
+//! Linear fitting with errors
+/*!
+  Linear fit with errors of y(k) as the sum of \f$ c(i)f(i,x(k)), i=0..nf-1 \f$.
+  \param x : vector of X values
+  \param y : vector of datas
+  \param errY2 : vector of errors square on Y
+  \param nf: number of functions
+  \param f : f(i,x(k)), i=0..nf-1
+  \return c : vector of solutions
+  \return errC : vector of errors on solutions C
+  \return return chisquare
+ */
+template <class T>
+r_8 LinFitter<T>::LinFit(const TVector<T>& x, const TVector<T>& y,
+             const TVector<T>& errY2, uint_4 nf, T (*f)(uint_4,T),
+                                  TVector<T>& c, TVector<T>& errC)
+{
+uint_4 n = x.NElts();
+if (n != y.NElts()) THROW(sizeMismatchErr);
+
+TMatrix<T> fx(nf, n);
+for(uint_4 i=0; i<nf; i++)
+  for(uint_4 j=0; j<n; j++)
+    fx(i,j) = f(i,x(j));
+
+return LinFit(fx,y,errY2,c,errC);
+}
+
+// fit lineaire des y(k) en tant que somme de c(i)f(i,x(k)), i=0..nf-1,
+// la matrice fx contient les valeurs des f:
+// fx(i,j) = f(i, x(j)).
+// errY2 contient les carres des erreurs sur les Y.
+// au retour, errC contient les erreurs sur les coefs.
+//! Linear fitting with errors
+/*!
+  \param fx : matrix which contains \f$ fx(i,j) = f(i, x(j)) \f$.
+  \param y : vector of datas
+  \param errY2 : vector of errors square on Y
+  \return c : vector of solutions
+  \return errC : vector of errors on solutions on C
+  \return return chisquare
+ */
+template <class T>
+r_8 LinFitter<T>::LinFit(const TMatrix<T>& fx, const TVector<T>& y,
+           const TVector<T>& errY2,TVector<T>& c, TVector<T>& errC)
+{
+uint_4 n = y.NElts();
+if( n != errY2.NElts()) THROW(sizeMismatchErr);
+if( n != fx.NCol()) THROW(sizeMismatchErr);
+
+uint_4 nf = fx.NRows();
+
+TMatrix<T> a(nf,nf);
+
+c.Realloc(nf);
+errC.Realloc(nf);
+
+for(uint_4 j=0; j<nf; j++)
+  for(uint_4 k=j; k<nf; k++) {
+    T x=0;
+    // Matrice a inverser
+    for(uint_4 l=0; l<n; l++) x += fx(j,l)*fx(k,l)/errY2(l);
+    a(j,k) = a(k,j) = x;
+  }
+ 
+TMatrix<T> d(nf,nf+1);
+for(uint_4 k=0; k<nf; k++) {
+  T x = (T) 0;
+  // Second membre 1ere colonne
+  for(uint_4 l=0; l<n; l++) x += fx(k,l)*y(l)/errY2(l);             
+  d(k,0) = x;
+  // Reste second membre = Id.
+  for(uint_4 m=1; m<=nf; m++) d(k,m) = (k==m)?1:0;
+}
+
+if(SimpleMatrixOperation<T>::GausPiv(a,d)==(T) 0) THROW(singMatxErr);
+                                            /* $CHECK$ Reza 10/3/2000 */
+
+for(uint_4 l=0; l<nf; l++) {
+  c(l) = d(l,0);        // Parametres = 1ere colonne
+  errC(l) = d(l,l+1);   // Erreurs = diag inverse.
+}
+
+r_8 xi2 = 0., ax;
+T x;
+for(uint_4 jj=0; jj<n; jj++) {
+  x = (T) 0;
+  for(uint_4 ii=0; ii<nf; ii++) x += c(ii) * fx(ii,jj);
+    x -= y(jj);
+    ax = TMatrixRC<T>::Abs_Value(x);
+    xi2 += ax*ax/TMatrixRC<T>::Abs_Value(errY2(jj));
   }
   return xi2;
 }
 
-
-
-double LinFitter::LinFit(const Vector& x, const Vector& y, const Vector& errY2, int nf,
-	                 double (*f)(int, double), Vector& c, Vector& errC)
-{
-  int n = x.NElts();
-  if (n != y.NElts()) THROW(sizeMismatchErr);
-  
-  Matrix fx(nf, n);
-  for (int i=0; i<nf; i++)
-    for (int j=0; j<n; j++)
-      fx(i,j) = f(i,x(j));
-
-  return LinFit(fx,y,errY2,c,errC);
-}
-
-
-double LinFitter::LinFit(const Matrix& fx, const Vector& y, const Vector& errY2,
-	   Vector& c, Vector& errC)
-{
-  int n = y.NElts();
-  if ( n != errY2.NElts()) THROW(sizeMismatchErr);
-  if ( n != fx.NCol()) THROW(sizeMismatchErr);
-
-  int nf = fx.NRows();
-
-  Matrix a(nf,nf);
-
-  c.Realloc(nf);
-  errC.Realloc(nf);
-
-  for (int j=0; j<nf; j++)
-    for (int k=j; k<nf; k++) {
-      double x=0;
-      for (int l=0; l<n; l++)
-	x += fx(j,l) * fx(k,l) / errY2(l);        // Matrice a inverser
-      a(j,k) = a(k,j) = x;
-    }
- 
-  Matrix d(nf,nf+1);
-  for (int k=0; k<nf; k++) {
-    double x=0;
-    for (int l=0; l<n; l++)
-      x += fx(k,l) * y(l) / errY2(l);             // Second membre 1ere colonne
-    d(k,0) = x;
-    for (int m=1; m<=nf; m++)
-      d(k,m) = (k==m) ? 1 : 0;                // Reste second membre = Id.
-  }
-
-  if (SimpleMatrixOperation<r_8>::GausPiv(a,d) == 0) THROW(singMatxErr); /* $CHECK$ Reza 10/3/2000 */
-
-  for (int l=0; l<nf; l++) {
-    c(l) = d(l,0);                            // Parametres = 1ere colonne
-    errC(l) = d(l,l+1);                        // Erreurs = diag inverse.
-  }
-
-  double xi2 = 0;
-
-  for (int jj=0; jj<n; jj++) {
-    double x = 0;
-    for (int ii=0; ii<nf; ii++)
-      x += c(ii) * fx(ii,jj);
-    x -= y(jj);
-    xi2 += x*x/errY2(jj);
-  }
-  return xi2;
-}
-
-
+///////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////
 void _ZZ_TestTMatrixRC_YY_(TMatrix<r_8> & m)
 {
@@ -712,4 +777,8 @@
 #pragma define_template SimpleMatrixOperation< complex<r_4> >
 #pragma define_template SimpleMatrixOperation< complex<r_8> >
+#pragma define_template LinFitter<r_4>
+#pragma define_template LinFitter<r_8>
+#pragma define_template LinFitter< complex<r_4> >
+#pragma define_template LinFitter< complex<r_8> >
 #endif
 
@@ -726,4 +795,8 @@
 template class SimpleMatrixOperation< complex<r_4> >;
 template class SimpleMatrixOperation< complex<r_8> >;
+template class LinFitter<r_4>;
+template class LinFitter<r_8>;
+template class LinFitter< complex<r_4> >;
+template class LinFitter< complex<r_8> >;
 #endif
 
