// 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 {



////////////////////////////////////////////////////////////////
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);
}

inline r_4 LinSolve(const TMatrix<r_4>& a, const TVector<r_4>& b, TVector<r_4>& c)
{ 
if(a.NCols() != b.NRows() || a.NCols() != a.NRows())
  throw(SzMismatchError("LinSolve(TMatrix<r_4>,TVector<r_4>) size mismatch"));
c = b;
TMatrix<r_4> a1(a);
return SimpleMatrixOperation<r_4>::GausPiv(a1,c);
}

// Inverse d'une matrice 
inline TMatrix<r_8> Inverse(TMatrix<r_8> const & A)
{
  return SimpleMatrixOperation<r_8>::Inverse(A);
}

inline TMatrix<r_4> Inverse(TMatrix<r_4> const & A)
{
  return SimpleMatrixOperation<r_4>::Inverse(A);
}


//--------------------------------------
//        Linear fitting 
//--------------------------------------

class LinFitter {
public :
		LinFitter();
  virtual	~LinFitter();

  double	LinFit(const Vector& x, const Vector& y, int nf,
		       double (*f)(int, double), Vector& c);
// fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1;

  double	LinFit(const Matrix& fx, const Vector& y, Vector& c);
// fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1,
// la matrice fx contient les valeurs des f:
// fx(i,j) = f(i, x(j)).
		      
  double	LinFit(const Vector& x, const Vector& y, const Vector& errY2, int nf,
		       double (*f)(int, double), Vector& c, Vector& errC);
// fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1,
// errY2 contient les carres des erreurs sur les Y.
// au retour, errC contient les erreurs sur les coefs.

  double	LinFit(const Matrix& fx, const Vector& y, const Vector& errY2,
		       Vector& c, Vector& errC);
// fit lineaire des y en tant que somme de c(i)f(i,x), 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.
};


} // Fin du namespace

#endif