source: Sophya/trunk/SophyaLib/TArray/triangmtx.h@ 3572

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

#ifndef TRIANGMTX_H_SEEN
#define TRIANGMTX_H_SEEN

#include "ndatablock.h"
#include "pexceptions.h"

// doit etre mis en dehors du namespace
/*!
  \class SOPHYA::TriangularMatrix
  \ingroup TArray
  \brief Class for inferior triangular matrix (base class for the class Alm)
  The inferior triangular matrix is represented in memory as column packed,
  as illustrated below for a 5x5 triangular matrix.
  \verbatim
  5x5 Inf.Triang.Matrix, Size= 15 elements (0 ... 14) 
  | 0                  |
  | 1   5              |
  | 2   6   9          |
  | 3   7   10  12     |
  | 4   8   11  13  14 |
  \endverbatim 
*/

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

//! Default constructor
TriangularMatrix()
  : long_diag_(0)
{
}

//! instanciate a triangular matrix from the number of rows
TriangularMatrix(sa_size_t rowSize)
  : long_diag_(rowSize) 
{
  elem_.ReSize((rowSize*(rowSize+1)/2) ); 
}

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

//! resize the matrix with a new number of rows
inline void ReSizeRow(sa_size_t rowSize) 
{
  long_diag_=(uint_4)rowSize;
  elem_.ReSize(long_diag_*(long_diag_+1)/2);
}

TriangularMatrix<T>& SetT(T a)
  {
    if (long_diag_ < 1)
    throw RangeCheckError("TriangularMatrix<T>::SetT(T )  - TriangularMatrix not dimensionned ! ");
    elem_ = a;
    return (*this);
  }

//! () operator : access to elements row \b l and column \b m
inline T& operator()(sa_size_t l, sa_size_t m) 
  {
      return  elem_(indexOfElement(l,m));
  }

inline T& operator()(sa_size_t index) 
  {
      return  elem_(index);
  }


//! () operator : access to elements row \b l and column \b m
inline T const& operator()(sa_size_t l, sa_size_t m) const 
  {
      return *(elem_.Begin()+ indexOfElement(l,m));
  }

inline T const& operator()(sa_size_t index) const 
  {
      return *(elem_.Begin()+ index);
  }

TriangularMatrix<T>& Set(const TriangularMatrix<T>& a)
  {
    if (this != &a) 
      {
        if (a.Size() < 1) 
          throw RangeCheckError(" TriangularMatrix<T>::Set()- Array a not allocated ! ");
      }
    if (Size() < 1)  CloneOrShare(a);
    else CopyElt(a);
    return(*this);
  }

inline TriangularMatrix<T>& operator = (const TriangularMatrix<T>& a)
                                                           {return Set(a);}

TriangularMatrix<T>& CopyElt(const  TriangularMatrix<T>& a)
  {
    if (Size() < 1) 
      throw RangeCheckError("TriangularMatrix<T>::CopyElt(const TriangularMatrix<T>& )  - Not Allocated Array ! ");
    if (Size() != a.Size() )
      throw(SzMismatchError("TriangularMatrix<T>::CopyElt(const TriangularMatrix<T>&) SizeMismatch")) ;
    long_diag_ = a.long_diag_;
    sa_size_t k;
    for (k=0; k< Size(); k++) elem_(k) = a.elem_(k);
    return(*this);
  }

void CloneOrShare(const  TriangularMatrix<T>& a)
  {
    long_diag_ = a.long_diag_;
    elem_.CloneOrShare(a.elem_);
  }


//! Return number of rows
inline  sa_size_t rowNumber() const {return (int_4)long_diag_;}

//! Return size of the total array
  inline sa_size_t Size() const {return elem_.Size();}

  inline bool CheckRelativeIndices(sa_size_t l, sa_size_t m) const 
   {
     if ( l < m ) 
       {
         throw RangeCheckError("TriangularMatrix<T>::CheckRelativeIndices: indices out of range " );
       }
     return true;
   }
  inline bool CheckAbsoluteIndice(sa_size_t l, sa_size_t m) const 
   {
     if ( indexOfElement(l,m) >= elem_.Size() )
       {
         throw RangeCheckError("TriangularMatrix<T>::CheckAbsoluteIndice: indices out of range " );
       }
   }
  inline bool CheckAbsoluteIndice(sa_size_t ind) const 
   {
     if ( ind >= elem_.Size() )
       {
         throw RangeCheckError("TriangularMatrix<T>::CheckAbsoluteIndice: indices out of range " );
       }
   }

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

inline void Print(sa_size_t nbLignes=0) const  { Print(cout, nbLignes); }


//! 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 elem_.Begin()+(long_diag_*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 elem_.Begin()+(long_diag_*j-j*(j-1)/2) ;}

//! compute the address of an element in the single array representing the matrix
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+ long_diag_*j-j*(j+1)/2);
}

private: 

sa_size_t long_diag_;    //!< size of the square matrix
NDataBlock<T> elem_;  //!< Data block

  };

template <class T>
inline ostream& operator << (ostream& os, const TriangularMatrix<T>& a)
                            { a.Print(os, 0);    return(os);    }
  
}   // namespace SOPHYA

#endif