mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-03-13 18:37:27 +08:00
update sparse*sparse product: the default is now a conservative algorithm preserving symbolic non zeros. The previous with auto pruning of the small value is avaible doing: (A*B).pruned() or (A*B).pruned(ref) or (A*B).pruned(ref,eps)
This commit is contained in:
Gael Guennebaud
parent
a997dacc67
commit
1ddf88060b
@ -51,8 +51,9 @@ struct Sparse {};
|
||||
#include "src/Sparse/SparseAssign.h"
|
||||
#include "src/Sparse/SparseRedux.h"
|
||||
#include "src/Sparse/SparseFuzzy.h"
|
||||
#include "src/Sparse/ConservativeSparseSparseProduct.h"
|
||||
#include "src/Sparse/SparseSparseProductWithPruning.h"
|
||||
#include "src/Sparse/SparseProduct.h"
|
||||
#include "src/Sparse/SparseSparseProduct.h"
|
||||
#include "src/Sparse/SparseDenseProduct.h"
|
||||
#include "src/Sparse/SparseDiagonalProduct.h"
|
||||
#include "src/Sparse/SparseTriangularView.h"
|
||||
|
||||
@ -299,7 +299,7 @@ class AmbiVector<_Scalar,_Index>::Iterator
|
||||
* In practice, all coefficients having a magnitude smaller than \a epsilon
|
||||
* are skipped.
|
||||
*/
|
||||
Iterator(const AmbiVector& vec, RealScalar epsilon = RealScalar(0.1)*NumTraits<RealScalar>::dummy_precision())
|
||||
Iterator(const AmbiVector& vec, RealScalar epsilon = 0)
|
||||
: m_vector(vec)
|
||||
{
|
||||
m_epsilon = epsilon;
|
||||
@ -315,7 +315,7 @@ class AmbiVector<_Scalar,_Index>::Iterator
|
||||
{
|
||||
ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_vector.m_buffer);
|
||||
m_currentEl = m_vector.m_llStart;
|
||||
while (m_currentEl>=0 && internal::abs(llElements[m_currentEl].value)<m_epsilon)
|
||||
while (m_currentEl>=0 && internal::abs(llElements[m_currentEl].value)<=m_epsilon)
|
||||
m_currentEl = llElements[m_currentEl].next;
|
||||
if (m_currentEl<0)
|
||||
{
|
||||
|
||||
253
Eigen/src/Sparse/ConservativeSparseSparseProduct.h
Normal file
253
Eigen/src/Sparse/ConservativeSparseSparseProduct.h
Normal file
@ -0,0 +1,253 @@
|
||||
// This file is part of Eigen, a lightweight C++ template library
|
||||
// for linear algebra.
|
||||
//
|
||||
// Copyright (C) 2008-2011 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_CONSERVATIVESPARSESPARSEPRODUCT_H
|
||||
#define EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H
|
||||
|
||||
namespace internal {
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef typename remove_all<Lhs>::type::Scalar Scalar;
|
||||
typedef typename remove_all<Lhs>::type::Index Index;
|
||||
|
||||
// make sure to call innerSize/outerSize since we fake the storage order.
|
||||
Index rows = lhs.innerSize();
|
||||
Index cols = rhs.outerSize();
|
||||
eigen_assert(lhs.outerSize() == rhs.innerSize());
|
||||
|
||||
std::vector<bool> mask(rows,false);
|
||||
Matrix<Scalar,Dynamic,1> values(rows);
|
||||
Matrix<Index,Dynamic,1> indices(rows);
|
||||
|
||||
// estimate the number of non zero entries
|
||||
float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols()));
|
||||
float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
|
||||
float ratioRes = (std::min)(ratioLhs * avgNnzPerRhsColumn, 1.f);
|
||||
|
||||
res.setZero();
|
||||
res.reserve(Index(ratioRes*rows*cols));
|
||||
// we compute each column of the result, one after the other
|
||||
for (Index j=0; j<cols; ++j)
|
||||
{
|
||||
|
||||
res.startVec(j);
|
||||
Index nnz = 0;
|
||||
for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
|
||||
{
|
||||
Scalar y = rhsIt.value();
|
||||
Index k = rhsIt.index();
|
||||
for (typename Lhs::InnerIterator lhsIt(lhs, k); lhsIt; ++lhsIt)
|
||||
{
|
||||
Index i = lhsIt.index();
|
||||
Scalar x = lhsIt.value();
|
||||
if(!mask[i])
|
||||
{
|
||||
mask[i] = true;
|
||||
values[i] = x * y;
|
||||
indices[nnz] = i;
|
||||
++nnz;
|
||||
}
|
||||
else
|
||||
values[i] += x * y;
|
||||
}
|
||||
}
|
||||
|
||||
// unordered insertion
|
||||
for(int k=0; k<nnz; ++k)
|
||||
{
|
||||
int i = indices[k];
|
||||
res.insertBackByOuterInnerUnordered(j,i) = values[i];
|
||||
mask[i] = false;
|
||||
}
|
||||
|
||||
#if 0
|
||||
// alternative ordered insertion code:
|
||||
|
||||
int t200 = rows/(log2(200)*1.39);
|
||||
int t = (rows*100)/139;
|
||||
|
||||
// FIXME reserve nnz non zeros
|
||||
// FIXME implement fast sort algorithms for very small nnz
|
||||
// if the result is sparse enough => 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 && nnz<t200) || nnz * log2(nnz) < t)
|
||||
//res.startVec(j);
|
||||
if(true)
|
||||
{
|
||||
if(nnz>1) std::sort(indices.data(),indices.data()+nnz);
|
||||
for(int k=0; k<nnz; ++k)
|
||||
{
|
||||
int i = indices[k];
|
||||
res.insertBackByOuterInner(j,i) = values[i];
|
||||
mask[i] = false;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// dense path
|
||||
for(int i=0; i<rows; ++i)
|
||||
{
|
||||
if(mask[i])
|
||||
{
|
||||
mask[i] = false;
|
||||
res.insertBackByOuterInner(j,i) = values[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
}
|
||||
res.finalize();
|
||||
}
|
||||
|
||||
|
||||
} // end namespace internal
|
||||
|
||||
namespace internal {
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType,
|
||||
int LhsStorageOrder = traits<Lhs>::Flags&RowMajorBit,
|
||||
int RhsStorageOrder = traits<Rhs>::Flags&RowMajorBit,
|
||||
int ResStorageOrder = traits<ResultType>::Flags&RowMajorBit>
|
||||
struct conservative_sparse_sparse_product_selector;
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
|
||||
{
|
||||
typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
|
||||
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
|
||||
typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
|
||||
ColMajorMatrix resCol(lhs.rows(),rhs.cols());
|
||||
conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol);
|
||||
// sort the non zeros:
|
||||
RowMajorMatrix resRow(resCol);
|
||||
res = resRow;
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor,ColMajor>
|
||||
{
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
|
||||
RowMajorMatrix rhsRow = rhs;
|
||||
RowMajorMatrix resRow(lhs.rows(), rhs.cols());
|
||||
conservative_sparse_sparse_product_impl<RowMajorMatrix,Lhs,RowMajorMatrix>(rhsRow, lhs, resRow);
|
||||
res = resRow;
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,RowMajor,ColMajor>
|
||||
{
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
|
||||
RowMajorMatrix lhsRow = lhs;
|
||||
RowMajorMatrix resRow(lhs.rows(), rhs.cols());
|
||||
conservative_sparse_sparse_product_impl<Rhs,RowMajorMatrix,RowMajorMatrix>(rhs, lhsRow, resRow);
|
||||
res = resRow;
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
|
||||
{
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
|
||||
RowMajorMatrix resRow(lhs.rows(), rhs.cols());
|
||||
conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
|
||||
res = resRow;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
|
||||
{
|
||||
typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
|
||||
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
|
||||
ColMajorMatrix resCol(lhs.rows(), rhs.cols());
|
||||
conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol);
|
||||
res = resCol;
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor,RowMajor>
|
||||
{
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
|
||||
ColMajorMatrix lhsCol = lhs;
|
||||
ColMajorMatrix resCol(lhs.rows(), rhs.cols());
|
||||
conservative_sparse_sparse_product_impl<ColMajorMatrix,Rhs,ColMajorMatrix>(lhsCol, rhs, resCol);
|
||||
res = resCol;
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,RowMajor,RowMajor>
|
||||
{
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
|
||||
ColMajorMatrix rhsCol = rhs;
|
||||
ColMajorMatrix resCol(lhs.rows(), rhs.cols());
|
||||
conservative_sparse_sparse_product_impl<Lhs,ColMajorMatrix,ColMajorMatrix>(lhs, rhsCol, resCol);
|
||||
res = resCol;
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
|
||||
{
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
|
||||
typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
|
||||
RowMajorMatrix resRow(lhs.rows(),rhs.cols());
|
||||
conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
|
||||
// sort the non zeros:
|
||||
ColMajorMatrix resCol(resRow);
|
||||
res = resCol;
|
||||
}
|
||||
};
|
||||
|
||||
} // end namespace internal
|
||||
|
||||
|
||||
#endif // EIGEN_CONSERVATIVESPARSESPARSEPRODUCT_H
|
||||
@ -106,9 +106,42 @@ class SparseSparseProduct : internal::no_assignment_operator,
|
||||
|
||||
template<typename Lhs, typename Rhs>
|
||||
EIGEN_STRONG_INLINE SparseSparseProduct(const Lhs& lhs, const Rhs& rhs)
|
||||
: m_lhs(lhs), m_rhs(rhs)
|
||||
: m_lhs(lhs), m_rhs(rhs), m_tolerance(0), m_conservative(true)
|
||||
{
|
||||
eigen_assert(lhs.cols() == rhs.rows());
|
||||
init();
|
||||
}
|
||||
|
||||
template<typename Lhs, typename Rhs>
|
||||
EIGEN_STRONG_INLINE SparseSparseProduct(const Lhs& lhs, const Rhs& rhs, RealScalar tolerance)
|
||||
: m_lhs(lhs), m_rhs(rhs), m_tolerance(tolerance), m_conservative(false)
|
||||
{
|
||||
init();
|
||||
}
|
||||
|
||||
SparseSparseProduct pruned(Scalar reference = 0, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision()) const
|
||||
{
|
||||
return SparseSparseProduct(m_lhs,m_rhs,internal::abs(reference)*epsilon);
|
||||
}
|
||||
|
||||
template<typename Dest>
|
||||
void evalTo(Dest& result) const
|
||||
{
|
||||
if(m_conservative)
|
||||
internal::conservative_sparse_sparse_product_selector<_LhsNested, _RhsNested, Dest>::run(lhs(),rhs(),result);
|
||||
else
|
||||
internal::sparse_sparse_product_with_pruning_selector<_LhsNested, _RhsNested, Dest>::run(lhs(),rhs(),result,m_tolerance);
|
||||
}
|
||||
|
||||
EIGEN_STRONG_INLINE Index rows() const { return m_lhs.rows(); }
|
||||
EIGEN_STRONG_INLINE Index cols() const { return m_rhs.cols(); }
|
||||
|
||||
EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; }
|
||||
EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; }
|
||||
|
||||
protected:
|
||||
void init()
|
||||
{
|
||||
eigen_assert(m_lhs.cols() == m_rhs.rows());
|
||||
|
||||
enum {
|
||||
ProductIsValid = _LhsNested::ColsAtCompileTime==Dynamic
|
||||
@ -127,15 +160,28 @@ class SparseSparseProduct : internal::no_assignment_operator,
|
||||
EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
|
||||
}
|
||||
|
||||
EIGEN_STRONG_INLINE Index rows() const { return m_lhs.rows(); }
|
||||
EIGEN_STRONG_INLINE Index cols() const { return m_rhs.cols(); }
|
||||
|
||||
EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; }
|
||||
EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; }
|
||||
|
||||
protected:
|
||||
LhsNested m_lhs;
|
||||
RhsNested m_rhs;
|
||||
RealScalar m_tolerance;
|
||||
bool m_conservative;
|
||||
};
|
||||
|
||||
// sparse = sparse * sparse
|
||||
template<typename Derived>
|
||||
template<typename Lhs, typename Rhs>
|
||||
inline Derived& SparseMatrixBase<Derived>::operator=(const SparseSparseProduct<Lhs,Rhs>& product)
|
||||
{
|
||||
product.evalTo(derived());
|
||||
return derived();
|
||||
}
|
||||
|
||||
// sparse * sparse
|
||||
template<typename Derived>
|
||||
template<typename OtherDerived>
|
||||
inline const typename SparseSparseProductReturnType<Derived,OtherDerived>::Type
|
||||
SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other) const
|
||||
{
|
||||
return typename SparseSparseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
|
||||
}
|
||||
|
||||
#endif // EIGEN_SPARSEPRODUCT_H
|
||||
|
||||
@ -1,401 +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_SPARSESPARSEPRODUCT_H
|
||||
#define EIGEN_SPARSESPARSEPRODUCT_H
|
||||
|
||||
namespace internal {
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
static void sparse_product_impl2(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef typename remove_all<Lhs>::type::Scalar Scalar;
|
||||
typedef typename remove_all<Lhs>::type::Index Index;
|
||||
|
||||
// make sure to call innerSize/outerSize since we fake the storage order.
|
||||
Index rows = lhs.innerSize();
|
||||
Index cols = rhs.outerSize();
|
||||
eigen_assert(lhs.outerSize() == rhs.innerSize());
|
||||
|
||||
std::vector<bool> mask(rows,false);
|
||||
Matrix<Scalar,Dynamic,1> values(rows);
|
||||
Matrix<Index,Dynamic,1> indices(rows);
|
||||
|
||||
// estimate the number of non zero entries
|
||||
float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols()));
|
||||
float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
|
||||
float ratioRes = (std::min)(ratioLhs * avgNnzPerRhsColumn, 1.f);
|
||||
|
||||
int t200 = rows/(log2(200)*1.39);
|
||||
int t = (rows*100)/139;
|
||||
|
||||
res.resize(rows, cols);
|
||||
res.reserve(Index(ratioRes*rows*cols));
|
||||
// we compute each column of the result, one after the other
|
||||
for (Index j=0; j<cols; ++j)
|
||||
{
|
||||
|
||||
res.startVec(j);
|
||||
Index nnz = 0;
|
||||
for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
|
||||
{
|
||||
Scalar y = rhsIt.value();
|
||||
Index k = rhsIt.index();
|
||||
for (typename Lhs::InnerIterator lhsIt(lhs, k); lhsIt; ++lhsIt)
|
||||
{
|
||||
Index i = lhsIt.index();
|
||||
Scalar x = lhsIt.value();
|
||||
if(!mask[i])
|
||||
{
|
||||
mask[i] = true;
|
||||
// values[i] = x * y;
|
||||
// indices[nnz] = i;
|
||||
++nnz;
|
||||
}
|
||||
else
|
||||
values[i] += x * y;
|
||||
}
|
||||
}
|
||||
// FIXME reserve nnz non zeros
|
||||
// FIXME implement fast sort algorithms for very small nnz
|
||||
// if the result is sparse enough => 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 && nnz<t200) || nnz * log2(nnz) < t)
|
||||
{
|
||||
if(nnz>1) std::sort(indices.data(),indices.data()+nnz);
|
||||
for(int k=0; k<nnz; ++k)
|
||||
{
|
||||
int i = indices[k];
|
||||
res.insertBackNoCheck(j,i) = values[i];
|
||||
mask[i] = false;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// dense path
|
||||
for(int i=0; i<rows; ++i)
|
||||
{
|
||||
if(mask[i])
|
||||
{
|
||||
mask[i] = false;
|
||||
res.insertBackNoCheck(j,i) = values[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
res.finalize();
|
||||
}
|
||||
|
||||
// perform a pseudo in-place sparse * sparse product assuming all matrices are col major
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
static void sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
// return sparse_product_impl2(lhs,rhs,res);
|
||||
|
||||
typedef typename remove_all<Lhs>::type::Scalar Scalar;
|
||||
typedef typename remove_all<Lhs>::type::Index Index;
|
||||
|
||||
// make sure to call innerSize/outerSize since we fake the storage order.
|
||||
Index rows = lhs.innerSize();
|
||||
Index cols = rhs.outerSize();
|
||||
//int size = lhs.outerSize();
|
||||
eigen_assert(lhs.outerSize() == rhs.innerSize());
|
||||
|
||||
// allocate a temporary buffer
|
||||
AmbiVector<Scalar,Index> tempVector(rows);
|
||||
|
||||
// estimate the number of non zero entries
|
||||
float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols()));
|
||||
float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
|
||||
float ratioRes = (std::min)(ratioLhs * avgNnzPerRhsColumn, 1.f);
|
||||
|
||||
// mimics a resizeByInnerOuter:
|
||||
if(ResultType::IsRowMajor)
|
||||
res.resize(cols, rows);
|
||||
else
|
||||
res.resize(rows, cols);
|
||||
|
||||
res.reserve(Index(ratioRes*rows*cols));
|
||||
for (Index j=0; j<cols; ++j)
|
||||
{
|
||||
// let's do a more accurate determination of the nnz ratio for the current column j of res
|
||||
//float ratioColRes = (std::min)(ratioLhs * rhs.innerNonZeros(j), 1.f);
|
||||
// FIXME find a nice way to get the number of nonzeros of a sub matrix (here an inner vector)
|
||||
float ratioColRes = ratioRes;
|
||||
tempVector.init(ratioColRes);
|
||||
tempVector.setZero();
|
||||
for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
|
||||
{
|
||||
// FIXME should be written like this: tmp += rhsIt.value() * lhs.col(rhsIt.index())
|
||||
tempVector.restart();
|
||||
Scalar x = rhsIt.value();
|
||||
for (typename Lhs::InnerIterator lhsIt(lhs, rhsIt.index()); lhsIt; ++lhsIt)
|
||||
{
|
||||
tempVector.coeffRef(lhsIt.index()) += lhsIt.value() * x;
|
||||
}
|
||||
}
|
||||
res.startVec(j);
|
||||
for (typename AmbiVector<Scalar,Index>::Iterator it(tempVector); it; ++it)
|
||||
res.insertBackByOuterInner(j,it.index()) = it.value();
|
||||
}
|
||||
res.finalize();
|
||||
}
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType,
|
||||
int LhsStorageOrder = traits<Lhs>::Flags&RowMajorBit,
|
||||
int RhsStorageOrder = traits<Rhs>::Flags&RowMajorBit,
|
||||
int ResStorageOrder = traits<ResultType>::Flags&RowMajorBit>
|
||||
struct sparse_product_selector;
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
|
||||
{
|
||||
typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
|
||||
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
// std::cerr << __LINE__ << "\n";
|
||||
typename remove_all<ResultType>::type _res(res.rows(), res.cols());
|
||||
sparse_product_impl<Lhs,Rhs,ResultType>(lhs, rhs, _res);
|
||||
res.swap(_res);
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
|
||||
{
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
// std::cerr << __LINE__ << "\n";
|
||||
// we need a col-major matrix to hold the result
|
||||
typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
|
||||
SparseTemporaryType _res(res.rows(), res.cols());
|
||||
sparse_product_impl<Lhs,Rhs,SparseTemporaryType>(lhs, rhs, _res);
|
||||
res = _res;
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
|
||||
{
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
// std::cerr << __LINE__ << "\n";
|
||||
// let's transpose the product to get a column x column product
|
||||
typename remove_all<ResultType>::type _res(res.rows(), res.cols());
|
||||
sparse_product_impl<Rhs,Lhs,ResultType>(rhs, lhs, _res);
|
||||
res.swap(_res);
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
|
||||
{
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
// std::cerr << "here...\n";
|
||||
typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
|
||||
ColMajorMatrix colLhs(lhs);
|
||||
ColMajorMatrix colRhs(rhs);
|
||||
// std::cerr << "more...\n";
|
||||
sparse_product_impl<ColMajorMatrix,ColMajorMatrix,ResultType>(colLhs, colRhs, res);
|
||||
// std::cerr << "OK.\n";
|
||||
|
||||
// let's transpose the product to get a column x column product
|
||||
|
||||
// typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
|
||||
// SparseTemporaryType _res(res.cols(), res.rows());
|
||||
// sparse_product_impl<Rhs,Lhs,SparseTemporaryType>(rhs, lhs, _res);
|
||||
// res = _res.transpose();
|
||||
}
|
||||
};
|
||||
|
||||
// NOTE the 2 others cases (col row *) must never occur since they are caught
|
||||
// by ProductReturnType which transforms it to (col col *) by evaluating rhs.
|
||||
|
||||
} // end namespace internal
|
||||
|
||||
// sparse = sparse * sparse
|
||||
template<typename Derived>
|
||||
template<typename Lhs, typename Rhs>
|
||||
inline Derived& SparseMatrixBase<Derived>::operator=(const SparseSparseProduct<Lhs,Rhs>& product)
|
||||
{
|
||||
// std::cerr << "there..." << typeid(Lhs).name() << " " << typeid(Lhs).name() << " " << (Derived::Flags&&RowMajorBit) << "\n";
|
||||
internal::sparse_product_selector<
|
||||
typename internal::remove_all<Lhs>::type,
|
||||
typename internal::remove_all<Rhs>::type,
|
||||
Derived>::run(product.lhs(),product.rhs(),derived());
|
||||
return derived();
|
||||
}
|
||||
|
||||
namespace internal {
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType,
|
||||
int LhsStorageOrder = traits<Lhs>::Flags&RowMajorBit,
|
||||
int RhsStorageOrder = traits<Rhs>::Flags&RowMajorBit,
|
||||
int ResStorageOrder = traits<ResultType>::Flags&RowMajorBit>
|
||||
struct sparse_product_selector2;
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_product_selector2<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
|
||||
{
|
||||
typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
|
||||
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
sparse_product_impl2<Lhs,Rhs,ResultType>(lhs, rhs, res);
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_product_selector2<Lhs,Rhs,ResultType,RowMajor,ColMajor,ColMajor>
|
||||
{
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
// prevent warnings until the code is fixed
|
||||
EIGEN_UNUSED_VARIABLE(lhs);
|
||||
EIGEN_UNUSED_VARIABLE(rhs);
|
||||
EIGEN_UNUSED_VARIABLE(res);
|
||||
|
||||
// typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
|
||||
// RowMajorMatrix rhsRow = rhs;
|
||||
// RowMajorMatrix resRow(res.rows(), res.cols());
|
||||
// sparse_product_impl2<RowMajorMatrix,Lhs,RowMajorMatrix>(rhsRow, lhs, resRow);
|
||||
// res = resRow;
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_product_selector2<Lhs,Rhs,ResultType,ColMajor,RowMajor,ColMajor>
|
||||
{
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
|
||||
RowMajorMatrix lhsRow = lhs;
|
||||
RowMajorMatrix resRow(res.rows(), res.cols());
|
||||
sparse_product_impl2<Rhs,RowMajorMatrix,RowMajorMatrix>(rhs, lhsRow, resRow);
|
||||
res = resRow;
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_product_selector2<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
|
||||
{
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix;
|
||||
RowMajorMatrix resRow(res.rows(), res.cols());
|
||||
sparse_product_impl2<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow);
|
||||
res = resRow;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_product_selector2<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
|
||||
{
|
||||
typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
|
||||
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
|
||||
ColMajorMatrix resCol(res.rows(), res.cols());
|
||||
sparse_product_impl2<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol);
|
||||
res = resCol;
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_product_selector2<Lhs,Rhs,ResultType,RowMajor,ColMajor,RowMajor>
|
||||
{
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
|
||||
ColMajorMatrix lhsCol = lhs;
|
||||
ColMajorMatrix resCol(res.rows(), res.cols());
|
||||
sparse_product_impl2<ColMajorMatrix,Rhs,ColMajorMatrix>(lhsCol, rhs, resCol);
|
||||
res = resCol;
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_product_selector2<Lhs,Rhs,ResultType,ColMajor,RowMajor,RowMajor>
|
||||
{
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
|
||||
ColMajorMatrix rhsCol = rhs;
|
||||
ColMajorMatrix resCol(res.rows(), res.cols());
|
||||
sparse_product_impl2<Lhs,ColMajorMatrix,ColMajorMatrix>(lhs, rhsCol, resCol);
|
||||
res = resCol;
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_product_selector2<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
|
||||
{
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
|
||||
// ColMajorMatrix lhsTr(lhs);
|
||||
// ColMajorMatrix rhsTr(rhs);
|
||||
// ColMajorMatrix aux(res.rows(), res.cols());
|
||||
// sparse_product_impl2<Rhs,Lhs,ColMajorMatrix>(rhs, lhs, aux);
|
||||
// // ColMajorMatrix aux2 = aux.transpose();
|
||||
// res = aux;
|
||||
typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
|
||||
ColMajorMatrix lhsCol(lhs);
|
||||
ColMajorMatrix rhsCol(rhs);
|
||||
ColMajorMatrix resCol(res.rows(), res.cols());
|
||||
sparse_product_impl2<ColMajorMatrix,ColMajorMatrix,ColMajorMatrix>(lhsCol, rhsCol, resCol);
|
||||
res = resCol;
|
||||
}
|
||||
};
|
||||
|
||||
} // end namespace internal
|
||||
|
||||
template<typename Derived>
|
||||
template<typename Lhs, typename Rhs>
|
||||
inline void SparseMatrixBase<Derived>::_experimentalNewProduct(const Lhs& lhs, const Rhs& rhs)
|
||||
{
|
||||
//derived().resize(lhs.rows(), rhs.cols());
|
||||
internal::sparse_product_selector2<
|
||||
typename internal::remove_all<Lhs>::type,
|
||||
typename internal::remove_all<Rhs>::type,
|
||||
Derived>::run(lhs,rhs,derived());
|
||||
}
|
||||
|
||||
// sparse * sparse
|
||||
template<typename Derived>
|
||||
template<typename OtherDerived>
|
||||
inline const typename SparseSparseProductReturnType<Derived,OtherDerived>::Type
|
||||
SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other) const
|
||||
{
|
||||
return typename SparseSparseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
|
||||
}
|
||||
|
||||
#endif // EIGEN_SPARSESPARSEPRODUCT_H
|
||||
157
Eigen/src/Sparse/SparseSparseProductWithPruning.h
Normal file
157
Eigen/src/Sparse/SparseSparseProductWithPruning.h
Normal file
@ -0,0 +1,157 @@
|
||||
// This file is part of Eigen, a lightweight C++ template library
|
||||
// for linear algebra.
|
||||
//
|
||||
// Copyright (C) 2008-2011 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_SPARSESPARSEPRODUCTWITHPRUNING_H
|
||||
#define EIGEN_SPARSESPARSEPRODUCTWITHPRUNING_H
|
||||
|
||||
namespace internal {
|
||||
|
||||
|
||||
// perform a pseudo in-place sparse * sparse product assuming all matrices are col major
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
static void sparse_sparse_product_with_pruning_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res, typename ResultType::RealScalar tolerance)
|
||||
{
|
||||
// return sparse_sparse_product_with_pruning_impl2(lhs,rhs,res);
|
||||
|
||||
typedef typename remove_all<Lhs>::type::Scalar Scalar;
|
||||
typedef typename remove_all<Lhs>::type::Index Index;
|
||||
|
||||
// make sure to call innerSize/outerSize since we fake the storage order.
|
||||
Index rows = lhs.innerSize();
|
||||
Index cols = rhs.outerSize();
|
||||
//int size = lhs.outerSize();
|
||||
eigen_assert(lhs.outerSize() == rhs.innerSize());
|
||||
|
||||
// allocate a temporary buffer
|
||||
AmbiVector<Scalar,Index> tempVector(rows);
|
||||
|
||||
// estimate the number of non zero entries
|
||||
float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols()));
|
||||
float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
|
||||
float ratioRes = (std::min)(ratioLhs * avgNnzPerRhsColumn, 1.f);
|
||||
|
||||
// mimics a resizeByInnerOuter:
|
||||
if(ResultType::IsRowMajor)
|
||||
res.resize(cols, rows);
|
||||
else
|
||||
res.resize(rows, cols);
|
||||
|
||||
res.reserve(Index(ratioRes*rows*cols));
|
||||
for (Index j=0; j<cols; ++j)
|
||||
{
|
||||
// let's do a more accurate determination of the nnz ratio for the current column j of res
|
||||
//float ratioColRes = (std::min)(ratioLhs * rhs.innerNonZeros(j), 1.f);
|
||||
// FIXME find a nice way to get the number of nonzeros of a sub matrix (here an inner vector)
|
||||
float ratioColRes = ratioRes;
|
||||
tempVector.init(ratioColRes);
|
||||
tempVector.setZero();
|
||||
for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
|
||||
{
|
||||
// FIXME should be written like this: tmp += rhsIt.value() * lhs.col(rhsIt.index())
|
||||
tempVector.restart();
|
||||
Scalar x = rhsIt.value();
|
||||
for (typename Lhs::InnerIterator lhsIt(lhs, rhsIt.index()); lhsIt; ++lhsIt)
|
||||
{
|
||||
tempVector.coeffRef(lhsIt.index()) += lhsIt.value() * x;
|
||||
}
|
||||
}
|
||||
res.startVec(j);
|
||||
for (typename AmbiVector<Scalar,Index>::Iterator it(tempVector,tolerance); it; ++it)
|
||||
res.insertBackByOuterInner(j,it.index()) = it.value();
|
||||
}
|
||||
res.finalize();
|
||||
}
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType,
|
||||
int LhsStorageOrder = traits<Lhs>::Flags&RowMajorBit,
|
||||
int RhsStorageOrder = traits<Rhs>::Flags&RowMajorBit,
|
||||
int ResStorageOrder = traits<ResultType>::Flags&RowMajorBit>
|
||||
struct sparse_sparse_product_with_pruning_selector;
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_sparse_product_with_pruning_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
|
||||
{
|
||||
typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar;
|
||||
typedef typename ResultType::RealScalar RealScalar;
|
||||
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res, RealScalar tolerance)
|
||||
{
|
||||
typename remove_all<ResultType>::type _res(res.rows(), res.cols());
|
||||
sparse_sparse_product_with_pruning_impl<Lhs,Rhs,ResultType>(lhs, rhs, _res, tolerance);
|
||||
res.swap(_res);
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_sparse_product_with_pruning_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
|
||||
{
|
||||
typedef typename ResultType::RealScalar RealScalar;
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res, RealScalar tolerance)
|
||||
{
|
||||
// we need a col-major matrix to hold the result
|
||||
typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
|
||||
SparseTemporaryType _res(res.rows(), res.cols());
|
||||
sparse_sparse_product_with_pruning_impl<Lhs,Rhs,SparseTemporaryType>(lhs, rhs, _res, tolerance);
|
||||
res = _res;
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_sparse_product_with_pruning_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
|
||||
{
|
||||
typedef typename ResultType::RealScalar RealScalar;
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res, RealScalar tolerance)
|
||||
{
|
||||
// let's transpose the product to get a column x column product
|
||||
typename remove_all<ResultType>::type _res(res.rows(), res.cols());
|
||||
sparse_sparse_product_with_pruning_impl<Rhs,Lhs,ResultType>(rhs, lhs, _res, tolerance);
|
||||
res.swap(_res);
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct sparse_sparse_product_with_pruning_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
|
||||
{
|
||||
typedef typename ResultType::RealScalar RealScalar;
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res, RealScalar tolerance)
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix;
|
||||
ColMajorMatrix colLhs(lhs);
|
||||
ColMajorMatrix colRhs(rhs);
|
||||
sparse_sparse_product_with_pruning_impl<ColMajorMatrix,ColMajorMatrix,ResultType>(colLhs, colRhs, res, tolerance);
|
||||
|
||||
// let's transpose the product to get a column x column product
|
||||
// typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
|
||||
// SparseTemporaryType _res(res.cols(), res.rows());
|
||||
// sparse_sparse_product_with_pruning_impl<Rhs,Lhs,SparseTemporaryType>(rhs, lhs, _res);
|
||||
// res = _res.transpose();
|
||||
}
|
||||
};
|
||||
|
||||
// NOTE the 2 others cases (col row *) must never occur since they are caught
|
||||
// by ProductReturnType which transforms it to (col col *) by evaluating rhs.
|
||||
|
||||
} // end namespace internal
|
||||
|
||||
#endif // EIGEN_SPARSESPARSEPRODUCTWITHPRUNING_H
|
||||
@ -94,6 +94,7 @@ template<typename SparseMatrixType> void sparse_product()
|
||||
|
||||
// int c = internal::random<int>(0,depth-1);
|
||||
|
||||
// sparse * sparse
|
||||
VERIFY_IS_APPROX(m4=m2*m3, refMat4=refMat2*refMat3);
|
||||
VERIFY_IS_APPROX(m4=m2t.transpose()*m3, refMat4=refMat2t.transpose()*refMat3);
|
||||
VERIFY_IS_APPROX(m4=m2t.transpose()*m3t.transpose(), refMat4=refMat2t.transpose()*refMat3t.transpose());
|
||||
@ -103,6 +104,11 @@ template<typename SparseMatrixType> void sparse_product()
|
||||
VERIFY_IS_APPROX(m4 = m2*m3*s1, refMat4 = refMat2*refMat3*s1);
|
||||
VERIFY_IS_APPROX(m4 = s2*m2*m3*s1, refMat4 = s2*refMat2*refMat3*s1);
|
||||
|
||||
VERIFY_IS_APPROX(m4=(m2*m3).pruned(0), refMat4=refMat2*refMat3);
|
||||
VERIFY_IS_APPROX(m4=(m2t.transpose()*m3).pruned(0), refMat4=refMat2t.transpose()*refMat3);
|
||||
VERIFY_IS_APPROX(m4=(m2t.transpose()*m3t.transpose()).pruned(0), refMat4=refMat2t.transpose()*refMat3t.transpose());
|
||||
VERIFY_IS_APPROX(m4=(m2*m3t.transpose()).pruned(0), refMat4=refMat2*refMat3t.transpose());
|
||||
|
||||
// sparse * dense
|
||||
VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3);
|
||||
VERIFY_IS_APPROX(dm4=m2*refMat3t.transpose(), refMat4=refMat2*refMat3t.transpose());
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user