mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-07-08 22:21:49 +08:00
this is essentially backporting all the changes made in the Sparse module up to KDE SVN revision r945600, aka changeset:
df9dfa145547529bf71afd4c6e8f3af947acaad0 This is what is needed to make Step (in KDE/kdeedu) build. The rest of Eigen (outside of Sparse) is unaffected except for a few trivial changes that were needed. calling this 2.0.3, will tag if no problem.
This commit is contained in:
Benoit Jacob
parent
12570d97ce
commit
5a18f7545d
@ -1,5 +1,5 @@
|
||||
project(Eigen)
|
||||
set(EIGEN_VERSION_NUMBER "2.0.2")
|
||||
set(EIGEN_VERSION_NUMBER "2.0.3")
|
||||
|
||||
#if the svnversion program is absent, this will leave the SVN_REVISION string empty,
|
||||
#but won't stop CMake.
|
||||
|
||||
@ -22,7 +22,6 @@
|
||||
#endif
|
||||
|
||||
#ifdef EIGEN_TAUCS_SUPPORT
|
||||
|
||||
// taucs.h declares a lot of mess
|
||||
#define isnan
|
||||
#define finite
|
||||
@ -40,7 +39,9 @@
|
||||
#ifdef max
|
||||
#undef max
|
||||
#endif
|
||||
|
||||
#ifdef complex
|
||||
#undef complex
|
||||
#endif
|
||||
#endif
|
||||
|
||||
#ifdef EIGEN_SUPERLU_SUPPORT
|
||||
@ -102,6 +103,7 @@ namespace Eigen {
|
||||
#include "src/Sparse/SparseFuzzy.h"
|
||||
#include "src/Sparse/SparseFlagged.h"
|
||||
#include "src/Sparse/SparseProduct.h"
|
||||
#include "src/Sparse/SparseDiagonalProduct.h"
|
||||
#include "src/Sparse/TriangularSolver.h"
|
||||
#include "src/Sparse/SparseLLT.h"
|
||||
#include "src/Sparse/SparseLDLT.h"
|
||||
|
||||
@ -62,6 +62,7 @@ class DiagonalMatrix : ei_no_assignment_operator,
|
||||
public:
|
||||
|
||||
EIGEN_GENERIC_PUBLIC_INTERFACE(DiagonalMatrix)
|
||||
typedef CoeffsVectorType _CoeffsVectorType;
|
||||
|
||||
// needed to evaluate a DiagonalMatrix<Xpr> to a DiagonalMatrix<NestByValue<Vector> >
|
||||
template<typename OtherCoeffsVectorType>
|
||||
|
||||
@ -30,7 +30,7 @@
|
||||
|
||||
#define EIGEN_WORLD_VERSION 2
|
||||
#define EIGEN_MAJOR_VERSION 0
|
||||
#define EIGEN_MINOR_VERSION 2
|
||||
#define EIGEN_MINOR_VERSION 3
|
||||
|
||||
#define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
|
||||
(EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \
|
||||
|
||||
@ -73,7 +73,9 @@
|
||||
YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY,
|
||||
THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES,
|
||||
THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES,
|
||||
INVALID_MATRIX_TEMPLATE_PARAMETERS
|
||||
INVALID_MATRIX_TEMPLATE_PARAMETERS,
|
||||
BOTH_MATRICES_MUST_HAVE_THE_SAME_STORAGE_ORDER,
|
||||
THIS_METHOD_IS_ONLY_FOR_DIAGONAL_MATRIX
|
||||
};
|
||||
};
|
||||
|
||||
|
||||
@ -2,5 +2,5 @@ FILE(GLOB Eigen_Sparse_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_Sparse_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Sparse
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Sparse COMPONENT Devel
|
||||
)
|
||||
|
||||
@ -1,7 +1,7 @@
|
||||
// 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>
|
||||
// Copyright (C) 2008-2009 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
|
||||
@ -54,16 +54,17 @@ void ei_cholmod_configure_matrix(CholmodType& mat)
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Scalar, int Flags>
|
||||
cholmod_sparse SparseMatrixBase<Scalar,Flags>::asCholmodMatrix()
|
||||
template<typename Derived>
|
||||
cholmod_sparse SparseMatrixBase<Derived>::asCholmodMatrix()
|
||||
{
|
||||
typedef typename Derived::Scalar Scalar;
|
||||
cholmod_sparse res;
|
||||
res.nzmax = nonZeros();
|
||||
res.nrow = rows();;
|
||||
res.ncol = cols();
|
||||
res.p = _outerIndexPtr();
|
||||
res.i = _innerIndexPtr();
|
||||
res.x = _valuePtr();
|
||||
res.p = derived()._outerIndexPtr();
|
||||
res.i = derived()._innerIndexPtr();
|
||||
res.x = derived()._valuePtr();
|
||||
res.xtype = CHOLMOD_REAL;
|
||||
res.itype = CHOLMOD_INT;
|
||||
res.sorted = 1;
|
||||
@ -73,11 +74,11 @@ cholmod_sparse SparseMatrixBase<Scalar,Flags>::asCholmodMatrix()
|
||||
|
||||
ei_cholmod_configure_matrix<Scalar>(res);
|
||||
|
||||
if (Flags & SelfAdjoint)
|
||||
if (Derived::Flags & SelfAdjoint)
|
||||
{
|
||||
if (Flags & UpperTriangular)
|
||||
if (Derived::Flags & UpperTriangular)
|
||||
res.stype = 1;
|
||||
else if (Flags & LowerTriangular)
|
||||
else if (Derived::Flags & LowerTriangular)
|
||||
res.stype = -1;
|
||||
else
|
||||
res.stype = 0;
|
||||
@ -108,14 +109,14 @@ cholmod_dense ei_cholmod_map_eigen_to_dense(MatrixBase<Derived>& mat)
|
||||
}
|
||||
|
||||
template<typename Scalar, int Flags>
|
||||
MappedSparseMatrix<Scalar,Flags>::MappedSparseMatrix(taucs_ccs_matrix& taucsMat)
|
||||
MappedSparseMatrix<Scalar,Flags>::MappedSparseMatrix(cholmod_sparse& cm)
|
||||
{
|
||||
m_innerSize = cm.nrow;
|
||||
m_outerSize = cm.ncol;
|
||||
m_outerIndex = reinterpret_cast<int*>(cm.p);
|
||||
m_innerIndices = reinterpret_cast<int*>(cm.i);
|
||||
m_values = reinterpret_cast<Scalar*>(cm.x);
|
||||
m_nnz = res.m_outerIndex[cm.ncol]);
|
||||
m_nnz = m_outerIndex[cm.ncol];
|
||||
}
|
||||
|
||||
template<typename MatrixType>
|
||||
@ -123,8 +124,8 @@ class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType>
|
||||
{
|
||||
protected:
|
||||
typedef SparseLLT<MatrixType> Base;
|
||||
using typename Base::Scalar;
|
||||
using Base::RealScalar;
|
||||
typedef typename Base::Scalar Scalar;
|
||||
typedef typename Base::RealScalar RealScalar;
|
||||
using Base::MatrixLIsDirty;
|
||||
using Base::SupernodalFactorIsDirty;
|
||||
using Base::m_flags;
|
||||
@ -205,7 +206,7 @@ SparseLLT<MatrixType,Cholmod>::matrixL() const
|
||||
ei_assert(!(m_status & SupernodalFactorIsDirty));
|
||||
|
||||
cholmod_sparse* cmRes = cholmod_factor_to_sparse(m_cholmodFactor, &m_cholmod);
|
||||
const_cast<typename Base::CholMatrixType&>(m_matrix) = Base::CholMatrixType::Map(*cmRes);
|
||||
const_cast<typename Base::CholMatrixType&>(m_matrix) = MappedSparseMatrix<Scalar>(*cmRes);
|
||||
free(cmRes);
|
||||
|
||||
m_status = (m_status & ~MatrixLIsDirty);
|
||||
|
||||
@ -1,7 +1,7 @@
|
||||
// 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>
|
||||
// Copyright (C) 2008-2009 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
|
||||
@ -84,6 +84,9 @@ class DynamicSparseMatrix
|
||||
inline int outerSize() const { return m_data.size(); }
|
||||
inline int innerNonZeros(int j) const { return m_data[j].size(); }
|
||||
|
||||
std::vector<CompressedStorage<Scalar> >& _data() { return m_data; }
|
||||
const std::vector<CompressedStorage<Scalar> >& _data() const { return m_data; }
|
||||
|
||||
/** \returns the coefficient value at given position \a row, \a col
|
||||
* This operation involes a log(rho*outer_size) binary search.
|
||||
*/
|
||||
@ -124,6 +127,8 @@ class DynamicSparseMatrix
|
||||
|
||||
/** Set the matrix to zero and reserve the memory for \a reserveSize nonzero coefficients. */
|
||||
inline void startFill(int reserveSize = 1000)
|
||||
{
|
||||
if (outerSize()>0)
|
||||
{
|
||||
int reserveSizePerVector = std::max(reserveSize/outerSize(),4);
|
||||
for (int j=0; j<outerSize(); ++j)
|
||||
@ -132,6 +137,7 @@ class DynamicSparseMatrix
|
||||
m_data[j].reserve(reserveSizePerVector);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/** inserts a nonzero coefficient at given coordinates \a row, \a col and returns its reference assuming that:
|
||||
* 1 - the coefficient does not exist yet
|
||||
@ -215,7 +221,7 @@ class DynamicSparseMatrix
|
||||
}
|
||||
|
||||
inline DynamicSparseMatrix()
|
||||
: m_innerSize(0)
|
||||
: m_innerSize(0), m_data(0)
|
||||
{
|
||||
ei_assert(innerSize()==0 && outerSize()==0);
|
||||
}
|
||||
|
||||
@ -65,10 +65,10 @@ class MappedSparseMatrix
|
||||
|
||||
//----------------------------------------
|
||||
// direct access interface
|
||||
inline const Scalar* _valuePtr() const { return &m_values; }
|
||||
inline Scalar* _valuePtr() { return &m_values; }
|
||||
inline const Scalar* _valuePtr() const { return m_values; }
|
||||
inline Scalar* _valuePtr() { return m_values; }
|
||||
|
||||
inline const int* _innerIndexPtr() const { return &m_innerIndices; }
|
||||
inline const int* _innerIndexPtr() const { return m_innerIndices; }
|
||||
inline int* _innerIndexPtr() { return m_innerIndices; }
|
||||
|
||||
inline const int* _outerIndexPtr() const { return m_outerIndex; }
|
||||
@ -140,21 +140,25 @@ class MappedSparseMatrix<Scalar,_Flags>::InnerIterator
|
||||
{
|
||||
public:
|
||||
InnerIterator(const MappedSparseMatrix& mat, int outer)
|
||||
: m_matrix(mat), m_outer(outer), m_id(mat._outerIndexPtr[outer]), m_start(m_id), m_end(mat._outerIndexPtr[outer+1])
|
||||
: m_matrix(mat),
|
||||
m_outer(outer),
|
||||
m_id(mat._outerIndexPtr()[outer]),
|
||||
m_start(m_id),
|
||||
m_end(mat._outerIndexPtr()[outer+1])
|
||||
{}
|
||||
|
||||
template<unsigned int Added, unsigned int Removed>
|
||||
InnerIterator(const Flagged<MappedSparseMatrix,Added,Removed>& mat, int outer)
|
||||
: m_matrix(mat._expression()), m_id(m_matrix._outerIndexPtr[outer]),
|
||||
m_start(m_id), m_end(m_matrix._outerIndexPtr[outer+1])
|
||||
: m_matrix(mat._expression()), m_id(m_matrix._outerIndexPtr()[outer]),
|
||||
m_start(m_id), m_end(m_matrix._outerIndexPtr()[outer+1])
|
||||
{}
|
||||
|
||||
inline InnerIterator& operator++() { m_id++; return *this; }
|
||||
|
||||
inline Scalar value() const { return m_matrix.m_valuePtr[m_id]; }
|
||||
inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix._valuePtr[m_id]); }
|
||||
inline Scalar value() const { return m_matrix._valuePtr()[m_id]; }
|
||||
inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix._valuePtr()[m_id]); }
|
||||
|
||||
inline int index() const { return m_matrix._innerIndexPtr(m_id); }
|
||||
inline int index() const { return m_matrix._innerIndexPtr()[m_id]; }
|
||||
inline int row() const { return IsRowMajor ? m_outer : index(); }
|
||||
inline int col() const { return IsRowMajor ? index() : m_outer; }
|
||||
|
||||
|
||||
@ -26,61 +26,237 @@
|
||||
#ifndef EIGEN_SPARSE_BLOCK_H
|
||||
#define EIGEN_SPARSE_BLOCK_H
|
||||
|
||||
template<typename MatrixType>
|
||||
struct ei_traits<SparseInnerVector<MatrixType> >
|
||||
template<typename MatrixType, int Size>
|
||||
struct ei_traits<SparseInnerVectorSet<MatrixType, Size> >
|
||||
{
|
||||
typedef typename ei_traits<MatrixType>::Scalar Scalar;
|
||||
enum {
|
||||
IsRowMajor = (int(MatrixType::Flags)&RowMajorBit)==RowMajorBit,
|
||||
Flags = MatrixType::Flags,
|
||||
RowsAtCompileTime = IsRowMajor ? 1 : MatrixType::RowsAtCompileTime,
|
||||
ColsAtCompileTime = IsRowMajor ? MatrixType::ColsAtCompileTime : 1,
|
||||
RowsAtCompileTime = IsRowMajor ? Size : MatrixType::RowsAtCompileTime,
|
||||
ColsAtCompileTime = IsRowMajor ? MatrixType::ColsAtCompileTime : Size,
|
||||
CoeffReadCost = MatrixType::CoeffReadCost
|
||||
};
|
||||
};
|
||||
|
||||
template<typename MatrixType>
|
||||
class SparseInnerVector : ei_no_assignment_operator,
|
||||
public SparseMatrixBase<SparseInnerVector<MatrixType> >
|
||||
template<typename MatrixType, int Size>
|
||||
class SparseInnerVectorSet : ei_no_assignment_operator,
|
||||
public SparseMatrixBase<SparseInnerVectorSet<MatrixType, Size> >
|
||||
{
|
||||
enum {
|
||||
IsRowMajor = ei_traits<SparseInnerVector>::IsRowMajor
|
||||
};
|
||||
enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor };
|
||||
public:
|
||||
|
||||
EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVector)
|
||||
class InnerIterator;
|
||||
|
||||
inline SparseInnerVector(const MatrixType& matrix, int outer)
|
||||
: m_matrix(matrix), m_outer(outer)
|
||||
EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet)
|
||||
class InnerIterator: public MatrixType::InnerIterator
|
||||
{
|
||||
public:
|
||||
inline InnerIterator(const SparseInnerVectorSet& xpr, int outer)
|
||||
: MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer)
|
||||
{}
|
||||
};
|
||||
|
||||
inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize)
|
||||
: m_matrix(matrix), m_outerStart(outerStart), m_outerSize(outerSize)
|
||||
{
|
||||
ei_assert( (outerStart>=0) && ((outerStart+outerSize)<=matrix.outerSize()) );
|
||||
}
|
||||
|
||||
inline SparseInnerVectorSet(const MatrixType& matrix, int outer)
|
||||
: m_matrix(matrix), m_outerStart(outer), m_outerSize(Size)
|
||||
{
|
||||
ei_assert(Size!=Dynamic);
|
||||
ei_assert( (outer>=0) && (outer<matrix.outerSize()) );
|
||||
}
|
||||
|
||||
EIGEN_STRONG_INLINE int rows() const { return IsRowMajor ? 1 : m_matrix.rows(); }
|
||||
EIGEN_STRONG_INLINE int cols() const { return IsRowMajor ? m_matrix.cols() : 1; }
|
||||
// template<typename OtherDerived>
|
||||
// inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
|
||||
// {
|
||||
// return *this;
|
||||
// }
|
||||
|
||||
// template<typename Sparse>
|
||||
// inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
|
||||
// {
|
||||
// return *this;
|
||||
// }
|
||||
|
||||
EIGEN_STRONG_INLINE int rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); }
|
||||
EIGEN_STRONG_INLINE int cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); }
|
||||
|
||||
protected:
|
||||
|
||||
const typename MatrixType::Nested m_matrix;
|
||||
int m_outer;
|
||||
int m_outerStart;
|
||||
const ei_int_if_dynamic<Size> m_outerSize;
|
||||
|
||||
};
|
||||
|
||||
template<typename MatrixType>
|
||||
class SparseInnerVector<MatrixType>::InnerIterator : public MatrixType::InnerIterator
|
||||
/***************************************************************************
|
||||
* specialisation for DynamicSparseMatrix
|
||||
***************************************************************************/
|
||||
|
||||
template<typename _Scalar, int _Options, int Size>
|
||||
class SparseInnerVectorSet<DynamicSparseMatrix<_Scalar, _Options>, Size>
|
||||
: public SparseMatrixBase<SparseInnerVectorSet<DynamicSparseMatrix<_Scalar, _Options>, Size> >
|
||||
{
|
||||
typedef DynamicSparseMatrix<_Scalar, _Options> MatrixType;
|
||||
enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor };
|
||||
public:
|
||||
|
||||
EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet)
|
||||
class InnerIterator: public MatrixType::InnerIterator
|
||||
{
|
||||
public:
|
||||
inline InnerIterator(const SparseInnerVector& xpr, int outer=0)
|
||||
: MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outer)
|
||||
{
|
||||
ei_assert(outer==0);
|
||||
}
|
||||
inline InnerIterator(const SparseInnerVectorSet& xpr, int outer)
|
||||
: MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer)
|
||||
{}
|
||||
};
|
||||
|
||||
inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize)
|
||||
: m_matrix(matrix), m_outerStart(outerStart), m_outerSize(outerSize)
|
||||
{
|
||||
ei_assert( (outerStart>=0) && ((outerStart+outerSize)<=matrix.outerSize()) );
|
||||
}
|
||||
|
||||
inline SparseInnerVectorSet(const MatrixType& matrix, int outer)
|
||||
: m_matrix(matrix), m_outerStart(outer), m_outerSize(Size)
|
||||
{
|
||||
ei_assert(Size!=Dynamic);
|
||||
ei_assert( (outer>=0) && (outer<matrix.outerSize()) );
|
||||
}
|
||||
|
||||
template<typename OtherDerived>
|
||||
inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
|
||||
{
|
||||
if (IsRowMajor != ((OtherDerived::Flags&RowMajorBit)==RowMajorBit))
|
||||
{
|
||||
// need to transpose => perform a block evaluation followed by a big swap
|
||||
DynamicSparseMatrix<Scalar,IsRowMajor?RowMajorBit:0> aux(other);
|
||||
*this = aux.markAsRValue();
|
||||
}
|
||||
else
|
||||
{
|
||||
// evaluate/copy vector per vector
|
||||
for (int j=0; j<m_outerSize.value(); ++j)
|
||||
{
|
||||
SparseVector<Scalar,IsRowMajor ? RowMajorBit : 0> aux(other.innerVector(j));
|
||||
m_matrix.const_cast_derived()._data()[m_outerStart+j].swap(aux._data());
|
||||
}
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
inline SparseInnerVectorSet& operator=(const SparseInnerVectorSet& other)
|
||||
{
|
||||
return operator=<SparseInnerVectorSet>(other);
|
||||
}
|
||||
|
||||
// template<typename Sparse>
|
||||
// inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
|
||||
// {
|
||||
// return *this;
|
||||
// }
|
||||
|
||||
EIGEN_STRONG_INLINE int rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); }
|
||||
EIGEN_STRONG_INLINE int cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); }
|
||||
|
||||
protected:
|
||||
|
||||
const typename MatrixType::Nested m_matrix;
|
||||
int m_outerStart;
|
||||
const ei_int_if_dynamic<Size> m_outerSize;
|
||||
|
||||
};
|
||||
|
||||
|
||||
/***************************************************************************
|
||||
* specialisation for SparseMatrix
|
||||
***************************************************************************/
|
||||
/*
|
||||
template<typename _Scalar, int _Options, int Size>
|
||||
class SparseInnerVectorSet<SparseMatrix<_Scalar, _Options>, Size>
|
||||
: public SparseMatrixBase<SparseInnerVectorSet<SparseMatrix<_Scalar, _Options>, Size> >
|
||||
{
|
||||
typedef DynamicSparseMatrix<_Scalar, _Options> MatrixType;
|
||||
enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor };
|
||||
public:
|
||||
|
||||
EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet)
|
||||
class InnerIterator: public MatrixType::InnerIterator
|
||||
{
|
||||
public:
|
||||
inline InnerIterator(const SparseInnerVectorSet& xpr, int outer)
|
||||
: MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer)
|
||||
{}
|
||||
};
|
||||
|
||||
inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize)
|
||||
: m_matrix(matrix), m_outerStart(outerStart), m_outerSize(outerSize)
|
||||
{
|
||||
ei_assert( (outerStart>=0) && ((outerStart+outerSize)<=matrix.outerSize()) );
|
||||
}
|
||||
|
||||
inline SparseInnerVectorSet(const MatrixType& matrix, int outer)
|
||||
: m_matrix(matrix), m_outerStart(outer)
|
||||
{
|
||||
ei_assert(Size==1);
|
||||
ei_assert( (outer>=0) && (outer<matrix.outerSize()) );
|
||||
}
|
||||
|
||||
template<typename OtherDerived>
|
||||
inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
|
||||
{
|
||||
if (IsRowMajor != ((OtherDerived::Flags&RowMajorBit)==RowMajorBit))
|
||||
{
|
||||
// need to transpose => perform a block evaluation followed by a big swap
|
||||
DynamicSparseMatrix<Scalar,IsRowMajor?RowMajorBit:0> aux(other);
|
||||
*this = aux.markAsRValue();
|
||||
}
|
||||
else
|
||||
{
|
||||
// evaluate/copy vector per vector
|
||||
for (int j=0; j<m_outerSize.value(); ++j)
|
||||
{
|
||||
SparseVector<Scalar,IsRowMajor ? RowMajorBit : 0> aux(other.innerVector(j));
|
||||
m_matrix.const_cast_derived()._data()[m_outerStart+j].swap(aux._data());
|
||||
}
|
||||
}
|
||||
return *this;
|
||||
}
|
||||
|
||||
inline SparseInnerVectorSet& operator=(const SparseInnerVectorSet& other)
|
||||
{
|
||||
return operator=<SparseInnerVectorSet>(other);
|
||||
}
|
||||
|
||||
inline const Scalar* _valuePtr() const
|
||||
{ return m_matrix._valuePtr() + m_matrix._outerIndexPtr()[m_outerStart]; }
|
||||
inline const int* _innerIndexPtr() const
|
||||
{ return m_matrix._innerIndexPtr() + m_matrix._outerIndexPtr()[m_outerStart]; }
|
||||
inline const int* _outerIndexPtr() const { return m_matrix._outerIndexPtr() + m_outerStart; }
|
||||
|
||||
// template<typename Sparse>
|
||||
// inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
|
||||
// {
|
||||
// return *this;
|
||||
// }
|
||||
|
||||
EIGEN_STRONG_INLINE int rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); }
|
||||
EIGEN_STRONG_INLINE int cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); }
|
||||
|
||||
protected:
|
||||
|
||||
const typename MatrixType::Nested m_matrix;
|
||||
int m_outerStart;
|
||||
const ei_int_if_dynamic<Size> m_outerSize;
|
||||
|
||||
};
|
||||
*/
|
||||
//----------
|
||||
|
||||
/** \returns the i-th row of the matrix \c *this. For row-major matrix only. */
|
||||
template<typename Derived>
|
||||
SparseInnerVector<Derived> SparseMatrixBase<Derived>::row(int i)
|
||||
SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::row(int i)
|
||||
{
|
||||
EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
|
||||
return innerVector(i);
|
||||
@ -89,7 +265,7 @@ SparseInnerVector<Derived> SparseMatrixBase<Derived>::row(int i)
|
||||
/** \returns the i-th row of the matrix \c *this. For row-major matrix only.
|
||||
* (read-only version) */
|
||||
template<typename Derived>
|
||||
const SparseInnerVector<Derived> SparseMatrixBase<Derived>::row(int i) const
|
||||
const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::row(int i) const
|
||||
{
|
||||
EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
|
||||
return innerVector(i);
|
||||
@ -97,18 +273,18 @@ const SparseInnerVector<Derived> SparseMatrixBase<Derived>::row(int i) const
|
||||
|
||||
/** \returns the i-th column of the matrix \c *this. For column-major matrix only. */
|
||||
template<typename Derived>
|
||||
SparseInnerVector<Derived> SparseMatrixBase<Derived>::col(int i)
|
||||
SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::col(int i)
|
||||
{
|
||||
EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
|
||||
EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
|
||||
return innerVector(i);
|
||||
}
|
||||
|
||||
/** \returns the i-th column of the matrix \c *this. For column-major matrix only.
|
||||
* (read-only version) */
|
||||
template<typename Derived>
|
||||
const SparseInnerVector<Derived> SparseMatrixBase<Derived>::col(int i) const
|
||||
const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::col(int i) const
|
||||
{
|
||||
EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
|
||||
EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
|
||||
return innerVector(i);
|
||||
}
|
||||
|
||||
@ -116,15 +292,65 @@ const SparseInnerVector<Derived> SparseMatrixBase<Derived>::col(int i) const
|
||||
* is col-major (resp. row-major).
|
||||
*/
|
||||
template<typename Derived>
|
||||
SparseInnerVector<Derived> SparseMatrixBase<Derived>::innerVector(int outer)
|
||||
{ return SparseInnerVector<Derived>(derived(), outer); }
|
||||
SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::innerVector(int outer)
|
||||
{ return SparseInnerVectorSet<Derived,1>(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<typename Derived>
|
||||
const SparseInnerVector<Derived> SparseMatrixBase<Derived>::innerVector(int outer) const
|
||||
{ return SparseInnerVector<Derived>(derived(), outer); }
|
||||
const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::innerVector(int outer) const
|
||||
{ return SparseInnerVectorSet<Derived,1>(derived(), outer); }
|
||||
|
||||
//----------
|
||||
|
||||
/** \returns the i-th row of the matrix \c *this. For row-major matrix only. */
|
||||
template<typename Derived>
|
||||
SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::subrows(int start, int size)
|
||||
{
|
||||
EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
|
||||
return innerVectors(start, size);
|
||||
}
|
||||
|
||||
/** \returns the i-th row of the matrix \c *this. For row-major matrix only.
|
||||
* (read-only version) */
|
||||
template<typename Derived>
|
||||
const SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::subrows(int start, int size) const
|
||||
{
|
||||
EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
|
||||
return innerVectors(start, size);
|
||||
}
|
||||
|
||||
/** \returns the i-th column of the matrix \c *this. For column-major matrix only. */
|
||||
template<typename Derived>
|
||||
SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::subcols(int start, int size)
|
||||
{
|
||||
EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
|
||||
return innerVectors(start, size);
|
||||
}
|
||||
|
||||
/** \returns the i-th column of the matrix \c *this. For column-major matrix only.
|
||||
* (read-only version) */
|
||||
template<typename Derived>
|
||||
const SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::subcols(int start, int size) const
|
||||
{
|
||||
EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
|
||||
return innerVectors(start, size);
|
||||
}
|
||||
|
||||
/** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this
|
||||
* is col-major (resp. row-major).
|
||||
*/
|
||||
template<typename Derived>
|
||||
SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::innerVectors(int outerStart, int outerSize)
|
||||
{ return SparseInnerVectorSet<Derived,Dynamic>(derived(), outerStart, 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<typename Derived>
|
||||
const SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::innerVectors(int outerStart, int outerSize) const
|
||||
{ return SparseInnerVectorSet<Derived,Dynamic>(derived(), outerStart, outerSize); }
|
||||
|
||||
# if 0
|
||||
template<typename MatrixType, int BlockRows, int BlockCols, int PacketAccess>
|
||||
|
||||
@ -86,6 +86,8 @@ class SparseCwiseBinaryOp : ei_no_assignment_operator,
|
||||
EIGEN_STRONG_INLINE SparseCwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp())
|
||||
: m_lhs(lhs), m_rhs(rhs), m_functor(func)
|
||||
{
|
||||
EIGEN_STATIC_ASSERT((_LhsNested::Flags&RowMajorBit)==(_RhsNested::Flags&RowMajorBit),
|
||||
BOTH_MATRICES_MUST_HAVE_THE_SAME_STORAGE_ORDER)
|
||||
EIGEN_STATIC_ASSERT((ei_functor_allows_mixing_real_and_complex<BinaryOp>::ret
|
||||
? int(ei_is_same_type<typename Lhs::RealScalar, typename Rhs::RealScalar>::ret)
|
||||
: int(ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret)),
|
||||
@ -130,11 +132,10 @@ class SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator
|
||||
* Implementation of inner-iterators
|
||||
***************************************************************************/
|
||||
|
||||
// template<typename T> struct ei_is_scalar_product { enum { ret = false }; };
|
||||
// template<typename T> struct ei_is_scalar_product<ei_scalar_product_op<T> > { enum { ret = true }; };
|
||||
|
||||
// helper class
|
||||
// template<typename T> struct ei_func_is_conjunction { enum { ret = false }; };
|
||||
// template<typename T> struct ei_func_is_conjunction<ei_scalar_product_op<T> > { enum { ret = true }; };
|
||||
|
||||
// TODO generalize the ei_scalar_product_op specialization to all conjunctions if any !
|
||||
|
||||
// sparse - sparse (generic)
|
||||
template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived>
|
||||
@ -259,12 +260,13 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>,
|
||||
typedef SparseCwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr;
|
||||
typedef typename CwiseBinaryXpr::Scalar Scalar;
|
||||
typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
|
||||
typedef typename ei_traits<CwiseBinaryXpr>::RhsNested RhsNested;
|
||||
typedef typename _LhsNested::InnerIterator LhsIterator;
|
||||
enum { IsRowMajor = (int(Lhs::Flags)&RowMajorBit)==RowMajorBit };
|
||||
public:
|
||||
|
||||
EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
|
||||
: m_xpr(xpr), m_lhsIter(xpr.lhs(),outer), m_functor(xpr.functor()), m_outer(outer)
|
||||
: m_rhs(xpr.rhs()), m_lhsIter(xpr.lhs(),outer), m_functor(xpr.functor()), m_outer(outer)
|
||||
{}
|
||||
|
||||
EIGEN_STRONG_INLINE Derived& operator++()
|
||||
@ -275,7 +277,7 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>,
|
||||
|
||||
EIGEN_STRONG_INLINE Scalar value() const
|
||||
{ return m_functor(m_lhsIter.value(),
|
||||
m_xpr.rhs().coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); }
|
||||
m_rhs.coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); }
|
||||
|
||||
EIGEN_STRONG_INLINE int index() const { return m_lhsIter.index(); }
|
||||
EIGEN_STRONG_INLINE int row() const { return m_lhsIter.row(); }
|
||||
@ -284,9 +286,9 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>,
|
||||
EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter; }
|
||||
|
||||
protected:
|
||||
const CwiseBinaryXpr& m_xpr;
|
||||
const RhsNested m_rhs;
|
||||
LhsIterator m_lhsIter;
|
||||
const BinaryFunc& m_functor;
|
||||
const BinaryFunc m_functor;
|
||||
const int m_outer;
|
||||
};
|
||||
|
||||
@ -329,6 +331,10 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>,
|
||||
};
|
||||
|
||||
|
||||
/***************************************************************************
|
||||
* Implementation of SparseMatrixBase and SparseCwise functions/operators
|
||||
***************************************************************************/
|
||||
|
||||
template<typename Derived>
|
||||
template<typename OtherDerived>
|
||||
EIGEN_STRONG_INLINE const SparseCwiseBinaryOp<ei_scalar_difference_op<typename ei_traits<Derived>::Scalar>,
|
||||
|
||||
@ -89,7 +89,7 @@ class SparseCwiseUnaryOp<UnaryOp,MatrixType>::InnerIterator
|
||||
|
||||
protected:
|
||||
MatrixTypeIterator m_iter;
|
||||
const UnaryOp& m_functor;
|
||||
const UnaryOp m_functor;
|
||||
};
|
||||
|
||||
template<typename Derived>
|
||||
|
||||
157
Eigen/src/Sparse/SparseDiagonalProduct.h
Normal file
157
Eigen/src/Sparse/SparseDiagonalProduct.h
Normal file
@ -0,0 +1,157 @@
|
||||
// This file is part of Eigen, a lightweight C++ template library
|
||||
// for linear algebra. Eigen itself is part of the KDE project.
|
||||
//
|
||||
// Copyright (C) 2009 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_SPARSE_DIAGONAL_PRODUCT_H
|
||||
#define EIGEN_SPARSE_DIAGONAL_PRODUCT_H
|
||||
|
||||
// the product a diagonal matrix with a sparse matrix can be easily
|
||||
// implemented using expression template. We have two 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.
|
||||
|
||||
template<typename Lhs, typename Rhs>
|
||||
struct ei_traits<SparseDiagonalProduct<Lhs, Rhs> > : ei_traits<SparseProduct<Lhs, Rhs, DiagonalProduct> >
|
||||
{
|
||||
typedef typename ei_cleantype<Lhs>::type _Lhs;
|
||||
typedef typename ei_cleantype<Rhs>::type _Rhs;
|
||||
enum {
|
||||
SparseFlags = ((int(_Lhs::Flags)&Diagonal)==Diagonal) ? int(_Rhs::Flags) : int(_Lhs::Flags),
|
||||
Flags = SparseBit | (SparseFlags&RowMajorBit)
|
||||
};
|
||||
};
|
||||
|
||||
enum {SDP_IsDiagonal, SDP_IsSparseRowMajor, SDP_IsSparseColMajor};
|
||||
template<typename Lhs, typename Rhs, typename SparseDiagonalProductType, int RhsMode, int LhsMode>
|
||||
class ei_sparse_diagonal_product_inner_iterator_selector;
|
||||
|
||||
template<typename LhsNested, typename RhsNested>
|
||||
class SparseDiagonalProduct : public SparseMatrixBase<SparseDiagonalProduct<LhsNested,RhsNested> >, ei_no_assignment_operator
|
||||
{
|
||||
typedef typename ei_traits<SparseDiagonalProduct>::_LhsNested _LhsNested;
|
||||
typedef typename ei_traits<SparseDiagonalProduct>::_RhsNested _RhsNested;
|
||||
|
||||
enum {
|
||||
LhsMode = (_LhsNested::Flags&Diagonal)==Diagonal ? SDP_IsDiagonal
|
||||
: (_LhsNested::Flags&RowMajorBit) ? SDP_IsSparseRowMajor : SDP_IsSparseColMajor,
|
||||
RhsMode = (_RhsNested::Flags&Diagonal)==Diagonal ? SDP_IsDiagonal
|
||||
: (_RhsNested::Flags&RowMajorBit) ? SDP_IsSparseRowMajor : SDP_IsSparseColMajor
|
||||
};
|
||||
|
||||
public:
|
||||
|
||||
EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseDiagonalProduct)
|
||||
|
||||
typedef ei_sparse_diagonal_product_inner_iterator_selector
|
||||
<_LhsNested,_RhsNested,SparseDiagonalProduct,LhsMode,RhsMode> InnerIterator;
|
||||
|
||||
template<typename Lhs, typename Rhs>
|
||||
EIGEN_STRONG_INLINE SparseDiagonalProduct(const Lhs& lhs, const Rhs& rhs)
|
||||
: m_lhs(lhs), m_rhs(rhs)
|
||||
{
|
||||
ei_assert(lhs.cols() == rhs.rows() && "invalid sparse matrix * diagonal matrix product");
|
||||
}
|
||||
|
||||
EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); }
|
||||
EIGEN_STRONG_INLINE int 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;
|
||||
};
|
||||
|
||||
|
||||
template<typename Lhs, typename Rhs, typename SparseDiagonalProductType>
|
||||
class ei_sparse_diagonal_product_inner_iterator_selector
|
||||
<Lhs,Rhs,SparseDiagonalProductType,SDP_IsDiagonal,SDP_IsSparseRowMajor>
|
||||
: public SparseCwiseUnaryOp<ei_scalar_multiple_op<typename Lhs::Scalar>,Rhs>::InnerIterator
|
||||
{
|
||||
typedef typename SparseCwiseUnaryOp<ei_scalar_multiple_op<typename Lhs::Scalar>,Rhs>::InnerIterator Base;
|
||||
public:
|
||||
inline ei_sparse_diagonal_product_inner_iterator_selector(
|
||||
const SparseDiagonalProductType& expr, int outer)
|
||||
: Base(expr.rhs()*(expr.lhs().diagonal().coeff(outer)), outer)
|
||||
{}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename SparseDiagonalProductType>
|
||||
class ei_sparse_diagonal_product_inner_iterator_selector
|
||||
<Lhs,Rhs,SparseDiagonalProductType,SDP_IsDiagonal,SDP_IsSparseColMajor>
|
||||
: public SparseCwiseBinaryOp<
|
||||
ei_scalar_product_op<typename Lhs::Scalar>,
|
||||
SparseInnerVectorSet<Rhs,1>,
|
||||
typename Lhs::_CoeffsVectorType>::InnerIterator
|
||||
{
|
||||
typedef typename SparseCwiseBinaryOp<
|
||||
ei_scalar_product_op<typename Lhs::Scalar>,
|
||||
SparseInnerVectorSet<Rhs,1>,
|
||||
typename Lhs::_CoeffsVectorType>::InnerIterator Base;
|
||||
public:
|
||||
inline ei_sparse_diagonal_product_inner_iterator_selector(
|
||||
const SparseDiagonalProductType& expr, int outer)
|
||||
: Base(expr.rhs().innerVector(outer) .cwise()* expr.lhs().diagonal(), 0)
|
||||
{}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename SparseDiagonalProductType>
|
||||
class ei_sparse_diagonal_product_inner_iterator_selector
|
||||
<Lhs,Rhs,SparseDiagonalProductType,SDP_IsSparseColMajor,SDP_IsDiagonal>
|
||||
: public SparseCwiseUnaryOp<ei_scalar_multiple_op<typename Rhs::Scalar>,Lhs>::InnerIterator
|
||||
{
|
||||
typedef typename SparseCwiseUnaryOp<ei_scalar_multiple_op<typename Rhs::Scalar>,Lhs>::InnerIterator Base;
|
||||
public:
|
||||
inline ei_sparse_diagonal_product_inner_iterator_selector(
|
||||
const SparseDiagonalProductType& expr, int outer)
|
||||
: Base(expr.lhs()*expr.rhs().diagonal().coeff(outer), outer)
|
||||
{}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename SparseDiagonalProductType>
|
||||
class ei_sparse_diagonal_product_inner_iterator_selector
|
||||
<Lhs,Rhs,SparseDiagonalProductType,SDP_IsSparseRowMajor,SDP_IsDiagonal>
|
||||
: public SparseCwiseBinaryOp<
|
||||
ei_scalar_product_op<typename Rhs::Scalar>,
|
||||
SparseInnerVectorSet<Lhs,1>,
|
||||
NestByValue<Transpose<typename Rhs::_CoeffsVectorType> > >::InnerIterator
|
||||
{
|
||||
typedef typename SparseCwiseBinaryOp<
|
||||
ei_scalar_product_op<typename Rhs::Scalar>,
|
||||
SparseInnerVectorSet<Lhs,1>,
|
||||
NestByValue<Transpose<typename Rhs::_CoeffsVectorType> > >::InnerIterator Base;
|
||||
public:
|
||||
inline ei_sparse_diagonal_product_inner_iterator_selector(
|
||||
const SparseDiagonalProductType& expr, int outer)
|
||||
: Base(expr.lhs().innerVector(outer) .cwise()* expr.rhs().diagonal().transpose().nestByValue(), 0)
|
||||
{}
|
||||
};
|
||||
|
||||
#endif // EIGEN_SPARSE_DIAGONAL_PRODUCT_H
|
||||
@ -1,7 +1,7 @@
|
||||
// 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>
|
||||
// Copyright (C) 2008-2009 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
|
||||
@ -138,7 +138,6 @@ class SparseMatrix
|
||||
*/
|
||||
inline void startFill(int reserveSize = 1000)
|
||||
{
|
||||
// std::cerr << this << " startFill\n";
|
||||
setZero();
|
||||
m_data.reserve(reserveSize);
|
||||
}
|
||||
@ -161,6 +160,10 @@ class SparseMatrix
|
||||
}
|
||||
m_outerIndex[outer+1] = m_outerIndex[outer];
|
||||
}
|
||||
else
|
||||
{
|
||||
ei_assert(m_data.index(m_data.size()-1)<inner && "wrong sorted insertion");
|
||||
}
|
||||
assert(size_t(m_outerIndex[outer+1]) == m_data.size());
|
||||
int id = m_outerIndex[outer+1];
|
||||
++m_outerIndex[outer+1];
|
||||
@ -390,11 +393,11 @@ class SparseMatrix
|
||||
s << std::endl;
|
||||
s << std::endl;
|
||||
s << "Column pointers:\n";
|
||||
for (int i=0; i<m.cols(); ++i)
|
||||
for (int i=0; i<m.outerSize(); ++i)
|
||||
{
|
||||
s << m.m_outerIndex[i] << " ";
|
||||
}
|
||||
s << std::endl;
|
||||
s << " $" << std::endl;
|
||||
s << std::endl;
|
||||
);
|
||||
s << static_cast<const SparseMatrixBase<SparseMatrix>&>(m);
|
||||
|
||||
@ -327,18 +327,21 @@ template<typename Derived> class SparseMatrixBase
|
||||
// void transposeInPlace();
|
||||
const AdjointReturnType adjoint() const { return conjugate()/*.nestByValue()*/; }
|
||||
|
||||
SparseInnerVector<Derived> row(int i);
|
||||
const SparseInnerVector<Derived> row(int i) const;
|
||||
SparseInnerVector<Derived> col(int j);
|
||||
const SparseInnerVector<Derived> col(int j) const;
|
||||
SparseInnerVector<Derived> innerVector(int outer);
|
||||
const SparseInnerVector<Derived> innerVector(int outer) const;
|
||||
// sub-vector
|
||||
SparseInnerVectorSet<Derived,1> row(int i);
|
||||
const SparseInnerVectorSet<Derived,1> row(int i) const;
|
||||
SparseInnerVectorSet<Derived,1> col(int j);
|
||||
const SparseInnerVectorSet<Derived,1> col(int j) const;
|
||||
SparseInnerVectorSet<Derived,1> innerVector(int outer);
|
||||
const SparseInnerVectorSet<Derived,1> innerVector(int outer) const;
|
||||
|
||||
// RowXpr row(int i);
|
||||
// const RowXpr row(int i) const;
|
||||
|
||||
// ColXpr col(int i);
|
||||
// const ColXpr col(int i) const;
|
||||
// set of sub-vectors
|
||||
SparseInnerVectorSet<Derived,Dynamic> subrows(int start, int size);
|
||||
const SparseInnerVectorSet<Derived,Dynamic> subrows(int start, int size) const;
|
||||
SparseInnerVectorSet<Derived,Dynamic> subcols(int start, int size);
|
||||
const SparseInnerVectorSet<Derived,Dynamic> subcols(int start, int size) const;
|
||||
SparseInnerVectorSet<Derived,Dynamic> innerVectors(int outerStart, int outerSize);
|
||||
const SparseInnerVectorSet<Derived,Dynamic> innerVectors(int outerStart, int outerSize) const;
|
||||
|
||||
// typename BlockReturnType<Derived>::Type block(int startRow, int startCol, int blockRows, int blockCols);
|
||||
// const typename BlockReturnType<Derived>::Type
|
||||
|
||||
@ -29,7 +29,9 @@ template<typename Lhs, typename Rhs> struct ei_sparse_product_mode
|
||||
{
|
||||
enum {
|
||||
|
||||
value = (Rhs::Flags&Lhs::Flags&SparseBit)==SparseBit
|
||||
value = ((Lhs::Flags&Diagonal)==Diagonal || (Rhs::Flags&Diagonal)==Diagonal)
|
||||
? DiagonalProduct
|
||||
: (Rhs::Flags&Lhs::Flags&SparseBit)==SparseBit
|
||||
? SparseTimeSparseProduct
|
||||
: (Lhs::Flags&SparseBit)==SparseBit
|
||||
? SparseTimeDenseProduct
|
||||
@ -45,6 +47,15 @@ struct SparseProductReturnType
|
||||
typedef SparseProduct<LhsNested, RhsNested, ProductMode> Type;
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs>
|
||||
struct SparseProductReturnType<Lhs,Rhs,DiagonalProduct>
|
||||
{
|
||||
typedef const typename ei_nested<Lhs,Rhs::RowsAtCompileTime>::type LhsNested;
|
||||
typedef const typename ei_nested<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
|
||||
|
||||
typedef SparseDiagonalProduct<LhsNested, RhsNested> Type;
|
||||
};
|
||||
|
||||
// sparse product return type specialization
|
||||
template<typename Lhs, typename Rhs>
|
||||
struct SparseProductReturnType<Lhs,Rhs,SparseTimeSparseProduct>
|
||||
@ -95,7 +106,7 @@ struct ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >
|
||||
// RhsIsRowMajor = (RhsFlags & RowMajorBit)==RowMajorBit,
|
||||
|
||||
EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit),
|
||||
ResultIsSparse = ProductMode==SparseTimeSparseProduct,
|
||||
ResultIsSparse = ProductMode==SparseTimeSparseProduct || ProductMode==DiagonalProduct,
|
||||
|
||||
RemovedBits = ~( (EvalToRowMajor ? 0 : RowMajorBit) | (ResultIsSparse ? 0 : SparseBit) ),
|
||||
|
||||
@ -112,7 +123,8 @@ struct ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >
|
||||
};
|
||||
|
||||
template<typename LhsNested, typename RhsNested, int ProductMode>
|
||||
class SparseProduct : ei_no_assignment_operator, public ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >::Base
|
||||
class SparseProduct : ei_no_assignment_operator,
|
||||
public ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >::Base
|
||||
{
|
||||
public:
|
||||
|
||||
@ -159,19 +171,12 @@ class SparseProduct : ei_no_assignment_operator, public ei_traits<SparseProduct<
|
||||
RhsNested m_rhs;
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType,
|
||||
int LhsStorageOrder = ei_traits<Lhs>::Flags&RowMajorBit,
|
||||
int RhsStorageOrder = ei_traits<Rhs>::Flags&RowMajorBit,
|
||||
int ResStorageOrder = ei_traits<ResultType>::Flags&RowMajorBit>
|
||||
struct ei_sparse_product_selector;
|
||||
|
||||
// perform a pseudo in-place sparse * sparse product assuming all matrices are col major
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct ei_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
|
||||
static void ei_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typedef typename ei_traits<typename ei_cleantype<Lhs>::type>::Scalar Scalar;
|
||||
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
// make sure to call innerSize/outerSize since we fake the storage order.
|
||||
int rows = lhs.innerSize();
|
||||
int cols = rhs.outerSize();
|
||||
@ -182,7 +187,7 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
|
||||
AmbiVector<Scalar> tempVector(rows);
|
||||
|
||||
// estimate the number of non zero entries
|
||||
float ratioLhs = float(lhs.nonZeros())/float(lhs.rows()*lhs.cols());
|
||||
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);
|
||||
|
||||
@ -214,16 +219,35 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
|
||||
}
|
||||
res.endFill();
|
||||
}
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType,
|
||||
int LhsStorageOrder = ei_traits<Lhs>::Flags&RowMajorBit,
|
||||
int RhsStorageOrder = ei_traits<Rhs>::Flags&RowMajorBit,
|
||||
int ResStorageOrder = ei_traits<ResultType>::Flags&RowMajorBit>
|
||||
struct ei_sparse_product_selector;
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct ei_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
|
||||
{
|
||||
typedef typename ei_traits<typename ei_cleantype<Lhs>::type>::Scalar Scalar;
|
||||
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
typename ei_cleantype<ResultType>::type _res(res.rows(), res.cols());
|
||||
ei_sparse_product_impl<Lhs,Rhs,ResultType>(lhs, rhs, _res);
|
||||
res.swap(_res);
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct ei_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
// we need a col-major matrix to hold the result
|
||||
typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
|
||||
SparseTemporaryType _res(res.rows(), res.cols());
|
||||
ei_sparse_product_selector<Lhs,Rhs,SparseTemporaryType,ColMajor,ColMajor,ColMajor>::run(lhs, rhs, _res);
|
||||
ei_sparse_product_impl<Lhs,Rhs,SparseTemporaryType>(lhs, rhs, _res);
|
||||
res = _res;
|
||||
}
|
||||
};
|
||||
@ -234,20 +258,21 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
// let's transpose the product to get a column x column product
|
||||
ei_sparse_product_selector<Rhs,Lhs,ResultType,ColMajor,ColMajor,ColMajor>::run(rhs, lhs, res);
|
||||
typename ei_cleantype<ResultType>::type _res(res.rows(), res.cols());
|
||||
ei_sparse_product_impl<Rhs,Lhs,ResultType>(rhs, lhs, _res);
|
||||
res.swap(_res);
|
||||
}
|
||||
};
|
||||
|
||||
template<typename Lhs, typename Rhs, typename ResultType>
|
||||
struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
|
||||
{
|
||||
typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
|
||||
static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
|
||||
{
|
||||
// let's transpose the product to get a column x column product
|
||||
typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
|
||||
SparseTemporaryType _res(res.cols(), res.rows());
|
||||
ei_sparse_product_selector<Rhs,Lhs,SparseTemporaryType,ColMajor,ColMajor,ColMajor>
|
||||
::run(rhs, lhs, _res);
|
||||
ei_sparse_product_impl<Rhs,Lhs,SparseTemporaryType>(rhs, lhs, _res);
|
||||
res = _res.transpose();
|
||||
}
|
||||
};
|
||||
@ -285,7 +310,6 @@ template<typename Derived>
|
||||
template<typename Lhs, typename Rhs>
|
||||
inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs,SparseTimeSparseProduct>& product)
|
||||
{
|
||||
// std::cout << "sparse product to sparse\n";
|
||||
ei_sparse_product_selector<
|
||||
typename ei_cleantype<Lhs>::type,
|
||||
typename ei_cleantype<Rhs>::type,
|
||||
@ -333,7 +357,7 @@ Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,SparseTimeD
|
||||
derived().row(j) += i.value() * product.rhs().row(j);
|
||||
++i;
|
||||
}
|
||||
Block<Derived,1,Derived::ColsAtCompileTime> foo = derived().row(j);
|
||||
Block<Derived,1,Derived::ColsAtCompileTime> res(derived().row(LhsIsRowMajor ? j : 0));
|
||||
for (; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
|
||||
{
|
||||
if (LhsIsSelfAdjoint)
|
||||
@ -345,7 +369,7 @@ Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,SparseTimeD
|
||||
derived().row(b) += ei_conj(v) * product.rhs().row(a);
|
||||
}
|
||||
else if (LhsIsRowMajor)
|
||||
foo += i.value() * product.rhs().row(i.index());
|
||||
res += i.value() * product.rhs().row(i.index());
|
||||
else
|
||||
derived().row(i.index()) += i.value() * product.rhs().row(j);
|
||||
}
|
||||
|
||||
@ -107,12 +107,13 @@ template<typename _Scalar, int _Flags = 0> class SparseVector;
|
||||
template<typename _Scalar, int _Flags = 0> class MappedSparseMatrix;
|
||||
|
||||
template<typename MatrixType> class SparseTranspose;
|
||||
template<typename MatrixType> class SparseInnerVector;
|
||||
template<typename MatrixType, int Size> class SparseInnerVectorSet;
|
||||
template<typename Derived> class SparseCwise;
|
||||
template<typename UnaryOp, typename MatrixType> class SparseCwiseUnaryOp;
|
||||
template<typename BinaryOp, typename Lhs, typename Rhs> class SparseCwiseBinaryOp;
|
||||
template<typename ExpressionType,
|
||||
unsigned int Added, unsigned int Removed> class SparseFlagged;
|
||||
template<typename Lhs, typename Rhs> class SparseDiagonalProduct;
|
||||
|
||||
template<typename Lhs, typename Rhs> struct ei_sparse_product_mode;
|
||||
template<typename Lhs, typename Rhs, int ProductMode = ei_sparse_product_mode<Lhs,Rhs>::value> struct SparseProductReturnType;
|
||||
|
||||
@ -1,7 +1,7 @@
|
||||
// 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>
|
||||
// Copyright (C) 2008-2009 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
|
||||
@ -59,6 +59,7 @@ class SparseVector
|
||||
EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseVector)
|
||||
EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseVector, +=)
|
||||
EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseVector, -=)
|
||||
// EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseVector, =)
|
||||
|
||||
protected:
|
||||
public:
|
||||
@ -69,6 +70,9 @@ class SparseVector
|
||||
CompressedStorage<Scalar> m_data;
|
||||
int m_size;
|
||||
|
||||
CompressedStorage<Scalar>& _data() { return m_data; }
|
||||
CompressedStorage<Scalar>& _data() const { return m_data; }
|
||||
|
||||
public:
|
||||
|
||||
EIGEN_STRONG_INLINE int rows() const { return IsColVector ? m_size : 1; }
|
||||
@ -199,6 +203,13 @@ class SparseVector
|
||||
*this = other.derived();
|
||||
}
|
||||
|
||||
template<typename OtherDerived>
|
||||
inline SparseVector(const SparseMatrixBase<OtherDerived>& other)
|
||||
: m_size(0)
|
||||
{
|
||||
*this = other.derived();
|
||||
}
|
||||
|
||||
inline SparseVector(const SparseVector& other)
|
||||
: m_size(0)
|
||||
{
|
||||
@ -225,9 +236,12 @@ class SparseVector
|
||||
return *this;
|
||||
}
|
||||
|
||||
// template<typename OtherDerived>
|
||||
// inline SparseVector& operator=(const MatrixBase<OtherDerived>& other)
|
||||
// {
|
||||
template<typename OtherDerived>
|
||||
inline SparseVector& operator=(const SparseMatrixBase<OtherDerived>& other)
|
||||
{
|
||||
return Base::operator=(other);
|
||||
}
|
||||
|
||||
// const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
|
||||
// if (needToTranspose)
|
||||
// {
|
||||
|
||||
@ -1,7 +1,7 @@
|
||||
// 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>
|
||||
// Copyright (C) 2008-2009 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
|
||||
@ -35,7 +35,8 @@
|
||||
FLOATTYPE *recip_pivot_growth, \
|
||||
FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \
|
||||
SuperLUStat_t *stats, int *info, KEYTYPE) { \
|
||||
NAMESPACE::mem_usage_t mem_usage; \
|
||||
using namespace NAMESPACE; \
|
||||
mem_usage_t mem_usage; \
|
||||
NAMESPACE::FNAME(options, A, perm_c, perm_r, etree, equed, R, C, L, \
|
||||
U, work, lwork, B, X, recip_pivot_growth, rcond, \
|
||||
ferr, berr, &mem_usage, stats, info); \
|
||||
@ -59,7 +60,10 @@ struct SluMatrixMapHelper;
|
||||
*/
|
||||
struct SluMatrix : SuperMatrix
|
||||
{
|
||||
SluMatrix() {}
|
||||
SluMatrix()
|
||||
{
|
||||
Store = &storage;
|
||||
}
|
||||
|
||||
SluMatrix(const SluMatrix& other)
|
||||
: SuperMatrix(other)
|
||||
@ -68,6 +72,14 @@ struct SluMatrix : SuperMatrix
|
||||
storage = other.storage;
|
||||
}
|
||||
|
||||
SluMatrix& operator=(const SluMatrix& other)
|
||||
{
|
||||
SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
|
||||
Store = &storage;
|
||||
storage = other.storage;
|
||||
return *this;
|
||||
}
|
||||
|
||||
struct
|
||||
{
|
||||
union {int nnz;int lda;};
|
||||
@ -223,6 +235,7 @@ SluMatrix SparseMatrixBase<Derived>::asSluMatrix()
|
||||
return SluMatrix::Map(derived());
|
||||
}
|
||||
|
||||
/** View a Super LU matrix as an Eigen expression */
|
||||
template<typename Scalar, int Flags>
|
||||
MappedSparseMatrix<Scalar,Flags>::MappedSparseMatrix(SluMatrix& sluMat)
|
||||
{
|
||||
|
||||
@ -1,7 +1,7 @@
|
||||
// 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>
|
||||
// Copyright (C) 2008-2009 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
|
||||
@ -32,9 +32,9 @@ taucs_ccs_matrix SparseMatrixBase<Derived>::asTaucsMatrix()
|
||||
res.n = cols();
|
||||
res.m = rows();
|
||||
res.flags = 0;
|
||||
res.colptr = _outerIndexPtr();
|
||||
res.rowind = _innerIndexPtr();
|
||||
res.values.v = _valuePtr();
|
||||
res.colptr = derived()._outerIndexPtr();
|
||||
res.rowind = derived()._innerIndexPtr();
|
||||
res.values.v = derived()._valuePtr();
|
||||
if (ei_is_same_type<Scalar,int>::ret)
|
||||
res.flags |= TAUCS_INT;
|
||||
else if (ei_is_same_type<Scalar,float>::ret)
|
||||
@ -78,8 +78,8 @@ class SparseLLT<MatrixType,Taucs> : public SparseLLT<MatrixType>
|
||||
{
|
||||
protected:
|
||||
typedef SparseLLT<MatrixType> Base;
|
||||
using Base::Scalar;
|
||||
using Base::RealScalar;
|
||||
typedef typename Base::Scalar Scalar;
|
||||
typedef typename Base::RealScalar RealScalar;
|
||||
using Base::MatrixLIsDirty;
|
||||
using Base::SupernodalFactorIsDirty;
|
||||
using Base::m_flags;
|
||||
@ -129,7 +129,10 @@ void SparseLLT<MatrixType,Taucs>::compute(const MatrixType& a)
|
||||
{
|
||||
taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix();
|
||||
taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, Base::m_precision, 0);
|
||||
m_matrix = Base::CholMatrixType::Map(*taucsRes);
|
||||
// the matrix returned by Taucs is not necessarily sorted,
|
||||
// so let's copy it in two steps
|
||||
DynamicSparseMatrix<Scalar,RowMajor> tmp = MappedSparseMatrix<Scalar>(*taucsRes);
|
||||
m_matrix = tmp;
|
||||
free(taucsRes);
|
||||
m_status = (m_status & ~(CompleteFactorization|MatrixLIsDirty))
|
||||
| IncompleteFactorization
|
||||
@ -161,7 +164,11 @@ SparseLLT<MatrixType,Taucs>::matrixL() const
|
||||
ei_assert(!(m_status & SupernodalFactorIsDirty));
|
||||
|
||||
taucs_ccs_matrix* taucsL = taucs_supernodal_factor_to_ccs(m_taucsSupernodalFactor);
|
||||
const_cast<typename Base::CholMatrixType&>(m_matrix) = Base::CholMatrixType::Map(*taucsL);
|
||||
|
||||
// the matrix returned by Taucs is not necessarily sorted,
|
||||
// so let's copy it in two steps
|
||||
DynamicSparseMatrix<Scalar,RowMajor> tmp = MappedSparseMatrix<Scalar>(*taucsL);
|
||||
const_cast<typename Base::CholMatrixType&>(m_matrix) = tmp;
|
||||
free(taucsL);
|
||||
m_status = (m_status & ~MatrixLIsDirty);
|
||||
}
|
||||
@ -172,22 +179,32 @@ template<typename MatrixType>
|
||||
template<typename Derived>
|
||||
void SparseLLT<MatrixType,Taucs>::solveInPlace(MatrixBase<Derived> &b) const
|
||||
{
|
||||
if (m_status & MatrixLIsDirty)
|
||||
bool inputIsCompatibleWithTaucs = (Derived::Flags&RowMajorBit)==0;
|
||||
|
||||
if (!inputIsCompatibleWithTaucs)
|
||||
{
|
||||
// TODO use taucs's supernodal solver, in particular check types, storage order, etc.
|
||||
// VectorXb x(b.rows());
|
||||
// for (int j=0; j<b.cols(); ++j)
|
||||
// {
|
||||
// taucs_supernodal_solve_llt(m_taucsSupernodalFactor,x.data(),&b.col(j).coeffRef(0));
|
||||
// b.col(j) = x;
|
||||
// }
|
||||
matrixL();
|
||||
}
|
||||
|
||||
|
||||
{
|
||||
Base::solveInPlace(b);
|
||||
}
|
||||
else if (m_flags & IncompleteFactorization)
|
||||
{
|
||||
taucs_ccs_matrix taucsLLT = const_cast<typename Base::CholMatrixType&>(m_matrix).asTaucsMatrix();
|
||||
typename ei_plain_matrix_type<Derived>::type x(b.rows());
|
||||
for (int j=0; j<b.cols(); ++j)
|
||||
{
|
||||
taucs_ccs_solve_llt(&taucsLLT,x.data(),&b.col(j).coeffRef(0));
|
||||
b.col(j) = x;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
typename ei_plain_matrix_type<Derived>::type x(b.rows());
|
||||
for (int j=0; j<b.cols(); ++j)
|
||||
{
|
||||
taucs_supernodal_solve_llt(m_taucsSupernodalFactor,x.data(),&b.col(j).coeffRef(0));
|
||||
b.col(j) = x;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#endif // EIGEN_TAUCSSUPPORT_H
|
||||
|
||||
@ -1,7 +1,7 @@
|
||||
// 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>
|
||||
// Copyright (C) 2008-2009 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
|
||||
|
||||
@ -217,6 +217,7 @@ endif(QT4_FOUND)
|
||||
ei_add_test(sparse_vector)
|
||||
ei_add_test(sparse_basic)
|
||||
ei_add_test(sparse_solvers " " "${SPARSE_LIBS}")
|
||||
ei_add_test(sparse_product)
|
||||
|
||||
# print a summary of the different options
|
||||
message("************************************************************")
|
||||
|
||||
@ -193,7 +193,6 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
|
||||
}
|
||||
}
|
||||
m2.endFill();
|
||||
//std::cerr << m1 << "\n\n" << m2 << "\n";
|
||||
VERIFY_IS_APPROX(m2,m1);
|
||||
}
|
||||
|
||||
@ -239,6 +238,8 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
|
||||
VERIFY_IS_APPROX(m1+=m2, refM1+=refM2);
|
||||
VERIFY_IS_APPROX(m1-=m2, refM1-=refM2);
|
||||
|
||||
VERIFY_IS_APPROX(m1.col(0).dot(refM2.row(0)), refM1.col(0).dot(refM2.row(0)));
|
||||
|
||||
refM4.setRandom();
|
||||
// sparse cwise* dense
|
||||
VERIFY_IS_APPROX(m3.cwise()*refM4, refM3.cwise()*refM4);
|
||||
@ -254,6 +255,24 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
|
||||
int j1 = ei_random(0,rows-1);
|
||||
VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.col(j0));
|
||||
VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.col(j0)+refMat2.col(j1));
|
||||
//m2.innerVector(j0) = 2*m2.innerVector(j1);
|
||||
//refMat2.col(j0) = 2*refMat2.col(j1);
|
||||
//VERIFY_IS_APPROX(m2, refMat2);
|
||||
}
|
||||
|
||||
// test innerVectors()
|
||||
{
|
||||
DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
|
||||
SparseMatrixType m2(rows, rows);
|
||||
initSparse<Scalar>(density, refMat2, m2);
|
||||
int j0 = ei_random(0,rows-2);
|
||||
int j1 = ei_random(0,rows-2);
|
||||
int n0 = ei_random<int>(1,rows-std::max(j0,j1));
|
||||
VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(0,j0,rows,n0));
|
||||
VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0),
|
||||
refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0));
|
||||
//m2.innerVectors(j0,n0) = m2.innerVectors(j0,n0) + m2.innerVectors(j1,n0);
|
||||
//refMat2.block(0,j0,rows,n0) = refMat2.block(0,j0,rows,n0) + refMat2.block(0,j1,rows,n0);
|
||||
}
|
||||
|
||||
// test transpose
|
||||
@ -265,69 +284,6 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
|
||||
VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose());
|
||||
}
|
||||
|
||||
// test matrix product
|
||||
{
|
||||
DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
|
||||
DenseMatrix refMat3 = DenseMatrix::Zero(rows, rows);
|
||||
DenseMatrix refMat4 = DenseMatrix::Zero(rows, rows);
|
||||
DenseMatrix dm4 = DenseMatrix::Zero(rows, rows);
|
||||
SparseMatrixType m2(rows, rows);
|
||||
SparseMatrixType m3(rows, rows);
|
||||
SparseMatrixType m4(rows, rows);
|
||||
initSparse<Scalar>(density, refMat2, m2);
|
||||
initSparse<Scalar>(density, refMat3, m3);
|
||||
initSparse<Scalar>(density, refMat4, m4);
|
||||
VERIFY_IS_APPROX(m4=m2*m3, refMat4=refMat2*refMat3);
|
||||
VERIFY_IS_APPROX(m4=m2.transpose()*m3, refMat4=refMat2.transpose()*refMat3);
|
||||
VERIFY_IS_APPROX(m4=m2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
|
||||
VERIFY_IS_APPROX(m4=m2*m3.transpose(), refMat4=refMat2*refMat3.transpose());
|
||||
|
||||
// sparse * dense
|
||||
VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3);
|
||||
VERIFY_IS_APPROX(dm4=m2*refMat3.transpose(), refMat4=refMat2*refMat3.transpose());
|
||||
VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3, refMat4=refMat2.transpose()*refMat3);
|
||||
VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
|
||||
|
||||
// dense * sparse
|
||||
VERIFY_IS_APPROX(dm4=refMat2*m3, refMat4=refMat2*refMat3);
|
||||
VERIFY_IS_APPROX(dm4=refMat2*m3.transpose(), refMat4=refMat2*refMat3.transpose());
|
||||
VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3, refMat4=refMat2.transpose()*refMat3);
|
||||
VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
|
||||
}
|
||||
|
||||
// test self adjoint products
|
||||
{
|
||||
DenseMatrix b = DenseMatrix::Random(rows, rows);
|
||||
DenseMatrix x = DenseMatrix::Random(rows, rows);
|
||||
DenseMatrix refX = DenseMatrix::Random(rows, rows);
|
||||
DenseMatrix refUp = DenseMatrix::Zero(rows, rows);
|
||||
DenseMatrix refLo = DenseMatrix::Zero(rows, rows);
|
||||
DenseMatrix refS = DenseMatrix::Zero(rows, rows);
|
||||
SparseMatrixType mUp(rows, rows);
|
||||
SparseMatrixType mLo(rows, rows);
|
||||
SparseMatrixType mS(rows, rows);
|
||||
do {
|
||||
initSparse<Scalar>(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular);
|
||||
} while (refUp.isZero());
|
||||
refLo = refUp.transpose().conjugate();
|
||||
mLo = mUp.transpose().conjugate();
|
||||
refS = refUp + refLo;
|
||||
refS.diagonal() *= 0.5;
|
||||
mS = mUp + mLo;
|
||||
for (int k=0; k<mS.outerSize(); ++k)
|
||||
for (typename SparseMatrixType::InnerIterator it(mS,k); it; ++it)
|
||||
if (it.index() == k)
|
||||
it.valueRef() *= 0.5;
|
||||
|
||||
VERIFY_IS_APPROX(refS.adjoint(), refS);
|
||||
VERIFY_IS_APPROX(mS.transpose().conjugate(), mS);
|
||||
VERIFY_IS_APPROX(mS, refS);
|
||||
VERIFY_IS_APPROX(x=mS*b, refX=refS*b);
|
||||
VERIFY_IS_APPROX(x=mUp.template marked<UpperTriangular|SelfAdjoint>()*b, refX=refS*b);
|
||||
VERIFY_IS_APPROX(x=mLo.template marked<LowerTriangular|SelfAdjoint>()*b, refX=refS*b);
|
||||
VERIFY_IS_APPROX(x=mS.template marked<SelfAdjoint>()*b, refX=refS*b);
|
||||
}
|
||||
|
||||
// test prune
|
||||
{
|
||||
SparseMatrixType m2(rows, rows);
|
||||
|
||||
130
test/sparse_product.cpp
Normal file
130
test/sparse_product.cpp
Normal file
@ -0,0 +1,130 @@
|
||||
// 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 Daniel Gomez Ferro <dgomezferro@gmail.com>
|
||||
//
|
||||
// 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/>.
|
||||
|
||||
#include "sparse.h"
|
||||
|
||||
template<typename SparseMatrixType> void sparse_product(const SparseMatrixType& ref)
|
||||
{
|
||||
const int rows = ref.rows();
|
||||
const int cols = ref.cols();
|
||||
typedef typename SparseMatrixType::Scalar Scalar;
|
||||
enum { Flags = SparseMatrixType::Flags };
|
||||
|
||||
double density = std::max(8./(rows*cols), 0.01);
|
||||
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
|
||||
typedef Matrix<Scalar,Dynamic,1> DenseVector;
|
||||
|
||||
// test matrix-matrix product
|
||||
{
|
||||
DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
|
||||
DenseMatrix refMat3 = DenseMatrix::Zero(rows, rows);
|
||||
DenseMatrix refMat4 = DenseMatrix::Zero(rows, rows);
|
||||
DenseMatrix dm4 = DenseMatrix::Zero(rows, rows);
|
||||
SparseMatrixType m2(rows, rows);
|
||||
SparseMatrixType m3(rows, rows);
|
||||
SparseMatrixType m4(rows, rows);
|
||||
initSparse<Scalar>(density, refMat2, m2);
|
||||
initSparse<Scalar>(density, refMat3, m3);
|
||||
initSparse<Scalar>(density, refMat4, m4);
|
||||
VERIFY_IS_APPROX(m4=m2*m3, refMat4=refMat2*refMat3);
|
||||
VERIFY_IS_APPROX(m4=m2.transpose()*m3, refMat4=refMat2.transpose()*refMat3);
|
||||
VERIFY_IS_APPROX(m4=m2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
|
||||
VERIFY_IS_APPROX(m4=m2*m3.transpose(), refMat4=refMat2*refMat3.transpose());
|
||||
|
||||
// sparse * dense
|
||||
VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3);
|
||||
VERIFY_IS_APPROX(dm4=m2*refMat3.transpose(), refMat4=refMat2*refMat3.transpose());
|
||||
VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3, refMat4=refMat2.transpose()*refMat3);
|
||||
VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
|
||||
|
||||
// dense * sparse
|
||||
VERIFY_IS_APPROX(dm4=refMat2*m3, refMat4=refMat2*refMat3);
|
||||
VERIFY_IS_APPROX(dm4=refMat2*m3.transpose(), refMat4=refMat2*refMat3.transpose());
|
||||
VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3, refMat4=refMat2.transpose()*refMat3);
|
||||
VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
|
||||
|
||||
VERIFY_IS_APPROX(m3=m3*m3, refMat3=refMat3*refMat3);
|
||||
}
|
||||
|
||||
// test matrix - diagonal product
|
||||
if(false) // it compiles, but the precision is terrible. probably doesn't matter in this branch....
|
||||
{
|
||||
DenseMatrix refM2 = DenseMatrix::Zero(rows, rows);
|
||||
DenseMatrix refM3 = DenseMatrix::Zero(rows, rows);
|
||||
DiagonalMatrix<DenseVector> d1(DenseVector::Random(rows));
|
||||
SparseMatrixType m2(rows, rows);
|
||||
SparseMatrixType m3(rows, rows);
|
||||
initSparse<Scalar>(density, refM2, m2);
|
||||
initSparse<Scalar>(density, refM3, m3);
|
||||
VERIFY_IS_APPROX(m3=m2*d1, refM3=refM2*d1);
|
||||
VERIFY_IS_APPROX(m3=m2.transpose()*d1, refM3=refM2.transpose()*d1);
|
||||
VERIFY_IS_APPROX(m3=d1*m2, refM3=d1*refM2);
|
||||
VERIFY_IS_APPROX(m3=d1*m2.transpose(), refM3=d1 * refM2.transpose());
|
||||
}
|
||||
|
||||
// test self adjoint products
|
||||
{
|
||||
DenseMatrix b = DenseMatrix::Random(rows, rows);
|
||||
DenseMatrix x = DenseMatrix::Random(rows, rows);
|
||||
DenseMatrix refX = DenseMatrix::Random(rows, rows);
|
||||
DenseMatrix refUp = DenseMatrix::Zero(rows, rows);
|
||||
DenseMatrix refLo = DenseMatrix::Zero(rows, rows);
|
||||
DenseMatrix refS = DenseMatrix::Zero(rows, rows);
|
||||
SparseMatrixType mUp(rows, rows);
|
||||
SparseMatrixType mLo(rows, rows);
|
||||
SparseMatrixType mS(rows, rows);
|
||||
do {
|
||||
initSparse<Scalar>(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular);
|
||||
} while (refUp.isZero());
|
||||
refLo = refUp.transpose().conjugate();
|
||||
mLo = mUp.transpose().conjugate();
|
||||
refS = refUp + refLo;
|
||||
refS.diagonal() *= 0.5;
|
||||
mS = mUp + mLo;
|
||||
for (int k=0; k<mS.outerSize(); ++k)
|
||||
for (typename SparseMatrixType::InnerIterator it(mS,k); it; ++it)
|
||||
if (it.index() == k)
|
||||
it.valueRef() *= 0.5;
|
||||
|
||||
VERIFY_IS_APPROX(refS.adjoint(), refS);
|
||||
VERIFY_IS_APPROX(mS.transpose().conjugate(), mS);
|
||||
VERIFY_IS_APPROX(mS, refS);
|
||||
VERIFY_IS_APPROX(x=mS*b, refX=refS*b);
|
||||
VERIFY_IS_APPROX(x=mUp.template marked<UpperTriangular|SelfAdjoint>()*b, refX=refS*b);
|
||||
VERIFY_IS_APPROX(x=mLo.template marked<LowerTriangular|SelfAdjoint>()*b, refX=refS*b);
|
||||
VERIFY_IS_APPROX(x=mS.template marked<SelfAdjoint>()*b, refX=refS*b);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void test_sparse_product()
|
||||
{
|
||||
for(int i = 0; i < g_repeat; i++) {
|
||||
CALL_SUBTEST( sparse_product(SparseMatrix<double>(8, 8)) );
|
||||
CALL_SUBTEST( sparse_product(SparseMatrix<std::complex<double> >(16, 16)) );
|
||||
CALL_SUBTEST( sparse_product(SparseMatrix<double>(33, 33)) );
|
||||
|
||||
CALL_SUBTEST( sparse_product(DynamicSparseMatrix<double>(8, 8)) );
|
||||
}
|
||||
}
|
||||
@ -84,6 +84,7 @@ template<typename Scalar> void sparse_vector(int rows, int cols)
|
||||
VERIFY_IS_APPROX(v1-=v2, refV1-=refV2);
|
||||
|
||||
VERIFY_IS_APPROX(v1.dot(v2), refV1.dot(refV2));
|
||||
VERIFY_IS_APPROX(v1.dot(refV2), refV1.dot(refV2));
|
||||
|
||||
}
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user