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

////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////
//------------------------------------------------------------//
//                     Classe TMatrixRC                       //
//------------------------------------------------------------//
////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////

namespace SOPHYA {

/*!
  \class SimpleMatrixOperation
  \ingroup TArray
  Class for simple operation on TMatrix
  \sa TMatrix TArray
  */

//! Class for simple operation on TMatrix
template <class T>
class SimpleMatrixOperation {
public:
  static TMatrix<T> Inverse(TMatrix<T> const & A);
  static T GausPiv(TMatrix<T>& A, TMatrix<T>& B);
};

} // Fin du namespace

////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////
//------------------------------------------------------------//
//              Resolution de systemes lineaires              //
//------------------------------------------------------------//
////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////

namespace SOPHYA {

//------------------------------------------------------------
// Resolution du systeme A*C = B
//------------------------------------------------------------

//! Solve A*C = B for C in place and return determinant
/*! \ingroup TArray \fn LinSolveInPlace */
inline r_4 LinSolveInPlace(TMatrix<r_4>& a, TVector<r_4>& b)
{
if(a.NCols() != b.NRows() || a.NCols() != a.NRows())
  throw(SzMismatchError("LinSolveInPlace(TMatrix<r_4>,TVector<r_4>) size mismatch"));
return SimpleMatrixOperation<r_4>::GausPiv(a,b);
}

//! Solve A*X = B in place and return determinant
/*! \ingroup TArray \fn LinSolveInPlace */
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);
}

//! Solve A*X = B in place and return determinant
/*! \ingroup TArray \fn LinSolveInPlace */
inline complex<r_4> LinSolveInPlace(TMatrix< complex<r_4> >& a, TVector< complex<r_4> >& b)
{
if(a.NCols() != b.NRows() || a.NCols() != a.NRows())
  throw(SzMismatchError("LinSolveInPlace(TMatrix< complex<r_4> >,TVector< complex<r_4> >) size mismatch"));
return SimpleMatrixOperation< complex<r_4> >::GausPiv(a,b);
}

//! Solve A*X = B in place and return determinant
/*! \ingroup TArray \fn LinSolveInPlace */
inline complex<r_8> LinSolveInPlace(TMatrix< complex<r_8> >& a, TVector< complex<r_8> >& b)
{
if(a.NCols() != b.NRows() || a.NCols() != a.NRows())
  throw(SzMismatchError("LinSolveInPlace(TMatrix< complex<r_8> >,TVector< complex<r_8> >) size mismatch"));
return SimpleMatrixOperation< complex<r_8> >::GausPiv(a,b);
}

//------------------------------------------------------------
// Resolution du systeme A*C = B, avec C retourne dans B
//------------------------------------------------------------

//! Solve A*C = B and return C and determinant
/*! \ingroup TArray \fn LinSolve */
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);
}

//! Solve A*C = B and return C and determinant
/*! \ingroup TArray \fn LinSolve */
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);
}

//! Solve A*C = B and return C and determinant
/*! \ingroup TArray \fn LinSolve */
inline complex<r_4> LinSolve(const TMatrix< complex<r_4> >& a, const TVector< complex<r_4> >& b, TVector< complex<r_4> >& c) {
  if(a.NCols()!=b.NRows() || a.NCols()!=a.NRows())
    throw(SzMismatchError("LinSolve(TMatrix< complex<r_4> >,TVector< complex<r_4> >) size mismatch"));
    c = b; TMatrix< complex<r_4> > a1(a);
  return SimpleMatrixOperation< complex<r_4> >::GausPiv(a1,c);
}

//! Solve A*C = B and return C and determinant
/*! \ingroup TArray \fn LinSolve */
inline complex<r_8> LinSolve(const TMatrix< complex<r_8> >& a, const TVector< complex<r_8> >& b, TVector< complex<r_8> >& c) {
  if(a.NCols()!=b.NRows() || a.NCols()!=a.NRows())
    throw(SzMismatchError("LinSolve(TMatrix< complex<r_8> >,TVector< complex<r_8> >) size mismatch"));
    c = b; TMatrix< complex<r_8> > a1(a);
  return SimpleMatrixOperation< complex<r_8> >::GausPiv(a1,c);
}

} // Fin du namespace

////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////
//------------------------------------------------------------//
//                   Inverse d'une matrice                    //
//------------------------------------------------------------//
////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////

namespace SOPHYA {

//! To inverse a TMatrix
/*! \ingroup TArray \fn Inverse */
inline TMatrix<r_4> Inverse(TMatrix<r_4> const & A)
  {return SimpleMatrixOperation<r_4>::Inverse(A);}
//! To inverse a TMatrix
/*! \ingroup TArray \fn Inverse */
inline TMatrix<r_8> Inverse(TMatrix<r_8> const & A)
  {return SimpleMatrixOperation<r_8>::Inverse(A);}
//! To inverse a TMatrix
/*! \ingroup TArray \fn Inverse */
inline TMatrix< complex<r_4> > Inverse(TMatrix< complex<r_4> > const & A)
  {return SimpleMatrixOperation< complex<r_4> >::Inverse(A);}
//! To inverse a TMatrix
/*! \ingroup TArray \fn Inverse */
inline TMatrix< complex<r_8> > Inverse(TMatrix< complex<r_8> > const & A)
  {return SimpleMatrixOperation< complex<r_8> >::Inverse(A);}

} // Fin du namespace


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

namespace SOPHYA {

/*!
  \class LinFitter
  \ingroup TArray
  Class for linear fitting
  \sa TMatrix TArray
  */

//!  Class for 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