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

////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////
//------------------------------------------------------------//
//      La classe de calcul simple sur les TMatrix            //
//------------------------------------------------------------//
////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////

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

/*! \ingroup TArray 
    \fn LinSolveInPlace(TMatrix<r_4>&, TVector<r_4>&)
    \brief  Solve A*C = B for C in place and return determinant
*/
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);
}

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

/*! \ingroup TArray 
    \fn LinSolveInPlace(TMatrix< complex<r_4> >&, TVector< complex<r_4> >&)
    \brief Solve A*X = B in place and return determinant */
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);
}

/*! \ingroup TArray 
    \fn LinSolveInPlace(TMatrix< complex<r_8> >&, TVector< complex<r_8> >&)
    \brief  Solve A*X = B in place and return determinant
*/
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
//------------------------------------------------------------

/*! \ingroup TArray 
    \fn LinSolve(const TMatrix<r_4>&, const TVector<r_4>&, TVector<r_4>&)
    \brief Solve A*C = B and return C and determinant
*/

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

/*! \ingroup TArray 
    \fn LinSolve(const TMatrix<r_8>&, const TVector<r_8>&, TVector<r_8>&)
    \brief Solve A*C = B and return C and determinant
*/
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);
}

/*! \ingroup TArray 
    \fn LinSolve(const TMatrix< complex<r_4> >&, const TVector< complex<r_4> >&, TVector< complex<r_4> >&)
    \brief Solve A*C = B and return C and determinant
*/
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);
}

/*! \ingroup TArray 
    \fn LinSolve(const TMatrix< complex<r_8> >&, const TVector< complex<r_8> >&, TVector< complex<r_8> >&)
    \brief  Solve A*C = B and return C and determinant
*/
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 {

/*! \ingroup TArray 
    \fn Inverse(TMatrix<r_4> const &)
    \brief To inverse a TMatrix
*/
inline TMatrix<r_4> Inverse(TMatrix<r_4> const & A)
  {return SimpleMatrixOperation<r_4>::Inverse(A);}
/*! \ingroup TArray 
    \fn Inverse(TMatrix<r_8> const &)
    \brief To inverse a TMatrix
*/
inline TMatrix<r_8> Inverse(TMatrix<r_8> const & A)
  {return SimpleMatrixOperation<r_8>::Inverse(A);}
/*! \ingroup TArray 
    \fn Inverse(TMatrix< complex<r_4> > const &)
    \brief To inverse a TMatrix (complex numbers)
*/
inline TMatrix< complex<r_4> > Inverse(TMatrix< complex<r_4> > const & A)
  {return SimpleMatrixOperation< complex<r_4> >::Inverse(A);}
/*! \ingroup TArray 
    \fn Inverse(TMatrix< complex<r_8> > const &)
    \brief To inverse a TMatrix (complex numbers)
*/
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
template <class T>
class LinFitter {

public :

         LinFitter();
virtual	~LinFitter();

//! Linear fitting
r_8 LinFit(const TVector<T>& x, const TVector<T>& y,
	   uint_4 nf, T (*f)(uint_4,T), TVector<T>& c);

//! Linear fitting
r_8 LinFit(const TMatrix<T>& fx, const TVector<T>& y, TVector<T>& c);
		      
//! Linear fitting with errors
r_8 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);

//! Linear fitting with errors
r_8 LinFit(const TMatrix<T>& fx, const TVector<T>& y,
           const TVector<T>& errY2, TVector<T>& c, TVector<T>& errC);
};

} // Fin du namespace

#endif