source: Sophya/trunk/SophyaLib/TArray/trngmtx.h@ 3809

// This may look like C code, but it is really -*- C++ -*-

#ifndef TRNGMTX_H_SEEN
#define TRNGMTX_H_SEEN

#include "spesqmtx.h"

// doit etre mis en dehors du namespace
/*!
  \class SOPHYA::LowerTriangularMatrix
  \ingroup TArray
  \brief Class representing a lower (inferior) triangular matrix. This is the 
  base class for the Alm<T> class (in module Samba), representing complex 
  coefficients of spherical harmonic decomposition.

  The lower triangular matrix is represented in memory as column packed,
  as illustrated below for a 5x5 triangular matrix.
  \verbatim
  5x5 Inf.Triang.Matrix, Size= 5*(5+1)/2 = 15 elements (0 ... 14) 
  | 0                  |
  | 1   5              |
  | 2   6   9          |
  | 3   7   10  12     |
  | 4   8   11  13  14 |
  \endverbatim 

  This class offers similar functionalities to the TArray<T> / TMatrix<T> classes, like 
  reference sharing and counting, arithmetic operators ... However, this class has no 
  sub matrix extraction method.
*/

namespace SOPHYA {
  
//! Class for inferior triangular matrix (base class for the class Alm)
template <class T>
class LowerTriangularMatrix : public SpecialSquareMatrix<T> {
public :

#include "spesqmtx_tsnl.h"

//! Default constructor - TriangMatrix of size 0, SetSize() should be called before the object is used
explicit LowerTriangularMatrix()
  : SpecialSquareMatrix<T>(C_LowerTriangularMatrix)
{
  mZeros = T(0);
}

//! Instanciate a triangular matrix from the number of rows (rowSize must be > 0)
explicit LowerTriangularMatrix(sa_size_t rowSize)
  : SpecialSquareMatrix<T>(rowSize, C_LowerTriangularMatrix)
{
  if (rowSize < 1) 
    throw ParmError("LowerTriangularMatrix<T>::LowerTriangularMatrix(rsz) rsz <= 0");
  mElems.ReSize((rowSize*(rowSize+1)/2) ); 
  mInfo = NULL;
  mZeros = T(0);
}

//! Copy constructor (possibility of sharing datas)
LowerTriangularMatrix(LowerTriangularMatrix<T> const & a,  bool share=false)  
  : SpecialSquareMatrix<T>(a, share) 
{
}

//! Copy constructor (possibility of sharing datas)
LowerTriangularMatrix(SpecialSquareMatrix<T> const & a,  bool share=false)  
  : SpecialSquareMatrix<T>(a, share) 
{
  if (a.MtxType() != C_LowerTriangularMatrix) 
    throw TypeMismatchExc("LowerTriangularMatrix(a) a NOT a LowerTriangularMatrix");
  mZeros = T(0);
}

/*!
 \brief Create a lower triangular matrix from a square matrix. 
 Elements above the diagonal are ignored.
*/
explicit LowerTriangularMatrix(TMatrix<T> const & mx)  
  : SpecialSquareMatrix<T>(C_LowerTriangularMatrix)
{
  if ((mx.NRows() != mx.NCols()) || (mx.NRows() < 1)) 
    throw ParmError("LowerTriangularMatrix<T>::(TMatrix<T> const & mx) mx not allocated OR NOT a square matrix");
  SetSize(mx.NRows());
  for(sa_size_t l=0; l<NRows(); l++) 
    for(sa_size_t m=0; m<=l; m++) (*this)(l,m) = mx(l,m);
  mZeros = T(0);
}

//! Sets or change the triangular matrix size, specifying the new number of rows
virtual void SetSize(sa_size_t rowSize) 
{
  if (rowSize < 1) 
    throw ParmError("LowerTriangularMatrix<T>::SetSize(rsz) rsz <= 0");
  if (rowSize == mNrows)  return;
  mNrows=rowSize;
  mElems.ReSize(mNrows*(mNrows+1)/2);  
}

//! Return number of rows (for compatibility with the old TriangularMatrix interface)
inline  sa_size_t rowNumber() const {return (int_4)mNrows;}

//! Return the object (triangular matrix) as a standard square matrix
virtual TMatrix<T> ConvertToStdMatrix() const
{
  if (mNrows < 1)
    throw SzMismatchError("LowerTriangularMatrix<T>::ConvertToStdMatrix() (this) not allocated !");
  SOPHYA::TMatrix<T> mx(NRows(), NRows());
  for(sa_size_t l=0; l<NRows(); l++) 
    for(sa_size_t m=0; m<=l; m++) mx(l,m) = (*this)(l,m);
  return mx;
}
 
//--- Operateurs = (T b) , = (LowerTriangularMatrix<T> b)
//! operator = a , to set all elements to the value \b a
inline LowerTriangularMatrix<T>& operator = (T a)
  {  SetCst(a);  return (*this);  }
//! operator = LowerTriangularMatrix<T> a , element by element copy operator
inline LowerTriangularMatrix<T>& operator = (LowerTriangularMatrix<T> const & a)
  {  Set(a); return (*this); }
//! operator = Sequence seq  
inline LowerTriangularMatrix<T>& operator = (Sequence const & seq)
  {  SetSeq(seq); return (*this); }


//--- Operateurs d'acces aux elements 
//! Element access operator (R/W): access to elements row \b r and column \b c
inline T& operator()(sa_size_t r, sa_size_t c) 
{ 
  if ((r<0)||(r>=mNrows))  
    throw RangeCheckError("DiagonalMatrix<T>::operator()(r,c) (r<0)||(r>=NRows())");
  if (c>r) { mZeros = T(0); return mZeros; }
  // the (inferior triangular )matrix is stored column by column
  return(mElems(r+ mNrows*c-c*(c+1)/2));
}
//! Element access operator (RO): access to elements row \b l and column \b m
inline T operator()(sa_size_t r, sa_size_t c) const
{ 
  if ((r<0)||(r>=mNrows))  
    throw RangeCheckError("DiagonalMatrix<T>::operator()(r,c) (r<0)||(r>=NRows())");
  if (c>r) { mZeros = T(0); return mZeros; }
  // the (inferior triangular )matrix is stored column by column
  return(mElems(r+ mNrows*c-c*(c+1)/2));
}

//! Return the pointer to the first non zero element in column \b j = &(tmmtx(j,j)) 
inline const T* columnData(sa_size_t j)  const {return mElems.Begin()+(mNrows*j-j*(j-1)/2) ;}

//! Return the pointer to the first non zero element in column \b j = &(tmmtx(j,j)) 
inline T* columnData(sa_size_t j) {return mElems.Begin()+(mNrows*j-j*(j-1)/2) ;}

//! compute the position of the element \b tm(i,j) relative to the first element
inline sa_size_t indexOfElement(sa_size_t i,sa_size_t j) const 
{
  //  return(i*(i+1)/2+j);
  // the (inferior triangular )matrix is stored column by column
  return(i+ mNrows*j-j*(j+1)/2);
}

//! Triangular Matrix product (multiplication) : ret_matrix = (*this) * tmx
LowerTriangularMatrix<T> Multiply(LowerTriangularMatrix<T> const & tmx) const
{
  if (NRows() != tmx.NRows()) 
    throw SzMismatchError("Matrix<T>::Multiply(LowerTriangularMatrix<T> tmx): different sizes");
// codage peu efficace : on utilise la multiplication de matrices generales ...
  TMatrix<T> a = ConvertToStdMatrix();
  TMatrix<T> b = tmx.ConvertToStdMatrix();
  LowerTriangularMatrix<T> ret(a.Multiply(b));
  ret.SetTemp(true);
  return ret;
}

//! Matrix product (multiplication) : ret_matrix = (*this) * mx
TMatrix<T> MultiplyTG(TMatrix<T> const & mx) const 
{
  if (NCols() != mx.NRows())
    throw SzMismatchError("LowerTriangularMatrix<T>::MultiplyTG(TMatrix<T> mx): NCols()!=mx.NRows()");
  TMatrix<T> a = ConvertToStdMatrix();
  return a.Multiply(mx); 
}

//! Matrix product (multiplication) : ret_matrix = mx * (*this)
TMatrix<T> MultiplyGT(TMatrix<T> const & mx) const 
{
  if (NRows() != mx.NCols())
    throw SzMismatchError("LowerTriangularMatrix<T>::MultiplyGT(TMatrix<T> mx): NRows()!=mx.NCols()");
  TMatrix<T> a = ConvertToStdMatrix();
  return mx.Multiply(a); 
}

//! ASCII dump/print of the triangular matrix object (set nbLignes=-1 for dumping the complete matrix)
ostream& Print(ostream& os, sa_size_t nbLignes=0) const
{
  os << "LowerTriangularMatrix< " << typeid(T).name() 
     << " > NRow=" << mNrows << " NbElem<>0 : " << Size() << endl;
  if (nbLignes == 0) return os;
  if (nbLignes < 0 ) nbLignes = mNrows;
  if (nbLignes > mNrows ) nbLignes = mNrows;
  for (sa_size_t k=0; k < nbLignes; k++)  {
    os << "L[" << k << "]: " ;
    for (sa_size_t kc = 0; kc <= k ; kc++) 
      os << " " << mElems(indexOfElement(k,kc));
    os << endl;
  }
  if (nbLignes < mNrows)  os << " ... ... ... " << endl;
  return os;
}

protected: 
mutable T mZeros; 
};

//----- Surcharge d'operateurs C = A * B (multiplication matricielle)
/*! \ingroup TArray \fn operator*(const LowerTriangularMatrix<T>&,const LowerTriangularMatrix<T>&)
  \brief * : LowerTriangularMatrix multiplication \b a and \b b */
template <class T> 
inline LowerTriangularMatrix<T> operator * (const LowerTriangularMatrix<T>& a, const LowerTriangularMatrix<T>& b)
   { return(a.Multiply(b)); }

/*! \ingroup TArray \fn operator*(const LowerTriangularMatrix<T>&,const TMatrix<T>&)
  \brief * : Matrix multiplication LowerTriangularMatrix (\b a ) *  TMatrix<T> ( \b b ) */
template <class T> 
inline TMatrix<T> operator * (const LowerTriangularMatrix<T>& a, const TMatrix<T>& b)
   { return(a.MultiplyTG(b)); }

/*! \ingroup TArray \fn operator*(const TMatrix<T>&,const LowerTriangularMatrix<T>&)
  \brief * : Matrix multiplication TMatrix (\b a ) *  LowerTriangularMatrix<T> ( \b b ) */
template <class T> 
inline TMatrix<T> operator * (const TMatrix<T>& a, const LowerTriangularMatrix<T>& b)
   { return(b.MultiplyGT(a)); }

  
}   // namespace SOPHYA

#endif