// This may look like C code, but it is really -*- C++ -*- // 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-2010 #ifndef TRNGMTX_H_SEEN #define TRNGMTX_H_SEEN #include "spesqmtx.h" namespace SOPHYA { /*! \class LowerTriangularMatrix \ingroup TArray \brief Class representing a lower (inferior) triangular matrix. This is the base class for the Alm 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 / TMatrix classes, like reference sharing and counting, arithmetic operators ... However, this class has no sub matrix extraction method. */ template class LowerTriangularMatrix : public SpecialSquareMatrix { public : #include "spesqmtx_tsnl.h" //! Default constructor - TriangMatrix of size 0, SetSize() should be called before the object is used explicit LowerTriangularMatrix() : SpecialSquareMatrix(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(rowSize, C_LowerTriangularMatrix) { if (rowSize < 1) throw ParmError("LowerTriangularMatrix::LowerTriangularMatrix(rsz) rsz <= 0"); mElems.ReSize((rowSize*(rowSize+1)/2) ); mInfo = NULL; mZeros = T(0); } //! Copy constructor (possibility of sharing datas) LowerTriangularMatrix(LowerTriangularMatrix const & a, bool share=false) : SpecialSquareMatrix(a, share) { } //! Copy constructor (possibility of sharing datas) LowerTriangularMatrix(SpecialSquareMatrix const & a, bool share=false) : SpecialSquareMatrix(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 const & mx) : SpecialSquareMatrix(C_LowerTriangularMatrix) { if ((mx.NRows() != mx.NCols()) || (mx.NRows() < 1)) throw ParmError("LowerTriangularMatrix::(TMatrix const & mx) mx not allocated OR NOT a square matrix"); SetSize(mx.NRows()); for(sa_size_t l=0; l::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 ConvertToStdMatrix() const { if (mNrows < 1) throw SzMismatchError("LowerTriangularMatrix::ConvertToStdMatrix() (this) not allocated !"); SOPHYA::TMatrix mx(NRows(), NRows()); for(sa_size_t l=0; l b) //! operator = a , to set all elements to the value \b a inline LowerTriangularMatrix& operator = (T a) { SetCst(a); return (*this); } //! operator = LowerTriangularMatrix a , element by element copy operator inline LowerTriangularMatrix& operator = (LowerTriangularMatrix const & a) { Set(a); return (*this); } //! operator = Sequence seq inline LowerTriangularMatrix& 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::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::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 Multiply(LowerTriangularMatrix const & tmx) const { if (NRows() != tmx.NRows()) throw SzMismatchError("Matrix::Multiply(LowerTriangularMatrix tmx): different sizes"); // codage peu efficace : on utilise la multiplication de matrices generales ... TMatrix a = ConvertToStdMatrix(); TMatrix b = tmx.ConvertToStdMatrix(); LowerTriangularMatrix ret(a.Multiply(b)); ret.SetTemp(true); return ret; } //! Matrix product (multiplication) : ret_matrix = (*this) * mx TMatrix MultiplyTG(TMatrix const & mx) const { if (NCols() != mx.NRows()) throw SzMismatchError("LowerTriangularMatrix::MultiplyTG(TMatrix mx): NCols()!=mx.NRows()"); TMatrix a = ConvertToStdMatrix(); return a.Multiply(mx); } //! Matrix product (multiplication) : ret_matrix = mx * (*this) TMatrix MultiplyGT(TMatrix const & mx) const { if (NRows() != mx.NCols()) throw SzMismatchError("LowerTriangularMatrix::MultiplyGT(TMatrix mx): NRows()!=mx.NCols()"); TMatrix 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&,const LowerTriangularMatrix&) \brief * : LowerTriangularMatrix multiplication \b a and \b b */ template inline LowerTriangularMatrix operator * (const LowerTriangularMatrix& a, const LowerTriangularMatrix& b) { return(a.Multiply(b)); } /*! \ingroup TArray \fn operator*(const LowerTriangularMatrix&,const TMatrix&) \brief * : Matrix multiplication LowerTriangularMatrix (\b a ) * TMatrix ( \b b ) */ template inline TMatrix operator * (const LowerTriangularMatrix& a, const TMatrix& b) { return(a.MultiplyTG(b)); } /*! \ingroup TArray \fn operator*(const TMatrix&,const LowerTriangularMatrix&) \brief * : Matrix multiplication TMatrix (\b a ) * LowerTriangularMatrix ( \b b ) */ template inline TMatrix operator * (const TMatrix& a, const LowerTriangularMatrix& b) { return(b.MultiplyGT(a)); } } // namespace SOPHYA #endif