// 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 TriangularMatrix { public : //! Default constructor TriangularMatrix() {;}; //! 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& a, bool share=false) : elem_(a.elem_, share), long_diag_(a.long_diag_) {;} //! 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& SetT(T a) { if (long_diag_ < 1) throw RangeCheckError("TriangularMatrix::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& Set(const TriangularMatrix& a) { if (this != &a) { if (a.Size() < 1) throw RangeCheckError(" TriangularMatrix::Set()- Array a not allocated ! "); } if (Size() < 1) CloneOrShare(a); else CopyElt(a); return(*this); } inline TriangularMatrix& operator = (const TriangularMatrix& a) {return Set(a);} TriangularMatrix& CopyElt(const TriangularMatrix& a) { if (Size() < 1) throw RangeCheckError("TriangularMatrix::CopyElt(const TriangularMatrix& ) - Not Allocated Array ! "); if (Size() != a.Size() ) throw(SzMismatchError("TriangularMatrix::CopyElt(const TriangularMatrix&) 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& 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::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::CheckAbsoluteIndice: indices out of range " ); } } inline bool CheckAbsoluteIndice(sa_size_t ind) const { if ( ind >= elem_.Size() ) { throw RangeCheckError("TriangularMatrix::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 elem_; //!< Data block }; template inline ostream& operator << (ostream& os, const TriangularMatrix& a) { a.Print(os, 0); return(os); } } // namespace SOPHYA #endif