// 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 SYMMTX_H_SEEN #define SYMMTX_H_SEEN #include "spesqmtx.h" namespace SOPHYA { /*! \class SymmetricMatrix \ingroup TArray \brief Class representing a symmetric matrix. The symmetric matrix is represented in memory as column packed, corresponding to the lower triangular part, as illustrated below for a 5x5 matrix. \verbatim 5x5 symmetric.Matrix, Size= 5*(5+1)/2 = 15 independent 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 SymmetricMatrix : public SpecialSquareMatrix { public : #include "spesqmtx_tsnl.h" //! Default constructor - TriangMatrix of size 0, SetSize() should be called before the object is used explicit SymmetricMatrix() : SpecialSquareMatrix(C_SymmetricMatrix) { } //! Instanciate a triangular matrix from the number of rows (rowSize must be > 0) explicit SymmetricMatrix(sa_size_t rowSize) : SpecialSquareMatrix(rowSize, C_SymmetricMatrix) { if (rowSize < 1) throw ParmError("SymmetricMatrix::SymmetricMatrix(rsz) rsz <= 0"); mElems.ReSize((rowSize*(rowSize+1)/2) ); mInfo = NULL; } //! Copy constructor (possibility of sharing datas) SymmetricMatrix(SymmetricMatrix const & a, bool share=false) : SpecialSquareMatrix(a, share) { } //! Copy constructor (possibility of sharing datas) SymmetricMatrix(SpecialSquareMatrix const & a, bool share=false) : SpecialSquareMatrix(a, share) { if (a.MtxType() != C_SymmetricMatrix) throw TypeMismatchExc("SymmetricMatrix(a) a NOT a SymmetricMatrix"); } /*! \brief Create a lower triangular matrix from a square matrix. Elements above the diagonal are ignored. */ explicit SymmetricMatrix(TMatrix const & mx) : SpecialSquareMatrix(C_SymmetricMatrix) { if ((mx.NRows() != mx.NCols()) || (mx.NRows() < 1)) throw ParmError("SymmetricMatrix::(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("SymmetricMatrix::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 SymmetricMatrix& operator = (T a) { SetCst(a); return (*this); } //! operator = SymmetricMatrix a , element by element copy operator inline SymmetricMatrix& operator = (SymmetricMatrix const & a) { Set(a); return (*this); } //! operator = Sequence seq inline SymmetricMatrix& 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) { sa_size_t rc = r; r=c; c=rc; } // the inferior triangular part of the 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) { sa_size_t rc = r; r=c; c=rc; } // the inferior triangular part of the 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 TMatrix Multiply(SymmetricMatrix const & tmx) const { if (NRows() != tmx.NRows()) throw SzMismatchError("Matrix::Multiply(SymmetricMatrix tmx): different sizes"); // codage peu efficace : on utilise la multiplication de matrices generales ... TMatrix a = ConvertToStdMatrix(); TMatrix b = tmx.ConvertToStdMatrix(); return (a.Multiply(b)); } //! Matrix product (multiplication) : ret_matrix = (*this) * mx TMatrix MultiplySG(TMatrix const & mx) const { if (NCols() != mx.NRows()) throw SzMismatchError("SymmetricMatrix::MultiplySG(TMatrix mx): NCols()!=mx.NRows()"); TMatrix a = ConvertToStdMatrix(); return a.Multiply(mx); } //! Matrix product (multiplication) : ret_matrix = mx * (*this) TMatrix MultiplyGS(TMatrix const & mx) const { if (NRows() != mx.NCols()) throw SzMismatchError("SymmetricMatrix::MultiplyGS(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 << "SymmetricMatrix< " << 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 r=0; r&,const SymmetricMatrix&) \brief * : SymmetricMatrix multiplication \b a and \b b */ template inline TMatrix operator * (const SymmetricMatrix& a, const SymmetricMatrix& b) { return(a.Multiply(b)); } /*! \ingroup TArray \fn operator*(const SymmetricMatrix&,const TMatrix&) \brief * : Matrix multiplication SymmetricMatrix (\b a ) * TMatrix ( \b b ) */ template inline TMatrix operator * (const SymmetricMatrix& a, const TMatrix& b) { return(a.MultiplySG(b)); } /*! \ingroup TArray \fn operator*(const TMatrix&,const SymmetricMatrix&) \brief * : Matrix multiplication TMatrix (\b a ) * SymmetricMatrix ( \b b ) */ template inline TMatrix operator * (const TMatrix& a, const SymmetricMatrix& b) { return(b.MultiplyGS(a)); } } // namespace SOPHYA #endif