Index: trunk/SophyaLib/NTools/tmatrix.cc
===================================================================
--- trunk/SophyaLib/NTools/tmatrix.cc	(revision 298)
+++ trunk/SophyaLib/NTools/tmatrix.cc	(revision 299)
@@ -1,3 +1,3 @@
-// $Id: tmatrix.cc,v 1.4 1999-05-05 05:51:54 ansari Exp $
+// $Id: tmatrix.cc,v 1.5 1999-05-17 17:16:43 ansari Exp $
 //                         C.Magneville          04/99
 #include "machdefs.h"
@@ -9,4 +9,5 @@
 #include "tmatrix.h"
 #include "objfio.h"
+#include "generalfit.h"
 
 using namespace PlanckDPC;
@@ -130,5 +131,5 @@
 {
 if(mNr!=b.mNr || mNc!=b.mNc)
-    throw(SzMismatchError("NDataBlock operator C=A+B size mismatch\n"));
+    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;
@@ -139,5 +140,5 @@
 {
 if(mNr!=b.mNr || mNc!=b.mNc)
-    throw(SzMismatchError("NDataBlock operator C=A-B size mismatch\n"));
+    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;
@@ -149,5 +150,5 @@
 {
 if(mNr==0 || mNc==0 || b.mNc==0 || mNc!=b.mNr)
-  throw(SzMismatchError("NDataBlock operator C=A*B size mismatch\n"));
+  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;
@@ -160,4 +161,74 @@
 }
 return r;
+}
+
+#include "matrix.h"
+////////////////////////////////////////////////////////////////
+//**** Pour inversion
+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 Matrix"));
+return b;
+}
+
+double TMatrix<r_8>::GausPiv(TMatrix<r_8>& a, TMatrix<r_8>& b)
+// Pivot de Gauss
+{
+Matrix A(a);
+Matrix B(b);
+return Matrix::GausPiv(A,B);
+}
+
+//////////////////////////////////////////////////////////
+//**** Residus des fits
+TMatrix<r_8> TMatrix<r_8>::FitResidus(GeneralFit& gfit)
+//	Retourne une classe contenant les residus du fit ``gfit''.
+//	On suppose que x=j (colonnes) et y=i (lignes) pour m(i,j).
+{
+if(NCols()<=0||NRows()<=0)
+  throw(SzMismatchError("TMatrix::FitResidus size mismatch\n"));
+GeneralFunction* f = gfit.GetFunction();
+if(f==NULL)
+  throw(NullPtrError("TMatrix::FitResidus GeneraFit==NULL\n"));
+int npar =  gfit.GetNPar();
+if(npar==0)
+  throw(SzMismatchError("TMatrix::FitResidus GeneraFit 0 parametre\n"));
+double* par = new double[npar];
+{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
+TMatrix<r_8> m(*this);
+for(int i=0;i<NRows();i++) for(int j=0;j<NCols();j++) {
+  double x[2] = {(double)j,(double)i}; // attention x=col y=row!!
+  m(i,j) -= f->Value(x,par);
+}
+delete [] par;
+return m;
+}
+
+TMatrix<r_8> TMatrix<r_8>::FitFunction(GeneralFit& gfit)
+//	Retourne une classe contenant la fonction du fit ``gfit''.
+//	On suppose que x=j (colonnes) et y=i (lignes) pour m(i,j).
+{
+if(NCols()<=0||NRows()<=0)
+  throw(SzMismatchError("TMatrix::FitFunction size mismatch\n"));
+GeneralFunction* f = gfit.GetFunction();
+if(f==NULL)
+  throw(NullPtrError("TMatrix::FitFunction GeneraFit==NULL\n"));
+int npar =  gfit.GetNPar();
+if(npar==0)
+  throw(SzMismatchError("TMatrix::FitFunction GeneraFit 0 parametre\n"));
+double* par = new double[npar];
+{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
+TMatrix<r_8> m(*this);
+for(int i=0;i<NRows();i++) for(int j=0;j<NCols();j++) {
+  double x[2] = {(double)j,(double)i};
+  m(i,j) = f->Value(x,par);
+}
+delete [] par;
+return m;
 }
 
@@ -235,5 +306,4 @@
 os.Put(itab,3);
 // On ecrit le NDataBlock
-//cmv   encore des problemes avec les consteries
 FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
 fio_nd.Write(os);
Index: trunk/SophyaLib/NTools/tmatrix.h
===================================================================
--- trunk/SophyaLib/NTools/tmatrix.h	(revision 298)
+++ trunk/SophyaLib/NTools/tmatrix.h	(revision 299)
@@ -10,4 +10,5 @@
 #include "anydataobj.h"
 #include "ndatablock.h"
+class GeneralFit;
 
 namespace PlanckDPC {
@@ -93,4 +94,13 @@
   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);
+
+  // Residus et fonction fittees.
+  TMatrix<T> FitResidus(GeneralFit& gfit);
+  TMatrix<T> FitFunction(GeneralFit& gfit);
+
 protected:
   // partage les donnees si "a" temporaire, clone sinon.
Index: trunk/SophyaLib/NTools/tvector.cc
===================================================================
--- trunk/SophyaLib/NTools/tvector.cc	(revision 299)
+++ trunk/SophyaLib/NTools/tvector.cc	(revision 299)
@@ -0,0 +1,128 @@
+// $Id: tvector.cc,v 1.1 1999-05-17 17:16:44 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"
+#include "generalfit.h"
+
+using namespace PlanckDPC;
+
+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"));
+}
+
+//////////////////////////////////////////////////////////
+//**** Residus des fits
+TVector<r_8> TVector<r_8>::FitResidus(GeneralFit& gfit)
+//	Retourne une classe contenant les residus du fit ``gfit''.
+//	On suppose que x=i.
+{
+if(NElts()<=0)
+  throw(SzMismatchError("TVector::FitResidus size mismatch\n"));
+GeneralFunction* f = gfit.GetFunction();
+if(f==NULL)
+  throw(NullPtrError("TVector::FitResidus GeneraFit==NULL\n"));
+int npar =  gfit.GetNPar();
+if(npar==0)
+  throw(SzMismatchError("TVector::FitResidus GeneraFit 0 parametre\n"));
+double* par = new double[npar];
+{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
+TVector<r_8> v(*this);
+for(int i=0;i<NElts();i++) {
+  double x = (double)i;
+  v(i) -= f->Value(&x,par);
+}
+delete [] par;
+return v;
+}
+
+TVector<r_8> TVector<r_8>::FitFunction(GeneralFit& gfit)
+//	Retourne une classe contenant la fonction du fit ``gfit''.
+//	On suppose que x=i.
+{
+if(NElts()<=0)
+  throw(SzMismatchError("TVector::FitFunction size mismatch\n"));
+GeneralFunction* f = gfit.GetFunction();
+if(f==NULL)
+  throw(NullPtrError("TVector::FitFunction GeneraFit==NULL\n"));
+int npar =  gfit.GetNPar();
+if(npar==0)
+  throw(SzMismatchError("TVector::FitFunction GeneraFit 0 parametre\n"));
+double* par = new double[npar];
+{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
+TVector<r_8> v(*this);
+for(int i=0;i<NElts();i++) {
+  double x = (double)i;
+  v(i) = f->Value(&x,par);
+}
+delete [] par;
+return v;
+}
+
+///////////////////////////////////////////////////////////////
+#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<float> > 
+#pragma define_template TVector< complex<double> > 
+#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<float> >;
+template class TVector< complex<double> >;
+#endif
Index: trunk/SophyaLib/NTools/tvector.h
===================================================================
--- trunk/SophyaLib/NTools/tvector.h	(revision 299)
+++ trunk/SophyaLib/NTools/tvector.h	(revision 299)
@@ -0,0 +1,76 @@
+// This may look like C code, but it is really -*- C++ -*-
+//                         C.Magneville          05/99
+#ifndef TVECTOR_SEEN
+#define TVECTOR_SEEN
+
+#include "tmatrix.h"
+class GeneralFit;
+
+namespace PlanckDPC {
+
+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
+
+  // 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];}
+
+  // 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;}
+
+  // Residus et fonction fittees.
+  TVector<T> FitResidus(GeneralFit& gfit);
+  TVector<T> FitFunction(GeneralFit& gfit);
+
+};
+
+// 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);
+}
+
+} // Fin du namespace
+
+#endif
