// --- Special Square matrices base class 
// This code is part of the SOPHYA library
//  (C) Univ. Paris-Sud   (C) LAL-IN2P3/CNRS   (C) IRFU-CEA
//  (C) R. Ansari, C.Magneville    2009  

#define SPESQMTX_CC_BFILE  // avoid extern template declarations 
#include "spesqmtx.h"
#include "pexceptions.h"
#include "tmatrix.h"

namespace SOPHYA {
/*!
  \class SpecialSquareMatrix
  \ingroup TArray
  \brief The base class for representing and manipulating special square matrices
  Diagonal matrix, Symmetric matrix, Triangular matrix
  This class implements the common operations for special square matrices. Only objects 

  from the sub-classes (DiagonalMatrix<T>, LowerTriangularMatrix<T>, SymmetricMatrix<T>) 
  can be instanciated. Theses classes offer similar functionalities to the TArray<T> / TMatrix<T> 
  classes, like reference sharing and arithmetic operations. However, no sub-matrix extraction 
  method is provided. 

  \code 
  // Create and initialize a diagonal matrix
  DiagonalMatrix<double> diag(3);
  diag(0,0)=1.;  diag(1,1)=2.;  diag(2,2)=3.;
  // use of multiplication by a constant operator
  diag *= 10.;
  // check the diagonal matrix
  cout << diag << diag.ConvertToStdMatrix() << endl;
  // define and initialize a general matrix
  TMatrix<double> mx(3,3);
  mx=RegularSequence();
  cout << mx;
  // check the result of a matrix mutiplication
  cout << mx*diag;
  \endcode
*/

//! Default constructor - SpecialSquareMatrix of size 0, SetSize() should be called before the object is used
template <class T>
SpecialSquareMatrix<T>::SpecialSquareMatrix(SpSqMtxType typ)
  : mType(typ) , mNrows(0) , mInfo(NULL) 
{
}

//! Instanciate a  SpecialSquareMatrix matrix from the number of rows (rowSize must be > 0)
template <class T>
SpecialSquareMatrix<T>::SpecialSquareMatrix(sa_size_t rowSize, SpSqMtxType typ)
  : mType(typ) , mNrows(rowSize) , mInfo(NULL)
{
  if (rowSize < 1) 
    throw ParmError("SpecialSquareMatrix<T>::SpecialSquareMatrix(rsz) rsz <= 0");
}

//! Copy constructor (possibility of sharing datas)
template <class T>
SpecialSquareMatrix<T>::SpecialSquareMatrix(SpecialSquareMatrix<T> const & a,  bool share)  
  : mType(a.mType) , mNrows(a.mNrows) , mElems(a.mElems, share)  , mInfo(NULL)
{
  if (a.mInfo)  mInfo = new DVList(*(a.mInfo)); 
}

//! Sets or change the matrix size, specifying the new number of rows
template <class T>
void SpecialSquareMatrix<T>::SetSize(sa_size_t rowSize)
{
  throw  ForbiddenError("SpecialSquareMatrix<T>::SetSize() Only sub-classes should be instanciated");
}

//! Return the object  as a standard (TMatrix<T>) object
template <class T>
TMatrix<T> SpecialSquareMatrix<T>::ConvertToStdMatrix() const
{
  throw  ForbiddenError("SpecialSquareMatrix<T>::ConvertToStdMatrix() Only sub-classes should be instanciated");
}




//! Clone (duplicate) the SpecialSquareMatrix<T>  matrix \b a
template <class T>
void SpecialSquareMatrix<T>::Clone(SpecialSquareMatrix<T> const & a)
{
  if ((mType != C_UndefinedMatrix) && (a.mType != mType)) 
    throw TypeMismatchExc("SpecialSquareMatrix<T>::Clone()- Different type matrices");
  mType = a.mType;
  mNrows = a.mNrows;
  mElems.Clone(a.mElems);
  if (a.mInfo)  mInfo = new DVList(*(a.mInfo)); 
}

//! Share the data with \b a 
template <class T>
void SpecialSquareMatrix<T>::Share(SpecialSquareMatrix<T> const & a)
{
  if ((mType != C_UndefinedMatrix) && (a.mType != mType)) 
    throw TypeMismatchExc("SpecialSquareMatrix<T>::Share()- Different type matrices");
  mType = a.mType;
  mNrows= a.mNrows; 
  mElems.Share(a.mElems);
  if (a.mInfo)  mInfo = new DVList(*(a.mInfo)); 
}

//! Clone or share the data. Data is shared if \b is flagged as temporary
template <class T>
void SpecialSquareMatrix<T>::CloneOrShare(SpecialSquareMatrix<T> const & a)
{
  if ((mType != C_UndefinedMatrix) && (a.mType != mType)) 
    throw TypeMismatchExc("SpecialSquareMatrix<T>::CloneOrShare()- Different type matrices");
  mType = a.mType;
  mNrows = a.mNrows;
  mElems.CloneOrShare(a.mElems);
  if (a.mInfo)  mInfo = new DVList(*(a.mInfo)); 
}



//! Set all matrix elements to the value \b a
template <class T>
SpecialSquareMatrix<T>& SpecialSquareMatrix<T>::SetCst(T a)
{
  if (mElems.Size() < 1)
    throw SzMismatchError("SpecialSquareMatrix<T>::SetCst(T )- (this) not allocated !");
  mElems = a;
  return (*this);
}


//! Element by element copy method 
template <class T>
SpecialSquareMatrix<T>& SpecialSquareMatrix<T>::Set(SpecialSquareMatrix<T> const & a)
{
  if (a.mType != mType) 
    throw TypeMismatchExc("SpecialSquareMatrix<T>::Set(a)- Different type matrices");
  if (a.mElems.Size() < 1)
    throw SzMismatchError("SpecialSquareMatrix<T>::Set(a)- Object a not allocated !");  
  if (mElems.Size() < 1) {
    mNrows = a.mNrows;
    mElems.CloneOrShare(a.mElems);
    return (*this);
  }
  if (mNrows != a.mNrows)  
    throw SzMismatchError("SpecialSquareMatrix<T>::Set(a)- NRows() != a.NRows() ! ");
  mElems = a.mElems;
  return (*this);
}

//! Set the matrix elements using the values from the sequence \b seq
template <class T>
SpecialSquareMatrix<T>& SpecialSquareMatrix<T>::SetSeq(Sequence const & seq)
{
  if (mElems.Size() < 1)
    throw SzMismatchError("SpecialSquareMatrix<T>::SetSeq()- (this) not allocated !");
  for(size_t k=0; k<mElems.Size(); k++)  seq(k).Convert(mElems(k));
  return (*this);
}

//! Add a constant value (\b b ) to all matrix elements 
template <class T>
SpecialSquareMatrix<T>& SpecialSquareMatrix<T>::AddCst(T b)
{
  if (mElems.Size() < 1)
    throw SzMismatchError("SpecialSquareMatrix<T>::AddCst(T )- (this) not allocated !");
  mElems += b;
  return (*this);
}

//! Subtract a constant value (\b b ) from all matrix elements. fginv==true -> elem = b-elem
template <class T>
SpecialSquareMatrix<T>& SpecialSquareMatrix<T>::SubCst(T b, bool fginv)
{
  if (mElems.Size() < 1)
    throw SzMismatchError("SpecialSquareMatrix<T>::Subst(T )- (this) not allocated !");
  if (fginv) 
    mElems = mElems.Sub(b, true);
  else   mElems -= b;
  return (*this);
}

//! Multiply all matrix elements by a constant value \b b . 
template <class T>
SpecialSquareMatrix<T>& SpecialSquareMatrix<T>::MulCst(T b)
{
  if (mElems.Size() < 1)
    throw SzMismatchError("SpecialSquareMatrix<T>::MulCst(T )- (this) not allocated !");
  mElems *= b; 
  return (*this);
}

//! Divide all matrix elements by a constant value \b b . if fginv==true elem = b / elem 
template <class T>
SpecialSquareMatrix<T>& SpecialSquareMatrix<T>::DivCst(T b, bool fginv)
{
  if (mElems.Size() < 1)
    throw SzMismatchError("SpecialSquareMatrix<T>::DivCst(T )- (this) not allocated !");
  if (fginv)  mElems = mElems.Div(b, true);
  else mElems /= b;
  return (*this);
}

//! Element by element addition : elem(l,m) = elem(l,m) + b.elem(l,m)
template <class T>
SpecialSquareMatrix<T>& SpecialSquareMatrix<T>::AddElt(SpecialSquareMatrix<T> const& b)
{
  if (b.mType != mType) 
    throw TypeMismatchExc("SpecialSquareMatrix<T>::AddElt() Different type matrices");
  if ((mNrows < 1)||(mNrows != b.mNrows)||(mElems.Size() < 1)||(mElems.Size() != b.mElems.Size()))
    throw SzMismatchError("SpecialSquareMatrix<T>::AddElt()- Not allocated or different sizes !");
  mElems += b.mElems;
  return (*this);
}

//! Element by element subtraction : elem(l,m) = elem(l,m) - b.elem(l,m)
template <class T>
SpecialSquareMatrix<T>& SpecialSquareMatrix<T>::SubElt(SpecialSquareMatrix<T> const& b)
{
  if (b.mType != mType) 
    throw TypeMismatchExc("SpecialSquareMatrix<T>::SubElt() Different type matrices");
  if ((mNrows < 1)||(mNrows != b.mNrows)||(mElems.Size() < 1)||(mElems.Size() != b.mElems.Size()))
    throw SzMismatchError("SpecialSquareMatrix<T>::SubElt()- Not allocated or different sizes !");
  mElems -= b.mElems;
  return (*this);
}
//! Element by element multiplication : elem(l,m) = elem(l,m) * b.elem(l,m)
template <class T>
SpecialSquareMatrix<T>& SpecialSquareMatrix<T>::MulElt(SpecialSquareMatrix<T> const& b)
{
  if (b.mType != mType) 
    throw TypeMismatchExc("SpecialSquareMatrix<T>::MulElt() Different type matrices");
  if ((mNrows < 1)||(mNrows != b.mNrows)||(mElems.Size() < 1)||(mElems.Size() != b.mElems.Size()))
    throw SzMismatchError("SpecialSquareMatrix<T>::MulElt()- Not allocated or different sizes !");
  mElems *= b.mElems;
  return (*this);
}
//! Element by element divison : elem(l,m) = elem(l,m) / b.elem(l,m)
template <class T>
SpecialSquareMatrix<T>& SpecialSquareMatrix<T>::DivElt(SpecialSquareMatrix<T> const& b)
{
  if (b.mType != mType) 
    throw TypeMismatchExc("SpecialSquareMatrix<T>::DivElt() Different type matrices");
  if ((mNrows < 1)||(mNrows != b.mNrows)||(mElems.Size() < 1)||(mElems.Size() != b.mElems.Size()))
    throw SzMismatchError("SpecialSquareMatrix<T>::DivElt()- Not allocated or different sizes !");
  mElems /= b.mElems;
  return (*this);
}


//! Return the minimum and the maximum values of the SpecialSquare Matrix elements 
/*!
  This method generates an exception (\c MathExc) if called for complex matrices
*/
template <class T>
void SpecialSquareMatrix<T>::MinMax(T& min, T& max) const 
{
  const T * pe = mElems.Begin();
  min = pe[0];
  max = pe[0];
  for(size_t k=0; k<mElems.Size(); k++) {
    if (pe[k]<min)  min = pe[k];
    else if (pe[k]>max)  max = pe[k];      
  }
  return;
}

DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void SpecialSquareMatrix< complex<r_4> >::MinMax(complex<r_4>& min, complex<r_4>& max) const
{
  throw MathExc("SpecialSquareMatrix< complex<r_4> >::MinMax(...) - No order in complex");
}
DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void SpecialSquareMatrix< complex<r_8> >::MinMax(complex<r_8>& min, complex<r_8>& max) const
{
  throw MathExc("SpecialSquareMatrix< complex<r_4> >::MinMax(...) - No order in complex");
}

//! Short ASCII representatio of the object pushed to stream \b os 
template <class T>
ostream& SpecialSquareMatrix<T>::Show(ostream& os) const
{
  os << typeid(*this).name() << " < " << typeid(T).name() 
     << " > NRow=" << mNrows << " NbElem<>0 : " << Size() << endl;
  return os;
}


///////////////////////////////////////////////////////////////
#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template SpecialSquareMatrix<uint_1>
#pragma define_template SpecialSquareMatrix<uint_2>
#pragma define_template SpecialSquareMatrix<uint_4>
#pragma define_template SpecialSquareMatrix<uint_8>
#pragma define_template SpecialSquareMatrix<int_1>
#pragma define_template SpecialSquareMatrix<int_2>
#pragma define_template SpecialSquareMatrix<int_4>
#pragma define_template SpecialSquareMatrix<int_8>
#pragma define_template SpecialSquareMatrix<r_4>
#pragma define_template SpecialSquareMatrix<r_8>
#pragma define_template SpecialSquareMatrix< complex<r_4> > 
#pragma define_template SpecialSquareMatrix< complex<r_8> > 
#ifdef SO_LDBLE128
#pragma define_template SpecialSquareMatrix<r_16>
#pragma define_template SpecialSquareMatrix< complex<r_16> > 
#endif
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class SpecialSquareMatrix<uint_1>;
template class SpecialSquareMatrix<uint_2>;
template class SpecialSquareMatrix<uint_4>;
template class SpecialSquareMatrix<uint_8>;
template class SpecialSquareMatrix<int_1>;
template class SpecialSquareMatrix<int_2>;
template class SpecialSquareMatrix<int_4>;
template class SpecialSquareMatrix<int_8>;
template class SpecialSquareMatrix<r_4>;
template class SpecialSquareMatrix<r_8>;
template class SpecialSquareMatrix< complex<r_4> >;
template class SpecialSquareMatrix< complex<r_8> >;
#ifdef SO_LDBLE128
template class SpecialSquareMatrix<r_16>;
template class SpecialSquareMatrix< complex<r_16> >;
#endif
#endif

}   // namespace SOPHYA