From f0238cfb6997c4acfc2bd200de7295f3fa36968f Mon Sep 17 00:00:00 2001 From: Stanislaw Halik Date: Sun, 3 Mar 2019 21:09:10 +0100 Subject: don't index Eigen --- eigen/Eigen/src/SparseCore/AmbiVector.h | 377 ------ eigen/Eigen/src/SparseCore/CompressedStorage.h | 258 ---- .../SparseCore/ConservativeSparseSparseProduct.h | 352 ----- eigen/Eigen/src/SparseCore/MappedSparseMatrix.h | 67 - eigen/Eigen/src/SparseCore/SparseAssign.h | 216 --- eigen/Eigen/src/SparseCore/SparseBlock.h | 603 --------- eigen/Eigen/src/SparseCore/SparseColEtree.h | 206 --- eigen/Eigen/src/SparseCore/SparseCompressedBase.h | 341 ----- eigen/Eigen/src/SparseCore/SparseCwiseBinaryOp.h | 726 ---------- eigen/Eigen/src/SparseCore/SparseCwiseUnaryOp.h | 148 --- eigen/Eigen/src/SparseCore/SparseDenseProduct.h | 320 ----- eigen/Eigen/src/SparseCore/SparseDiagonalProduct.h | 138 -- eigen/Eigen/src/SparseCore/SparseDot.h | 98 -- eigen/Eigen/src/SparseCore/SparseFuzzy.h | 29 - eigen/Eigen/src/SparseCore/SparseMap.h | 305 ----- eigen/Eigen/src/SparseCore/SparseMatrix.h | 1403 -------------------- eigen/Eigen/src/SparseCore/SparseMatrixBase.h | 405 ------ eigen/Eigen/src/SparseCore/SparsePermutation.h | 178 --- eigen/Eigen/src/SparseCore/SparseProduct.h | 169 --- eigen/Eigen/src/SparseCore/SparseRedux.h | 49 - eigen/Eigen/src/SparseCore/SparseRef.h | 397 ------ eigen/Eigen/src/SparseCore/SparseSelfAdjointView.h | 656 --------- eigen/Eigen/src/SparseCore/SparseSolverBase.h | 124 -- .../SparseCore/SparseSparseProductWithPruning.h | 198 --- eigen/Eigen/src/SparseCore/SparseTranspose.h | 92 -- eigen/Eigen/src/SparseCore/SparseTriangularView.h | 189 --- eigen/Eigen/src/SparseCore/SparseUtil.h | 178 --- eigen/Eigen/src/SparseCore/SparseVector.h | 478 ------- eigen/Eigen/src/SparseCore/SparseView.h | 253 ---- eigen/Eigen/src/SparseCore/TriangularSolver.h | 315 ----- 30 files changed, 9268 deletions(-) delete mode 100644 eigen/Eigen/src/SparseCore/AmbiVector.h delete mode 100644 eigen/Eigen/src/SparseCore/CompressedStorage.h delete mode 100644 eigen/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h delete mode 100644 eigen/Eigen/src/SparseCore/MappedSparseMatrix.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseAssign.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseBlock.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseColEtree.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseCompressedBase.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseCwiseBinaryOp.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseCwiseUnaryOp.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseDenseProduct.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseDiagonalProduct.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseDot.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseFuzzy.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseMap.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseMatrix.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseMatrixBase.h delete mode 100644 eigen/Eigen/src/SparseCore/SparsePermutation.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseProduct.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseRedux.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseRef.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseSelfAdjointView.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseSolverBase.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseSparseProductWithPruning.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseTranspose.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseTriangularView.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseUtil.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseVector.h delete mode 100644 eigen/Eigen/src/SparseCore/SparseView.h delete 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 deleted file mode 100644 index e0295f2..0000000 --- a/eigen/Eigen/src/SparseCore/AmbiVector.h +++ /dev/null @@ -1,377 +0,0 @@ -// 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 _StorageIndex StorageIndex; - typedef typename NumTraits::Real RealScalar; - - explicit 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 = convert_index(start); m_end = convert_index(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 = convert_index(size); - } - - StorageIndex size() const { return m_size; } - - protected: - StorageIndex convert_index(Index idx) - { - return internal::convert_index(idx); - } - - 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 = convert_index((allocSize*sizeof(Scalar))/sizeof(ListEl)); - m_buffer = new Scalar[allocSize]; - } - else - { - m_allocatedElements = convert_index((size*sizeof(Scalar))/sizeof(ListEl)); - m_buffer = new Scalar[size]; - } - m_size = convert_index(size); - m_start = 0; - m_end = m_size; - } - - void reallocateSparse() - { - Index copyElements = m_allocatedElements; - m_allocatedElements = (std::min)(StorageIndex(m_allocatedElements*1.5),m_size); - Index allocSize = m_allocatedElements * sizeof(ListEl); - allocSize = (allocSize + sizeof(Scalar) - 1)/sizeof(Scalar); - Scalar* newBuffer = new Scalar[allocSize]; - std::memcpy(newBuffer, m_buffer, copyElements * sizeof(ListEl)); - delete[] m_buffer; - m_buffer = newBuffer; - } - - protected: - // element type of the linked list - struct ListEl - { - StorageIndex next; - StorageIndex index; - Scalar value; - }; - - // used to store data in both mode - Scalar* m_buffer; - Scalar m_zero; - StorageIndex m_size; - StorageIndex m_start; - StorageIndex m_end; - StorageIndex m_allocatedSize; - StorageIndex m_allocatedElements; - StorageIndex m_mode; - - // linked list mode - StorageIndex m_llStart; - StorageIndex m_llCurrent; - StorageIndex m_llSize; -}; - -/** \returns the number of non zeros in the current sub vector */ -template -Index AmbiVector<_Scalar,_StorageIndex>::nonZeros() const -{ - if (m_mode==IsSparse) - return m_llSize; - else - return m_end - m_start; -} - -template -void AmbiVector<_Scalar,_StorageIndex>::init(double estimatedDensity) -{ - if (estimatedDensity>0.1) - init(IsDense); - else - init(IsSparse); -} - -template -void AmbiVector<_Scalar,_StorageIndex>::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,_StorageIndex>::restart() -{ - m_llCurrent = m_llStart; -} - -/** Set all coefficients of current subvector to zero */ -template -void AmbiVector<_Scalar,_StorageIndex>::setZero() -{ - if (m_mode==IsDense) - { - for (Index i=m_start; i -_Scalar& AmbiVector<_Scalar,_StorageIndex>::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 = convert_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,_StorageIndex>::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,_StorageIndex>::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. - */ - explicit 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; - } - } - } - - StorageIndex 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)<=m_epsilon); - if (m_currentEl<0) - { - m_cachedIndex = -1; - } - else - { - m_cachedIndex = llElements[m_currentEl].index; - m_cachedValue = llElements[m_currentEl].value; - } - } - return *this; - } - - protected: - const AmbiVector& m_vector; // the target vector - StorageIndex m_currentEl; // the current element in sparse/linked-list mode - RealScalar m_epsilon; // epsilon used to prune zero coefficients - StorageIndex m_cachedIndex; // current coordinate - Scalar m_cachedValue; // current value - bool m_isDense; // mode of the vector -}; - -} // end namespace internal - -} // end namespace Eigen - -#endif // EIGEN_AMBIVECTOR_H diff --git a/eigen/Eigen/src/SparseCore/CompressedStorage.h b/eigen/Eigen/src/SparseCore/CompressedStorage.h deleted file mode 100644 index d89fa0d..0000000 --- a/eigen/Eigen/src/SparseCore/CompressedStorage.h +++ /dev/null @@ -1,258 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2014 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_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 _StorageIndex StorageIndex; - - protected: - - typedef typename NumTraits::Real RealScalar; - - public: - - CompressedStorage() - : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0) - {} - - explicit CompressedStorage(Index 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()); - if(other.size()>0) - { - internal::smart_copy(other.m_values, other.m_values + m_size, m_values); - internal::smart_copy(other.m_indices, other.m_indices + m_size, m_indices); - } - 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(Index size) - { - Index newAllocatedSize = m_size + size; - if (newAllocatedSize > m_allocatedSize) - reallocate(newAllocatedSize); - } - - void squeeze() - { - if (m_allocatedSize>m_size) - reallocate(m_size); - } - - void resize(Index size, double reserveSizeFactor = 0) - { - if (m_allocatedSize)(NumTraits::highest(), size + Index(reserveSizeFactor*double(size))); - if(realloc_size(i); - } - - inline Index size() const { return m_size; } - inline Index allocatedSize() const { return m_allocatedSize; } - inline void clear() { m_size = 0; } - - const Scalar* valuePtr() const { return m_values; } - Scalar* valuePtr() { return m_values; } - const StorageIndex* indexPtr() const { return m_indices; } - StorageIndex* indexPtr() { return m_indices; } - - inline Scalar& value(Index i) { eigen_internal_assert(m_values!=0); return m_values[i]; } - inline const Scalar& value(Index i) const { eigen_internal_assert(m_values!=0); return m_values[i]; } - - inline StorageIndex& index(Index i) { eigen_internal_assert(m_indices!=0); return m_indices[i]; } - inline const StorageIndex& index(Index i) const { eigen_internal_assert(m_indices!=0); return m_indices[i]; } - - /** \returns the largest \c k such that for all \c j in [0,k) index[\c j]\<\a key */ - inline Index searchLowerIndex(Index key) const - { - return searchLowerIndex(0, m_size, key); - } - - /** \returns the largest \c k in [start,end) such that for all \c j in [start,k) index[\c j]\<\a key */ - inline Index searchLowerIndex(Index start, Index end, Index key) const - { - while(end>start) - { - Index mid = (end+start)>>1; - if (m_indices[mid]=end) - return defaultValue; - 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 Index id = searchLowerIndex(start,end-1,key); - return ((id=m_size || m_indices[id]!=key) - { - if (m_allocatedSize newValues(m_allocatedSize); - internal::scoped_array newIndices(m_allocatedSize); - - // copy first chunk - internal::smart_copy(m_values, m_values +id, newValues.ptr()); - internal::smart_copy(m_indices, m_indices+id, newIndices.ptr()); - - // copy the rest - if(m_size>id) - { - internal::smart_copy(m_values +id, m_values +m_size, newValues.ptr() +id+1); - internal::smart_copy(m_indices+id, m_indices+m_size, newIndices.ptr()+id+1); - } - std::swap(m_values,newValues.ptr()); - std::swap(m_indices,newIndices.ptr()); - } - else if(m_size>id) - { - internal::smart_memmove(m_values +id, m_values +m_size, m_values +id+1); - internal::smart_memmove(m_indices+id, m_indices+m_size, m_indices+id+1); - } - m_size++; - m_indices[id] = internal::convert_index(key); - m_values[id] = defaultValue; - } - return m_values[id]; - } - - void prune(const Scalar& reference, const RealScalar& epsilon = NumTraits::dummy_precision()) - { - Index k = 0; - Index n = size(); - for (Index i=0; i newValues(size); - internal::scoped_array newIndices(size); - Index copySize = (std::min)(size, m_size); - if (copySize>0) { - internal::smart_copy(m_values, m_values+copySize, newValues.ptr()); - internal::smart_copy(m_indices, m_indices+copySize, newIndices.ptr()); - } - std::swap(m_values,newValues.ptr()); - std::swap(m_indices,newIndices.ptr()); - m_allocatedSize = size; - } - - protected: - Scalar* m_values; - StorageIndex* m_indices; - Index m_size; - Index 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 deleted file mode 100644 index 9db119b..0000000 --- a/eigen/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h +++ /dev/null @@ -1,352 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2015 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, bool sortedInsertion = false) -{ - typedef typename remove_all::type::Scalar LhsScalar; - typedef typename remove_all::type::Scalar RhsScalar; - typedef typename remove_all::type::Scalar ResScalar; - - // 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()); - - ei_declare_aligned_stack_constructed_variable(bool, mask, rows, 0); - ei_declare_aligned_stack_constructed_variable(ResScalar, values, rows, 0); - ei_declare_aligned_stack_constructed_variable(Index, indices, rows, 0); - - std::memset(mask,0,sizeof(bool)*rows); - - evaluator lhsEval(lhs); - evaluator rhsEval(rhs); - - // 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 = lhsEval.nonZerosEstimate() + rhsEval.nonZerosEstimate(); - - res.setZero(); - res.reserve(Index(estimated_nnz_prod)); - // we compute each column of the result, one after the other - for (Index j=0; j::InnerIterator rhsIt(rhsEval, j); rhsIt; ++rhsIt) - { - RhsScalar y = rhsIt.value(); - Index k = rhsIt.index(); - for (typename evaluator::InnerIterator lhsIt(lhsEval, k); lhsIt; ++lhsIt) - { - Index i = lhsIt.index(); - LhsScalar x = lhsIt.value(); - if(!mask[i]) - { - mask[i] = true; - values[i] = x * y; - indices[nnz] = i; - ++nnz; - } - else - values[i] += x * y; - } - } - if(!sortedInsertion) - { - // unordered insertion - for(Index k=0; k 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,indices+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 ColMajorMatrixAux; - typedef typename sparse_eval::type ColMajorMatrix; - - // If the result is tall and thin (in the extreme case a column vector) - // then it is faster to sort the coefficients inplace instead of transposing twice. - // FIXME, the following heuristic is probably not very good. - if(lhs.rows()>rhs.cols()) - { - ColMajorMatrix resCol(lhs.rows(),rhs.cols()); - // perform sorted insertion - internal::conservative_sparse_sparse_product_impl(lhs, rhs, resCol, true); - res = resCol.markAsRValue(); - } - else - { - ColMajorMatrixAux resCol(lhs.rows(),rhs.cols()); - // ressort to transpose to sort the entries - internal::conservative_sparse_sparse_product_impl(lhs, rhs, resCol, false); - RowMajorMatrix resRow(resCol); - res = resRow.markAsRValue(); - } - } -}; - -template -struct conservative_sparse_sparse_product_selector -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix RowMajorRhs; - typedef SparseMatrix RowMajorRes; - RowMajorRhs rhsRow = rhs; - RowMajorRes 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 RowMajorLhs; - typedef SparseMatrix RowMajorRes; - RowMajorLhs lhsRow = lhs; - RowMajorRes 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 ColMajorLhs; - typedef SparseMatrix ColMajorRes; - ColMajorLhs lhsCol = lhs; - ColMajorRes 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 ColMajorRhs; - typedef SparseMatrix ColMajorRes; - ColMajorRhs rhsCol = rhs; - ColMajorRes 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 - - -namespace internal { - -template -static void sparse_sparse_to_dense_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res) -{ - typedef typename remove_all::type::Scalar LhsScalar; - typedef typename remove_all::type::Scalar RhsScalar; - Index cols = rhs.outerSize(); - eigen_assert(lhs.outerSize() == rhs.innerSize()); - - evaluator lhsEval(lhs); - evaluator rhsEval(rhs); - - for (Index j=0; j::InnerIterator rhsIt(rhsEval, j); rhsIt; ++rhsIt) - { - RhsScalar y = rhsIt.value(); - Index k = rhsIt.index(); - for (typename evaluator::InnerIterator lhsIt(lhsEval, k); lhsIt; ++lhsIt) - { - Index i = lhsIt.index(); - LhsScalar x = lhsIt.value(); - res.coeffRef(i,j) += x * y; - } - } - } -} - - -} // end namespace internal - -namespace internal { - -template::Flags&RowMajorBit) ? RowMajor : ColMajor, - int RhsStorageOrder = (traits::Flags&RowMajorBit) ? RowMajor : ColMajor> -struct sparse_sparse_to_dense_product_selector; - -template -struct sparse_sparse_to_dense_product_selector -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - internal::sparse_sparse_to_dense_product_impl(lhs, rhs, res); - } -}; - -template -struct sparse_sparse_to_dense_product_selector -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix ColMajorLhs; - ColMajorLhs lhsCol(lhs); - internal::sparse_sparse_to_dense_product_impl(lhsCol, rhs, res); - } -}; - -template -struct sparse_sparse_to_dense_product_selector -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - typedef SparseMatrix ColMajorRhs; - ColMajorRhs rhsCol(rhs); - internal::sparse_sparse_to_dense_product_impl(lhs, rhsCol, res); - } -}; - -template -struct sparse_sparse_to_dense_product_selector -{ - static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) - { - Transpose trRes(res); - internal::sparse_sparse_to_dense_product_impl >(rhs, lhs, trRes); - } -}; - - -} // 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 deleted file mode 100644 index 67718c8..0000000 --- a/eigen/Eigen/src/SparseCore/MappedSparseMatrix.h +++ /dev/null @@ -1,67 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2014 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 { - -/** \deprecated Use Map > - * \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 > -{}; -} // end namespace internal - -template -class MappedSparseMatrix - : public Map > -{ - typedef Map > Base; - - public: - - typedef typename Base::StorageIndex StorageIndex; - typedef typename Base::Scalar Scalar; - - inline MappedSparseMatrix(Index rows, Index cols, Index nnz, StorageIndex* outerIndexPtr, StorageIndex* innerIndexPtr, Scalar* valuePtr, StorageIndex* innerNonZeroPtr = 0) - : Base(rows, cols, nnz, outerIndexPtr, innerIndexPtr, valuePtr, innerNonZeroPtr) - {} - - /** Empty destructor */ - inline ~MappedSparseMatrix() {} -}; - -namespace internal { - -template -struct evaluator > - : evaluator > > -{ - typedef MappedSparseMatrix<_Scalar,_Options,_StorageIndex> XprType; - typedef evaluator > Base; - - evaluator() : Base() {} - explicit evaluator(const XprType &mat) : Base(mat) {} -}; - -} - -} // end namespace Eigen - -#endif // EIGEN_MAPPED_SPARSEMATRIX_H diff --git a/eigen/Eigen/src/SparseCore/SparseAssign.h b/eigen/Eigen/src/SparseCore/SparseAssign.h deleted file mode 100644 index 18352a8..0000000 --- a/eigen/Eigen/src/SparseCore/SparseAssign.h +++ /dev/null @@ -1,216 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2014 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_SPARSEASSIGN_H -#define EIGEN_SPARSEASSIGN_H - -namespace Eigen { - -template -template -Derived& SparseMatrixBase::operator=(const EigenBase &other) -{ - internal::call_assignment_no_alias(derived(), other.derived()); - return derived(); -} - -template -template -Derived& SparseMatrixBase::operator=(const ReturnByValue& other) -{ - // TODO use the evaluator mechanism - other.evalTo(derived()); - return derived(); -} - -template -template -inline Derived& SparseMatrixBase::operator=(const SparseMatrixBase& other) -{ - // by default sparse evaluation do not alias, so we can safely bypass the generic call_assignment routine - internal::Assignment > - ::run(derived(), other.derived(), internal::assign_op()); - return derived(); -} - -template -inline Derived& SparseMatrixBase::operator=(const Derived& other) -{ - internal::call_assignment_no_alias(derived(), other.derived()); - return derived(); -} - -namespace internal { - -template<> -struct storage_kind_to_evaluator_kind { - typedef IteratorBased Kind; -}; - -template<> -struct storage_kind_to_shape { - typedef SparseShape Shape; -}; - -struct Sparse2Sparse {}; -struct Sparse2Dense {}; - -template<> struct AssignmentKind { typedef Sparse2Sparse Kind; }; -template<> struct AssignmentKind { typedef Sparse2Sparse Kind; }; -template<> struct AssignmentKind { typedef Sparse2Dense Kind; }; -template<> struct AssignmentKind { typedef Sparse2Dense Kind; }; - - -template -void assign_sparse_to_sparse(DstXprType &dst, const SrcXprType &src) -{ - typedef typename DstXprType::Scalar Scalar; - typedef internal::evaluator DstEvaluatorType; - typedef internal::evaluator SrcEvaluatorType; - - SrcEvaluatorType srcEvaluator(src); - - const bool transpose = (DstEvaluatorType::Flags & RowMajorBit) != (SrcEvaluatorType::Flags & RowMajorBit); - const Index outerEvaluationSize = (SrcEvaluatorType::Flags&RowMajorBit) ? src.rows() : src.cols(); - if ((!transpose) && src.isRValue()) - { - // eval without temporary - dst.resize(src.rows(), src.cols()); - dst.setZero(); - dst.reserve((std::max)(src.rows(),src.cols())*2); - for (Index j=0; j::SupportedAccessPatterns & OuterRandomAccessPattern)==OuterRandomAccessPattern) || - (!((DstEvaluatorType::Flags & RowMajorBit) != (SrcEvaluatorType::Flags & RowMajorBit)))) && - "the transpose operation is supposed to be handled in SparseMatrix::operator="); - - enum { Flip = (DstEvaluatorType::Flags & RowMajorBit) != (SrcEvaluatorType::Flags & RowMajorBit) }; - - - DstXprType temp(src.rows(), src.cols()); - - temp.reserve((std::max)(src.rows(),src.cols())*2); - for (Index j=0; j -struct Assignment -{ - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) - { - assign_sparse_to_sparse(dst.derived(), src.derived()); - } -}; - -// Generic Sparse to Dense assignment -template< typename DstXprType, typename SrcXprType, typename Functor> -struct Assignment -{ - static void run(DstXprType &dst, const SrcXprType &src, const Functor &func) - { - if(internal::is_same >::value) - dst.setZero(); - - internal::evaluator srcEval(src); - resize_if_allowed(dst, src, func); - internal::evaluator dstEval(dst); - - const Index outerEvaluationSize = (internal::evaluator::Flags&RowMajorBit) ? src.rows() : src.cols(); - for (Index j=0; j::InnerIterator i(srcEval,j); i; ++i) - func.assignCoeff(dstEval.coeffRef(i.row(),i.col()), i.value()); - } -}; - -// Specialization for "dst = dec.solve(rhs)" -// NOTE we need to specialize it for Sparse2Sparse to avoid ambiguous specialization error -template -struct Assignment, internal::assign_op, Sparse2Sparse> -{ - typedef Solve SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) - { - Index dstRows = src.rows(); - Index dstCols = src.cols(); - if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) - dst.resize(dstRows, dstCols); - - src.dec()._solve_impl(src.rhs(), dst); - } -}; - -struct Diagonal2Sparse {}; - -template<> struct AssignmentKind { typedef Diagonal2Sparse Kind; }; - -template< typename DstXprType, typename SrcXprType, typename Functor> -struct Assignment -{ - typedef typename DstXprType::StorageIndex StorageIndex; - typedef typename DstXprType::Scalar Scalar; - typedef Array ArrayXI; - typedef Array ArrayXS; - template - static void run(SparseMatrix &dst, const SrcXprType &src, const internal::assign_op &/*func*/) - { - Index dstRows = src.rows(); - Index dstCols = src.cols(); - if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) - dst.resize(dstRows, dstCols); - - Index size = src.diagonal().size(); - dst.makeCompressed(); - dst.resizeNonZeros(size); - Map(dst.innerIndexPtr(), size).setLinSpaced(0,StorageIndex(size)-1); - Map(dst.outerIndexPtr(), size+1).setLinSpaced(0,StorageIndex(size)); - Map(dst.valuePtr(), size) = src.diagonal(); - } - - template - static void run(SparseMatrixBase &dst, const SrcXprType &src, const internal::assign_op &/*func*/) - { - dst.diagonal() = src.diagonal(); - } - - static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op &/*func*/) - { dst.diagonal() += src.diagonal(); } - - static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op &/*func*/) - { dst.diagonal() -= src.diagonal(); } -}; -} // end namespace internal - -} // end namespace Eigen - -#endif // EIGEN_SPARSEASSIGN_H diff --git a/eigen/Eigen/src/SparseCore/SparseBlock.h b/eigen/Eigen/src/SparseCore/SparseBlock.h deleted file mode 100644 index 511e92b..0000000 --- a/eigen/Eigen/src/SparseCore/SparseBlock.h +++ /dev/null @@ -1,603 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2014 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 { - -// Subset of columns or rows -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 }; - typedef SparseMatrixBase Base; - using Base::convert_index; -public: - EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) - - inline BlockImpl(XprType& xpr, Index i) - : m_matrix(xpr), m_outerStart(convert_index(i)), m_outerSize(OuterSize) - {} - - inline BlockImpl(XprType& xpr, Index startRow, Index startCol, Index blockRows, Index blockCols) - : m_matrix(xpr), m_outerStart(convert_index(IsRowMajor ? startRow : startCol)), m_outerSize(convert_index(IsRowMajor ? blockRows : blockCols)) - {} - - 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(); } - - Index nonZeros() const - { - typedef internal::evaluator EvaluatorType; - EvaluatorType matEval(m_matrix); - Index nnz = 0; - Index end = m_outerStart + m_outerSize.value(); - for(Index j=m_outerStart; j::non_const_type m_matrix; - Index m_outerStart; - const internal::variable_if_dynamic m_outerSize; - - protected: - // Disable assignment with clear error message. - // Note that simply removing operator= yields compilation errors with ICC+MSVC - template - BlockImpl& operator=(const T&) - { - EIGEN_STATIC_ASSERT(sizeof(T)==0, THIS_SPARSE_BLOCK_SUBEXPRESSION_IS_READ_ONLY); - return *this; - } -}; - - -/*************************************************************************** -* specialization for SparseMatrix -***************************************************************************/ - -namespace internal { - -template -class sparse_matrix_block_impl - : public SparseCompressedBase > -{ - typedef typename internal::remove_all::type _MatrixTypeNested; - typedef Block BlockType; - typedef SparseCompressedBase > Base; - using Base::convert_index; -public: - enum { IsRowMajor = internal::traits::IsRowMajor }; - EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) -protected: - typedef typename Base::IndexVector IndexVector; - enum { OuterSize = IsRowMajor ? BlockRows : BlockCols }; -public: - - inline sparse_matrix_block_impl(SparseMatrixType& xpr, Index i) - : m_matrix(xpr), m_outerStart(convert_index(i)), m_outerSize(OuterSize) - {} - - inline sparse_matrix_block_impl(SparseMatrixType& xpr, Index startRow, Index startCol, Index blockRows, Index blockCols) - : m_matrix(xpr), m_outerStart(convert_index(IsRowMajor ? startRow : startCol)), m_outerSize(convert_index(IsRowMajor ? blockRows : blockCols)) - {} - - template - inline BlockType& operator=(const SparseMatrixBase& other) - { - typedef typename internal::remove_all::type _NestedMatrixType; - _NestedMatrixType& matrix = m_matrix; - // This assignment 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 - Ref > tmp(other.derived()); - eigen_internal_assert(tmp.outerSize()==m_outerSize.value()); - - // 2 - let's check whether there is enough allocated memory - Index nnz = tmp.nonZeros(); - Index start = m_outerStart==0 ? 0 : m_matrix.outerIndexPtr()[m_outerStart]; // starting position of the current block - Index end = m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()]; // ending position 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; - - Index tmp_start = tmp.outerIndexPtr()[0]; - - bool update_trailing_pointers = false; - if(nnz>free_size) - { - // realloc manually to reduce copies - typename SparseMatrixType::Storage newdata(m_matrix.data().allocatedSize() - block_size + nnz); - - internal::smart_copy(m_matrix.valuePtr(), m_matrix.valuePtr() + start, newdata.valuePtr()); - internal::smart_copy(m_matrix.innerIndexPtr(), m_matrix.innerIndexPtr() + start, newdata.indexPtr()); - - internal::smart_copy(tmp.valuePtr() + tmp_start, tmp.valuePtr() + tmp_start + nnz, newdata.valuePtr() + start); - internal::smart_copy(tmp.innerIndexPtr() + tmp_start, tmp.innerIndexPtr() + tmp_start + nnz, newdata.indexPtr() + start); - - internal::smart_copy(matrix.valuePtr()+end, matrix.valuePtr()+end + tail_size, newdata.valuePtr()+start+nnz); - internal::smart_copy(matrix.innerIndexPtr()+end, matrix.innerIndexPtr()+end + tail_size, newdata.indexPtr()+start+nnz); - - newdata.resize(m_matrix.outerIndexPtr()[m_matrix.outerSize()] - block_size + nnz); - - matrix.data().swap(newdata); - - update_trailing_pointers = true; - } - else - { - if(m_matrix.isCompressed()) - { - // no need to realloc, simply copy the tail at its respective position and insert tmp - matrix.data().resize(start + nnz + tail_size); - - internal::smart_memmove(matrix.valuePtr()+end, matrix.valuePtr() + end+tail_size, matrix.valuePtr() + start+nnz); - internal::smart_memmove(matrix.innerIndexPtr()+end, matrix.innerIndexPtr() + end+tail_size, matrix.innerIndexPtr() + start+nnz); - - update_trailing_pointers = true; - } - - internal::smart_copy(tmp.valuePtr() + tmp_start, tmp.valuePtr() + tmp_start + nnz, matrix.valuePtr() + start); - internal::smart_copy(tmp.innerIndexPtr() + tmp_start, tmp.innerIndexPtr() + tmp_start + nnz, matrix.innerIndexPtr() + start); - } - - // update outer index pointers and innerNonZeros - if(IsVectorAtCompileTime) - { - if(!m_matrix.isCompressed()) - matrix.innerNonZeroPtr()[m_outerStart] = StorageIndex(nnz); - matrix.outerIndexPtr()[m_outerStart] = StorageIndex(start); - } - else - { - StorageIndex p = StorageIndex(start); - for(Index k=0; k(tmp.innerVector(k).nonZeros()); - if(!m_matrix.isCompressed()) - matrix.innerNonZeroPtr()[m_outerStart+k] = nnz_k; - matrix.outerIndexPtr()[m_outerStart+k] = p; - p += nnz_k; - } - } - - if(update_trailing_pointers) - { - StorageIndex offset = internal::convert_index(nnz - block_size); - for(Index k = m_outerStart + m_outerSize.value(); k<=matrix.outerSize(); ++k) - { - matrix.outerIndexPtr()[k] += offset; - } - } - - return derived(); - } - - inline BlockType& operator=(const BlockType& other) - { - return operator=(other); - } - - inline const Scalar* valuePtr() const - { return m_matrix.valuePtr(); } - inline Scalar* valuePtr() - { return m_matrix.valuePtr(); } - - inline const StorageIndex* innerIndexPtr() const - { return m_matrix.innerIndexPtr(); } - inline StorageIndex* innerIndexPtr() - { return m_matrix.innerIndexPtr(); } - - inline const StorageIndex* outerIndexPtr() const - { return m_matrix.outerIndexPtr() + m_outerStart; } - inline StorageIndex* outerIndexPtr() - { return m_matrix.outerIndexPtr() + m_outerStart; } - - inline const StorageIndex* innerNonZeroPtr() const - { return isCompressed() ? 0 : (m_matrix.innerNonZeroPtr()+m_outerStart); } - inline StorageIndex* innerNonZeroPtr() - { return isCompressed() ? 0 : (m_matrix.innerNonZeroPtr()+m_outerStart); } - - bool isCompressed() const { return m_matrix.innerNonZeroPtr()==0; } - - inline Scalar& coeffRef(Index row, Index col) - { - return m_matrix.coeffRef(row + (IsRowMajor ? m_outerStart : 0), col + (IsRowMajor ? 0 : m_outerStart)); - } - - inline const Scalar coeff(Index row, Index col) const - { - return m_matrix.coeff(row + (IsRowMajor ? m_outerStart : 0), col + (IsRowMajor ? 0 : m_outerStart)); - } - - inline const Scalar coeff(Index index) const - { - return m_matrix.coeff(IsRowMajor ? m_outerStart : index, IsRowMajor ? index : m_outerStart); - } - - const Scalar& lastCoeff() const - { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(sparse_matrix_block_impl); - eigen_assert(Base::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(); } - - inline const SparseMatrixType& nestedExpression() const { return m_matrix; } - inline SparseMatrixType& nestedExpression() { return m_matrix; } - Index startRow() const { return IsRowMajor ? m_outerStart : 0; } - Index startCol() const { return IsRowMajor ? 0 : m_outerStart; } - Index blockRows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); } - Index blockCols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); } - - protected: - - typename internal::ref_selector::non_const_type m_matrix; - Index m_outerStart; - const internal::variable_if_dynamic m_outerSize; - -}; - -} // namespace internal - -template -class BlockImpl,BlockRows,BlockCols,true,Sparse> - : public internal::sparse_matrix_block_impl,BlockRows,BlockCols> -{ -public: - typedef _StorageIndex StorageIndex; - typedef SparseMatrix<_Scalar, _Options, _StorageIndex> SparseMatrixType; - typedef internal::sparse_matrix_block_impl Base; - inline BlockImpl(SparseMatrixType& xpr, Index i) - : Base(xpr, i) - {} - - inline BlockImpl(SparseMatrixType& xpr, Index startRow, Index startCol, Index blockRows, Index blockCols) - : Base(xpr, startRow, startCol, blockRows, blockCols) - {} - - using Base::operator=; -}; - -template -class BlockImpl,BlockRows,BlockCols,true,Sparse> - : public internal::sparse_matrix_block_impl,BlockRows,BlockCols> -{ -public: - typedef _StorageIndex StorageIndex; - typedef const SparseMatrix<_Scalar, _Options, _StorageIndex> SparseMatrixType; - typedef internal::sparse_matrix_block_impl Base; - inline BlockImpl(SparseMatrixType& xpr, Index i) - : Base(xpr, i) - {} - - inline BlockImpl(SparseMatrixType& xpr, Index startRow, Index startCol, Index blockRows, Index blockCols) - : Base(xpr, startRow, startCol, blockRows, blockCols) - {} - - using Base::operator=; -private: - template BlockImpl(const SparseMatrixBase& xpr, Index i); - template BlockImpl(const SparseMatrixBase& xpr); -}; - -//---------- - -/** \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 Block BlockType; - typedef SparseMatrixBase Base; - using Base::convert_index; -public: - enum { IsRowMajor = internal::traits::IsRowMajor }; - EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) - - typedef typename internal::remove_all::type _MatrixTypeNested; - - /** Column or Row constructor - */ - inline BlockImpl(XprType& xpr, Index i) - : m_matrix(xpr), - m_startRow( (BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) ? convert_index(i) : 0), - m_startCol( (BlockRows==XprType::RowsAtCompileTime) && (BlockCols==1) ? convert_index(i) : 0), - m_blockRows(BlockRows==1 ? 1 : xpr.rows()), - m_blockCols(BlockCols==1 ? 1 : xpr.cols()) - {} - - /** Dynamic-size constructor - */ - inline BlockImpl(XprType& xpr, Index startRow, Index startCol, Index blockRows, Index blockCols) - : m_matrix(xpr), m_startRow(convert_index(startRow)), m_startCol(convert_index(startCol)), m_blockRows(convert_index(blockRows)), m_blockCols(convert_index(blockCols)) - {} - - inline Index rows() const { return m_blockRows.value(); } - inline Index cols() const { return m_blockCols.value(); } - - inline Scalar& coeffRef(Index row, Index col) - { - return m_matrix.coeffRef(row + m_startRow.value(), col + m_startCol.value()); - } - - inline const Scalar coeff(Index row, Index col) const - { - return m_matrix.coeff(row + m_startRow.value(), col + m_startCol.value()); - } - - inline Scalar& coeffRef(Index index) - { - return m_matrix.coeffRef(m_startRow.value() + (RowsAtCompileTime == 1 ? 0 : index), - m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0)); - } - - inline const Scalar coeff(Index index) const - { - return m_matrix.coeff(m_startRow.value() + (RowsAtCompileTime == 1 ? 0 : index), - m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0)); - } - - inline const XprType& nestedExpression() const { return m_matrix; } - inline XprType& nestedExpression() { return m_matrix; } - Index startRow() const { return m_startRow.value(); } - Index startCol() const { return m_startCol.value(); } - Index blockRows() const { return m_blockRows.value(); } - Index blockCols() const { return m_blockCols.value(); } - - protected: -// friend class internal::GenericSparseBlockInnerIteratorImpl; - friend struct internal::unary_evaluator, internal::IteratorBased, Scalar >; - - Index nonZeros() const { return Dynamic; } - - typename internal::ref_selector::non_const_type 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; - - protected: - // Disable assignment with clear error message. - // Note that simply removing operator= yields compilation errors with ICC+MSVC - template - BlockImpl& operator=(const T&) - { - EIGEN_STATIC_ASSERT(sizeof(T)==0, THIS_SPARSE_BLOCK_SUBEXPRESSION_IS_READ_ONLY); - return *this; - } - -}; - -namespace internal { - -template -struct unary_evaluator, IteratorBased > - : public evaluator_base > -{ - class InnerVectorInnerIterator; - class OuterVectorInnerIterator; - public: - typedef Block XprType; - typedef typename XprType::StorageIndex StorageIndex; - typedef typename XprType::Scalar Scalar; - - enum { - IsRowMajor = XprType::IsRowMajor, - - OuterVector = (BlockCols==1 && ArgType::IsRowMajor) - | // FIXME | instead of || to please GCC 4.4.0 stupid warning "suggest parentheses around &&". - // revert to || as soon as not needed anymore. - (BlockRows==1 && !ArgType::IsRowMajor), - - CoeffReadCost = evaluator::CoeffReadCost, - Flags = XprType::Flags - }; - - typedef typename internal::conditional::type InnerIterator; - - explicit unary_evaluator(const XprType& op) - : m_argImpl(op.nestedExpression()), m_block(op) - {} - - inline Index nonZerosEstimate() const { - Index nnz = m_block.nonZeros(); - if(nnz<0) - return m_argImpl.nonZerosEstimate() * m_block.size() / m_block.nestedExpression().size(); - return nnz; - } - - protected: - typedef typename evaluator::InnerIterator EvalIterator; - - evaluator m_argImpl; - const XprType &m_block; -}; - -template -class unary_evaluator, IteratorBased>::InnerVectorInnerIterator - : public EvalIterator -{ - enum { IsRowMajor = unary_evaluator::IsRowMajor }; - const XprType& m_block; - Index m_end; -public: - - EIGEN_STRONG_INLINE InnerVectorInnerIterator(const unary_evaluator& aEval, Index outer) - : EvalIterator(aEval.m_argImpl, outer + (IsRowMajor ? aEval.m_block.startRow() : aEval.m_block.startCol())), - m_block(aEval.m_block), - m_end(IsRowMajor ? aEval.m_block.startCol()+aEval.m_block.blockCols() : aEval.m_block.startRow()+aEval.m_block.blockRows()) - { - while( (EvalIterator::operator bool()) && (EvalIterator::index() < (IsRowMajor ? m_block.startCol() : m_block.startRow())) ) - EvalIterator::operator++(); - } - - inline StorageIndex index() const { return EvalIterator::index() - convert_index(IsRowMajor ? m_block.startCol() : m_block.startRow()); } - inline Index outer() const { return EvalIterator::outer() - (IsRowMajor ? m_block.startRow() : m_block.startCol()); } - inline Index row() const { return EvalIterator::row() - m_block.startRow(); } - inline Index col() const { return EvalIterator::col() - m_block.startCol(); } - - inline operator bool() const { return EvalIterator::operator bool() && EvalIterator::index() < m_end; } -}; - -template -class unary_evaluator, IteratorBased>::OuterVectorInnerIterator -{ - enum { IsRowMajor = unary_evaluator::IsRowMajor }; - const unary_evaluator& m_eval; - Index m_outerPos; - const Index m_innerIndex; - Index m_end; - EvalIterator m_it; -public: - - EIGEN_STRONG_INLINE OuterVectorInnerIterator(const unary_evaluator& aEval, Index outer) - : m_eval(aEval), - m_outerPos( (IsRowMajor ? aEval.m_block.startCol() : aEval.m_block.startRow()) ), - m_innerIndex(IsRowMajor ? aEval.m_block.startRow() : aEval.m_block.startCol()), - m_end(IsRowMajor ? aEval.m_block.startCol()+aEval.m_block.blockCols() : aEval.m_block.startRow()+aEval.m_block.blockRows()), - m_it(m_eval.m_argImpl, m_outerPos) - { - EIGEN_UNUSED_VARIABLE(outer); - eigen_assert(outer==0); - - while(m_it && m_it.index() < m_innerIndex) ++m_it; - if((!m_it) || (m_it.index()!=m_innerIndex)) - ++(*this); - } - - inline StorageIndex index() const { return convert_index(m_outerPos - (IsRowMajor ? m_eval.m_block.startCol() : m_eval.m_block.startRow())); } - inline Index outer() const { return 0; } - inline Index row() const { return IsRowMajor ? 0 : index(); } - inline Index col() const { return IsRowMajor ? index() : 0; } - - inline Scalar value() const { return m_it.value(); } - inline Scalar& valueRef() { return m_it.valueRef(); } - - inline OuterVectorInnerIterator& operator++() - { - // search next non-zero entry - while(++m_outerPos -struct unary_evaluator,BlockRows,BlockCols,true>, IteratorBased> - : evaluator,BlockRows,BlockCols,true> > > -{ - typedef Block,BlockRows,BlockCols,true> XprType; - typedef evaluator > Base; - explicit unary_evaluator(const XprType &xpr) : Base(xpr) {} -}; - -template -struct unary_evaluator,BlockRows,BlockCols,true>, IteratorBased> - : evaluator,BlockRows,BlockCols,true> > > -{ - typedef Block,BlockRows,BlockCols,true> XprType; - typedef evaluator > Base; - explicit unary_evaluator(const XprType &xpr) : Base(xpr) {} -}; - -} // end namespace internal - - -} // end namespace Eigen - -#endif // EIGEN_SPARSE_BLOCK_H diff --git a/eigen/Eigen/src/SparseCore/SparseColEtree.h b/eigen/Eigen/src/SparseCore/SparseColEtree.h deleted file mode 100644 index ebe02d1..0000000 --- a/eigen/Eigen/src/SparseCore/SparseColEtree.h +++ /dev/null @@ -1,206 +0,0 @@ -// 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::StorageIndex *perm=0) -{ - typedef typename MatrixType::StorageIndex StorageIndex; - StorageIndex nc = convert_index(mat.cols()); // Number of columns - StorageIndex m = convert_index(mat.rows()); - StorageIndex 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 - firstRowElt.resize(m); - firstRowElt.setConstant(nc); - firstRowElt.segment(0, diagSize).setLinSpaced(diagSize, 0, diagSize-1); - bool found_diag; - for (StorageIndex col = 0; col < nc; col++) - { - StorageIndex pcol = col; - if(perm) pcol = perm[col]; - for (typename MatrixType::InnerIterator it(mat, pcol); it; ++it) - { - Index 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. */ - StorageIndex rset, cset, rroot; - for (StorageIndex 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 */ - StorageIndex 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; - - StorageIndex 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 (typename IndexVector::Scalar n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, typename IndexVector::Scalar postnum) -{ - typedef typename IndexVector::Scalar StorageIndex; - StorageIndex 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(typename IndexVector::Scalar n, IndexVector& parent, IndexVector& post) -{ - typedef typename IndexVector::Scalar StorageIndex; - IndexVector first_kid, next_kid; // Linked list of children - StorageIndex 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 - first_kid.setConstant(-1); - for (StorageIndex v = n-1; v >= 0; v--) - { - StorageIndex 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/SparseCompressedBase.h b/eigen/Eigen/src/SparseCore/SparseCompressedBase.h deleted file mode 100644 index 5ccb466..0000000 --- a/eigen/Eigen/src/SparseCore/SparseCompressedBase.h +++ /dev/null @@ -1,341 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2015 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_COMPRESSED_BASE_H -#define EIGEN_SPARSE_COMPRESSED_BASE_H - -namespace Eigen { - -template class SparseCompressedBase; - -namespace internal { - -template -struct traits > : traits -{}; - -} // end namespace internal - -/** \ingroup SparseCore_Module - * \class SparseCompressedBase - * \brief Common base class for sparse [compressed]-{row|column}-storage format. - * - * This class defines the common interface for all derived classes implementing the compressed sparse storage format, such as: - * - SparseMatrix - * - Ref - * - Map - * - */ -template -class SparseCompressedBase - : public SparseMatrixBase -{ - public: - typedef SparseMatrixBase Base; - EIGEN_SPARSE_PUBLIC_INTERFACE(SparseCompressedBase) - using Base::operator=; - using Base::IsRowMajor; - - class InnerIterator; - class ReverseInnerIterator; - - protected: - typedef typename Base::IndexVector IndexVector; - Eigen::Map innerNonZeros() { return Eigen::Map(innerNonZeroPtr(), isCompressed()?0:derived().outerSize()); } - const Eigen::Map innerNonZeros() const { return Eigen::Map(innerNonZeroPtr(), isCompressed()?0:derived().outerSize()); } - - public: - - /** \returns the number of non zero coefficients */ - inline Index nonZeros() const - { - if(Derived::IsVectorAtCompileTime && outerIndexPtr()==0) - return derived().nonZeros(); - else if(isCompressed()) - return outerIndexPtr()[derived().outerSize()]-outerIndexPtr()[0]; - else if(derived().outerSize()==0) - return 0; - else - return innerNonZeros().sum(); - } - - /** \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 derived().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 derived().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 StorageIndex* innerIndexPtr() const { return derived().innerIndexPtr(); } - /** \returns a non-const pointer to the array of inner indices. - * This function is aimed at interoperability with other libraries. - * \sa valuePtr(), outerIndexPtr() */ - inline StorageIndex* innerIndexPtr() { return derived().innerIndexPtr(); } - - /** \returns a const pointer to the array of the starting positions of the inner vectors. - * This function is aimed at interoperability with other libraries. - * \warning it returns the null pointer 0 for SparseVector - * \sa valuePtr(), innerIndexPtr() */ - inline const StorageIndex* outerIndexPtr() const { return derived().outerIndexPtr(); } - /** \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. - * \warning it returns the null pointer 0 for SparseVector - * \sa valuePtr(), innerIndexPtr() */ - inline StorageIndex* outerIndexPtr() { return derived().outerIndexPtr(); } - - /** \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 StorageIndex* innerNonZeroPtr() const { return derived().innerNonZeroPtr(); } - /** \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 StorageIndex* innerNonZeroPtr() { return derived().innerNonZeroPtr(); } - - /** \returns whether \c *this is in compressed form. */ - inline bool isCompressed() const { return innerNonZeroPtr()==0; } - - /** \returns a read-only view of the stored coefficients as a 1D array expression. - * - * \warning this method is for \b compressed \b storage \b only, and it will trigger an assertion otherwise. - * - * \sa valuePtr(), isCompressed() */ - const Map > coeffs() const { eigen_assert(isCompressed()); return Array::Map(valuePtr(),nonZeros()); } - - /** \returns a read-write view of the stored coefficients as a 1D array expression - * - * \warning this method is for \b compressed \b storage \b only, and it will trigger an assertion otherwise. - * - * Here is an example: - * \include SparseMatrix_coeffs.cpp - * and the output is: - * \include SparseMatrix_coeffs.out - * - * \sa valuePtr(), isCompressed() */ - Map > coeffs() { eigen_assert(isCompressed()); return Array::Map(valuePtr(),nonZeros()); } - - protected: - /** Default constructor. Do nothing. */ - SparseCompressedBase() {} - private: - template explicit SparseCompressedBase(const SparseCompressedBase&); -}; - -template -class SparseCompressedBase::InnerIterator -{ - public: - InnerIterator() - : m_values(0), m_indices(0), m_outer(0), m_id(0), m_end(0) - {} - - InnerIterator(const InnerIterator& other) - : m_values(other.m_values), m_indices(other.m_indices), m_outer(other.m_outer), m_id(other.m_id), m_end(other.m_end) - {} - - InnerIterator& operator=(const InnerIterator& other) - { - m_values = other.m_values; - m_indices = other.m_indices; - const_cast(m_outer).setValue(other.m_outer.value()); - m_id = other.m_id; - m_end = other.m_end; - return *this; - } - - InnerIterator(const SparseCompressedBase& mat, Index outer) - : m_values(mat.valuePtr()), m_indices(mat.innerIndexPtr()), m_outer(outer) - { - if(Derived::IsVectorAtCompileTime && mat.outerIndexPtr()==0) - { - m_id = 0; - m_end = mat.nonZeros(); - } - else - { - m_id = mat.outerIndexPtr()[outer]; - if(mat.isCompressed()) - m_end = mat.outerIndexPtr()[outer+1]; - else - m_end = m_id + mat.innerNonZeroPtr()[outer]; - } - } - - explicit InnerIterator(const SparseCompressedBase& mat) - : m_values(mat.valuePtr()), m_indices(mat.innerIndexPtr()), m_outer(0), m_id(0), m_end(mat.nonZeros()) - { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived); - } - - explicit InnerIterator(const internal::CompressedStorage& data) - : m_values(data.valuePtr()), m_indices(data.indexPtr()), m_outer(0), m_id(0), m_end(data.size()) - { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived); - } - - 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 StorageIndex index() const { return m_indices[m_id]; } - inline Index outer() const { return m_outer.value(); } - inline Index row() const { return IsRowMajor ? m_outer.value() : index(); } - inline Index col() const { return IsRowMajor ? index() : m_outer.value(); } - - inline operator bool() const { return (m_id < m_end); } - - protected: - const Scalar* m_values; - const StorageIndex* m_indices; - typedef internal::variable_if_dynamic OuterType; - const OuterType m_outer; - Index m_id; - Index m_end; - private: - // If you get here, then you're not using the right InnerIterator type, e.g.: - // SparseMatrix A; - // SparseMatrix::InnerIterator it(A,0); - template InnerIterator(const SparseMatrixBase&, Index outer); -}; - -template -class SparseCompressedBase::ReverseInnerIterator -{ - public: - ReverseInnerIterator(const SparseCompressedBase& mat, Index outer) - : m_values(mat.valuePtr()), m_indices(mat.innerIndexPtr()), m_outer(outer) - { - if(Derived::IsVectorAtCompileTime && mat.outerIndexPtr()==0) - { - m_start = 0; - m_id = mat.nonZeros(); - } - else - { - m_start = mat.outerIndexPtr()[outer]; - if(mat.isCompressed()) - m_id = mat.outerIndexPtr()[outer+1]; - else - m_id = m_start + mat.innerNonZeroPtr()[outer]; - } - } - - explicit ReverseInnerIterator(const SparseCompressedBase& mat) - : m_values(mat.valuePtr()), m_indices(mat.innerIndexPtr()), m_outer(0), m_start(0), m_id(mat.nonZeros()) - { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived); - } - - explicit ReverseInnerIterator(const internal::CompressedStorage& data) - : m_values(data.valuePtr()), m_indices(data.indexPtr()), m_outer(0), m_start(0), m_id(data.size()) - { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived); - } - - 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 StorageIndex index() const { return m_indices[m_id-1]; } - inline Index outer() const { return m_outer.value(); } - inline Index row() const { return IsRowMajor ? m_outer.value() : index(); } - inline Index col() const { return IsRowMajor ? index() : m_outer.value(); } - - inline operator bool() const { return (m_id > m_start); } - - protected: - const Scalar* m_values; - const StorageIndex* m_indices; - typedef internal::variable_if_dynamic OuterType; - const OuterType m_outer; - Index m_start; - Index m_id; -}; - -namespace internal { - -template -struct evaluator > - : evaluator_base -{ - typedef typename Derived::Scalar Scalar; - typedef typename Derived::InnerIterator InnerIterator; - - enum { - CoeffReadCost = NumTraits::ReadCost, - Flags = Derived::Flags - }; - - evaluator() : m_matrix(0), m_zero(0) - { - EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); - } - explicit evaluator(const Derived &mat) : m_matrix(&mat), m_zero(0) - { - EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); - } - - inline Index nonZerosEstimate() const { - return m_matrix->nonZeros(); - } - - operator Derived&() { return m_matrix->const_cast_derived(); } - operator const Derived&() const { return *m_matrix; } - - typedef typename DenseCoeffsBase::CoeffReturnType CoeffReturnType; - const Scalar& coeff(Index row, Index col) const - { - Index p = find(row,col); - - if(p==Dynamic) - return m_zero; - else - return m_matrix->const_cast_derived().valuePtr()[p]; - } - - Scalar& coeffRef(Index row, Index col) - { - Index p = find(row,col); - eigen_assert(p!=Dynamic && "written coefficient does not exist"); - return m_matrix->const_cast_derived().valuePtr()[p]; - } - -protected: - - Index find(Index row, Index col) const - { - eigen_internal_assert(row>=0 && rowrows() && col>=0 && colcols()); - - const Index outer = Derived::IsRowMajor ? row : col; - const Index inner = Derived::IsRowMajor ? col : row; - - Index start = m_matrix->outerIndexPtr()[outer]; - Index end = m_matrix->isCompressed() ? m_matrix->outerIndexPtr()[outer+1] : m_matrix->outerIndexPtr()[outer] + m_matrix->innerNonZeroPtr()[outer]; - eigen_assert(end>=start && "you are using a non finalized sparse matrix or written coefficient does not exist"); - const Index p = std::lower_bound(m_matrix->innerIndexPtr()+start, m_matrix->innerIndexPtr()+end,inner) - m_matrix->innerIndexPtr(); - - return ((pinnerIndexPtr()[p]==inner)) ? p : Dynamic; - } - - const Derived *m_matrix; - const Scalar m_zero; -}; - -} - -} // end namespace Eigen - -#endif // EIGEN_SPARSE_COMPRESSED_BASE_H diff --git a/eigen/Eigen/src/SparseCore/SparseCwiseBinaryOp.h b/eigen/Eigen/src/SparseCore/SparseCwiseBinaryOp.h deleted file mode 100644 index e315e35..0000000 --- a/eigen/Eigen/src/SparseCore/SparseCwiseBinaryOp.h +++ /dev/null @@ -1,726 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2014 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 -// -// TODO to ease compiler job, we could specialize product/quotient with a scalar -// and fallback to cwise-unary evaluator using bind1st_op and bind2nd_op. - -template -class CwiseBinaryOpImpl - : public SparseMatrixBase > -{ - public: - typedef CwiseBinaryOp Derived; - typedef SparseMatrixBase Base; - EIGEN_SPARSE_PUBLIC_INTERFACE(Derived) - CwiseBinaryOpImpl() - { - EIGEN_STATIC_ASSERT(( - (!internal::is_same::StorageKind, - typename internal::traits::StorageKind>::value) - || ((internal::evaluator::Flags&RowMajorBit) == (internal::evaluator::Flags&RowMajorBit))), - THE_STORAGE_ORDER_OF_BOTH_SIDES_MUST_MATCH); - } -}; - -namespace internal { - - -// Generic "sparse OP sparse" -template struct binary_sparse_evaluator; - -template -struct binary_evaluator, IteratorBased, IteratorBased> - : evaluator_base > -{ -protected: - typedef typename evaluator::InnerIterator LhsIterator; - typedef typename evaluator::InnerIterator RhsIterator; - typedef CwiseBinaryOp XprType; - typedef typename traits::Scalar Scalar; - typedef typename XprType::StorageIndex StorageIndex; -public: - - class InnerIterator - { - public: - - EIGEN_STRONG_INLINE InnerIterator(const binary_evaluator& aEval, Index outer) - : m_lhsIter(aEval.m_lhsImpl,outer), m_rhsIter(aEval.m_rhsImpl,outer), m_functor(aEval.m_functor) - { - this->operator++(); - } - - EIGEN_STRONG_INLINE InnerIterator& 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 *this; - } - - EIGEN_STRONG_INLINE Scalar value() const { return m_value; } - - EIGEN_STRONG_INLINE StorageIndex index() const { return m_id; } - EIGEN_STRONG_INLINE Index outer() const { return m_lhsIter.outer(); } - 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; - StorageIndex m_id; - }; - - - enum { - CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, - Flags = XprType::Flags - }; - - explicit binary_evaluator(const XprType& xpr) - : m_functor(xpr.functor()), - m_lhsImpl(xpr.lhs()), - m_rhsImpl(xpr.rhs()) - { - EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); - EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); - } - - inline Index nonZerosEstimate() const { - return m_lhsImpl.nonZerosEstimate() + m_rhsImpl.nonZerosEstimate(); - } - -protected: - const BinaryOp m_functor; - evaluator m_lhsImpl; - evaluator m_rhsImpl; -}; - -// dense op sparse -template -struct binary_evaluator, IndexBased, IteratorBased> - : evaluator_base > -{ -protected: - typedef typename evaluator::InnerIterator RhsIterator; - typedef CwiseBinaryOp XprType; - typedef typename traits::Scalar Scalar; - typedef typename XprType::StorageIndex StorageIndex; -public: - - class InnerIterator - { - enum { IsRowMajor = (int(Rhs::Flags)&RowMajorBit)==RowMajorBit }; - public: - - EIGEN_STRONG_INLINE InnerIterator(const binary_evaluator& aEval, Index outer) - : m_lhsEval(aEval.m_lhsImpl), m_rhsIter(aEval.m_rhsImpl,outer), m_functor(aEval.m_functor), m_value(0), m_id(-1), m_innerSize(aEval.m_expr.rhs().innerSize()) - { - this->operator++(); - } - - EIGEN_STRONG_INLINE InnerIterator& operator++() - { - ++m_id; - if(m_id &m_lhsEval; - RhsIterator m_rhsIter; - const BinaryOp& m_functor; - Scalar m_value; - StorageIndex m_id; - StorageIndex m_innerSize; - }; - - - enum { - CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, - // Expose storage order of the sparse expression - Flags = (XprType::Flags & ~RowMajorBit) | (int(Rhs::Flags)&RowMajorBit) - }; - - explicit binary_evaluator(const XprType& xpr) - : m_functor(xpr.functor()), - m_lhsImpl(xpr.lhs()), - m_rhsImpl(xpr.rhs()), - m_expr(xpr) - { - EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); - EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); - } - - inline Index nonZerosEstimate() const { - return m_expr.size(); - } - -protected: - const BinaryOp m_functor; - evaluator m_lhsImpl; - evaluator m_rhsImpl; - const XprType &m_expr; -}; - -// sparse op dense -template -struct binary_evaluator, IteratorBased, IndexBased> - : evaluator_base > -{ -protected: - typedef typename evaluator::InnerIterator LhsIterator; - typedef CwiseBinaryOp XprType; - typedef typename traits::Scalar Scalar; - typedef typename XprType::StorageIndex StorageIndex; -public: - - class InnerIterator - { - enum { IsRowMajor = (int(Lhs::Flags)&RowMajorBit)==RowMajorBit }; - public: - - EIGEN_STRONG_INLINE InnerIterator(const binary_evaluator& aEval, Index outer) - : m_lhsIter(aEval.m_lhsImpl,outer), m_rhsEval(aEval.m_rhsImpl), m_functor(aEval.m_functor), m_value(0), m_id(-1), m_innerSize(aEval.m_expr.lhs().innerSize()) - { - this->operator++(); - } - - EIGEN_STRONG_INLINE InnerIterator& operator++() - { - ++m_id; - if(m_id &m_rhsEval; - const BinaryOp& m_functor; - Scalar m_value; - StorageIndex m_id; - StorageIndex m_innerSize; - }; - - - enum { - CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, - // Expose storage order of the sparse expression - Flags = (XprType::Flags & ~RowMajorBit) | (int(Lhs::Flags)&RowMajorBit) - }; - - explicit binary_evaluator(const XprType& xpr) - : m_functor(xpr.functor()), - m_lhsImpl(xpr.lhs()), - m_rhsImpl(xpr.rhs()), - m_expr(xpr) - { - EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); - EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); - } - - inline Index nonZerosEstimate() const { - return m_expr.size(); - } - -protected: - const BinaryOp m_functor; - evaluator m_lhsImpl; - evaluator m_rhsImpl; - const XprType &m_expr; -}; - -template::Kind, - typename RhsKind = typename evaluator_traits::Kind, - typename LhsScalar = typename traits::Scalar, - typename RhsScalar = typename traits::Scalar> struct sparse_conjunction_evaluator; - -// "sparse .* sparse" -template -struct binary_evaluator, Lhs, Rhs>, IteratorBased, IteratorBased> - : sparse_conjunction_evaluator, Lhs, Rhs> > -{ - typedef CwiseBinaryOp, Lhs, Rhs> XprType; - typedef sparse_conjunction_evaluator Base; - explicit binary_evaluator(const XprType& xpr) : Base(xpr) {} -}; -// "dense .* sparse" -template -struct binary_evaluator, Lhs, Rhs>, IndexBased, IteratorBased> - : sparse_conjunction_evaluator, Lhs, Rhs> > -{ - typedef CwiseBinaryOp, Lhs, Rhs> XprType; - typedef sparse_conjunction_evaluator Base; - explicit binary_evaluator(const XprType& xpr) : Base(xpr) {} -}; -// "sparse .* dense" -template -struct binary_evaluator, Lhs, Rhs>, IteratorBased, IndexBased> - : sparse_conjunction_evaluator, Lhs, Rhs> > -{ - typedef CwiseBinaryOp, Lhs, Rhs> XprType; - typedef sparse_conjunction_evaluator Base; - explicit binary_evaluator(const XprType& xpr) : Base(xpr) {} -}; - -// "sparse ./ dense" -template -struct binary_evaluator, Lhs, Rhs>, IteratorBased, IndexBased> - : sparse_conjunction_evaluator, Lhs, Rhs> > -{ - typedef CwiseBinaryOp, Lhs, Rhs> XprType; - typedef sparse_conjunction_evaluator Base; - explicit binary_evaluator(const XprType& xpr) : Base(xpr) {} -}; - -// "sparse && sparse" -template -struct binary_evaluator, IteratorBased, IteratorBased> - : sparse_conjunction_evaluator > -{ - typedef CwiseBinaryOp XprType; - typedef sparse_conjunction_evaluator Base; - explicit binary_evaluator(const XprType& xpr) : Base(xpr) {} -}; -// "dense && sparse" -template -struct binary_evaluator, IndexBased, IteratorBased> - : sparse_conjunction_evaluator > -{ - typedef CwiseBinaryOp XprType; - typedef sparse_conjunction_evaluator Base; - explicit binary_evaluator(const XprType& xpr) : Base(xpr) {} -}; -// "sparse && dense" -template -struct binary_evaluator, IteratorBased, IndexBased> - : sparse_conjunction_evaluator > -{ - typedef CwiseBinaryOp XprType; - typedef sparse_conjunction_evaluator Base; - explicit binary_evaluator(const XprType& xpr) : Base(xpr) {} -}; - -// "sparse ^ sparse" -template -struct sparse_conjunction_evaluator - : evaluator_base -{ -protected: - typedef typename XprType::Functor BinaryOp; - typedef typename XprType::Lhs LhsArg; - typedef typename XprType::Rhs RhsArg; - typedef typename evaluator::InnerIterator LhsIterator; - typedef typename evaluator::InnerIterator RhsIterator; - typedef typename XprType::StorageIndex StorageIndex; - typedef typename traits::Scalar Scalar; -public: - - class InnerIterator - { - public: - - EIGEN_STRONG_INLINE InnerIterator(const sparse_conjunction_evaluator& aEval, Index outer) - : m_lhsIter(aEval.m_lhsImpl,outer), m_rhsIter(aEval.m_rhsImpl,outer), m_functor(aEval.m_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 InnerIterator& 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 *this; - } - - EIGEN_STRONG_INLINE Scalar value() const { return m_functor(m_lhsIter.value(), m_rhsIter.value()); } - - EIGEN_STRONG_INLINE StorageIndex index() const { return m_lhsIter.index(); } - EIGEN_STRONG_INLINE Index outer() const { return m_lhsIter.outer(); } - 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 BinaryOp& m_functor; - }; - - - enum { - CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, - Flags = XprType::Flags - }; - - explicit sparse_conjunction_evaluator(const XprType& xpr) - : m_functor(xpr.functor()), - m_lhsImpl(xpr.lhs()), - m_rhsImpl(xpr.rhs()) - { - EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); - EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); - } - - inline Index nonZerosEstimate() const { - return (std::min)(m_lhsImpl.nonZerosEstimate(), m_rhsImpl.nonZerosEstimate()); - } - -protected: - const BinaryOp m_functor; - evaluator m_lhsImpl; - evaluator m_rhsImpl; -}; - -// "dense ^ sparse" -template -struct sparse_conjunction_evaluator - : evaluator_base -{ -protected: - typedef typename XprType::Functor BinaryOp; - typedef typename XprType::Lhs LhsArg; - typedef typename XprType::Rhs RhsArg; - typedef evaluator LhsEvaluator; - typedef typename evaluator::InnerIterator RhsIterator; - typedef typename XprType::StorageIndex StorageIndex; - typedef typename traits::Scalar Scalar; -public: - - class InnerIterator - { - enum { IsRowMajor = (int(RhsArg::Flags)&RowMajorBit)==RowMajorBit }; - - public: - - EIGEN_STRONG_INLINE InnerIterator(const sparse_conjunction_evaluator& aEval, Index outer) - : m_lhsEval(aEval.m_lhsImpl), m_rhsIter(aEval.m_rhsImpl,outer), m_functor(aEval.m_functor), m_outer(outer) - {} - - EIGEN_STRONG_INLINE InnerIterator& operator++() - { - ++m_rhsIter; - return *this; - } - - EIGEN_STRONG_INLINE Scalar value() const - { return m_functor(m_lhsEval.coeff(IsRowMajor?m_outer:m_rhsIter.index(),IsRowMajor?m_rhsIter.index():m_outer), m_rhsIter.value()); } - - EIGEN_STRONG_INLINE StorageIndex index() const { return m_rhsIter.index(); } - EIGEN_STRONG_INLINE Index outer() const { return m_rhsIter.outer(); } - 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 LhsEvaluator &m_lhsEval; - RhsIterator m_rhsIter; - const BinaryOp& m_functor; - const Index m_outer; - }; - - - enum { - CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, - // Expose storage order of the sparse expression - Flags = (XprType::Flags & ~RowMajorBit) | (int(RhsArg::Flags)&RowMajorBit) - }; - - explicit sparse_conjunction_evaluator(const XprType& xpr) - : m_functor(xpr.functor()), - m_lhsImpl(xpr.lhs()), - m_rhsImpl(xpr.rhs()) - { - EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); - EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); - } - - inline Index nonZerosEstimate() const { - return m_rhsImpl.nonZerosEstimate(); - } - -protected: - const BinaryOp m_functor; - evaluator m_lhsImpl; - evaluator m_rhsImpl; -}; - -// "sparse ^ dense" -template -struct sparse_conjunction_evaluator - : evaluator_base -{ -protected: - typedef typename XprType::Functor BinaryOp; - typedef typename XprType::Lhs LhsArg; - typedef typename XprType::Rhs RhsArg; - typedef typename evaluator::InnerIterator LhsIterator; - typedef evaluator RhsEvaluator; - typedef typename XprType::StorageIndex StorageIndex; - typedef typename traits::Scalar Scalar; -public: - - class InnerIterator - { - enum { IsRowMajor = (int(LhsArg::Flags)&RowMajorBit)==RowMajorBit }; - - public: - - EIGEN_STRONG_INLINE InnerIterator(const sparse_conjunction_evaluator& aEval, Index outer) - : m_lhsIter(aEval.m_lhsImpl,outer), m_rhsEval(aEval.m_rhsImpl), m_functor(aEval.m_functor), m_outer(outer) - {} - - EIGEN_STRONG_INLINE InnerIterator& operator++() - { - ++m_lhsIter; - return *this; - } - - EIGEN_STRONG_INLINE Scalar value() const - { return m_functor(m_lhsIter.value(), - m_rhsEval.coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); } - - EIGEN_STRONG_INLINE StorageIndex index() const { return m_lhsIter.index(); } - EIGEN_STRONG_INLINE Index outer() const { return m_lhsIter.outer(); } - 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: - LhsIterator m_lhsIter; - const evaluator &m_rhsEval; - const BinaryOp& m_functor; - const Index m_outer; - }; - - - enum { - CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, - // Expose storage order of the sparse expression - Flags = (XprType::Flags & ~RowMajorBit) | (int(LhsArg::Flags)&RowMajorBit) - }; - - explicit sparse_conjunction_evaluator(const XprType& xpr) - : m_functor(xpr.functor()), - m_lhsImpl(xpr.lhs()), - m_rhsImpl(xpr.rhs()) - { - EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); - EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); - } - - inline Index nonZerosEstimate() const { - return m_lhsImpl.nonZerosEstimate(); - } - -protected: - const BinaryOp m_functor; - evaluator m_lhsImpl; - evaluator m_rhsImpl; -}; - -} - -/*************************************************************************** -* Implementation of SparseMatrixBase and SparseCwise functions/operators -***************************************************************************/ - -template -template -Derived& SparseMatrixBase::operator+=(const EigenBase &other) -{ - call_assignment(derived(), other.derived(), internal::add_assign_op()); - return derived(); -} - -template -template -Derived& SparseMatrixBase::operator-=(const EigenBase &other) -{ - call_assignment(derived(), other.derived(), internal::assign_op()); - return derived(); -} - -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 -Derived& SparseMatrixBase::operator+=(const DiagonalBase& other) -{ - call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op()); - return derived(); -} - -template -template -Derived& SparseMatrixBase::operator-=(const DiagonalBase& other) -{ - call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op()); - return 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()); -} - -template -EIGEN_STRONG_INLINE const CwiseBinaryOp, const DenseDerived, const SparseDerived> -operator+(const MatrixBase &a, const SparseMatrixBase &b) -{ - return CwiseBinaryOp, const DenseDerived, const SparseDerived>(a.derived(), b.derived()); -} - -template -EIGEN_STRONG_INLINE const CwiseBinaryOp, const SparseDerived, const DenseDerived> -operator+(const SparseMatrixBase &a, const MatrixBase &b) -{ - return CwiseBinaryOp, const SparseDerived, const DenseDerived>(a.derived(), b.derived()); -} - -template -EIGEN_STRONG_INLINE const CwiseBinaryOp, const DenseDerived, const SparseDerived> -operator-(const MatrixBase &a, const SparseMatrixBase &b) -{ - return CwiseBinaryOp, const DenseDerived, const SparseDerived>(a.derived(), b.derived()); -} - -template -EIGEN_STRONG_INLINE const CwiseBinaryOp, const SparseDerived, const DenseDerived> -operator-(const SparseMatrixBase &a, const MatrixBase &b) -{ - return CwiseBinaryOp, const SparseDerived, const DenseDerived>(a.derived(), b.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 deleted file mode 100644 index ea79737..0000000 --- a/eigen/Eigen/src/SparseCore/SparseCwiseUnaryOp.h +++ /dev/null @@ -1,148 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2015 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 { - -namespace internal { - -template -struct unary_evaluator, IteratorBased> - : public evaluator_base > -{ - public: - typedef CwiseUnaryOp XprType; - - class InnerIterator; - - enum { - CoeffReadCost = evaluator::CoeffReadCost + functor_traits::Cost, - Flags = XprType::Flags - }; - - explicit unary_evaluator(const XprType& op) : m_functor(op.functor()), m_argImpl(op.nestedExpression()) - { - EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); - EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); - } - - inline Index nonZerosEstimate() const { - return m_argImpl.nonZerosEstimate(); - } - - protected: - typedef typename evaluator::InnerIterator EvalIterator; - - const UnaryOp m_functor; - evaluator m_argImpl; -}; - -template -class unary_evaluator, IteratorBased>::InnerIterator - : public unary_evaluator, IteratorBased>::EvalIterator -{ - typedef typename XprType::Scalar Scalar; - typedef typename unary_evaluator, IteratorBased>::EvalIterator Base; - public: - - EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& unaryOp, Index outer) - : Base(unaryOp.m_argImpl,outer), m_functor(unaryOp.m_functor) - {} - - EIGEN_STRONG_INLINE InnerIterator& operator++() - { Base::operator++(); return *this; } - - EIGEN_STRONG_INLINE Scalar value() const { return m_functor(Base::value()); } - - protected: - const UnaryOp m_functor; - private: - Scalar& valueRef(); -}; - -template -struct unary_evaluator, IteratorBased> - : public evaluator_base > -{ - public: - typedef CwiseUnaryView XprType; - - class InnerIterator; - - enum { - CoeffReadCost = evaluator::CoeffReadCost + functor_traits::Cost, - Flags = XprType::Flags - }; - - explicit unary_evaluator(const XprType& op) : m_functor(op.functor()), m_argImpl(op.nestedExpression()) - { - EIGEN_INTERNAL_CHECK_COST_VALUE(functor_traits::Cost); - EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); - } - - protected: - typedef typename evaluator::InnerIterator EvalIterator; - - const ViewOp m_functor; - evaluator m_argImpl; -}; - -template -class unary_evaluator, IteratorBased>::InnerIterator - : public unary_evaluator, IteratorBased>::EvalIterator -{ - typedef typename XprType::Scalar Scalar; - typedef typename unary_evaluator, IteratorBased>::EvalIterator Base; - public: - - EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& unaryOp, Index outer) - : Base(unaryOp.m_argImpl,outer), m_functor(unaryOp.m_functor) - {} - - EIGEN_STRONG_INLINE InnerIterator& operator++() - { Base::operator++(); return *this; } - - EIGEN_STRONG_INLINE Scalar value() const { return m_functor(Base::value()); } - EIGEN_STRONG_INLINE Scalar& valueRef() { return m_functor(Base::valueRef()); } - - protected: - const ViewOp m_functor; -}; - -} // end namespace internal - -template -EIGEN_STRONG_INLINE Derived& -SparseMatrixBase::operator*=(const Scalar& other) -{ - typedef typename internal::evaluator::InnerIterator EvalIterator; - internal::evaluator thisEval(derived()); - for (Index j=0; j -EIGEN_STRONG_INLINE Derived& -SparseMatrixBase::operator/=(const Scalar& other) -{ - typedef typename internal::evaluator::InnerIterator EvalIterator; - internal::evaluator thisEval(derived()); - 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 { - -namespace internal { - -template <> struct product_promote_storage_type { typedef Sparse ret; }; -template <> struct product_promote_storage_type { typedef Sparse ret; }; - -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 evaluator::InnerIterator LhsInnerIterator; - typedef evaluator LhsEval; - static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha) - { - LhsEval lhsEval(lhs); - - Index n = lhs.outerSize(); -#ifdef EIGEN_HAS_OPENMP - Eigen::initParallel(); - Index threads = Eigen::nbThreads(); -#endif - - for(Index c=0; c1 && lhsEval.nonZerosEstimate() > 20000) - { - #pragma omp parallel for schedule(dynamic,(n+threads*4-1)/(threads*4)) num_threads(threads) - for(Index i=0; i let's disable it for now as it is conflicting with generic scalar*matrix and matrix*scalar operators -// template -// struct ScalarBinaryOpTraits > -// { -// enum { -// Defined = 1 -// }; -// typedef typename CwiseUnaryOp, T2>::PlainObject ReturnType; -// }; - -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 evaluator::InnerIterator LhsInnerIterator; - static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha) - { - evaluator lhsEval(lhs); - for(Index c=0; c::ReturnType rhs_j(alpha * rhs.coeff(j,c)); - for(LhsInnerIterator it(lhsEval,j); it ;++it) - res.coeffRef(it.index(),c) += it.value() * rhs_j; - } - } - } -}; - -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 evaluator::InnerIterator LhsInnerIterator; - static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha) - { - evaluator lhsEval(lhs); - 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 evaluator::InnerIterator LhsInnerIterator; - static void run(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const typename Res::Scalar& alpha) - { - evaluator lhsEval(lhs); - 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 - -namespace internal { - -template -struct generic_product_impl - : generic_product_impl_base > -{ - typedef typename Product::Scalar Scalar; - - template - static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha) - { - typedef typename nested_eval::type LhsNested; - typedef typename nested_eval::type RhsNested; - LhsNested lhsNested(lhs); - RhsNested rhsNested(rhs); - internal::sparse_time_dense_product(lhsNested, rhsNested, dst, alpha); - } -}; - -template -struct generic_product_impl - : generic_product_impl -{}; - -template -struct generic_product_impl - : generic_product_impl_base > -{ - typedef typename Product::Scalar Scalar; - - template - static void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha) - { - typedef typename nested_eval::type LhsNested; - typedef typename nested_eval::type RhsNested; - LhsNested lhsNested(lhs); - RhsNested rhsNested(rhs); - - // transpose everything - Transpose dstT(dst); - internal::sparse_time_dense_product(rhsNested.transpose(), lhsNested.transpose(), dstT, alpha); - } -}; - -template -struct generic_product_impl - : generic_product_impl -{}; - -template -struct sparse_dense_outer_product_evaluator -{ -protected: - typedef typename conditional::type Lhs1; - typedef typename conditional::type ActualRhs; - typedef Product ProdXprType; - - // if the actual left-hand side is a dense vector, - // then build a sparse-view so that we can seamlessly iterate over it. - typedef typename conditional::StorageKind,Sparse>::value, - Lhs1, SparseView >::type ActualLhs; - typedef typename conditional::StorageKind,Sparse>::value, - Lhs1 const&, SparseView >::type LhsArg; - - typedef evaluator LhsEval; - typedef evaluator RhsEval; - typedef typename evaluator::InnerIterator LhsIterator; - typedef typename ProdXprType::Scalar Scalar; - -public: - enum { - Flags = NeedToTranspose ? RowMajorBit : 0, - CoeffReadCost = HugeCost - }; - - class InnerIterator : public LhsIterator - { - public: - InnerIterator(const sparse_dense_outer_product_evaluator &xprEval, Index outer) - : LhsIterator(xprEval.m_lhsXprImpl, 0), - m_outer(outer), - m_empty(false), - m_factor(get(xprEval.m_rhsXprImpl, outer, typename internal::traits::StorageKind() )) - {} - - EIGEN_STRONG_INLINE Index outer() const { return m_outer; } - EIGEN_STRONG_INLINE Index row() const { return NeedToTranspose ? m_outer : LhsIterator::index(); } - EIGEN_STRONG_INLINE Index col() const { return NeedToTranspose ? LhsIterator::index() : m_outer; } - - EIGEN_STRONG_INLINE Scalar value() const { return LhsIterator::value() * m_factor; } - EIGEN_STRONG_INLINE operator bool() const { return LhsIterator::operator bool() && (!m_empty); } - - protected: - Scalar get(const RhsEval &rhs, Index outer, Dense = Dense()) const - { - return rhs.coeff(outer); - } - - Scalar get(const RhsEval &rhs, Index outer, Sparse = Sparse()) - { - typename RhsEval::InnerIterator it(rhs, outer); - if (it && it.index()==0 && it.value()!=Scalar(0)) - return it.value(); - m_empty = true; - return Scalar(0); - } - - Index m_outer; - bool m_empty; - Scalar m_factor; - }; - - sparse_dense_outer_product_evaluator(const Lhs1 &lhs, const ActualRhs &rhs) - : m_lhs(lhs), m_lhsXprImpl(m_lhs), m_rhsXprImpl(rhs) - { - EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); - } - - // transpose case - sparse_dense_outer_product_evaluator(const ActualRhs &rhs, const Lhs1 &lhs) - : m_lhs(lhs), m_lhsXprImpl(m_lhs), m_rhsXprImpl(rhs) - { - EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); - } - -protected: - const LhsArg m_lhs; - evaluator m_lhsXprImpl; - evaluator m_rhsXprImpl; -}; - -// sparse * dense outer product -template -struct product_evaluator, OuterProduct, SparseShape, DenseShape> - : sparse_dense_outer_product_evaluator -{ - typedef sparse_dense_outer_product_evaluator Base; - - typedef Product XprType; - typedef typename XprType::PlainObject PlainObject; - - explicit product_evaluator(const XprType& xpr) - : Base(xpr.lhs(), xpr.rhs()) - {} - -}; - -template -struct product_evaluator, OuterProduct, DenseShape, SparseShape> - : sparse_dense_outer_product_evaluator -{ - typedef sparse_dense_outer_product_evaluator Base; - - typedef Product XprType; - typedef typename XprType::PlainObject PlainObject; - - explicit product_evaluator(const XprType& xpr) - : Base(xpr.lhs(), xpr.rhs()) - {} - -}; - -} // end namespace internal - -} // end namespace Eigen - -#endif // EIGEN_SPARSEDENSEPRODUCT_H diff --git a/eigen/Eigen/src/SparseCore/SparseDiagonalProduct.h b/eigen/Eigen/src/SparseCore/SparseDiagonalProduct.h deleted file mode 100644 index 941c03b..0000000 --- a/eigen/Eigen/src/SparseCore/SparseDiagonalProduct.h +++ /dev/null @@ -1,138 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2009-2015 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 { - -enum { - SDP_AsScalarProduct, - SDP_AsCwiseProduct -}; - -template -struct sparse_diagonal_product_evaluator; - -template -struct product_evaluator, ProductTag, DiagonalShape, SparseShape> - : public sparse_diagonal_product_evaluator -{ - typedef Product XprType; - enum { CoeffReadCost = HugeCost, Flags = Rhs::Flags&RowMajorBit, Alignment = 0 }; // FIXME CoeffReadCost & Flags - - typedef sparse_diagonal_product_evaluator Base; - explicit product_evaluator(const XprType& xpr) : Base(xpr.rhs(), xpr.lhs().diagonal()) {} -}; - -template -struct product_evaluator, ProductTag, SparseShape, DiagonalShape> - : public sparse_diagonal_product_evaluator, Lhs::Flags&RowMajorBit?SDP_AsCwiseProduct:SDP_AsScalarProduct> -{ - typedef Product XprType; - enum { CoeffReadCost = HugeCost, Flags = Lhs::Flags&RowMajorBit, Alignment = 0 }; // FIXME CoeffReadCost & Flags - - typedef sparse_diagonal_product_evaluator, Lhs::Flags&RowMajorBit?SDP_AsCwiseProduct:SDP_AsScalarProduct> Base; - explicit product_evaluator(const XprType& xpr) : Base(xpr.lhs(), xpr.rhs().diagonal().transpose()) {} -}; - -template -struct sparse_diagonal_product_evaluator -{ -protected: - typedef typename evaluator::InnerIterator SparseXprInnerIterator; - typedef typename SparseXprType::Scalar Scalar; - -public: - class InnerIterator : public SparseXprInnerIterator - { - public: - InnerIterator(const sparse_diagonal_product_evaluator &xprEval, Index outer) - : SparseXprInnerIterator(xprEval.m_sparseXprImpl, outer), - m_coeff(xprEval.m_diagCoeffImpl.coeff(outer)) - {} - - EIGEN_STRONG_INLINE Scalar value() const { return m_coeff * SparseXprInnerIterator::value(); } - protected: - typename DiagonalCoeffType::Scalar m_coeff; - }; - - sparse_diagonal_product_evaluator(const SparseXprType &sparseXpr, const DiagonalCoeffType &diagCoeff) - : m_sparseXprImpl(sparseXpr), m_diagCoeffImpl(diagCoeff) - {} - - Index nonZerosEstimate() const { return m_sparseXprImpl.nonZerosEstimate(); } - -protected: - evaluator m_sparseXprImpl; - evaluator m_diagCoeffImpl; -}; - - -template -struct sparse_diagonal_product_evaluator -{ - typedef typename SparseXprType::Scalar Scalar; - typedef typename SparseXprType::StorageIndex StorageIndex; - - typedef typename nested_eval::type DiagCoeffNested; - - class InnerIterator - { - typedef typename evaluator::InnerIterator SparseXprIter; - public: - InnerIterator(const sparse_diagonal_product_evaluator &xprEval, Index outer) - : m_sparseIter(xprEval.m_sparseXprEval, outer), m_diagCoeffNested(xprEval.m_diagCoeffNested) - {} - - inline Scalar value() const { return m_sparseIter.value() * m_diagCoeffNested.coeff(index()); } - inline StorageIndex index() const { return m_sparseIter.index(); } - inline Index outer() const { return m_sparseIter.outer(); } - inline Index col() const { return SparseXprType::IsRowMajor ? m_sparseIter.index() : m_sparseIter.outer(); } - inline Index row() const { return SparseXprType::IsRowMajor ? m_sparseIter.outer() : m_sparseIter.index(); } - - EIGEN_STRONG_INLINE InnerIterator& operator++() { ++m_sparseIter; return *this; } - inline operator bool() const { return m_sparseIter; } - - protected: - SparseXprIter m_sparseIter; - DiagCoeffNested m_diagCoeffNested; - }; - - sparse_diagonal_product_evaluator(const SparseXprType &sparseXpr, const DiagCoeffType &diagCoeff) - : m_sparseXprEval(sparseXpr), m_diagCoeffNested(diagCoeff) - {} - - Index nonZerosEstimate() const { return m_sparseXprEval.nonZerosEstimate(); } - -protected: - evaluator m_sparseXprEval; - DiagCoeffNested m_diagCoeffNested; -}; - -} // end namespace internal - -} // 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 deleted file mode 100644 index 38bc4aa..0000000 --- a/eigen/Eigen/src/SparseCore/SparseDot.h +++ /dev/null @@ -1,98 +0,0 @@ -// 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"); - - internal::evaluator thisEval(derived()); - typename internal::evaluator::InnerIterator i(thisEval, 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()); - - internal::evaluator thisEval(derived()); - typename internal::evaluator::InnerIterator i(thisEval, 0); - - internal::evaluator otherEval(other.derived()); - typename internal::evaluator::InnerIterator j(otherEval, 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 deleted file mode 100644 index 7d47eb9..0000000 --- a/eigen/Eigen/src/SparseCore/SparseFuzzy.h +++ /dev/null @@ -1,29 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2014 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 - -namespace Eigen { - -template -template -bool SparseMatrixBase::isApprox(const SparseMatrixBase& other, const RealScalar &prec) const -{ - const typename internal::nested_eval::type actualA(derived()); - typename internal::conditional::type, - const PlainObject>::type actualB(other.derived()); - - return (actualA - actualB).squaredNorm() <= prec * prec * numext::mini(actualA.squaredNorm(), actualB.squaredNorm()); -} - -} // end namespace Eigen - -#endif // EIGEN_SPARSE_FUZZY_H diff --git a/eigen/Eigen/src/SparseCore/SparseMap.h b/eigen/Eigen/src/SparseCore/SparseMap.h deleted file mode 100644 index f99be33..0000000 --- a/eigen/Eigen/src/SparseCore/SparseMap.h +++ /dev/null @@ -1,305 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2015 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_MAP_H -#define EIGEN_SPARSE_MAP_H - -namespace Eigen { - -namespace internal { - -template -struct traits, Options, StrideType> > - : public traits > -{ - typedef SparseMatrix PlainObjectType; - typedef traits TraitsBase; - enum { - Flags = TraitsBase::Flags & (~NestByRefBit) - }; -}; - -template -struct traits, Options, StrideType> > - : public traits > -{ - typedef SparseMatrix PlainObjectType; - typedef traits TraitsBase; - enum { - Flags = TraitsBase::Flags & (~ (NestByRefBit | LvalueBit)) - }; -}; - -} // end namespace internal - -template::has_write_access ? WriteAccessors : ReadOnlyAccessors -> class SparseMapBase; - -/** \ingroup SparseCore_Module - * class SparseMapBase - * \brief Common base class for Map and Ref instance of sparse matrix and vector. - */ -template -class SparseMapBase - : public SparseCompressedBase -{ - public: - typedef SparseCompressedBase Base; - typedef typename Base::Scalar Scalar; - typedef typename Base::StorageIndex StorageIndex; - enum { IsRowMajor = Base::IsRowMajor }; - using Base::operator=; - protected: - - typedef typename internal::conditional< - bool(internal::is_lvalue::value), - Scalar *, const Scalar *>::type ScalarPointer; - typedef typename internal::conditional< - bool(internal::is_lvalue::value), - StorageIndex *, const StorageIndex *>::type IndexPointer; - - Index m_outerSize; - Index m_innerSize; - Array m_zero_nnz; - IndexPointer m_outerIndex; - IndexPointer m_innerIndices; - ScalarPointer m_values; - IndexPointer m_innerNonZeros; - - public: - - /** \copydoc SparseMatrixBase::rows() */ - inline Index rows() const { return IsRowMajor ? m_outerSize : m_innerSize; } - /** \copydoc SparseMatrixBase::cols() */ - inline Index cols() const { return IsRowMajor ? m_innerSize : m_outerSize; } - /** \copydoc SparseMatrixBase::innerSize() */ - inline Index innerSize() const { return m_innerSize; } - /** \copydoc SparseMatrixBase::outerSize() */ - inline Index outerSize() const { return m_outerSize; } - /** \copydoc SparseCompressedBase::nonZeros */ - inline Index nonZeros() const { return m_zero_nnz[1]; } - - /** \copydoc SparseCompressedBase::isCompressed */ - bool isCompressed() const { return m_innerNonZeros==0; } - - //---------------------------------------- - // direct access interface - /** \copydoc SparseMatrix::valuePtr */ - inline const Scalar* valuePtr() const { return m_values; } - /** \copydoc SparseMatrix::innerIndexPtr */ - inline const StorageIndex* innerIndexPtr() const { return m_innerIndices; } - /** \copydoc SparseMatrix::outerIndexPtr */ - inline const StorageIndex* outerIndexPtr() const { return m_outerIndex; } - /** \copydoc SparseMatrix::innerNonZeroPtr */ - inline const StorageIndex* innerNonZeroPtr() const { return m_innerNonZeros; } - //---------------------------------------- - - /** \copydoc SparseMatrix::coeff */ - 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 = isCompressed() ? m_outerIndex[outer+1] : start + m_innerNonZeros[outer]; - 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 StorageIndex* r = std::lower_bound(&m_innerIndices[start],&m_innerIndices[end-1],inner); - const Index id = r-&m_innerIndices[0]; - return ((*r==inner) && (id(nnz)), m_outerIndex(outerIndexPtr), - m_innerIndices(innerIndexPtr), m_values(valuePtr), m_innerNonZeros(innerNonZerosPtr) - {} - - // for vectors - inline SparseMapBase(Index size, Index nnz, IndexPointer innerIndexPtr, ScalarPointer valuePtr) - : m_outerSize(1), m_innerSize(size), m_zero_nnz(0,internal::convert_index(nnz)), m_outerIndex(m_zero_nnz.data()), - m_innerIndices(innerIndexPtr), m_values(valuePtr), m_innerNonZeros(0) - {} - - /** Empty destructor */ - inline ~SparseMapBase() {} - - protected: - inline SparseMapBase() {} -}; - -/** \ingroup SparseCore_Module - * class SparseMapBase - * \brief Common base class for writable Map and Ref instance of sparse matrix and vector. - */ -template -class SparseMapBase - : public SparseMapBase -{ - typedef MapBase ReadOnlyMapBase; - - public: - typedef SparseMapBase Base; - typedef typename Base::Scalar Scalar; - typedef typename Base::StorageIndex StorageIndex; - enum { IsRowMajor = Base::IsRowMajor }; - - using Base::operator=; - - public: - - //---------------------------------------- - // direct access interface - using Base::valuePtr; - using Base::innerIndexPtr; - using Base::outerIndexPtr; - using Base::innerNonZeroPtr; - /** \copydoc SparseMatrix::valuePtr */ - inline Scalar* valuePtr() { return Base::m_values; } - /** \copydoc SparseMatrix::innerIndexPtr */ - inline StorageIndex* innerIndexPtr() { return Base::m_innerIndices; } - /** \copydoc SparseMatrix::outerIndexPtr */ - inline StorageIndex* outerIndexPtr() { return Base::m_outerIndex; } - /** \copydoc SparseMatrix::innerNonZeroPtr */ - inline StorageIndex* innerNonZeroPtr() { return Base::m_innerNonZeros; } - //---------------------------------------- - - /** \copydoc SparseMatrix::coeffRef */ - inline Scalar& coeffRef(Index row, Index col) - { - const Index outer = IsRowMajor ? row : col; - const Index inner = IsRowMajor ? col : row; - - Index start = Base::m_outerIndex[outer]; - Index end = Base::isCompressed() ? Base::m_outerIndex[outer+1] : start + Base::m_innerNonZeros[outer]; - eigen_assert(end>=start && "you probably called coeffRef on a non finalized matrix"); - eigen_assert(end>start && "coeffRef cannot be called on a zero coefficient"); - StorageIndex* r = std::lower_bound(&Base::m_innerIndices[start],&Base::m_innerIndices[end],inner); - const Index id = r - &Base::m_innerIndices[0]; - eigen_assert((*r==inner) && (id(Base::m_values)[id]; - } - - inline SparseMapBase(Index rows, Index cols, Index nnz, StorageIndex* outerIndexPtr, StorageIndex* innerIndexPtr, - Scalar* valuePtr, StorageIndex* innerNonZerosPtr = 0) - : Base(rows, cols, nnz, outerIndexPtr, innerIndexPtr, valuePtr, innerNonZerosPtr) - {} - - // for vectors - inline SparseMapBase(Index size, Index nnz, StorageIndex* innerIndexPtr, Scalar* valuePtr) - : Base(size, nnz, innerIndexPtr, valuePtr) - {} - - /** Empty destructor */ - inline ~SparseMapBase() {} - - protected: - inline SparseMapBase() {} -}; - -/** \ingroup SparseCore_Module - * - * \brief Specialization of class Map for SparseMatrix-like storage. - * - * \tparam SparseMatrixType the equivalent sparse matrix type of the referenced data, it must be a template instance of class SparseMatrix. - * - * \sa class Map, class SparseMatrix, class Ref - */ -#ifndef EIGEN_PARSED_BY_DOXYGEN -template -class Map, Options, StrideType> - : public SparseMapBase, Options, StrideType> > -#else -template -class Map - : public SparseMapBase -#endif -{ - public: - typedef SparseMapBase Base; - EIGEN_SPARSE_PUBLIC_INTERFACE(Map) - enum { IsRowMajor = Base::IsRowMajor }; - - public: - - /** Constructs a read-write Map to a sparse matrix of size \a rows x \a cols, containing \a nnz non-zero coefficients, - * stored as a sparse format as defined by the pointers \a outerIndexPtr, \a innerIndexPtr, and \a valuePtr. - * If the optional parameter \a innerNonZerosPtr is the null pointer, then a standard compressed format is assumed. - * - * This constructor is available only if \c SparseMatrixType is non-const. - * - * More details on the expected storage schemes are given in the \ref TutorialSparse "manual pages". - */ - inline Map(Index rows, Index cols, Index nnz, StorageIndex* outerIndexPtr, - StorageIndex* innerIndexPtr, Scalar* valuePtr, StorageIndex* innerNonZerosPtr = 0) - : Base(rows, cols, nnz, outerIndexPtr, innerIndexPtr, valuePtr, innerNonZerosPtr) - {} -#ifndef EIGEN_PARSED_BY_DOXYGEN - /** Empty destructor */ - inline ~Map() {} -}; - -template -class Map, Options, StrideType> - : public SparseMapBase, Options, StrideType> > -{ - public: - typedef SparseMapBase Base; - EIGEN_SPARSE_PUBLIC_INTERFACE(Map) - enum { IsRowMajor = Base::IsRowMajor }; - - public: -#endif - /** This is the const version of the above constructor. - * - * This constructor is available only if \c SparseMatrixType is const, e.g.: - * \code Map > \endcode - */ - inline Map(Index rows, Index cols, Index nnz, const StorageIndex* outerIndexPtr, - const StorageIndex* innerIndexPtr, const Scalar* valuePtr, const StorageIndex* innerNonZerosPtr = 0) - : Base(rows, cols, nnz, outerIndexPtr, innerIndexPtr, valuePtr, innerNonZerosPtr) - {} - - /** Empty destructor */ - inline ~Map() {} -}; - -namespace internal { - -template -struct evaluator, Options, StrideType> > - : evaluator, Options, StrideType> > > -{ - typedef evaluator, Options, StrideType> > > Base; - typedef Map, Options, StrideType> XprType; - evaluator() : Base() {} - explicit evaluator(const XprType &mat) : Base(mat) {} -}; - -template -struct evaluator, Options, StrideType> > - : evaluator, Options, StrideType> > > -{ - typedef evaluator, Options, StrideType> > > Base; - typedef Map, Options, StrideType> XprType; - evaluator() : Base() {} - explicit evaluator(const XprType &mat) : Base(mat) {} -}; - -} - -} // end namespace Eigen - -#endif // EIGEN_SPARSE_MAP_H diff --git a/eigen/Eigen/src/SparseCore/SparseMatrix.h b/eigen/Eigen/src/SparseCore/SparseMatrix.h deleted file mode 100644 index 0a2490b..0000000 --- a/eigen/Eigen/src/SparseCore/SparseMatrix.h +++ /dev/null @@ -1,1403 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2014 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 _StorageIndex the type of the indices. It has to be a \b signed type (e.g., short, int, std::ptrdiff_t). Default is \c int. - * - * \warning In %Eigen 3.2, the undocumented type \c SparseMatrix::Index was improperly defined as the storage index type (e.g., int), - * whereas it is now (starting from %Eigen 3.3) deprecated and always defined as Eigen::Index. - * Codes making use of \c SparseMatrix::Index, might thus likely have to be changed to use \c SparseMatrix::StorageIndex instead. - * - * This class can be extended with the help of the plugin mechanism described on the page - * \ref TopicCustomizing_Plugins by defining the preprocessor symbol \c EIGEN_SPARSEMATRIX_PLUGIN. - */ - -namespace internal { -template -struct traits > -{ - typedef _Scalar Scalar; - typedef _StorageIndex StorageIndex; - typedef Sparse StorageKind; - typedef MatrixXpr XprKind; - enum { - RowsAtCompileTime = Dynamic, - ColsAtCompileTime = Dynamic, - MaxRowsAtCompileTime = Dynamic, - MaxColsAtCompileTime = Dynamic, - Flags = _Options | NestByRefBit | LvalueBit | CompressedAccessBit, - SupportedAccessPatterns = InnerRandomAccessPattern - }; -}; - -template -struct traits, DiagIndex> > -{ - typedef SparseMatrix<_Scalar, _Options, _StorageIndex> MatrixType; - typedef typename ref_selector::type MatrixTypeNested; - typedef typename remove_reference::type _MatrixTypeNested; - - typedef _Scalar Scalar; - typedef Dense StorageKind; - typedef _StorageIndex StorageIndex; - typedef MatrixXpr XprKind; - - enum { - RowsAtCompileTime = Dynamic, - ColsAtCompileTime = 1, - MaxRowsAtCompileTime = Dynamic, - MaxColsAtCompileTime = 1, - Flags = LvalueBit - }; -}; - -template -struct traits, DiagIndex> > - : public traits, DiagIndex> > -{ - enum { - Flags = 0 - }; -}; - -} // end namespace internal - -template -class SparseMatrix - : public SparseCompressedBase > -{ - typedef SparseCompressedBase Base; - using Base::convert_index; - friend class SparseVector<_Scalar,0,_StorageIndex>; - public: - using Base::isCompressed; - using Base::nonZeros; - EIGEN_SPARSE_PUBLIC_INTERFACE(SparseMatrix) - using Base::operator+=; - using Base::operator-=; - - typedef MappedSparseMatrix Map; - typedef Diagonal DiagonalReturnType; - typedef Diagonal ConstDiagonalReturnType; - typedef typename Base::InnerIterator InnerIterator; - typedef typename Base::ReverseInnerIterator ReverseInnerIterator; - - - using Base::IsRowMajor; - typedef internal::CompressedStorage Storage; - enum { - Options = _Options - }; - - typedef typename Base::IndexVector IndexVector; - typedef typename Base::ScalarVector ScalarVector; - protected: - typedef SparseMatrix TransposedSparseMatrix; - - Index m_outerSize; - Index m_innerSize; - StorageIndex* m_outerIndex; - StorageIndex* m_innerNonZeros; // optional, if null then the data is compressed - Storage m_data; - - public: - - /** \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 StorageIndex* 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 StorageIndex* 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 StorageIndex* 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 StorageIndex* 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 StorageIndex* 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 StorageIndex* 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,StorageIndex(inner)); - if((pinnerSize() non zeros if reserve(Index) has not been called earlier. - * In this case, the insertion procedure is optimized for a \e sequential insertion mode where elements are assumed to be - * inserted by increasing outer-indices. - * - * If that's not the case, then it is strongly recommended to either use a triplet-list to assemble the matrix, or to first - * call reserve(const SizesType &) to reserve the appropriate number of non-zero elements per inner vector. - * - * Assuming memory has been appropriately reserved, this function performs a sorted insertion in O(1) - * if the elements of each inner vector are inserted in increasing inner index order, and in O(nnz_j) for a random insertion. - * - */ - Scalar& insert(Index row, Index col); - - public: - - /** Removes all non zeros but keep allocated memory - * - * This function does not free the currently allocated memory. To release as much as memory as possible, - * call \code mat.data().squeeze(); \endcode after resizing it. - * - * \sa resize(Index,Index), data() - */ - inline void setZero() - { - m_data.clear(); - memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(StorageIndex)); - if(m_innerNonZeros) - memset(m_innerNonZeros, 0, (m_outerSize)*sizeof(StorageIndex)); - } - - /** 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. - * - * The type \c SizesType must expose the following interface: - \code - typedef value_type; - const value_type& operator[](i) const; - \endcode - * for \c i in the [0,this->outerSize()[ range. - * Typical choices include std::vector, Eigen::VectorXi, Eigen::VectorXi::Constant, etc. - */ - template - inline void reserve(const SizesType& reserveSizes); - #else - template - inline void reserve(const SizesType& reserveSizes, const typename SizesType::value_type& enableif = - #if (!EIGEN_COMP_MSVC) || (EIGEN_COMP_MSVC>=1500) // MSVC 2005 fails to compile with this typename - typename - #endif - SizesType::value_type()) - { - EIGEN_UNUSED_VARIABLE(enableif); - reserveInnerVectors(reserveSizes); - } - #endif // EIGEN_PARSED_BY_DOXYGEN - protected: - template - inline void reserveInnerVectors(const SizesType& reserveSizes) - { - if(isCompressed()) - { - Index totalReserveSize = 0; - // turn the matrix into non-compressed mode - m_innerNonZeros = static_cast(std::malloc(m_outerSize * sizeof(StorageIndex))); - if (!m_innerNonZeros) internal::throw_std_bad_alloc(); - - // temporarily use m_innerSizes to hold the new starting points. - StorageIndex* newOuterIndex = m_innerNonZeros; - - StorageIndex count = 0; - for(Index j=0; j=0; --j) - { - StorageIndex 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 - { - StorageIndex* newOuterIndex = static_cast(std::malloc((m_outerSize+1)*sizeof(StorageIndex))); - if (!newOuterIndex) internal::throw_std_bad_alloc(); - - StorageIndex 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) - { - StorageIndex 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(Index(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); - - template - void setFromTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func); - - void sumupDuplicates() { collapseDuplicates(internal::scalar_sum_op()); } - - template - void collapseDuplicates(DupFunctor dup_func = DupFunctor()); - - //--- - - /** \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; - - eigen_internal_assert(m_outerIndex!=0 && m_outerSize>0); - - 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(StorageIndex))); - 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 - makeCompressed(); - - StorageIndex 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(); - StorageIndex newInnerSize = convert_index(IsRowMajor ? cols : rows); - - // Deals with inner non zeros - if (m_innerNonZeros) - { - // Resize m_innerNonZeros - StorageIndex *newInnerNonZeros = static_cast(std::realloc(m_innerNonZeros, (m_outerSize + outerChange) * sizeof(StorageIndex))); - if (!newInnerNonZeros) internal::throw_std_bad_alloc(); - m_innerNonZeros = newInnerNonZeros; - - for(Index i=m_outerSize; i(std::malloc((m_outerSize+outerChange+1) * sizeof(StorageIndex))); - 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++) - { - StorageIndex &n = m_innerNonZeros[i]; - StorageIndex 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; - - StorageIndex *newOuterIndex = static_cast(std::realloc(m_outerIndex, (m_outerSize + outerChange + 1) * sizeof(StorageIndex))); - if (!newOuterIndex) internal::throw_std_bad_alloc(); - m_outerIndex = newOuterIndex; - if (outerChange > 0) - { - StorageIndex last = m_outerSize == 0 ? 0 : m_outerIndex[m_outerSize]; - for(Index i=m_outerSize; i(std::malloc((outerSize + 1) * sizeof(StorageIndex))); - 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(StorageIndex)); - } - - /** \internal - * Resize the nonzero vector to \a size */ - void resizeNonZeros(Index size) - { - m_data.resize(size); - } - - /** \returns a const expression of the diagonal coefficients. */ - const ConstDiagonalReturnType diagonal() const { return ConstDiagonalReturnType(*this); } - - /** \returns a read-write expression of the diagonal coefficients. - * \warning If the diagonal entries are written, then all diagonal - * entries \b must already exist, otherwise an assertion will be raised. - */ - DiagonalReturnType diagonal() { return DiagonalReturnType(*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(); - const bool needToTranspose = (Flags & RowMajorBit) != (internal::evaluator::Flags & RowMajorBit); - if (needToTranspose) - *this = other.derived(); - else - { - #ifdef EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN - EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN - #endif - internal::call_assignment_no_alias(*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(); - Base::operator=(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); - } - - /** \brief Copy constructor with in-place evaluation */ - template - explicit SparseMatrix(const DiagonalBase& other) - : Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0) - { - check_template_parameters(); - *this = other.derived(); - } - - /** 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, StorageIndex(rows()-1)); - Eigen::Map(this->m_data.valuePtr(), rows()).setOnes(); - Eigen::Map(this->m_outerIndex, rows()+1).setLinSpaced(0, StorageIndex(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) - { - #ifdef EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN - EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN - #endif - initAssignment(other); - if(other.isCompressed()) - { - internal::smart_copy(other.m_outerIndex, other.m_outerIndex + m_outerSize + 1, m_outerIndex); - m_data = other.m_data; - } - else - { - Base::operator=(other); - } - } - return *this; - } - -#ifndef EIGEN_PARSED_BY_DOXYGEN - template - inline SparseMatrix& operator=(const EigenBase& other) - { return Base::operator=(other.derived()); } -#endif // EIGEN_PARSED_BY_DOXYGEN - - 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); - } - - /** Overloaded for performance */ - Scalar sum() const; - -# 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 - { - StorageIndex m_index; - StorageIndex m_value; - public: - typedef StorageIndex value_type; - SingletonVector(Index i, Index v) - : m_index(convert_index(i)), m_value(convert_index(v)) - {} - - StorageIndex 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) = convert_index(inner); - return (m_data.value(p) = Scalar(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; - }; -}; - -namespace internal { - -template -void set_from_triplets(const InputIterator& begin, const InputIterator& end, SparseMatrixType& mat, DupFunctor dup_func) -{ - enum { IsRowMajor = SparseMatrixType::IsRowMajor }; - typedef typename SparseMatrixType::Scalar Scalar; - typedef typename SparseMatrixType::StorageIndex StorageIndex; - SparseMatrix trMat(mat.rows(),mat.cols()); - - if(begin!=end) - { - // pass 1: count the nnz per inner-vector - typename SparseMatrixType::IndexVector 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.collapseDuplicates(dup_func); - } - - // 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::scalar_sum_op()); -} - -/** The same as setFromTriplets but when duplicates are met the functor \a dup_func is applied: - * \code - * value = dup_func(OldValue, NewValue) - * \endcode - * Here is a C++11 example keeping the latest entry only: - * \code - * mat.setFromTriplets(triplets.begin(), triplets.end(), [] (const Scalar&,const Scalar &b) { return b; }); - * \endcode - */ -template -template -void SparseMatrix::setFromTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func) -{ - internal::set_from_triplets, DupFunctor>(begin, end, *this, dup_func); -} - -/** \internal */ -template -template -void SparseMatrix::collapseDuplicates(DupFunctor dup_func) -{ - eigen_assert(!isCompressed()); - // TODO, in practice we should be able to use m_innerNonZeros for that task - IndexVector wi(innerSize()); - wi.fill(-1); - StorageIndex 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)) = dup_func(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) - - #ifdef EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN - EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN - #endif - - const bool needToTranspose = (Flags & RowMajorBit) != (internal::evaluator::Flags & RowMajorBit); - if (needToTranspose) - { - #ifdef EIGEN_SPARSE_TRANSPOSED_COPY_PLUGIN - EIGEN_SPARSE_TRANSPOSED_COPY_PLUGIN - #endif - // 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_eval::type >::type OtherCopy; - typedef typename internal::remove_all::type _OtherCopy; - typedef internal::evaluator<_OtherCopy> OtherCopyEval; - OtherCopy otherCopy(other.derived()); - OtherCopyEval otherCopyEval(otherCopy); - - 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; jswap(dest); - return *this; - } - else - { - if(other.isRValue()) - { - initAssignment(other.derived()); - } - // there is no special optimization - return Base::operator=(other.derived()); - } -} - -template -typename SparseMatrix<_Scalar,_Options,_StorageIndex>::Scalar& SparseMatrix<_Scalar,_Options,_StorageIndex>::insert(Index row, Index col) -{ - eigen_assert(row>=0 && row=0 && col(std::malloc(m_outerSize * sizeof(StorageIndex))); - if(!m_innerNonZeros) internal::throw_std_bad_alloc(); - - memset(m_innerNonZeros, 0, (m_outerSize)*sizeof(StorageIndex)); - - // pack all inner-vectors to the end of the pre-allocated space - // and allocate the entire free-space to the first inner-vector - StorageIndex end = convert_index(m_data.allocatedSize()); - for(Index j=1; j<=m_outerSize; ++j) - m_outerIndex[j] = end; - } - else - { - // turn the matrix into non-compressed mode - m_innerNonZeros = static_cast(std::malloc(m_outerSize * sizeof(StorageIndex))); - if(!m_innerNonZeros) internal::throw_std_bad_alloc(); - for(Index j=0; j=0 && m_innerNonZeros[j]==0) - m_outerIndex[j--] = p; - - // push back the new element - ++m_innerNonZeros[outer]; - m_data.append(Scalar(0), inner); - - // check for reallocation - if(data_end != m_data.allocatedSize()) - { - // m_data has been reallocated - // -> move remaining inner-vectors back to the end of the free-space - // so that the entire free-space is allocated to the current inner-vector. - eigen_internal_assert(data_end < m_data.allocatedSize()); - StorageIndex new_end = convert_index(m_data.allocatedSize()); - for(Index k=outer+1; k<=m_outerSize; ++k) - if(m_outerIndex[k]==data_end) - m_outerIndex[k] = new_end; - } - return m_data.value(p); - } - - // Second case: the next inner-vector is packed to the end - // and the current inner-vector end match the used-space. - if(m_outerIndex[outer+1]==data_end && m_outerIndex[outer]+m_innerNonZeros[outer]==m_data.size()) - { - eigen_internal_assert(outer+1==m_outerSize || m_innerNonZeros[outer+1]==0); - - // add space for the new element - ++m_innerNonZeros[outer]; - m_data.resize(m_data.size()+1); - - // check for reallocation - if(data_end != m_data.allocatedSize()) - { - // m_data has been reallocated - // -> move remaining inner-vectors back to the end of the free-space - // so that the entire free-space is allocated to the current inner-vector. - eigen_internal_assert(data_end < m_data.allocatedSize()); - StorageIndex new_end = convert_index(m_data.allocatedSize()); - for(Index k=outer+1; k<=m_outerSize; ++k) - if(m_outerIndex[k]==data_end) - m_outerIndex[k] = new_end; - } - - // and insert it at the right position (sorted insertion) - Index startId = m_outerIndex[outer]; - Index p = m_outerIndex[outer]+m_innerNonZeros[outer]-1; - 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) = convert_index(inner); - return (m_data.value(p) = 0); - } - - if(m_data.size() != m_data.allocatedSize()) - { - // make sure the matrix is compatible to random un-compressed insertion: - m_data.resize(m_data.allocatedSize()); - this->reserveInnerVectors(Array::Constant(m_outerSize, 2)); - } - - return insertUncompressed(row,col); -} - -template -EIGEN_DONT_INLINE typename SparseMatrix<_Scalar,_Options,_StorageIndex>::Scalar& SparseMatrix<_Scalar,_Options,_StorageIndex>::insertUncompressed(Index row, Index col) -{ - eigen_assert(!isCompressed()); - - const Index outer = IsRowMajor ? row : col; - const StorageIndex inner = convert_index(IsRowMajor ? col : row); - - Index room = m_outerIndex[outer+1] - m_outerIndex[outer]; - StorageIndex 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 exists, you must call coeffRef to this end"); - - m_innerNonZeros[outer]++; - - m_data.index(p) = inner; - return (m_data.value(p) = Scalar(0)); -} - -template -EIGEN_DONT_INLINE typename SparseMatrix<_Scalar,_Options,_StorageIndex>::Scalar& SparseMatrix<_Scalar,_Options,_StorageIndex>::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] = convert_index(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)) - && (std::size_t(m_outerIndex[outer+1]) == m_data.size()); - - std::size_t startId = m_outerIndex[outer]; - // FIXME let's make sure sizeof(long int) == sizeof(std::size_t) - std::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) = Scalar(0)); -} - -namespace internal { - -template -struct evaluator > - : evaluator > > -{ - typedef evaluator > > Base; - typedef SparseMatrix<_Scalar,_Options,_StorageIndex> SparseMatrixType; - evaluator() : Base() {} - explicit evaluator(const SparseMatrixType &mat) : Base(mat) {} -}; - -} - -} // end namespace Eigen - -#endif // EIGEN_SPARSEMATRIX_H diff --git a/eigen/Eigen/src/SparseCore/SparseMatrixBase.h b/eigen/Eigen/src/SparseCore/SparseMatrixBase.h deleted file mode 100644 index c6b548f..0000000 --- a/eigen/Eigen/src/SparseCore/SparseMatrixBase.h +++ /dev/null @@ -1,405 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2014 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 is the derived type, e.g. a sparse matrix type, or an expression, etc. - * - * This class can be extended with the help of the plugin mechanism described on the page - * \ref TopicCustomizing_Plugins by defining the preprocessor symbol \c EIGEN_SPARSEMATRIXBASE_PLUGIN. - */ -template class SparseMatrixBase - : public EigenBase -{ - public: - - typedef typename internal::traits::Scalar Scalar; - - /** The numeric type of the expression' coefficients, e.g. float, double, int or std::complex, etc. - * - * It is an alias for the Scalar type */ - typedef Scalar value_type; - - typedef typename internal::packet_traits::type PacketScalar; - typedef typename internal::traits::StorageKind StorageKind; - - /** The integer type used to \b store indices within a SparseMatrix. - * For a \c SparseMatrix it an alias of the third template parameter \c IndexType. */ - typedef typename internal::traits::StorageIndex StorageIndex; - - typedef typename internal::add_const_on_value_type_if_arithmetic< - typename internal::packet_traits::type - >::type PacketReturnType; - - typedef SparseMatrixBase StorageBaseType; - - typedef Matrix IndexVector; - typedef Matrix ScalarVector; - - template - Derived& operator=(const EigenBase &other); - - 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". - */ - - 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 Transpose TransposeReturnType; - typedef typename internal::add_const >::type ConstTransposeReturnType; - - // FIXME storage order do not match evaluator storage order - 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 dense matrix */ - typedef Matrix DenseMatrixType; - /** 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 EigenBase Base; - -#endif // not EIGEN_PARSED_BY_DOXYGEN - -#define EIGEN_CURRENT_STORAGE_BASE_CLASS Eigen::SparseMatrixBase -#ifdef EIGEN_PARSED_BY_DOXYGEN -#define EIGEN_DOC_UNARY_ADDONS(METHOD,OP) /**

