mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-06-03 18:24:02 +08:00
332 lines
10 KiB
C++
332 lines
10 KiB
C++
// This file is part of Eigen, a lightweight C++ template library
|
|
// for linear algebra. Eigen itself is part of the KDE project.
|
|
//
|
|
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.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_SPARSEVECTOR_H
|
|
#define EIGEN_SPARSEVECTOR_H
|
|
|
|
/** \class SparseVector
|
|
*
|
|
* \brief a sparse vector class
|
|
*
|
|
* \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.
|
|
*
|
|
*/
|
|
template<typename _Scalar, int _Flags>
|
|
struct ei_traits<SparseVector<_Scalar, _Flags> >
|
|
{
|
|
typedef _Scalar Scalar;
|
|
enum {
|
|
IsColVector = _Flags & RowMajorBit ? 0 : 1,
|
|
|
|
RowsAtCompileTime = IsColVector ? Dynamic : 1,
|
|
ColsAtCompileTime = IsColVector ? 1 : Dynamic,
|
|
MaxRowsAtCompileTime = RowsAtCompileTime,
|
|
MaxColsAtCompileTime = ColsAtCompileTime,
|
|
Flags = SparseBit | _Flags,
|
|
CoeffReadCost = NumTraits<Scalar>::ReadCost,
|
|
SupportedAccessPatterns = FullyCoherentAccessPattern
|
|
};
|
|
};
|
|
|
|
|
|
|
|
template<typename _Scalar, int _Flags>
|
|
class SparseVector
|
|
: public SparseMatrixBase<SparseVector<_Scalar, _Flags> >
|
|
{
|
|
public:
|
|
EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseVector)
|
|
|
|
protected:
|
|
public:
|
|
|
|
typedef SparseMatrixBase<SparseVector> SparseBase;
|
|
enum { IsColVector = ei_traits<SparseVector>::IsColVector };
|
|
|
|
CompressedStorage<Scalar> m_data;
|
|
int m_size;
|
|
|
|
|
|
public:
|
|
|
|
inline int rows() const { return IsColVector ? m_size : 1; }
|
|
inline int cols() const { return IsColVector ? 1 : m_size; }
|
|
inline int innerSize() const { return m_size; }
|
|
inline int outerSize() const { return 1; }
|
|
inline int innerNonZeros(int j) const { ei_assert(j==0); return m_size; }
|
|
|
|
inline const Scalar* _valuePtr() const { return &m_data.value(0); }
|
|
inline Scalar* _valuePtr() { return &m_data.value(0); }
|
|
|
|
inline const int* _innerIndexPtr() const { return &m_data.index(0); }
|
|
inline int* _innerIndexPtr() { return &m_data.index(0); }
|
|
|
|
inline Scalar coeff(int row, int col) const
|
|
{
|
|
ei_assert((IsColVector ? col : row)==0);
|
|
return coeff(IsColVector ? row : col);
|
|
}
|
|
inline Scalar coeff(int i) const
|
|
{
|
|
int start = 0;
|
|
int end = m_data.size();
|
|
if (start==end)
|
|
return Scalar(0);
|
|
else if (end>0 && i==m_data.index(end-1))
|
|
return m_data.value(end-1);
|
|
// ^^ optimization: let's first check if it is the last coefficient
|
|
// (very common in high level algorithms)
|
|
|
|
// TODO move this search to ScalarArray
|
|
const int* r = std::lower_bound(&m_data.index(start),&m_data.index(end-1),i);
|
|
const int id = r-&m_data.index(0);
|
|
return ((*r==i) && (id<end)) ? m_data.value(id) : Scalar(0);
|
|
}
|
|
|
|
inline Scalar& coeffRef(int row, int col)
|
|
{
|
|
ei_assert((IsColVector ? col : row)==0);
|
|
return coeff(IsColVector ? row : col);
|
|
}
|
|
|
|
inline Scalar& coeffRef(int i)
|
|
{
|
|
int start = 0;
|
|
int end = m_data.size();
|
|
ei_assert(end>=start && "you probably called coeffRef on a non finalized vector");
|
|
ei_assert(end>start && "coeffRef cannot be called on a zero coefficient");
|
|
int* r = std::lower_bound(&m_data.index(start),&m_data.index(end),i);
|
|
const int id = r-&m_data.index(0);
|
|
ei_assert((*r==i) && (id<end) && "coeffRef cannot be called on a zero coefficient");
|
|
return m_data.value(id);
|
|
}
|
|
|
|
public:
|
|
|
|
class InnerIterator;
|
|
|
|
inline void setZero() { m_data.clear(); }
|
|
|
|
/** \returns the number of non zero coefficients */
|
|
inline int nonZeros() const { return m_data.size(); }
|
|
|
|
/**
|
|
*/
|
|
inline void reserve(int reserveSize) { m_data.reserve(reserveSize); }
|
|
|
|
/**
|
|
*/
|
|
inline Scalar& fill(int i)
|
|
{
|
|
m_data.append(0, i);
|
|
return m_data.value(m_data.size()-1);
|
|
}
|
|
|
|
/** Like fill() but with random coordinates.
|
|
*/
|
|
inline Scalar& fillrand(int i)
|
|
{
|
|
int startId = 0;
|
|
int id = m_data.size() - 1;
|
|
m_data.resize(id+2);
|
|
|
|
while ( (id >= startId) && (m_data.index(id) > i) )
|
|
{
|
|
m_data.index(id+1) = m_data.index(id);
|
|
m_data.value(id+1) = m_data.value(id);
|
|
--id;
|
|
}
|
|
m_data.index(id+1) = i;
|
|
m_data.value(id+1) = 0;
|
|
return m_data.value(id+1);
|
|
}
|
|
|
|
void resize(int newSize)
|
|
{
|
|
m_size = newSize;
|
|
m_data.clear();
|
|
}
|
|
|
|
void resizeNonZeros(int size) { m_data.resize(size); }
|
|
|
|
inline SparseVector() : m_size(0) { resize(0, 0); }
|
|
|
|
inline SparseVector(int size) : m_size(0) { resize(size); }
|
|
|
|
template<typename OtherDerived>
|
|
inline SparseVector(const MatrixBase<OtherDerived>& other)
|
|
: m_size(0)
|
|
{
|
|
*this = other.derived();
|
|
}
|
|
|
|
inline SparseVector(const SparseVector& other)
|
|
: m_size(0)
|
|
{
|
|
*this = other.derived();
|
|
}
|
|
|
|
inline void swap(SparseVector& other)
|
|
{
|
|
std::swap(m_size, other.m_size);
|
|
m_data.swap(other.m_data);
|
|
}
|
|
|
|
inline SparseVector& operator=(const SparseVector& other)
|
|
{
|
|
if (other.isRValue())
|
|
{
|
|
swap(other.const_cast_derived());
|
|
}
|
|
else
|
|
{
|
|
resize(other.size());
|
|
m_data = other.m_data;
|
|
}
|
|
return *this;
|
|
}
|
|
|
|
// template<typename OtherDerived>
|
|
// inline SparseVector& operator=(const MatrixBase<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 evauluate it if needed
|
|
// typedef typename ei_nested<OtherDerived,2>::type OtherCopy;
|
|
// OtherCopy otherCopy(other.derived());
|
|
// typedef typename ei_cleantype<OtherCopy>::type _OtherCopy;
|
|
//
|
|
// resize(other.rows(), other.cols());
|
|
// Eigen::Map<VectorXi>(m_outerIndex,outerSize()).setZero();
|
|
// // pass 1
|
|
// // FIXME the above copy could be merged with that pass
|
|
// for (int j=0; j<otherCopy.outerSize(); ++j)
|
|
// for (typename _OtherCopy::InnerIterator it(otherCopy, j); it; ++it)
|
|
// ++m_outerIndex[it.index()];
|
|
//
|
|
// // prefix sum
|
|
// int count = 0;
|
|
// VectorXi positions(outerSize());
|
|
// for (int j=0; j<outerSize(); ++j)
|
|
// {
|
|
// int 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 (int j=0; j<otherCopy.outerSize(); ++j)
|
|
// for (typename _OtherCopy::InnerIterator it(otherCopy, j); it; ++it)
|
|
// {
|
|
// int 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 SparseVector& m)
|
|
{
|
|
for (unsigned int i=0; i<m.nonZeros(); ++i)
|
|
s << "(" << m.m_data.value(i) << "," << m.m_data.index(i) << ") ";
|
|
s << std::endl;
|
|
return s;
|
|
}
|
|
|
|
// this specialized version does not seems to be faster
|
|
// Scalar dot(const SparseVector& other) const
|
|
// {
|
|
// int i=0, j=0;
|
|
// Scalar res = 0;
|
|
// asm("#begindot");
|
|
// while (i<nonZeros() && j<other.nonZeros())
|
|
// {
|
|
// if (m_data.index(i)==other.m_data.index(j))
|
|
// {
|
|
// res += m_data.value(i) * ei_conj(other.m_data.value(j));
|
|
// ++i; ++j;
|
|
// }
|
|
// else if (m_data.index(i)<other.m_data.index(j))
|
|
// ++i;
|
|
// else
|
|
// ++j;
|
|
// }
|
|
// asm("#enddot");
|
|
// return res;
|
|
// }
|
|
|
|
/** Destructor */
|
|
inline ~SparseVector() {}
|
|
};
|
|
|
|
template<typename Scalar, int _Flags>
|
|
class SparseVector<Scalar,_Flags>::InnerIterator
|
|
{
|
|
public:
|
|
InnerIterator(const SparseVector& vec, int outer=0)
|
|
: m_vector(vec), m_id(0), m_end(vec.nonZeros())
|
|
{
|
|
ei_assert(outer==0);
|
|
}
|
|
|
|
template<unsigned int Added, unsigned int Removed>
|
|
InnerIterator(const Flagged<SparseVector,Added,Removed>& vec, int outer)
|
|
: m_vector(vec._expression()), m_id(0), m_end(m_vector.nonZeros())
|
|
{}
|
|
|
|
inline InnerIterator& operator++() { m_id++; return *this; }
|
|
|
|
inline Scalar value() const { return m_vector.m_data.value(m_id); }
|
|
inline Scalar& valueRef() { return const_cast<Scalar&>(m_vector.m_data.value(m_id)); }
|
|
|
|
inline int index() const { return m_vector.m_data.index(m_id); }
|
|
inline int row() const { return IsColVector ? index() : 0; }
|
|
inline int col() const { return IsColVector ? 0 : index(); }
|
|
|
|
inline operator bool() const { return (m_id < m_end); }
|
|
|
|
protected:
|
|
const SparseVector& m_vector;
|
|
int m_id;
|
|
const int m_end;
|
|
};
|
|
|
|
#endif // EIGEN_SPARSEVECTOR_H
|