diff options
Diffstat (limited to 'extern/Eigen3/Eigen/src/Sparse/SparseMatrix.h')
| -rw-r--r-- | extern/Eigen3/Eigen/src/Sparse/SparseMatrix.h | 651 |
1 files changed, 0 insertions, 651 deletions
diff --git a/extern/Eigen3/Eigen/src/Sparse/SparseMatrix.h b/extern/Eigen3/Eigen/src/Sparse/SparseMatrix.h deleted file mode 100644 index 0e175ec6e71..00000000000 --- a/extern/Eigen3/Eigen/src/Sparse/SparseMatrix.h +++ /dev/null @@ -1,651 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr> -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. - -#ifndef EIGEN_SPARSEMATRIX_H -#define EIGEN_SPARSEMATRIX_H - -/** \ingroup Sparse_Module - * - * \class SparseMatrix - * - * \brief The main sparse matrix class - * - * This class implements a sparse matrix using the very common compressed row/column storage - * scheme. - * - * \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 RowMajor. The default is 0 which means column-major. - * \tparam _Index the type of the indices. Default is \c int. - * - * See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme. - * - * This class can be extended with the help of the plugin mechanism described on the page - * \ref TopicCustomizingEigen by defining the preprocessor symbol \c EIGEN_SPARSEMATRIX_PLUGIN. - */ - -namespace internal { -template<typename _Scalar, int _Options, typename _Index> -struct traits<SparseMatrix<_Scalar, _Options, _Index> > -{ - typedef _Scalar Scalar; - typedef _Index Index; - typedef Sparse StorageKind; - typedef MatrixXpr XprKind; - enum { - RowsAtCompileTime = Dynamic, - ColsAtCompileTime = Dynamic, - MaxRowsAtCompileTime = Dynamic, - MaxColsAtCompileTime = Dynamic, - Flags = _Options | NestByRefBit | LvalueBit, - CoeffReadCost = NumTraits<Scalar>::ReadCost, - SupportedAccessPatterns = InnerRandomAccessPattern - }; -}; - -} // end namespace internal - -template<typename _Scalar, int _Options, typename _Index> -class SparseMatrix - : public SparseMatrixBase<SparseMatrix<_Scalar, _Options, _Index> > -{ - public: - EIGEN_SPARSE_PUBLIC_INTERFACE(SparseMatrix) -// using Base::operator=; - EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseMatrix, +=) - EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseMatrix, -=) - // FIXME: why are these operator already alvailable ??? - // EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(SparseMatrix, *=) - // EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(SparseMatrix, /=) - - typedef MappedSparseMatrix<Scalar,Flags> Map; - using Base::IsRowMajor; - typedef CompressedStorage<Scalar,Index> Storage; - enum { - Options = _Options - }; - - protected: - - typedef SparseMatrix<Scalar,(Flags&~RowMajorBit)|(IsRowMajor?RowMajorBit:0)> TransposedSparseMatrix; - - Index m_outerSize; - Index m_innerSize; - Index* m_outerIndex; - CompressedStorage<Scalar,Index> m_data; - - public: - - inline Index rows() const { return IsRowMajor ? m_outerSize : m_innerSize; } - inline Index cols() const { return IsRowMajor ? m_innerSize : m_outerSize; } - - inline Index innerSize() const { return m_innerSize; } - inline Index outerSize() const { return m_outerSize; } - inline Index innerNonZeros(Index j) const { return m_outerIndex[j+1]-m_outerIndex[j]; } - - inline const Scalar* _valuePtr() const { return &m_data.value(0); } - inline Scalar* _valuePtr() { return &m_data.value(0); } - - inline const Index* _innerIndexPtr() const { return &m_data.index(0); } - inline Index* _innerIndexPtr() { return &m_data.index(0); } - - inline const Index* _outerIndexPtr() const { return m_outerIndex; } - inline Index* _outerIndexPtr() { return m_outerIndex; } - - inline Storage& data() { return m_data; } - inline const Storage& data() const { return m_data; } - - inline Scalar coeff(Index row, Index col) const - { - const Index outer = IsRowMajor ? row : col; - const Index inner = IsRowMajor ? col : row; - return m_data.atInRange(m_outerIndex[outer], m_outerIndex[outer+1], inner); - } - - inline Scalar& coeffRef(Index row, Index col) - { - const Index outer = IsRowMajor ? row : col; - const Index inner = IsRowMajor ? col : row; - - Index start = m_outerIndex[outer]; - Index end = m_outerIndex[outer+1]; - 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"); - const Index p = m_data.searchLowerIndex(start,end-1,inner); - eigen_assert((p<end) && (m_data.index(p)==inner) && "coeffRef cannot be called on a zero coefficient"); - return m_data.value(p); - } - - public: - - class InnerIterator; - - /** Removes all non zeros */ - inline void setZero() - { - m_data.clear(); - memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(Index)); - } - - /** \returns the number of non zero coefficients */ - inline Index nonZeros() const { return static_cast<Index>(m_data.size()); } - - /** Preallocates \a reserveSize non zeros */ - inline void reserve(Index reserveSize) - { - m_data.reserve(reserveSize); - } - - //--- low level purely coherent filling --- - - /** \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); - } - - /** \sa insertBack, startVec */ - inline Scalar& insertBackByOuterInner(Index outer, Index inner) - { - eigen_assert(size_t(m_outerIndex[outer+1]) == m_data.size() && "Invalid ordered insertion (invalid outer index)"); - eigen_assert( (m_outerIndex[outer+1]-m_outerIndex[outer]==0 || m_data.index(m_data.size()-1)<inner) && "Invalid ordered insertion (invalid inner index)"); - Index p = m_outerIndex[outer+1]; - ++m_outerIndex[outer+1]; - m_data.append(0, inner); - return m_data.value(p); - } - - /** \warning use it only if you know what you are doing */ - inline Scalar& insertBackByOuterInnerUnordered(Index outer, Index inner) - { - Index p = m_outerIndex[outer+1]; - ++m_outerIndex[outer+1]; - m_data.append(0, inner); - return m_data.value(p); - } - - /** \sa insertBack, insertBackByOuterInner */ - inline void startVec(Index outer) - { - eigen_assert(m_outerIndex[outer]==int(m_data.size()) && "You must call startVec for each inner vector sequentially"); - eigen_assert(m_outerIndex[outer+1]==0 && "You must call startVec for each inner vector sequentially"); - m_outerIndex[outer+1] = m_outerIndex[outer]; - } - - //--- - - /** \returns a reference to a novel non zero coefficient with coordinates \a row x \a col. - * The non zero coefficient must \b not already exist. - * - * \warning This function can be extremely slow if the non zero coefficients - * are not inserted in a coherent order. - * - * After an insertion session, you should call the finalize() function. - */ - EIGEN_DONT_INLINE Scalar& insert(Index row, Index col) - { - const Index outer = IsRowMajor ? row : col; - const Index inner = IsRowMajor ? col : row; - - Index previousOuter = outer; - if (m_outerIndex[outer+1]==0) - { - // we start a new inner vector - while (previousOuter>=0 && m_outerIndex[previousOuter]==0) - { - m_outerIndex[previousOuter] = static_cast<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 inserting in, e.g., - // the 2nd inner vector... - bool isLastVec = (!(previousOuter==-1 && m_data.size()!=0)) - && (size_t(m_outerIndex[outer+1]) == m_data.size()); - - size_t startId = m_outerIndex[outer]; - // FIXME let's make sure sizeof(long int) == sizeof(size_t) - size_t p = m_outerIndex[outer+1]; - ++m_outerIndex[outer+1]; - - float 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 float to avoid integers overflows - float nnzEstimate = float(m_outerIndex[outer])*float(m_outerSize)/float(outer+1); - reallocRatio = (nnzEstimate-float(m_data.size()))/float(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.5f),8.f); - } - } - m_data.resize(m_data.size()+1,reallocRatio); - - if (!isLastVec) - { - if (previousOuter==-1) - { - // oops wrong guess. - // let's correct the outer offsets - for (Index k=0; k<=(outer+1); ++k) - m_outerIndex[k] = 0; - Index k=outer+1; - while(m_outerIndex[k]==0) - m_outerIndex[k++] = 1; - while (k<=m_outerSize && m_outerIndex[k]!=0) - m_outerIndex[k++]++; - p = 0; - --k; - k = m_outerIndex[k]-1; - while (k>0) - { - m_data.index(k) = m_data.index(k-1); - m_data.value(k) = m_data.value(k-1); - k--; - } - } - else - { - // we are not inserting into the last inner vec - // update outer indices: - Index j = outer+2; - while (j<=m_outerSize && m_outerIndex[j]!=0) - m_outerIndex[j++]++; - --j; - // shift data of last vecs: - Index k = m_outerIndex[j]-1; - while (k>=Index(p)) - { - m_data.index(k) = m_data.index(k-1); - m_data.value(k) = m_data.value(k-1); - k--; - } - } - } - - while ( (p > startId) && (m_data.index(p-1) > inner) ) - { - m_data.index(p) = m_data.index(p-1); - m_data.value(p) = m_data.value(p-1); - --p; - } - - m_data.index(p) = inner; - return (m_data.value(p) = 0); - } - - - - - /** Must be called after inserting a set of non zero entries. - */ - inline void finalize() - { - Index size = static_cast<Index>(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; - } - } - - /** Suppress all nonzeros which are smaller than \a reference under the tolerence \a epsilon */ - void prune(Scalar reference, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision()) - { - prune(default_prunning_func(reference,epsilon)); - } - - /** Suppress 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<typename KeepFunc> - void prune(const KeepFunc& keep = KeepFunc()) - { - Index k = 0; - for(Index j=0; j<m_outerSize; ++j) - { - Index previousStart = m_outerIndex[j]; - m_outerIndex[j] = k; - Index end = m_outerIndex[j+1]; - for(Index i=previousStart; i<end; ++i) - { - if(keep(IsRowMajor?j:m_data.index(i), IsRowMajor?m_data.index(i):j, m_data.value(i))) - { - m_data.value(k) = m_data.value(i); - m_data.index(k) = m_data.index(i); - ++k; - } - } - } - m_outerIndex[m_outerSize] = k; - m_data.resize(k,0); - } - - /** Resizes the matrix to a \a rows x \a cols matrix and initializes it to zero - * \sa resizeNonZeros(Index), reserve(), setZero() - */ - void resize(Index rows, Index cols) - { - const Index outerSize = IsRowMajor ? rows : cols; - m_innerSize = IsRowMajor ? cols : rows; - m_data.clear(); - if (m_outerSize != outerSize || m_outerSize==0) - { - delete[] m_outerIndex; - m_outerIndex = new Index [outerSize+1]; - m_outerSize = outerSize; - } - memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(Index)); - } - - /** Low level API - * Resize the nonzero vector to \a size */ - void resizeNonZeros(Index size) - { - m_data.resize(size); - } - - /** Default constructor yielding an empty \c 0 \c x \c 0 matrix */ - inline SparseMatrix() - : m_outerSize(-1), m_innerSize(0), m_outerIndex(0) - { - 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) - { - resize(rows, cols); - } - - /** Constructs a sparse matrix from the sparse expression \a other */ - template<typename OtherDerived> - inline SparseMatrix(const SparseMatrixBase<OtherDerived>& other) - : m_outerSize(0), m_innerSize(0), m_outerIndex(0) - { - *this = other.derived(); - } - - /** Copy constructor */ - inline SparseMatrix(const SparseMatrix& other) - : Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0) - { - *this = other.derived(); - } - - /** Swap the content of two sparse matrices of same type (optimization) */ - 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); - m_data.swap(other.m_data); - } - - inline SparseMatrix& operator=(const SparseMatrix& other) - { -// std::cout << "SparseMatrix& operator=(const SparseMatrix& other)\n"; - if (other.isRValue()) - { - swap(other.const_cast_derived()); - } - else - { - resize(other.rows(), other.cols()); - memcpy(m_outerIndex, other.m_outerIndex, (m_outerSize+1)*sizeof(Index)); - m_data = other.m_data; - } - return *this; - } - - #ifndef EIGEN_PARSED_BY_DOXYGEN - template<typename Lhs, typename Rhs> - inline SparseMatrix& operator=(const SparseSparseProduct<Lhs,Rhs>& product) - { return Base::operator=(product); } - - template<typename OtherDerived> - inline SparseMatrix& operator=(const ReturnByValue<OtherDerived>& other) - { return Base::operator=(other); } - - template<typename OtherDerived> - inline SparseMatrix& operator=(const EigenBase<OtherDerived>& other) - { return Base::operator=(other); } - #endif - - template<typename OtherDerived> - EIGEN_DONT_INLINE SparseMatrix& operator=(const SparseMatrixBase<OtherDerived>& other) - { - const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); - if (needToTranspose) - { - // two passes algorithm: - // 1 - compute the number of coeffs per dest inner vector - // 2 - do the actual copy/eval - // Since each coeff of the rhs has to be evaluated twice, let's evaluate it if needed - typedef typename internal::nested<OtherDerived,2>::type OtherCopy; - typedef typename internal::remove_all<OtherCopy>::type _OtherCopy; - OtherCopy otherCopy(other.derived()); - - resize(other.rows(), other.cols()); - Eigen::Map<Matrix<Index, Dynamic, 1> > (m_outerIndex,outerSize()).setZero(); - // pass 1 - // FIXME the above copy could be merged with that pass - for (Index j=0; j<otherCopy.outerSize(); ++j) - for (typename _OtherCopy::InnerIterator it(otherCopy, j); it; ++it) - ++m_outerIndex[it.index()]; - - // prefix sum - Index count = 0; - VectorXi positions(outerSize()); - for (Index j=0; j<outerSize(); ++j) - { - Index tmp = m_outerIndex[j]; - m_outerIndex[j] = count; - positions[j] = count; - count += tmp; - } - m_outerIndex[outerSize()] = count; - // alloc - m_data.resize(count); - // pass 2 - for (Index j=0; j<otherCopy.outerSize(); ++j) - { - for (typename _OtherCopy::InnerIterator it(otherCopy, j); it; ++it) - { - Index pos = positions[it.index()]++; - m_data.index(pos) = j; - m_data.value(pos) = it.value(); - } - } - return *this; - } - else - { - // there is no special optimization - return SparseMatrixBase<SparseMatrix>::operator=(other.derived()); - } - } - - friend std::ostream & operator << (std::ostream & s, const SparseMatrix& m) - { - EIGEN_DBG_SPARSE( - s << "Nonzero entries:\n"; - for (Index i=0; i<m.nonZeros(); ++i) - { - s << "(" << m.m_data.value(i) << "," << m.m_data.index(i) << ") "; - } - s << std::endl; - s << std::endl; - s << "Column pointers:\n"; - for (Index i=0; i<m.outerSize(); ++i) - { - s << m.m_outerIndex[i] << " "; - } - s << " $" << std::endl; - s << std::endl; - ); - s << static_cast<const SparseMatrixBase<SparseMatrix>&>(m); - return s; - } - - /** Destructor */ - inline ~SparseMatrix() - { - delete[] m_outerIndex; - } - - /** Overloaded for performance */ - Scalar sum() const; - - public: - - /** \deprecated use setZero() and reserve() - * Initializes the filling process of \c *this. - * \param reserveSize approximate number of nonzeros - * Note that the matrix \c *this is zero-ed. - */ - EIGEN_DEPRECATED void startFill(Index reserveSize = 1000) - { - setZero(); - m_data.reserve(reserveSize); - } - - /** \deprecated use insert() - * Like fill() but with random inner coordinates. - */ - EIGEN_DEPRECATED Scalar& fillrand(Index row, Index col) - { - return insert(row,col); - } - - /** \deprecated use insert() - */ - EIGEN_DEPRECATED Scalar& fill(Index row, Index col) - { - const Index outer = IsRowMajor ? row : col; - const Index inner = IsRowMajor ? col : row; - - if (m_outerIndex[outer+1]==0) - { - // we start a new inner vector - Index i = outer; - while (i>=0 && m_outerIndex[i]==0) - { - m_outerIndex[i] = m_data.size(); - --i; - } - m_outerIndex[outer+1] = m_outerIndex[outer]; - } - else - { - eigen_assert(m_data.index(m_data.size()-1)<inner && "wrong sorted insertion"); - } -// std::cerr << size_t(m_outerIndex[outer+1]) << " == " << m_data.size() << "\n"; - assert(size_t(m_outerIndex[outer+1]) == m_data.size()); - Index p = m_outerIndex[outer+1]; - ++m_outerIndex[outer+1]; - - m_data.append(0, inner); - return m_data.value(p); - } - - /** \deprecated use finalize */ - EIGEN_DEPRECATED void endFill() { finalize(); } - -# ifdef EIGEN_SPARSEMATRIX_PLUGIN -# include EIGEN_SPARSEMATRIX_PLUGIN -# endif - -private: - struct default_prunning_func { - default_prunning_func(Scalar ref, RealScalar eps) : reference(ref), epsilon(eps) {} - inline bool operator() (const Index&, const Index&, const Scalar& value) const - { - return !internal::isMuchSmallerThan(value, reference, epsilon); - } - Scalar reference; - RealScalar epsilon; - }; -}; - -template<typename Scalar, int _Options, typename _Index> -class SparseMatrix<Scalar,_Options,_Index>::InnerIterator -{ - public: - InnerIterator(const SparseMatrix& mat, Index outer) - : m_values(mat._valuePtr()), m_indices(mat._innerIndexPtr()), m_outer(outer), m_id(mat.m_outerIndex[outer]), m_end(mat.m_outerIndex[outer+1]) - {} - - inline InnerIterator& operator++() { m_id++; return *this; } - - inline const Scalar& value() const { return m_values[m_id]; } - inline Scalar& valueRef() { return const_cast<Scalar&>(m_values[m_id]); } - - inline Index index() const { return m_indices[m_id]; } - inline Index outer() const { return m_outer; } - inline Index row() const { return IsRowMajor ? m_outer : index(); } - inline Index col() const { return IsRowMajor ? index() : m_outer; } - - inline operator bool() const { return (m_id < m_end); } - - protected: - const Scalar* m_values; - const Index* m_indices; - const Index m_outer; - Index m_id; - const Index m_end; -}; - -#endif // EIGEN_SPARSEMATRIX_H |