Insert from triplets

This commit is contained in:
Charles Schlosser
2023-04-12 20:01:48 +00:00
committed by Rasmus Munk Larsen
parent 990a282fc4
commit 0d12fcc34e
5 changed files with 679 additions and 110 deletions

View File

@@ -151,7 +151,7 @@ class SparseMatrix
typedef typename Base::IndexVector IndexVector;
typedef typename Base::ScalarVector ScalarVector;
protected:
typedef SparseMatrix<Scalar,(Flags&~RowMajorBit)|(IsRowMajor?RowMajorBit:0),StorageIndex> TransposedSparseMatrix;
typedef SparseMatrix<Scalar, IsRowMajor ? ColMajor : RowMajor, StorageIndex> TransposedSparseMatrix;
Index m_outerSize;
Index m_innerSize;
@@ -489,6 +489,18 @@ class SparseMatrix
template<typename InputIterators, typename DupFunctor>
void setFromSortedTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func);
template<typename InputIterators>
void insertFromTriplets(const InputIterators& begin, const InputIterators& end);
template<typename InputIterators, typename DupFunctor>
void insertFromTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func);
template<typename InputIterators>
void insertFromSortedTriplets(const InputIterators& begin, const InputIterators& end);
template<typename InputIterators, typename DupFunctor>
void insertFromSortedTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func);
//---
/** \internal
@@ -1095,49 +1107,51 @@ namespace internal {
template <typename InputIterator, typename SparseMatrixType, typename DupFunctor>
void set_from_triplets(const InputIterator& begin, const InputIterator& end, SparseMatrixType& mat,
DupFunctor dup_func) {
constexpr bool IsRowMajor = SparseMatrixType::IsRowMajor;
typedef typename SparseMatrixType::StorageIndex StorageIndex;
typedef typename VectorX<StorageIndex>::AlignedMapType IndexMap;
using StorageIndex = typename SparseMatrixType::StorageIndex;
using IndexMap = typename VectorX<StorageIndex>::AlignedMapType;
using TransposedSparseMatrix = SparseMatrix<typename SparseMatrixType::Scalar, IsRowMajor ? ColMajor : RowMajor, StorageIndex>;
if (begin == end) return;
// free innerNonZeroPtr (if present) and zero outerIndexPtr
mat.resize(mat.rows(), mat.cols());
// allocate temporary storage for nonzero insertion (outer size) and duplicate removal (inner size)
ei_declare_aligned_stack_constructed_variable(StorageIndex, tmp, numext::maxi(mat.innerSize(), mat.outerSize()), 0);
// There are two strategies to consider for constructing a matrix from unordered triplets:
// A) construct the 'mat' in its native storage order and sort in-place (less memory); or,
// B) construct the transposed matrix and use an implicit sort upon assignment to `mat` (less time).
// This routine uses B) for faster execution time.
TransposedSparseMatrix trmat(mat.rows(), mat.cols());
// scan triplets to determine allocation size before constructing matrix
IndexMap outerIndexMap(mat.outerIndexPtr(), mat.outerSize() + 1);
Index nonZeros = 0;
for (InputIterator it(begin); it != end; ++it) {
eigen_assert(it->row() >= 0 && it->row() < mat.rows() && it->col() >= 0 && it->col() < mat.cols());
StorageIndex j = convert_index<StorageIndex>(IsRowMajor ? it->row() : it->col());
outerIndexMap.coeffRef(j + 1)++;
StorageIndex j = convert_index<StorageIndex>(IsRowMajor ? it->col() : it->row());
if (nonZeros == NumTraits<StorageIndex>::highest()) internal::throw_std_bad_alloc();
trmat.outerIndexPtr()[j + 1]++;
nonZeros++;
}
// finalize outer indices and allocate memory
std::partial_sum(outerIndexMap.begin(), outerIndexMap.end(), outerIndexMap.begin());
eigen_assert(nonZeros == mat.outerIndexPtr()[mat.outerSize()]);
mat.resizeNonZeros(nonZeros);
std::partial_sum(trmat.outerIndexPtr(), trmat.outerIndexPtr() + trmat.outerSize() + 1, trmat.outerIndexPtr());
eigen_assert(nonZeros == trmat.outerIndexPtr()[trmat.outerSize()]);
trmat.resizeNonZeros(nonZeros);
// use tmp to track nonzero insertions
IndexMap back(tmp, mat.outerSize());
back = outerIndexMap.head(mat.outerSize());
// construct temporary array to track insertions (outersize) and collapse duplicates (innersize)
ei_declare_aligned_stack_constructed_variable(StorageIndex, tmp, numext::maxi(mat.innerSize(), mat.outerSize()), 0);
smart_copy(trmat.outerIndexPtr(), trmat.outerIndexPtr() + trmat.outerSize(), tmp);
// push triplets to back of each inner vector
// push triplets to back of each vector
for (InputIterator it(begin); it != end; ++it) {
StorageIndex j = convert_index<StorageIndex>(IsRowMajor ? it->row() : it->col());
StorageIndex i = convert_index<StorageIndex>(IsRowMajor ? it->col() : it->row());
mat.data().index(back.coeff(j)) = i;
mat.data().value(back.coeff(j)) = it->value();
back.coeffRef(j)++;
StorageIndex j = convert_index<StorageIndex>(IsRowMajor ? it->col() : it->row());
StorageIndex i = convert_index<StorageIndex>(IsRowMajor ? it->row() : it->col());
StorageIndex k = tmp[j];
trmat.data().index(k) = i;
trmat.data().value(k) = it->value();
tmp[j]++;
}
// use tmp to collapse duplicates
IndexMap wi(tmp, mat.innerSize());
mat.collapseDuplicates(wi, dup_func);
mat.sortInnerIndices();
IndexMap wi(tmp, trmat.innerSize());
trmat.collapseDuplicates(wi, dup_func);
// implicit sorting
mat = trmat;
}
// Creates a compressed sparse matrix from a sorted range of triplets
@@ -1145,8 +1159,8 @@ template <typename InputIterator, typename SparseMatrixType, typename DupFunctor
void set_from_triplets_sorted(const InputIterator& begin, const InputIterator& end, SparseMatrixType& mat,
DupFunctor dup_func) {
constexpr bool IsRowMajor = SparseMatrixType::IsRowMajor;
typedef typename SparseMatrixType::StorageIndex StorageIndex;
typedef typename VectorX<StorageIndex>::AlignedMapType IndexMap;
using StorageIndex = typename SparseMatrixType::StorageIndex;
if (begin == end) return;
constexpr StorageIndex kEmptyIndexValue(-1);
@@ -1156,7 +1170,6 @@ void set_from_triplets_sorted(const InputIterator& begin, const InputIterator& e
StorageIndex previous_j = kEmptyIndexValue;
StorageIndex previous_i = kEmptyIndexValue;
// scan triplets to determine allocation size before constructing matrix
IndexMap outerIndexMap(mat.outerIndexPtr(), mat.outerSize() + 1);
Index nonZeros = 0;
for (InputIterator it(begin); it != end; ++it) {
eigen_assert(it->row() >= 0 && it->row() < mat.rows() && it->col() >= 0 && it->col() < mat.cols());
@@ -1166,16 +1179,16 @@ void set_from_triplets_sorted(const InputIterator& begin, const InputIterator& e
// identify duplicates by examining previous location
bool duplicate = (previous_j == j) && (previous_i == i);
if (!duplicate) {
outerIndexMap.coeffRef(j + 1)++;
if (nonZeros == NumTraits<StorageIndex>::highest()) internal::throw_std_bad_alloc();
nonZeros++;
mat.outerIndexPtr()[j + 1]++;
previous_j = j;
previous_i = i;
}
previous_j = j;
previous_i = i;
}
// finalize outer indices and allocate memory
std::partial_sum(outerIndexMap.begin(), outerIndexMap.end(), outerIndexMap.begin());
std::partial_sum(mat.outerIndexPtr(), mat.outerIndexPtr() + mat.outerSize() + 1, mat.outerIndexPtr());
eigen_assert(nonZeros == mat.outerIndexPtr()[mat.outerSize()]);
mat.resizeNonZeros(nonZeros);
@@ -1197,27 +1210,73 @@ void set_from_triplets_sorted(const InputIterator& begin, const InputIterator& e
back++;
}
}
eigen_assert(back == nonZeros);
// matrix is finalized
}
// thin wrapper around a generic binary functor to use the sparse disjunction evaulator instead of the default "arithmetic" evaulator
template<typename DupFunctor, typename LhsScalar, typename RhsScalar = LhsScalar>
struct scalar_disjunction_op
{
using result_type = typename result_of<DupFunctor(LhsScalar, RhsScalar)>::type;
scalar_disjunction_op(const DupFunctor& op) : m_functor(op) {}
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE result_type operator() (const LhsScalar& a, const RhsScalar& b) const { return m_functor(a, b); }
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const DupFunctor& functor() const { return m_functor; }
const DupFunctor& m_functor;
};
template <typename DupFunctor, typename LhsScalar, typename RhsScalar>
struct functor_traits<scalar_disjunction_op<DupFunctor, LhsScalar, RhsScalar>> : public functor_traits<DupFunctor> {};
// Creates a compressed sparse matrix from its existing entries and those from an unsorted range of triplets
template <typename InputIterator, typename SparseMatrixType, typename DupFunctor>
void insert_from_triplets(const InputIterator& begin, const InputIterator& end, SparseMatrixType& mat,
DupFunctor dup_func) {
using Scalar = typename SparseMatrixType::Scalar;
using SrcXprType = CwiseBinaryOp<scalar_disjunction_op<DupFunctor, Scalar>, const SparseMatrixType, const SparseMatrixType>;
// set_from_triplets is necessary to sort the inner indices and remove the duplicate entries
SparseMatrixType trips(mat.rows(), mat.cols());
set_from_triplets(begin, end, trips, dup_func);
SrcXprType src = mat.binaryExpr(trips, scalar_disjunction_op<DupFunctor, Scalar>(dup_func));
// the sparse assignment procedure creates a temporary matrix and swaps the final result
assign_sparse_to_sparse<SparseMatrixType, SrcXprType>(mat, src);
}
// Creates a compressed sparse matrix from its existing entries and those from an sorted range of triplets
template <typename InputIterator, typename SparseMatrixType, typename DupFunctor>
void insert_from_triplets_sorted(const InputIterator& begin, const InputIterator& end, SparseMatrixType& mat,
DupFunctor dup_func) {
using Scalar = typename SparseMatrixType::Scalar;
using SrcXprType = CwiseBinaryOp<scalar_disjunction_op<DupFunctor, Scalar>, const SparseMatrixType, const SparseMatrixType>;
/** Fill the matrix \c *this with the list of \em triplets defined by the iterator range \a begin - \a end.
// TODO: process triplets without making a copy
SparseMatrixType trips(mat.rows(), mat.cols());
set_from_triplets_sorted(begin, end, trips, dup_func);
SrcXprType src = mat.binaryExpr(trips, scalar_disjunction_op<DupFunctor, Scalar>(dup_func));
// the sparse assignment procedure creates a temporary matrix and swaps the final result
assign_sparse_to_sparse<SparseMatrixType, SrcXprType>(mat, src);
}
} // namespace internal
/** Fill the matrix \c *this with the list of \em triplets defined in the half-open range from \a begin to \a end.
*
* 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.
* The input list of triplets does not have to be sorted, and may contain 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 initial contents of \c *this are 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
* IndexType row() const; // the row index i
* IndexType col() const; // the column index j
* \endcode
* See for instance the Eigen::Triplet template class.
*
@@ -1247,20 +1306,6 @@ void SparseMatrix<Scalar,Options_,StorageIndex_>::setFromTriplets(const InputIte
internal::set_from_triplets<InputIterators, SparseMatrix<Scalar,Options_,StorageIndex_> >(begin, end, *this, internal::scalar_sum_op<Scalar,Scalar>());
}
/** The same as setFromTriplets but triplets are assumed to be pre-sorted. This is faster and requires less temporary storage.
* Two triplets `a` and `b` are appropriately ordered if:
* \code
* ColMajor: ((a.col() != b.col()) ? (a.col() < b.col()) : (a.row() < b.row())
* RowMajor: ((a.row() != b.row()) ? (a.row() < b.row()) : (a.col() < b.col())
* \endcode
*/
template<typename Scalar, int Options_, typename StorageIndex_>
template<typename InputIterators>
void SparseMatrix<Scalar, Options_, StorageIndex_>::setFromSortedTriplets(const InputIterators& begin, const InputIterators& end)
{
internal::set_from_triplets_sorted<InputIterators, SparseMatrix<Scalar, Options_, StorageIndex_> >(begin, end, *this, internal::scalar_sum_op<Scalar, Scalar>());
}
/** The same as setFromTriplets but when duplicates are met the functor \a dup_func is applied:
* \code
* value = dup_func(OldValue, NewValue)
@@ -1274,7 +1319,21 @@ template<typename Scalar, int Options_, typename StorageIndex_>
template<typename InputIterators, typename DupFunctor>
void SparseMatrix<Scalar, Options_, StorageIndex_>::setFromTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func)
{
internal::set_from_triplets<InputIterators, SparseMatrix<Scalar, Options_, StorageIndex_>, DupFunctor>(begin, end, *this, dup_func);
internal::set_from_triplets<InputIterators, SparseMatrix<Scalar, Options_, StorageIndex_>, DupFunctor>(begin, end, *this, dup_func);
}
/** The same as setFromTriplets but triplets are assumed to be pre-sorted. This is faster and requires less temporary storage.
* Two triplets `a` and `b` are appropriately ordered if:
* \code
* ColMajor: ((a.col() != b.col()) ? (a.col() < b.col()) : (a.row() < b.row())
* RowMajor: ((a.row() != b.row()) ? (a.row() < b.row()) : (a.col() < b.col())
* \endcode
*/
template<typename Scalar, int Options_, typename StorageIndex_>
template<typename InputIterators>
void SparseMatrix<Scalar, Options_, StorageIndex_>::setFromSortedTriplets(const InputIterators& begin, const InputIterators& end)
{
internal::set_from_triplets_sorted<InputIterators, SparseMatrix<Scalar, Options_, StorageIndex_> >(begin, end, *this, internal::scalar_sum_op<Scalar, Scalar>());
}
/** The same as setFromSortedTriplets but when duplicates are met the functor \a dup_func is applied:
@@ -1283,52 +1342,147 @@ void SparseMatrix<Scalar, Options_, StorageIndex_>::setFromTriplets(const InputI
* \endcode
* Here is a C++11 example keeping the latest entry only:
* \code
* mat.setFromTriplets(triplets.begin(), triplets.end(), [] (const Scalar&,const Scalar &b) { return b; });
* mat.setFromSortedTriplets(triplets.begin(), triplets.end(), [] (const Scalar&,const Scalar &b) { return b; });
* \endcode
*/
template<typename Scalar, int Options_, typename StorageIndex_>
template<typename InputIterators, typename DupFunctor>
void SparseMatrix<Scalar, Options_, StorageIndex_>::setFromSortedTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func)
{
internal::set_from_triplets_sorted<InputIterators, SparseMatrix<Scalar, Options_, StorageIndex_>, DupFunctor>(begin, end, *this, dup_func);
internal::set_from_triplets_sorted<InputIterators, SparseMatrix<Scalar, Options_, StorageIndex_>, DupFunctor>(begin, end, *this, dup_func);
}
/** Insert a batch of elements into the matrix \c *this with the list of \em triplets defined in the half-open range from \a begin to \a end.
*
* 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 may contain 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 are preserved (except for the summation of duplicate elements).
* The matrix \c *this must be properly sized beforehand. 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
* IndexType row() const; // the row index i
* IndexType col() const; // the column index j
* \endcode
* See for instance the Eigen::Triplet template class.
*
* Here is a typical usage example:
* \code
SparseMatrixType m(rows,cols); // m contains nonzero entries
typedef Triplet<double> T;
std::vector<T> tripletList;
tripletList.reserve(estimation_of_entries);
for(...)
{
// ...
tripletList.push_back(T(i,j,v_ij));
}
m.insertFromTriplets(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 explicitly stored into a std::vector for instance.
*/
template<typename Scalar, int Options_, typename StorageIndex_>
template<typename InputIterators>
void SparseMatrix<Scalar, Options_, StorageIndex_>::insertFromTriplets(const InputIterators& begin, const InputIterators& end)
{
internal::insert_from_triplets<InputIterators, SparseMatrix<Scalar, Options_, StorageIndex_> >(begin, end, *this, internal::scalar_sum_op<Scalar, Scalar>());
}
/** The same as insertFromTriplets but when duplicates are met the functor \a dup_func is applied:
* \code
* value = dup_func(OldValue, NewValue)
* \endcode
* Here is a C++11 example keeping the latest entry only:
* \code
* mat.insertFromTriplets(triplets.begin(), triplets.end(), [] (const Scalar&,const Scalar &b) { return b; });
* \endcode
*/
template<typename Scalar, int Options_, typename StorageIndex_>
template<typename InputIterators, typename DupFunctor>
void SparseMatrix<Scalar, Options_, StorageIndex_>::insertFromTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func)
{
internal::insert_from_triplets<InputIterators, SparseMatrix<Scalar, Options_, StorageIndex_>, DupFunctor>(begin, end, *this, dup_func);
}
/** The same as insertFromTriplets but triplets are assumed to be pre-sorted. This is faster and requires less temporary storage.
* Two triplets `a` and `b` are appropriately ordered if:
* \code
* ColMajor: ((a.col() != b.col()) ? (a.col() < b.col()) : (a.row() < b.row())
* RowMajor: ((a.row() != b.row()) ? (a.row() < b.row()) : (a.col() < b.col())
* \endcode
*/
template<typename Scalar, int Options_, typename StorageIndex_>
template<typename InputIterators>
void SparseMatrix<Scalar, Options_, StorageIndex_>::insertFromSortedTriplets(const InputIterators& begin, const InputIterators& end)
{
internal::insert_from_triplets_sorted<InputIterators, SparseMatrix<Scalar, Options_, StorageIndex_> >(begin, end, *this, internal::scalar_sum_op<Scalar, Scalar>());
}
/** The same as insertFromSortedTriplets but when duplicates are met the functor \a dup_func is applied:
* \code
* value = dup_func(OldValue, NewValue)
* \endcode
* Here is a C++11 example keeping the latest entry only:
* \code
* mat.insertFromSortedTriplets(triplets.begin(), triplets.end(), [] (const Scalar&,const Scalar &b) { return b; });
* \endcode
*/
template<typename Scalar, int Options_, typename StorageIndex_>
template<typename InputIterators, typename DupFunctor>
void SparseMatrix<Scalar, Options_, StorageIndex_>::insertFromSortedTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func)
{
internal::insert_from_triplets_sorted<InputIterators, SparseMatrix<Scalar, Options_, StorageIndex_>, DupFunctor>(begin, end, *this, dup_func);
}
/** \internal */
template<typename Scalar, int Options_, typename StorageIndex_>
template<typename Derived, typename DupFunctor>
void SparseMatrix<Scalar, Options_, StorageIndex_>::collapseDuplicates(DenseBase<Derived>& wi, DupFunctor dup_func)
{
eigen_assert(wi.size() >= m_innerSize);
template <typename Scalar_, int Options_, typename StorageIndex_>
template <typename Derived, typename DupFunctor>
void SparseMatrix<Scalar_, Options_, StorageIndex_>::collapseDuplicates(DenseBase<Derived>& wi, DupFunctor dup_func) {
// removes duplicate entries and compresses the matrix
// the excess allocated memory is not released
// the inner indices do not need to be sorted, nor is the matrix returned in a sorted state
eigen_assert(wi.size() == m_innerSize);
constexpr StorageIndex kEmptyIndexValue(-1);
wi.setConstant(kEmptyIndexValue);
StorageIndex count = 0;
const bool is_compressed = isCompressed();
// for each inner-vector, wi[inner_index] will hold the position of first element into the index/value buffers
for (Index j = 0; j < m_outerSize; ++j) {
StorageIndex start = count;
StorageIndex oldEnd = isCompressed() ? m_outerIndex[j + 1] : m_outerIndex[j] + m_innerNonZeros[j];
for (StorageIndex k = m_outerIndex[j]; k < oldEnd; ++k) {
const StorageIndex newBegin = count;
const StorageIndex end = is_compressed ? m_outerIndex[j + 1] : m_outerIndex[j] + m_innerNonZeros[j];
for (StorageIndex k = m_outerIndex[j]; k < end; ++k) {
StorageIndex i = m_data.index(k);
if (wi(i) >= start) {
// we already meet this entry => accumulate it
if (wi(i) >= newBegin) {
// entry at k is a duplicate
// accumulate it into the primary entry located at wi(i)
m_data.value(wi(i)) = dup_func(m_data.value(wi(i)), m_data.value(k));
} else {
// k is the primary entry in j with inner index i
// shift it to the left and record its location at wi(i)
m_data.index(count) = i;
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[j] = newBegin;
}
m_outerIndex[m_outerSize] = count;
// turn the matrix into compressed form
m_data.resize(count);
// turn the matrix into compressed form (if it is not already)
internal::conditional_aligned_delete_auto<StorageIndex, true>(m_innerNonZeros, m_outerSize);
m_innerNonZeros = 0;
m_data.resize(m_outerIndex[m_outerSize]);
}
/** \internal */
template<typename Scalar, int Options_, typename StorageIndex_>
template<typename OtherDerived>