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

//! Class of matrixes
/*!
  \sa TArray
 */

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);
  //! Operator = between matrices
  inline  TMatrix<T>& operator = (const TMatrix<T>& a)
                     { Set(a);  return(*this); }

  // Size - Changing the Size 
  //! return number of rows
  inline uint_4 NRows() const {return Size(marowi_); }
  //! return number of columns
  inline uint_4 NCols() const {return Size(macoli_); }
  //! return number of columns
  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 ces methode const ?
  TMatrix<T> SubMatrix(Range rline, Range rcol) const ;
  //! () : Return submatrix define by \b Range \b rline and \b rcol
  inline TMatrix<T> operator () (Range rline, Range rcol) const 
                    { return SubMatrix(rline, rcol); }
  // Lignes et colonnes de la matrice
  //! Return submatrix define by line \b ir (line vector)
  inline TMatrix<T> Row(uint_4 ir) const
                    { return SubMatrix(Range(ir,ir), Range(0,NCols()-1)); }
  //! Return submatrix define by column \b ic (column vector)
  inline TMatrix<T> Column(uint_4 ic) const
                    { return SubMatrix(Range(0,NRows()-1), Range(ic,ic)); }

  // 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);
  // = : fill matrix with an identity matrix \b imx
  inline  TMatrix<T>& operator = (IdentityMatrix imx) { return SetIdentity(imx); }

  // = : fill matrix with a Sequence \b seq
  inline  TMatrix<T>&  operator = (Sequence seq)    { SetSeq(seq); return(*this); }

  // Operations diverses  avec une constante
  //! = : fill matrix with constant value \b x
  inline  TMatrix<T>&  operator = (T x)             { SetT(x); return(*this); }
  //! += : add constant value \b x to matrix
  inline  TMatrix<T>&  operator += (T x)            { Add(x); return(*this); }
  //! -= : substract constant value \b x to matrix
  inline  TMatrix<T>&  operator -= (T x)            { Sub(x); return(*this); }
  //! *= : multiply matrix by constant value \b x
  inline  TMatrix<T>&  operator *= (T x)            { Mul(x); return(*this); }
  //! /= : divide matrix by constant value \b x
  inline  TMatrix<T>&  operator /= (T x)            { Div(x); return(*this); }

  //  operations avec matrices 
  //! += : add a matrix
  inline  TMatrix<T>&  operator += (const TMatrix<T>& a)  { AddElt(a); return(*this); }
  //! -= : substract a matrix
  inline  TMatrix<T>&  operator -= (const TMatrix<T>& a)  { SubElt(a); return(*this); }
  TMatrix<T>  Multiply(const TMatrix<T>& b, short mm=SameMemoryMapping) const;
  //! *= : matrix product : C = (*this)*B
  inline  TMatrix<T>&  operator *= (const TMatrix<T>& b)
          { this->Set(Multiply(b));  return(*this); }

  // I/O print, ...
  virtual string InfoString() const;
  virtual void  Print(ostream& os, int_4 maxprt=-1, bool si=false) const ;

protected:
};

//  ---- inline acces methods ------
 //! () : return element for line \b r and column \b c
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_] ) );
}

//! () : return element for line \b r and column \b c
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
// $CHECK$ Reza 3/4/2000 Pas necessaire  de redefinir les operateurs
// Defini au niveau de TArray<T> - Pour ameliorer l'efficacite
// Doit-on le faire aussi pour les constantes ? - Fin de $CHECK$ Reza 3/4/2000

//! + : add matrixes \b a and \b b
template <class T>
inline TMatrix<T> operator + (const TMatrix<T>& a,const TMatrix<T>& b)
    {TMatrix<T> result(a); result.SetTemp(true); result.AddElt(b); return result;}

//! - : substract matrixes \b a and \b b
template <class T>
inline TMatrix<T> operator - (const TMatrix<T>& a,const TMatrix<T>& b)
    {TMatrix<T> result(a); result.SetTemp(true); result.SubElt(b); return result;}

// Surcharge d'operateurs C = A * B
//! - : multiply matrixes \b a and \b b
template <class T> inline TMatrix<T> operator * (const TMatrix<T>& a, const TMatrix<T>& b)
{ TMatrix<T> result(a); result.SetTemp(true); return(result.Multiply(b)); }

//! Define Matrix to be TMatrix<r_8>
typedef TMatrix<r_8> Matrix;

} // Fin du namespace

#endif