diff --git a/Eigen/SparseCore b/Eigen/SparseCore index 56a9401af34a61f2bdfb2d646ddc71d9f2aa1b7d..e0b24147b4ba068dff20b2574bd891231e264f92 100644 --- a/Eigen/SparseCore +++ b/Eigen/SparseCore @@ -35,6 +35,7 @@ */ // IWYU pragma: begin_exports +#include "src/SparseCore/IndexStorage.h" #include "src/SparseCore/SparseUtil.h" #include "src/SparseCore/SparseMatrixBase.h" #include "src/SparseCore/SparseAssign.h" diff --git a/Eigen/src/CholmodSupport/CholmodSupport.h b/Eigen/src/CholmodSupport/CholmodSupport.h index 7e3c881aec9891435d82b082bce449a96bbb9944..4fe63d6cd3626c09bfc6cbcebec23277c163cd84 100644 --- a/Eigen/src/CholmodSupport/CholmodSupport.h +++ b/Eigen/src/CholmodSupport/CholmodSupport.h @@ -60,8 +60,8 @@ struct cholmod_configure_matrix > { /** Wraps the Eigen sparse matrix \a mat into a Cholmod sparse matrix object. * Note that the data are shared. */ -template -cholmod_sparse viewAsCholmod(Ref > mat) { +template +cholmod_sparse viewAsCholmod(Ref > mat) { cholmod_sparse res; res.nzmax = mat.nonZeros(); res.nrow = mat.rows(); @@ -98,23 +98,27 @@ cholmod_sparse viewAsCholmod(Ref return res; } -template -const cholmod_sparse viewAsCholmod(const SparseMatrix& mat) { - cholmod_sparse res = viewAsCholmod(Ref >(mat.const_cast_derived())); +template +const cholmod_sparse viewAsCholmod(const SparseMatrix& mat) { + cholmod_sparse res = + viewAsCholmod(Ref >(mat.const_cast_derived())); return res; } -template -const cholmod_sparse viewAsCholmod(const SparseVector& mat) { - cholmod_sparse res = viewAsCholmod(Ref >(mat.const_cast_derived())); +template +const cholmod_sparse viewAsCholmod(const SparseVector& mat) { + cholmod_sparse res = + viewAsCholmod(Ref >(mat.const_cast_derived())); return res; } /** Returns a view of the Eigen sparse matrix \a mat as Cholmod sparse matrix. * The data are not copied but shared. */ template -cholmod_sparse viewAsCholmod(const SparseSelfAdjointView, UpLo>& mat) { - cholmod_sparse res = viewAsCholmod(Ref >(mat.matrix().const_cast_derived())); +cholmod_sparse viewAsCholmod( + const SparseSelfAdjointView, UpLo>& mat) { + cholmod_sparse res = viewAsCholmod( + Ref >(mat.matrix().const_cast_derived())); if (UpLo == Upper) res.stype = 1; if (UpLo == Lower) res.stype = -1; diff --git a/Eigen/src/SparseCore/CompressedStorage.h b/Eigen/src/SparseCore/CompressedStorage.h index 8f8a6963a983e49dbc2e913c0b271924b13db6b1..a8f04f0e688d762837a0d314943c07d63829e799 100644 --- a/Eigen/src/SparseCore/CompressedStorage.h +++ b/Eigen/src/SparseCore/CompressedStorage.h @@ -17,11 +17,52 @@ namespace Eigen { namespace internal { +template +struct InnerStorage { + Scalar values[MaxSize_]; + StorageIndex indices[MaxSize_]; + Index size = MaxSize_; + Index allocatedSize = MaxSize_; + + InnerStorage() = default; + ~InnerStorage() = default; +}; + +template +struct InnerStorage { + Scalar* values; + StorageIndex* indices; + Index size; + Index allocatedSize; + + InnerStorage() : values(0), indices(0), size(0), allocatedSize(0) {} + + ~InnerStorage() { + conditional_aligned_delete_auto(values, allocatedSize); + conditional_aligned_delete_auto(indices, allocatedSize); + } +}; + +template +void inner_storage_swap(InnerStorage& lhs, + InnerStorage& rhs) { + std::swap(lhs.values, rhs.values); + std::swap(lhs.indices, rhs.indices); + std::swap(lhs.size, rhs.size); + std::swap(lhs.allocatedSize, rhs.allocatedSize); +} + +template +void inner_storage_assign_impl(InnerStorage_& lhs, const InnerStorage_& rhs) { + internal::smart_copy(rhs.values, rhs.values + lhs.size, lhs.values); + internal::smart_copy(rhs.indices, rhs.indices + lhs.size, lhs.indices); +} + /** \internal * Stores a sparse set of values as a list of values and a list of indices. * */ -template +template class CompressedStorage { public: typedef Scalar_ Scalar; @@ -30,47 +71,45 @@ class CompressedStorage { protected: typedef typename NumTraits::Real RealScalar; + template ::type = 0> + void before_equal(const CompressedStorage& other) { + resize(other.size()); + } + + template ::type = 0> + void before_equal(const CompressedStorage& other) { + eigen_assert(other.size() == MaxSize); + } + public: - CompressedStorage() : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0) {} + CompressedStorage() : m_storage() {} - explicit CompressedStorage(Index size) : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0) { resize(size); } + explicit CompressedStorage(Index size) : m_storage() { resize(size); } - CompressedStorage(const CompressedStorage& other) : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0) { - *this = other; - } + CompressedStorage(const CompressedStorage& other) : m_storage() { *this = other; } CompressedStorage& operator=(const CompressedStorage& other) { - resize(other.size()); - if (other.size() > 0) { - internal::smart_copy(other.m_values, other.m_values + m_size, m_values); - internal::smart_copy(other.m_indices, other.m_indices + m_size, m_indices); - } + before_equal(other); + inner_storage_assign_impl(m_storage, other.m_storage); return *this; } - void swap(CompressedStorage& other) { - std::swap(m_values, other.m_values); - std::swap(m_indices, other.m_indices); - std::swap(m_size, other.m_size); - std::swap(m_allocatedSize, other.m_allocatedSize); - } + void swap(CompressedStorage& other) { inner_storage_swap(m_storage, other.m_storage); } - ~CompressedStorage() { - conditional_aligned_delete_auto(m_values, m_allocatedSize); - conditional_aligned_delete_auto(m_indices, m_allocatedSize); - } + ~CompressedStorage() = default; void reserve(Index size) { - Index newAllocatedSize = m_size + size; - if (newAllocatedSize > m_allocatedSize) reallocate(newAllocatedSize); + Index newAllocatedSize = m_storage.size + size; + if (newAllocatedSize > m_storage.allocatedSize) reallocate(newAllocatedSize); } void squeeze() { - if (m_allocatedSize > m_size) reallocate(m_size); + if (m_storage.allocatedSize > m_storage.size) reallocate(m_storage.size); } void resize(Index size, double reserveSizeFactor = 0) { - if (m_allocatedSize < size) { + eigen_assert(MaxSize == Dynamic && "CompressedStorage can be resized only for Dynamic sizes"); + if (m_storage.allocatedSize < size) { // Avoid underflow on the std::min call by choosing the smaller index type. using SmallerIndexType = typename std::conditional((std::numeric_limits::max)()) < @@ -81,104 +120,106 @@ class CompressedStorage { if (realloc_size < size) internal::throw_std_bad_alloc(); reallocate(realloc_size); } - m_size = size; + m_storage.size = size; } void append(const Scalar& v, Index i) { - Index id = m_size; - resize(m_size + 1, 1); - m_values[id] = v; - m_indices[id] = internal::convert_index(i); + Index id = m_storage.size; + resize(m_storage.size + 1, 1); + m_storage.values[id] = v; + m_storage.indices[id] = internal::convert_index(i); } - inline Index size() const { return m_size; } - inline Index allocatedSize() const { return m_allocatedSize; } - inline void clear() { m_size = 0; } + inline Index size() const { return m_storage.size; } + inline Index allocatedSize() const { return m_storage.allocatedSize; } + inline void clear() { m_storage.size = 0; } - const Scalar* valuePtr() const { return m_values; } - Scalar* valuePtr() { return m_values; } - const StorageIndex* indexPtr() const { return m_indices; } - StorageIndex* indexPtr() { return m_indices; } + const Scalar* valuePtr() const { return m_storage.values; } + Scalar* valuePtr() { return m_storage.values; } + const StorageIndex* indexPtr() const { return m_storage.indices; } + StorageIndex* indexPtr() { return m_storage.indices; } inline Scalar& value(Index i) { - eigen_internal_assert(m_values != 0); - return m_values[i]; + eigen_internal_assert(m_storage.values != 0); + return m_storage.values[i]; } inline const Scalar& value(Index i) const { - eigen_internal_assert(m_values != 0); - return m_values[i]; + eigen_internal_assert(m_storage.values != 0); + return m_storage.values[i]; } inline StorageIndex& index(Index i) { - eigen_internal_assert(m_indices != 0); - return m_indices[i]; + eigen_internal_assert(m_storage.indices != 0); + return m_storage.indices[i]; } inline const StorageIndex& index(Index i) const { - eigen_internal_assert(m_indices != 0); - return m_indices[i]; + eigen_internal_assert(m_storage.indices != 0); + return m_storage.indices[i]; } /** \returns the largest \c k such that for all \c j in [0,k) index[\c j]\<\a key */ - inline Index searchLowerIndex(Index key) const { return searchLowerIndex(0, m_size, key); } + inline Index searchLowerIndex(Index key) const { return searchLowerIndex(0, m_storage.size, key); } /** \returns the largest \c k in [start,end) such that for all \c j in [start,k) index[\c j]\<\a key */ inline Index searchLowerIndex(Index start, Index end, Index key) const { - return static_cast(std::distance(m_indices, std::lower_bound(m_indices + start, m_indices + end, key))); + return static_cast( + std::distance(m_storage.indices, std::lower_bound(m_storage.indices + start, m_storage.indices + end, key))); } /** \returns the stored value at index \a key * If the value does not exist, then the value \a defaultValue is returned without any insertion. */ inline Scalar at(Index key, const Scalar& defaultValue = Scalar(0)) const { - if (m_size == 0) + if (m_storage.size == 0) return defaultValue; - else if (key == m_indices[m_size - 1]) - return m_values[m_size - 1]; + else if (key == m_storage.indices[m_storage.size - 1]) + return m_storage.values[m_storage.size - 1]; // ^^ optimization: let's first check if it is the last coefficient // (very common in high level algorithms) - const Index id = searchLowerIndex(0, m_size - 1, key); - return ((id < m_size) && (m_indices[id] == key)) ? m_values[id] : defaultValue; + const Index id = searchLowerIndex(0, m_storage.size - 1, key); + return ((id < m_storage.size) && (m_storage.indices[id] == key)) ? m_storage.values[id] : defaultValue; } /** Like at(), but the search is performed in the range [start,end) */ inline Scalar atInRange(Index start, Index end, Index key, const Scalar& defaultValue = Scalar(0)) const { if (start >= end) return defaultValue; - else if (end > start && key == m_indices[end - 1]) - return m_values[end - 1]; + else if (end > start && key == m_storage.indices[end - 1]) + return m_storage.values[end - 1]; // ^^ optimization: let's first check if it is the last coefficient // (very common in high level algorithms) const Index id = searchLowerIndex(start, end - 1, key); - return ((id < end) && (m_indices[id] == key)) ? m_values[id] : defaultValue; + return ((id < end) && (m_storage.indices[id] == key)) ? m_storage.values[id] : defaultValue; } /** \returns a reference to the value at index \a key * If the value does not exist, then the value \a defaultValue is inserted * such that the keys are sorted. */ inline Scalar& atWithInsertion(Index key, const Scalar& defaultValue = Scalar(0)) { - Index id = searchLowerIndex(0, m_size, key); - if (id >= m_size || m_indices[id] != key) { - if (m_allocatedSize < m_size + 1) { - Index newAllocatedSize = 2 * (m_size + 1); - m_values = conditional_aligned_realloc_new_auto(m_values, newAllocatedSize, m_allocatedSize); - m_indices = - conditional_aligned_realloc_new_auto(m_indices, newAllocatedSize, m_allocatedSize); - m_allocatedSize = newAllocatedSize; + Index id = searchLowerIndex(0, m_storage.size, key); + if (id >= m_storage.size || m_storage.indices[id] != key) { + if (m_storage.allocatedSize < m_storage.size + 1) { + Index newAllocatedSize = 2 * (m_storage.size + 1); + m_storage.values = conditional_aligned_realloc_new_auto(m_storage.values, newAllocatedSize, + m_storage.allocatedSize); + m_storage.indices = conditional_aligned_realloc_new_auto( + m_storage.indices, newAllocatedSize, m_storage.allocatedSize); + m_storage.allocatedSize = newAllocatedSize; } - if (m_size > id) { - internal::smart_memmove(m_values + id, m_values + m_size, m_values + id + 1); - internal::smart_memmove(m_indices + id, m_indices + m_size, m_indices + id + 1); + if (m_storage.size > id) { + internal::smart_memmove(m_storage.values + id, m_storage.values + m_storage.size, m_storage.values + id + 1); + internal::smart_memmove(m_storage.indices + id, m_storage.indices + m_storage.size, m_storage.indices + id + 1); } - m_size++; - m_indices[id] = internal::convert_index(key); - m_values[id] = defaultValue; + m_storage.size++; + m_storage.indices[id] = internal::convert_index(key); + m_storage.values[id] = defaultValue; } - return m_values[id]; + return m_storage.values[id]; } inline void moveChunk(Index from, Index to, Index chunkSize) { - eigen_internal_assert(chunkSize >= 0 && to + chunkSize <= m_size); - internal::smart_memmove(m_values + from, m_values + from + chunkSize, m_values + to); - internal::smart_memmove(m_indices + from, m_indices + from + chunkSize, m_indices + to); + eigen_internal_assert(chunkSize >= 0 && to + chunkSize <= m_storage.size); + internal::smart_memmove(m_storage.values + from, m_storage.values + from + chunkSize, m_storage.values + to); + internal::smart_memmove(m_storage.indices + from, m_storage.indices + from + chunkSize, m_storage.indices + to); } protected: @@ -186,21 +227,21 @@ class CompressedStorage { #ifdef EIGEN_SPARSE_COMPRESSED_STORAGE_REALLOCATE_PLUGIN EIGEN_SPARSE_COMPRESSED_STORAGE_REALLOCATE_PLUGIN #endif - eigen_internal_assert(size != m_allocatedSize); - m_values = conditional_aligned_realloc_new_auto(m_values, size, m_allocatedSize); - m_indices = conditional_aligned_realloc_new_auto(m_indices, size, m_allocatedSize); - m_allocatedSize = size; + eigen_internal_assert(size != m_storage.allocatedSize); + m_storage.values = + conditional_aligned_realloc_new_auto(m_storage.values, size, m_storage.allocatedSize); + m_storage.indices = + conditional_aligned_realloc_new_auto(m_storage.indices, size, m_storage.allocatedSize); + m_storage.allocatedSize = size; } protected: - Scalar* m_values; - StorageIndex* m_indices; - Index m_size; - Index m_allocatedSize; + using Storage = InnerStorage; + Storage m_storage; }; } // end namespace internal } // end namespace Eigen -#endif // EIGEN_COMPRESSED_STORAGE_H +#endif // EIGEN_COMPRESSED_STORAGE_H \ No newline at end of file diff --git a/Eigen/src/SparseCore/IndexStorage.h b/Eigen/src/SparseCore/IndexStorage.h new file mode 100644 index 0000000000000000000000000000000000000000..c0d717237f01322d0b606c4aefbbcf7bc498ffcd --- /dev/null +++ b/Eigen/src/SparseCore/IndexStorage.h @@ -0,0 +1,46 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_INDEXSTORAGE_STORAGE_H +#define EIGEN_INDEXSTORAGE_STORAGE_H + +// IWYU pragma: private +#include "./InternalHeaderCheck.h" + +namespace Eigen { + +namespace internal { + +template +struct IndexStorage { + StorageIndex outerIndex[MaxSize_ + 1] = {0}; + StorageIndex* innerNonZeros; // not used in static size + StorageIndex* outerPtr() { return outerIndex; } + const StorageIndex* outerPtr() const { return outerIndex; } + StorageIndex* innerNZPtr() { return innerNonZeros; } + const StorageIndex* innerNZPtr() const { return innerNonZeros; } + IndexStorage() : innerNonZeros(0) {} + ~IndexStorage() = default; +}; + +template +struct IndexStorage { + StorageIndex* outerIndex; + StorageIndex* innerNonZeros; // optional, if null then the data is compressed + StorageIndex* outerPtr() { return outerIndex; } + const StorageIndex* outerPtr() const { return outerIndex; } + StorageIndex* innerNZPtr() { return innerNonZeros; } + const StorageIndex* innerNZPtr() const { return innerNonZeros; } + IndexStorage() : outerIndex(0), innerNonZeros(0) {} + ~IndexStorage() = default; +}; + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_INDEXSTORAGE_STORAGE_H \ No newline at end of file diff --git a/Eigen/src/SparseCore/SparseBlock.h b/Eigen/src/SparseCore/SparseBlock.h index 1342f4e7bd7a6c483b9b308543d17ce5a3f49896..512c471c3f892cd71c55edb4f1609d135b1746a8 100644 --- a/Eigen/src/SparseCore/SparseBlock.h +++ b/Eigen/src/SparseCore/SparseBlock.h @@ -257,12 +257,15 @@ class sparse_matrix_block_impl : public SparseCompressedBase -class BlockImpl, BlockRows, BlockCols, true, Sparse> - : public internal::sparse_matrix_block_impl, BlockRows, BlockCols> { +template +class BlockImpl, BlockRows, BlockCols, true, + Sparse> + : public internal::sparse_matrix_block_impl, + BlockRows, BlockCols> { public: typedef StorageIndex_ StorageIndex; - typedef SparseMatrix SparseMatrixType; + typedef SparseMatrix SparseMatrixType; typedef internal::sparse_matrix_block_impl Base; inline BlockImpl(SparseMatrixType& xpr, Index i) : Base(xpr, i) {} @@ -272,13 +275,15 @@ class BlockImpl, BlockRows, Block using Base::operator=; }; -template -class BlockImpl, BlockRows, BlockCols, true, Sparse> - : public internal::sparse_matrix_block_impl, BlockRows, - BlockCols> { +template +class BlockImpl, BlockRows, BlockCols, true, + Sparse> + : public internal::sparse_matrix_block_impl< + const SparseMatrix, BlockRows, BlockCols> { public: typedef StorageIndex_ StorageIndex; - typedef const SparseMatrix SparseMatrixType; + typedef const SparseMatrix SparseMatrixType; typedef internal::sparse_matrix_block_impl Base; inline BlockImpl(SparseMatrixType& xpr, Index i) : Base(xpr, i) {} @@ -508,21 +513,28 @@ class unary_evaluator, Iterator inline operator bool() const { return m_outerPos < m_end; } }; -template -struct unary_evaluator, BlockRows, BlockCols, true>, IteratorBased> - : evaluator< - SparseCompressedBase, BlockRows, BlockCols, true> > > { - typedef Block, BlockRows, BlockCols, true> XprType; +template +struct unary_evaluator< + Block, BlockRows, BlockCols, true>, + IteratorBased> + : evaluator, BlockRows, BlockCols, true> > > { + typedef Block, BlockRows, BlockCols, true> + XprType; typedef evaluator > Base; explicit unary_evaluator(const XprType& xpr) : Base(xpr) {} }; -template -struct unary_evaluator, BlockRows, BlockCols, true>, - IteratorBased> - : evaluator, BlockRows, BlockCols, true> > > { - typedef Block, BlockRows, BlockCols, true> XprType; +template +struct unary_evaluator< + Block, BlockRows, BlockCols, true>, + IteratorBased> + : evaluator, + BlockRows, BlockCols, true> > > { + typedef Block, BlockRows, BlockCols, true> + XprType; typedef evaluator > Base; explicit unary_evaluator(const XprType& xpr) : Base(xpr) {} }; diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index 7ddb6fab901621f15ff026d4ec405509f13abf68..7daa66e420e06b15dfdfec33f361baa4de99cd15 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -47,26 +47,27 @@ namespace Eigen { */ namespace internal { -template -struct traits> { +template +struct traits> { typedef Scalar_ Scalar; typedef StorageIndex_ StorageIndex; typedef Sparse StorageKind; typedef MatrixXpr XprKind; enum { - RowsAtCompileTime = Dynamic, - ColsAtCompileTime = Dynamic, - MaxRowsAtCompileTime = Dynamic, - MaxColsAtCompileTime = Dynamic, + RowsAtCompileTime = Rows_, + ColsAtCompileTime = Cols_, + MaxRowsAtCompileTime = Rows_, + MaxColsAtCompileTime = Cols_, + MaxNZ = MaxNZ_, Options = Options_, Flags = Options_ | NestByRefBit | LvalueBit | CompressedAccessBit, SupportedAccessPatterns = InnerRandomAccessPattern }; }; -template -struct traits, DiagIndex>> { - typedef SparseMatrix MatrixType; +template +struct traits, DiagIndex>> { + typedef SparseMatrix MatrixType; typedef typename ref_selector::type MatrixTypeNested; typedef std::remove_reference_t MatrixTypeNested_; @@ -76,17 +77,18 @@ struct traits, DiagIndex typedef MatrixXpr XprKind; enum { - RowsAtCompileTime = Dynamic, + RowsAtCompileTime = Rows_, ColsAtCompileTime = 1, - MaxRowsAtCompileTime = Dynamic, + MaxRowsAtCompileTime = Rows_, MaxColsAtCompileTime = 1, + MaxNZ = MaxNZ_, Flags = LvalueBit }; }; -template -struct traits, DiagIndex>> - : public traits, DiagIndex>> { +template +struct traits, DiagIndex>> + : public traits, DiagIndex>> { enum { Flags = 0 }; }; @@ -117,8 +119,8 @@ struct functor_traits> { } // end namespace internal -template -class SparseMatrix : public SparseCompressedBase> { +template +class SparseMatrix : public SparseCompressedBase> { typedef SparseCompressedBase Base; using Base::convert_index; friend class SparseVector; @@ -139,8 +141,8 @@ class SparseMatrix : public SparseCompressedBase Storage; - enum { Options = Options_ }; + typedef internal::CompressedStorage Storage; + enum { Options = Options_, isStatic = MaxNZ_ != Dynamic }; typedef typename Base::IndexVector IndexVector; typedef typename Base::ScalarVector ScalarVector; @@ -150,10 +152,10 @@ class SparseMatrix : public SparseCompressedBase m_idx_data; + public: /** \returns the number of rows of the matrix */ inline Index rows() const { return IsRowMajor ? m_outerSize : m_innerSize; } @@ -186,20 +188,20 @@ class SparseMatrix : public SparseCompressedBase= 0 && row < rows() && col >= 0 && col < cols()); const Index outer = IsRowMajor ? row : col; const Index inner = IsRowMajor ? col : row; - Index start = m_outerIndex[outer]; - Index end = isCompressed() ? m_outerIndex[outer + 1] : m_outerIndex[outer] + m_innerNonZeros[outer]; + Index start = m_idx_data.outerIndex[outer]; + Index end = isCompressed() ? m_idx_data.outerIndex[outer + 1] + : m_idx_data.outerIndex[outer] + m_idx_data.innerNonZeros[outer]; eigen_assert(end >= start && "you probably called coeffRef on a non finalized matrix"); Index dst = start == end ? end : m_data.searchLowerIndex(start, end, inner); if (dst == end) { - Index capacity = m_outerIndex[outer + 1] - end; + Index capacity = m_idx_data.outerIndex[outer + 1] - end; if (capacity > 0) { // implies uncompressed: push to back of vector - m_innerNonZeros[outer]++; + m_idx_data.innerNonZeros[outer]++; m_data.index(end) = StorageIndex(inner); m_data.value(end) = Scalar(0); if (inserted != nullptr) { @@ -303,9 +307,9 @@ class SparseMatrix : public SparseCompressedBase + template ::type = 0> inline void reserveInnerVectors(const SizesType& reserveSizes) { if (isCompressed()) { Index totalReserveSize = 0; @@ -351,63 +355,69 @@ class SparseMatrix : public SparseCompressedBase(m_outerSize); + m_idx_data.innerNonZeros = internal::conditional_aligned_new_auto(m_outerSize); // temporarily use m_innerSizes to hold the new starting points. - StorageIndex* newOuterIndex = m_innerNonZeros; + StorageIndex* newOuterIndex = m_idx_data.innerNonZeros; Index count = 0; for (Index j = 0; j < m_outerSize; ++j) { newOuterIndex[j] = internal::convert_index(count); Index reserveSize = internal::convert_index(reserveSizes[j]); - count += reserveSize + internal::convert_index(m_outerIndex[j + 1] - m_outerIndex[j]); + count += reserveSize + internal::convert_index(m_idx_data.outerIndex[j + 1] - m_idx_data.outerIndex[j]); } m_data.reserve(totalReserveSize); - StorageIndex previousOuterIndex = m_outerIndex[m_outerSize]; + StorageIndex previousOuterIndex = m_idx_data.outerIndex[m_outerSize]; for (Index j = m_outerSize - 1; j >= 0; --j) { - StorageIndex innerNNZ = previousOuterIndex - m_outerIndex[j]; - StorageIndex begin = m_outerIndex[j]; + StorageIndex innerNNZ = previousOuterIndex - m_idx_data.outerIndex[j]; + StorageIndex begin = m_idx_data.outerIndex[j]; StorageIndex end = begin + innerNNZ; StorageIndex target = newOuterIndex[j]; internal::smart_memmove(innerIndexPtr() + begin, innerIndexPtr() + end, innerIndexPtr() + target); internal::smart_memmove(valuePtr() + begin, valuePtr() + end, valuePtr() + target); - previousOuterIndex = m_outerIndex[j]; - m_outerIndex[j] = newOuterIndex[j]; - m_innerNonZeros[j] = innerNNZ; + previousOuterIndex = m_idx_data.outerIndex[j]; + m_idx_data.outerIndex[j] = newOuterIndex[j]; + m_idx_data.innerNonZeros[j] = innerNNZ; } if (m_outerSize > 0) - m_outerIndex[m_outerSize] = m_outerIndex[m_outerSize - 1] + m_innerNonZeros[m_outerSize - 1] + - internal::convert_index(reserveSizes[m_outerSize - 1]); + m_idx_data.outerIndex[m_outerSize] = m_idx_data.outerIndex[m_outerSize - 1] + + m_idx_data.innerNonZeros[m_outerSize - 1] + + internal::convert_index(reserveSizes[m_outerSize - 1]); - m_data.resize(m_outerIndex[m_outerSize]); + m_data.resize(m_idx_data.outerIndex[m_outerSize]); } else { StorageIndex* newOuterIndex = internal::conditional_aligned_new_auto(m_outerSize + 1); Index count = 0; for (Index j = 0; j < m_outerSize; ++j) { newOuterIndex[j] = internal::convert_index(count); - Index alreadyReserved = - internal::convert_index(m_outerIndex[j + 1] - m_outerIndex[j] - m_innerNonZeros[j]); + Index alreadyReserved = internal::convert_index(m_idx_data.outerIndex[j + 1] - m_idx_data.outerIndex[j] - + m_idx_data.innerNonZeros[j]); Index reserveSize = internal::convert_index(reserveSizes[j]); Index toReserve = numext::maxi(reserveSize, alreadyReserved); - count += toReserve + internal::convert_index(m_innerNonZeros[j]); + count += toReserve + internal::convert_index(m_idx_data.innerNonZeros[j]); } newOuterIndex[m_outerSize] = internal::convert_index(count); m_data.resize(count); for (Index j = m_outerSize - 1; j >= 0; --j) { - StorageIndex innerNNZ = m_innerNonZeros[j]; - StorageIndex begin = m_outerIndex[j]; + StorageIndex innerNNZ = m_idx_data.innerNonZeros[j]; + StorageIndex begin = m_idx_data.outerIndex[j]; StorageIndex target = newOuterIndex[j]; m_data.moveChunk(begin, target, innerNNZ); } - std::swap(m_outerIndex, newOuterIndex); + std::swap(m_idx_data.outerIndex, newOuterIndex); internal::conditional_aligned_delete_auto(newOuterIndex, m_outerSize + 1); } } + template ::type = 0> + inline void reserveInnerVectors(const SizesType&) { + eigen_assert(false && "You cannot call dynamic size methods for static size matrices"); + } + public: //--- low level purely coherent filling --- @@ -428,11 +438,13 @@ class SparseMatrix : public SparseCompressedBase(m_data.size()); Index i = m_outerSize; // find the last filled column - while (i >= 0 && m_outerIndex[i] == 0) --i; + while (i >= 0 && m_idx_data.outerIndex[i] == 0) --i; ++i; while (i <= m_outerSize) { - m_outerIndex[i] = size; + m_idx_data.outerIndex[i] = size; ++i; } } @@ -484,21 +496,23 @@ class SparseMatrix : public SparseCompressedBase m_outerIndex[j]) uncompress(); + if (m_idx_data.outerIndex[j + num] > m_idx_data.outerIndex[j]) uncompress(); - // shift m_outerIndex and m_innerNonZeros [num] to the left - internal::smart_memmove(m_outerIndex + begin, m_outerIndex + end + 1, m_outerIndex + target); + // shift m_idx_data.outerIndex and m_idx_data.innerNonZeros [num] to the left + internal::smart_memmove(m_idx_data.outerIndex + begin, m_idx_data.outerIndex + end + 1, + m_idx_data.outerIndex + target); if (!isCompressed()) - internal::smart_memmove(m_innerNonZeros + begin, m_innerNonZeros + end, m_innerNonZeros + target); + internal::smart_memmove(m_idx_data.innerNonZeros + begin, m_idx_data.innerNonZeros + end, + m_idx_data.innerNonZeros + target); - // if m_outerIndex[0] > 0, shift the data within the first vector while it is easy to do so - if (m_outerIndex[0] > StorageIndex(0)) { + // if m_idx_data.outerIndex[0] > 0, shift the data within the first vector while it is easy to do so + if (m_idx_data.outerIndex[0] > StorageIndex(0)) { uncompress(); - const Index from = internal::convert_index(m_outerIndex[0]); + const Index from = internal::convert_index(m_idx_data.outerIndex[0]); const Index to = Index(0); - const Index chunkSize = internal::convert_index(m_innerNonZeros[0]); + const Index chunkSize = internal::convert_index(m_idx_data.innerNonZeros[0]); m_data.moveChunk(from, to, chunkSize); - m_outerIndex[0] = StorageIndex(0); + m_idx_data.outerIndex[0] = StorageIndex(0); } // truncate the matrix to the smaller size @@ -520,15 +534,17 @@ class SparseMatrix : public SparseCompressedBase void setFromSortedTriplets(const InputIterators& begin, const InputIterators& end); + template + void assignFromSortedTriplets(const InputIterators& begin, const InputIterators& end); + template void setFromSortedTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func); @@ -559,6 +578,9 @@ class SparseMatrix : public SparseCompressedBase void insertFromSortedTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func); + // template + // void setFromTripletsStatic(const InputIterators& begin, const InputIterators& end); + //--- /** \internal @@ -566,14 +588,14 @@ class SparseMatrix : public SparseCompressedBase= 0 && j < m_outerSize && "invalid outer index"); eigen_assert(i >= 0 && i < m_innerSize && "invalid inner index"); - Index start = m_outerIndex[j]; - Index end = isCompressed() ? m_outerIndex[j + 1] : start + m_innerNonZeros[j]; + Index start = m_idx_data.outerIndex[j]; + Index end = isCompressed() ? m_idx_data.outerIndex[j + 1] : start + m_idx_data.innerNonZeros[j]; Index dst = start == end ? end : m_data.searchLowerIndex(start, end, i); if (dst == end) { - Index capacity = m_outerIndex[j + 1] - end; + Index capacity = m_idx_data.outerIndex[j + 1] - end; if (capacity > 0) { // implies uncompressed: push to back of vector - m_innerNonZeros[j]++; + m_idx_data.innerNonZeros[j]++; m_data.index(end) = StorageIndex(i); m_data.value(end) = Scalar(0); return m_data.value(end); @@ -589,16 +611,16 @@ class SparseMatrix : public SparseCompressedBase 0); + eigen_internal_assert(m_idx_data.outerIndex != 0 && m_outerSize > 0); - StorageIndex start = m_outerIndex[1]; - m_outerIndex[1] = m_innerNonZeros[0]; + StorageIndex start = m_idx_data.outerIndex[1]; + m_idx_data.outerIndex[1] = m_idx_data.innerNonZeros[0]; // try to move fewer, larger contiguous chunks Index copyStart = start; - Index copyTarget = m_innerNonZeros[0]; + Index copyTarget = m_idx_data.innerNonZeros[0]; for (Index j = 1; j < m_outerSize; j++) { - StorageIndex end = start + m_innerNonZeros[j]; - StorageIndex nextStart = m_outerIndex[j + 1]; + StorageIndex end = start + m_idx_data.innerNonZeros[j]; + StorageIndex nextStart = m_idx_data.outerIndex[j + 1]; // dont forget to move the last chunk! bool breakUpCopy = (end != nextStart) || (j == m_outerSize - 1); if (breakUpCopy) { @@ -608,25 +630,26 @@ class SparseMatrix : public SparseCompressedBase(m_innerNonZeros, m_outerSize); - m_innerNonZeros = 0; + internal::conditional_aligned_delete_auto(m_idx_data.innerNonZeros, m_outerSize); + m_idx_data.innerNonZeros = 0; m_data.squeeze(); } /** Turns the matrix into the uncompressed mode */ void uncompress() { if (!isCompressed()) return; - m_innerNonZeros = internal::conditional_aligned_new_auto(m_outerSize); - if (m_outerIndex[m_outerSize] == 0) { + m_idx_data.innerNonZeros = internal::conditional_aligned_new_auto(m_outerSize); + if (m_idx_data.outerIndex[m_outerSize] == 0) { using std::fill_n; - fill_n(m_innerNonZeros, m_outerSize, StorageIndex(0)); + fill_n(m_idx_data.innerNonZeros, m_outerSize, StorageIndex(0)); } else { - for (Index j = 0; j < m_outerSize; j++) m_innerNonZeros[j] = m_outerIndex[j + 1] - m_outerIndex[j]; + for (Index j = 0; j < m_outerSize; j++) + m_idx_data.innerNonZeros[j] = m_idx_data.outerIndex[j + 1] - m_idx_data.outerIndex[j]; } } @@ -646,12 +669,12 @@ class SparseMatrix : public SparseCompressedBase(m_outerIndex, newOuterSize + 1, - m_outerSize + 1); + m_idx_data.outerIndex = internal::conditional_aligned_realloc_new_auto( + m_idx_data.outerIndex, newOuterSize + 1, m_outerSize + 1); if (!isCompressed()) - m_innerNonZeros = internal::conditional_aligned_realloc_new_auto(m_innerNonZeros, - newOuterSize, m_outerSize); + m_idx_data.innerNonZeros = internal::conditional_aligned_realloc_new_auto( + m_idx_data.innerNonZeros, newOuterSize, m_outerSize); if (outerChange > 0) { - StorageIndex lastIdx = m_outerSize == 0 ? StorageIndex(0) : m_outerIndex[m_outerSize]; + StorageIndex lastIdx = m_outerSize == 0 ? StorageIndex(0) : m_idx_data.outerIndex[m_outerSize]; using std::fill_n; - fill_n(m_outerIndex + m_outerSize, outerChange + 1, lastIdx); + fill_n(m_idx_data.outerIndex + m_outerSize, outerChange + 1, lastIdx); - if (!isCompressed()) fill_n(m_innerNonZeros + m_outerSize, outerChange, StorageIndex(0)); + if (!isCompressed()) fill_n(m_idx_data.innerNonZeros + m_outerSize, outerChange, StorageIndex(0)); } } m_outerSize = newOuterSize; if (innerChange < 0) { for (Index j = 0; j < m_outerSize; j++) { - Index start = m_outerIndex[j]; - Index end = isCompressed() ? m_outerIndex[j + 1] : start + m_innerNonZeros[j]; + Index start = m_idx_data.outerIndex[j]; + Index end = isCompressed() ? m_idx_data.outerIndex[j + 1] : start + m_idx_data.innerNonZeros[j]; Index lb = m_data.searchLowerIndex(start, end, newInnerSize); if (lb != end) { uncompress(); - m_innerNonZeros[j] = StorageIndex(lb - start); + m_idx_data.innerNonZeros[j] = StorageIndex(lb - start); } } } m_innerSize = newInnerSize; - Index newSize = m_outerIndex[m_outerSize]; + Index newSize = m_idx_data.outerIndex[m_outerSize]; eigen_assert(newSize <= m_data.size()); m_data.resize(newSize); } @@ -731,22 +754,29 @@ class SparseMatrix : public SparseCompressedBase::type = 0> void resize(Index rows, Index cols) { const Index outerSize = IsRowMajor ? rows : cols; m_innerSize = IsRowMajor ? cols : rows; m_data.clear(); - if ((m_outerIndex == 0) || (m_outerSize != outerSize)) { - m_outerIndex = internal::conditional_aligned_realloc_new_auto(m_outerIndex, outerSize + 1, - m_outerSize + 1); + if ((m_idx_data.outerIndex == 0) || (m_outerSize != outerSize)) { + m_idx_data.outerIndex = internal::conditional_aligned_realloc_new_auto( + m_idx_data.outerIndex, outerSize + 1, m_outerSize + 1); m_outerSize = outerSize; } - internal::conditional_aligned_delete_auto(m_innerNonZeros, m_outerSize); - m_innerNonZeros = 0; + internal::conditional_aligned_delete_auto(m_idx_data.innerNonZeros, m_outerSize); + m_idx_data.innerNonZeros = 0; using std::fill_n; - fill_n(m_outerIndex, m_outerSize + 1, StorageIndex(0)); + fill_n(m_idx_data.outerIndex, m_outerSize + 1, StorageIndex(0)); + } + + template ::type = 0> + void resize(Index, Index) { + eigen_assert(false && "You cannot call dynamic size methods for static size matrices"); } /** \internal @@ -763,17 +793,23 @@ class SparseMatrix : public SparseCompressedBase - inline SparseMatrix(const SparseMatrixBase& other) - : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0) { + inline SparseMatrix(const SparseMatrixBase& other) : m_outerSize(0), m_innerSize(0), m_idx_data() { EIGEN_STATIC_ASSERT( (internal::is_same::value), YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) @@ -791,7 +827,7 @@ class SparseMatrix : public SparseCompressedBase inline SparseMatrix(const SparseSelfAdjointView& other) - : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0) { + : m_outerSize(0), m_innerSize(0), m_idx_data() { Base::operator=(other); } @@ -804,15 +840,13 @@ class SparseMatrix : public SparseCompressedBase - SparseMatrix(const ReturnByValue& other) - : Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0) { + SparseMatrix(const ReturnByValue& other) : Base(), m_outerSize(0), m_innerSize(0), m_idx_data() { initAssignment(other); other.evalTo(*this); } @@ -820,7 +854,7 @@ class SparseMatrix : public SparseCompressedBase explicit SparseMatrix(const DiagonalBase& other) - : Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0) { + : Base(), m_outerSize(0), m_innerSize(0), m_idx_data() { *this = other.derived(); } @@ -828,10 +862,10 @@ class SparseMatrix : public SparseCompressedBase(m_innerNonZeros, m_outerSize); - m_innerNonZeros = 0; + internal::conditional_aligned_delete_auto(m_idx_data.innerNonZeros, m_outerSize); + m_idx_data.innerNonZeros = 0; m_data.resize(m_outerSize); // is it necessary to squeeze? m_data.squeeze(); - std::iota(m_outerIndex, m_outerIndex + m_outerSize + 1, StorageIndex(0)); + std::iota(m_idx_data.outerIndex, m_idx_data.outerIndex + m_outerSize + 1, StorageIndex(0)); std::iota(innerIndexPtr(), innerIndexPtr() + m_outerSize, StorageIndex(0)); using std::fill_n; fill_n(valuePtr(), m_outerSize, Scalar(1)); @@ -861,7 +895,8 @@ class SparseMatrix : public SparseCompressedBase(m_outerIndex, m_outerSize + 1); - internal::conditional_aligned_delete_auto(m_innerNonZeros, m_outerSize); + EIGEN_IF_CONSTEXPR(!isStatic) { + internal::conditional_aligned_delete_auto(m_idx_data.outerIndex, m_outerSize + 1); + internal::conditional_aligned_delete_auto(m_idx_data.innerNonZeros, m_outerSize); + } } /** Overloaded for performance */ @@ -942,8 +979,8 @@ class SparseMatrix : public SparseCompressedBase void initAssignment(const Other& other) { resize(other.rows(), other.cols()); - internal::conditional_aligned_delete_auto(m_innerNonZeros, m_outerSize); - m_innerNonZeros = 0; + internal::conditional_aligned_delete_auto(m_idx_data.innerNonZeros, m_outerSize); + m_idx_data.innerNonZeros = 0; } /** \internal @@ -975,9 +1012,9 @@ class SparseMatrix : public SparseCompressedBase(m_innerNonZeros, m_outerSize); - m_innerNonZeros = 0; + internal::conditional_aligned_delete_auto(m_idx_data.innerNonZeros, m_outerSize); + m_idx_data.innerNonZeros = 0; resizeNonZeros(n); ValueMap valueMap(valuePtr(), n); - std::iota(m_outerIndex, m_outerIndex + n + 1, StorageIndex(0)); + std::iota(m_idx_data.outerIndex, m_idx_data.outerIndex + n + 1, StorageIndex(0)); std::iota(innerIndexPtr(), innerIndexPtr() + n, StorageIndex(0)); valueMap.setZero(); internal::call_assignment_no_alias(valueMap, diagXpr, assignFunc); @@ -1031,9 +1068,9 @@ class SparseMatrix : public SparseCompressedBase 0) { m_data.resize(m_data.size() + shift); - Index copyEnd = isCompressed() ? m_outerIndex[m_outerSize] - : m_outerIndex[m_outerSize - 1] + m_innerNonZeros[m_outerSize - 1]; + Index copyEnd = isCompressed() + ? m_idx_data.outerIndex[m_outerSize] + : m_idx_data.outerIndex[m_outerSize - 1] + m_idx_data.innerNonZeros[m_outerSize - 1]; for (Index j = m_outerSize - 1; deferredInsertions > 0; j--) { - Index begin = m_outerIndex[j]; - Index end = isCompressed() ? m_outerIndex[j + 1] : begin + m_innerNonZeros[j]; - Index capacity = m_outerIndex[j + 1] - end; + Index begin = m_idx_data.outerIndex[j]; + Index end = isCompressed() ? m_idx_data.outerIndex[j + 1] : begin + m_idx_data.innerNonZeros[j]; + Index capacity = m_idx_data.outerIndex[j + 1] - end; bool doInsertion = insertionLocations(j) >= 0; bool breakUpCopy = doInsertion && (capacity > 0); @@ -1066,14 +1104,14 @@ class SparseMatrix : public SparseCompressedBase= 0` to skip inactive nonzeros in each vector // this reduces the total number of copied elements, but requires more moveChunk calls if (breakUpCopy) { - Index copyBegin = m_outerIndex[j + 1]; + Index copyBegin = m_idx_data.outerIndex[j + 1]; Index to = copyBegin + shift; Index chunkSize = copyEnd - copyBegin; m_data.moveChunk(copyBegin, to, chunkSize); copyEnd = end; } - m_outerIndex[j + 1] += shift; + m_idx_data.outerIndex[j + 1] += shift; if (doInsertion) { // if there is capacity, shift into the inactive nonzeros @@ -1086,7 +1124,7 @@ class SparseMatrix : public SparseCompressedBase +void assign_from_sorted_triplets(const InputIterator& begin, const InputIterator& end, SparseMatrixType& mat, + DupFunctor dup_func) { + constexpr bool IsRowMajor = SparseMatrixType::IsRowMajor; + using StorageIndex = typename SparseMatrixType::StorageIndex; + + if (begin == end) return; + + constexpr StorageIndex kEmptyIndexValue(-1); + + StorageIndex previous_i = kEmptyIndexValue; + StorageIndex previous_j = kEmptyIndexValue; + Index back = 0; + + typename SparseMatrixType::Scalar* values = mat.data().valuePtr(); + typename SparseMatrixType::StorageIndex* inner_indices = mat.data().indexPtr(); + typename SparseMatrixType::StorageIndex* outer_indices = mat.outerIndexPtr(); + + for (InputIterator it(begin); it != end; ++it) { + StorageIndex j = convert_index(IsRowMajor ? it->row() : it->col()); + StorageIndex i = convert_index(IsRowMajor ? it->col() : it->row()); + bool duplicate = (previous_j == j) && (previous_i == i); + + if (duplicate) { + values[back - 1] = dup_func(values[back - 1], it->value()); + } else { + // push triplets to back + typename SparseMatrixType::Scalar val = it->value(); + values[back] = val; + + inner_indices[back] = i; + outer_indices[j + 1]++; + + back++; + previous_j = j; + previous_i = i; + } + } + + for (Index j = 1; j < mat.outerSize() + 1; ++j) { + outer_indices[j] += outer_indices[j - 1]; + } + // eigen_assert(back == SparseMatrixType::MaxNz); +} + // Creates a compressed sparse matrix from a sorted range of triplets template void set_from_triplets_sorted(const InputIterator& begin, const InputIterator& end, SparseMatrixType& mat, @@ -1321,11 +1404,20 @@ void insert_from_triplets_sorted(const InputIterator& begin, const InputIterator * 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 +template template -void SparseMatrix::setFromTriplets(const InputIterators& begin, - const InputIterators& end) { - internal::set_from_triplets>( +void SparseMatrix::setFromTriplets(const InputIterators& begin, + const InputIterators& end) { + internal::set_from_triplets>( + begin, end, *this, internal::scalar_sum_op()); +} + +template +template +void SparseMatrix::assignFromSortedTriplets( + const InputIterators& begin, const InputIterators& end) { + internal::assign_from_sorted_triplets>( begin, end, *this, internal::scalar_sum_op()); } @@ -1338,23 +1430,25 @@ void SparseMatrix::setFromTriplets(const InputI * mat.setFromTriplets(triplets.begin(), triplets.end(), [] (const Scalar&,const Scalar &b) { return b; }); * \endcode */ -template +template template -void SparseMatrix::setFromTriplets(const InputIterators& begin, - const InputIterators& end, DupFunctor dup_func) { - internal::set_from_triplets, DupFunctor>( - begin, end, *this, dup_func); +void SparseMatrix::setFromTriplets(const InputIterators& begin, + const InputIterators& end, + DupFunctor dup_func) { + internal::set_from_triplets, + 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 +template template -void SparseMatrix::setFromSortedTriplets(const InputIterators& begin, - const InputIterators& end) { - internal::set_from_triplets_sorted>( +void SparseMatrix::setFromSortedTriplets( + const InputIterators& begin, const InputIterators& end) { + internal::set_from_triplets_sorted>( begin, end, *this, internal::scalar_sum_op()); } @@ -1367,12 +1461,12 @@ void SparseMatrix::setFromSortedTriplets(const * mat.setFromSortedTriplets(triplets.begin(), triplets.end(), [] (const Scalar&,const Scalar &b) { return b; }); * \endcode */ -template +template template -void SparseMatrix::setFromSortedTriplets(const InputIterators& begin, - const InputIterators& end, - DupFunctor dup_func) { - internal::set_from_triplets_sorted, DupFunctor>( +void SparseMatrix::setFromSortedTriplets( + const InputIterators& begin, const InputIterators& end, DupFunctor dup_func) { + internal::set_from_triplets_sorted, DupFunctor>( begin, end, *this, dup_func); } @@ -1414,11 +1508,11 @@ void SparseMatrix::setFromSortedTriplets(const * 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 +template template -void SparseMatrix::insertFromTriplets(const InputIterators& begin, - const InputIterators& end) { - internal::insert_from_triplets>( +void SparseMatrix::insertFromTriplets( + const InputIterators& begin, const InputIterators& end) { + internal::insert_from_triplets>( begin, end, *this, internal::scalar_sum_op()); } @@ -1431,23 +1525,24 @@ void SparseMatrix::insertFromTriplets(const Inp * mat.insertFromTriplets(triplets.begin(), triplets.end(), [] (const Scalar&,const Scalar &b) { return b; }); * \endcode */ -template +template template -void SparseMatrix::insertFromTriplets(const InputIterators& begin, - const InputIterators& end, DupFunctor dup_func) { - internal::insert_from_triplets, DupFunctor>( - begin, end, *this, dup_func); +void SparseMatrix::insertFromTriplets( + const InputIterators& begin, const InputIterators& end, DupFunctor dup_func) { + internal::insert_from_triplets, + 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 +template template -void SparseMatrix::insertFromSortedTriplets(const InputIterators& begin, - const InputIterators& end) { - internal::insert_from_triplets_sorted>( +void SparseMatrix::insertFromSortedTriplets( + const InputIterators& begin, const InputIterators& end) { + internal::insert_from_triplets_sorted>( begin, end, *this, internal::scalar_sum_op()); } @@ -1460,19 +1555,20 @@ void SparseMatrix::insertFromSortedTriplets(con * mat.insertFromSortedTriplets(triplets.begin(), triplets.end(), [] (const Scalar&,const Scalar &b) { return b; }); * \endcode */ -template +template template -void SparseMatrix::insertFromSortedTriplets(const InputIterators& begin, - const InputIterators& end, - DupFunctor dup_func) { - internal::insert_from_triplets_sorted, DupFunctor>( - begin, end, *this, dup_func); +void SparseMatrix::insertFromSortedTriplets( + const InputIterators& begin, const InputIterators& end, DupFunctor dup_func) { + internal::insert_from_triplets_sorted< + InputIterators, SparseMatrix, DupFunctor>(begin, end, + *this, dup_func); } /** \internal */ -template +template template -void SparseMatrix::collapseDuplicates(DenseBase& wi, DupFunctor dup_func) { +void SparseMatrix::collapseDuplicates(DenseBase& 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 @@ -1484,8 +1580,9 @@ void SparseMatrix::collapseDuplicates(DenseBas // 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) { 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) { + const StorageIndex end = + is_compressed ? m_idx_data.outerIndex[j + 1] : m_idx_data.outerIndex[j] + m_idx_data.innerNonZeros[j]; + for (StorageIndex k = m_idx_data.outerIndex[j]; k < end; ++k) { StorageIndex i = m_data.index(k); if (wi(i) >= newBegin) { // entry at k is a duplicate @@ -1500,21 +1597,22 @@ void SparseMatrix::collapseDuplicates(DenseBas ++count; } } - m_outerIndex[j] = newBegin; + m_idx_data.outerIndex[j] = newBegin; } - m_outerIndex[m_outerSize] = count; + m_idx_data.outerIndex[m_outerSize] = count; m_data.resize(count); // turn the matrix into compressed form (if it is not already) - internal::conditional_aligned_delete_auto(m_innerNonZeros, m_outerSize); - m_innerNonZeros = 0; + internal::conditional_aligned_delete_auto(m_idx_data.innerNonZeros, m_outerSize); + m_idx_data.innerNonZeros = 0; } /** \internal */ -template +template template -EIGEN_DONT_INLINE SparseMatrix& -SparseMatrix::operator=(const SparseMatrixBase& other) { +EIGEN_DONT_INLINE SparseMatrix& +SparseMatrix::operator=( + const SparseMatrixBase& other) { EIGEN_STATIC_ASSERT( (internal::is_same::value), YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) @@ -1541,23 +1639,24 @@ SparseMatrix::operator=(const SparseMatrixBase< OtherCopyEval otherCopyEval(otherCopy); SparseMatrix dest(other.rows(), other.cols()); - Eigen::Map(dest.m_outerIndex, dest.outerSize()).setZero(); + Eigen::Map(dest.m_idx_data.outerIndex, dest.outerSize()).setZero(); // pass 1 // FIXME the above copy could be merged with that pass for (Index j = 0; j < otherCopy.outerSize(); ++j) - for (typename OtherCopyEval::InnerIterator it(otherCopyEval, j); it; ++it) ++dest.m_outerIndex[it.index()]; + for (typename OtherCopyEval::InnerIterator it(otherCopyEval, j); it; ++it) + ++dest.m_idx_data.outerIndex[it.index()]; // prefix sum StorageIndex count = 0; IndexVector positions(dest.outerSize()); for (Index j = 0; j < dest.outerSize(); ++j) { - StorageIndex tmp = dest.m_outerIndex[j]; - dest.m_outerIndex[j] = count; + StorageIndex tmp = dest.m_idx_data.outerIndex[j]; + dest.m_idx_data.outerIndex[j] = count; positions[j] = count; count += tmp; } - dest.m_outerIndex[dest.outerSize()] = count; + dest.m_idx_data.outerIndex[dest.outerSize()] = count; // alloc dest.m_data.resize(count); // pass 2 @@ -1579,34 +1678,36 @@ SparseMatrix::operator=(const SparseMatrixBase< } } -template -inline typename SparseMatrix::Scalar& -SparseMatrix::insert(Index row, Index col) { +template +inline typename SparseMatrix::Scalar& +SparseMatrix::insert(Index row, Index col) { return insertByOuterInner(IsRowMajor ? row : col, IsRowMajor ? col : row); } -template -EIGEN_STRONG_INLINE typename SparseMatrix::Scalar& -SparseMatrix::insertAtByOuterInner(Index outer, Index inner, Index dst) { +template +EIGEN_STRONG_INLINE typename SparseMatrix::Scalar& +SparseMatrix::insertAtByOuterInner(Index outer, Index inner, + Index dst) { // random insertion into compressed matrix is very slow uncompress(); return insertUncompressedAtByOuterInner(outer, inner, dst); } -template -EIGEN_DEPRECATED EIGEN_DONT_INLINE typename SparseMatrix::Scalar& -SparseMatrix::insertUncompressed(Index row, Index col) { +template +EIGEN_DEPRECATED EIGEN_DONT_INLINE +typename SparseMatrix::Scalar& +SparseMatrix::insertUncompressed(Index row, Index col) { eigen_assert(!isCompressed()); Index outer = IsRowMajor ? row : col; Index inner = IsRowMajor ? col : row; - Index start = m_outerIndex[outer]; - Index end = start + m_innerNonZeros[outer]; + Index start = m_idx_data.outerIndex[outer]; + Index end = start + m_idx_data.innerNonZeros[outer]; Index dst = start == end ? end : m_data.searchLowerIndex(start, end, inner); if (dst == end) { - Index capacity = m_outerIndex[outer + 1] - end; + Index capacity = m_idx_data.outerIndex[outer + 1] - end; if (capacity > 0) { // implies uncompressed: push to back of vector - m_innerNonZeros[outer]++; + m_idx_data.innerNonZeros[outer]++; m_data.index(end) = StorageIndex(inner); m_data.value(end) = Scalar(0); return m_data.value(end); @@ -1617,23 +1718,26 @@ SparseMatrix::insertUncompressed(Index row, In return insertUncompressedAtByOuterInner(outer, inner, dst); } -template -EIGEN_DEPRECATED EIGEN_DONT_INLINE typename SparseMatrix::Scalar& -SparseMatrix::insertCompressed(Index row, Index col) { +template +EIGEN_DEPRECATED EIGEN_DONT_INLINE +typename SparseMatrix::Scalar& +SparseMatrix::insertCompressed(Index row, Index col) { eigen_assert(isCompressed()); Index outer = IsRowMajor ? row : col; Index inner = IsRowMajor ? col : row; - Index start = m_outerIndex[outer]; - Index end = m_outerIndex[outer + 1]; + Index start = m_idx_data.outerIndex[outer]; + Index end = m_idx_data.outerIndex[outer + 1]; Index dst = start == end ? end : m_data.searchLowerIndex(start, end, inner); eigen_assert((dst == end || m_data.index(dst) != inner) && "you cannot insert an element that already exists, you must call coeffRef to this end"); return insertCompressedAtByOuterInner(outer, inner, dst); } -template -typename SparseMatrix::Scalar& -SparseMatrix::insertCompressedAtByOuterInner(Index outer, Index inner, Index dst) { +template +typename SparseMatrix::Scalar& +SparseMatrix::insertCompressedAtByOuterInner(Index outer, + Index inner, + Index dst) { eigen_assert(isCompressed()); // compressed insertion always requires expanding the buffer // first, check if there is adequate allocated memory @@ -1645,12 +1749,12 @@ SparseMatrix::insertCompressedAtByOuterInner(I m_data.reserve(reserveSize); } m_data.resize(m_data.size() + 1); - Index chunkSize = m_outerIndex[m_outerSize] - dst; + Index chunkSize = m_idx_data.outerIndex[m_outerSize] - dst; // shift the existing data to the right if necessary m_data.moveChunk(dst, dst + 1, chunkSize); // update nonzero counts // potentially O(outerSize) bottleneck! - for (Index j = outer; j < m_outerSize; j++) m_outerIndex[j + 1]++; + for (Index j = outer; j < m_outerSize; j++) m_idx_data.outerIndex[j + 1]++; // initialize the coefficient m_data.index(dst) = StorageIndex(inner); m_data.value(dst) = Scalar(0); @@ -1658,23 +1762,25 @@ SparseMatrix::insertCompressedAtByOuterInner(I return m_data.value(dst); } -template -typename SparseMatrix::Scalar& -SparseMatrix::insertUncompressedAtByOuterInner(Index outer, Index inner, Index dst) { +template +typename SparseMatrix::Scalar& +SparseMatrix::insertUncompressedAtByOuterInner(Index outer, + Index inner, + Index dst) { eigen_assert(!isCompressed()); // find a vector with capacity, starting at `outer` and searching to the left and right for (Index leftTarget = outer - 1, rightTarget = outer; (leftTarget >= 0) || (rightTarget < m_outerSize);) { if (rightTarget < m_outerSize) { - Index start = m_outerIndex[rightTarget]; - Index end = start + m_innerNonZeros[rightTarget]; - Index nextStart = m_outerIndex[rightTarget + 1]; + Index start = m_idx_data.outerIndex[rightTarget]; + Index end = start + m_idx_data.innerNonZeros[rightTarget]; + Index nextStart = m_idx_data.outerIndex[rightTarget + 1]; Index capacity = nextStart - end; if (capacity > 0) { // move [dst, end) to dst+1 and insert at dst Index chunkSize = end - dst; if (chunkSize > 0) m_data.moveChunk(dst, dst + 1, chunkSize); - m_innerNonZeros[outer]++; - for (Index j = outer; j < rightTarget; j++) m_outerIndex[j + 1]++; + m_idx_data.innerNonZeros[outer]++; + for (Index j = outer; j < rightTarget; j++) m_idx_data.outerIndex[j + 1]++; m_data.index(dst) = StorageIndex(inner); m_data.value(dst) = Scalar(0); return m_data.value(dst); @@ -1682,17 +1788,17 @@ SparseMatrix::insertUncompressedAtByOuterInner rightTarget++; } if (leftTarget >= 0) { - Index start = m_outerIndex[leftTarget]; - Index end = start + m_innerNonZeros[leftTarget]; - Index nextStart = m_outerIndex[leftTarget + 1]; + Index start = m_idx_data.outerIndex[leftTarget]; + Index end = start + m_idx_data.innerNonZeros[leftTarget]; + Index nextStart = m_idx_data.outerIndex[leftTarget + 1]; Index capacity = nextStart - end; if (capacity > 0) { // tricky: dst is a lower bound, so we must insert at dst-1 when shifting left // move [nextStart, dst) to nextStart-1 and insert at dst-1 Index chunkSize = dst - nextStart; if (chunkSize > 0) m_data.moveChunk(nextStart, nextStart - 1, chunkSize); - m_innerNonZeros[outer]++; - for (Index j = leftTarget; j < outer; j++) m_outerIndex[j + 1]--; + m_idx_data.innerNonZeros[outer]++; + for (Index j = leftTarget; j < outer; j++) m_idx_data.outerIndex[j + 1]--; m_data.index(dst - 1) = StorageIndex(inner); m_data.value(dst - 1) = Scalar(0); return m_data.value(dst - 1); @@ -1704,12 +1810,12 @@ SparseMatrix::insertUncompressedAtByOuterInner // no room for interior insertion // nonZeros() == m_data.size() // record offset as outerIndxPtr will change - Index dst_offset = dst - m_outerIndex[outer]; + Index dst_offset = dst - m_idx_data.outerIndex[outer]; // allocate space for random insertion if (m_data.allocatedSize() == 0) { // fast method to allocate space for one element per vector in empty matrix m_data.resize(m_outerSize); - std::iota(m_outerIndex, m_outerIndex + m_outerSize + 1, StorageIndex(0)); + std::iota(m_idx_data.outerIndex, m_idx_data.outerIndex + m_outerSize + 1, StorageIndex(0)); } else { // check for integer overflow: if maxReserveSize == 0, insertion is not possible Index maxReserveSize = static_cast(NumTraits::highest()) - m_data.allocatedSize(); @@ -1727,12 +1833,12 @@ SparseMatrix::insertUncompressedAtByOuterInner } } // insert element at `dst` with new outer indices - Index start = m_outerIndex[outer]; - Index end = start + m_innerNonZeros[outer]; + Index start = m_idx_data.outerIndex[outer]; + Index end = start + m_idx_data.innerNonZeros[outer]; Index new_dst = start + dst_offset; Index chunkSize = end - new_dst; if (chunkSize > 0) m_data.moveChunk(new_dst, new_dst + 1, chunkSize); - m_innerNonZeros[outer]++; + m_idx_data.innerNonZeros[outer]++; m_data.index(new_dst) = StorageIndex(inner); m_data.value(new_dst) = Scalar(0); return m_data.value(new_dst); @@ -1740,11 +1846,11 @@ SparseMatrix::insertUncompressedAtByOuterInner namespace internal { -template -struct evaluator> - : evaluator>> { - typedef evaluator>> Base; - typedef SparseMatrix SparseMatrixType; +template +struct evaluator> + : evaluator>> { + typedef evaluator>> Base; + typedef SparseMatrix SparseMatrixType; evaluator() : Base() {} explicit evaluator(const SparseMatrixType& mat) : Base(mat) {} }; diff --git a/Eigen/src/SparseCore/SparseProduct.h b/Eigen/src/SparseCore/SparseProduct.h index 249dabc14588fe33190b039ec034e32461302f2f..2d9d3ad52eb6f70623ac10221803076065074fb1 100644 --- a/Eigen/src/SparseCore/SparseProduct.h +++ b/Eigen/src/SparseCore/SparseProduct.h @@ -162,9 +162,10 @@ struct unary_evaluator >, IteratorBased> } // end namespace internal // sparse matrix = sparse-product (can be sparse*sparse, sparse*perm, etc.) -template +template template -SparseMatrix& SparseMatrix::operator=( +SparseMatrix& +SparseMatrix::operator=( const Product& src) { // std::cout << "in Assignment : " << DstOptions << "\n"; SparseMatrix dst(src.rows(), src.cols()); diff --git a/Eigen/src/SparseCore/SparseRedux.h b/Eigen/src/SparseCore/SparseRedux.h index 732e4f77a7c213e1504d28c24431441e3c97c96b..92fe68259ec2c975959b49f1000a9c723e6959db 100644 --- a/Eigen/src/SparseCore/SparseRedux.h +++ b/Eigen/src/SparseCore/SparseRedux.h @@ -25,9 +25,9 @@ typename internal::traits::Scalar SparseMatrixBase::sum() cons return res; } -template -typename internal::traits >::Scalar -SparseMatrix::sum() const { +template +typename internal::traits >::Scalar +SparseMatrix::sum() const { eigen_assert(rows() > 0 && cols() > 0 && "you are using a non initialized matrix"); if (this->isCompressed()) return Matrix::Map(m_data.valuePtr(), m_data.size()).sum(); @@ -35,9 +35,9 @@ SparseMatrix::sum() const { return Base::sum(); } -template -typename internal::traits >::Scalar -SparseVector::sum() const { +template +typename internal::traits >::Scalar +SparseVector::sum() const { eigen_assert(rows() > 0 && cols() > 0 && "you are using a non initialized matrix"); return Matrix::Map(m_data.valuePtr(), m_data.size()).sum(); } diff --git a/Eigen/src/SparseCore/SparseUtil.h b/Eigen/src/SparseCore/SparseUtil.h index 33cedaf3eb3b2089694194834d64cfdec0135f7a..e02f4c99cd5d72ccd98eaac0013f59bd32ab2ab4 100644 --- a/Eigen/src/SparseCore/SparseUtil.h +++ b/Eigen/src/SparseCore/SparseUtil.h @@ -43,9 +43,11 @@ const int InnerRandomAccessPattern = 0x2 | CoherentAccessPattern; const int OuterRandomAccessPattern = 0x4 | CoherentAccessPattern; const int RandomAccessPattern = 0x8 | OuterRandomAccessPattern | InnerRandomAccessPattern; -template +template class SparseMatrix; -template +template class SparseVector; template @@ -91,18 +93,30 @@ template struct sparse_eval { typedef typename traits::Scalar Scalar_; typedef typename traits::StorageIndex StorageIndex_; + enum { + Rows_ = traits::RowsAtCompileTime, + Cols_ = traits::ColsAtCompileTime, + // MaxNZ_ = traits::MaxNZ + }; public: - typedef SparseVector type; + // TODO Find a clean way to forward MaxNZ (enum does not work if T is e.g. Product<...>) + typedef SparseVector type; }; template struct sparse_eval { typedef typename traits::Scalar Scalar_; typedef typename traits::StorageIndex StorageIndex_; + enum { + Rows_ = traits::RowsAtCompileTime, + Cols_ = traits::ColsAtCompileTime, + // MaxNZ_ = traits::MaxNZ + }; public: - typedef SparseVector type; + // TODO Find a clean way to forward MaxNZ (enum does not work if T is e.g. Product<...>) + typedef SparseVector type; }; // TODO this seems almost identical to plain_matrix_type @@ -112,8 +126,15 @@ struct sparse_eval { typedef typename traits::StorageIndex StorageIndex_; enum { Options_ = ((Flags & RowMajorBit) == RowMajorBit) ? RowMajor : ColMajor }; + enum { + Rows_ = traits::RowsAtCompileTime, + Cols_ = traits::ColsAtCompileTime, + // MaxNZ_ = traits::MaxNZ + }; + public: - typedef SparseMatrix type; + // TODO Find a clean way to forward MaxNZ (enum does not work if T is e.g. Product<...>) + typedef SparseMatrix type; }; template @@ -128,10 +149,12 @@ template struct plain_matrix_type { typedef typename traits::Scalar Scalar_; typedef typename traits::StorageIndex StorageIndex_; - enum { Options_ = ((evaluator::Flags & RowMajorBit) == RowMajorBit) ? RowMajor : ColMajor }; + enum { + Options_ = ((evaluator::Flags & RowMajorBit) == RowMajorBit) ? RowMajor : ColMajor, + }; public: - typedef SparseMatrix type; + typedef SparseMatrix type; }; template diff --git a/Eigen/src/SparseCore/SparseVector.h b/Eigen/src/SparseCore/SparseVector.h index 3f72a34dabd9510a72b45e63a618ae22c769bf10..193c84478e4d59bc15e3fcc764fb7de6c1a4601c 100644 --- a/Eigen/src/SparseCore/SparseVector.h +++ b/Eigen/src/SparseCore/SparseVector.h @@ -29,8 +29,8 @@ namespace Eigen { */ namespace internal { -template -struct traits > { +template +struct traits> { typedef Scalar_ Scalar; typedef StorageIndex_ StorageIndex; typedef Sparse StorageKind; @@ -38,10 +38,12 @@ struct traits > { enum { IsColVector = (Options_ & RowMajorBit) ? 0 : 1, - RowsAtCompileTime = IsColVector ? Dynamic : 1, - ColsAtCompileTime = IsColVector ? 1 : Dynamic, + IsStatic = MaxNZ_ != Dynamic, + RowsAtCompileTime = IsColVector ? Rows_ : 1, + ColsAtCompileTime = IsColVector ? 1 : Cols_, MaxRowsAtCompileTime = RowsAtCompileTime, MaxColsAtCompileTime = ColsAtCompileTime, + MaxNZ = MaxNZ_, Flags = Options_ | NestByRefBit | LvalueBit | (IsColVector ? 0 : RowMajorBit) | CompressedAccessBit, SupportedAccessPatterns = InnerRandomAccessPattern }; @@ -58,8 +60,8 @@ struct sparse_vector_assign_selector; } // namespace internal -template -class SparseVector : public SparseCompressedBase > { +template +class SparseVector : public SparseCompressedBase> { typedef SparseCompressedBase Base; using Base::convert_index; @@ -68,8 +70,14 @@ class SparseVector : public SparseCompressedBase Storage; - enum { IsColVector = internal::traits::IsColVector }; + template + void setFromSortedPairs(InputIterators begin, InputIterators end); + + typedef internal::CompressedStorage Storage; + enum { + IsColVector = internal::traits::IsColVector, + IsStatic = internal::traits::IsStatic + }; enum { Options = Options_ }; @@ -254,7 +262,12 @@ class SparseVector : public SparseCompressedBase +void set_from_sorted_pairs(const InputIterator& begin, const InputIterator& end, SparseVectorType& vec) { + using StorageIndex = typename SparseVectorType::StorageIndex; + using Scalar = typename SparseVectorType::Scalar; + Scalar* values = vec.valuePtr(); + StorageIndex* inner_indices = vec.innerIndexPtr(); + + Eigen::Index i = 0; + for (InputIterator it(begin); it != end; ++it, ++i) { + inner_indices[i] = it->first; + values[i] = it->second; + } +} +} // namespace internal + +template +template +void SparseVector::setFromSortedPairs(InputIterator begin, + InputIterator end) { + internal::set_from_sorted_pairs>( + begin, end, *this); +} + namespace internal { -template -struct evaluator > : evaluator_base > { - typedef SparseVector SparseVectorType; +template +struct evaluator> + : evaluator_base> { + typedef SparseVector SparseVectorType; typedef evaluator_base Base; typedef typename SparseVectorType::InnerIterator InnerIterator; typedef typename SparseVectorType::ReverseInnerIterator ReverseInnerIterator; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e62ec4529a4ddabdbf8c649e698d339ac708f342..0b9dd29171adbbedf530b0d74a65f3eebcee19d4 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -277,8 +277,10 @@ ei_add_test(stdlist_overload) ei_add_test(stddeque) ei_add_test(stddeque_overload) ei_add_test(sparse_basic) +ei_add_test(sparse_nomalloc) ei_add_test(sparse_block) ei_add_test(sparse_vector) +ei_add_test(sparse_vector_nomalloc) ei_add_test(sparse_product) ei_add_test(sparse_ref) ei_add_test(sparse_solvers) diff --git a/test/sparse_nomalloc.cpp b/test/sparse_nomalloc.cpp new file mode 100644 index 0000000000000000000000000000000000000000..2e00d18256b3e5292c8bd5bce9a41c274665b79f --- /dev/null +++ b/test/sparse_nomalloc.cpp @@ -0,0 +1,127 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#define EIGEN_NO_MALLOC + +#ifndef EIGEN_SPARSE_TEST_INCLUDED_FROM_SPARSE_EXTRA +static long g_dense_op_sparse_count = 0; +#endif + +#include "sparse.h" + +template +void sparse_fixed_nnz_test(const SparseMatrixType&) { + typedef typename SparseMatrixType::Scalar Scalar; + + Scalar eps = Scalar(1e-6); + + // maxNZ should not be Dynamic + VERIFY(internal::traits::MaxNZ != Dynamic); + + // Create a sparse matrix with dimensions matching the reference + SparseMatrixType m; + + // maximum number of non-zeros + const int max_nz = internal::traits::MaxNZ; + + std::vector> triplets(max_nz); + // Generate non-duplicate triplets + for (int i_count = 0; i_count < max_nz; ++i_count) { + int i = (i_count * 2 + 1) % SparseMatrixType::RowsAtCompileTime; + int j = (i_count * 2 + 1) / SparseMatrixType::RowsAtCompileTime; + Scalar tmp = Scalar(i * 10 + j + 1 + 3 * i * j); + triplets[i_count] = Triplet(i, j, tmp); + } + + // Triplets must be sorted + std::sort(triplets.begin(), triplets.end(), [&](const Triplet& a, const Triplet& b) { + if (SparseMatrixType::IsRowMajor) { + if (a.row() != b.row()) { + return a.row() < b.row(); + } + return a.col() < b.col(); + } else { + if (a.col() != b.col()) { + return a.col() < b.col(); + } + return a.row() < b.row(); + } + }); + m.assignFromSortedTriplets(triplets.begin(), triplets.end()); + + typedef Matrix + DenseMatrixType; + DenseMatrixType a; + a.setZero(); + for (size_t i = 0; i < triplets.size(); ++i) { + a(triplets[i].row(), triplets[i].col()) = triplets[i].value(); + } + VERIFY_IS_MUCH_SMALLER_THAN((a - m).norm(), eps); + VERIFY_IS_APPROX(m, a); +} + +EIGEN_DECLARE_TEST(sparse_nomalloc) { + g_dense_op_sparse_count = 0; // Suppresses compiler warning + for (int i = 0; i < g_repeat; i++) { + // Square matrices, ColMajor, + { + CALL_SUBTEST_1((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_1((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_1((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_1((sparse_fixed_nnz_test(SparseMatrix()))); + } + + // Square matrices, RowMajor + { + CALL_SUBTEST_2((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_2((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_2((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_2((sparse_fixed_nnz_test(SparseMatrix()))); + } + + // complex numbers, 8 non-zeros + { + CALL_SUBTEST_3((sparse_fixed_nnz_test(SparseMatrix, RowMajor, int, 10, 10, 8>()))); + CALL_SUBTEST_3((sparse_fixed_nnz_test(SparseMatrix, RowMajor, int, 2, 16, 8>()))); + CALL_SUBTEST_3((sparse_fixed_nnz_test(SparseMatrix, RowMajor, int, 8, 3, 8>()))); + CALL_SUBTEST_3((sparse_fixed_nnz_test(SparseMatrix, RowMajor, int, 10, 2, 8>()))); + CALL_SUBTEST_3((sparse_fixed_nnz_test(SparseMatrix, ColMajor, int, 10, 10, 8>()))); + CALL_SUBTEST_3((sparse_fixed_nnz_test(SparseMatrix, ColMajor, int, 2, 16, 8>()))); + CALL_SUBTEST_3((sparse_fixed_nnz_test(SparseMatrix, ColMajor, int, 8, 3, 8>()))); + CALL_SUBTEST_3((sparse_fixed_nnz_test(SparseMatrix, ColMajor, int, 10, 2, 8>()))); + } + + // Test with non-square matrices, double + { + CALL_SUBTEST_4((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_4((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_4((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_4((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_4((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_4((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_4((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_4((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_4((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_4((sparse_fixed_nnz_test(SparseMatrix()))); + } + + // Test with non-square matrices, float + { + CALL_SUBTEST_5((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_5((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_5((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_5((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_5((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_5((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_5((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_5((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_5((sparse_fixed_nnz_test(SparseMatrix()))); + CALL_SUBTEST_5((sparse_fixed_nnz_test(SparseMatrix()))); + } + } +} \ No newline at end of file diff --git a/test/sparse_vector_nomalloc.cpp b/test/sparse_vector_nomalloc.cpp new file mode 100644 index 0000000000000000000000000000000000000000..19d72e580fd7a88be42209a47333cfaef69d9952 --- /dev/null +++ b/test/sparse_vector_nomalloc.cpp @@ -0,0 +1,90 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#define EIGEN_NO_MALLOC + +#ifndef EIGEN_SPARSE_TEST_INCLUDED_FROM_SPARSE_EXTRA +static long g_dense_op_sparse_count = 0; +#endif + +#include "sparse.h" + +template +void sparsevector_fixed_nnz_test(const SparseVectorType&) { + typedef typename SparseVectorType::Scalar Scalar; + + Scalar eps = Scalar(1e-6); + + // maxNZ should not be Dynamic + VERIFY(internal::traits::MaxNZ != Dynamic); + + // Create a sparse matrix with dimensions matching the reference + SparseVectorType vec; + using StorageIndex = typename SparseVectorType::StorageIndex; + + // maximum number of non-zeros + const int max_nz = internal::traits::MaxNZ; + + std::vector> pairs(max_nz); + // Generate non-duplicate triplets + for (int i_count = 0; i_count < max_nz; ++i_count) { + int i = (i_count * 2 + 1); + Scalar tmp = Scalar(i * 10 + 1); + pairs[i_count] = std::make_pair(i, tmp); + } + // Pairs must be sorted + vec.setFromSortedPairs(pairs.begin(), pairs.end()); + + typedef Matrix + DenseMatrixType; + DenseMatrixType vec_full; + + const Scalar* values = vec.valuePtr(); + const StorageIndex* indices = vec.innerIndexPtr(); + + vec_full.setZero(); + for (size_t i = 0; i < max_nz; ++i) { + vec_full(indices[i]) = values[i]; + } + VERIFY_IS_MUCH_SMALLER_THAN((vec - vec_full).norm(), eps); + VERIFY_IS_APPROX(vec, vec_full); +} + +EIGEN_DECLARE_TEST(sparse_vector_nomalloc) { + g_dense_op_sparse_count = 0; // Suppresses compiler warning + for (int i = 0; i < g_repeat; i++) { + // ColMajor, + { + CALL_SUBTEST_1((sparsevector_fixed_nnz_test(SparseVector()))); + CALL_SUBTEST_1((sparsevector_fixed_nnz_test(SparseVector()))); + CALL_SUBTEST_1((sparsevector_fixed_nnz_test(SparseVector()))); + CALL_SUBTEST_1((sparsevector_fixed_nnz_test(SparseVector()))); + CALL_SUBTEST_1((sparsevector_fixed_nnz_test(SparseVector()))); + } + + // RowMajor + { + CALL_SUBTEST_2((sparsevector_fixed_nnz_test(SparseVector()))); + CALL_SUBTEST_2((sparsevector_fixed_nnz_test(SparseVector()))); + CALL_SUBTEST_2((sparsevector_fixed_nnz_test(SparseVector()))); + CALL_SUBTEST_2((sparsevector_fixed_nnz_test(SparseVector()))); + } + + // Complex values + { + CALL_SUBTEST_3((sparsevector_fixed_nnz_test(SparseVector, RowMajor, int, 1, 21, 8>()))); + CALL_SUBTEST_3((sparsevector_fixed_nnz_test(SparseVector, RowMajor, int, 1, 22, 8>()))); + CALL_SUBTEST_3((sparsevector_fixed_nnz_test(SparseVector, RowMajor, int, 1, 60, 8>()))); + CALL_SUBTEST_3((sparsevector_fixed_nnz_test(SparseVector, RowMajor, int, 1, 41, 8>()))); + CALL_SUBTEST_3((sparsevector_fixed_nnz_test(SparseVector, ColMajor, int, 21, 1, 8>()))); + CALL_SUBTEST_3((sparsevector_fixed_nnz_test(SparseVector, ColMajor, int, 22, 1, 8>()))); + CALL_SUBTEST_3((sparsevector_fixed_nnz_test(SparseVector, ColMajor, int, 60, 1, 8>()))); + CALL_SUBTEST_3((sparsevector_fixed_nnz_test(SparseVector, ColMajor, int, 41, 1, 8>()))); + } + } +} \ No newline at end of file diff --git a/unsupported/Eigen/src/SparseExtra/BlockSparseMatrix.h b/unsupported/Eigen/src/SparseExtra/BlockSparseMatrix.h index 6e8be8490d5b1c8cb18c695547c728668916e192..081b52294d6fc2fe59f75950be8488b38029396a 100644 --- a/unsupported/Eigen/src/SparseExtra/BlockSparseMatrix.h +++ b/unsupported/Eigen/src/SparseExtra/BlockSparseMatrix.h @@ -286,7 +286,7 @@ class BlockSparseMatrix m_values(0), m_blockPtr(0), m_indices(0), - m_outerIndex(0), + m_idx_data.outerIndex(0), m_blockSize(BlockSize) {} /** @@ -302,7 +302,7 @@ class BlockSparseMatrix m_values(0), m_blockPtr(0), m_indices(0), - m_outerIndex(0), + m_idx_data.outerIndex(0), m_blockSize(BlockSize) {} /** @@ -325,7 +325,7 @@ class BlockSparseMatrix if (m_blockSize != Dynamic) std::copy(other.m_blockPtr, other.m_blockPtr + m_nonzerosblocks, m_blockPtr); std::copy(other.m_indices, other.m_indices + m_nonzerosblocks, m_indices); - std::copy(other.m_outerIndex, other.m_outerIndex + m_outerBSize, m_outerIndex); + std::copy(other.m_idx_data.outerIndex, other.m_idx_data.outerIndex + m_outerBSize, m_idx_data.outerIndex); } friend void swap(BlockSparseMatrix& first, BlockSparseMatrix& second) { @@ -338,7 +338,7 @@ class BlockSparseMatrix std::swap(first.m_values, second.m_values); std::swap(first.m_blockPtr, second.m_blockPtr); std::swap(first.m_indices, second.m_indices); - std::swap(first.m_outerIndex, second.m_outerIndex); + std::swap(first.m_idx_data.outerIndex, second.m_idx_data.outerIndex); std::swap(first.m_BlockSize, second.m_blockSize); } @@ -350,7 +350,7 @@ class BlockSparseMatrix // Destructor ~BlockSparseMatrix() { - delete[] m_outerIndex; + delete[] m_idx_data.outerIndex; delete[] m_innerOffset; delete[] m_outerOffset; delete[] m_indices; @@ -420,24 +420,24 @@ class BlockSparseMatrix StorageIndex idx = 0; // Position of this block in the column block StorageIndex bi = innerToBlock(it_spmat.index()); // Index of the current nonzero block // Go to the inner block where this element belongs to - while (bi > m_indices[m_outerIndex[bj] + idx]) ++idx; // Not expensive for ordered blocks + while (bi > m_indices[m_idx_data.outerIndex[bj] + idx]) ++idx; // Not expensive for ordered blocks StorageIndex idxVal; // Get the right position in the array of values for this element if (m_blockSize == Dynamic) { // Offset from all blocks before ... - idxVal = m_blockPtr[m_outerIndex[bj] + idx]; + idxVal = m_blockPtr[m_idx_data.outerIndex[bj] + idx]; // ... and offset inside the block idxVal += (j - blockOuterIndex(bj)) * blockOuterSize(bj) + it_spmat.index() - m_innerOffset[bi]; } else { // All blocks before - idxVal = (m_outerIndex[bj] + idx) * m_blockSize * m_blockSize; + idxVal = (m_idx_data.outerIndex[bj] + idx) * m_blockSize * m_blockSize; // inside the block idxVal += (j - blockOuterIndex(bj)) * m_blockSize + (it_spmat.index() % m_blockSize); } // Insert the value m_values[idxVal] = it_spmat.value(); } // end of this column - } // end of this block - } // end of this outer block + } // end of this block + } // end of this outer block return *this; } @@ -465,7 +465,7 @@ class BlockSparseMatrix reserve(blockPattern.nonZeros()); // Browse the block pattern and set up the various pointers - m_outerIndex[0] = 0; + m_idx_data.outerIndex[0] = 0; if (m_blockSize == Dynamic) m_blockPtr[0] = 0; for (StorageIndex nz = 0; nz < m_nonzeros; ++nz) m_values[nz] = Scalar(0); for (StorageIndex bj = 0; bj < m_outerBSize; ++bj) { @@ -482,14 +482,14 @@ class BlockSparseMatrix // Now, fill block indices and (eventually) pointers to blocks for (StorageIndex idx = 0; idx < nzBlockIdx.size(); ++idx) { - StorageIndex offset = m_outerIndex[bj] + idx; // offset in m_indices + StorageIndex offset = m_idx_data.outerIndex[bj] + idx; // offset in m_indices m_indices[offset] = nzBlockIdx[idx]; if (m_blockSize == Dynamic) m_blockPtr[offset] = m_blockPtr[offset - 1] + blockInnerSize(nzBlockIdx[idx]) * blockOuterSize(bj); // There is no blockPtr for fixed-size blocks... not needed !??? } // Save the pointer to the next outer block - m_outerIndex[bj + 1] = m_outerIndex[bj] + nzBlockIdx.size(); + m_idx_data.outerIndex[bj + 1] = m_idx_data.outerIndex[bj] + nzBlockIdx.size(); } } @@ -552,7 +552,7 @@ class BlockSparseMatrix "TRYING TO RESERVE ZERO-SIZE MATRICES, CALL resize() first"); // FIXME Should free if already allocated - m_outerIndex = new StorageIndex[m_outerBSize + 1]; + m_idx_data.outerIndex = new StorageIndex[m_outerBSize + 1]; m_nonzerosblocks = nonzerosblocks; if (m_blockSize != Dynamic) { @@ -623,13 +623,13 @@ class BlockSparseMatrix VectorXi block_id(m_outerBSize); // To be used as a block marker during insertion // Setup outer index pointers and markers - m_outerIndex[0] = 0; + m_idx_data.outerIndex[0] = 0; if (m_blockSize == Dynamic) m_blockPtr[0] = 0; for (StorageIndex bj = 0; bj < m_outerBSize; ++bj) { - m_outerIndex[bj + 1] = m_outerIndex[bj] + nzblock_outer(bj); - block_id(bj) = m_outerIndex[bj]; + m_idx_data.outerIndex[bj + 1] = m_idx_data.outerIndex[bj] + nzblock_outer(bj); + block_id(bj) = m_idx_data.outerIndex[bj]; if (m_blockSize == Dynamic) { - m_blockPtr[m_outerIndex[bj + 1]] = m_blockPtr[m_outerIndex[bj]] + nz_outer(bj); + m_blockPtr[m_idx_data.outerIndex[bj + 1]] = m_blockPtr[m_idx_data.outerIndex[bj]] + nz_outer(bj); } } @@ -656,9 +656,9 @@ class BlockSparseMatrix // while (idvalue().rows()*it->value().cols(); // m_blockPtr[id_nzblock+1] = m_blockPtr[id_nzblock] + block_size; @@ -669,7 +669,7 @@ class BlockSparseMatrix // while(id < m_outerBSize-1) // Empty columns at the end // { // id++; - // m_outerIndex[id+1]=m_outerIndex[id]; + // m_idx_data.outerIndex[id+1]=m_idx_data.outerIndex[id]; // } // } } @@ -743,8 +743,8 @@ class BlockSparseMatrix StorageIndex csize = IsColMajor ? blockOuterSize(bcol) : blockInnerSize(brow); StorageIndex inner = IsColMajor ? brow : bcol; StorageIndex outer = IsColMajor ? bcol : brow; - StorageIndex offset = m_outerIndex[outer]; - while (offset < m_outerIndex[outer + 1] && m_indices[offset] != inner) offset++; + StorageIndex offset = m_idx_data.outerIndex[outer]; + while (offset < m_idx_data.outerIndex[outer + 1] && m_indices[offset] != inner) offset++; if (m_indices[offset] == inner) { return Map(&(m_values[blockPtr(offset)]), rsize, csize); } else { @@ -764,8 +764,8 @@ class BlockSparseMatrix StorageIndex csize = IsColMajor ? blockOuterSize(bcol) : blockInnerSize(brow); StorageIndex inner = IsColMajor ? brow : bcol; StorageIndex outer = IsColMajor ? bcol : brow; - StorageIndex offset = m_outerIndex[outer]; - while (offset < m_outerIndex[outer + 1] && m_indices[offset] != inner) offset++; + StorageIndex offset = m_idx_data.outerIndex[outer]; + while (offset < m_idx_data.outerIndex[outer + 1] && m_indices[offset] != inner) offset++; if (m_indices[offset] == inner) { return Map(&(m_values[blockPtr(offset)]), rsize, csize); } else @@ -788,8 +788,8 @@ class BlockSparseMatrix // inline Scalar *valuePtr(){ return m_values; } inline StorageIndex* innerIndexPtr() { return m_indices; } inline const StorageIndex* innerIndexPtr() const { return m_indices; } - inline StorageIndex* outerIndexPtr() { return m_outerIndex; } - inline const StorageIndex* outerIndexPtr() const { return m_outerIndex; } + inline StorageIndex* outerIndexPtr() { return m_idx_data.outerIndex; } + inline const StorageIndex* outerIndexPtr() const { return m_idx_data.outerIndex; } /** \brief for compatibility purposes with the SparseMatrix class */ inline bool isCompressed() const { return true; } @@ -875,8 +875,8 @@ class BlockSparseMatrix StorageIndex* m_blockPtr; // Pointer to the beginning of each block in m_values, size m_nonzeroblocks ... null for // fixed-size blocks StorageIndex* m_indices; // Inner block indices, size m_nonzerosblocks ... OK - StorageIndex* m_outerIndex; // Starting pointer of each block column in m_indices (size m_outerBSize)... OK - Index m_blockSize; // Size of a block for fixed-size blocks, otherwise -1 + StorageIndex* m_idx_data.outerIndex; // Starting pointer of each block column in m_indices (size m_outerBSize)... OK + Index m_blockSize; // Size of a block for fixed-size blocks, otherwise -1 }; template @@ -885,7 +885,10 @@ class BlockSparseMatrix:: enum { Flags = Options_ }; BlockInnerIterator(const BlockSparseMatrix& mat, const Index outer) - : m_mat(mat), m_outer(outer), m_id(mat.m_outerIndex[outer]), m_end(mat.m_outerIndex[outer + 1]) {} + : m_mat(mat), + m_outer(outer), + m_id(mat.m_idx_data.outerIndex[outer]), + m_end(mat.m_idx_data.outerIndex[outer + 1]) {} inline BlockInnerIterator& operator++() { m_id++;