From 44861dcbfeee041223c4aac1ee075e92fa4daa01 Mon Sep 17 00:00:00 2001 From: Stanislaw Halik Date: Sun, 18 Sep 2016 12:42:15 +0200 Subject: update --- eigen/Eigen/src/SparseCore/AmbiVector.h | 373 ++++++ eigen/Eigen/src/SparseCore/CMakeLists.txt | 6 + eigen/Eigen/src/SparseCore/CompressedStorage.h | 240 ++++ .../SparseCore/ConservativeSparseSparseProduct.h | 245 ++++ eigen/Eigen/src/SparseCore/MappedSparseMatrix.h | 181 +++ eigen/Eigen/src/SparseCore/SparseBlock.h | 538 +++++++++ eigen/Eigen/src/SparseCore/SparseColEtree.h | 206 ++++ eigen/Eigen/src/SparseCore/SparseCwiseBinaryOp.h | 324 +++++ eigen/Eigen/src/SparseCore/SparseCwiseUnaryOp.h | 163 +++ eigen/Eigen/src/SparseCore/SparseDenseProduct.h | 311 +++++ eigen/Eigen/src/SparseCore/SparseDiagonalProduct.h | 196 +++ eigen/Eigen/src/SparseCore/SparseDot.h | 101 ++ eigen/Eigen/src/SparseCore/SparseFuzzy.h | 26 + eigen/Eigen/src/SparseCore/SparseMatrix.h | 1262 ++++++++++++++++++++ eigen/Eigen/src/SparseCore/SparseMatrixBase.h | 462 +++++++ eigen/Eigen/src/SparseCore/SparsePermutation.h | 148 +++ eigen/Eigen/src/SparseCore/SparseProduct.h | 188 +++ eigen/Eigen/src/SparseCore/SparseRedux.h | 48 + eigen/Eigen/src/SparseCore/SparseSelfAdjointView.h | 507 ++++++++ .../SparseCore/SparseSparseProductWithPruning.h | 150 +++ eigen/Eigen/src/SparseCore/SparseTranspose.h | 63 + eigen/Eigen/src/SparseCore/SparseTriangularView.h | 179 +++ eigen/Eigen/src/SparseCore/SparseUtil.h | 172 +++ eigen/Eigen/src/SparseCore/SparseVector.h | 448 +++++++ eigen/Eigen/src/SparseCore/SparseView.h | 99 ++ eigen/Eigen/src/SparseCore/TriangularSolver.h | 334 ++++++ 26 files changed, 6970 insertions(+) create mode 100644 eigen/Eigen/src/SparseCore/AmbiVector.h create mode 100644 eigen/Eigen/src/SparseCore/CMakeLists.txt create mode 100644 eigen/Eigen/src/SparseCore/CompressedStorage.h create mode 100644 eigen/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h create mode 100644 eigen/Eigen/src/SparseCore/MappedSparseMatrix.h create mode 100644 eigen/Eigen/src/SparseCore/SparseBlock.h create mode 100644 eigen/Eigen/src/SparseCore/SparseColEtree.h create mode 100644 eigen/Eigen/src/SparseCore/SparseCwiseBinaryOp.h create mode 100644 eigen/Eigen/src/SparseCore/SparseCwiseUnaryOp.h create mode 100644 eigen/Eigen/src/SparseCore/SparseDenseProduct.h create mode 100644 eigen/Eigen/src/SparseCore/SparseDiagonalProduct.h create mode 100644 eigen/Eigen/src/SparseCore/SparseDot.h create mode 100644 eigen/Eigen/src/SparseCore/SparseFuzzy.h create mode 100644 eigen/Eigen/src/SparseCore/SparseMatrix.h create mode 100644 eigen/Eigen/src/SparseCore/SparseMatrixBase.h create mode 100644 eigen/Eigen/src/SparseCore/SparsePermutation.h create mode 100644 eigen/Eigen/src/SparseCore/SparseProduct.h create mode 100644 eigen/Eigen/src/SparseCore/SparseRedux.h create mode 100644 eigen/Eigen/src/SparseCore/SparseSelfAdjointView.h create mode 100644 eigen/Eigen/src/SparseCore/SparseSparseProductWithPruning.h create mode 100644 eigen/Eigen/src/SparseCore/SparseTranspose.h create mode 100644 eigen/Eigen/src/SparseCore/SparseTriangularView.h create mode 100644 eigen/Eigen/src/SparseCore/SparseUtil.h create mode 100644 eigen/Eigen/src/SparseCore/SparseVector.h create mode 100644 eigen/Eigen/src/SparseCore/SparseView.h create mode 100644 eigen/Eigen/src/SparseCore/TriangularSolver.h (limited to 'eigen/Eigen/src/SparseCore') diff --git a/eigen/Eigen/src/SparseCore/AmbiVector.h b/eigen/Eigen/src/SparseCore/AmbiVector.h new file mode 100644 index 0000000..220c645 --- /dev/null +++ b/eigen/Eigen/src/SparseCore/AmbiVector.h @@ -0,0 +1,373 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// +// 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_AMBIVECTOR_H +#define EIGEN_AMBIVECTOR_H + +namespace Eigen { + +namespace internal { + +/** \internal + * Hybrid sparse/dense vector class designed for intensive read-write operations. + * + * See BasicSparseLLT and SparseProduct for usage examples. + */ +template +class AmbiVector +{ + public: + typedef _Scalar Scalar; + typedef _Index Index; + typedef typename NumTraits::Real RealScalar; + + AmbiVector(Index size) + : m_buffer(0), m_zero(0), m_size(0), m_allocatedSize(0), m_allocatedElements(0), m_mode(-1) + { + resize(size); + } + + void init(double estimatedDensity); + void init(int mode); + + Index nonZeros() const; + + /** Specifies a sub-vector to work on */ + void setBounds(Index start, Index end) { m_start = start; m_end = end; } + + void setZero(); + + void restart(); + Scalar& coeffRef(Index i); + Scalar& coeff(Index i); + + class Iterator; + + ~AmbiVector() { delete[] m_buffer; } + + void resize(Index size) + { + if (m_allocatedSize < size) + reallocate(size); + m_size = size; + } + + Index size() const { return m_size; } + + protected: + + void reallocate(Index size) + { + // if the size of the matrix is not too large, let's allocate a bit more than needed such + // that we can handle dense vector even in sparse mode. + delete[] m_buffer; + if (size<1000) + { + Index allocSize = (size * sizeof(ListEl) + sizeof(Scalar) - 1)/sizeof(Scalar); + m_allocatedElements = (allocSize*sizeof(Scalar))/sizeof(ListEl); + m_buffer = new Scalar[allocSize]; + } + else + { + m_allocatedElements = (size*sizeof(Scalar))/sizeof(ListEl); + m_buffer = new Scalar[size]; + } + m_size = size; + m_start = 0; + m_end = m_size; + } + + void reallocateSparse() + { + Index copyElements = m_allocatedElements; + m_allocatedElements = (std::min)(Index(m_allocatedElements*1.5),m_size); + Index allocSize = m_allocatedElements * sizeof(ListEl); + allocSize = (allocSize + sizeof(Scalar) - 1)/sizeof(Scalar); + Scalar* newBuffer = new Scalar[allocSize]; + memcpy(newBuffer, m_buffer, copyElements * sizeof(ListEl)); + delete[] m_buffer; + m_buffer = newBuffer; + } + + protected: + // element type of the linked list + struct ListEl + { + Index next; + Index index; + Scalar value; + }; + + // used to store data in both mode + Scalar* m_buffer; + Scalar m_zero; + Index m_size; + Index m_start; + Index m_end; + Index m_allocatedSize; + Index m_allocatedElements; + Index m_mode; + + // linked list mode + Index m_llStart; + Index m_llCurrent; + Index m_llSize; +}; + +/** \returns the number of non zeros in the current sub vector */ +template +_Index AmbiVector<_Scalar,_Index>::nonZeros() const +{ + if (m_mode==IsSparse) + return m_llSize; + else + return m_end - m_start; +} + +template +void AmbiVector<_Scalar,_Index>::init(double estimatedDensity) +{ + if (estimatedDensity>0.1) + init(IsDense); + else + init(IsSparse); +} + +template +void AmbiVector<_Scalar,_Index>::init(int mode) +{ + m_mode = mode; + if (m_mode==IsSparse) + { + m_llSize = 0; + m_llStart = -1; + } +} + +/** Must be called whenever we might perform a write access + * with an index smaller than the previous one. + * + * Don't worry, this function is extremely cheap. + */ +template +void AmbiVector<_Scalar,_Index>::restart() +{ + m_llCurrent = m_llStart; +} + +/** Set all coefficients of current subvector to zero */ +template +void AmbiVector<_Scalar,_Index>::setZero() +{ + if (m_mode==IsDense) + { + for (Index i=m_start; i +_Scalar& AmbiVector<_Scalar,_Index>::coeffRef(_Index i) +{ + if (m_mode==IsDense) + return m_buffer[i]; + else + { + ListEl* EIGEN_RESTRICT llElements = reinterpret_cast(m_buffer); + // TODO factorize the following code to reduce code generation + eigen_assert(m_mode==IsSparse); + if (m_llSize==0) + { + // this is the first element + m_llStart = 0; + m_llCurrent = 0; + ++m_llSize; + llElements[0].value = Scalar(0); + llElements[0].index = i; + llElements[0].next = -1; + return llElements[0].value; + } + else if (i=llElements[m_llCurrent].index && "you must call restart() before inserting an element with lower or equal index"); + while (nextel >= 0 && llElements[nextel].index<=i) + { + m_llCurrent = nextel; + nextel = llElements[nextel].next; + } + + if (llElements[m_llCurrent].index==i) + { + // the coefficient already exists and we found it ! + return llElements[m_llCurrent].value; + } + else + { + if (m_llSize>=m_allocatedElements) + { + reallocateSparse(); + llElements = reinterpret_cast(m_buffer); + } + eigen_internal_assert(m_llSize +_Scalar& AmbiVector<_Scalar,_Index>::coeff(_Index i) +{ + if (m_mode==IsDense) + return m_buffer[i]; + else + { + ListEl* EIGEN_RESTRICT llElements = reinterpret_cast(m_buffer); + eigen_assert(m_mode==IsSparse); + if ((m_llSize==0) || (i= 0 && llElements[elid].index +class AmbiVector<_Scalar,_Index>::Iterator +{ + public: + typedef _Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + + /** Default constructor + * \param vec the vector on which we iterate + * \param epsilon the minimal value used to prune zero coefficients. + * In practice, all coefficients having a magnitude smaller than \a epsilon + * are skipped. + */ + Iterator(const AmbiVector& vec, const RealScalar& epsilon = 0) + : m_vector(vec) + { + using std::abs; + m_epsilon = epsilon; + m_isDense = m_vector.m_mode==IsDense; + if (m_isDense) + { + m_currentEl = 0; // this is to avoid a compilation warning + m_cachedValue = 0; // this is to avoid a compilation warning + m_cachedIndex = m_vector.m_start-1; + ++(*this); + } + else + { + ListEl* EIGEN_RESTRICT llElements = reinterpret_cast(m_vector.m_buffer); + m_currentEl = m_vector.m_llStart; + while (m_currentEl>=0 && abs(llElements[m_currentEl].value)<=m_epsilon) + m_currentEl = llElements[m_currentEl].next; + if (m_currentEl<0) + { + m_cachedValue = 0; // this is to avoid a compilation warning + m_cachedIndex = -1; + } + else + { + m_cachedIndex = llElements[m_currentEl].index; + m_cachedValue = llElements[m_currentEl].value; + } + } + } + + Index index() const { return m_cachedIndex; } + Scalar value() const { return m_cachedValue; } + + operator bool() const { return m_cachedIndex>=0; } + + Iterator& operator++() + { + using std::abs; + if (m_isDense) + { + do { + ++m_cachedIndex; + } while (m_cachedIndex(m_vector.m_buffer); + do { + m_currentEl = llElements[m_currentEl].next; + } while (m_currentEl>=0 && abs(llElements[m_currentEl].value) +// +// 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_COMPRESSED_STORAGE_H +#define EIGEN_COMPRESSED_STORAGE_H + +namespace Eigen { + +namespace internal { + +/** \internal + * Stores a sparse set of values as a list of values and a list of indices. + * + */ +template +class CompressedStorage +{ + public: + + typedef _Scalar Scalar; + typedef _Index Index; + + protected: + + typedef typename NumTraits::Real RealScalar; + + public: + + CompressedStorage() + : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0) + {} + + CompressedStorage(size_t size) + : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0) + { + resize(size); + } + + CompressedStorage(const CompressedStorage& other) + : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0) + { + *this = other; + } + + CompressedStorage& operator=(const CompressedStorage& other) + { + resize(other.size()); + 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); + 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); + } + + ~CompressedStorage() + { + delete[] m_values; + delete[] m_indices; + } + + void reserve(size_t size) + { + size_t newAllocatedSize = m_size + size; + if (newAllocatedSize > m_allocatedSize) + reallocate(newAllocatedSize); + } + + void squeeze() + { + if (m_allocatedSize>m_size) + reallocate(m_size); + } + + void resize(size_t size, double reserveSizeFactor = 0) + { + if (m_allocatedSize(m_size); + resize(m_size+1, 1); + m_values[id] = v; + m_indices[id] = i; + } + + inline size_t size() const { return m_size; } + inline size_t allocatedSize() const { return m_allocatedSize; } + inline void clear() { m_size = 0; } + + const Scalar* valuePtr() const { return m_values; } + Scalar* valuePtr() { return m_values; } + const Index* indexPtr() const { return m_indices; } + Index* indexPtr() { return m_indices; } + + inline Scalar& value(size_t i) { return m_values[i]; } + inline const Scalar& value(size_t i) const { return m_values[i]; } + + inline Index& index(size_t i) { return m_indices[i]; } + inline const Index& index(size_t i) const { return m_indices[i]; } + + static CompressedStorage Map(Index* indices, Scalar* values, size_t size) + { + CompressedStorage res; + res.m_indices = indices; + res.m_values = values; + res.m_allocatedSize = res.m_size = size; + return res; + } + + /** \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); + } + + /** \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(size_t start, size_t end, Index key) const + { + while(end>start) + { + size_t mid = (end+start)>>1; + if (m_indices[mid](start); + } + + /** \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) + return defaultValue; + else if (key==m_indices[m_size-1]) + return m_values[m_size-1]; + // ^^ optimization: let's first check if it is the last coefficient + // (very common in high level algorithms) + const size_t id = searchLowerIndex(0,m_size-1,key); + return ((id=end) + return Scalar(0); + else if (end>start && key==m_indices[end-1]) + return m_values[end-1]; + // ^^ optimization: let's first check if it is the last coefficient + // (very common in high level algorithms) + const size_t id = searchLowerIndex(start,end-1,key); + return ((id=m_size || m_indices[id]!=key) + { + resize(m_size+1,1); + for (size_t j=m_size-1; j>id; --j) + { + m_indices[j] = m_indices[j-1]; + m_values[j] = m_values[j-1]; + } + m_indices[id] = key; + m_values[id] = defaultValue; + } + return m_values[id]; + } + + void prune(const Scalar& reference, const RealScalar& epsilon = NumTraits::dummy_precision()) + { + size_t k = 0; + size_t n = size(); + for (size_t i=0; i0) { + internal::smart_copy(m_values, m_values+copySize, newValues); + internal::smart_copy(m_indices, m_indices+copySize, newIndices); + } + // delete old stuff + delete[] m_values; + delete[] m_indices; + m_values = newValues; + m_indices = newIndices; + m_allocatedSize = size; + } + + protected: + Scalar* m_values; + Index* m_indices; + size_t m_size; + size_t m_allocatedSize; + +}; + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_COMPRESSED_STORAGE_H diff --git a/eigen/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h b/eigen/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h new file mode 100644 index 0000000..5c320e2 --- /dev/null +++ b/eigen/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h @@ -0,0 +1,245 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2011 Gael Guennebaud +// +// 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_CONSERVATIVESPARSESPARSEPRODUCT_H +#define EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H + +namespace Eigen { + +namespace internal { + +template +static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res) +{ + typedef typename remove_all::type::Scalar Scalar; + typedef typename remove_all::type::Index Index; + + // make sure to call innerSize/outerSize since we fake the storage order. + Index rows = lhs.innerSize(); + Index cols = rhs.outerSize(); + eigen_assert(lhs.outerSize() == rhs.innerSize()); + + std::vector mask(rows,false); + Matrix values(rows); + Matrix indices(rows); + + // estimate the number of non zero entries + // given a rhs column containing Y non zeros, we assume that the respective Y columns + // of the lhs differs in average of one non zeros, thus the number of non zeros for + // the product of a rhs column with the lhs is X+Y where X is the average number of non zero + // per column of the lhs. + // Therefore, we have nnz(lhs*rhs) = nnz(lhs) + nnz(rhs) + Index estimated_nnz_prod = lhs.nonZeros() + rhs.nonZeros(); + + res.setZero(); + res.reserve(Index(estimated_nnz_prod)); + // we compute each column of the result, one after the other + for (Index j=0; j use a quick sort + // otherwise => loop through the entire vector + // In order to avoid to perform an expensive log2 when the + // result is clearly very sparse we use a linear bound up to 200. + //if((nnz<200 && nnz1) std::sort(indices.data(),indices.data()+nnz); + for(Index k=0; k::Flags&RowMajorBit) ? RowMajor : ColMajor, + int RhsStorageOrder = (traits::Flags&RowMajorBit) ? RowMajor : ColMajor, + int ResStorageOrder = (traits::Flags&RowMajorBit) ? RowMajor : ColMajor> +struct conservative_sparse_sparse_product_selector; + +template +struct conservative_sparse_sparse_product_selector +{ + typedef typename remove_all::type LhsCleaned; + typedef typename LhsCleaned::Scalar Scalar; + + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix RowMajorMatrix; + typedef SparseMatrix ColMajorMatrix; + ColMajorMatrix resCol(lhs.rows(),rhs.cols()); + internal::conservative_sparse_sparse_product_impl(lhs, rhs, resCol); + // sort the non zeros: + RowMajorMatrix resRow(resCol); + res = resRow; + } +}; + +template +struct conservative_sparse_sparse_product_selector +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix RowMajorMatrix; + RowMajorMatrix rhsRow = rhs; + RowMajorMatrix resRow(lhs.rows(), rhs.cols()); + internal::conservative_sparse_sparse_product_impl(rhsRow, lhs, resRow); + res = resRow; + } +}; + +template +struct conservative_sparse_sparse_product_selector +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix RowMajorMatrix; + RowMajorMatrix lhsRow = lhs; + RowMajorMatrix resRow(lhs.rows(), rhs.cols()); + internal::conservative_sparse_sparse_product_impl(rhs, lhsRow, resRow); + res = resRow; + } +}; + +template +struct conservative_sparse_sparse_product_selector +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix RowMajorMatrix; + RowMajorMatrix resRow(lhs.rows(), rhs.cols()); + internal::conservative_sparse_sparse_product_impl(rhs, lhs, resRow); + res = resRow; + } +}; + + +template +struct conservative_sparse_sparse_product_selector +{ + typedef typename traits::type>::Scalar Scalar; + + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix ColMajorMatrix; + ColMajorMatrix resCol(lhs.rows(), rhs.cols()); + internal::conservative_sparse_sparse_product_impl(lhs, rhs, resCol); + res = resCol; + } +}; + +template +struct conservative_sparse_sparse_product_selector +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix ColMajorMatrix; + ColMajorMatrix lhsCol = lhs; + ColMajorMatrix resCol(lhs.rows(), rhs.cols()); + internal::conservative_sparse_sparse_product_impl(lhsCol, rhs, resCol); + res = resCol; + } +}; + +template +struct conservative_sparse_sparse_product_selector +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix ColMajorMatrix; + ColMajorMatrix rhsCol = rhs; + ColMajorMatrix resCol(lhs.rows(), rhs.cols()); + internal::conservative_sparse_sparse_product_impl(lhs, rhsCol, resCol); + res = resCol; + } +}; + +template +struct conservative_sparse_sparse_product_selector +{ + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) + { + typedef SparseMatrix RowMajorMatrix; + typedef SparseMatrix ColMajorMatrix; + RowMajorMatrix resRow(lhs.rows(),rhs.cols()); + internal::conservative_sparse_sparse_product_impl(rhs, lhs, resRow); + // sort the non zeros: + ColMajorMatrix resCol(resRow); + res = resCol; + } +}; + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H diff --git a/eigen/Eigen/src/SparseCore/MappedSparseMatrix.h b/eigen/Eigen/src/SparseCore/MappedSparseMatrix.h new file mode 100644 index 0000000..ab1a266 --- /dev/null +++ b/eigen/Eigen/src/SparseCore/MappedSparseMatrix.h @@ -0,0 +1,181 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// +// 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_MAPPED_SPARSEMATRIX_H +#define EIGEN_MAPPED_SPARSEMATRIX_H + +namespace Eigen { + +/** \class MappedSparseMatrix + * + * \brief Sparse matrix + * + * \param _Scalar the scalar type, i.e. the type of the coefficients + * + * See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme. + * + */ +namespace internal { +template +struct traits > : traits > +{}; +} + +template +class MappedSparseMatrix + : public SparseMatrixBase > +{ + public: + EIGEN_SPARSE_PUBLIC_INTERFACE(MappedSparseMatrix) + enum { IsRowMajor = Base::IsRowMajor }; + + protected: + + Index m_outerSize; + Index m_innerSize; + Index m_nnz; + Index* m_outerIndex; + Index* m_innerIndices; + Scalar* m_values; + + public: + + inline Index rows() const { return IsRowMajor ? m_outerSize : m_innerSize; } + inline Index cols() const { return IsRowMajor ? m_innerSize : m_outerSize; } + inline Index innerSize() const { return m_innerSize; } + inline Index outerSize() const { return m_outerSize; } + + bool isCompressed() const { return true; } + + //---------------------------------------- + // direct access interface + inline const Scalar* valuePtr() const { return m_values; } + inline Scalar* valuePtr() { return m_values; } + + inline const Index* innerIndexPtr() const { return m_innerIndices; } + inline Index* innerIndexPtr() { return m_innerIndices; } + + inline const Index* outerIndexPtr() const { return m_outerIndex; } + inline Index* outerIndexPtr() { return m_outerIndex; } + //---------------------------------------- + + inline Scalar coeff(Index row, Index col) const + { + const Index outer = IsRowMajor ? row : col; + const Index inner = IsRowMajor ? col : row; + + Index start = m_outerIndex[outer]; + Index end = m_outerIndex[outer+1]; + if (start==end) + return Scalar(0); + else if (end>0 && inner==m_innerIndices[end-1]) + return m_values[end-1]; + // ^^ optimization: let's first check if it is the last coefficient + // (very common in high level algorithms) + + const Index* r = std::lower_bound(&m_innerIndices[start],&m_innerIndices[end-1],inner); + const Index id = r-&m_innerIndices[0]; + return ((*r==inner) && (id=start && "you probably called coeffRef on a non finalized matrix"); + eigen_assert(end>start && "coeffRef cannot be called on a zero coefficient"); + Index* r = std::lower_bound(&m_innerIndices[start],&m_innerIndices[end],inner); + const Index id = r-&m_innerIndices[0]; + eigen_assert((*r==inner) && (id +class MappedSparseMatrix::InnerIterator +{ + public: + InnerIterator(const MappedSparseMatrix& mat, Index outer) + : m_matrix(mat), + m_outer(outer), + m_id(mat.outerIndexPtr()[outer]), + m_start(m_id), + m_end(mat.outerIndexPtr()[outer+1]) + {} + + inline InnerIterator& operator++() { m_id++; return *this; } + + inline Scalar value() const { return m_matrix.valuePtr()[m_id]; } + inline Scalar& valueRef() { return const_cast(m_matrix.valuePtr()[m_id]); } + + inline Index index() const { return m_matrix.innerIndexPtr()[m_id]; } + inline Index row() const { return IsRowMajor ? m_outer : index(); } + inline Index col() const { return IsRowMajor ? index() : m_outer; } + + inline operator bool() const { return (m_id < m_end) && (m_id>=m_start); } + + protected: + const MappedSparseMatrix& m_matrix; + const Index m_outer; + Index m_id; + const Index m_start; + const Index m_end; +}; + +template +class MappedSparseMatrix::ReverseInnerIterator +{ + public: + ReverseInnerIterator(const MappedSparseMatrix& mat, Index outer) + : m_matrix(mat), + m_outer(outer), + m_id(mat.outerIndexPtr()[outer+1]), + m_start(mat.outerIndexPtr()[outer]), + m_end(m_id) + {} + + inline ReverseInnerIterator& operator--() { m_id--; return *this; } + + inline Scalar value() const { return m_matrix.valuePtr()[m_id-1]; } + inline Scalar& valueRef() { return const_cast(m_matrix.valuePtr()[m_id-1]); } + + inline Index index() const { return m_matrix.innerIndexPtr()[m_id-1]; } + inline Index row() const { return IsRowMajor ? m_outer : index(); } + inline Index col() const { return IsRowMajor ? index() : m_outer; } + + inline operator bool() const { return (m_id <= m_end) && (m_id>m_start); } + + protected: + const MappedSparseMatrix& m_matrix; + const Index m_outer; + Index m_id; + const Index m_start; + const Index m_end; +}; + +} // end namespace Eigen + +#endif // EIGEN_MAPPED_SPARSEMATRIX_H diff --git a/eigen/Eigen/src/SparseCore/SparseBlock.h b/eigen/Eigen/src/SparseCore/SparseBlock.h new file mode 100644 index 0000000..ac4124b --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseBlock.h @@ -0,0 +1,538 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 Gael Guennebaud +// +// 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_SPARSE_BLOCK_H +#define EIGEN_SPARSE_BLOCK_H + +namespace Eigen { + +template +class BlockImpl + : public SparseMatrixBase > +{ + typedef typename internal::remove_all::type _MatrixTypeNested; + typedef Block BlockType; +public: + enum { IsRowMajor = internal::traits::IsRowMajor }; +protected: + enum { OuterSize = IsRowMajor ? BlockRows : BlockCols }; +public: + EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) + + class InnerIterator: public XprType::InnerIterator + { + typedef typename BlockImpl::Index Index; + public: + inline InnerIterator(const BlockType& xpr, Index outer) + : XprType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) + {} + inline Index row() const { return IsRowMajor ? m_outer : this->index(); } + inline Index col() const { return IsRowMajor ? this->index() : m_outer; } + protected: + Index m_outer; + }; + class ReverseInnerIterator: public XprType::ReverseInnerIterator + { + typedef typename BlockImpl::Index Index; + public: + inline ReverseInnerIterator(const BlockType& xpr, Index outer) + : XprType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) + {} + inline Index row() const { return IsRowMajor ? m_outer : this->index(); } + inline Index col() const { return IsRowMajor ? this->index() : m_outer; } + protected: + Index m_outer; + }; + + inline BlockImpl(const XprType& xpr, int i) + : m_matrix(xpr), m_outerStart(i), m_outerSize(OuterSize) + {} + + inline BlockImpl(const XprType& xpr, int startRow, int startCol, int blockRows, int blockCols) + : m_matrix(xpr), m_outerStart(IsRowMajor ? startRow : startCol), m_outerSize(IsRowMajor ? blockRows : blockCols) + {} + + inline const Scalar coeff(int row, int col) const + { + return m_matrix.coeff(row + IsRowMajor ? m_outerStart : 0, col +IsRowMajor ? 0 : m_outerStart); + } + + inline const Scalar coeff(int index) const + { + return m_matrix.coeff(IsRowMajor ? m_outerStart : index, IsRowMajor ? index : m_outerStart); + } + + EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); } + EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); } + + protected: + + typename XprType::Nested m_matrix; + Index m_outerStart; + const internal::variable_if_dynamic m_outerSize; + + EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl) + private: + Index nonZeros() const; +}; + + +/*************************************************************************** +* specialisation for SparseMatrix +***************************************************************************/ + +template +class BlockImpl,BlockRows,BlockCols,true,Sparse> + : public SparseMatrixBase,BlockRows,BlockCols,true> > +{ + typedef SparseMatrix<_Scalar, _Options, _Index> SparseMatrixType; + typedef typename internal::remove_all::type _MatrixTypeNested; + typedef Block BlockType; + typedef Block ConstBlockType; +public: + enum { IsRowMajor = internal::traits::IsRowMajor }; + EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) +protected: + enum { OuterSize = IsRowMajor ? BlockRows : BlockCols }; +public: + + class InnerIterator: public SparseMatrixType::InnerIterator + { + public: + inline InnerIterator(const BlockType& xpr, Index outer) + : SparseMatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) + {} + inline Index row() const { return IsRowMajor ? m_outer : this->index(); } + inline Index col() const { return IsRowMajor ? this->index() : m_outer; } + protected: + Index m_outer; + }; + class ReverseInnerIterator: public SparseMatrixType::ReverseInnerIterator + { + public: + inline ReverseInnerIterator(const BlockType& xpr, Index outer) + : SparseMatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) + {} + inline Index row() const { return IsRowMajor ? m_outer : this->index(); } + inline Index col() const { return IsRowMajor ? this->index() : m_outer; } + protected: + Index m_outer; + }; + + inline BlockImpl(const SparseMatrixType& xpr, int i) + : m_matrix(xpr), m_outerStart(i), m_outerSize(OuterSize) + {} + + inline BlockImpl(const SparseMatrixType& xpr, int startRow, int startCol, int blockRows, int blockCols) + : m_matrix(xpr), m_outerStart(IsRowMajor ? startRow : startCol), m_outerSize(IsRowMajor ? blockRows : blockCols) + {} + + template + inline BlockType& operator=(const SparseMatrixBase& other) + { + typedef typename internal::remove_all::type _NestedMatrixType; + _NestedMatrixType& matrix = const_cast<_NestedMatrixType&>(m_matrix);; + // This assignement is slow if this vector set is not empty + // and/or it is not at the end of the nonzeros of the underlying matrix. + + // 1 - eval to a temporary to avoid transposition and/or aliasing issues + SparseMatrix tmp(other); + + // 2 - let's check whether there is enough allocated memory + Index nnz = tmp.nonZeros(); + Index start = m_outerStart==0 ? 0 : matrix.outerIndexPtr()[m_outerStart]; // starting position of the current block + Index end = m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()]; // ending posiiton of the current block + Index block_size = end - start; // available room in the current block + Index tail_size = m_matrix.outerIndexPtr()[m_matrix.outerSize()] - end; + + Index free_size = m_matrix.isCompressed() + ? Index(matrix.data().allocatedSize()) + block_size + : block_size; + + if(nnz>free_size) + { + // realloc manually to reduce copies + typename SparseMatrixType::Storage newdata(m_matrix.data().allocatedSize() - block_size + nnz); + + std::memcpy(newdata.valuePtr(), m_matrix.data().valuePtr(), start*sizeof(Scalar)); + std::memcpy(newdata.indexPtr(), m_matrix.data().indexPtr(), start*sizeof(Index)); + + std::memcpy(newdata.valuePtr() + start, tmp.data().valuePtr(), nnz*sizeof(Scalar)); + std::memcpy(newdata.indexPtr() + start, tmp.data().indexPtr(), nnz*sizeof(Index)); + + std::memcpy(newdata.valuePtr()+start+nnz, matrix.data().valuePtr()+end, tail_size*sizeof(Scalar)); + std::memcpy(newdata.indexPtr()+start+nnz, matrix.data().indexPtr()+end, tail_size*sizeof(Index)); + + newdata.resize(m_matrix.outerIndexPtr()[m_matrix.outerSize()] - block_size + nnz); + + matrix.data().swap(newdata); + } + else + { + // no need to realloc, simply copy the tail at its respective position and insert tmp + matrix.data().resize(start + nnz + tail_size); + + std::memmove(matrix.data().valuePtr()+start+nnz, matrix.data().valuePtr()+end, tail_size*sizeof(Scalar)); + std::memmove(matrix.data().indexPtr()+start+nnz, matrix.data().indexPtr()+end, tail_size*sizeof(Index)); + + std::memcpy(matrix.data().valuePtr()+start, tmp.data().valuePtr(), nnz*sizeof(Scalar)); + std::memcpy(matrix.data().indexPtr()+start, tmp.data().indexPtr(), nnz*sizeof(Index)); + } + + // update innerNonZeros + if(!m_matrix.isCompressed()) + for(Index j=0; j(other); + } + + inline const Scalar* valuePtr() const + { return m_matrix.valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; } + inline Scalar* valuePtr() + { return m_matrix.const_cast_derived().valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; } + + inline const Index* innerIndexPtr() const + { return m_matrix.innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; } + inline Index* innerIndexPtr() + { return m_matrix.const_cast_derived().innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; } + + inline const Index* outerIndexPtr() const + { return m_matrix.outerIndexPtr() + m_outerStart; } + inline Index* outerIndexPtr() + { return m_matrix.const_cast_derived().outerIndexPtr() + m_outerStart; } + + Index nonZeros() const + { + if(m_matrix.isCompressed()) + return std::size_t(m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()]) + - std::size_t(m_matrix.outerIndexPtr()[m_outerStart]); + else if(m_outerSize.value()==0) + return 0; + else + return Map >(m_matrix.innerNonZeroPtr()+m_outerStart, m_outerSize.value()).sum(); + } + + inline Scalar& coeffRef(int row, int col) + { + return m_matrix.const_cast_derived().coeffRef(row + (IsRowMajor ? m_outerStart : 0), col + (IsRowMajor ? 0 : m_outerStart)); + } + + inline const Scalar coeff(int row, int col) const + { + return m_matrix.coeff(row + (IsRowMajor ? m_outerStart : 0), col + (IsRowMajor ? 0 : m_outerStart)); + } + + inline const Scalar coeff(int index) const + { + return m_matrix.coeff(IsRowMajor ? m_outerStart : index, IsRowMajor ? index : m_outerStart); + } + + const Scalar& lastCoeff() const + { + EIGEN_STATIC_ASSERT_VECTOR_ONLY(BlockImpl); + eigen_assert(nonZeros()>0); + if(m_matrix.isCompressed()) + return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart+1]-1]; + else + return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart]+m_matrix.innerNonZeroPtr()[m_outerStart]-1]; + } + + EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); } + EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); } + + protected: + + typename SparseMatrixType::Nested m_matrix; + Index m_outerStart; + const internal::variable_if_dynamic m_outerSize; + +}; + + +template +class BlockImpl,BlockRows,BlockCols,true,Sparse> + : public SparseMatrixBase,BlockRows,BlockCols,true> > +{ + typedef SparseMatrix<_Scalar, _Options, _Index> SparseMatrixType; + typedef typename internal::remove_all::type _MatrixTypeNested; + typedef Block BlockType; +public: + enum { IsRowMajor = internal::traits::IsRowMajor }; + EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) +protected: + enum { OuterSize = IsRowMajor ? BlockRows : BlockCols }; +public: + + class InnerIterator: public SparseMatrixType::InnerIterator + { + public: + inline InnerIterator(const BlockType& xpr, Index outer) + : SparseMatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) + {} + inline Index row() const { return IsRowMajor ? m_outer : this->index(); } + inline Index col() const { return IsRowMajor ? this->index() : m_outer; } + protected: + Index m_outer; + }; + class ReverseInnerIterator: public SparseMatrixType::ReverseInnerIterator + { + public: + inline ReverseInnerIterator(const BlockType& xpr, Index outer) + : SparseMatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) + {} + inline Index row() const { return IsRowMajor ? m_outer : this->index(); } + inline Index col() const { return IsRowMajor ? this->index() : m_outer; } + protected: + Index m_outer; + }; + + inline BlockImpl(const SparseMatrixType& xpr, int i) + : m_matrix(xpr), m_outerStart(i), m_outerSize(OuterSize) + {} + + inline BlockImpl(const SparseMatrixType& xpr, int startRow, int startCol, int blockRows, int blockCols) + : m_matrix(xpr), m_outerStart(IsRowMajor ? startRow : startCol), m_outerSize(IsRowMajor ? blockRows : blockCols) + {} + + inline const Scalar* valuePtr() const + { return m_matrix.valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; } + + inline const Index* innerIndexPtr() const + { return m_matrix.innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; } + + inline const Index* outerIndexPtr() const + { return m_matrix.outerIndexPtr() + m_outerStart; } + + Index nonZeros() const + { + if(m_matrix.isCompressed()) + return std::size_t(m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()]) + - std::size_t(m_matrix.outerIndexPtr()[m_outerStart]); + else if(m_outerSize.value()==0) + return 0; + else + return Map >(m_matrix.innerNonZeroPtr()+m_outerStart, m_outerSize.value()).sum(); + } + + inline const Scalar coeff(int row, int col) const + { + return m_matrix.coeff(row + (IsRowMajor ? m_outerStart : 0), col + (IsRowMajor ? 0 : m_outerStart)); + } + + inline const Scalar coeff(int index) const + { + return m_matrix.coeff(IsRowMajor ? m_outerStart : index, IsRowMajor ? index : m_outerStart); + } + + const Scalar& lastCoeff() const + { + EIGEN_STATIC_ASSERT_VECTOR_ONLY(BlockImpl); + eigen_assert(nonZeros()>0); + if(m_matrix.isCompressed()) + return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart+1]-1]; + else + return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart]+m_matrix.innerNonZeroPtr()[m_outerStart]-1]; + } + + EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); } + EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); } + + protected: + + EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl) + + typename SparseMatrixType::Nested m_matrix; + Index m_outerStart; + const internal::variable_if_dynamic m_outerSize; +}; + +//---------- + +/** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this + * is col-major (resp. row-major). + */ +template +typename SparseMatrixBase::InnerVectorReturnType SparseMatrixBase::innerVector(Index outer) +{ return InnerVectorReturnType(derived(), outer); } + +/** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this + * is col-major (resp. row-major). Read-only. + */ +template +const typename SparseMatrixBase::ConstInnerVectorReturnType SparseMatrixBase::innerVector(Index outer) const +{ return ConstInnerVectorReturnType(derived(), outer); } + +/** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this + * is col-major (resp. row-major). + */ +template +typename SparseMatrixBase::InnerVectorsReturnType +SparseMatrixBase::innerVectors(Index outerStart, Index outerSize) +{ + return Block(derived(), + IsRowMajor ? outerStart : 0, IsRowMajor ? 0 : outerStart, + IsRowMajor ? outerSize : rows(), IsRowMajor ? cols() : outerSize); + +} + +/** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this + * is col-major (resp. row-major). Read-only. + */ +template +const typename SparseMatrixBase::ConstInnerVectorsReturnType +SparseMatrixBase::innerVectors(Index outerStart, Index outerSize) const +{ + return Block(derived(), + IsRowMajor ? outerStart : 0, IsRowMajor ? 0 : outerStart, + IsRowMajor ? outerSize : rows(), IsRowMajor ? cols() : outerSize); + +} + +/** Generic implementation of sparse Block expression. + * Real-only. + */ +template +class BlockImpl + : public SparseMatrixBase >, internal::no_assignment_operator +{ + typedef typename internal::remove_all::type _MatrixTypeNested; + typedef Block BlockType; +public: + enum { IsRowMajor = internal::traits::IsRowMajor }; + EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) + + /** Column or Row constructor + */ + inline BlockImpl(const XprType& xpr, int i) + : m_matrix(xpr), + m_startRow( (BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) ? i : 0), + m_startCol( (BlockRows==XprType::RowsAtCompileTime) && (BlockCols==1) ? i : 0), + m_blockRows(BlockRows==1 ? 1 : xpr.rows()), + m_blockCols(BlockCols==1 ? 1 : xpr.cols()) + {} + + /** Dynamic-size constructor + */ + inline BlockImpl(const XprType& xpr, int startRow, int startCol, int blockRows, int blockCols) + : m_matrix(xpr), m_startRow(startRow), m_startCol(startCol), m_blockRows(blockRows), m_blockCols(blockCols) + {} + + inline int rows() const { return m_blockRows.value(); } + inline int cols() const { return m_blockCols.value(); } + + inline Scalar& coeffRef(int row, int col) + { + return m_matrix.const_cast_derived() + .coeffRef(row + m_startRow.value(), col + m_startCol.value()); + } + + inline const Scalar coeff(int row, int col) const + { + return m_matrix.coeff(row + m_startRow.value(), col + m_startCol.value()); + } + + inline Scalar& coeffRef(int index) + { + return m_matrix.const_cast_derived() + .coeffRef(m_startRow.value() + (RowsAtCompileTime == 1 ? 0 : index), + m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0)); + } + + inline const Scalar coeff(int index) const + { + return m_matrix + .coeff(m_startRow.value() + (RowsAtCompileTime == 1 ? 0 : index), + m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0)); + } + + inline const _MatrixTypeNested& nestedExpression() const { return m_matrix; } + + class InnerIterator : public _MatrixTypeNested::InnerIterator + { + typedef typename _MatrixTypeNested::InnerIterator Base; + const BlockType& m_block; + Index m_end; + public: + + EIGEN_STRONG_INLINE InnerIterator(const BlockType& block, Index outer) + : Base(block.derived().nestedExpression(), outer + (IsRowMajor ? block.m_startRow.value() : block.m_startCol.value())), + m_block(block), + m_end(IsRowMajor ? block.m_startCol.value()+block.m_blockCols.value() : block.m_startRow.value()+block.m_blockRows.value()) + { + while( (Base::operator bool()) && (Base::index() < (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value())) ) + Base::operator++(); + } + + inline Index index() const { return Base::index() - (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value()); } + inline Index outer() const { return Base::outer() - (IsRowMajor ? m_block.m_startRow.value() : m_block.m_startCol.value()); } + inline Index row() const { return Base::row() - m_block.m_startRow.value(); } + inline Index col() const { return Base::col() - m_block.m_startCol.value(); } + + inline operator bool() const { return Base::operator bool() && Base::index() < m_end; } + }; + class ReverseInnerIterator : public _MatrixTypeNested::ReverseInnerIterator + { + typedef typename _MatrixTypeNested::ReverseInnerIterator Base; + const BlockType& m_block; + Index m_begin; + public: + + EIGEN_STRONG_INLINE ReverseInnerIterator(const BlockType& block, Index outer) + : Base(block.derived().nestedExpression(), outer + (IsRowMajor ? block.m_startRow.value() : block.m_startCol.value())), + m_block(block), + m_begin(IsRowMajor ? block.m_startCol.value() : block.m_startRow.value()) + { + while( (Base::operator bool()) && (Base::index() >= (IsRowMajor ? m_block.m_startCol.value()+block.m_blockCols.value() : m_block.m_startRow.value()+block.m_blockRows.value())) ) + Base::operator--(); + } + + inline Index index() const { return Base::index() - (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value()); } + inline Index outer() const { return Base::outer() - (IsRowMajor ? m_block.m_startRow.value() : m_block.m_startCol.value()); } + inline Index row() const { return Base::row() - m_block.m_startRow.value(); } + inline Index col() const { return Base::col() - m_block.m_startCol.value(); } + + inline operator bool() const { return Base::operator bool() && Base::index() >= m_begin; } + }; + protected: + friend class InnerIterator; + friend class ReverseInnerIterator; + + EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl) + + typename XprType::Nested m_matrix; + const internal::variable_if_dynamic m_startRow; + const internal::variable_if_dynamic m_startCol; + const internal::variable_if_dynamic m_blockRows; + const internal::variable_if_dynamic m_blockCols; + private: + Index nonZeros() const; +}; + +} // end namespace Eigen + +#endif // EIGEN_SPARSE_BLOCK_H diff --git a/eigen/Eigen/src/SparseCore/SparseColEtree.h b/eigen/Eigen/src/SparseCore/SparseColEtree.h new file mode 100644 index 0000000..f8745f4 --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseColEtree.h @@ -0,0 +1,206 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Désiré Nuentsa-Wakam +// +// 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/. + + +/* + + * NOTE: This file is the modified version of sp_coletree.c file in SuperLU + + * -- SuperLU routine (version 3.1) -- + * Univ. of California Berkeley, Xerox Palo Alto Research Center, + * and Lawrence Berkeley National Lab. + * August 1, 2008 + * + * Copyright (c) 1994 by Xerox Corporation. All rights reserved. + * + * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY + * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. + * + * Permission is hereby granted to use or copy this program for any + * purpose, provided the above notices are retained on all copies. + * Permission to modify the code and to distribute modified code is + * granted, provided the above notices are retained, and a notice that + * the code was modified is included with the above copyright notice. + */ +#ifndef SPARSE_COLETREE_H +#define SPARSE_COLETREE_H + +namespace Eigen { + +namespace internal { + +/** Find the root of the tree/set containing the vertex i : Use Path halving */ +template +Index etree_find (Index i, IndexVector& pp) +{ + Index p = pp(i); // Parent + Index gp = pp(p); // Grand parent + while (gp != p) + { + pp(i) = gp; // Parent pointer on find path is changed to former grand parent + i = gp; + p = pp(i); + gp = pp(p); + } + return p; +} + +/** Compute the column elimination tree of a sparse matrix + * \param mat The matrix in column-major format. + * \param parent The elimination tree + * \param firstRowElt The column index of the first element in each row + * \param perm The permutation to apply to the column of \b mat + */ +template +int coletree(const MatrixType& mat, IndexVector& parent, IndexVector& firstRowElt, typename MatrixType::Index *perm=0) +{ + typedef typename MatrixType::Index Index; + Index nc = mat.cols(); // Number of columns + Index m = mat.rows(); + Index diagSize = (std::min)(nc,m); + IndexVector root(nc); // root of subtree of etree + root.setZero(); + IndexVector pp(nc); // disjoint sets + pp.setZero(); // Initialize disjoint sets + parent.resize(mat.cols()); + //Compute first nonzero column in each row + Index row,col; + firstRowElt.resize(m); + firstRowElt.setConstant(nc); + firstRowElt.segment(0, diagSize).setLinSpaced(diagSize, 0, diagSize-1); + bool found_diag; + for (col = 0; col < nc; col++) + { + Index pcol = col; + if(perm) pcol = perm[col]; + for (typename MatrixType::InnerIterator it(mat, pcol); it; ++it) + { + row = it.row(); + firstRowElt(row) = (std::min)(firstRowElt(row), col); + } + } + /* Compute etree by Liu's algorithm for symmetric matrices, + except use (firstRowElt[r],c) in place of an edge (r,c) of A. + Thus each row clique in A'*A is replaced by a star + centered at its first vertex, which has the same fill. */ + Index rset, cset, rroot; + for (col = 0; col < nc; col++) + { + found_diag = col>=m; + pp(col) = col; + cset = col; + root(cset) = col; + parent(col) = nc; + /* The diagonal element is treated here even if it does not exist in the matrix + * hence the loop is executed once more */ + Index pcol = col; + if(perm) pcol = perm[col]; + for (typename MatrixType::InnerIterator it(mat, pcol); it||!found_diag; ++it) + { // A sequence of interleaved find and union is performed + Index i = col; + if(it) i = it.index(); + if (i == col) found_diag = true; + + row = firstRowElt(i); + if (row >= col) continue; + rset = internal::etree_find(row, pp); // Find the name of the set containing row + rroot = root(rset); + if (rroot != col) + { + parent(rroot) = col; + pp(cset) = rset; + cset = rset; + root(cset) = col; + } + } + } + return 0; +} + +/** + * Depth-first search from vertex n. No recursion. + * This routine was contributed by Cédric Doucet, CEDRAT Group, Meylan, France. +*/ +template +void nr_etdfs (Index n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, Index postnum) +{ + Index current = n, first, next; + while (postnum != n) + { + // No kid for the current node + first = first_kid(current); + + // no kid for the current node + if (first == -1) + { + // Numbering this node because it has no kid + post(current) = postnum++; + + // looking for the next kid + next = next_kid(current); + while (next == -1) + { + // No more kids : back to the parent node + current = parent(current); + // numbering the parent node + post(current) = postnum++; + + // Get the next kid + next = next_kid(current); + } + // stopping criterion + if (postnum == n+1) return; + + // Updating current node + current = next; + } + else + { + current = first; + } + } +} + + +/** + * \brief Post order a tree + * \param n the number of nodes + * \param parent Input tree + * \param post postordered tree + */ +template +void treePostorder(Index n, IndexVector& parent, IndexVector& post) +{ + IndexVector first_kid, next_kid; // Linked list of children + Index postnum; + // Allocate storage for working arrays and results + first_kid.resize(n+1); + next_kid.setZero(n+1); + post.setZero(n+1); + + // Set up structure describing children + Index v, dad; + first_kid.setConstant(-1); + for (v = n-1; v >= 0; v--) + { + dad = parent(v); + next_kid(v) = first_kid(dad); + first_kid(dad) = v; + } + + // Depth-first search from dummy root vertex #n + postnum = 0; + internal::nr_etdfs(n, parent, first_kid, next_kid, post, postnum); +} + +} // end namespace internal + +} // end namespace Eigen + +#endif // SPARSE_COLETREE_H diff --git a/eigen/Eigen/src/SparseCore/SparseCwiseBinaryOp.h b/eigen/Eigen/src/SparseCore/SparseCwiseBinaryOp.h new file mode 100644 index 0000000..5462737 --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseCwiseBinaryOp.h @@ -0,0 +1,324 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// +// 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_SPARSE_CWISE_BINARY_OP_H +#define EIGEN_SPARSE_CWISE_BINARY_OP_H + +namespace Eigen { + +// Here we have to handle 3 cases: +// 1 - sparse op dense +// 2 - dense op sparse +// 3 - sparse op sparse +// We also need to implement a 4th iterator for: +// 4 - dense op dense +// Finally, we also need to distinguish between the product and other operations : +// configuration returned mode +// 1 - sparse op dense product sparse +// generic dense +// 2 - dense op sparse product sparse +// generic dense +// 3 - sparse op sparse product sparse +// generic sparse +// 4 - dense op dense product dense +// generic dense + +namespace internal { + +template<> struct promote_storage_type +{ typedef Sparse ret; }; + +template<> struct promote_storage_type +{ typedef Sparse ret; }; + +template::StorageKind, + typename _RhsStorageMode = typename traits::StorageKind> +class sparse_cwise_binary_op_inner_iterator_selector; + +} // end namespace internal + +template +class CwiseBinaryOpImpl + : public SparseMatrixBase > +{ + public: + class InnerIterator; + class ReverseInnerIterator; + typedef CwiseBinaryOp Derived; + EIGEN_SPARSE_PUBLIC_INTERFACE(Derived) + CwiseBinaryOpImpl() + { + EIGEN_STATIC_ASSERT(( + (!internal::is_same::StorageKind, + typename internal::traits::StorageKind>::value) + || ((Lhs::Flags&RowMajorBit) == (Rhs::Flags&RowMajorBit))), + THE_STORAGE_ORDER_OF_BOTH_SIDES_MUST_MATCH); + } +}; + +template +class CwiseBinaryOpImpl::InnerIterator + : public internal::sparse_cwise_binary_op_inner_iterator_selector::InnerIterator> +{ + public: + typedef typename Lhs::Index Index; + typedef internal::sparse_cwise_binary_op_inner_iterator_selector< + BinaryOp,Lhs,Rhs, InnerIterator> Base; + + // NOTE: we have to prefix Index by "typename Lhs::" to avoid an ICE with VC11 + EIGEN_STRONG_INLINE InnerIterator(const CwiseBinaryOpImpl& binOp, typename Lhs::Index outer) + : Base(binOp.derived(),outer) + {} +}; + +/*************************************************************************** +* Implementation of inner-iterators +***************************************************************************/ + +// template struct internal::func_is_conjunction { enum { ret = false }; }; +// template struct internal::func_is_conjunction > { enum { ret = true }; }; + +// TODO generalize the internal::scalar_product_op specialization to all conjunctions if any ! + +namespace internal { + +// sparse - sparse (generic) +template +class sparse_cwise_binary_op_inner_iterator_selector +{ + typedef CwiseBinaryOp CwiseBinaryXpr; + typedef typename traits::Scalar Scalar; + typedef typename traits::_LhsNested _LhsNested; + typedef typename traits::_RhsNested _RhsNested; + typedef typename _LhsNested::InnerIterator LhsIterator; + typedef typename _RhsNested::InnerIterator RhsIterator; + typedef typename Lhs::Index Index; + + public: + + EIGEN_STRONG_INLINE sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, Index outer) + : m_lhsIter(xpr.lhs(),outer), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor()) + { + this->operator++(); + } + + EIGEN_STRONG_INLINE Derived& operator++() + { + if (m_lhsIter && m_rhsIter && (m_lhsIter.index() == m_rhsIter.index())) + { + m_id = m_lhsIter.index(); + m_value = m_functor(m_lhsIter.value(), m_rhsIter.value()); + ++m_lhsIter; + ++m_rhsIter; + } + else if (m_lhsIter && (!m_rhsIter || (m_lhsIter.index() < m_rhsIter.index()))) + { + m_id = m_lhsIter.index(); + m_value = m_functor(m_lhsIter.value(), Scalar(0)); + ++m_lhsIter; + } + else if (m_rhsIter && (!m_lhsIter || (m_lhsIter.index() > m_rhsIter.index()))) + { + m_id = m_rhsIter.index(); + m_value = m_functor(Scalar(0), m_rhsIter.value()); + ++m_rhsIter; + } + else + { + m_value = 0; // this is to avoid a compilation warning + m_id = -1; + } + return *static_cast(this); + } + + EIGEN_STRONG_INLINE Scalar value() const { return m_value; } + + EIGEN_STRONG_INLINE Index index() const { return m_id; } + EIGEN_STRONG_INLINE Index row() const { return Lhs::IsRowMajor ? m_lhsIter.row() : index(); } + EIGEN_STRONG_INLINE Index col() const { return Lhs::IsRowMajor ? index() : m_lhsIter.col(); } + + EIGEN_STRONG_INLINE operator bool() const { return m_id>=0; } + + protected: + LhsIterator m_lhsIter; + RhsIterator m_rhsIter; + const BinaryOp& m_functor; + Scalar m_value; + Index m_id; +}; + +// sparse - sparse (product) +template +class sparse_cwise_binary_op_inner_iterator_selector, Lhs, Rhs, Derived, Sparse, Sparse> +{ + typedef scalar_product_op BinaryFunc; + typedef CwiseBinaryOp CwiseBinaryXpr; + typedef typename CwiseBinaryXpr::Scalar Scalar; + typedef typename traits::_LhsNested _LhsNested; + typedef typename _LhsNested::InnerIterator LhsIterator; + typedef typename traits::_RhsNested _RhsNested; + typedef typename _RhsNested::InnerIterator RhsIterator; + typedef typename Lhs::Index Index; + public: + + EIGEN_STRONG_INLINE sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, Index outer) + : m_lhsIter(xpr.lhs(),outer), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor()) + { + while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index())) + { + if (m_lhsIter.index() < m_rhsIter.index()) + ++m_lhsIter; + else + ++m_rhsIter; + } + } + + EIGEN_STRONG_INLINE Derived& operator++() + { + ++m_lhsIter; + ++m_rhsIter; + while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index())) + { + if (m_lhsIter.index() < m_rhsIter.index()) + ++m_lhsIter; + else + ++m_rhsIter; + } + return *static_cast(this); + } + + EIGEN_STRONG_INLINE Scalar value() const { return m_functor(m_lhsIter.value(), m_rhsIter.value()); } + + EIGEN_STRONG_INLINE Index index() const { return m_lhsIter.index(); } + EIGEN_STRONG_INLINE Index row() const { return m_lhsIter.row(); } + EIGEN_STRONG_INLINE Index col() const { return m_lhsIter.col(); } + + EIGEN_STRONG_INLINE operator bool() const { return (m_lhsIter && m_rhsIter); } + + protected: + LhsIterator m_lhsIter; + RhsIterator m_rhsIter; + const BinaryFunc& m_functor; +}; + +// sparse - dense (product) +template +class sparse_cwise_binary_op_inner_iterator_selector, Lhs, Rhs, Derived, Sparse, Dense> +{ + typedef scalar_product_op BinaryFunc; + typedef CwiseBinaryOp CwiseBinaryXpr; + typedef typename CwiseBinaryXpr::Scalar Scalar; + typedef typename traits::_LhsNested _LhsNested; + typedef typename traits::RhsNested RhsNested; + typedef typename _LhsNested::InnerIterator LhsIterator; + typedef typename Lhs::Index Index; + enum { IsRowMajor = (int(Lhs::Flags)&RowMajorBit)==RowMajorBit }; + public: + + EIGEN_STRONG_INLINE sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, Index outer) + : m_rhs(xpr.rhs()), m_lhsIter(xpr.lhs(),outer), m_functor(xpr.functor()), m_outer(outer) + {} + + EIGEN_STRONG_INLINE Derived& operator++() + { + ++m_lhsIter; + return *static_cast(this); + } + + EIGEN_STRONG_INLINE Scalar value() const + { return m_functor(m_lhsIter.value(), + m_rhs.coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); } + + EIGEN_STRONG_INLINE Index index() const { return m_lhsIter.index(); } + EIGEN_STRONG_INLINE Index row() const { return m_lhsIter.row(); } + EIGEN_STRONG_INLINE Index col() const { return m_lhsIter.col(); } + + EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter; } + + protected: + RhsNested m_rhs; + LhsIterator m_lhsIter; + const BinaryFunc m_functor; + const Index m_outer; +}; + +// sparse - dense (product) +template +class sparse_cwise_binary_op_inner_iterator_selector, Lhs, Rhs, Derived, Dense, Sparse> +{ + typedef scalar_product_op BinaryFunc; + typedef CwiseBinaryOp CwiseBinaryXpr; + typedef typename CwiseBinaryXpr::Scalar Scalar; + typedef typename traits::_RhsNested _RhsNested; + typedef typename _RhsNested::InnerIterator RhsIterator; + typedef typename Lhs::Index Index; + + enum { IsRowMajor = (int(Rhs::Flags)&RowMajorBit)==RowMajorBit }; + public: + + EIGEN_STRONG_INLINE sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, Index outer) + : m_xpr(xpr), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor()), m_outer(outer) + {} + + EIGEN_STRONG_INLINE Derived& operator++() + { + ++m_rhsIter; + return *static_cast(this); + } + + EIGEN_STRONG_INLINE Scalar value() const + { return m_functor(m_xpr.lhs().coeff(IsRowMajor?m_outer:m_rhsIter.index(),IsRowMajor?m_rhsIter.index():m_outer), m_rhsIter.value()); } + + EIGEN_STRONG_INLINE Index index() const { return m_rhsIter.index(); } + EIGEN_STRONG_INLINE Index row() const { return m_rhsIter.row(); } + EIGEN_STRONG_INLINE Index col() const { return m_rhsIter.col(); } + + EIGEN_STRONG_INLINE operator bool() const { return m_rhsIter; } + + protected: + const CwiseBinaryXpr& m_xpr; + RhsIterator m_rhsIter; + const BinaryFunc& m_functor; + const Index m_outer; +}; + +} // end namespace internal + +/*************************************************************************** +* Implementation of SparseMatrixBase and SparseCwise functions/operators +***************************************************************************/ + +template +template +EIGEN_STRONG_INLINE Derived & +SparseMatrixBase::operator-=(const SparseMatrixBase &other) +{ + return derived() = derived() - other.derived(); +} + +template +template +EIGEN_STRONG_INLINE Derived & +SparseMatrixBase::operator+=(const SparseMatrixBase& other) +{ + return derived() = derived() + other.derived(); +} + +template +template +EIGEN_STRONG_INLINE const typename SparseMatrixBase::template CwiseProductDenseReturnType::Type +SparseMatrixBase::cwiseProduct(const MatrixBase &other) const +{ + return typename CwiseProductDenseReturnType::Type(derived(), other.derived()); +} + +} // end namespace Eigen + +#endif // EIGEN_SPARSE_CWISE_BINARY_OP_H diff --git a/eigen/Eigen/src/SparseCore/SparseCwiseUnaryOp.h b/eigen/Eigen/src/SparseCore/SparseCwiseUnaryOp.h new file mode 100644 index 0000000..5a50c78 --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseCwiseUnaryOp.h @@ -0,0 +1,163 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2010 Gael Guennebaud +// +// 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_SPARSE_CWISE_UNARY_OP_H +#define EIGEN_SPARSE_CWISE_UNARY_OP_H + +namespace Eigen { + +template +class CwiseUnaryOpImpl + : public SparseMatrixBase > +{ + public: + + class InnerIterator; + class ReverseInnerIterator; + + typedef CwiseUnaryOp Derived; + EIGEN_SPARSE_PUBLIC_INTERFACE(Derived) + + protected: + typedef typename internal::traits::_XprTypeNested _MatrixTypeNested; + typedef typename _MatrixTypeNested::InnerIterator MatrixTypeIterator; + typedef typename _MatrixTypeNested::ReverseInnerIterator MatrixTypeReverseIterator; +}; + +template +class CwiseUnaryOpImpl::InnerIterator + : public CwiseUnaryOpImpl::MatrixTypeIterator +{ + typedef typename CwiseUnaryOpImpl::Scalar Scalar; + typedef typename CwiseUnaryOpImpl::MatrixTypeIterator Base; + public: + + EIGEN_STRONG_INLINE InnerIterator(const CwiseUnaryOpImpl& unaryOp, typename CwiseUnaryOpImpl::Index outer) + : Base(unaryOp.derived().nestedExpression(),outer), m_functor(unaryOp.derived().functor()) + {} + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { Base::operator++(); return *this; } + + EIGEN_STRONG_INLINE typename CwiseUnaryOpImpl::Scalar value() const { return m_functor(Base::value()); } + + protected: + const UnaryOp m_functor; + private: + typename CwiseUnaryOpImpl::Scalar& valueRef(); +}; + +template +class CwiseUnaryOpImpl::ReverseInnerIterator + : public CwiseUnaryOpImpl::MatrixTypeReverseIterator +{ + typedef typename CwiseUnaryOpImpl::Scalar Scalar; + typedef typename CwiseUnaryOpImpl::MatrixTypeReverseIterator Base; + public: + + EIGEN_STRONG_INLINE ReverseInnerIterator(const CwiseUnaryOpImpl& unaryOp, typename CwiseUnaryOpImpl::Index outer) + : Base(unaryOp.derived().nestedExpression(),outer), m_functor(unaryOp.derived().functor()) + {} + + EIGEN_STRONG_INLINE ReverseInnerIterator& operator--() + { Base::operator--(); return *this; } + + EIGEN_STRONG_INLINE typename CwiseUnaryOpImpl::Scalar value() const { return m_functor(Base::value()); } + + protected: + const UnaryOp m_functor; + private: + typename CwiseUnaryOpImpl::Scalar& valueRef(); +}; + +template +class CwiseUnaryViewImpl + : public SparseMatrixBase > +{ + public: + + class InnerIterator; + class ReverseInnerIterator; + + typedef CwiseUnaryView Derived; + EIGEN_SPARSE_PUBLIC_INTERFACE(Derived) + + protected: + typedef typename internal::traits::_MatrixTypeNested _MatrixTypeNested; + typedef typename _MatrixTypeNested::InnerIterator MatrixTypeIterator; + typedef typename _MatrixTypeNested::ReverseInnerIterator MatrixTypeReverseIterator; +}; + +template +class CwiseUnaryViewImpl::InnerIterator + : public CwiseUnaryViewImpl::MatrixTypeIterator +{ + typedef typename CwiseUnaryViewImpl::Scalar Scalar; + typedef typename CwiseUnaryViewImpl::MatrixTypeIterator Base; + public: + + EIGEN_STRONG_INLINE InnerIterator(const CwiseUnaryViewImpl& unaryOp, typename CwiseUnaryViewImpl::Index outer) + : Base(unaryOp.derived().nestedExpression(),outer), m_functor(unaryOp.derived().functor()) + {} + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { Base::operator++(); return *this; } + + EIGEN_STRONG_INLINE typename CwiseUnaryViewImpl::Scalar value() const { return m_functor(Base::value()); } + EIGEN_STRONG_INLINE typename CwiseUnaryViewImpl::Scalar& valueRef() { return m_functor(Base::valueRef()); } + + protected: + const ViewOp m_functor; +}; + +template +class CwiseUnaryViewImpl::ReverseInnerIterator + : public CwiseUnaryViewImpl::MatrixTypeReverseIterator +{ + typedef typename CwiseUnaryViewImpl::Scalar Scalar; + typedef typename CwiseUnaryViewImpl::MatrixTypeReverseIterator Base; + public: + + EIGEN_STRONG_INLINE ReverseInnerIterator(const CwiseUnaryViewImpl& unaryOp, typename CwiseUnaryViewImpl::Index outer) + : Base(unaryOp.derived().nestedExpression(),outer), m_functor(unaryOp.derived().functor()) + {} + + EIGEN_STRONG_INLINE ReverseInnerIterator& operator--() + { Base::operator--(); return *this; } + + EIGEN_STRONG_INLINE typename CwiseUnaryViewImpl::Scalar value() const { return m_functor(Base::value()); } + EIGEN_STRONG_INLINE typename CwiseUnaryViewImpl::Scalar& valueRef() { return m_functor(Base::valueRef()); } + + protected: + const ViewOp m_functor; +}; + +template +EIGEN_STRONG_INLINE Derived& +SparseMatrixBase::operator*=(const Scalar& other) +{ + for (Index j=0; j +EIGEN_STRONG_INLINE Derived& +SparseMatrixBase::operator/=(const Scalar& other) +{ + for (Index j=0; j +// +// 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_SPARSEDENSEPRODUCT_H +#define EIGEN_SPARSEDENSEPRODUCT_H + +namespace Eigen { + +template struct SparseDenseProductReturnType +{ + typedef SparseTimeDenseProduct Type; +}; + +template struct SparseDenseProductReturnType +{ + typedef typename internal::conditional< + Lhs::IsRowMajor, + SparseDenseOuterProduct, + SparseDenseOuterProduct >::type Type; +}; + +template struct DenseSparseProductReturnType +{ + typedef DenseTimeSparseProduct Type; +}; + +template struct DenseSparseProductReturnType +{ + typedef typename internal::conditional< + Rhs::IsRowMajor, + SparseDenseOuterProduct, + SparseDenseOuterProduct >::type Type; +}; + +namespace internal { + +template +struct traits > +{ + typedef Sparse StorageKind; + typedef typename scalar_product_traits::Scalar, + typename traits::Scalar>::ReturnType Scalar; + typedef typename Lhs::Index Index; + typedef typename Lhs::Nested LhsNested; + typedef typename Rhs::Nested RhsNested; + typedef typename remove_all::type _LhsNested; + typedef typename remove_all::type _RhsNested; + + enum { + LhsCoeffReadCost = traits<_LhsNested>::CoeffReadCost, + RhsCoeffReadCost = traits<_RhsNested>::CoeffReadCost, + + RowsAtCompileTime = Tr ? int(traits::RowsAtCompileTime) : int(traits::RowsAtCompileTime), + ColsAtCompileTime = Tr ? int(traits::ColsAtCompileTime) : int(traits::ColsAtCompileTime), + MaxRowsAtCompileTime = Tr ? int(traits::MaxRowsAtCompileTime) : int(traits::MaxRowsAtCompileTime), + MaxColsAtCompileTime = Tr ? int(traits::MaxColsAtCompileTime) : int(traits::MaxColsAtCompileTime), + + Flags = Tr ? RowMajorBit : 0, + + CoeffReadCost = LhsCoeffReadCost + RhsCoeffReadCost + NumTraits::MulCost + }; +}; + +} // end namespace internal + +template +class SparseDenseOuterProduct + : public SparseMatrixBase > +{ + public: + + typedef SparseMatrixBase Base; + EIGEN_DENSE_PUBLIC_INTERFACE(SparseDenseOuterProduct) + typedef internal::traits Traits; + + private: + + typedef typename Traits::LhsNested LhsNested; + typedef typename Traits::RhsNested RhsNested; + typedef typename Traits::_LhsNested _LhsNested; + typedef typename Traits::_RhsNested _RhsNested; + + public: + + class InnerIterator; + + EIGEN_STRONG_INLINE SparseDenseOuterProduct(const Lhs& lhs, const Rhs& rhs) + : m_lhs(lhs), m_rhs(rhs) + { + EIGEN_STATIC_ASSERT(!Tr,YOU_MADE_A_PROGRAMMING_MISTAKE); + } + + EIGEN_STRONG_INLINE SparseDenseOuterProduct(const Rhs& rhs, const Lhs& lhs) + : m_lhs(lhs), m_rhs(rhs) + { + EIGEN_STATIC_ASSERT(Tr,YOU_MADE_A_PROGRAMMING_MISTAKE); + } + + EIGEN_STRONG_INLINE Index rows() const { return Tr ? m_rhs.rows() : m_lhs.rows(); } + EIGEN_STRONG_INLINE Index cols() const { return Tr ? m_lhs.cols() : m_rhs.cols(); } + + EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; } + EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; } + + protected: + LhsNested m_lhs; + RhsNested m_rhs; +}; + +template +class SparseDenseOuterProduct::InnerIterator : public _LhsNested::InnerIterator +{ + typedef typename _LhsNested::InnerIterator Base; + typedef typename SparseDenseOuterProduct::Index Index; + public: + EIGEN_STRONG_INLINE InnerIterator(const SparseDenseOuterProduct& prod, Index outer) + : Base(prod.lhs(), 0), m_outer(outer), m_factor(get(prod.rhs(), outer, typename internal::traits::StorageKind() )) + { } + + inline Index outer() const { return m_outer; } + inline Index row() const { return Transpose ? m_outer : Base::index(); } + inline Index col() const { return Transpose ? Base::index() : m_outer; } + + inline Scalar value() const { return Base::value() * m_factor; } + + protected: + static Scalar get(const _RhsNested &rhs, Index outer, Dense = Dense()) + { + return rhs.coeff(outer); + } + + static Scalar get(const _RhsNested &rhs, Index outer, Sparse = Sparse()) + { + typename Traits::_RhsNested::InnerIterator it(rhs, outer); + if (it && it.index()==0) + return it.value(); + + return Scalar(0); + } + + Index m_outer; + Scalar m_factor; +}; + +namespace internal { +template +struct traits > + : traits, Lhs, Rhs> > +{ + typedef Dense StorageKind; + typedef MatrixXpr XprKind; +}; + +template +struct sparse_time_dense_product_impl; + +template +struct sparse_time_dense_product_impl +{ + typedef typename internal::remove_all::type Lhs; + typedef typename internal::remove_all::type Rhs; + typedef typename internal::remove_all::type Res; + typedef typename Lhs::Index Index; + typedef typename Lhs::InnerIterator LhsInnerIterator; + static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha) + { + for(Index c=0; c +struct sparse_time_dense_product_impl +{ + typedef typename internal::remove_all::type Lhs; + typedef typename internal::remove_all::type Rhs; + typedef typename internal::remove_all::type Res; + typedef typename Lhs::InnerIterator LhsInnerIterator; + typedef typename Lhs::Index Index; + static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha) + { + for(Index c=0; c +struct sparse_time_dense_product_impl +{ + typedef typename internal::remove_all::type Lhs; + typedef typename internal::remove_all::type Rhs; + typedef typename internal::remove_all::type Res; + typedef typename Lhs::InnerIterator LhsInnerIterator; + typedef typename Lhs::Index Index; + static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha) + { + for(Index j=0; j +struct sparse_time_dense_product_impl +{ + typedef typename internal::remove_all::type Lhs; + typedef typename internal::remove_all::type Rhs; + typedef typename internal::remove_all::type Res; + typedef typename Lhs::InnerIterator LhsInnerIterator; + typedef typename Lhs::Index Index; + static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha) + { + for(Index j=0; j +inline void sparse_time_dense_product(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha) +{ + sparse_time_dense_product_impl::run(lhs, rhs, res, alpha); +} + +} // end namespace internal + +template +class SparseTimeDenseProduct + : public ProductBase, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseTimeDenseProduct) + + SparseTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + {} + + template void scaleAndAddTo(Dest& dest, const Scalar& alpha) const + { + internal::sparse_time_dense_product(m_lhs, m_rhs, dest, alpha); + } + + private: + SparseTimeDenseProduct& operator=(const SparseTimeDenseProduct&); +}; + + +// dense = dense * sparse +namespace internal { +template +struct traits > + : traits, Lhs, Rhs> > +{ + typedef Dense StorageKind; +}; +} // end namespace internal + +template +class DenseTimeSparseProduct + : public ProductBase, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseProduct) + + DenseTimeSparseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + {} + + template void scaleAndAddTo(Dest& dest, const Scalar& alpha) const + { + Transpose lhs_t(m_lhs); + Transpose rhs_t(m_rhs); + Transpose dest_t(dest); + internal::sparse_time_dense_product(rhs_t, lhs_t, dest_t, alpha); + } + + private: + DenseTimeSparseProduct& operator=(const DenseTimeSparseProduct&); +}; + +} // end namespace Eigen + +#endif // EIGEN_SPARSEDENSEPRODUCT_H diff --git a/eigen/Eigen/src/SparseCore/SparseDiagonalProduct.h b/eigen/Eigen/src/SparseCore/SparseDiagonalProduct.h new file mode 100644 index 0000000..1bb590e --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseDiagonalProduct.h @@ -0,0 +1,196 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud +// +// 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_SPARSE_DIAGONAL_PRODUCT_H +#define EIGEN_SPARSE_DIAGONAL_PRODUCT_H + +namespace Eigen { + +// The product of a diagonal matrix with a sparse matrix can be easily +// implemented using expression template. +// We have two consider very different cases: +// 1 - diag * row-major sparse +// => each inner vector <=> scalar * sparse vector product +// => so we can reuse CwiseUnaryOp::InnerIterator +// 2 - diag * col-major sparse +// => each inner vector <=> densevector * sparse vector cwise product +// => again, we can reuse specialization of CwiseBinaryOp::InnerIterator +// for that particular case +// The two other cases are symmetric. + +namespace internal { + +template +struct traits > +{ + typedef typename remove_all::type _Lhs; + typedef typename remove_all::type _Rhs; + typedef typename _Lhs::Scalar Scalar; + typedef typename promote_index_type::Index, + typename traits::Index>::type Index; + typedef Sparse StorageKind; + typedef MatrixXpr XprKind; + enum { + RowsAtCompileTime = _Lhs::RowsAtCompileTime, + ColsAtCompileTime = _Rhs::ColsAtCompileTime, + + MaxRowsAtCompileTime = _Lhs::MaxRowsAtCompileTime, + MaxColsAtCompileTime = _Rhs::MaxColsAtCompileTime, + + SparseFlags = is_diagonal<_Lhs>::ret ? int(_Rhs::Flags) : int(_Lhs::Flags), + Flags = (SparseFlags&RowMajorBit), + CoeffReadCost = Dynamic + }; +}; + +enum {SDP_IsDiagonal, SDP_IsSparseRowMajor, SDP_IsSparseColMajor}; +template +class sparse_diagonal_product_inner_iterator_selector; + +} // end namespace internal + +template +class SparseDiagonalProduct + : public SparseMatrixBase >, + internal::no_assignment_operator +{ + typedef typename Lhs::Nested LhsNested; + typedef typename Rhs::Nested RhsNested; + + typedef typename internal::remove_all::type _LhsNested; + typedef typename internal::remove_all::type _RhsNested; + + enum { + LhsMode = internal::is_diagonal<_LhsNested>::ret ? internal::SDP_IsDiagonal + : (_LhsNested::Flags&RowMajorBit) ? internal::SDP_IsSparseRowMajor : internal::SDP_IsSparseColMajor, + RhsMode = internal::is_diagonal<_RhsNested>::ret ? internal::SDP_IsDiagonal + : (_RhsNested::Flags&RowMajorBit) ? internal::SDP_IsSparseRowMajor : internal::SDP_IsSparseColMajor + }; + + public: + + EIGEN_SPARSE_PUBLIC_INTERFACE(SparseDiagonalProduct) + + typedef internal::sparse_diagonal_product_inner_iterator_selector + <_LhsNested,_RhsNested,SparseDiagonalProduct,LhsMode,RhsMode> InnerIterator; + + // We do not want ReverseInnerIterator for diagonal-sparse products, + // but this dummy declaration is needed to make diag * sparse * diag compile. + class ReverseInnerIterator; + + EIGEN_STRONG_INLINE SparseDiagonalProduct(const Lhs& lhs, const Rhs& rhs) + : m_lhs(lhs), m_rhs(rhs) + { + eigen_assert(lhs.cols() == rhs.rows() && "invalid sparse matrix * diagonal matrix product"); + } + + EIGEN_STRONG_INLINE Index rows() const { return m_lhs.rows(); } + EIGEN_STRONG_INLINE Index cols() const { return m_rhs.cols(); } + + EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; } + EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; } + + protected: + LhsNested m_lhs; + RhsNested m_rhs; +}; + +namespace internal { + +template +class sparse_diagonal_product_inner_iterator_selector + + : public CwiseUnaryOp,const Rhs>::InnerIterator +{ + typedef typename CwiseUnaryOp,const Rhs>::InnerIterator Base; + typedef typename Lhs::Index Index; + public: + inline sparse_diagonal_product_inner_iterator_selector( + const SparseDiagonalProductType& expr, Index outer) + : Base(expr.rhs()*(expr.lhs().diagonal().coeff(outer)), outer) + {} +}; + +template +class sparse_diagonal_product_inner_iterator_selector + + : public CwiseBinaryOp< + scalar_product_op, + const typename Rhs::ConstInnerVectorReturnType, + const typename Lhs::DiagonalVectorType>::InnerIterator +{ + typedef typename CwiseBinaryOp< + scalar_product_op, + const typename Rhs::ConstInnerVectorReturnType, + const typename Lhs::DiagonalVectorType>::InnerIterator Base; + typedef typename Lhs::Index Index; + Index m_outer; + public: + inline sparse_diagonal_product_inner_iterator_selector( + const SparseDiagonalProductType& expr, Index outer) + : Base(expr.rhs().innerVector(outer) .cwiseProduct(expr.lhs().diagonal()), 0), m_outer(outer) + {} + + inline Index outer() const { return m_outer; } + inline Index col() const { return m_outer; } +}; + +template +class sparse_diagonal_product_inner_iterator_selector + + : public CwiseUnaryOp,const Lhs>::InnerIterator +{ + typedef typename CwiseUnaryOp,const Lhs>::InnerIterator Base; + typedef typename Lhs::Index Index; + public: + inline sparse_diagonal_product_inner_iterator_selector( + const SparseDiagonalProductType& expr, Index outer) + : Base(expr.lhs()*expr.rhs().diagonal().coeff(outer), outer) + {} +}; + +template +class sparse_diagonal_product_inner_iterator_selector + + : public CwiseBinaryOp< + scalar_product_op, + const typename Lhs::ConstInnerVectorReturnType, + const Transpose >::InnerIterator +{ + typedef typename CwiseBinaryOp< + scalar_product_op, + const typename Lhs::ConstInnerVectorReturnType, + const Transpose >::InnerIterator Base; + typedef typename Lhs::Index Index; + Index m_outer; + public: + inline sparse_diagonal_product_inner_iterator_selector( + const SparseDiagonalProductType& expr, Index outer) + : Base(expr.lhs().innerVector(outer) .cwiseProduct(expr.rhs().diagonal().transpose()), 0), m_outer(outer) + {} + + inline Index outer() const { return m_outer; } + inline Index row() const { return m_outer; } +}; + +} // end namespace internal + +// SparseMatrixBase functions + +template +template +const SparseDiagonalProduct +SparseMatrixBase::operator*(const DiagonalBase &other) const +{ + return SparseDiagonalProduct(this->derived(), other.derived()); +} + +} // end namespace Eigen + +#endif // EIGEN_SPARSE_DIAGONAL_PRODUCT_H diff --git a/eigen/Eigen/src/SparseCore/SparseDot.h b/eigen/Eigen/src/SparseCore/SparseDot.h new file mode 100644 index 0000000..db39c9a --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseDot.h @@ -0,0 +1,101 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// +// 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_SPARSE_DOT_H +#define EIGEN_SPARSE_DOT_H + +namespace Eigen { + +template +template +typename internal::traits::Scalar +SparseMatrixBase::dot(const MatrixBase& other) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) + EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived) + 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) + + eigen_assert(size() == other.size()); + eigen_assert(other.size()>0 && "you are using a non initialized vector"); + + typename Derived::InnerIterator i(derived(),0); + Scalar res(0); + while (i) + { + res += numext::conj(i.value()) * other.coeff(i.index()); + ++i; + } + return res; +} + +template +template +typename internal::traits::Scalar +SparseMatrixBase::dot(const SparseMatrixBase& other) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived) + EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived) + 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) + + eigen_assert(size() == other.size()); + + typedef typename Derived::Nested Nested; + typedef typename OtherDerived::Nested OtherNested; + typedef typename internal::remove_all::type NestedCleaned; + typedef typename internal::remove_all::type OtherNestedCleaned; + + Nested nthis(derived()); + OtherNested nother(other.derived()); + + typename NestedCleaned::InnerIterator i(nthis,0); + typename OtherNestedCleaned::InnerIterator j(nother,0); + Scalar res(0); + while (i && j) + { + if (i.index()==j.index()) + { + res += numext::conj(i.value()) * j.value(); + ++i; ++j; + } + else if (i.index() +inline typename NumTraits::Scalar>::Real +SparseMatrixBase::squaredNorm() const +{ + return numext::real((*this).cwiseAbs2().sum()); +} + +template +inline typename NumTraits::Scalar>::Real +SparseMatrixBase::norm() const +{ + using std::sqrt; + return sqrt(squaredNorm()); +} + +template +inline typename NumTraits::Scalar>::Real +SparseMatrixBase::blueNorm() const +{ + return internal::blueNorm_impl(*this); +} +} // end namespace Eigen + +#endif // EIGEN_SPARSE_DOT_H diff --git a/eigen/Eigen/src/SparseCore/SparseFuzzy.h b/eigen/Eigen/src/SparseCore/SparseFuzzy.h new file mode 100644 index 0000000..45f36e9 --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseFuzzy.h @@ -0,0 +1,26 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// +// 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_SPARSE_FUZZY_H +#define EIGEN_SPARSE_FUZZY_H + +// template +// template +// bool SparseMatrixBase::isApprox( +// const OtherDerived& other, +// typename NumTraits::Real prec +// ) const +// { +// const typename internal::nested::type nested(derived()); +// const typename internal::nested::type otherNested(other.derived()); +// return (nested - otherNested).cwise().abs2().sum() +// <= prec * prec * (std::min)(nested.cwise().abs2().sum(), otherNested.cwise().abs2().sum()); +// } + +#endif // EIGEN_SPARSE_FUZZY_H diff --git a/eigen/Eigen/src/SparseCore/SparseMatrix.h b/eigen/Eigen/src/SparseCore/SparseMatrix.h new file mode 100644 index 0000000..3b8946a --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseMatrix.h @@ -0,0 +1,1262 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2010 Gael Guennebaud +// +// 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_SPARSEMATRIX_H +#define EIGEN_SPARSEMATRIX_H + +namespace Eigen { + +/** \ingroup SparseCore_Module + * + * \class SparseMatrix + * + * \brief A versatible sparse matrix representation + * + * This class implements a more versatile variants of the common \em compressed row/column storage format. + * Each colmun's (resp. row) non zeros are stored as a pair of value with associated row (resp. colmiun) index. + * All the non zeros are stored in a single large buffer. Unlike the \em compressed format, there might be extra + * space inbetween the nonzeros of two successive colmuns (resp. rows) such that insertion of new non-zero + * can be done with limited memory reallocation and copies. + * + * A call to the function makeCompressed() turns the matrix into the standard \em compressed format + * compatible with many library. + * + * More details on this storage sceheme are given in the \ref TutorialSparse "manual pages". + * + * \tparam _Scalar the scalar type, i.e. the type of the coefficients + * \tparam _Options Union of bit flags controlling the storage scheme. Currently the only possibility + * is ColMajor or RowMajor. The default is 0 which means column-major. + * \tparam _Index the type of the indices. It has to be a \b signed type (e.g., short, int, std::ptrdiff_t). Default is \c int. + * + * This class can be extended with the help of the plugin mechanism described on the page + * \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_SPARSEMATRIX_PLUGIN. + */ + +namespace internal { +template +struct traits > +{ + typedef _Scalar Scalar; + typedef _Index Index; + typedef Sparse StorageKind; + typedef MatrixXpr XprKind; + enum { + RowsAtCompileTime = Dynamic, + ColsAtCompileTime = Dynamic, + MaxRowsAtCompileTime = Dynamic, + MaxColsAtCompileTime = Dynamic, + Flags = _Options | NestByRefBit | LvalueBit, + CoeffReadCost = NumTraits::ReadCost, + SupportedAccessPatterns = InnerRandomAccessPattern + }; +}; + +template +struct traits, DiagIndex> > +{ + typedef SparseMatrix<_Scalar, _Options, _Index> MatrixType; + typedef typename nested::type MatrixTypeNested; + typedef typename remove_reference::type _MatrixTypeNested; + + typedef _Scalar Scalar; + typedef Dense StorageKind; + typedef _Index Index; + typedef MatrixXpr XprKind; + + enum { + RowsAtCompileTime = Dynamic, + ColsAtCompileTime = 1, + MaxRowsAtCompileTime = Dynamic, + MaxColsAtCompileTime = 1, + Flags = 0, + CoeffReadCost = _MatrixTypeNested::CoeffReadCost*10 + }; +}; + +} // end namespace internal + +template +class SparseMatrix + : public SparseMatrixBase > +{ + public: + EIGEN_SPARSE_PUBLIC_INTERFACE(SparseMatrix) + EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseMatrix, +=) + EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseMatrix, -=) + + typedef MappedSparseMatrix Map; + using Base::IsRowMajor; + typedef internal::CompressedStorage Storage; + enum { + Options = _Options + }; + + protected: + + typedef SparseMatrix TransposedSparseMatrix; + + Index m_outerSize; + Index m_innerSize; + Index* m_outerIndex; + Index* m_innerNonZeros; // optional, if null then the data is compressed + Storage m_data; + + Eigen::Map > innerNonZeros() { return Eigen::Map >(m_innerNonZeros, m_innerNonZeros?m_outerSize:0); } + const Eigen::Map > innerNonZeros() const { return Eigen::Map >(m_innerNonZeros, m_innerNonZeros?m_outerSize:0); } + + public: + + /** \returns whether \c *this is in compressed form. */ + inline bool isCompressed() const { return m_innerNonZeros==0; } + + /** \returns the number of rows of the matrix */ + inline Index rows() const { return IsRowMajor ? m_outerSize : m_innerSize; } + /** \returns the number of columns of the matrix */ + inline Index cols() const { return IsRowMajor ? m_innerSize : m_outerSize; } + + /** \returns the number of rows (resp. columns) of the matrix if the storage order column major (resp. row major) */ + inline Index innerSize() const { return m_innerSize; } + /** \returns the number of columns (resp. rows) of the matrix if the storage order column major (resp. row major) */ + inline Index outerSize() const { return m_outerSize; } + + /** \returns a const pointer to the array of values. + * This function is aimed at interoperability with other libraries. + * \sa innerIndexPtr(), outerIndexPtr() */ + inline const Scalar* valuePtr() const { return m_data.valuePtr(); } + /** \returns a non-const pointer to the array of values. + * This function is aimed at interoperability with other libraries. + * \sa innerIndexPtr(), outerIndexPtr() */ + inline Scalar* valuePtr() { return m_data.valuePtr(); } + + /** \returns a const pointer to the array of inner indices. + * This function is aimed at interoperability with other libraries. + * \sa valuePtr(), outerIndexPtr() */ + inline const Index* innerIndexPtr() const { return m_data.indexPtr(); } + /** \returns a non-const pointer to the array of inner indices. + * This function is aimed at interoperability with other libraries. + * \sa valuePtr(), outerIndexPtr() */ + inline Index* innerIndexPtr() { return m_data.indexPtr(); } + + /** \returns a const pointer to the array of the starting positions of the inner vectors. + * This function is aimed at interoperability with other libraries. + * \sa valuePtr(), innerIndexPtr() */ + inline const Index* outerIndexPtr() const { return m_outerIndex; } + /** \returns a non-const pointer to the array of the starting positions of the inner vectors. + * This function is aimed at interoperability with other libraries. + * \sa valuePtr(), innerIndexPtr() */ + inline Index* outerIndexPtr() { return m_outerIndex; } + + /** \returns a const pointer to the array of the number of non zeros of the inner vectors. + * This function is aimed at interoperability with other libraries. + * \warning it returns the null pointer 0 in compressed mode */ + inline const Index* innerNonZeroPtr() const { return m_innerNonZeros; } + /** \returns a non-const pointer to the array of the number of non zeros of the inner vectors. + * This function is aimed at interoperability with other libraries. + * \warning it returns the null pointer 0 in compressed mode */ + inline Index* innerNonZeroPtr() { return m_innerNonZeros; } + + /** \internal */ + inline Storage& data() { return m_data; } + /** \internal */ + inline const Storage& data() const { return m_data; } + + /** \returns the value of the matrix at position \a i, \a j + * This function returns Scalar(0) if the element is an explicit \em zero */ + inline Scalar coeff(Index row, Index col) const + { + eigen_assert(row>=0 && row=0 && col=0 && row=0 && col=start && "you probably called coeffRef on a non finalized matrix"); + if(end<=start) + return insert(row,col); + const Index p = m_data.searchLowerIndex(start,end-1,inner); + if((p=0 && row=0 && col::Constant(outerSize(), 2)); + } + return insertUncompressed(row,col); + } + + public: + + class InnerIterator; + class ReverseInnerIterator; + + /** Removes all non zeros but keep allocated memory */ + inline void setZero() + { + m_data.clear(); + memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(Index)); + if(m_innerNonZeros) + memset(m_innerNonZeros, 0, (m_outerSize)*sizeof(Index)); + } + + /** \returns the number of non zero coefficients */ + inline Index nonZeros() const + { + if(m_innerNonZeros) + return innerNonZeros().sum(); + return static_cast(m_data.size()); + } + + /** Preallocates \a reserveSize non zeros. + * + * Precondition: the matrix must be in compressed mode. */ + inline void reserve(Index reserveSize) + { + eigen_assert(isCompressed() && "This function does not make sense in non compressed mode."); + m_data.reserve(reserveSize); + } + + #ifdef EIGEN_PARSED_BY_DOXYGEN + /** Preallocates \a reserveSize[\c j] non zeros for each column (resp. row) \c j. + * + * This function turns the matrix in non-compressed mode */ + template + inline void reserve(const SizesType& reserveSizes); + #else + template + inline void reserve(const SizesType& reserveSizes, const typename SizesType::value_type& enableif = typename SizesType::value_type()) + { + EIGEN_UNUSED_VARIABLE(enableif); + reserveInnerVectors(reserveSizes); + } + template + inline void reserve(const SizesType& reserveSizes, const typename SizesType::Scalar& enableif = + #if (!defined(_MSC_VER)) || (_MSC_VER>=1500) // MSVC 2005 fails to compile with this typename + typename + #endif + SizesType::Scalar()) + { + EIGEN_UNUSED_VARIABLE(enableif); + reserveInnerVectors(reserveSizes); + } + #endif // EIGEN_PARSED_BY_DOXYGEN + protected: + template + inline void reserveInnerVectors(const SizesType& reserveSizes) + { + if(isCompressed()) + { + std::size_t totalReserveSize = 0; + // turn the matrix into non-compressed mode + m_innerNonZeros = static_cast(std::malloc(m_outerSize * sizeof(Index))); + if (!m_innerNonZeros) internal::throw_std_bad_alloc(); + + // temporarily use m_innerSizes to hold the new starting points. + Index* newOuterIndex = m_innerNonZeros; + + Index count = 0; + for(Index j=0; j=0; --j) + { + Index innerNNZ = previousOuterIndex - m_outerIndex[j]; + for(Index i=innerNNZ-1; i>=0; --i) + { + m_data.index(newOuterIndex[j]+i) = m_data.index(m_outerIndex[j]+i); + m_data.value(newOuterIndex[j]+i) = m_data.value(m_outerIndex[j]+i); + } + previousOuterIndex = m_outerIndex[j]; + m_outerIndex[j] = newOuterIndex[j]; + m_innerNonZeros[j] = innerNNZ; + } + m_outerIndex[m_outerSize] = m_outerIndex[m_outerSize-1] + m_innerNonZeros[m_outerSize-1] + reserveSizes[m_outerSize-1]; + + m_data.resize(m_outerIndex[m_outerSize]); + } + else + { + Index* newOuterIndex = static_cast(std::malloc((m_outerSize+1)*sizeof(Index))); + if (!newOuterIndex) internal::throw_std_bad_alloc(); + + Index count = 0; + for(Index j=0; j(reserveSizes[j], alreadyReserved); + count += toReserve + m_innerNonZeros[j]; + } + newOuterIndex[m_outerSize] = count; + + m_data.resize(count); + for(Index j=m_outerSize-1; j>=0; --j) + { + Index offset = newOuterIndex[j] - m_outerIndex[j]; + if(offset>0) + { + Index innerNNZ = m_innerNonZeros[j]; + for(Index i=innerNNZ-1; i>=0; --i) + { + m_data.index(newOuterIndex[j]+i) = m_data.index(m_outerIndex[j]+i); + m_data.value(newOuterIndex[j]+i) = m_data.value(m_outerIndex[j]+i); + } + } + } + + std::swap(m_outerIndex, newOuterIndex); + std::free(newOuterIndex); + } + + } + public: + + //--- low level purely coherent filling --- + + /** \internal + * \returns a reference to the non zero coefficient at position \a row, \a col assuming that: + * - the nonzero does not already exist + * - the new coefficient is the last one according to the storage order + * + * Before filling a given inner vector you must call the statVec(Index) function. + * + * After an insertion session, you should call the finalize() function. + * + * \sa insert, insertBackByOuterInner, startVec */ + inline Scalar& insertBack(Index row, Index col) + { + return insertBackByOuterInner(IsRowMajor?row:col, IsRowMajor?col:row); + } + + /** \internal + * \sa insertBack, startVec */ + inline Scalar& insertBackByOuterInner(Index outer, Index inner) + { + eigen_assert(size_t(m_outerIndex[outer+1]) == m_data.size() && "Invalid ordered insertion (invalid outer index)"); + eigen_assert( (m_outerIndex[outer+1]-m_outerIndex[outer]==0 || m_data.index(m_data.size()-1)(m_data.size()); + Index i = m_outerSize; + // find the last filled column + while (i>=0 && m_outerIndex[i]==0) + --i; + ++i; + while (i<=m_outerSize) + { + m_outerIndex[i] = size; + ++i; + } + } + } + + //--- + + template + void setFromTriplets(const InputIterators& begin, const InputIterators& end); + + void sumupDuplicates(); + + //--- + + /** \internal + * same as insert(Index,Index) except that the indices are given relative to the storage order */ + Scalar& insertByOuterInner(Index j, Index i) + { + return insert(IsRowMajor ? j : i, IsRowMajor ? i : j); + } + + /** Turns the matrix into the \em compressed format. + */ + void makeCompressed() + { + if(isCompressed()) + return; + + Index oldStart = m_outerIndex[1]; + m_outerIndex[1] = m_innerNonZeros[0]; + for(Index j=1; j0) + { + for(Index k=0; k(std::malloc(m_outerSize * sizeof(Index))); + for (Index i = 0; i < m_outerSize; i++) + { + m_innerNonZeros[i] = m_outerIndex[i+1] - m_outerIndex[i]; + } + } + + /** Suppresses all nonzeros which are \b much \b smaller \b than \a reference under the tolerence \a epsilon */ + void prune(const Scalar& reference, const RealScalar& epsilon = NumTraits::dummy_precision()) + { + prune(default_prunning_func(reference,epsilon)); + } + + /** Turns the matrix into compressed format, and suppresses all nonzeros which do not satisfy the predicate \a keep. + * The functor type \a KeepFunc must implement the following function: + * \code + * bool operator() (const Index& row, const Index& col, const Scalar& value) const; + * \endcode + * \sa prune(Scalar,RealScalar) + */ + template + void prune(const KeepFunc& keep = KeepFunc()) + { + // TODO optimize the uncompressed mode to avoid moving and allocating the data twice + // TODO also implement a unit test + makeCompressed(); + + Index k = 0; + for(Index j=0; jrows() == rows && this->cols() == cols) return; + + // If one dimension is null, then there is nothing to be preserved + if(rows==0 || cols==0) return resize(rows,cols); + + Index innerChange = IsRowMajor ? cols - this->cols() : rows - this->rows(); + Index outerChange = IsRowMajor ? rows - this->rows() : cols - this->cols(); + Index newInnerSize = IsRowMajor ? cols : rows; + + // Deals with inner non zeros + if (m_innerNonZeros) + { + // Resize m_innerNonZeros + Index *newInnerNonZeros = static_cast(std::realloc(m_innerNonZeros, (m_outerSize + outerChange) * sizeof(Index))); + if (!newInnerNonZeros) internal::throw_std_bad_alloc(); + m_innerNonZeros = newInnerNonZeros; + + for(Index i=m_outerSize; i(std::malloc((m_outerSize+outerChange+1) * sizeof(Index))); + if (!m_innerNonZeros) internal::throw_std_bad_alloc(); + for(Index i = 0; i < m_outerSize; i++) + m_innerNonZeros[i] = m_outerIndex[i+1] - m_outerIndex[i]; + } + + // Change the m_innerNonZeros in case of a decrease of inner size + if (m_innerNonZeros && innerChange < 0) + { + for(Index i = 0; i < m_outerSize + (std::min)(outerChange, Index(0)); i++) + { + Index &n = m_innerNonZeros[i]; + Index start = m_outerIndex[i]; + while (n > 0 && m_data.index(start+n-1) >= newInnerSize) --n; + } + } + + m_innerSize = newInnerSize; + + // Re-allocate outer index structure if necessary + if (outerChange == 0) + return; + + Index *newOuterIndex = static_cast(std::realloc(m_outerIndex, (m_outerSize + outerChange + 1) * sizeof(Index))); + if (!newOuterIndex) internal::throw_std_bad_alloc(); + m_outerIndex = newOuterIndex; + if (outerChange > 0) + { + Index last = m_outerSize == 0 ? 0 : m_outerIndex[m_outerSize]; + for(Index i=m_outerSize; i(std::malloc((outerSize + 1) * sizeof(Index))); + if (!m_outerIndex) internal::throw_std_bad_alloc(); + + m_outerSize = outerSize; + } + if(m_innerNonZeros) + { + std::free(m_innerNonZeros); + m_innerNonZeros = 0; + } + memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(Index)); + } + + /** \internal + * Resize the nonzero vector to \a size */ + void resizeNonZeros(Index size) + { + // TODO remove this function + m_data.resize(size); + } + + /** \returns a const expression of the diagonal coefficients */ + const Diagonal diagonal() const { return *this; } + + /** Default constructor yielding an empty \c 0 \c x \c 0 matrix */ + inline SparseMatrix() + : m_outerSize(-1), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0) + { + check_template_parameters(); + resize(0, 0); + } + + /** Constructs a \a rows \c x \a cols empty matrix */ + inline SparseMatrix(Index rows, Index cols) + : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0) + { + check_template_parameters(); + resize(rows, cols); + } + + /** Constructs a sparse matrix from the sparse expression \a other */ + template + inline SparseMatrix(const SparseMatrixBase& other) + : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0) + { + 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) + check_template_parameters(); + *this = other.derived(); + } + + /** Constructs a sparse matrix from the sparse selfadjoint view \a other */ + template + inline SparseMatrix(const SparseSelfAdjointView& other) + : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0) + { + check_template_parameters(); + *this = other; + } + + /** Copy constructor (it performs a deep copy) */ + inline SparseMatrix(const SparseMatrix& other) + : Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0) + { + check_template_parameters(); + *this = other.derived(); + } + + /** \brief Copy constructor with in-place evaluation */ + template + SparseMatrix(const ReturnByValue& other) + : Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0) + { + check_template_parameters(); + initAssignment(other); + other.evalTo(*this); + } + + /** Swaps the content of two sparse matrices of the same type. + * This is a fast operation that simply swaps the underlying pointers and parameters. */ + inline void swap(SparseMatrix& other) + { + //EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: swap\n"); + std::swap(m_outerIndex, other.m_outerIndex); + std::swap(m_innerSize, other.m_innerSize); + std::swap(m_outerSize, other.m_outerSize); + std::swap(m_innerNonZeros, other.m_innerNonZeros); + m_data.swap(other.m_data); + } + + /** Sets *this to the identity matrix. + * This function also turns the matrix into compressed mode, and drop any reserved memory. */ + inline void setIdentity() + { + eigen_assert(rows() == cols() && "ONLY FOR SQUARED MATRICES"); + this->m_data.resize(rows()); + Eigen::Map >(this->m_data.indexPtr(), rows()).setLinSpaced(0, rows()-1); + Eigen::Map >(this->m_data.valuePtr(), rows()).setOnes(); + Eigen::Map >(this->m_outerIndex, rows()+1).setLinSpaced(0, rows()); + std::free(m_innerNonZeros); + m_innerNonZeros = 0; + } + inline SparseMatrix& operator=(const SparseMatrix& other) + { + if (other.isRValue()) + { + swap(other.const_cast_derived()); + } + else if(this!=&other) + { + initAssignment(other); + if(other.isCompressed()) + { + memcpy(m_outerIndex, other.m_outerIndex, (m_outerSize+1)*sizeof(Index)); + m_data = other.m_data; + } + else + { + Base::operator=(other); + } + } + return *this; + } + + #ifndef EIGEN_PARSED_BY_DOXYGEN + template + inline SparseMatrix& operator=(const SparseSparseProduct& product) + { return Base::operator=(product); } + + template + inline SparseMatrix& operator=(const ReturnByValue& other) + { + initAssignment(other); + return Base::operator=(other.derived()); + } + + template + inline SparseMatrix& operator=(const EigenBase& other) + { return Base::operator=(other.derived()); } + #endif + + template + EIGEN_DONT_INLINE SparseMatrix& operator=(const SparseMatrixBase& other); + + friend std::ostream & operator << (std::ostream & s, const SparseMatrix& m) + { + EIGEN_DBG_SPARSE( + s << "Nonzero entries:\n"; + if(m.isCompressed()) + for (Index i=0; i&>(m); + return s; + } + + /** Destructor */ + inline ~SparseMatrix() + { + std::free(m_outerIndex); + std::free(m_innerNonZeros); + } + +#ifndef EIGEN_PARSED_BY_DOXYGEN + /** Overloaded for performance */ + Scalar sum() const; +#endif + +# ifdef EIGEN_SPARSEMATRIX_PLUGIN +# include EIGEN_SPARSEMATRIX_PLUGIN +# endif + +protected: + + template + void initAssignment(const Other& other) + { + resize(other.rows(), other.cols()); + if(m_innerNonZeros) + { + std::free(m_innerNonZeros); + m_innerNonZeros = 0; + } + } + + /** \internal + * \sa insert(Index,Index) */ + EIGEN_DONT_INLINE Scalar& insertCompressed(Index row, Index col); + + /** \internal + * A vector object that is equal to 0 everywhere but v at the position i */ + class SingletonVector + { + Index m_index; + Index m_value; + public: + typedef Index value_type; + SingletonVector(Index i, Index v) + : m_index(i), m_value(v) + {} + + Index operator[](Index i) const { return i==m_index ? m_value : 0; } + }; + + /** \internal + * \sa insert(Index,Index) */ + EIGEN_DONT_INLINE Scalar& insertUncompressed(Index row, Index col); + +public: + /** \internal + * \sa insert(Index,Index) */ + EIGEN_STRONG_INLINE Scalar& insertBackUncompressed(Index row, Index col) + { + const Index outer = IsRowMajor ? row : col; + const Index inner = IsRowMajor ? col : row; + + eigen_assert(!isCompressed()); + eigen_assert(m_innerNonZeros[outer]<=(m_outerIndex[outer+1] - m_outerIndex[outer])); + + Index p = m_outerIndex[outer] + m_innerNonZeros[outer]++; + m_data.index(p) = inner; + return (m_data.value(p) = 0); + } + +private: + static void check_template_parameters() + { + EIGEN_STATIC_ASSERT(NumTraits::IsSigned,THE_INDEX_TYPE_MUST_BE_A_SIGNED_TYPE); + EIGEN_STATIC_ASSERT((Options&(ColMajor|RowMajor))==Options,INVALID_MATRIX_TEMPLATE_PARAMETERS); + } + + struct default_prunning_func { + default_prunning_func(const Scalar& ref, const RealScalar& eps) : reference(ref), epsilon(eps) {} + inline bool operator() (const Index&, const Index&, const Scalar& value) const + { + return !internal::isMuchSmallerThan(value, reference, epsilon); + } + Scalar reference; + RealScalar epsilon; + }; +}; + +template +class SparseMatrix::InnerIterator +{ + public: + InnerIterator(const SparseMatrix& mat, Index outer) + : m_values(mat.valuePtr()), m_indices(mat.innerIndexPtr()), m_outer(outer), m_id(mat.m_outerIndex[outer]) + { + if(mat.isCompressed()) + m_end = mat.m_outerIndex[outer+1]; + else + m_end = m_id + mat.m_innerNonZeros[outer]; + } + + inline InnerIterator& operator++() { m_id++; return *this; } + + inline const Scalar& value() const { return m_values[m_id]; } + inline Scalar& valueRef() { return const_cast(m_values[m_id]); } + + inline Index index() const { return m_indices[m_id]; } + inline Index outer() const { return m_outer; } + inline Index row() const { return IsRowMajor ? m_outer : index(); } + inline Index col() const { return IsRowMajor ? index() : m_outer; } + + inline operator bool() const { return (m_id < m_end); } + + protected: + const Scalar* m_values; + const Index* m_indices; + const Index m_outer; + Index m_id; + Index m_end; +}; + +template +class SparseMatrix::ReverseInnerIterator +{ + public: + ReverseInnerIterator(const SparseMatrix& mat, Index outer) + : m_values(mat.valuePtr()), m_indices(mat.innerIndexPtr()), m_outer(outer), m_start(mat.m_outerIndex[outer]) + { + if(mat.isCompressed()) + m_id = mat.m_outerIndex[outer+1]; + else + m_id = m_start + mat.m_innerNonZeros[outer]; + } + + inline ReverseInnerIterator& operator--() { --m_id; return *this; } + + inline const Scalar& value() const { return m_values[m_id-1]; } + inline Scalar& valueRef() { return const_cast(m_values[m_id-1]); } + + inline Index index() const { return m_indices[m_id-1]; } + inline Index outer() const { return m_outer; } + inline Index row() const { return IsRowMajor ? m_outer : index(); } + inline Index col() const { return IsRowMajor ? index() : m_outer; } + + inline operator bool() const { return (m_id > m_start); } + + protected: + const Scalar* m_values; + const Index* m_indices; + const Index m_outer; + Index m_id; + const Index m_start; +}; + +namespace internal { + +template +void set_from_triplets(const InputIterator& begin, const InputIterator& end, SparseMatrixType& mat, int Options = 0) +{ + EIGEN_UNUSED_VARIABLE(Options); + enum { IsRowMajor = SparseMatrixType::IsRowMajor }; + typedef typename SparseMatrixType::Scalar Scalar; + typedef typename SparseMatrixType::Index Index; + SparseMatrix trMat(mat.rows(),mat.cols()); + + if(begin!=end) + { + // pass 1: count the nnz per inner-vector + Matrix wi(trMat.outerSize()); + wi.setZero(); + for(InputIterator it(begin); it!=end; ++it) + { + eigen_assert(it->row()>=0 && it->row()col()>=0 && it->col()col() : it->row())++; + } + + // pass 2: insert all the elements into trMat + trMat.reserve(wi); + for(InputIterator it(begin); it!=end; ++it) + trMat.insertBackUncompressed(it->row(),it->col()) = it->value(); + + // pass 3: + trMat.sumupDuplicates(); + } + + // pass 4: transposed copy -> implicit sorting + mat = trMat; +} + +} + + +/** Fill the matrix \c *this with the list of \em triplets defined by the iterator range \a begin - \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. + * In any case, the result is a \b sorted and \b compressed sparse matrix where the duplicates have been summed up. + * This is a \em O(n) operation, with \em n the number of triplet elements. + * The initial contents of \c *this is destroyed. + * The matrix \c *this must be properly resized beforehand using the SparseMatrix(Index,Index) constructor, + * or the resize(Index,Index) method. The sizes are not extracted from the triplet list. + * + * The \a InputIterators value_type must provide the following interface: + * \code + * Scalar value() const; // the value + * Scalar row() const; // the row index i + * Scalar col() const; // the column index j + * \endcode + * See for instance the Eigen::Triplet template class. + * + * Here is a typical usage example: + * \code + typedef Triplet T; + std::vector tripletList; + triplets.reserve(estimation_of_entries); + for(...) + { + // ... + tripletList.push_back(T(i,j,v_ij)); + } + SparseMatrixType m(rows,cols); + m.setFromTriplets(tripletList.begin(), tripletList.end()); + // m is ready to go! + * \endcode + * + * \warning The list of triplets is read multiple times (at least twice). Therefore, it is not recommended to define + * an abstract iterator over a complex data-structure that would be expensive to evaluate. The triplets should rather + * be explicitely stored into a std::vector for instance. + */ +template +template +void SparseMatrix::setFromTriplets(const InputIterators& begin, const InputIterators& end) +{ + internal::set_from_triplets(begin, end, *this); +} + +/** \internal */ +template +void SparseMatrix::sumupDuplicates() +{ + eigen_assert(!isCompressed()); + // TODO, in practice we should be able to use m_innerNonZeros for that task + Matrix wi(innerSize()); + wi.fill(-1); + Index count = 0; + // for each inner-vector, wi[inner_index] will hold the position of first element into the index/value buffers + for(Index j=0; j=start) + { + // we already meet this entry => accumulate it + m_data.value(wi(i)) += m_data.value(k); + } + else + { + m_data.value(count) = m_data.value(k); + m_data.index(count) = m_data.index(k); + wi(i) = count; + ++count; + } + } + m_outerIndex[j] = start; + } + m_outerIndex[m_outerSize] = count; + + // turn the matrix into compressed form + std::free(m_innerNonZeros); + m_innerNonZeros = 0; + m_data.resize(m_outerIndex[m_outerSize]); +} + +template +template +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) + + const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); + if (needToTranspose) + { + // two passes algorithm: + // 1 - compute the number of coeffs per dest inner vector + // 2 - do the actual copy/eval + // Since each coeff of the rhs has to be evaluated twice, let's evaluate it if needed + typedef typename internal::nested::type OtherCopy; + typedef typename internal::remove_all::type _OtherCopy; + OtherCopy otherCopy(other.derived()); + + SparseMatrix dest(other.rows(),other.cols()); + Eigen::Map > (dest.m_outerIndex,dest.outerSize()).setZero(); + + // pass 1 + // FIXME the above copy could be merged with that pass + for (Index j=0; j positions(dest.outerSize()); + for (Index j=0; jswap(dest); + return *this; + } + else + { + if(other.isRValue()) + initAssignment(other.derived()); + // there is no special optimization + return Base::operator=(other.derived()); + } +} + +template +EIGEN_DONT_INLINE typename SparseMatrix<_Scalar,_Options,_Index>::Scalar& SparseMatrix<_Scalar,_Options,_Index>::insertUncompressed(Index row, Index col) +{ + eigen_assert(!isCompressed()); + + const Index outer = IsRowMajor ? row : col; + const Index inner = IsRowMajor ? col : row; + + Index room = m_outerIndex[outer+1] - m_outerIndex[outer]; + Index innerNNZ = m_innerNonZeros[outer]; + if(innerNNZ>=room) + { + // this inner vector is full, we need to reallocate the whole buffer :( + reserve(SingletonVector(outer,std::max(2,innerNNZ))); + } + + Index startId = m_outerIndex[outer]; + Index p = startId + m_innerNonZeros[outer]; + while ( (p > startId) && (m_data.index(p-1) > inner) ) + { + m_data.index(p) = m_data.index(p-1); + m_data.value(p) = m_data.value(p-1); + --p; + } + eigen_assert((p<=startId || m_data.index(p-1)!=inner) && "you cannot insert an element that already exist, you must call coeffRef to this end"); + + m_innerNonZeros[outer]++; + + m_data.index(p) = inner; + return (m_data.value(p) = 0); +} + +template +EIGEN_DONT_INLINE typename SparseMatrix<_Scalar,_Options,_Index>::Scalar& SparseMatrix<_Scalar,_Options,_Index>::insertCompressed(Index row, Index col) +{ + eigen_assert(isCompressed()); + + const Index outer = IsRowMajor ? row : col; + const Index inner = IsRowMajor ? col : row; + + Index previousOuter = outer; + if (m_outerIndex[outer+1]==0) + { + // we start a new inner vector + while (previousOuter>=0 && m_outerIndex[previousOuter]==0) + { + m_outerIndex[previousOuter] = static_cast(m_data.size()); + --previousOuter; + } + m_outerIndex[outer+1] = m_outerIndex[outer]; + } + + // here we have to handle the tricky case where the outerIndex array + // starts with: [ 0 0 0 0 0 1 ...] and we are inserted in, e.g., + // the 2nd inner vector... + bool isLastVec = (!(previousOuter==-1 && m_data.size()!=0)) + && (size_t(m_outerIndex[outer+1]) == m_data.size()); + + size_t startId = m_outerIndex[outer]; + // FIXME let's make sure sizeof(long int) == sizeof(size_t) + size_t p = m_outerIndex[outer+1]; + ++m_outerIndex[outer+1]; + + double reallocRatio = 1; + if (m_data.allocatedSize()<=m_data.size()) + { + // if there is no preallocated memory, let's reserve a minimum of 32 elements + if (m_data.size()==0) + { + m_data.reserve(32); + } + else + { + // we need to reallocate the data, to reduce multiple reallocations + // we use a smart resize algorithm based on the current filling ratio + // in addition, we use double to avoid integers overflows + double nnzEstimate = double(m_outerIndex[outer])*double(m_outerSize)/double(outer+1); + reallocRatio = (nnzEstimate-double(m_data.size()))/double(m_data.size()); + // furthermore we bound the realloc ratio to: + // 1) reduce multiple minor realloc when the matrix is almost filled + // 2) avoid to allocate too much memory when the matrix is almost empty + reallocRatio = (std::min)((std::max)(reallocRatio,1.5),8.); + } + } + m_data.resize(m_data.size()+1,reallocRatio); + + if (!isLastVec) + { + if (previousOuter==-1) + { + // oops wrong guess. + // let's correct the outer offsets + for (Index k=0; k<=(outer+1); ++k) + m_outerIndex[k] = 0; + Index k=outer+1; + while(m_outerIndex[k]==0) + m_outerIndex[k++] = 1; + while (k<=m_outerSize && m_outerIndex[k]!=0) + m_outerIndex[k++]++; + p = 0; + --k; + k = m_outerIndex[k]-1; + while (k>0) + { + m_data.index(k) = m_data.index(k-1); + m_data.value(k) = m_data.value(k-1); + k--; + } + } + else + { + // we are not inserting into the last inner vec + // update outer indices: + Index j = outer+2; + while (j<=m_outerSize && m_outerIndex[j]!=0) + m_outerIndex[j++]++; + --j; + // shift data of last vecs: + Index k = m_outerIndex[j]-1; + while (k>=Index(p)) + { + m_data.index(k) = m_data.index(k-1); + m_data.value(k) = m_data.value(k-1); + k--; + } + } + } + + while ( (p > startId) && (m_data.index(p-1) > inner) ) + { + m_data.index(p) = m_data.index(p-1); + m_data.value(p) = m_data.value(p-1); + --p; + } + + m_data.index(p) = inner; + return (m_data.value(p) = 0); +} + +} // end namespace Eigen + +#endif // EIGEN_SPARSEMATRIX_H diff --git a/eigen/Eigen/src/SparseCore/SparseMatrixBase.h b/eigen/Eigen/src/SparseCore/SparseMatrixBase.h new file mode 100644 index 0000000..6f4a47c --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseMatrixBase.h @@ -0,0 +1,462 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2011 Gael Guennebaud +// +// 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_SPARSEMATRIXBASE_H +#define EIGEN_SPARSEMATRIXBASE_H + +namespace Eigen { + +/** \ingroup SparseCore_Module + * + * \class SparseMatrixBase + * + * \brief Base class of any sparse matrices or sparse expressions + * + * \tparam Derived + * + * This class can be extended with the help of the plugin mechanism described on the page + * \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_SPARSEMATRIXBASE_PLUGIN. + */ +template class SparseMatrixBase +#ifndef EIGEN_PARSED_BY_DOXYGEN + : public internal::special_scalar_op_base::Scalar, + typename NumTraits::Scalar>::Real, + EigenBase > +#else + : public EigenBase +#endif // not EIGEN_PARSED_BY_DOXYGEN +{ + public: + + typedef typename internal::traits::Scalar Scalar; + typedef typename internal::packet_traits::type PacketScalar; + typedef typename internal::traits::StorageKind StorageKind; + typedef typename internal::traits::Index Index; + typedef typename internal::traits::Index StorageIndex; + typedef typename internal::add_const_on_value_type_if_arithmetic< + typename internal::packet_traits::type + >::type PacketReturnType; + + typedef SparseMatrixBase StorageBaseType; + + template + Derived& operator=(const EigenBase &other) + { + other.derived().evalTo(derived()); + return derived(); + } + + enum { + + RowsAtCompileTime = internal::traits::RowsAtCompileTime, + /**< The number of rows at compile-time. This is just a copy of the value provided + * by the \a Derived type. If a value is not known at compile-time, + * it is set to the \a Dynamic constant. + * \sa MatrixBase::rows(), MatrixBase::cols(), ColsAtCompileTime, SizeAtCompileTime */ + + ColsAtCompileTime = internal::traits::ColsAtCompileTime, + /**< The number of columns at compile-time. This is just a copy of the value provided + * by the \a Derived type. If a value is not known at compile-time, + * it is set to the \a Dynamic constant. + * \sa MatrixBase::rows(), MatrixBase::cols(), RowsAtCompileTime, SizeAtCompileTime */ + + + SizeAtCompileTime = (internal::size_at_compile_time::RowsAtCompileTime, + internal::traits::ColsAtCompileTime>::ret), + /**< This is equal to the number of coefficients, i.e. the number of + * rows times the number of columns, or to \a Dynamic if this is not + * known at compile-time. \sa RowsAtCompileTime, ColsAtCompileTime */ + + MaxRowsAtCompileTime = RowsAtCompileTime, + MaxColsAtCompileTime = ColsAtCompileTime, + + MaxSizeAtCompileTime = (internal::size_at_compile_time::ret), + + IsVectorAtCompileTime = RowsAtCompileTime == 1 || ColsAtCompileTime == 1, + /**< This is set to true if either the number of rows or the number of + * columns is known at compile-time to be equal to 1. Indeed, in that case, + * we are dealing with a column-vector (if there is only one column) or with + * a row-vector (if there is only one row). */ + + Flags = internal::traits::Flags, + /**< This stores expression \ref flags flags which may or may not be inherited by new expressions + * constructed from this one. See the \ref flags "list of flags". + */ + + CoeffReadCost = internal::traits::CoeffReadCost, + /**< This is a rough measure of how expensive it is to read one coefficient from + * this expression. + */ + + IsRowMajor = Flags&RowMajorBit ? 1 : 0, + + InnerSizeAtCompileTime = int(IsVectorAtCompileTime) ? int(SizeAtCompileTime) + : int(IsRowMajor) ? int(ColsAtCompileTime) : int(RowsAtCompileTime), + + #ifndef EIGEN_PARSED_BY_DOXYGEN + _HasDirectAccess = (int(Flags)&DirectAccessBit) ? 1 : 0 // workaround sunCC + #endif + }; + + /** \internal the return type of MatrixBase::adjoint() */ + typedef typename internal::conditional::IsComplex, + CwiseUnaryOp, Eigen::Transpose >, + Transpose + >::type AdjointReturnType; + + + typedef SparseMatrix PlainObject; + + +#ifndef EIGEN_PARSED_BY_DOXYGEN + /** This is the "real scalar" type; if the \a Scalar type is already real numbers + * (e.g. int, float or double) then \a RealScalar is just the same as \a Scalar. If + * \a Scalar is \a std::complex then RealScalar is \a T. + * + * \sa class NumTraits + */ + typedef typename NumTraits::Real RealScalar; + + /** \internal the return type of coeff() + */ + typedef typename internal::conditional<_HasDirectAccess, const Scalar&, Scalar>::type CoeffReturnType; + + /** \internal Represents a matrix with all coefficients equal to one another*/ + typedef CwiseNullaryOp,Matrix > ConstantReturnType; + + /** type of the equivalent square matrix */ + typedef Matrix SquareMatrixType; + + inline const Derived& derived() const { return *static_cast(this); } + inline Derived& derived() { return *static_cast(this); } + inline Derived& const_cast_derived() const + { return *static_cast(const_cast(this)); } + + typedef internal::special_scalar_op_base > Base; + using Base::operator*; +#endif // not EIGEN_PARSED_BY_DOXYGEN + +#define EIGEN_CURRENT_STORAGE_BASE_CLASS Eigen::SparseMatrixBase +# include "../plugins/CommonCwiseUnaryOps.h" +# include "../plugins/CommonCwiseBinaryOps.h" +# include "../plugins/MatrixCwiseUnaryOps.h" +# include "../plugins/MatrixCwiseBinaryOps.h" +# include "../plugins/BlockMethods.h" +# ifdef EIGEN_SPARSEMATRIXBASE_PLUGIN +# include EIGEN_SPARSEMATRIXBASE_PLUGIN +# endif +# undef EIGEN_CURRENT_STORAGE_BASE_CLASS +#undef EIGEN_CURRENT_STORAGE_BASE_CLASS + + /** \returns the number of rows. \sa cols() */ + inline Index rows() const { return derived().rows(); } + /** \returns the number of columns. \sa rows() */ + inline Index cols() const { return derived().cols(); } + /** \returns the number of coefficients, which is \a rows()*cols(). + * \sa rows(), cols(). */ + inline Index size() const { return rows() * cols(); } + /** \returns the number of nonzero coefficients which is in practice the number + * of stored coefficients. */ + inline Index nonZeros() const { return derived().nonZeros(); } + /** \returns true if either the number of rows or the number of columns is equal to 1. + * In other words, this function returns + * \code rows()==1 || cols()==1 \endcode + * \sa rows(), cols(), IsVectorAtCompileTime. */ + inline bool isVector() const { return rows()==1 || cols()==1; } + /** \returns the size of the storage major dimension, + * i.e., the number of columns for a columns major matrix, and the number of rows otherwise */ + Index outerSize() const { return (int(Flags)&RowMajorBit) ? this->rows() : this->cols(); } + /** \returns the size of the inner dimension according to the storage order, + * i.e., the number of rows for a columns major matrix, and the number of cols otherwise */ + Index innerSize() const { return (int(Flags)&RowMajorBit) ? this->cols() : this->rows(); } + + bool isRValue() const { return m_isRValue; } + Derived& markAsRValue() { m_isRValue = true; return derived(); } + + SparseMatrixBase() : m_isRValue(false) { /* TODO check flags */ } + + + template + Derived& operator=(const ReturnByValue& other) + { + other.evalTo(derived()); + return derived(); + } + + + template + inline Derived& operator=(const SparseMatrixBase& other) + { + return assign(other.derived()); + } + + inline Derived& operator=(const Derived& other) + { +// if (other.isRValue()) +// derived().swap(other.const_cast_derived()); +// else + return assign(other.derived()); + } + + protected: + + template + inline Derived& assign(const OtherDerived& other) + { + const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); + const Index outerSize = (int(OtherDerived::Flags) & RowMajorBit) ? other.rows() : other.cols(); + if ((!transpose) && other.isRValue()) + { + // eval without temporary + derived().resize(other.rows(), other.cols()); + derived().setZero(); + derived().reserve((std::max)(this->rows(),this->cols())*2); + for (Index j=0; j + inline void assignGeneric(const OtherDerived& other) + { + //const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); + eigen_assert(( ((internal::traits::SupportedAccessPatterns&OuterRandomAccessPattern)==OuterRandomAccessPattern) || + (!((Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit)))) && + "the transpose operation is supposed to be handled in SparseMatrix::operator="); + + enum { Flip = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit) }; + + const Index outerSize = other.outerSize(); + //typedef typename internal::conditional, Derived>::type TempType; + // thanks to shallow copies, we always eval to a tempary + Derived temp(other.rows(), other.cols()); + + temp.reserve((std::max)(this->rows(),this->cols())*2); + for (Index j=0; j + inline Derived& operator=(const SparseSparseProduct& product); + + friend std::ostream & operator << (std::ostream & s, const SparseMatrixBase& m) + { + typedef typename Derived::Nested Nested; + typedef typename internal::remove_all::type NestedCleaned; + + if (Flags&RowMajorBit) + { + const Nested nm(m.derived()); + for (Index row=0; row trans = m; + s << static_cast >&>(trans); + } + } + return s; + } + + template + Derived& operator+=(const SparseMatrixBase& other); + template + Derived& operator-=(const SparseMatrixBase& other); + + Derived& operator*=(const Scalar& other); + Derived& operator/=(const Scalar& other); + + template struct CwiseProductDenseReturnType { + typedef CwiseBinaryOp::Scalar, + typename internal::traits::Scalar + >::ReturnType>, + const Derived, + const OtherDerived + > Type; + }; + + template + EIGEN_STRONG_INLINE const typename CwiseProductDenseReturnType::Type + cwiseProduct(const MatrixBase &other) const; + + // sparse * sparse + template + const typename SparseSparseProductReturnType::Type + operator*(const SparseMatrixBase &other) const; + + // sparse * diagonal + template + const SparseDiagonalProduct + operator*(const DiagonalBase &other) const; + + // diagonal * sparse + template friend + const SparseDiagonalProduct + operator*(const DiagonalBase &lhs, const SparseMatrixBase& rhs) + { return SparseDiagonalProduct(lhs.derived(), rhs.derived()); } + + /** dense * sparse (return a dense object unless it is an outer product) */ + template friend + const typename DenseSparseProductReturnType::Type + operator*(const MatrixBase& lhs, const Derived& rhs) + { return typename DenseSparseProductReturnType::Type(lhs.derived(),rhs); } + + /** sparse * dense (returns a dense object unless it is an outer product) */ + template + const typename SparseDenseProductReturnType::Type + operator*(const MatrixBase &other) const + { return typename SparseDenseProductReturnType::Type(derived(), other.derived()); } + + /** \returns an expression of P H P^-1 where H is the matrix represented by \c *this */ + SparseSymmetricPermutationProduct twistedBy(const PermutationMatrix& perm) const + { + return SparseSymmetricPermutationProduct(derived(), perm); + } + + template + Derived& operator*=(const SparseMatrixBase& other); + + #ifdef EIGEN2_SUPPORT + // deprecated + template + typename internal::plain_matrix_type_column_major::type + solveTriangular(const MatrixBase& other) const; + + // deprecated + template + void solveTriangularInPlace(MatrixBase& other) const; + #endif // EIGEN2_SUPPORT + + template + inline const SparseTriangularView triangularView() const; + + template inline const SparseSelfAdjointView selfadjointView() const; + template inline SparseSelfAdjointView selfadjointView(); + + template Scalar dot(const MatrixBase& other) const; + template Scalar dot(const SparseMatrixBase& other) const; + RealScalar squaredNorm() const; + RealScalar norm() const; + RealScalar blueNorm() const; + + Transpose transpose() { return derived(); } + const Transpose transpose() const { return derived(); } + const AdjointReturnType adjoint() const { return transpose(); } + + // inner-vector + typedef Block InnerVectorReturnType; + typedef Block ConstInnerVectorReturnType; + InnerVectorReturnType innerVector(Index outer); + const ConstInnerVectorReturnType innerVector(Index outer) const; + + // set of inner-vectors + typedef Block InnerVectorsReturnType; + typedef Block ConstInnerVectorsReturnType; + InnerVectorsReturnType innerVectors(Index outerStart, Index outerSize); + const ConstInnerVectorsReturnType innerVectors(Index outerStart, Index outerSize) const; + + /** \internal use operator= */ + template + void evalTo(MatrixBase& dst) const + { + dst.setZero(); + for (Index j=0; j toDense() const + { + return derived(); + } + + template + bool isApprox(const SparseMatrixBase& other, + const RealScalar& prec = NumTraits::dummy_precision()) const + { return toDense().isApprox(other.toDense(),prec); } + + template + bool isApprox(const MatrixBase& other, + const RealScalar& prec = NumTraits::dummy_precision()) const + { return toDense().isApprox(other,prec); } + + /** \returns the matrix or vector obtained by evaluating this expression. + * + * Notice that in the case of a plain matrix or vector (not an expression) this function just returns + * a const reference, in order to avoid a useless copy. + */ + inline const typename internal::eval::type eval() const + { return typename internal::eval::type(derived()); } + + Scalar sum() const; + + protected: + + bool m_isRValue; +}; + +} // end namespace Eigen + +#endif // EIGEN_SPARSEMATRIXBASE_H diff --git a/eigen/Eigen/src/SparseCore/SparsePermutation.h b/eigen/Eigen/src/SparseCore/SparsePermutation.h new file mode 100644 index 0000000..75e2100 --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparsePermutation.h @@ -0,0 +1,148 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Gael Guennebaud +// +// 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_SPARSE_PERMUTATION_H +#define EIGEN_SPARSE_PERMUTATION_H + +// This file implements sparse * permutation products + +namespace Eigen { + +namespace internal { + +template +struct traits > +{ + typedef typename remove_all::type MatrixTypeNestedCleaned; + typedef typename MatrixTypeNestedCleaned::Scalar Scalar; + typedef typename MatrixTypeNestedCleaned::Index Index; + enum { + SrcStorageOrder = MatrixTypeNestedCleaned::Flags&RowMajorBit ? RowMajor : ColMajor, + MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight + }; + + typedef typename internal::conditional, + SparseMatrix >::type ReturnType; +}; + +template +struct permut_sparsematrix_product_retval + : public ReturnByValue > +{ + typedef typename remove_all::type MatrixTypeNestedCleaned; + typedef typename MatrixTypeNestedCleaned::Scalar Scalar; + typedef typename MatrixTypeNestedCleaned::Index Index; + + enum { + SrcStorageOrder = MatrixTypeNestedCleaned::Flags&RowMajorBit ? RowMajor : ColMajor, + MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight + }; + + permut_sparsematrix_product_retval(const PermutationType& perm, const MatrixType& matrix) + : m_permutation(perm), m_matrix(matrix) + {} + + inline int rows() const { return m_matrix.rows(); } + inline int cols() const { return m_matrix.cols(); } + + template inline void evalTo(Dest& dst) const + { + if(MoveOuter) + { + SparseMatrix tmp(m_matrix.rows(), m_matrix.cols()); + Matrix sizes(m_matrix.outerSize()); + for(Index j=0; j tmp(m_matrix.rows(), m_matrix.cols()); + Matrix sizes(tmp.outerSize()); + sizes.setZero(); + PermutationMatrix perm; + if((Side==OnTheLeft) ^ Transposed) + perm = m_permutation; + else + perm = m_permutation.transpose(); + + for(Index j=0; j +inline const internal::permut_sparsematrix_product_retval, SparseDerived, OnTheRight, false> +operator*(const SparseMatrixBase& matrix, const PermutationBase& perm) +{ + return internal::permut_sparsematrix_product_retval, SparseDerived, OnTheRight, false>(perm, matrix.derived()); +} + +/** \returns the matrix with the permutation applied to the rows + */ +template +inline const internal::permut_sparsematrix_product_retval, SparseDerived, OnTheLeft, false> +operator*( const PermutationBase& perm, const SparseMatrixBase& matrix) +{ + return internal::permut_sparsematrix_product_retval, SparseDerived, OnTheLeft, false>(perm, matrix.derived()); +} + + + +/** \returns the matrix with the inverse permutation applied to the columns. + */ +template +inline const internal::permut_sparsematrix_product_retval, SparseDerived, OnTheRight, true> +operator*(const SparseMatrixBase& matrix, const Transpose >& tperm) +{ + return internal::permut_sparsematrix_product_retval, SparseDerived, OnTheRight, true>(tperm.nestedPermutation(), matrix.derived()); +} + +/** \returns the matrix with the inverse permutation applied to the rows. + */ +template +inline const internal::permut_sparsematrix_product_retval, SparseDerived, OnTheLeft, true> +operator*(const Transpose >& tperm, const SparseMatrixBase& matrix) +{ + return internal::permut_sparsematrix_product_retval, SparseDerived, OnTheLeft, true>(tperm.nestedPermutation(), matrix.derived()); +} + +} // end namespace Eigen + +#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H diff --git a/eigen/Eigen/src/SparseCore/SparseProduct.h b/eigen/Eigen/src/SparseCore/SparseProduct.h new file mode 100644 index 0000000..cf76630 --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseProduct.h @@ -0,0 +1,188 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2010 Gael Guennebaud +// +// 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_SPARSEPRODUCT_H +#define EIGEN_SPARSEPRODUCT_H + +namespace Eigen { + +template +struct SparseSparseProductReturnType +{ + typedef typename internal::traits::Scalar Scalar; + typedef typename internal::traits::Index Index; + enum { + LhsRowMajor = internal::traits::Flags & RowMajorBit, + RhsRowMajor = internal::traits::Flags & RowMajorBit, + TransposeRhs = (!LhsRowMajor) && RhsRowMajor, + TransposeLhs = LhsRowMajor && (!RhsRowMajor) + }; + + typedef typename internal::conditional, + typename internal::nested::type>::type LhsNested; + + typedef typename internal::conditional, + typename internal::nested::type>::type RhsNested; + + typedef SparseSparseProduct Type; +}; + +namespace internal { +template +struct traits > +{ + typedef MatrixXpr XprKind; + // clean the nested types: + typedef typename remove_all::type _LhsNested; + typedef typename remove_all::type _RhsNested; + typedef typename _LhsNested::Scalar Scalar; + typedef typename promote_index_type::Index, + typename traits<_RhsNested>::Index>::type Index; + + enum { + LhsCoeffReadCost = _LhsNested::CoeffReadCost, + RhsCoeffReadCost = _RhsNested::CoeffReadCost, + LhsFlags = _LhsNested::Flags, + RhsFlags = _RhsNested::Flags, + + RowsAtCompileTime = _LhsNested::RowsAtCompileTime, + ColsAtCompileTime = _RhsNested::ColsAtCompileTime, + MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime, + MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime, + + InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime), + + EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit), + + RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit), + + Flags = (int(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits) + | EvalBeforeAssigningBit + | EvalBeforeNestingBit, + + CoeffReadCost = Dynamic + }; + + typedef Sparse StorageKind; +}; + +} // end namespace internal + +template +class SparseSparseProduct : internal::no_assignment_operator, + public SparseMatrixBase > +{ + public: + + typedef SparseMatrixBase Base; + EIGEN_DENSE_PUBLIC_INTERFACE(SparseSparseProduct) + + private: + + typedef typename internal::traits::_LhsNested _LhsNested; + typedef typename internal::traits::_RhsNested _RhsNested; + + public: + + template + EIGEN_STRONG_INLINE SparseSparseProduct(const Lhs& lhs, const Rhs& rhs) + : m_lhs(lhs), m_rhs(rhs), m_tolerance(0), m_conservative(true) + { + init(); + } + + template + EIGEN_STRONG_INLINE SparseSparseProduct(const Lhs& lhs, const Rhs& rhs, const RealScalar& tolerance) + : m_lhs(lhs), m_rhs(rhs), m_tolerance(tolerance), m_conservative(false) + { + init(); + } + + SparseSparseProduct pruned(const Scalar& reference = 0, const RealScalar& epsilon = NumTraits::dummy_precision()) const + { + using std::abs; + return SparseSparseProduct(m_lhs,m_rhs,abs(reference)*epsilon); + } + + template + void evalTo(Dest& result) const + { + if(m_conservative) + internal::conservative_sparse_sparse_product_selector<_LhsNested, _RhsNested, Dest>::run(lhs(),rhs(),result); + else + internal::sparse_sparse_product_with_pruning_selector<_LhsNested, _RhsNested, Dest>::run(lhs(),rhs(),result,m_tolerance); + } + + EIGEN_STRONG_INLINE Index rows() const { return m_lhs.rows(); } + EIGEN_STRONG_INLINE Index cols() const { return m_rhs.cols(); } + + EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; } + EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; } + + protected: + void init() + { + eigen_assert(m_lhs.cols() == m_rhs.rows()); + + enum { + ProductIsValid = _LhsNested::ColsAtCompileTime==Dynamic + || _RhsNested::RowsAtCompileTime==Dynamic + || int(_LhsNested::ColsAtCompileTime)==int(_RhsNested::RowsAtCompileTime), + AreVectors = _LhsNested::IsVectorAtCompileTime && _RhsNested::IsVectorAtCompileTime, + SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(_LhsNested,_RhsNested) + }; + // note to the lost user: + // * for a dot product use: v1.dot(v2) + // * for a coeff-wise product use: v1.cwise()*v2 + EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes), + INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS) + EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors), + INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION) + EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT) + } + + LhsNested m_lhs; + RhsNested m_rhs; + RealScalar m_tolerance; + bool m_conservative; +}; + +// sparse = sparse * sparse +template +template +inline Derived& SparseMatrixBase::operator=(const SparseSparseProduct& product) +{ + product.evalTo(derived()); + return derived(); +} + +/** \returns an expression of the product of two sparse matrices. + * By default a conservative product preserving the symbolic non zeros is performed. + * The automatic pruning of the small values can be achieved by calling the pruned() function + * in which case a totally different product algorithm is employed: + * \code + * C = (A*B).pruned(); // supress numerical zeros (exact) + * C = (A*B).pruned(ref); + * C = (A*B).pruned(ref,epsilon); + * \endcode + * where \c ref is a meaningful non zero reference value. + * */ +template +template +inline const typename SparseSparseProductReturnType::Type +SparseMatrixBase::operator*(const SparseMatrixBase &other) const +{ + return typename SparseSparseProductReturnType::Type(derived(), other.derived()); +} + +} // end namespace Eigen + +#endif // EIGEN_SPARSEPRODUCT_H diff --git a/eigen/Eigen/src/SparseCore/SparseRedux.h b/eigen/Eigen/src/SparseCore/SparseRedux.h new file mode 100644 index 0000000..51ed9ae --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseRedux.h @@ -0,0 +1,48 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// +// 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_SPARSEREDUX_H +#define EIGEN_SPARSEREDUX_H + +namespace Eigen { + +template +typename internal::traits::Scalar +SparseMatrixBase::sum() const +{ + eigen_assert(rows()>0 && cols()>0 && "you are using a non initialized matrix"); + Scalar res(0); + for (Index j=0; j +typename internal::traits >::Scalar +SparseMatrix<_Scalar,_Options,_Index>::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(); + else + return Base::sum(); +} + +template +typename internal::traits >::Scalar +SparseVector<_Scalar,_Options,_Index>::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(); +} + +} // end namespace Eigen + +#endif // EIGEN_SPARSEREDUX_H diff --git a/eigen/Eigen/src/SparseCore/SparseSelfAdjointView.h b/eigen/Eigen/src/SparseCore/SparseSelfAdjointView.h new file mode 100644 index 0000000..0eda96b --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseSelfAdjointView.h @@ -0,0 +1,507 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud +// +// 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_SPARSE_SELFADJOINTVIEW_H +#define EIGEN_SPARSE_SELFADJOINTVIEW_H + +namespace Eigen { + +/** \ingroup SparseCore_Module + * \class SparseSelfAdjointView + * + * \brief Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix. + * + * \param MatrixType the type of the dense matrix storing the coefficients + * \param UpLo can be either \c #Lower or \c #Upper + * + * This class is an expression of a sefladjoint matrix from a triangular part of a matrix + * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView() + * and most of the time this is the only way that it is used. + * + * \sa SparseMatrixBase::selfadjointView() + */ +template +class SparseSelfAdjointTimeDenseProduct; + +template +class DenseTimeSparseSelfAdjointProduct; + +namespace internal { + +template +struct traits > : traits { +}; + +template +void permute_symm_to_symm(const MatrixType& mat, SparseMatrix& _dest, const typename MatrixType::Index* perm = 0); + +template +void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix& _dest, const typename MatrixType::Index* perm = 0); + +} + +template class SparseSelfAdjointView + : public EigenBase > +{ + public: + + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::Index Index; + typedef Matrix VectorI; + typedef typename MatrixType::Nested MatrixTypeNested; + typedef typename internal::remove_all::type _MatrixTypeNested; + + inline SparseSelfAdjointView(const MatrixType& matrix) : m_matrix(matrix) + { + eigen_assert(rows()==cols() && "SelfAdjointView is only for squared matrices"); + } + + inline Index rows() const { return m_matrix.rows(); } + inline Index cols() const { return m_matrix.cols(); } + + /** \internal \returns a reference to the nested matrix */ + const _MatrixTypeNested& matrix() const { return m_matrix; } + _MatrixTypeNested& matrix() { return m_matrix.const_cast_derived(); } + + /** \returns an expression of the matrix product between a sparse self-adjoint matrix \c *this and a sparse matrix \a rhs. + * + * Note that there is no algorithmic advantage of performing such a product compared to a general sparse-sparse matrix product. + * Indeed, the SparseSelfadjointView operand is first copied into a temporary SparseMatrix before computing the product. + */ + template + SparseSparseProduct + operator*(const SparseMatrixBase& rhs) const + { + return SparseSparseProduct(*this, rhs.derived()); + } + + /** \returns an expression of the matrix product between a sparse matrix \a lhs and a sparse self-adjoint matrix \a rhs. + * + * Note that there is no algorithmic advantage of performing such a product compared to a general sparse-sparse matrix product. + * Indeed, the SparseSelfadjointView operand is first copied into a temporary SparseMatrix before computing the product. + */ + template friend + SparseSparseProduct + operator*(const SparseMatrixBase& lhs, const SparseSelfAdjointView& rhs) + { + return SparseSparseProduct(lhs.derived(), rhs); + } + + /** Efficient sparse self-adjoint matrix times dense vector/matrix product */ + template + SparseSelfAdjointTimeDenseProduct + operator*(const MatrixBase& rhs) const + { + return SparseSelfAdjointTimeDenseProduct(m_matrix, rhs.derived()); + } + + /** Efficient dense vector/matrix times sparse self-adjoint matrix product */ + template friend + DenseTimeSparseSelfAdjointProduct + operator*(const MatrixBase& lhs, const SparseSelfAdjointView& rhs) + { + return DenseTimeSparseSelfAdjointProduct(lhs.derived(), rhs.m_matrix); + } + + /** Perform a symmetric rank K update of the selfadjoint matrix \c *this: + * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix. + * + * \returns a reference to \c *this + * + * To perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply + * call this function with u.adjoint(). + */ + template + SparseSelfAdjointView& rankUpdate(const SparseMatrixBase& u, const Scalar& alpha = Scalar(1)); + + /** \internal triggered by sparse_matrix = SparseSelfadjointView; */ + template void evalTo(SparseMatrix& _dest) const + { + internal::permute_symm_to_fullsymm(m_matrix, _dest); + } + + template void evalTo(DynamicSparseMatrix& _dest) const + { + // TODO directly evaluate into _dest; + SparseMatrix tmp(_dest.rows(),_dest.cols()); + internal::permute_symm_to_fullsymm(m_matrix, tmp); + _dest = tmp; + } + + /** \returns an expression of P H P^-1 */ + SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo> twistedBy(const PermutationMatrix& perm) const + { + return SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo>(m_matrix, perm); + } + + template + SparseSelfAdjointView& operator=(const SparseSymmetricPermutationProduct& permutedMatrix) + { + permutedMatrix.evalTo(*this); + return *this; + } + + + SparseSelfAdjointView& operator=(const SparseSelfAdjointView& src) + { + PermutationMatrix pnull; + return *this = src.twistedBy(pnull); + } + + template + SparseSelfAdjointView& operator=(const SparseSelfAdjointView& src) + { + PermutationMatrix pnull; + return *this = src.twistedBy(pnull); + } + + + // const SparseLLT llt() const; + // const SparseLDLT ldlt() const; + + protected: + + typename MatrixType::Nested m_matrix; + mutable VectorI m_countPerRow; + mutable VectorI m_countPerCol; +}; + +/*************************************************************************** +* Implementation of SparseMatrixBase methods +***************************************************************************/ + +template +template +const SparseSelfAdjointView SparseMatrixBase::selfadjointView() const +{ + return derived(); +} + +template +template +SparseSelfAdjointView SparseMatrixBase::selfadjointView() +{ + return derived(); +} + +/*************************************************************************** +* Implementation of SparseSelfAdjointView methods +***************************************************************************/ + +template +template +SparseSelfAdjointView& +SparseSelfAdjointView::rankUpdate(const SparseMatrixBase& u, const Scalar& alpha) +{ + SparseMatrix tmp = u * u.adjoint(); + if(alpha==Scalar(0)) + m_matrix.const_cast_derived() = tmp.template triangularView(); + else + m_matrix.const_cast_derived() += alpha * tmp.template triangularView(); + + return *this; +} + +/*************************************************************************** +* Implementation of sparse self-adjoint time dense matrix +***************************************************************************/ + +namespace internal { +template +struct traits > + : traits, Lhs, Rhs> > +{ + typedef Dense StorageKind; +}; +} + +template +class SparseSelfAdjointTimeDenseProduct + : public ProductBase, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseSelfAdjointTimeDenseProduct) + + SparseSelfAdjointTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + {} + + template void scaleAndAddTo(Dest& dest, const Scalar& alpha) const + { + EIGEN_ONLY_USED_FOR_DEBUG(alpha); + // TODO use alpha + eigen_assert(alpha==Scalar(1) && "alpha != 1 is not implemented yet, sorry"); + typedef typename internal::remove_all::type _Lhs; + typedef typename _Lhs::InnerIterator LhsInnerIterator; + enum { + LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit, + ProcessFirstHalf = + ((UpLo&(Upper|Lower))==(Upper|Lower)) + || ( (UpLo&Upper) && !LhsIsRowMajor) + || ( (UpLo&Lower) && LhsIsRowMajor), + ProcessSecondHalf = !ProcessFirstHalf + }; + for (Index j=0; j +struct traits > + : traits, Lhs, Rhs> > +{}; +} + +template +class DenseTimeSparseSelfAdjointProduct + : public ProductBase, Lhs, Rhs> +{ + public: + EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseSelfAdjointProduct) + + DenseTimeSparseSelfAdjointProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + {} + + template void scaleAndAddTo(Dest& /*dest*/, const Scalar& /*alpha*/) const + { + // TODO + } + + private: + DenseTimeSparseSelfAdjointProduct& operator=(const DenseTimeSparseSelfAdjointProduct&); +}; + +/*************************************************************************** +* Implementation of symmetric copies and permutations +***************************************************************************/ +namespace internal { + +template +struct traits > : traits { +}; + +template +void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix& _dest, const typename MatrixType::Index* perm) +{ + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + typedef SparseMatrix Dest; + typedef Matrix VectorI; + + Dest& dest(_dest.derived()); + enum { + StorageOrderMatch = int(Dest::IsRowMajor) == int(MatrixType::IsRowMajor) + }; + + Index size = mat.rows(); + VectorI count; + count.resize(size); + count.setZero(); + dest.resize(size,size); + for(Index j = 0; jc) || ( UpLo==Upper && rc) || ( (UpLo&Upper)==Upper && r +void permute_symm_to_symm(const MatrixType& mat, SparseMatrix& _dest, const typename MatrixType::Index* perm) +{ + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + SparseMatrix& dest(_dest.derived()); + typedef Matrix VectorI; + enum { + SrcOrder = MatrixType::IsRowMajor ? RowMajor : ColMajor, + StorageOrderMatch = int(SrcOrder) == int(DstOrder), + DstUpLo = DstOrder==RowMajor ? (_DstUpLo==Upper ? Lower : Upper) : _DstUpLo, + SrcUpLo = SrcOrder==RowMajor ? (_SrcUpLo==Upper ? Lower : Upper) : _SrcUpLo + }; + + Index size = mat.rows(); + VectorI count(size); + count.setZero(); + dest.resize(size,size); + for(Index j = 0; jj)) + continue; + + Index ip = perm ? perm[i] : i; + count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; + } + } + dest.outerIndexPtr()[0] = 0; + for(Index j=0; jj)) + continue; + + Index jp = perm ? perm[j] : j; + Index ip = perm? perm[i] : i; + + Index k = count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; + dest.innerIndexPtr()[k] = int(DstUpLo)==int(Lower) ? (std::max)(ip,jp) : (std::min)(ip,jp); + + if(!StorageOrderMatch) std::swap(ip,jp); + if( ((int(DstUpLo)==int(Lower) && ipjp))) + dest.valuePtr()[k] = numext::conj(it.value()); + else + dest.valuePtr()[k] = it.value(); + } + } +} + +} + +template +class SparseSymmetricPermutationProduct + : public EigenBase > +{ + public: + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::Index Index; + protected: + typedef PermutationMatrix Perm; + public: + typedef Matrix VectorI; + typedef typename MatrixType::Nested MatrixTypeNested; + typedef typename internal::remove_all::type _MatrixTypeNested; + + SparseSymmetricPermutationProduct(const MatrixType& mat, const Perm& perm) + : m_matrix(mat), m_perm(perm) + {} + + inline Index rows() const { return m_matrix.rows(); } + inline Index cols() const { return m_matrix.cols(); } + + template + void evalTo(SparseMatrix& _dest) const + { +// internal::permute_symm_to_fullsymm(m_matrix,_dest,m_perm.indices().data()); + SparseMatrix tmp; + internal::permute_symm_to_fullsymm(m_matrix,tmp,m_perm.indices().data()); + _dest = tmp; + } + + template void evalTo(SparseSelfAdjointView& dest) const + { + internal::permute_symm_to_symm(m_matrix,dest.matrix(),m_perm.indices().data()); + } + + protected: + MatrixTypeNested m_matrix; + const Perm& m_perm; + +}; + +} // end namespace Eigen + +#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H diff --git a/eigen/Eigen/src/SparseCore/SparseSparseProductWithPruning.h b/eigen/Eigen/src/SparseCore/SparseSparseProductWithPruning.h new file mode 100644 index 0000000..55b84a4 --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseSparseProductWithPruning.h @@ -0,0 +1,150 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2011 Gael Guennebaud +// +// 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_SPARSESPARSEPRODUCTWITHPRUNING_H +#define EIGEN_SPARSESPARSEPRODUCTWITHPRUNING_H + +namespace Eigen { + +namespace internal { + + +// perform a pseudo in-place sparse * sparse product assuming all matrices are col major +template +static void sparse_sparse_product_with_pruning_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res, const typename ResultType::RealScalar& tolerance) +{ + // return sparse_sparse_product_with_pruning_impl2(lhs,rhs,res); + + typedef typename remove_all::type::Scalar Scalar; + typedef typename remove_all::type::Index Index; + + // make sure to call innerSize/outerSize since we fake the storage order. + Index rows = lhs.innerSize(); + Index cols = rhs.outerSize(); + //Index size = lhs.outerSize(); + eigen_assert(lhs.outerSize() == rhs.innerSize()); + + // allocate a temporary buffer + AmbiVector tempVector(rows); + + // estimate the number of non zero entries + // given a rhs column containing Y non zeros, we assume that the respective Y columns + // of the lhs differs in average of one non zeros, thus the number of non zeros for + // the product of a rhs column with the lhs is X+Y where X is the average number of non zero + // per column of the lhs. + // Therefore, we have nnz(lhs*rhs) = nnz(lhs) + nnz(rhs) + Index estimated_nnz_prod = lhs.nonZeros() + rhs.nonZeros(); + + // mimics a resizeByInnerOuter: + if(ResultType::IsRowMajor) + res.resize(cols, rows); + else + res.resize(rows, cols); + + res.reserve(estimated_nnz_prod); + double ratioColRes = double(estimated_nnz_prod)/(double(lhs.rows())*double(rhs.cols())); + for (Index j=0; j::Iterator it(tempVector,tolerance); it; ++it) + res.insertBackByOuterInner(j,it.index()) = it.value(); + } + res.finalize(); +} + +template::Flags&RowMajorBit, + int RhsStorageOrder = traits::Flags&RowMajorBit, + int ResStorageOrder = traits::Flags&RowMajorBit> +struct sparse_sparse_product_with_pruning_selector; + +template +struct sparse_sparse_product_with_pruning_selector +{ + typedef typename traits::type>::Scalar Scalar; + typedef typename ResultType::RealScalar RealScalar; + + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res, const RealScalar& tolerance) + { + typename remove_all::type _res(res.rows(), res.cols()); + internal::sparse_sparse_product_with_pruning_impl(lhs, rhs, _res, tolerance); + res.swap(_res); + } +}; + +template +struct sparse_sparse_product_with_pruning_selector +{ + typedef typename ResultType::RealScalar RealScalar; + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res, const RealScalar& tolerance) + { + // we need a col-major matrix to hold the result + typedef SparseMatrix SparseTemporaryType; + SparseTemporaryType _res(res.rows(), res.cols()); + internal::sparse_sparse_product_with_pruning_impl(lhs, rhs, _res, tolerance); + res = _res; + } +}; + +template +struct sparse_sparse_product_with_pruning_selector +{ + typedef typename ResultType::RealScalar RealScalar; + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res, const RealScalar& tolerance) + { + // let's transpose the product to get a column x column product + typename remove_all::type _res(res.rows(), res.cols()); + internal::sparse_sparse_product_with_pruning_impl(rhs, lhs, _res, tolerance); + res.swap(_res); + } +}; + +template +struct sparse_sparse_product_with_pruning_selector +{ + typedef typename ResultType::RealScalar RealScalar; + static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res, const RealScalar& tolerance) + { + typedef SparseMatrix ColMajorMatrixLhs; + typedef SparseMatrix ColMajorMatrixRhs; + ColMajorMatrixLhs colLhs(lhs); + ColMajorMatrixRhs colRhs(rhs); + internal::sparse_sparse_product_with_pruning_impl(colLhs, colRhs, res, tolerance); + + // let's transpose the product to get a column x column product +// typedef SparseMatrix SparseTemporaryType; +// SparseTemporaryType _res(res.cols(), res.rows()); +// sparse_sparse_product_with_pruning_impl(rhs, lhs, _res); +// res = _res.transpose(); + } +}; + +// NOTE the 2 others cases (col row *) must never occur since they are caught +// by ProductReturnType which transforms it to (col col *) by evaluating rhs. + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_SPARSESPARSEPRODUCTWITHPRUNING_H diff --git a/eigen/Eigen/src/SparseCore/SparseTranspose.h b/eigen/Eigen/src/SparseCore/SparseTranspose.h new file mode 100644 index 0000000..76d031d --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseTranspose.h @@ -0,0 +1,63 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 Gael Guennebaud +// +// 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_SPARSETRANSPOSE_H +#define EIGEN_SPARSETRANSPOSE_H + +namespace Eigen { + +template class TransposeImpl + : public SparseMatrixBase > +{ + typedef typename internal::remove_all::type _MatrixTypeNested; + public: + + EIGEN_SPARSE_PUBLIC_INTERFACE(Transpose ) + + class InnerIterator; + class ReverseInnerIterator; + + inline Index nonZeros() const { return derived().nestedExpression().nonZeros(); } +}; + +// NOTE: VC10 and VC11 trigger an ICE if don't put typename TransposeImpl:: in front of Index, +// a typedef typename TransposeImpl::Index Index; +// does not fix the issue. +// An alternative is to define the nested class in the parent class itself. +template class TransposeImpl::InnerIterator + : public _MatrixTypeNested::InnerIterator +{ + typedef typename _MatrixTypeNested::InnerIterator Base; + typedef typename TransposeImpl::Index Index; + public: + + EIGEN_STRONG_INLINE InnerIterator(const TransposeImpl& trans, typename TransposeImpl::Index outer) + : Base(trans.derived().nestedExpression(), outer) + {} + typename TransposeImpl::Index row() const { return Base::col(); } + typename TransposeImpl::Index col() const { return Base::row(); } +}; + +template class TransposeImpl::ReverseInnerIterator + : public _MatrixTypeNested::ReverseInnerIterator +{ + typedef typename _MatrixTypeNested::ReverseInnerIterator Base; + typedef typename TransposeImpl::Index Index; + public: + + EIGEN_STRONG_INLINE ReverseInnerIterator(const TransposeImpl& xpr, typename TransposeImpl::Index outer) + : Base(xpr.derived().nestedExpression(), outer) + {} + typename TransposeImpl::Index row() const { return Base::col(); } + typename TransposeImpl::Index col() const { return Base::row(); } +}; + +} // end namespace Eigen + +#endif // EIGEN_SPARSETRANSPOSE_H diff --git a/eigen/Eigen/src/SparseCore/SparseTriangularView.h b/eigen/Eigen/src/SparseCore/SparseTriangularView.h new file mode 100644 index 0000000..333127b --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseTriangularView.h @@ -0,0 +1,179 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud +// Copyright (C) 2012 Désiré Nuentsa-Wakam +// +// 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_SPARSE_TRIANGULARVIEW_H +#define EIGEN_SPARSE_TRIANGULARVIEW_H + +namespace Eigen { + +namespace internal { + +template +struct traits > +: public traits +{}; + +} // namespace internal + +template class SparseTriangularView + : public SparseMatrixBase > +{ + enum { SkipFirst = ((Mode&Lower) && !(MatrixType::Flags&RowMajorBit)) + || ((Mode&Upper) && (MatrixType::Flags&RowMajorBit)), + SkipLast = !SkipFirst, + SkipDiag = (Mode&ZeroDiag) ? 1 : 0, + HasUnitDiag = (Mode&UnitDiag) ? 1 : 0 + }; + + public: + + EIGEN_SPARSE_PUBLIC_INTERFACE(SparseTriangularView) + + class InnerIterator; + class ReverseInnerIterator; + + inline Index rows() const { return m_matrix.rows(); } + inline Index cols() const { return m_matrix.cols(); } + + typedef typename MatrixType::Nested MatrixTypeNested; + typedef typename internal::remove_reference::type MatrixTypeNestedNonRef; + typedef typename internal::remove_all::type MatrixTypeNestedCleaned; + + inline SparseTriangularView(const MatrixType& matrix) : m_matrix(matrix) {} + + /** \internal */ + inline const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; } + + template + typename internal::plain_matrix_type_column_major::type + solve(const MatrixBase& other) const; + + template void solveInPlace(MatrixBase& other) const; + template void solveInPlace(SparseMatrixBase& other) const; + + protected: + MatrixTypeNested m_matrix; +}; + +template +class SparseTriangularView::InnerIterator : public MatrixTypeNestedCleaned::InnerIterator +{ + typedef typename MatrixTypeNestedCleaned::InnerIterator Base; + typedef typename SparseTriangularView::Index Index; + public: + + EIGEN_STRONG_INLINE InnerIterator(const SparseTriangularView& view, Index outer) + : Base(view.nestedExpression(), outer), m_returnOne(false) + { + if(SkipFirst) + { + while((*this) && ((HasUnitDiag||SkipDiag) ? this->index()<=outer : this->index()=Base::outer())) + { + if((!SkipFirst) && Base::operator bool()) + Base::operator++(); + m_returnOne = true; + } + } + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { + if(HasUnitDiag && m_returnOne) + m_returnOne = false; + else + { + Base::operator++(); + if(HasUnitDiag && (!SkipFirst) && ((!Base::operator bool()) || Base::index()>=Base::outer())) + { + if((!SkipFirst) && Base::operator bool()) + Base::operator++(); + m_returnOne = true; + } + } + return *this; + } + + inline Index row() const { return (MatrixType::Flags&RowMajorBit ? Base::outer() : this->index()); } + inline Index col() const { return (MatrixType::Flags&RowMajorBit ? this->index() : Base::outer()); } + inline Index index() const + { + if(HasUnitDiag && m_returnOne) return Base::outer(); + else return Base::index(); + } + inline Scalar value() const + { + if(HasUnitDiag && m_returnOne) return Scalar(1); + else return Base::value(); + } + + EIGEN_STRONG_INLINE operator bool() const + { + if(HasUnitDiag && m_returnOne) + return true; + if(SkipFirst) return Base::operator bool(); + else + { + if (SkipDiag) return (Base::operator bool() && this->index() < this->outer()); + else return (Base::operator bool() && this->index() <= this->outer()); + } + } + protected: + bool m_returnOne; +}; + +template +class SparseTriangularView::ReverseInnerIterator : public MatrixTypeNestedCleaned::ReverseInnerIterator +{ + typedef typename MatrixTypeNestedCleaned::ReverseInnerIterator Base; + typedef typename SparseTriangularView::Index Index; + public: + + EIGEN_STRONG_INLINE ReverseInnerIterator(const SparseTriangularView& view, Index outer) + : Base(view.nestedExpression(), outer) + { + eigen_assert((!HasUnitDiag) && "ReverseInnerIterator does not support yet triangular views with a unit diagonal"); + if(SkipLast) { + while((*this) && (SkipDiag ? this->index()>=outer : this->index()>outer)) + --(*this); + } + } + + EIGEN_STRONG_INLINE ReverseInnerIterator& operator--() + { Base::operator--(); return *this; } + + inline Index row() const { return Base::row(); } + inline Index col() const { return Base::col(); } + + EIGEN_STRONG_INLINE operator bool() const + { + if (SkipLast) return Base::operator bool() ; + else + { + if(SkipDiag) return (Base::operator bool() && this->index() > this->outer()); + else return (Base::operator bool() && this->index() >= this->outer()); + } + } +}; + +template +template +inline const SparseTriangularView +SparseMatrixBase::triangularView() const +{ + return derived(); +} + +} // end namespace Eigen + +#endif // EIGEN_SPARSE_TRIANGULARVIEW_H diff --git a/eigen/Eigen/src/SparseCore/SparseUtil.h b/eigen/Eigen/src/SparseCore/SparseUtil.h new file mode 100644 index 0000000..d627546 --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseUtil.h @@ -0,0 +1,172 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// +// 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_SPARSEUTIL_H +#define EIGEN_SPARSEUTIL_H + +namespace Eigen { + +#ifdef NDEBUG +#define EIGEN_DBG_SPARSE(X) +#else +#define EIGEN_DBG_SPARSE(X) X +#endif + +#define EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(Derived, Op) \ +template \ +EIGEN_STRONG_INLINE Derived& operator Op(const Eigen::SparseMatrixBase& other) \ +{ \ + return Base::operator Op(other.derived()); \ +} \ +EIGEN_STRONG_INLINE Derived& operator Op(const Derived& other) \ +{ \ + return Base::operator Op(other); \ +} + +#define EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, Op) \ +template \ +EIGEN_STRONG_INLINE Derived& operator Op(const Other& scalar) \ +{ \ + return Base::operator Op(scalar); \ +} + +#define EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATORS(Derived) \ +EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(Derived, =) \ +EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(Derived, +=) \ +EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(Derived, -=) \ +EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, *=) \ +EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, /=) + +#define _EIGEN_SPARSE_PUBLIC_INTERFACE(Derived, BaseClass) \ + typedef BaseClass Base; \ + typedef typename Eigen::internal::traits::Scalar Scalar; \ + typedef typename Eigen::NumTraits::Real RealScalar; \ + typedef typename Eigen::internal::nested::type Nested; \ + typedef typename Eigen::internal::traits::StorageKind StorageKind; \ + typedef typename Eigen::internal::traits::Index Index; \ + enum { RowsAtCompileTime = Eigen::internal::traits::RowsAtCompileTime, \ + ColsAtCompileTime = Eigen::internal::traits::ColsAtCompileTime, \ + Flags = Eigen::internal::traits::Flags, \ + CoeffReadCost = Eigen::internal::traits::CoeffReadCost, \ + SizeAtCompileTime = Base::SizeAtCompileTime, \ + IsVectorAtCompileTime = Base::IsVectorAtCompileTime }; \ + using Base::derived; \ + using Base::const_cast_derived; + +#define EIGEN_SPARSE_PUBLIC_INTERFACE(Derived) \ + _EIGEN_SPARSE_PUBLIC_INTERFACE(Derived, Eigen::SparseMatrixBase) + +const int CoherentAccessPattern = 0x1; +const int InnerRandomAccessPattern = 0x2 | CoherentAccessPattern; +const int OuterRandomAccessPattern = 0x4 | CoherentAccessPattern; +const int RandomAccessPattern = 0x8 | OuterRandomAccessPattern | InnerRandomAccessPattern; + +template class SparseMatrix; +template class DynamicSparseMatrix; +template class SparseVector; +template class MappedSparseMatrix; + +template class SparseTriangularView; +template class SparseSelfAdjointView; +template class SparseDiagonalProduct; +template class SparseView; + +template class SparseSparseProduct; +template class SparseTimeDenseProduct; +template class DenseTimeSparseProduct; +template class SparseDenseOuterProduct; + +template struct SparseSparseProductReturnType; +template::ColsAtCompileTime,internal::traits::RowsAtCompileTime)> struct DenseSparseProductReturnType; +template::ColsAtCompileTime,internal::traits::RowsAtCompileTime)> struct SparseDenseProductReturnType; +template class SparseSymmetricPermutationProduct; + +namespace internal { + +template struct sparse_eval; + +template struct eval + : public sparse_eval::RowsAtCompileTime,traits::ColsAtCompileTime> +{}; + +template struct sparse_eval { + typedef typename traits::Scalar _Scalar; + typedef typename traits::Index _Index; + public: + typedef SparseVector<_Scalar, RowMajor, _Index> type; +}; + +template struct sparse_eval { + typedef typename traits::Scalar _Scalar; + typedef typename traits::Index _Index; + public: + typedef SparseVector<_Scalar, ColMajor, _Index> type; +}; + +template struct sparse_eval { + typedef typename traits::Scalar _Scalar; + typedef typename traits::Index _Index; + enum { _Options = ((traits::Flags&RowMajorBit)==RowMajorBit) ? RowMajor : ColMajor }; + public: + typedef SparseMatrix<_Scalar, _Options, _Index> type; +}; + +template struct sparse_eval { + typedef typename traits::Scalar _Scalar; + public: + typedef Matrix<_Scalar, 1, 1> type; +}; + +template struct plain_matrix_type +{ + typedef typename traits::Scalar _Scalar; + typedef typename traits::Index _Index; + enum { _Options = ((traits::Flags&RowMajorBit)==RowMajorBit) ? RowMajor : ColMajor }; + public: + typedef SparseMatrix<_Scalar, _Options, _Index> type; +}; + +} // end namespace internal + +/** \ingroup SparseCore_Module + * + * \class Triplet + * + * \brief A small structure to hold a non zero as a triplet (i,j,value). + * + * \sa SparseMatrix::setFromTriplets() + */ +template::Index > +class Triplet +{ +public: + Triplet() : m_row(0), m_col(0), m_value(0) {} + + Triplet(const Index& i, const Index& j, const Scalar& v = Scalar(0)) + : m_row(i), m_col(j), m_value(v) + {} + + /** \returns the row index of the element */ + const Index& row() const { return m_row; } + + /** \returns the column index of the element */ + const Index& col() const { return m_col; } + + /** \returns the value of the element */ + const Scalar& value() const { return m_value; } +protected: + Index m_row, m_col; + Scalar m_value; +}; + +} // end namespace Eigen + +#endif // EIGEN_SPARSEUTIL_H diff --git a/eigen/Eigen/src/SparseCore/SparseVector.h b/eigen/Eigen/src/SparseCore/SparseVector.h new file mode 100644 index 0000000..c7ee89c --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseVector.h @@ -0,0 +1,448 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 Gael Guennebaud +// +// 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_SPARSEVECTOR_H +#define EIGEN_SPARSEVECTOR_H + +namespace Eigen { + +/** \ingroup SparseCore_Module + * \class SparseVector + * + * \brief a sparse vector class + * + * \tparam _Scalar the scalar type, i.e. the type of the coefficients + * + * See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme. + * + * This class can be extended with the help of the plugin mechanism described on the page + * \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_SPARSEVECTOR_PLUGIN. + */ + +namespace internal { +template +struct traits > +{ + typedef _Scalar Scalar; + typedef _Index Index; + typedef Sparse StorageKind; + typedef MatrixXpr XprKind; + enum { + IsColVector = (_Options & RowMajorBit) ? 0 : 1, + + RowsAtCompileTime = IsColVector ? Dynamic : 1, + ColsAtCompileTime = IsColVector ? 1 : Dynamic, + MaxRowsAtCompileTime = RowsAtCompileTime, + MaxColsAtCompileTime = ColsAtCompileTime, + Flags = _Options | NestByRefBit | LvalueBit | (IsColVector ? 0 : RowMajorBit), + CoeffReadCost = NumTraits::ReadCost, + SupportedAccessPatterns = InnerRandomAccessPattern + }; +}; + +// Sparse-Vector-Assignment kinds: +enum { + SVA_RuntimeSwitch, + SVA_Inner, + SVA_Outer +}; + +template< typename Dest, typename Src, + int AssignmentKind = !bool(Src::IsVectorAtCompileTime) ? SVA_RuntimeSwitch + : Src::InnerSizeAtCompileTime==1 ? SVA_Outer + : SVA_Inner> +struct sparse_vector_assign_selector; + +} + +template +class SparseVector + : public SparseMatrixBase > +{ + typedef SparseMatrixBase SparseBase; + + public: + EIGEN_SPARSE_PUBLIC_INTERFACE(SparseVector) + EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseVector, +=) + EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseVector, -=) + + typedef internal::CompressedStorage Storage; + enum { IsColVector = internal::traits::IsColVector }; + + enum { + Options = _Options + }; + + EIGEN_STRONG_INLINE Index rows() const { return IsColVector ? m_size : 1; } + EIGEN_STRONG_INLINE Index cols() const { return IsColVector ? 1 : m_size; } + EIGEN_STRONG_INLINE Index innerSize() const { return m_size; } + EIGEN_STRONG_INLINE Index outerSize() const { return 1; } + + EIGEN_STRONG_INLINE const Scalar* valuePtr() const { return m_data.valuePtr(); } + EIGEN_STRONG_INLINE Scalar* valuePtr() { return m_data.valuePtr(); } + + EIGEN_STRONG_INLINE const Index* innerIndexPtr() const { return m_data.indexPtr(); } + EIGEN_STRONG_INLINE Index* innerIndexPtr() { return m_data.indexPtr(); } + + /** \internal */ + inline Storage& data() { return m_data; } + /** \internal */ + inline const Storage& data() const { return m_data; } + + inline Scalar coeff(Index row, Index col) const + { + eigen_assert(IsColVector ? (col==0 && row>=0 && row=0 && col=0 && i=0 && row=0 && col=0 && i(m_data.size()); } + + inline void startVec(Index outer) + { + EIGEN_UNUSED_VARIABLE(outer); + eigen_assert(outer==0); + } + + inline Scalar& insertBackByOuterInner(Index outer, Index inner) + { + EIGEN_UNUSED_VARIABLE(outer); + eigen_assert(outer==0); + return insertBack(inner); + } + inline Scalar& insertBack(Index i) + { + m_data.append(0, i); + return m_data.value(m_data.size()-1); + } + + inline Scalar& insert(Index row, Index col) + { + eigen_assert(IsColVector ? (col==0 && row>=0 && row=0 && col=0 && i= startId) && (m_data.index(p) > i) ) + { + m_data.index(p+1) = m_data.index(p); + m_data.value(p+1) = m_data.value(p); + --p; + } + m_data.index(p+1) = i; + m_data.value(p+1) = 0; + return m_data.value(p+1); + } + + /** + */ + inline void reserve(Index reserveSize) { m_data.reserve(reserveSize); } + + + inline void finalize() {} + + void prune(const Scalar& reference, const RealScalar& epsilon = NumTraits::dummy_precision()) + { + m_data.prune(reference,epsilon); + } + + void resize(Index rows, Index cols) + { + eigen_assert(rows==1 || cols==1); + resize(IsColVector ? rows : cols); + } + + void resize(Index newSize) + { + m_size = newSize; + m_data.clear(); + } + + void resizeNonZeros(Index size) { m_data.resize(size); } + + inline SparseVector() : m_size(0) { check_template_parameters(); resize(0); } + + inline SparseVector(Index size) : m_size(0) { check_template_parameters(); resize(size); } + + inline SparseVector(Index rows, Index cols) : m_size(0) { check_template_parameters(); resize(rows,cols); } + + template + inline SparseVector(const SparseMatrixBase& other) + : m_size(0) + { + check_template_parameters(); + *this = other.derived(); + } + + inline SparseVector(const SparseVector& other) + : SparseBase(other), m_size(0) + { + check_template_parameters(); + *this = other.derived(); + } + + /** Swaps the values of \c *this and \a other. + * Overloaded for performance: this version performs a \em shallow swap by swaping pointers and attributes only. + * \sa SparseMatrixBase::swap() + */ + inline void swap(SparseVector& other) + { + std::swap(m_size, other.m_size); + m_data.swap(other.m_data); + } + + inline SparseVector& operator=(const SparseVector& other) + { + if (other.isRValue()) + { + swap(other.const_cast_derived()); + } + else + { + resize(other.size()); + m_data = other.m_data; + } + return *this; + } + + template + inline SparseVector& operator=(const SparseMatrixBase& other) + { + SparseVector tmp(other.size()); + internal::sparse_vector_assign_selector::run(tmp,other.derived()); + this->swap(tmp); + return *this; + } + + #ifndef EIGEN_PARSED_BY_DOXYGEN + template + inline SparseVector& operator=(const SparseSparseProduct& product) + { + return Base::operator=(product); + } + #endif + + friend std::ostream & operator << (std::ostream & s, const SparseVector& m) + { + for (Index i=0; i::IsSigned,THE_INDEX_TYPE_MUST_BE_A_SIGNED_TYPE); + EIGEN_STATIC_ASSERT((_Options&(ColMajor|RowMajor))==Options,INVALID_MATRIX_TEMPLATE_PARAMETERS); + } + + Storage m_data; + Index m_size; +}; + +template +class SparseVector::InnerIterator +{ + public: + InnerIterator(const SparseVector& vec, Index outer=0) + : m_data(vec.m_data), m_id(0), m_end(static_cast(m_data.size())) + { + EIGEN_UNUSED_VARIABLE(outer); + eigen_assert(outer==0); + } + + InnerIterator(const internal::CompressedStorage& data) + : m_data(data), m_id(0), m_end(static_cast(m_data.size())) + {} + + inline InnerIterator& operator++() { m_id++; return *this; } + + inline Scalar value() const { return m_data.value(m_id); } + inline Scalar& valueRef() { return const_cast(m_data.value(m_id)); } + + inline Index index() const { return m_data.index(m_id); } + inline Index row() const { return IsColVector ? index() : 0; } + inline Index col() const { return IsColVector ? 0 : index(); } + + inline operator bool() const { return (m_id < m_end); } + + protected: + const internal::CompressedStorage& m_data; + Index m_id; + const Index m_end; +}; + +template +class SparseVector::ReverseInnerIterator +{ + public: + ReverseInnerIterator(const SparseVector& vec, Index outer=0) + : m_data(vec.m_data), m_id(static_cast(m_data.size())), m_start(0) + { + EIGEN_UNUSED_VARIABLE(outer); + eigen_assert(outer==0); + } + + ReverseInnerIterator(const internal::CompressedStorage& data) + : m_data(data), m_id(static_cast(m_data.size())), m_start(0) + {} + + inline ReverseInnerIterator& operator--() { m_id--; return *this; } + + inline Scalar value() const { return m_data.value(m_id-1); } + inline Scalar& valueRef() { return const_cast(m_data.value(m_id-1)); } + + inline Index index() const { return m_data.index(m_id-1); } + inline Index row() const { return IsColVector ? index() : 0; } + inline Index col() const { return IsColVector ? 0 : index(); } + + inline operator bool() const { return (m_id > m_start); } + + protected: + const internal::CompressedStorage& m_data; + Index m_id; + const Index m_start; +}; + +namespace internal { + +template< typename Dest, typename Src> +struct sparse_vector_assign_selector { + static void run(Dest& dst, const Src& src) { + eigen_internal_assert(src.innerSize()==src.size()); + for(typename Src::InnerIterator it(src, 0); it; ++it) + dst.insert(it.index()) = it.value(); + } +}; + +template< typename Dest, typename Src> +struct sparse_vector_assign_selector { + static void run(Dest& dst, const Src& src) { + eigen_internal_assert(src.outerSize()==src.size()); + for(typename Dest::Index i=0; i +struct sparse_vector_assign_selector { + static void run(Dest& dst, const Src& src) { + if(src.outerSize()==1) sparse_vector_assign_selector::run(dst, src); + else sparse_vector_assign_selector::run(dst, src); + } +}; + +} + +} // end namespace Eigen + +#endif // EIGEN_SPARSEVECTOR_H diff --git a/eigen/Eigen/src/SparseCore/SparseView.h b/eigen/Eigen/src/SparseCore/SparseView.h new file mode 100644 index 0000000..2820b39 --- /dev/null +++ b/eigen/Eigen/src/SparseCore/SparseView.h @@ -0,0 +1,99 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2011 Gael Guennebaud +// Copyright (C) 2010 Daniel Lowengrub +// +// 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_SPARSEVIEW_H +#define EIGEN_SPARSEVIEW_H + +namespace Eigen { + +namespace internal { + +template +struct traits > : traits +{ + typedef typename MatrixType::Index Index; + typedef Sparse StorageKind; + enum { + Flags = int(traits::Flags) & (RowMajorBit) + }; +}; + +} // end namespace internal + +template +class SparseView : public SparseMatrixBase > +{ + typedef typename MatrixType::Nested MatrixTypeNested; + typedef typename internal::remove_all::type _MatrixTypeNested; +public: + EIGEN_SPARSE_PUBLIC_INTERFACE(SparseView) + + explicit SparseView(const MatrixType& mat, const Scalar& reference = Scalar(0), + const RealScalar &epsilon = NumTraits::dummy_precision()) + : m_matrix(mat), m_reference(reference), m_epsilon(epsilon) {} + + class InnerIterator; + + inline Index rows() const { return m_matrix.rows(); } + inline Index cols() const { return m_matrix.cols(); } + + inline Index innerSize() const { return m_matrix.innerSize(); } + inline Index outerSize() const { return m_matrix.outerSize(); } + +protected: + MatrixTypeNested m_matrix; + Scalar m_reference; + typename NumTraits::Real m_epsilon; +}; + +template +class SparseView::InnerIterator : public _MatrixTypeNested::InnerIterator +{ + typedef typename SparseView::Index Index; +public: + typedef typename _MatrixTypeNested::InnerIterator IterBase; + InnerIterator(const SparseView& view, Index outer) : + IterBase(view.m_matrix, outer), m_view(view) + { + incrementToNonZero(); + } + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { + IterBase::operator++(); + incrementToNonZero(); + return *this; + } + + using IterBase::value; + +protected: + const SparseView& m_view; + +private: + void incrementToNonZero() + { + while((bool(*this)) && internal::isMuchSmallerThan(value(), m_view.m_reference, m_view.m_epsilon)) + { + IterBase::operator++(); + } + } +}; + +template +const SparseView MatrixBase::sparseView(const Scalar& m_reference, + const typename NumTraits::Real& m_epsilon) const +{ + return SparseView(derived(), m_reference, m_epsilon); +} + +} // end namespace Eigen + +#endif diff --git a/eigen/Eigen/src/SparseCore/TriangularSolver.h b/eigen/Eigen/src/SparseCore/TriangularSolver.h new file mode 100644 index 0000000..ccc12af --- /dev/null +++ b/eigen/Eigen/src/SparseCore/TriangularSolver.h @@ -0,0 +1,334 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud +// +// 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_SPARSETRIANGULARSOLVER_H +#define EIGEN_SPARSETRIANGULARSOLVER_H + +namespace Eigen { + +namespace internal { + +template::Flags) & RowMajorBit> +struct sparse_solve_triangular_selector; + +// forward substitution, row-major +template +struct sparse_solve_triangular_selector +{ + typedef typename Rhs::Scalar Scalar; + static void run(const Lhs& lhs, Rhs& other) + { + for(int col=0 ; col +struct sparse_solve_triangular_selector +{ + typedef typename Rhs::Scalar Scalar; + static void run(const Lhs& lhs, Rhs& other) + { + for(int col=0 ; col=0 ; --i) + { + Scalar tmp = other.coeff(i,col); + Scalar l_ii(0); + typename Lhs::InnerIterator it(lhs, i); + while(it && it.index() +struct sparse_solve_triangular_selector +{ + typedef typename Rhs::Scalar Scalar; + static void run(const Lhs& lhs, Rhs& other) + { + for(int col=0 ; col +struct sparse_solve_triangular_selector +{ + typedef typename Rhs::Scalar Scalar; + static void run(const Lhs& lhs, Rhs& other) + { + for(int col=0 ; col=0; --i) + { + Scalar& tmp = other.coeffRef(i,col); + if (tmp!=Scalar(0)) // optimization when other is actually sparse + { + if(!(Mode & UnitDiag)) + { + // TODO replace this by a binary search. make sure the binary search is safe for partially sorted elements + typename Lhs::ReverseInnerIterator it(lhs, i); + while(it && it.index()!=i) + --it; + eigen_assert(it && it.index()==i); + other.coeffRef(i,col) /= it.value(); + } + typename Lhs::InnerIterator it(lhs, i); + for(; it && it.index() +template +void SparseTriangularView::solveInPlace(MatrixBase& other) const +{ + eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.rows()); + eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower))); + + enum { copy = internal::traits::Flags & RowMajorBit }; + + typedef typename internal::conditional::type, OtherDerived&>::type OtherCopy; + OtherCopy otherCopy(other.derived()); + + internal::sparse_solve_triangular_selector::type, Mode>::run(m_matrix, otherCopy); + + if (copy) + other = otherCopy; +} + +template +template +typename internal::plain_matrix_type_column_major::type +SparseTriangularView::solve(const MatrixBase& other) const +{ + typename internal::plain_matrix_type_column_major::type res(other); + solveInPlace(res); + return res; +} + +// pure sparse path + +namespace internal { + +template +struct sparse_solve_triangular_sparse_selector; + +// forward substitution, col-major +template +struct sparse_solve_triangular_sparse_selector +{ + typedef typename Rhs::Scalar Scalar; + typedef typename promote_index_type::Index, + typename traits::Index>::type Index; + static void run(const Lhs& lhs, Rhs& other) + { + const bool IsLower = (UpLo==Lower); + AmbiVector tempVector(other.rows()*2); + tempVector.setBounds(0,other.rows()); + + Rhs res(other.rows(), other.cols()); + res.reserve(other.nonZeros()); + + for(int col=0 ; col=0; + i+=IsLower?1:-1) + { + tempVector.restart(); + Scalar& ci = tempVector.coeffRef(i); + if (ci!=Scalar(0)) + { + // find + typename Lhs::InnerIterator it(lhs, i); + if(!(Mode & UnitDiag)) + { + if (IsLower) + { + eigen_assert(it.index()==i); + ci /= it.value(); + } + else + ci /= lhs.coeff(i,i); + } + tempVector.restart(); + if (IsLower) + { + if (it.index()==i) + ++it; + for(; it; ++it) + tempVector.coeffRef(it.index()) -= ci * it.value(); + } + else + { + for(; it && it.index()::Iterator it(tempVector/*,1e-12*/); it; ++it) + { + ++ count; +// std::cerr << "fill " << it.index() << ", " << col << "\n"; +// std::cout << it.value() << " "; + // FIXME use insertBack + res.insert(it.index(), col) = it.value(); + } +// std::cout << "tempVector.nonZeros() == " << int(count) << " / " << (other.rows()) << "\n"; + } + res.finalize(); + other = res.markAsRValue(); + } +}; + +} // end namespace internal + +template +template +void SparseTriangularView::solveInPlace(SparseMatrixBase& other) const +{ + eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.rows()); + eigen_assert( (!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower))); + +// enum { copy = internal::traits::Flags & RowMajorBit }; + +// typedef typename internal::conditional::type, OtherDerived&>::type OtherCopy; +// OtherCopy otherCopy(other.derived()); + + internal::sparse_solve_triangular_sparse_selector::run(m_matrix, other.derived()); + +// if (copy) +// other = otherCopy; +} + +#ifdef EIGEN2_SUPPORT + +// deprecated stuff: + +/** \deprecated */ +template +template +void SparseMatrixBase::solveTriangularInPlace(MatrixBase& other) const +{ + this->template triangular().solveInPlace(other); +} + +/** \deprecated */ +template +template +typename internal::plain_matrix_type_column_major::type +SparseMatrixBase::solveTriangular(const MatrixBase& other) const +{ + typename internal::plain_matrix_type_column_major::type res(other); + derived().solveTriangularInPlace(res); + return res; +} +#endif // EIGEN2_SUPPORT + +} // end namespace Eigen + +#endif // EIGEN_SPARSETRIANGULARSOLVER_H -- cgit v1.2.3