From 87138075dad1bcef1c496704f2cc3d6983156cbb Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sat, 28 Jan 2012 11:13:59 +0100 Subject: [PATCH] add the possibility to assemble a SparseMatrix object from a random list of triplets that may contain duplicated elements. It works in linear time, with O(1) re-allocations. --- Eigen/src/SparseCore/SparseMatrix.h | 138 ++++++++++++++++++++++++++++ Eigen/src/SparseCore/SparseUtil.h | 31 +++++++ test/sparse_basic.cpp | 21 +++++ 3 files changed, 190 insertions(+) diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index 587bb0fcd..eac79c826 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -431,7 +431,12 @@ class SparseMatrix //--- + template + void setFromTriplets(const InputIterators& begin, const InputIterators& end); + void sumupDuplicates(); + + //--- /** \internal * same as insert(Index,Index) except that the indices are given relative to the storage order */ @@ -899,6 +904,23 @@ protected: return (m_data.value(p) = 0); } +public: + /** \internal + * \sa insert(Index,Index) */ + inline Scalar& insertBackUncompressed(Index row, Index col) + { + const Index outer = IsRowMajor ? row : col; + const Index inner = IsRowMajor ? col : row; + + eigen_assert(!isCompressed()); + eigen_assert(m_innerNonZeros[outer]<=(m_outerIndex[outer+1] - m_outerIndex[outer])); + + Index p = m_outerIndex[outer] + m_innerNonZeros[outer]; + m_innerNonZeros[outer]++; + m_data.index(p) = inner; + return (m_data.value(p) = 0); + } + private: static void check_template_parameters() { @@ -982,4 +1004,120 @@ class SparseMatrix::ReverseInnerIterator const Index m_start; }; +namespace internal { + +template +void set_from_triplets(const InputIterator& begin, const InputIterator& end, SparseMatrixType& mat, int Options = 0) +{ + EIGEN_UNUSED_VARIABLE(Options); + enum { IsRowMajor = SparseMatrixType::IsRowMajor }; + typedef typename SparseMatrixType::Scalar Scalar; + typedef typename SparseMatrixType::Index Index; + SparseMatrix trMat(mat.rows(),mat.cols()); + + // pass 1: count the nnz per inner-vector + VectorXi wi(trMat.outerSize()); + wi.setZero(); + for(InputIterator it(begin); it!=end; ++it) + wi(IsRowMajor ? it->col() : it->row())++; + + // pass 2: insert all the elements into trMat + trMat.reserve(wi); + for(InputIterator it(begin); it!=end; ++it) + trMat.insertBackUncompressed(it->row(),it->col()) = it->value(); + + // pass 3: + trMat.sumupDuplicates(); + + // pass 4: transposed copy -> implicit sorting + mat = trMat; +} + +} + + +/** Fill the matrix \c *this with the list of \em triplets defined by the iterator range \a begin - \b. + * + * A \em triplet is a tuple (i,j,value) defining a non-zero element. + * The input list of triplets does not have to be sorted, and can contains duplicated elements. + * In any case, the result is a \b sorted and \b compressed sparse matrix where the duplicates have been summed up. + * This is a \em O(n) operation, with \em n the number of triplet elements. + * The initial contents of \c *this is destroyed. + * The matrix \c *this must be properly resized beforehand using the SparseMatrix(Index,Index) constructor, + * or the resize(Index,Index) method. The sizes are not extracted from the triplet list. + * + * The \a InputIterators value_type must provide the following interface: + * \code + * Scalar value() const; // the value + * Scalar row() const; // the row index i + * Scalar col() const; // the column index j + * \endcode + * See for instance the Eigen::Triplet template class. + * + * Here is a typical usage example: + * \code + typedef Triplet T; + std::vector tripletList; + triplets.reserve(estimation_of_entries); + for(...) + { + // ... + tripletList.push_back(T(i,j,v_ij)); + } + SparseMatrixType m(rows,cols); + m.setFromTriplets(tripletList.begin(), tripletList.end()); + // m is ready to go! + * \endcode + * + * \warning The list of triplets is read multiple times (at least twice). Therefore, it is not recommended to define + * an abstract iterator over a complex data-structure that would be expensive to evaluate. The triplets should rather + * be explicitely stored into a std::vector for instance. + */ +template +template +void SparseMatrix::setFromTriplets(const InputIterators& begin, const InputIterators& end) +{ + internal::set_from_triplets(begin, end, *this); +} + +/** \internal */ +template +void SparseMatrix::sumupDuplicates() +{ + eigen_assert(!isCompressed()); + // TODO, in practice we should be able to use m_innerNonZeros for that task + VectorXi wi(innerSize()); + wi.fill(-1); + Index count = 0; + // for each inner-vector, wi[inner_index] will hold the position of first element into the index/value buffers + for(int j=0; j=start) + { + // we already meet this entry => accumulate it + m_data.value(wi(i)) += m_data.value(k); + } + else + { + m_data.value(count) = m_data.value(k); + m_data.index(count) = m_data.index(k); + wi(i) = count; + ++count; + } + } + m_outerIndex[j] = start; + } + m_outerIndex[m_outerSize] = count; + + // turn the matrix into compressed form + delete[] m_innerNonZeros; + m_innerNonZeros = 0; + m_data.resize(m_outerIndex[m_outerSize]); +} + #endif // EIGEN_SPARSEMATRIX_H diff --git a/Eigen/src/SparseCore/SparseUtil.h b/Eigen/src/SparseCore/SparseUtil.h index a03cdbd5b..248ae7898 100644 --- a/Eigen/src/SparseCore/SparseUtil.h +++ b/Eigen/src/SparseCore/SparseUtil.h @@ -149,4 +149,35 @@ template struct plain_matrix_type } // end namespace internal +/** \ingroup SparseCore_Module + * + * \class Triplet + * + * \brief A small structure to hold a non zero as a triplet (i,j,value). + * + * \sa SparseMatrix::setFromTriplets() + */ +template +class Triplet +{ +public: + Triplet() : m_row(0), m_col(0), m_value(0) {} + + Triplet(const Index& i, const Index& j, const Scalar& v = Scalar(0)) + : m_row(i), m_col(j), m_value(v) + {} + + /** \returns the row index of the element */ + const Index& row() const { return m_row; } + + /** \returns the column index of the element */ + const Index& col() const { return m_col; } + + /** \returns the value of the element */ + const Scalar& value() const { return m_value; } +protected: + Index m_row, m_col; + Scalar m_value; +}; + #endif // EIGEN_SPARSEUTIL_H diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index 6a3e19105..5b22876b0 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -325,6 +325,27 @@ template void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(m2, refM2); } + // test setFromTriplets + { + typedef Triplet TripletType; + std::vector triplets; + int ntriplets = rows*cols; + triplets.reserve(ntriplets); + DenseMatrix refMat(rows,cols); + refMat.setZero(); + for(int i=0;i(0,rows-1); + int j = internal::random(0,cols-1); + Scalar v = internal::random(); + triplets.push_back(TripletType(i,j,v)); + refMat(i,j) += v; + } + SparseMatrixType m(rows,cols); + m.setFromTriplets(triplets.begin(), triplets.end()); + VERIFY_IS_APPROX(m, refMat); + } + // test triangularView { DenseMatrix refMat2(rows, rows), refMat3(rows, rows);