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

/////////////////////////////////////////////////////////////////////////
// Classe de lignes/colonnes de matrices
enum TRCKind {TmatrixRow=0, TmatrixCol=1, TmatrixDiag=2};
template <class T>
class TMatrixRC {
public:
  TMatrixRC();
  TMatrixRC(TMatrix<T>& m, TRCKind kind, uint_4 index=0);
  virtual ~TMatrixRC() {}

  // Acces aux rangees et colonnes de matrices
  static TMatrixRC<T> Row(TMatrix<T> & m, uint_4 r);
  static TMatrixRC<T> Col(TMatrix<T> & m, uint_4 c);
  static TMatrixRC<T> Diag(TMatrix<T> & m);

  int_4 Next();
  int_4 Prev();
  int_4 SetCol(int_4 c);
  int_4 SetRow(int_4 r);
  int_4 SetDiag();

  static uint_4 Step(const TMatrix<T>& m, TRCKind rckind);
  static T* Org(const TMatrix<T>&, TRCKind rckind, uint_4 ind=0);

  TRCKind Kind() const { return kind; }
  uint_4 NElts() const;
  T& operator()(uint_4 i);
  T  operator()(uint_4 i) const;

  TMatrixRC<T>& operator = (const TMatrixRC<T>& rc);
  TVector<T> GetVect() const;

  TMatrixRC<T>& operator += (const TMatrixRC<T>& rc);
  TMatrixRC<T>& operator -= (const TMatrixRC<T>& rc);

  TMatrixRC<T>& operator *= (T x);
  TMatrixRC<T>& operator /= (T x);
  TMatrixRC<T>& operator -= (T x);
  TMatrixRC<T>& operator += (T x);

  TMatrixRC<T>& LinComb(T a, T b, const TMatrixRC& rc, uint_4 first=0);
  TMatrixRC<T>& LinComb(T b, const TMatrixRC<T>& rc, uint_4 first=0);

  uint_4 IMaxAbs(uint_4 first=0);

  static void Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2);

  inline static double Abs_Value(uint_1 v) {return (double) v;}
  inline static double Abs_Value(uint_2 v) {return (double) v;}
  inline static double Abs_Value(int_2 v)  {return (v>0)? (double) v: (double) -v;}
  inline static double Abs_Value(int_4 v)  {return (v>0)? (double) v: (double) -v;}
  inline static double Abs_Value(int_8 v)  {return (v>0)? (double) v: (double) -v;}
  inline static double Abs_Value(uint_4 v) {return (double) v;}
  inline static double Abs_Value(uint_8 v) {return (double) v;}
  inline static double Abs_Value(r_4 v)    {return (double) fabsf(v);}
  inline static double Abs_Value(r_8 v)    {return fabs(v);}
  inline static double Abs_Value(complex<float> v)
                {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
  inline static double Abs_Value(complex<double> v)
                {return sqrt(v.real()*v.real()+v.imag()*v.imag());}

protected:
  TMatrix<T>* matrix;
  T*          data;
  int_4       index;
  uint_4      step;
  TRCKind     kind;
};


template <class T>
inline T operator * (const TMatrixRC<T>& a, const TMatrixRC<T>& b)
  {
  if ( a.NElts() != b.NElts() )
    throw(SzMismatchError("TMatrixRC::operator * size mismatch\n"));
  if ( a.Kind() != b.Kind() )
    throw(SzMismatchError("TMatrixRC::operator * type mismatch\n"));
  T sum = 0;
  for(uint_4 i=0; i<a.NElts(); i++) sum += a(i)*b(i);
  return sum;
  }

template <class T>
inline uint_4 TMatrixRC<T>::Step(const TMatrix<T>& m, TRCKind rckind)
  { switch (rckind) { case TmatrixRow  : return 1;
                      case TmatrixCol  : return m.NCols();
                      case TmatrixDiag : return m.NCols()+1; }
    return 0; }

template <class T>
inline T* TMatrixRC<T>::Org(const TMatrix<T>& m, TRCKind rckind, uint_4 index)
  { switch (rckind) { case TmatrixRow  : return const_cast<T *>(m.Data()) + index * m.NCols();
                      case TmatrixCol  : return const_cast<T *>(m.Data()) + index;
                      case TmatrixDiag : return const_cast<T *>(m.Data()); }
    return NULL; }

template <class T> inline uint_4 TMatrixRC<T>::NElts() const
  { if (!matrix) return 0;
    switch (kind) { case TmatrixRow  : return matrix->NCols();
                    case TmatrixCol  : return matrix->NRows();
                    case TmatrixDiag : return matrix->NCols(); }
    return 0; }

template <class T>
inline T& TMatrixRC<T>::operator()(uint_4 i) {return data[i*step];}
template <class T>
inline T  TMatrixRC<T>::operator()(uint_4 i) const {return data[i*step];}

////////////////////////////////////////////////////////////////
// Typedef pour simplifier et compatibilite Peida
typedef TMatrixRC<r_8> MatrixRC;

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


} // Fin du namespace

#endif