// This may look like C code, but it is really -*- C++ -*-
//                         C.Magneville          04/99
#ifndef TMatrix_SEEN
#define TMatrix_SEEN

#include "machdefs.h"
#include <stdio.h>
#include <iostream.h>
#include <complex>
#include "ppersist.h"
#include "anydataobj.h"
#include "ndatablock.h"

namespace PlanckDPC {

class GeneralFit;
template <class T> class TVector;
template <class T> class TMatrixRC;

template <class T>
class TMatrix : public AnyDataObj {
  friend class TMatrixRC<T>;
  friend class TVector<T>;
public:

  // Creation / destruction 
  TMatrix();
  TMatrix(uint_4 r,uint_4 c);
  TMatrix(uint_4 r,uint_4 c,T* values,Bridge* br=NULL);
  TMatrix(const TMatrix<T>& a);
  TMatrix(const TMatrix<T>& a,bool share);
  virtual ~TMatrix();

  // Temporaire?
  inline bool IsTemp(void) const {return mNDBlock.IsTemp();}
  inline void SetTemp(bool temp=false) const {mNDBlock.SetTemp(temp);}

  // Gestion taille/Remplissage
  inline void Clone(const TMatrix<T>& a)  // Clone: copie des donnees de "a"
              {mNDBlock.Clone(a.mNDBlock); mNr = a.mNr; mNc = a.mNc;}
  inline void Reset(T v=0) {mNDBlock.Reset(v);}
  inline void ReSize(uint_4 r,uint_4 c)  // Reallocation de place
       {if(r==0||c==0) throw(SzMismatchError("TMatrix::ReSize r ou c==0\n"));
        mNDBlock.ReSize(r*c); mNr = r; mNc = c;}
  inline void Realloc(uint_4 r,uint_4 c,bool force=false)
       {if(r==0||c==0) throw(SzMismatchError("TMatrix::Realloc r ou c==0\n"));
        mNDBlock.Realloc(r*c,force); mNr = r; mNc = c;}

  // Informations pointeur/data
  inline uint_4 NRows() const {return mNr;}
  inline uint_4 NCols() const {return mNc;}
  inline T const& operator()(uint_4 r,uint_4 c) const
                            {return *(mNDBlock.Begin()+r*mNc+c);}
  inline T&       operator()(uint_4 r,uint_4 c)
                            {return *(mNDBlock.Begin()+r*mNc+c);}
  inline T const& operator[](uint_4 ip) const
                            {return *(mNDBlock.Begin()+ip);}
  inline T&       operator[](uint_4 ip)
                            {return *(mNDBlock.Begin()+ip);}
  inline       T* Data()       {return mNDBlock.Begin();}
  inline const T* Data() const {return mNDBlock.Begin();}
  inline       NDataBlock<T>& DataBlock()       {return mNDBlock;}
  inline const NDataBlock<T>& DataBlock() const {return mNDBlock;}

  // Operations matricielles
  TMatrix<T> Transpose(void) const;

  // Operateur d'affectation
  // A = x (matrice diagonale x*Identite)
  inline TMatrix<T>& operator = (T x)
    {if(mNr!=mNc || mNr==0) throw(SzMismatchError("TMatrix::operator= mNc!=mNr ou ==0\n"));
     for(uint_4 r=0;r<mNr;r++) for(uint_4 c=0;c<mNc;c++) (*this)(r,c)=(r==c)?x:0;
     return *this;}
  // A = B : partage les donnees si "a" est temporaire, clone sinon.
  inline TMatrix<T>& operator = (const TMatrix<T>& a)
                {if(this == &a) return *this; CloneOrShare(a); return *this;}

  // Impression
  void Print(ostream& os,int lp=0,uint_4 i0=0,uint_4 ni=10,uint_4 j0=0,uint_4 nj=10) const;
  inline void Print(int lp=0,uint_4 i0=0,uint_4 ni=10,uint_4 j0=0,uint_4 nj=10) const
              {Print(cout,lp,i0,ni,j0,nj);}

  // Surcharge d'operateurs INPLACE: A (+=,-=,*=,/=) (T) x
  inline TMatrix<T>& operator += (T b) {mNDBlock += b; return *this;}
  inline TMatrix<T>& operator -= (T b) {mNDBlock -= b; return *this;}
  inline TMatrix<T>& operator *= (T b) {mNDBlock *= b; return *this;}
  inline TMatrix<T>& operator /= (T b) {mNDBlock /= b; return *this;}

  // Surcharge d'operateurs INPLACE: A (+=,-=,*=,/=) B
  inline TMatrix<T>& operator += (const TMatrix<T>& a)
              {if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
                  throw(SzMismatchError("TMatrix::operator+=A size mismatch"));
               mNDBlock += a.mNDBlock; return *this;}
  inline TMatrix<T>& operator -= (const TMatrix<T>& a)
              {if(mNr==0 || mNc==0 || mNr!=a.mNr || mNc!=a.mNc)
                  throw(SzMismatchError("TMatrix::operator-=A size mismatch"));
               mNDBlock -= a.mNDBlock; return *this;}
  TMatrix<T>& operator *= (const TMatrix<T>& a);