This method does not change the sparsity of \c *this: the OP is applied to explicitly stored coefficients only. \sa SparseCompressedBase::coeffs()

*/ -#define EIGEN_DOC_BLOCK_ADDONS_NOT_INNER_PANEL /**

\warning This method returns a read-only expression for any sparse matrices. \sa \ref TutorialSparse_SubMatrices "Sparse block operations"

*/ -#define EIGEN_DOC_BLOCK_ADDONS_INNER_PANEL_IF(COND) /**

\warning This method returns a read-write expression for COND sparse matrices only. Otherwise, the returned expression is read-only. \sa \ref TutorialSparse_SubMatrices "Sparse block operations"

*/ -#else -#define EIGEN_DOC_UNARY_ADDONS(X,Y) -#define EIGEN_DOC_BLOCK_ADDONS_NOT_INNER_PANEL -#define EIGEN_DOC_BLOCK_ADDONS_INNER_PANEL_IF(COND) -#endif -# 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_DOC_UNARY_ADDONS -#undef EIGEN_DOC_BLOCK_ADDONS_NOT_INNER_PANEL -#undef EIGEN_DOC_BLOCK_ADDONS_INNER_PANEL_IF - - /** \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 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); - - template - inline Derived& operator=(const SparseMatrixBase& other); - - inline Derived& operator=(const Derived& other); - - protected: - - template - inline Derived& assign(const OtherDerived& other); - - template - inline void assignGeneric(const OtherDerived& other); - - public: - - 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) - { - Nested nm(m.derived()); - internal::evaluator thisEval(nm); - for (Index row=0; row::InnerIterator it(thisEval, row); it; ++it) - { - for ( ; col thisEval(nm); - if (m.cols() == 1) { - Index row = 0; - for (typename internal::evaluator::InnerIterator it(thisEval, 0); it; ++it) - { - for ( ; row trans = m; - s << static_cast >&>(trans); - } - } - return s; - } - - template - Derived& operator+=(const SparseMatrixBase& other); - template - Derived& operator-=(const SparseMatrixBase& other); - - template - Derived& operator+=(const DiagonalBase& other); - template - Derived& operator-=(const DiagonalBase& other); - - template - Derived& operator+=(const EigenBase &other); - template - Derived& operator-=(const EigenBase &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 * diagonal - template - const Product - operator*(const DiagonalBase &other) const - { return Product(derived(), other.derived()); } - - // diagonal * sparse - template friend - const Product - operator*(const DiagonalBase &lhs, const SparseMatrixBase& rhs) - { return Product(lhs.derived(), rhs.derived()); } - - // sparse * sparse - template - const Product - operator*(const SparseMatrixBase &other) const; - - // sparse * dense - template - const Product - operator*(const MatrixBase &other) const - { return Product(derived(), other.derived()); } - - // dense * sparse - template friend - const Product - operator*(const MatrixBase &lhs, const SparseMatrixBase& rhs) - { return Product(lhs.derived(), rhs.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); - - template - inline const TriangularView triangularView() const; - - template struct SelfAdjointViewReturnType { typedef SparseSelfAdjointView Type; }; - template struct ConstSelfAdjointViewReturnType { typedef const SparseSelfAdjointView Type; }; - - template inline - typename ConstSelfAdjointViewReturnType::Type selfadjointView() const; - template inline - typename SelfAdjointViewReturnType::Type selfadjointView(); - - template Scalar dot(const MatrixBase& other) const; - template Scalar dot(const SparseMatrixBase& other) const; - RealScalar squaredNorm() const; - RealScalar norm() const; - RealScalar blueNorm() const; - - TransposeReturnType transpose() { return TransposeReturnType(derived()); } - const ConstTransposeReturnType transpose() const { return ConstTransposeReturnType(derived()); } - const AdjointReturnType adjoint() const { return AdjointReturnType(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; - - DenseMatrixType toDense() const - { - return DenseMatrixType(derived()); - } - - template - bool isApprox(const SparseMatrixBase& other, - const RealScalar& prec = NumTraits::dummy_precision()) const; - - 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; - - inline const SparseView - pruned(const Scalar& reference = Scalar(0), const RealScalar& epsilon = NumTraits::dummy_precision()) const; - - protected: - - bool m_isRValue; - - static inline StorageIndex convert_index(const Index idx) { - return internal::convert_index(idx); - } - private: - template void evalTo(Dest &) const; -}; - -} // end namespace Eigen - -#endif // EIGEN_SPARSEMATRIXBASE_H diff --git a/eigen/Eigen/src/SparseCore/SparsePermutation.h b/eigen/Eigen/src/SparseCore/SparsePermutation.h deleted file mode 100644 index ef38357..0000000 --- a/eigen/Eigen/src/SparseCore/SparsePermutation.h +++ /dev/null @@ -1,178 +0,0 @@ -// 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 permutation_matrix_product -{ - typedef typename nested_eval::type MatrixType; - typedef typename remove_all::type MatrixTypeCleaned; - - typedef typename MatrixTypeCleaned::Scalar Scalar; - typedef typename MatrixTypeCleaned::StorageIndex StorageIndex; - - enum { - SrcStorageOrder = MatrixTypeCleaned::Flags&RowMajorBit ? RowMajor : ColMajor, - MoveOuter = SrcStorageOrder==RowMajor ? Side==OnTheLeft : Side==OnTheRight - }; - - typedef typename internal::conditional, - SparseMatrix >::type ReturnType; - - template - static inline void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr) - { - MatrixType mat(xpr); - if(MoveOuter) - { - SparseMatrix tmp(mat.rows(), mat.cols()); - Matrix sizes(mat.outerSize()); - for(Index j=0; j tmp(mat.rows(), mat.cols()); - Matrix sizes(tmp.outerSize()); - sizes.setZero(); - PermutationMatrix perm_cpy; - if((Side==OnTheLeft) ^ Transposed) - perm_cpy = perm; - else - perm_cpy = perm.transpose(); - - for(Index j=0; j struct product_promote_storage_type { typedef Sparse ret; }; -template struct product_promote_storage_type { typedef Sparse ret; }; - -// TODO, the following two overloads are only needed to define the right temporary type through -// typename traits >::ReturnType -// whereas it should be correctly handled by traits >::PlainObject - -template -struct product_evaluator, ProductTag, PermutationShape, SparseShape> - : public evaluator::ReturnType> -{ - typedef Product XprType; - typedef typename permutation_matrix_product::ReturnType PlainObject; - typedef evaluator Base; - - enum { - Flags = Base::Flags | EvalBeforeNestingBit - }; - - explicit product_evaluator(const XprType& xpr) - : m_result(xpr.rows(), xpr.cols()) - { - ::new (static_cast(this)) Base(m_result); - generic_product_impl::evalTo(m_result, xpr.lhs(), xpr.rhs()); - } - -protected: - PlainObject m_result; -}; - -template -struct product_evaluator, ProductTag, SparseShape, PermutationShape > - : public evaluator::ReturnType> -{ - typedef Product XprType; - typedef typename permutation_matrix_product::ReturnType PlainObject; - typedef evaluator Base; - - enum { - Flags = Base::Flags | EvalBeforeNestingBit - }; - - explicit product_evaluator(const XprType& xpr) - : m_result(xpr.rows(), xpr.cols()) - { - ::new (static_cast(this)) Base(m_result); - generic_product_impl::evalTo(m_result, xpr.lhs(), xpr.rhs()); - } - -protected: - PlainObject m_result; -}; - -} // end namespace internal - -/** \returns the matrix with the permutation applied to the columns - */ -template -inline const Product -operator*(const SparseMatrixBase& matrix, const PermutationBase& perm) -{ return Product(matrix.derived(), perm.derived()); } - -/** \returns the matrix with the permutation applied to the rows - */ -template -inline const Product -operator*( const PermutationBase& perm, const SparseMatrixBase& matrix) -{ return Product(perm.derived(), matrix.derived()); } - - -/** \returns the matrix with the inverse permutation applied to the columns. - */ -template -inline const Product, AliasFreeProduct> -operator*(const SparseMatrixBase& matrix, const InverseImpl& tperm) -{ - return Product, AliasFreeProduct>(matrix.derived(), tperm.derived()); -} - -/** \returns the matrix with the inverse permutation applied to the rows. - */ -template -inline const Product, SparseDerived, AliasFreeProduct> -operator*(const InverseImpl& tperm, const SparseMatrixBase& matrix) -{ - return Product, SparseDerived, AliasFreeProduct>(tperm.derived(), 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 deleted file mode 100644 index 4cbf687..0000000 --- a/eigen/Eigen/src/SparseCore/SparseProduct.h +++ /dev/null @@ -1,169 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2015 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 { - -/** \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 Product -SparseMatrixBase::operator*(const SparseMatrixBase &other) const -{ - return Product(derived(), other.derived()); -} - -namespace internal { - -// sparse * sparse -template -struct generic_product_impl -{ - template - static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) - { - evalTo(dst, lhs, rhs, typename evaluator_traits::Shape()); - } - - // dense += sparse * sparse - template - static void addTo(Dest& dst, const ActualLhs& lhs, const Rhs& rhs, typename enable_if::Shape,DenseShape>::value,int*>::type* = 0) - { - typedef typename nested_eval::type LhsNested; - typedef typename nested_eval::type RhsNested; - LhsNested lhsNested(lhs); - RhsNested rhsNested(rhs); - internal::sparse_sparse_to_dense_product_selector::type, - typename remove_all::type, Dest>::run(lhsNested,rhsNested,dst); - } - - // dense -= sparse * sparse - template - static void subTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, typename enable_if::Shape,DenseShape>::value,int*>::type* = 0) - { - addTo(dst, -lhs, rhs); - } - -protected: - - // sparse = sparse * sparse - template - static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, SparseShape) - { - typedef typename nested_eval::type LhsNested; - typedef typename nested_eval::type RhsNested; - LhsNested lhsNested(lhs); - RhsNested rhsNested(rhs); - internal::conservative_sparse_sparse_product_selector::type, - typename remove_all::type, Dest>::run(lhsNested,rhsNested,dst); - } - - // dense = sparse * sparse - template - static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, DenseShape) - { - dst.setZero(); - addTo(dst, lhs, rhs); - } -}; - -// sparse * sparse-triangular -template -struct generic_product_impl - : public generic_product_impl -{}; - -// sparse-triangular * sparse -template -struct generic_product_impl - : public generic_product_impl -{}; - -// dense = sparse-product (can be sparse*sparse, sparse*perm, etc.) -template< typename DstXprType, typename Lhs, typename Rhs> -struct Assignment, internal::assign_op::Scalar>, Sparse2Dense> -{ - typedef Product SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) - { - Index dstRows = src.rows(); - Index dstCols = src.cols(); - if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) - dst.resize(dstRows, dstCols); - - generic_product_impl::evalTo(dst,src.lhs(),src.rhs()); - } -}; - -// dense += sparse-product (can be sparse*sparse, sparse*perm, etc.) -template< typename DstXprType, typename Lhs, typename Rhs> -struct Assignment, internal::add_assign_op::Scalar>, Sparse2Dense> -{ - typedef Product SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op &) - { - generic_product_impl::addTo(dst,src.lhs(),src.rhs()); - } -}; - -// dense -= sparse-product (can be sparse*sparse, sparse*perm, etc.) -template< typename DstXprType, typename Lhs, typename Rhs> -struct Assignment, internal::sub_assign_op::Scalar>, Sparse2Dense> -{ - typedef Product SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op &) - { - generic_product_impl::subTo(dst,src.lhs(),src.rhs()); - } -}; - -template -struct unary_evaluator >, IteratorBased> - : public evaluator::PlainObject> -{ - typedef SparseView > XprType; - typedef typename XprType::PlainObject PlainObject; - typedef evaluator Base; - - explicit unary_evaluator(const XprType& xpr) - : m_result(xpr.rows(), xpr.cols()) - { - using std::abs; - ::new (static_cast(this)) Base(m_result); - typedef typename nested_eval::type LhsNested; - typedef typename nested_eval::type RhsNested; - LhsNested lhsNested(xpr.nestedExpression().lhs()); - RhsNested rhsNested(xpr.nestedExpression().rhs()); - - internal::sparse_sparse_product_with_pruning_selector::type, - typename remove_all::type, PlainObject>::run(lhsNested,rhsNested,m_result, - abs(xpr.reference())*xpr.epsilon()); - } - -protected: - PlainObject m_result; -}; - -} // end namespace internal - -} // end namespace Eigen - -#endif // EIGEN_SPARSEPRODUCT_H diff --git a/eigen/Eigen/src/SparseCore/SparseRedux.h b/eigen/Eigen/src/SparseCore/SparseRedux.h deleted file mode 100644 index 4587749..0000000 --- a/eigen/Eigen/src/SparseCore/SparseRedux.h +++ /dev/null @@ -1,49 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2014 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); - internal::evaluator thisEval(derived()); - for (Index j=0; j::InnerIterator iter(thisEval,j); iter; ++iter) - res += iter.value(); - return res; -} - -template -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/SparseRef.h b/eigen/Eigen/src/SparseCore/SparseRef.h deleted file mode 100644 index d91f38f..0000000 --- a/eigen/Eigen/src/SparseCore/SparseRef.h +++ /dev/null @@ -1,397 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2015 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_REF_H -#define EIGEN_SPARSE_REF_H - -namespace Eigen { - -enum { - StandardCompressedFormat = 2 /**< used by Ref to specify whether the input storage must be in standard compressed form */ -}; - -namespace internal { - -template class SparseRefBase; - -template -struct traits, _Options, _StrideType> > - : public traits > -{ - typedef SparseMatrix PlainObjectType; - enum { - Options = _Options, - Flags = traits::Flags | CompressedAccessBit | NestByRefBit - }; - - template struct match { - enum { - StorageOrderMatch = PlainObjectType::IsVectorAtCompileTime || Derived::IsVectorAtCompileTime || ((PlainObjectType::Flags&RowMajorBit)==(Derived::Flags&RowMajorBit)), - MatchAtCompileTime = (Derived::Flags&CompressedAccessBit) && StorageOrderMatch - }; - typedef typename internal::conditional::type type; - }; - -}; - -template -struct traits, _Options, _StrideType> > - : public traits, _Options, _StrideType> > -{ - enum { - Flags = (traits >::Flags | CompressedAccessBit | NestByRefBit) & ~LvalueBit - }; -}; - -template -struct traits, _Options, _StrideType> > - : public traits > -{ - typedef SparseVector PlainObjectType; - enum { - Options = _Options, - Flags = traits::Flags | CompressedAccessBit | NestByRefBit - }; - - template struct match { - enum { - MatchAtCompileTime = (Derived::Flags&CompressedAccessBit) && Derived::IsVectorAtCompileTime - }; - typedef typename internal::conditional::type type; - }; - -}; - -template -struct traits, _Options, _StrideType> > - : public traits, _Options, _StrideType> > -{ - enum { - Flags = (traits >::Flags | CompressedAccessBit | NestByRefBit) & ~LvalueBit - }; -}; - -template -struct traits > : public traits {}; - -template class SparseRefBase - : public SparseMapBase -{ -public: - - typedef SparseMapBase Base; - EIGEN_SPARSE_PUBLIC_INTERFACE(SparseRefBase) - - SparseRefBase() - : Base(RowsAtCompileTime==Dynamic?0:RowsAtCompileTime,ColsAtCompileTime==Dynamic?0:ColsAtCompileTime, 0, 0, 0, 0, 0) - {} - -protected: - - template - void construct(Expression& expr) - { - if(expr.outerIndexPtr()==0) - ::new (static_cast(this)) Base(expr.size(), expr.nonZeros(), expr.innerIndexPtr(), expr.valuePtr()); - else - ::new (static_cast(this)) Base(expr.rows(), expr.cols(), expr.nonZeros(), expr.outerIndexPtr(), expr.innerIndexPtr(), expr.valuePtr(), expr.innerNonZeroPtr()); - } -}; - -} // namespace internal - - -/** - * \ingroup SparseCore_Module - * - * \brief A sparse matrix expression referencing an existing sparse expression - * - * \tparam SparseMatrixType the equivalent sparse matrix type of the referenced data, it must be a template instance of class SparseMatrix. - * \tparam Options specifies whether the a standard compressed format is required \c Options is \c #StandardCompressedFormat, or \c 0. - * The default is \c 0. - * - * \sa class Ref - */ -#ifndef EIGEN_PARSED_BY_DOXYGEN -template -class Ref, Options, StrideType > - : public internal::SparseRefBase, Options, StrideType > > -#else -template -class Ref - : public SparseMapBase // yes, that's weird to use Derived here, but that works! -#endif -{ - typedef SparseMatrix PlainObjectType; - typedef internal::traits Traits; - template - inline Ref(const SparseMatrix& expr); - template - inline Ref(const MappedSparseMatrix& expr); - public: - - typedef internal::SparseRefBase Base; - EIGEN_SPARSE_PUBLIC_INTERFACE(Ref) - - - #ifndef EIGEN_PARSED_BY_DOXYGEN - template - inline Ref(SparseMatrix& expr) - { - EIGEN_STATIC_ASSERT(bool(Traits::template match >::MatchAtCompileTime), STORAGE_LAYOUT_DOES_NOT_MATCH); - eigen_assert( ((Options & int(StandardCompressedFormat))==0) || (expr.isCompressed()) ); - Base::construct(expr.derived()); - } - - template - inline Ref(MappedSparseMatrix& expr) - { - EIGEN_STATIC_ASSERT(bool(Traits::template match >::MatchAtCompileTime), STORAGE_LAYOUT_DOES_NOT_MATCH); - eigen_assert( ((Options & int(StandardCompressedFormat))==0) || (expr.isCompressed()) ); - Base::construct(expr.derived()); - } - - template - inline Ref(const SparseCompressedBase& expr) - #else - /** Implicit constructor from any sparse expression (2D matrix or 1D vector) */ - template - inline Ref(SparseCompressedBase& expr) - #endif - { - EIGEN_STATIC_ASSERT(bool(internal::is_lvalue::value), THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY); - EIGEN_STATIC_ASSERT(bool(Traits::template match::MatchAtCompileTime), STORAGE_LAYOUT_DOES_NOT_MATCH); - eigen_assert( ((Options & int(StandardCompressedFormat))==0) || (expr.isCompressed()) ); - Base::construct(expr.const_cast_derived()); - } -}; - -// this is the const ref version -template -class Ref, Options, StrideType> - : public internal::SparseRefBase, Options, StrideType> > -{ - typedef SparseMatrix TPlainObjectType; - typedef internal::traits Traits; - public: - - typedef internal::SparseRefBase Base; - EIGEN_SPARSE_PUBLIC_INTERFACE(Ref) - - template - inline Ref(const SparseMatrixBase& expr) : m_hasCopy(false) - { - construct(expr.derived(), typename Traits::template match::type()); - } - - inline Ref(const Ref& other) : Base(other), m_hasCopy(false) { - // copy constructor shall not copy the m_object, to avoid unnecessary malloc and copy - } - - template - inline Ref(const RefBase& other) : m_hasCopy(false) { - construct(other.derived(), typename Traits::template match::type()); - } - - ~Ref() { - if(m_hasCopy) { - TPlainObjectType* obj = reinterpret_cast(m_object_bytes); - obj->~TPlainObjectType(); - } - } - - protected: - - template - void construct(const Expression& expr,internal::true_type) - { - if((Options & int(StandardCompressedFormat)) && (!expr.isCompressed())) - { - TPlainObjectType* obj = reinterpret_cast(m_object_bytes); - ::new (obj) TPlainObjectType(expr); - m_hasCopy = true; - Base::construct(*obj); - } - else - { - Base::construct(expr); - } - } - - template - void construct(const Expression& expr, internal::false_type) - { - TPlainObjectType* obj = reinterpret_cast(m_object_bytes); - ::new (obj) TPlainObjectType(expr); - m_hasCopy = true; - Base::construct(*obj); - } - - protected: - char m_object_bytes[sizeof(TPlainObjectType)]; - bool m_hasCopy; -}; - - - -/** - * \ingroup SparseCore_Module - * - * \brief A sparse vector expression referencing an existing sparse vector expression - * - * \tparam SparseVectorType the equivalent sparse vector type of the referenced data, it must be a template instance of class SparseVector. - * - * \sa class Ref - */ -#ifndef EIGEN_PARSED_BY_DOXYGEN -template -class Ref, Options, StrideType > - : public internal::SparseRefBase, Options, StrideType > > -#else -template -class Ref - : public SparseMapBase -#endif -{ - typedef SparseVector PlainObjectType; - typedef internal::traits Traits; - template - inline Ref(const SparseVector& expr); - public: - - typedef internal::SparseRefBase Base; - EIGEN_SPARSE_PUBLIC_INTERFACE(Ref) - - #ifndef EIGEN_PARSED_BY_DOXYGEN - template - inline Ref(SparseVector& expr) - { - EIGEN_STATIC_ASSERT(bool(Traits::template match >::MatchAtCompileTime), STORAGE_LAYOUT_DOES_NOT_MATCH); - Base::construct(expr.derived()); - } - - template - inline Ref(const SparseCompressedBase& expr) - #else - /** Implicit constructor from any 1D sparse vector expression */ - template - inline Ref(SparseCompressedBase& expr) - #endif - { - EIGEN_STATIC_ASSERT(bool(internal::is_lvalue::value), THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY); - EIGEN_STATIC_ASSERT(bool(Traits::template match::MatchAtCompileTime), STORAGE_LAYOUT_DOES_NOT_MATCH); - Base::construct(expr.const_cast_derived()); - } -}; - -// this is the const ref version -template -class Ref, Options, StrideType> - : public internal::SparseRefBase, Options, StrideType> > -{ - typedef SparseVector TPlainObjectType; - typedef internal::traits Traits; - public: - - typedef internal::SparseRefBase Base; - EIGEN_SPARSE_PUBLIC_INTERFACE(Ref) - - template - inline Ref(const SparseMatrixBase& expr) : m_hasCopy(false) - { - construct(expr.derived(), typename Traits::template match::type()); - } - - inline Ref(const Ref& other) : Base(other), m_hasCopy(false) { - // copy constructor shall not copy the m_object, to avoid unnecessary malloc and copy - } - - template - inline Ref(const RefBase& other) : m_hasCopy(false) { - construct(other.derived(), typename Traits::template match::type()); - } - - ~Ref() { - if(m_hasCopy) { - TPlainObjectType* obj = reinterpret_cast(m_object_bytes); - obj->~TPlainObjectType(); - } - } - - protected: - - template - void construct(const Expression& expr,internal::true_type) - { - Base::construct(expr); - } - - template - void construct(const Expression& expr, internal::false_type) - { - TPlainObjectType* obj = reinterpret_cast(m_object_bytes); - ::new (obj) TPlainObjectType(expr); - m_hasCopy = true; - Base::construct(*obj); - } - - protected: - char m_object_bytes[sizeof(TPlainObjectType)]; - bool m_hasCopy; -}; - -namespace internal { - -// FIXME shall we introduce a general evaluatior_ref that we can specialize for any sparse object once, and thus remove this copy-pasta thing... - -template -struct evaluator, Options, StrideType> > - : evaluator, Options, StrideType> > > -{ - typedef evaluator, Options, StrideType> > > Base; - typedef Ref, Options, StrideType> XprType; - evaluator() : Base() {} - explicit evaluator(const XprType &mat) : Base(mat) {} -}; - -template -struct evaluator, Options, StrideType> > - : evaluator, Options, StrideType> > > -{ - typedef evaluator, Options, StrideType> > > Base; - typedef Ref, Options, StrideType> XprType; - evaluator() : Base() {} - explicit evaluator(const XprType &mat) : Base(mat) {} -}; - -template -struct evaluator, Options, StrideType> > - : evaluator, Options, StrideType> > > -{ - typedef evaluator, Options, StrideType> > > Base; - typedef Ref, Options, StrideType> XprType; - evaluator() : Base() {} - explicit evaluator(const XprType &mat) : Base(mat) {} -}; - -template -struct evaluator, Options, StrideType> > - : evaluator, Options, StrideType> > > -{ - typedef evaluator, Options, StrideType> > > Base; - typedef Ref, Options, StrideType> XprType; - evaluator() : Base() {} - explicit evaluator(const XprType &mat) : Base(mat) {} -}; - -} - -} // end namespace Eigen - -#endif // EIGEN_SPARSE_REF_H diff --git a/eigen/Eigen/src/SparseCore/SparseSelfAdjointView.h b/eigen/Eigen/src/SparseCore/SparseSelfAdjointView.h deleted file mode 100644 index 65611b3..0000000 --- a/eigen/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ /dev/null @@ -1,656 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2009-2014 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 Mode 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() - */ -namespace internal { - -template -struct traits > : traits { -}; - -template -void permute_symm_to_symm(const MatrixType& mat, SparseMatrix& _dest, const typename MatrixType::StorageIndex* perm = 0); - -template -void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix& _dest, const typename MatrixType::StorageIndex* perm = 0); - -} - -template class SparseSelfAdjointView - : public EigenBase > -{ - public: - - enum { - Mode = _Mode, - TransposeMode = ((Mode & Upper) ? Lower : 0) | ((Mode & Lower) ? Upper : 0), - RowsAtCompileTime = internal::traits::RowsAtCompileTime, - ColsAtCompileTime = internal::traits::ColsAtCompileTime - }; - - typedef EigenBase Base; - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::StorageIndex StorageIndex; - typedef Matrix VectorI; - typedef typename internal::ref_selector::non_const_type MatrixTypeNested; - typedef typename internal::remove_all::type _MatrixTypeNested; - - explicit inline SparseSelfAdjointView(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; } - typename internal::remove_reference::type& matrix() { return m_matrix; } - - /** \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 - Product - operator*(const SparseMatrixBase& rhs) const - { - return Product(*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 - Product - operator*(const SparseMatrixBase& lhs, const SparseSelfAdjointView& rhs) - { - return Product(lhs.derived(), rhs); - } - - /** Efficient sparse self-adjoint matrix times dense vector/matrix product */ - template - Product - operator*(const MatrixBase& rhs) const - { - return Product(*this, rhs.derived()); - } - - /** Efficient dense vector/matrix times sparse self-adjoint matrix product */ - template friend - Product - operator*(const MatrixBase& lhs, const SparseSelfAdjointView& rhs) - { - return Product(lhs.derived(), rhs); - } - - /** 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)); - - /** \returns an expression of P H P^-1 */ - // TODO implement twists in a more evaluator friendly fashion - SparseSymmetricPermutationProduct<_MatrixTypeNested,Mode> twistedBy(const PermutationMatrix& perm) const - { - return SparseSymmetricPermutationProduct<_MatrixTypeNested,Mode>(m_matrix, perm); - } - - template - SparseSelfAdjointView& operator=(const SparseSymmetricPermutationProduct& permutedMatrix) - { - internal::call_assignment_no_alias_no_transpose(*this, permutedMatrix); - 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); - } - - void resize(Index rows, Index cols) - { - EIGEN_ONLY_USED_FOR_DEBUG(rows); - EIGEN_ONLY_USED_FOR_DEBUG(cols); - eigen_assert(rows == this->rows() && cols == this->cols() - && "SparseSelfadjointView::resize() does not actually allow to resize."); - } - - protected: - - MatrixTypeNested m_matrix; - //mutable VectorI m_countPerRow; - //mutable VectorI m_countPerCol; - private: - template void evalTo(Dest &) const; -}; - -/*************************************************************************** -* Implementation of SparseMatrixBase methods -***************************************************************************/ - -template -template -typename SparseMatrixBase::template ConstSelfAdjointViewReturnType::Type SparseMatrixBase::selfadjointView() const -{ - return SparseSelfAdjointView(derived()); -} - -template -template -typename SparseMatrixBase::template SelfAdjointViewReturnType::Type SparseMatrixBase::selfadjointView() -{ - return SparseSelfAdjointView(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 = tmp.template triangularView(); - else - m_matrix += alpha * tmp.template triangularView(); - - return *this; -} - -namespace internal { - -// TODO currently a selfadjoint expression has the form SelfAdjointView<.,.> -// in the future selfadjoint-ness should be defined by the expression traits -// such that Transpose > is valid. (currently TriangularBase::transpose() is overloaded to make it work) -template -struct evaluator_traits > -{ - typedef typename storage_kind_to_evaluator_kind::Kind Kind; - typedef SparseSelfAdjointShape Shape; -}; - -struct SparseSelfAdjoint2Sparse {}; - -template<> struct AssignmentKind { typedef SparseSelfAdjoint2Sparse Kind; }; -template<> struct AssignmentKind { typedef Sparse2Sparse Kind; }; - -template< typename DstXprType, typename SrcXprType, typename Functor> -struct Assignment -{ - typedef typename DstXprType::StorageIndex StorageIndex; - typedef internal::assign_op AssignOpType; - - template - static void run(SparseMatrix &dst, const SrcXprType &src, const AssignOpType&/*func*/) - { - internal::permute_symm_to_fullsymm(src.matrix(), dst); - } - - // FIXME: the handling of += and -= in sparse matrices should be cleanup so that next two overloads could be reduced to: - template - static void run(SparseMatrix &dst, const SrcXprType &src, const AssignFunc& func) - { - SparseMatrix tmp(src.rows(),src.cols()); - run(tmp, src, AssignOpType()); - call_assignment_no_alias_no_transpose(dst, tmp, func); - } - - template - static void run(SparseMatrix &dst, const SrcXprType &src, - const internal::add_assign_op& /* func */) - { - SparseMatrix tmp(src.rows(),src.cols()); - run(tmp, src, AssignOpType()); - dst += tmp; - } - - template - static void run(SparseMatrix &dst, const SrcXprType &src, - const internal::sub_assign_op& /* func */) - { - SparseMatrix tmp(src.rows(),src.cols()); - run(tmp, src, AssignOpType()); - dst -= tmp; - } - - template - static void run(DynamicSparseMatrix& dst, const SrcXprType &src, const AssignOpType&/*func*/) - { - // TODO directly evaluate into dst; - SparseMatrix tmp(dst.rows(),dst.cols()); - internal::permute_symm_to_fullsymm(src.matrix(), tmp); - dst = tmp; - } -}; - -} // end namespace internal - -/*************************************************************************** -* Implementation of sparse self-adjoint time dense matrix -***************************************************************************/ - -namespace internal { - -template -inline void sparse_selfadjoint_time_dense_product(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha) -{ - EIGEN_ONLY_USED_FOR_DEBUG(alpha); - - typedef typename internal::nested_eval::type SparseLhsTypeNested; - typedef typename internal::remove_all::type SparseLhsTypeNestedCleaned; - typedef evaluator LhsEval; - typedef typename LhsEval::InnerIterator LhsIterator; - typedef typename SparseLhsType::Scalar LhsScalar; - - enum { - LhsIsRowMajor = (LhsEval::Flags&RowMajorBit)==RowMajorBit, - ProcessFirstHalf = - ((Mode&(Upper|Lower))==(Upper|Lower)) - || ( (Mode&Upper) && !LhsIsRowMajor) - || ( (Mode&Lower) && LhsIsRowMajor), - ProcessSecondHalf = !ProcessFirstHalf - }; - - SparseLhsTypeNested lhs_nested(lhs); - LhsEval lhsEval(lhs_nested); - - // work on one column at once - for (Index k=0; k::ReturnType rhs_j(alpha*rhs(j,k)); - // accumulator for partial scalar product - typename DenseResType::Scalar res_j(0); - for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i) - { - LhsScalar lhs_ij = i.value(); - if(!LhsIsRowMajor) lhs_ij = numext::conj(lhs_ij); - res_j += lhs_ij * rhs.coeff(i.index(),k); - res(i.index(),k) += numext::conj(lhs_ij) * rhs_j; - } - res.coeffRef(j,k) += alpha * res_j; - - // handle diagonal coeff - if (ProcessFirstHalf && i && (i.index()==j)) - res.coeffRef(j,k) += alpha * i.value() * rhs.coeff(j,k); - } - } -} - - -template -struct generic_product_impl -: generic_product_impl_base > -{ - template - static void scaleAndAddTo(Dest& dst, const LhsView& lhsView, const Rhs& rhs, const typename Dest::Scalar& alpha) - { - typedef typename LhsView::_MatrixTypeNested Lhs; - typedef typename nested_eval::type LhsNested; - typedef typename nested_eval::type RhsNested; - LhsNested lhsNested(lhsView.matrix()); - RhsNested rhsNested(rhs); - - internal::sparse_selfadjoint_time_dense_product(lhsNested, rhsNested, dst, alpha); - } -}; - -template -struct generic_product_impl -: generic_product_impl_base > -{ - template - static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const RhsView& rhsView, const typename Dest::Scalar& alpha) - { - typedef typename RhsView::_MatrixTypeNested Rhs; - typedef typename nested_eval::type LhsNested; - typedef typename nested_eval::type RhsNested; - LhsNested lhsNested(lhs); - RhsNested rhsNested(rhsView.matrix()); - - // transpose everything - Transpose dstT(dst); - internal::sparse_selfadjoint_time_dense_product(rhsNested.transpose(), lhsNested.transpose(), dstT, alpha); - } -}; - -// NOTE: these two overloads are needed to evaluate the sparse selfadjoint view into a full sparse matrix -// TODO: maybe the copy could be handled by generic_product_impl so that these overloads would not be needed anymore - -template -struct product_evaluator, ProductTag, SparseSelfAdjointShape, SparseShape> - : public evaluator::PlainObject> -{ - typedef Product XprType; - typedef typename XprType::PlainObject PlainObject; - typedef evaluator Base; - - product_evaluator(const XprType& xpr) - : m_lhs(xpr.lhs()), m_result(xpr.rows(), xpr.cols()) - { - ::new (static_cast(this)) Base(m_result); - generic_product_impl::evalTo(m_result, m_lhs, xpr.rhs()); - } - -protected: - typename Rhs::PlainObject m_lhs; - PlainObject m_result; -}; - -template -struct product_evaluator, ProductTag, SparseShape, SparseSelfAdjointShape> - : public evaluator::PlainObject> -{ - typedef Product XprType; - typedef typename XprType::PlainObject PlainObject; - typedef evaluator Base; - - product_evaluator(const XprType& xpr) - : m_rhs(xpr.rhs()), m_result(xpr.rows(), xpr.cols()) - { - ::new (static_cast(this)) Base(m_result); - generic_product_impl::evalTo(m_result, xpr.lhs(), m_rhs); - } - -protected: - typename Lhs::PlainObject m_rhs; - PlainObject m_result; -}; - -} // namespace internal - -/*************************************************************************** -* Implementation of symmetric copies and permutations -***************************************************************************/ -namespace internal { - -template -void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix& _dest, const typename MatrixType::StorageIndex* perm) -{ - typedef typename MatrixType::StorageIndex StorageIndex; - typedef typename MatrixType::Scalar Scalar; - typedef SparseMatrix Dest; - typedef Matrix VectorI; - typedef evaluator MatEval; - typedef typename evaluator::InnerIterator MatIterator; - - MatEval matEval(mat); - 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) || ( Mode==Upper && r(it.index()); - Index r = it.row(); - Index c = it.col(); - - StorageIndex jp = perm ? perm[j] : j; - StorageIndex ip = perm ? perm[i] : i; - - if(Mode==(Upper|Lower)) - { - Index k = count[StorageOrderMatch ? jp : ip]++; - dest.innerIndexPtr()[k] = StorageOrderMatch ? ip : jp; - dest.valuePtr()[k] = it.value(); - } - else if(r==c) - { - Index k = count[ip]++; - dest.innerIndexPtr()[k] = ip; - dest.valuePtr()[k] = it.value(); - } - else if(( (Mode&Lower)==Lower && r>c) || ( (Mode&Upper)==Upper && r -void permute_symm_to_symm(const MatrixType& mat, SparseMatrix& _dest, const typename MatrixType::StorageIndex* perm) -{ - typedef typename MatrixType::StorageIndex StorageIndex; - typedef typename MatrixType::Scalar Scalar; - SparseMatrix& dest(_dest.derived()); - typedef Matrix VectorI; - typedef evaluator MatEval; - typedef typename evaluator::InnerIterator MatIterator; - - enum { - SrcOrder = MatrixType::IsRowMajor ? RowMajor : ColMajor, - StorageOrderMatch = int(SrcOrder) == int(DstOrder), - DstMode = DstOrder==RowMajor ? (_DstMode==Upper ? Lower : Upper) : _DstMode, - SrcMode = SrcOrder==RowMajor ? (_SrcMode==Upper ? Lower : Upper) : _SrcMode - }; - - MatEval matEval(mat); - - Index size = mat.rows(); - VectorI count(size); - count.setZero(); - dest.resize(size,size); - for(StorageIndex j = 0; jj)) - continue; - - StorageIndex ip = perm ? perm[i] : i; - count[int(DstMode)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; - } - } - dest.outerIndexPtr()[0] = 0; - for(Index j=0; jj)) - continue; - - StorageIndex jp = perm ? perm[j] : j; - StorageIndex ip = perm? perm[i] : i; - - Index k = count[int(DstMode)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; - dest.innerIndexPtr()[k] = int(DstMode)==int(Lower) ? (std::max)(ip,jp) : (std::min)(ip,jp); - - if(!StorageOrderMatch) std::swap(ip,jp); - if( ((int(DstMode)==int(Lower) && ipjp))) - dest.valuePtr()[k] = numext::conj(it.value()); - else - dest.valuePtr()[k] = it.value(); - } - } -} - -} - -// TODO implement twists in a more evaluator friendly fashion - -namespace internal { - -template -struct traits > : traits { -}; - -} - -template -class SparseSymmetricPermutationProduct - : public EigenBase > -{ - public: - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::StorageIndex StorageIndex; - enum { - RowsAtCompileTime = internal::traits::RowsAtCompileTime, - ColsAtCompileTime = internal::traits::ColsAtCompileTime - }; - protected: - typedef PermutationMatrix Perm; - public: - typedef Matrix VectorI; - typedef typename MatrixType::Nested MatrixTypeNested; - typedef typename internal::remove_all::type NestedExpression; - - 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(); } - - const NestedExpression& matrix() const { return m_matrix; } - const Perm& perm() const { return m_perm; } - - protected: - MatrixTypeNested m_matrix; - const Perm& m_perm; - -}; - -namespace internal { - -template -struct Assignment, internal::assign_op, Sparse2Sparse> -{ - typedef SparseSymmetricPermutationProduct SrcXprType; - typedef typename DstXprType::StorageIndex DstIndex; - template - static void run(SparseMatrix &dst, const SrcXprType &src, const internal::assign_op &) - { - // internal::permute_symm_to_fullsymm(m_matrix,_dest,m_perm.indices().data()); - SparseMatrix tmp; - internal::permute_symm_to_fullsymm(src.matrix(),tmp,src.perm().indices().data()); - dst = tmp; - } - - template - static void run(SparseSelfAdjointView& dst, const SrcXprType &src, const internal::assign_op &) - { - internal::permute_symm_to_symm(src.matrix(),dst.matrix(),src.perm().indices().data()); - } -}; - -} // end namespace internal - -} // end namespace Eigen - -#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H diff --git a/eigen/Eigen/src/SparseCore/SparseSolverBase.h b/eigen/Eigen/src/SparseCore/SparseSolverBase.h deleted file mode 100644 index b4c9a42..0000000 --- a/eigen/Eigen/src/SparseCore/SparseSolverBase.h +++ /dev/null @@ -1,124 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2014 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_SPARSESOLVERBASE_H -#define EIGEN_SPARSESOLVERBASE_H - -namespace Eigen { - -namespace internal { - - /** \internal - * Helper functions to solve with a sparse right-hand-side and result. - * The rhs is decomposed into small vertical panels which are solved through dense temporaries. - */ -template -typename enable_if::type -solve_sparse_through_dense_panels(const Decomposition &dec, const Rhs& rhs, Dest &dest) -{ - EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); - typedef typename Dest::Scalar DestScalar; - // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix. - static const Index NbColsAtOnce = 4; - Index rhsCols = rhs.cols(); - Index size = rhs.rows(); - // the temporary matrices do not need more columns than NbColsAtOnce: - Index tmpCols = (std::min)(rhsCols, NbColsAtOnce); - Eigen::Matrix tmp(size,tmpCols); - Eigen::Matrix tmpX(size,tmpCols); - for(Index k=0; k(rhsCols-k, NbColsAtOnce); - tmp.leftCols(actualCols) = rhs.middleCols(k,actualCols); - tmpX.leftCols(actualCols) = dec.solve(tmp.leftCols(actualCols)); - dest.middleCols(k,actualCols) = tmpX.leftCols(actualCols).sparseView(); - } -} - -// Overload for vector as rhs -template -typename enable_if::type -solve_sparse_through_dense_panels(const Decomposition &dec, const Rhs& rhs, Dest &dest) -{ - typedef typename Dest::Scalar DestScalar; - Index size = rhs.rows(); - Eigen::Matrix rhs_dense(rhs); - Eigen::Matrix dest_dense(size); - dest_dense = dec.solve(rhs_dense); - dest = dest_dense.sparseView(); -} - -} // end namespace internal - -/** \class SparseSolverBase - * \ingroup SparseCore_Module - * \brief A base class for sparse solvers - * - * \tparam Derived the actual type of the solver. - * - */ -template -class SparseSolverBase : internal::noncopyable -{ - public: - - /** Default constructor */ - SparseSolverBase() - : m_isInitialized(false) - {} - - ~SparseSolverBase() - {} - - Derived& derived() { return *static_cast(this); } - const Derived& derived() const { return *static_cast(this); } - - /** \returns an expression of the solution x of \f$ A x = b \f$ using the current decomposition of A. - * - * \sa compute() - */ - template - inline const Solve - solve(const MatrixBase& b) const - { - eigen_assert(m_isInitialized && "Solver is not initialized."); - eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b"); - return Solve(derived(), b.derived()); - } - - /** \returns an expression of the solution x of \f$ A x = b \f$ using the current decomposition of A. - * - * \sa compute() - */ - template - inline const Solve - solve(const SparseMatrixBase& b) const - { - eigen_assert(m_isInitialized && "Solver is not initialized."); - eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b"); - return Solve(derived(), b.derived()); - } - - #ifndef EIGEN_PARSED_BY_DOXYGEN - /** \internal default implementation of solving with a sparse rhs */ - template - void _solve_impl(const SparseMatrixBase &b, SparseMatrixBase &dest) const - { - internal::solve_sparse_through_dense_panels(derived(), b.derived(), dest.derived()); - } - #endif // EIGEN_PARSED_BY_DOXYGEN - - protected: - - mutable bool m_isInitialized; -}; - -} // end namespace Eigen - -#endif // EIGEN_SPARSESOLVERBASE_H diff --git a/eigen/Eigen/src/SparseCore/SparseSparseProductWithPruning.h b/eigen/Eigen/src/SparseCore/SparseSparseProductWithPruning.h deleted file mode 100644 index 88820a4..0000000 --- a/eigen/Eigen/src/SparseCore/SparseSparseProductWithPruning.h +++ /dev/null @@ -1,198 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2014 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 RhsScalar; - typedef typename remove_all::type::Scalar ResScalar; - typedef typename remove_all::type::StorageIndex StorageIndex; - - // 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); - - // mimics a resizeByInnerOuter: - if(ResultType::IsRowMajor) - res.resize(cols, rows); - else - res.resize(rows, cols); - - evaluator lhsEval(lhs); - evaluator rhsEval(rhs); - - // 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 = lhsEval.nonZerosEstimate() + rhsEval.nonZerosEstimate(); - - res.reserve(estimated_nnz_prod); - double ratioColRes = double(estimated_nnz_prod)/(double(lhs.rows())*double(rhs.cols())); - for (Index j=0; j::InnerIterator rhsIt(rhsEval, j); rhsIt; ++rhsIt) - { - // FIXME should be written like this: tmp += rhsIt.value() * lhs.col(rhsIt.index()) - tempVector.restart(); - RhsScalar x = rhsIt.value(); - for (typename evaluator::InnerIterator lhsIt(lhsEval, rhsIt.index()); lhsIt; ++lhsIt) - { - tempVector.coeffRef(lhsIt.index()) += lhsIt.value() * x; - } - } - res.startVec(j); - for (typename AmbiVector::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 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(); - } -}; - -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 RowMajorMatrixLhs; - RowMajorMatrixLhs rowLhs(lhs); - sparse_sparse_product_with_pruning_selector(rowLhs,rhs,res,tolerance); - } -}; - -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 RowMajorMatrixRhs; - RowMajorMatrixRhs rowRhs(rhs); - sparse_sparse_product_with_pruning_selector(lhs,rowRhs,res,tolerance); - } -}; - -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 ColMajorMatrixRhs; - ColMajorMatrixRhs colRhs(rhs); - internal::sparse_sparse_product_with_pruning_impl(lhs, colRhs, res, tolerance); - } -}; - -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; - ColMajorMatrixLhs colLhs(lhs); - internal::sparse_sparse_product_with_pruning_impl(colLhs, rhs, res, tolerance); - } -}; - -} // 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 deleted file mode 100644 index 3757d4c..0000000 --- a/eigen/Eigen/src/SparseCore/SparseTranspose.h +++ /dev/null @@ -1,92 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2015 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 { - -namespace internal { - template - class SparseTransposeImpl - : public SparseMatrixBase > - {}; - - template - class SparseTransposeImpl - : public SparseCompressedBase > - { - typedef SparseCompressedBase > Base; - public: - using Base::derived; - typedef typename Base::Scalar Scalar; - typedef typename Base::StorageIndex StorageIndex; - - inline Index nonZeros() const { return derived().nestedExpression().nonZeros(); } - - inline const Scalar* valuePtr() const { return derived().nestedExpression().valuePtr(); } - inline const StorageIndex* innerIndexPtr() const { return derived().nestedExpression().innerIndexPtr(); } - inline const StorageIndex* outerIndexPtr() const { return derived().nestedExpression().outerIndexPtr(); } - inline const StorageIndex* innerNonZeroPtr() const { return derived().nestedExpression().innerNonZeroPtr(); } - - inline Scalar* valuePtr() { return derived().nestedExpression().valuePtr(); } - inline StorageIndex* innerIndexPtr() { return derived().nestedExpression().innerIndexPtr(); } - inline StorageIndex* outerIndexPtr() { return derived().nestedExpression().outerIndexPtr(); } - inline StorageIndex* innerNonZeroPtr() { return derived().nestedExpression().innerNonZeroPtr(); } - }; -} - -template class TransposeImpl - : public internal::SparseTransposeImpl -{ - protected: - typedef internal::SparseTransposeImpl Base; -}; - -namespace internal { - -template -struct unary_evaluator, IteratorBased> - : public evaluator_base > -{ - typedef typename evaluator::InnerIterator EvalIterator; - public: - typedef Transpose XprType; - - inline Index nonZerosEstimate() const { - return m_argImpl.nonZerosEstimate(); - } - - class InnerIterator : public EvalIterator - { - public: - EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& unaryOp, Index outer) - : EvalIterator(unaryOp.m_argImpl,outer) - {} - - Index row() const { return EvalIterator::col(); } - Index col() const { return EvalIterator::row(); } - }; - - enum { - CoeffReadCost = evaluator::CoeffReadCost, - Flags = XprType::Flags - }; - - explicit unary_evaluator(const XprType& op) :m_argImpl(op.nestedExpression()) {} - - protected: - evaluator m_argImpl; -}; - -} // end namespace internal - -} // end namespace Eigen - -#endif // EIGEN_SPARSETRANSPOSE_H diff --git a/eigen/Eigen/src/SparseCore/SparseTriangularView.h b/eigen/Eigen/src/SparseCore/SparseTriangularView.h deleted file mode 100644 index 9ac1202..0000000 --- a/eigen/Eigen/src/SparseCore/SparseTriangularView.h +++ /dev/null @@ -1,189 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2009-2015 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 { - -/** \ingroup SparseCore_Module - * - * \brief Base class for a triangular part in a \b sparse matrix - * - * This class is an abstract base class of class TriangularView, and objects of type TriangularViewImpl cannot be instantiated. - * It extends class TriangularView with additional methods which are available for sparse expressions only. - * - * \sa class TriangularView, SparseMatrixBase::triangularView() - */ -template class TriangularViewImpl - : 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 - }; - - typedef TriangularView TriangularViewType; - - protected: - // dummy solve function to make TriangularView happy. - void solve() const; - - typedef SparseMatrixBase Base; - public: - - EIGEN_SPARSE_PUBLIC_INTERFACE(TriangularViewType) - - typedef typename MatrixType::Nested MatrixTypeNested; - typedef typename internal::remove_reference::type MatrixTypeNestedNonRef; - typedef typename internal::remove_all::type MatrixTypeNestedCleaned; - - template - EIGEN_DEVICE_FUNC - EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const { - if(!(internal::is_same::value && internal::extract_data(dst) == internal::extract_data(rhs))) - dst = rhs; - this->solveInPlace(dst); - } - - /** Applies the inverse of \c *this to the dense vector or matrix \a other, "in-place" */ - template void solveInPlace(MatrixBase& other) const; - - /** Applies the inverse of \c *this to the sparse vector or matrix \a other, "in-place" */ - template void solveInPlace(SparseMatrixBase& other) const; - -}; - -namespace internal { - -template -struct unary_evaluator, IteratorBased> - : evaluator_base > -{ - typedef TriangularView XprType; - -protected: - - typedef typename XprType::Scalar Scalar; - typedef typename XprType::StorageIndex StorageIndex; - typedef typename evaluator::InnerIterator EvalIterator; - - enum { SkipFirst = ((Mode&Lower) && !(ArgType::Flags&RowMajorBit)) - || ((Mode&Upper) && (ArgType::Flags&RowMajorBit)), - SkipLast = !SkipFirst, - SkipDiag = (Mode&ZeroDiag) ? 1 : 0, - HasUnitDiag = (Mode&UnitDiag) ? 1 : 0 - }; - -public: - - enum { - CoeffReadCost = evaluator::CoeffReadCost, - Flags = XprType::Flags - }; - - explicit unary_evaluator(const XprType &xpr) : m_argImpl(xpr.nestedExpression()), m_arg(xpr.nestedExpression()) {} - - inline Index nonZerosEstimate() const { - return m_argImpl.nonZerosEstimate(); - } - - class InnerIterator : public EvalIterator - { - typedef EvalIterator Base; - public: - - EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& xprEval, Index outer) - : Base(xprEval.m_argImpl,outer), m_returnOne(false), m_containsDiag(Base::outer()index()<=outer : this->index()=Base::outer())) - { - if((!SkipFirst) && Base::operator bool()) - Base::operator++(); - m_returnOne = m_containsDiag; - } - } - - 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 = m_containsDiag; - } - } - return *this; - } - - 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()); - } - } - -// inline Index row() const { return (ArgType::Flags&RowMajorBit ? Base::outer() : this->index()); } -// inline Index col() const { return (ArgType::Flags&RowMajorBit ? this->index() : Base::outer()); } - inline StorageIndex index() const - { - if(HasUnitDiag && m_returnOne) return internal::convert_index(Base::outer()); - else return Base::index(); - } - inline Scalar value() const - { - if(HasUnitDiag && m_returnOne) return Scalar(1); - else return Base::value(); - } - - protected: - bool m_returnOne; - bool m_containsDiag; - private: - Scalar& valueRef(); - }; - -protected: - evaluator m_argImpl; - const ArgType& m_arg; -}; - -} // end namespace internal - -template -template -inline const TriangularView -SparseMatrixBase::triangularView() const -{ - return TriangularView(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 deleted file mode 100644 index 74df0d4..0000000 --- a/eigen/Eigen/src/SparseCore/SparseUtil.h +++ /dev/null @@ -1,178 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2014 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, =) - - -#define EIGEN_SPARSE_PUBLIC_INTERFACE(Derived) \ - EIGEN_GENERIC_PUBLIC_INTERFACE(Derived) - - -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 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 - : sparse_eval::RowsAtCompileTime,traits::ColsAtCompileTime,traits::Flags> -{}; - -template struct sparse_eval { - typedef typename traits::Scalar _Scalar; - typedef typename traits::StorageIndex _StorageIndex; - public: - typedef SparseVector<_Scalar, RowMajor, _StorageIndex> type; -}; - -template struct sparse_eval { - typedef typename traits::Scalar _Scalar; - typedef typename traits::StorageIndex _StorageIndex; - public: - typedef SparseVector<_Scalar, ColMajor, _StorageIndex> type; -}; - -// TODO this seems almost identical to plain_matrix_type -template struct sparse_eval { - typedef typename traits::Scalar _Scalar; - typedef typename traits::StorageIndex _StorageIndex; - enum { _Options = ((Flags&RowMajorBit)==RowMajorBit) ? RowMajor : ColMajor }; - public: - typedef SparseMatrix<_Scalar, _Options, _StorageIndex> 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::StorageIndex _StorageIndex; - enum { _Options = ((evaluator::Flags&RowMajorBit)==RowMajorBit) ? RowMajor : ColMajor }; - public: - typedef SparseMatrix<_Scalar, _Options, _StorageIndex> type; -}; - -template -struct plain_object_eval - : sparse_eval::RowsAtCompileTime,traits::ColsAtCompileTime, evaluator::Flags> -{}; - -template -struct solve_traits -{ - typedef typename sparse_eval::Flags>::type PlainObject; -}; - -template -struct generic_xpr_base -{ - typedef SparseMatrixBase type; -}; - -struct SparseTriangularShape { static std::string debugName() { return "SparseTriangularShape"; } }; -struct SparseSelfAdjointShape { static std::string debugName() { return "SparseSelfAdjointShape"; } }; - -template<> struct glue_shapes { typedef SparseSelfAdjointShape type; }; -template<> struct glue_shapes { typedef SparseTriangularShape 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::StorageIndex > -class Triplet -{ -public: - Triplet() : m_row(0), m_col(0), m_value(0) {} - - Triplet(const StorageIndex& i, const StorageIndex& j, const Scalar& v = Scalar(0)) - : m_row(i), m_col(j), m_value(v) - {} - - /** \returns the row index of the element */ - const StorageIndex& row() const { return m_row; } - - /** \returns the column index of the element */ - const StorageIndex& col() const { return m_col; } - - /** \returns the value of the element */ - const Scalar& value() const { return m_value; } -protected: - StorageIndex 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 deleted file mode 100644 index 19b0fbc..0000000 --- a/eigen/Eigen/src/SparseCore/SparseVector.h +++ /dev/null @@ -1,478 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2015 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 TopicCustomizing_Plugins by defining the preprocessor symbol \c EIGEN_SPARSEVECTOR_PLUGIN. - */ - -namespace internal { -template -struct traits > -{ - typedef _Scalar Scalar; - typedef _StorageIndex StorageIndex; - 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) | CompressedAccessBit, - 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 SparseCompressedBase > -{ - typedef SparseCompressedBase Base; - using Base::convert_index; - 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 StorageIndex* innerIndexPtr() const { return m_data.indexPtr(); } - EIGEN_STRONG_INLINE StorageIndex* innerIndexPtr() { return m_data.indexPtr(); } - - inline const StorageIndex* outerIndexPtr() const { return 0; } - inline StorageIndex* outerIndexPtr() { return 0; } - inline const StorageIndex* innerNonZeroPtr() const { return 0; } - inline StorageIndex* innerNonZeroPtr() { return 0; } - - /** \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=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) = convert_index(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() {} - - /** \copydoc SparseMatrix::prune(const Scalar&,const RealScalar&) */ - void prune(const Scalar& reference, const RealScalar& epsilon = NumTraits::dummy_precision()) - { - m_data.prune(reference,epsilon); - } - - /** Resizes the sparse vector to \a rows x \a cols - * - * This method is provided for compatibility with matrices. - * For a column vector, \a cols must be equal to 1. - * For a row vector, \a rows must be equal to 1. - * - * \sa resize(Index) - */ - void resize(Index rows, Index cols) - { - eigen_assert((IsColVector ? cols : rows)==1 && "Outer dimension must equal 1"); - resize(IsColVector ? rows : cols); - } - - /** Resizes the sparse vector to \a newSize - * This method deletes all entries, thus leaving an empty sparse vector - * - * \sa conservativeResize(), setZero() */ - void resize(Index newSize) - { - m_size = newSize; - m_data.clear(); - } - - /** Resizes the sparse vector to \a newSize, while leaving old values untouched. - * - * If the size of the vector is decreased, then the storage of the out-of bounds coefficients is kept and reserved. - * Call .data().squeeze() to free extra memory. - * - * \sa reserve(), setZero() - */ - void conservativeResize(Index newSize) - { - if (newSize < m_size) - { - Index i = 0; - while (i - inline SparseVector(const SparseMatrixBase& other) - : m_size(0) - { - #ifdef EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN - EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN - #endif - check_template_parameters(); - *this = other.derived(); - } - - inline SparseVector(const SparseVector& other) - : Base(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); - } - - template - inline void swap(SparseMatrix& other) - { - eigen_assert(other.outerSize()==1); - std::swap(m_size, other.m_innerSize); - 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; -}; - -namespace internal { - -template -struct evaluator > - : evaluator_base > -{ - typedef SparseVector<_Scalar,_Options,_Index> SparseVectorType; - typedef evaluator_base Base; - typedef typename SparseVectorType::InnerIterator InnerIterator; - typedef typename SparseVectorType::ReverseInnerIterator ReverseInnerIterator; - - enum { - CoeffReadCost = NumTraits<_Scalar>::ReadCost, - Flags = SparseVectorType::Flags - }; - - evaluator() : Base() {} - - explicit evaluator(const SparseVectorType &mat) : m_matrix(&mat) - { - EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost); - } - - inline Index nonZerosEstimate() const { - return m_matrix->nonZeros(); - } - - operator SparseVectorType&() { return m_matrix->const_cast_derived(); } - operator const SparseVectorType&() const { return *m_matrix; } - - const SparseVectorType *m_matrix; -}; - -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()); - typedef internal::evaluator SrcEvaluatorType; - SrcEvaluatorType srcEval(src); - for(typename SrcEvaluatorType::InnerIterator it(srcEval, 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()); - typedef internal::evaluator SrcEvaluatorType; - SrcEvaluatorType srcEval(src); - for(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 deleted file mode 100644 index 7c4aea7..0000000 --- a/eigen/Eigen/src/SparseCore/SparseView.h +++ /dev/null @@ -1,253 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2011-2014 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::StorageIndex StorageIndex; - typedef Sparse StorageKind; - enum { - Flags = int(traits::Flags) & (RowMajorBit) - }; -}; - -} // end namespace internal - -/** \ingroup SparseCore_Module - * \class SparseView - * - * \brief Expression of a dense or sparse matrix with zero or too small values removed - * - * \tparam MatrixType the type of the object of which we are removing the small entries - * - * This class represents an expression of a given dense or sparse matrix with - * entries smaller than \c reference * \c epsilon are removed. - * It is the return type of MatrixBase::sparseView() and SparseMatrixBase::pruned() - * and most of the time this is the only way it is used. - * - * \sa MatrixBase::sparseView(), SparseMatrixBase::pruned() - */ -template -class SparseView : public SparseMatrixBase > -{ - typedef typename MatrixType::Nested MatrixTypeNested; - typedef typename internal::remove_all::type _MatrixTypeNested; - typedef SparseMatrixBase Base; -public: - EIGEN_SPARSE_PUBLIC_INTERFACE(SparseView) - typedef typename internal::remove_all::type NestedExpression; - - 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) {} - - 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(); } - - /** \returns the nested expression */ - const typename internal::remove_all::type& - nestedExpression() const { return m_matrix; } - - Scalar reference() const { return m_reference; } - RealScalar epsilon() const { return m_epsilon; } - -protected: - MatrixTypeNested m_matrix; - Scalar m_reference; - RealScalar m_epsilon; -}; - -namespace internal { - -// TODO find a way to unify the two following variants -// This is tricky because implementing an inner iterator on top of an IndexBased evaluator is -// not easy because the evaluators do not expose the sizes of the underlying expression. - -template -struct unary_evaluator, IteratorBased> - : public evaluator_base > -{ - typedef typename evaluator::InnerIterator EvalIterator; - public: - typedef SparseView XprType; - - class InnerIterator : public EvalIterator - { - typedef typename XprType::Scalar Scalar; - public: - - EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& sve, Index outer) - : EvalIterator(sve.m_argImpl,outer), m_view(sve.m_view) - { - incrementToNonZero(); - } - - EIGEN_STRONG_INLINE InnerIterator& operator++() - { - EvalIterator::operator++(); - incrementToNonZero(); - return *this; - } - - using EvalIterator::value; - - protected: - const XprType &m_view; - - private: - void incrementToNonZero() - { - while((bool(*this)) && internal::isMuchSmallerThan(value(), m_view.reference(), m_view.epsilon())) - { - EvalIterator::operator++(); - } - } - }; - - enum { - CoeffReadCost = evaluator::CoeffReadCost, - Flags = XprType::Flags - }; - - explicit unary_evaluator(const XprType& xpr) : m_argImpl(xpr.nestedExpression()), m_view(xpr) {} - - protected: - evaluator m_argImpl; - const XprType &m_view; -}; - -template -struct unary_evaluator, IndexBased> - : public evaluator_base > -{ - public: - typedef SparseView XprType; - protected: - enum { IsRowMajor = (XprType::Flags&RowMajorBit)==RowMajorBit }; - typedef typename XprType::Scalar Scalar; - typedef typename XprType::StorageIndex StorageIndex; - public: - - class InnerIterator - { - public: - - EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& sve, Index outer) - : m_sve(sve), m_inner(0), m_outer(outer), m_end(sve.m_view.innerSize()) - { - incrementToNonZero(); - } - - EIGEN_STRONG_INLINE InnerIterator& operator++() - { - m_inner++; - incrementToNonZero(); - return *this; - } - - EIGEN_STRONG_INLINE Scalar value() const - { - return (IsRowMajor) ? m_sve.m_argImpl.coeff(m_outer, m_inner) - : m_sve.m_argImpl.coeff(m_inner, m_outer); - } - - EIGEN_STRONG_INLINE StorageIndex index() const { return m_inner; } - inline Index row() const { return IsRowMajor ? m_outer : index(); } - inline Index col() const { return IsRowMajor ? index() : m_outer; } - - EIGEN_STRONG_INLINE operator bool() const { return m_inner < m_end && m_inner>=0; } - - protected: - const unary_evaluator &m_sve; - Index m_inner; - const Index m_outer; - const Index m_end; - - private: - void incrementToNonZero() - { - while((bool(*this)) && internal::isMuchSmallerThan(value(), m_sve.m_view.reference(), m_sve.m_view.epsilon())) - { - m_inner++; - } - } - }; - - enum { - CoeffReadCost = evaluator::CoeffReadCost, - Flags = XprType::Flags - }; - - explicit unary_evaluator(const XprType& xpr) : m_argImpl(xpr.nestedExpression()), m_view(xpr) {} - - protected: - evaluator m_argImpl; - const XprType &m_view; -}; - -} // end namespace internal - -/** \ingroup SparseCore_Module - * - * \returns a sparse expression of the dense expression \c *this with values smaller than - * \a reference * \a epsilon removed. - * - * This method is typically used when prototyping to convert a quickly assembled dense Matrix \c D to a SparseMatrix \c S: - * \code - * MatrixXd D(n,m); - * SparseMatrix S; - * S = D.sparseView(); // suppress numerical zeros (exact) - * S = D.sparseView(reference); - * S = D.sparseView(reference,epsilon); - * \endcode - * where \a reference is a meaningful non zero reference value, - * and \a epsilon is a tolerance factor defaulting to NumTraits::dummy_precision(). - * - * \sa SparseMatrixBase::pruned(), class SparseView */ -template -const SparseView MatrixBase::sparseView(const Scalar& reference, - const typename NumTraits::Real& epsilon) const -{ - return SparseView(derived(), reference, epsilon); -} - -/** \returns an expression of \c *this with values smaller than - * \a reference * \a epsilon removed. - * - * This method is typically used in conjunction with the product of two sparse matrices - * to automatically prune the smallest values as follows: - * \code - * C = (A*B).pruned(); // suppress 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 -const SparseView -SparseMatrixBase::pruned(const Scalar& reference, - const RealScalar& epsilon) const -{ - return SparseView(derived(), reference, epsilon); -} - -} // end namespace Eigen - -#endif diff --git a/eigen/Eigen/src/SparseCore/TriangularSolver.h b/eigen/Eigen/src/SparseCore/TriangularSolver.h deleted file mode 100644 index f9c56ba..0000000 --- a/eigen/Eigen/src/SparseCore/TriangularSolver.h +++ /dev/null @@ -1,315 +0,0 @@ -// 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; - typedef evaluator LhsEval; - typedef typename evaluator::InnerIterator LhsIterator; - static void run(const Lhs& lhs, Rhs& other) - { - LhsEval lhsEval(lhs); - for(Index col=0 ; col -struct sparse_solve_triangular_selector -{ - typedef typename Rhs::Scalar Scalar; - typedef evaluator LhsEval; - typedef typename evaluator::InnerIterator LhsIterator; - static void run(const Lhs& lhs, Rhs& other) - { - LhsEval lhsEval(lhs); - for(Index col=0 ; col=0 ; --i) - { - Scalar tmp = other.coeff(i,col); - Scalar l_ii(0); - LhsIterator it(lhsEval, i); - while(it && it.index() -struct sparse_solve_triangular_selector -{ - typedef typename Rhs::Scalar Scalar; - typedef evaluator LhsEval; - typedef typename evaluator::InnerIterator LhsIterator; - static void run(const Lhs& lhs, Rhs& other) - { - LhsEval lhsEval(lhs); - for(Index col=0 ; col -struct sparse_solve_triangular_selector -{ - typedef typename Rhs::Scalar Scalar; - typedef evaluator LhsEval; - typedef typename evaluator::InnerIterator LhsIterator; - static void run(const Lhs& lhs, Rhs& other) - { - LhsEval lhsEval(lhs); - for(Index 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 - LhsIterator it(lhsEval, i); - while(it && it.index()!=i) - ++it; - eigen_assert(it && it.index()==i); - other.coeffRef(i,col) /= it.value(); - } - LhsIterator it(lhsEval, i); - for(; it && it.index() -template -void TriangularViewImpl::solveInPlace(MatrixBase& other) const -{ - eigen_assert(derived().cols() == derived().rows() && derived().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(derived().nestedExpression(), otherCopy); - - if (copy) - other = otherCopy; -} -#endif - -// 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::StorageIndex, - typename traits::StorageIndex>::type StorageIndex; - 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(Index 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 - -#ifndef EIGEN_PARSED_BY_DOXYGEN -template -template -void TriangularViewImpl::solveInPlace(SparseMatrixBase& other) const -{ - eigen_assert(derived().cols() == derived().rows() && derived().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(derived().nestedExpression(), other.derived()); - -// if (copy) -// other = otherCopy; -} -#endif - -} // end namespace Eigen - -#endif // EIGEN_SPARSETRIANGULARSOLVER_H -- cgit v1.2.3