// 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 "tarray.h"

namespace SOPHYA {

template <class T>
class TMatrix : public TArray<T> {
public:
  // Creation / destruction 
  TMatrix();
  TMatrix(uint_4 r,uint_4 c, short mm=AutoMemoryMapping);
  TMatrix(const TMatrix<T>& a);
  TMatrix(const TMatrix<T>& a, bool share);
  TMatrix(const TArray<T>& a); 
  TMatrix(const TArray<T>& a,  bool share, short mm=CMemoryMapping);
  virtual ~TMatrix();

  //  Pour verifiez la compatibilite de dimensions lors de l'affectation
  virtual TArray<T>& Set(const TArray<T>& a);
  inline  TMatrix<T>& operator = (const TMatrix<T>& a)
                       { Set(a);  return(*this); }

  // Size - Changing the Size 
  inline uint_4 NRows() const {return Size(marowi_); }
  inline uint_4 NCols() const {return Size(macoli_); }
  inline uint_4 NCol()  const {return Size(macoli_); } // back-compat Peida

  void ReSize(uint_4 r,uint_4 c, short mm=SameMemoryMapping);  // Reallocation de place
  void Realloc(uint_4 r,uint_4 c, short mm=SameMemoryMapping, bool force=false);

  // Sub-matrix extraction $CHECK$ Reza 03/2000  Doit-on declarer cette methode const ?
  TMatrix<T> operator () (Range rline, Range rcol) const ;

  // Inline element acces methods 
  inline T const& operator()(uint_4 r,uint_4 c) const;
  inline T&       operator()(uint_4 r,uint_4 c);

  // Operations matricielles
  TMatrix<T>& Transpose();
  //mm = SameMemoryMapping or CMemoryMapping or FortranMemoryMapping 
  TMatrix<T>  Transpose(short mm); 
  // Rearranging Matrix Elements
  TMatrix<T>  Rearrange(short mm);

  // Operateur d'affectation
  // A = x (matrice diagonale Identite)
  virtual TMatrix<T>& SetIdentity(IdentityMatrix imx);
  inline  TMatrix<T>& operator = (IdentityMatrix imx) { return SetIdentity(imx); }

  // Operations diverses  avec une constante
  inline  TMatrix<T>&  operator = (T x)             { Set(x); return(*this); }
  inline  TMatrix<T>&  operator += (T x)            { Add(x); return(*this); }
  inline  TMatrix<T>&  operator -= (T x)            { Sub(x); return(*this); }
  inline  TMatrix<T>&  operator *= (T x)            { Mul(x); return(*this); }
  inline  TMatrix<T>&  operator /= (T x)            { Div(x); return(*this); }

  //  operations avec matrices 
  inline  TMatrix<T>&  operator += (const TMatrix<T>& a)  { return AddElt(a); }
  inline  TMatrix<T>&  operator -= (const TMatrix<T>& a)  { return SubElt(a); }

  // Input-Output
  virtual void  Print(ostream& os, int_4 maxprt=-1, bool si=false) const ;

  // Produit matriciel 
  TMatrix<T>  Multiply(const TMatrix<T>& b, short mm=SameMemoryMapping) const;
  //  inline TMatrix<T>& operator *= (const TMatrix<T>& b) 

protected:
};

//  ---- inline acces methods ------
template <class T>
inline T const& TMatrix<T>::operator()(uint_4 r, uint_4 c) const
{
#ifdef SO_BOUNDCHECKING
  if (marowi_ == 0)   CheckBound(r, c, 0, 0, 0, 4);
  else   CheckBound(c, r, 0, 0, 0, 4);
#endif
  return ( *( mNDBlock.Begin()+ offset_+ 
	      r*step_[marowi_] + c*step_[macoli_] ) );
}

template <class T>
inline T & TMatrix<T>::operator()(uint_4 r, uint_4 c) 
{
#ifdef SO_BOUNDCHECKING
  if (marowi_ == 0)   CheckBound(r, c, 0, 0, 0, 4);
  else   CheckBound(c, r, 0, 0, 0, 4);
#endif
  return ( *( mNDBlock.Begin()+ offset_+ 
	      r*step_[marowi_] + c*step_[macoli_] ) );
}


// Surcharge d'operateurs C = A * B

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

typedef TMatrix<r_8> Matrix;

} // Fin du namespace

#endif