  // Pour surcharge d'operateurs C = A (+,-,*) B
  TMatrix<T> Add(const TMatrix<T>& b) const;
  TMatrix<T> Sub(const TMatrix<T>& b) const;
  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
            ,double xorg=0.,double yorg=0.,double dx=1.,double dy=1.);
  TMatrix<T> FitFunction(GeneralFit& gfit
            ,double xorg=0.,double yorg=0.,double dx=1.,double dy=1.);

  // Acces aux rangees et colonnes
  TMatrixRC<T> Row(uint_4 r) const;
  TMatrixRC<T> Col(uint_4 c) const;
  TMatrixRC<T> Diag() const;

protected:
  // partage les donnees si "a" temporaire, clone sinon.
  inline void CloneOrShare(const TMatrix<T>& a)
              {mNDBlock.CloneOrShare(a.mNDBlock); mNr=a.mNr; mNc=a.mNc;}
  // Share: partage les donnees de "a"
  inline void Share(const TMatrix<T>& a)
              {mNDBlock.Share(a.mNDBlock); mNr=a.mNr; mNc=a.mNc;}

  uint_4 mNr,mNc;
  NDataBlock<T> mNDBlock;
};

////////////////////////////////////////////////////////////////
// Impression

template <class T>
inline ostream& operator << (ostream& os, const TMatrix<T>& a)
                            {a.Print(os); return(os);}

////////////////////////////////////////////////////////////////
// Surcharge d'operateurs A (+=,-=,*=,/=) (T) x

template <class T> inline TMatrix<T> operator + (const TMatrix<T>& a, T b)
    {TMatrix<T> result(a); result.SetTemp(true); result += b; return result;}

template <class T> inline TMatrix<T> operator + (T b,const TMatrix<T>& a)
    {TMatrix<T> result(a); result.SetTemp(true); result += b; return result;}

template <class T> inline TMatrix<T> operator - (const TMatrix<T>& a, T b)
    {TMatrix<T> result(a); result.SetTemp(true); result -= b; return result;}

template <class T> inline TMatrix<T> operator - (T b,const TMatrix<T>& a)
    {TMatrix<T> result(a); result.SetTemp(true);
     result.DataBlock() = b-result.DataBlock(); return result;}

template <class T> inline TMatrix<T> operator * (const TMatrix<T>& a, T b)
    {TMatrix<T> result(a); result.SetTemp(true); result *= b; return result;}

template <class T> inline TMatrix<T> operator * (T b,const TMatrix<T>& a)
    {TMatrix<T> result(a); result.SetTemp(true); result *= b; return result;}

template <class T> inline TMatrix<T> operator / (const TMatrix<T>& a, T b)
    {TMatrix<T> result(a); result.SetTemp(true); result /= b; return result;}

////////////////////////////////////////////////////////////////
// Surcharge d'operateurs C = A (+,-,*,/) B

template <class T>
inline TMatrix<T> operator + (const TMatrix<T>& a,const TMatrix<T>& b)
                  {return a.Add(b);}

template <class T>
inline TMatrix<T> operator - (const TMatrix<T>& a,const TMatrix<T>& b)
                  {return a.Sub(b);}

template <class T>
inline TMatrix<T> operator * (const TMatrix<T>& a,const TMatrix<T>& b)
                  {return a.Mul(b);}

////////////////////////////////////////////////////////////////
// Typedef pour simplifier et compatibilite Peida
typedef TMatrix<r_8> Matrix;

/////////////////////////////////////////////////////////////////////////
// Classe pour la gestion de persistance
template <class T>
class FIO_TMatrix : public  PPersist  {
public:
  FIO_TMatrix();
  FIO_TMatrix(string const & filename); 
  FIO_TMatrix(const TMatrix<T> & obj);
  FIO_TMatrix(TMatrix<T> * obj);
  virtual ~FIO_TMatrix();
  virtual AnyDataObj* DataObj();
  inline operator TMatrix<T>() { return(*dobj); }
protected :
  virtual void ReadSelf(PInPersist&);           
  virtual void WriteSelf(POutPersist&) const;  
  TMatrix<T> * dobj;
  bool ownobj;
};

template <class T>
inline POutPersist& operator << (POutPersist& os, TMatrix<T> & obj)
{ FIO_TMatrix<T> fio(&obj);  fio.Write(os);  return(os); }
template <class T>
inline PInPersist& operator >> (PInPersist& is, TMatrix<T> & obj)
{ FIO_TMatrix<T> fio(&obj);  fio.Read(is);  return(is); }

/////////////////////////////////////////////////////////////////////////
// Classe de lignes/colonnes de matrices
enum TRCKind {TmatrixRow=0, TmatrixCol=1, TmatrixDiag=2};
template <class T>
class TMatrixRC {
  friend class TVector<T>;
  friend class TMatrix<T>;
public:
  TMatrixRC();

  virtual ~TMatrixRC() {}

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

protected:
  TMatrixRC(TMatrix<T>& m, TRCKind kind, uint_4 index=0);
  TMatrix<T>* matrix;
  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());}

  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.mNc;
                      case TmatrixDiag : return m.mNc+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.mNc;
                      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->mNc;
                    case TmatrixCol  : return matrix->mNr;
                    case TmatrixDiag : return matrix->mNc; }
    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;

} // Fin du namespace

#endif