GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-09-25 07:43:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

* update sparse module wrt new diagonal matrix impl

Browse Source 
* fix a bug is SparseMatrix
This commit is contained in:
Gael Guennebaud 2009-07-04 11:16:27 +02:00
parent 2de9b7f537
commit 08e419dcb1
8 changed files with 75 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Eigen/src/Core/DiagonalMatrix.h
View File

@ -38,7 +38,8 @@ class DiagonalBase : public MultiplierBase<Derived>
ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime,
MaxRowsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime,
MaxColsAtCompileTime = DiagonalVectorType::MaxSizeAtCompileTime,
IsVectorAtCompileTime = 0
IsVectorAtCompileTime = 0,
Flags = 0
};
typedef Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime, 0, MaxRowsAtCompileTime, MaxColsAtCompileTime> DenseMatrixType;
@ -52,6 +53,9 @@ class DiagonalBase : public MultiplierBase<Derived>
inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); }
inline DiagonalVectorType& diagonal() { return derived().diagonal(); }
inline int rows() const { return diagonal().size(); }
inline int cols() const { return diagonal().size(); }
template<typename MatrixDerived>
const DiagonalProduct<MatrixDerived, Derived, DiagonalOnTheLeft>

12
Eigen/src/Core/util/Constants.h
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
@ -185,9 +185,6 @@ const unsigned int HereditaryBits = RowMajorBit
| EvalBeforeNestingBit
| EvalBeforeAssigningBit
| SparseBit;
// diagonal means both upper and lower triangular
const unsigned DiagonalBits = UpperTriangularBit | LowerTriangularBit;
// Possible values for the Mode parameter of part()
const unsigned int UpperTriangular = UpperTriangularBit;
@ -198,13 +195,6 @@ const unsigned int SelfAdjoint = SelfAdjointBit;
const unsigned int UnitUpperTriangular = UpperTriangularBit | UnitDiagBit;
const unsigned int UnitLowerTriangular = LowerTriangularBit | UnitDiagBit;
template<typename T> struct ei_is_diagonal
{
enum {
ret = ( (unsigned int)(T::Flags) & DiagonalBits ) == DiagonalBits
};
};
enum { DiagonalOnTheLeft, DiagonalOnTheRight };
enum { Aligned, Unaligned };

12
Eigen/src/Core/util/Meta.h
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
@ -186,6 +186,16 @@ struct ei_result_of<ei_scalar_product_op<Scalar>(ArgType0,ArgType1)> {
typedef typename ei_scalar_product_traits<typename ei_cleantype<ArgType0>::type, typename ei_cleantype<ArgType1>::type>::ReturnType type;
};
template<typename T> struct ei_is_diagonal
{ enum { ret = false }; };
template<typename T> struct ei_is_diagonal<DiagonalBase<T> >
{ enum { ret = true }; };
template<typename T> struct ei_is_diagonal<DiagonalWrapper<T> >
{ enum { ret = true }; };
template<typename T, int S> struct ei_is_diagonal<DiagonalMatrix<T,S> >
{ enum { ret = true }; };
#endif // EIGEN_META_H

52
Eigen/src/Sparse/SparseDiagonalProduct.h
View File

@ -25,8 +25,9 @@
#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:
// The product of a diagonal matrix with a sparse matrix can be easily
// implemented using expression template.
// We have two consider very different cases:
// 1 - diag * row-major sparse
// => each inner vector <=> scalar * sparse vector product
// => so we can reuse CwiseUnaryOp::InnerIterator
@ -37,13 +38,21 @@
// The two other cases are symmetric.
template<typename Lhs, typename Rhs>
struct ei_traits<SparseDiagonalProduct<Lhs, Rhs> > : ei_traits<SparseProduct<Lhs, Rhs, DiagonalProduct> >
struct ei_traits<SparseDiagonalProduct<Lhs, Rhs> >
{
typedef typename ei_cleantype<Lhs>::type _Lhs;
typedef typename ei_cleantype<Rhs>::type _Rhs;
typedef typename _Lhs::Scalar Scalar;
enum {
RowsAtCompileTime = _Lhs::RowsAtCompileTime,
ColsAtCompileTime = _Rhs::ColsAtCompileTime,
MaxRowsAtCompileTime = _Lhs::MaxRowsAtCompileTime,
MaxColsAtCompileTime = _Rhs::MaxColsAtCompileTime,
SparseFlags = ei_is_diagonal<_Lhs>::ret ? int(_Rhs::Flags) : int(_Lhs::Flags),
Flags = SparseBit | (SparseFlags&RowMajorBit)
Flags = SparseBit | (SparseFlags&RowMajorBit),
CoeffReadCost = Dynamic
};
};
@ -51,11 +60,16 @@ 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
template<typename Lhs, typename Rhs>
class SparseDiagonalProduct
: public SparseMatrixBase<SparseDiagonalProduct<Lhs,Rhs> >,
ei_no_assignment_operator
{
typedef typename ei_traits<SparseDiagonalProduct>::_LhsNested _LhsNested;
typedef typename ei_traits<SparseDiagonalProduct>::_RhsNested _RhsNested;
typedef typename Lhs::Nested LhsNested;
typedef typename Rhs::Nested RhsNested;
typedef typename ei_cleantype<LhsNested>::type _LhsNested;
typedef typename ei_cleantype<RhsNested>::type _RhsNested;
enum {
LhsMode = ei_is_diagonal<_LhsNested>::ret ? SDP_IsDiagonal
@ -71,8 +85,8 @@ class SparseDiagonalProduct : public SparseMatrixBase<SparseDiagonalProduct<LhsN
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)
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");
@ -109,12 +123,12 @@ class ei_sparse_diagonal_product_inner_iterator_selector
: public SparseCwiseBinaryOp<
ei_scalar_product_op<typename Lhs::Scalar>,
SparseInnerVectorSet<Rhs,1>,
typename Lhs::_CoeffsVectorType>::InnerIterator
typename Lhs::DiagonalVectorType>::InnerIterator
{
typedef typename SparseCwiseBinaryOp<
ei_scalar_product_op<typename Lhs::Scalar>,
SparseInnerVectorSet<Rhs,1>,
typename Lhs::_CoeffsVectorType>::InnerIterator Base;
typename Lhs::DiagonalVectorType>::InnerIterator Base;
public:
inline ei_sparse_diagonal_product_inner_iterator_selector(
const SparseDiagonalProductType& expr, int outer)
@ -141,12 +155,12 @@ class ei_sparse_diagonal_product_inner_iterator_selector
: public SparseCwiseBinaryOp<
ei_scalar_product_op<typename Rhs::Scalar>,
SparseInnerVectorSet<Lhs,1>,
NestByValue<Transpose<typename Rhs::_CoeffsVectorType> > >::InnerIterator
NestByValue<Transpose<typename Rhs::DiagonalVectorType> > >::InnerIterator
{
typedef typename SparseCwiseBinaryOp<
ei_scalar_product_op<typename Rhs::Scalar>,
SparseInnerVectorSet<Lhs,1>,
NestByValue<Transpose<typename Rhs::_CoeffsVectorType> > >::InnerIterator Base;
NestByValue<Transpose<typename Rhs::DiagonalVectorType> > >::InnerIterator Base;
public:
inline ei_sparse_diagonal_product_inner_iterator_selector(
const SparseDiagonalProductType& expr, int outer)
@ -154,4 +168,14 @@ class ei_sparse_diagonal_product_inner_iterator_selector
{}
};
// SparseMatrixBase functions
template<typename Derived>
template<typename OtherDerived>
const SparseDiagonalProduct<Derived,OtherDerived>
SparseMatrixBase<Derived>::operator*(const DiagonalBase<OtherDerived> &other) const
{
return SparseDiagonalProduct<Derived,OtherDerived>(this->derived(), other.derived());
}
#endif // EIGEN_SPARSE_DIAGONAL_PRODUCT_H

2
Eigen/src/Sparse/SparseMatrix.h
View File

@ -371,7 +371,7 @@ class SparseMatrix
const int outerSize = IsRowMajor ? rows : cols;
m_innerSize = IsRowMajor ? cols : rows;
m_data.clear();
if (m_outerSize != outerSize)
if (m_outerSize != outerSize || m_outerSize==0)
{
delete[] m_outerIndex;
m_outerIndex = new int [outerSize+1];

17
Eigen/src/Sparse/SparseMatrixBase.h
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// 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
@ -37,6 +37,9 @@
*
*/
template<typename Derived> class SparseMatrixBase
#ifndef EIGEN_PARSED_BY_DOXYGEN
: public MultiplierBase<Derived>
#endif // not EIGEN_PARSED_BY_DOXYGEN
{
public:
@ -301,10 +304,22 @@ template<typename Derived> class SparseMatrixBase
{ return matrix*scalar; }
// sparse * sparse
template<typename OtherDerived>
const typename SparseProductReturnType<Derived,OtherDerived>::Type
operator*(const SparseMatrixBase<OtherDerived> &other) const;
// sparse * diagonal
template<typename OtherDerived>
const SparseDiagonalProduct<Derived,OtherDerived>
operator*(const DiagonalBase<OtherDerived> &other) const;
// diagonal * sparse
template<typename OtherDerived> friend
const SparseDiagonalProduct<OtherDerived,Derived>
operator*(const DiagonalBase<OtherDerived> &lhs, const SparseMatrixBase& rhs)
{ return SparseDiagonalProduct<OtherDerived,Derived>(lhs.derived(), rhs.derived()); }
// dense * sparse (return a dense object)
template<typename OtherDerived> friend
const typename SparseProductReturnType<OtherDerived,Derived>::Type

17
Eigen/src/Sparse/SparseProduct.h
View File

@ -1,7 +1,7 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// 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
@ -29,9 +29,7 @@ template<typename Lhs, typename Rhs> struct ei_sparse_product_mode
{
enum {
value = ei_is_diagonal<Lhs>::ret || ei_is_diagonal<Rhs>::ret
? DiagonalProduct
: (Rhs::Flags&Lhs::Flags&SparseBit)==SparseBit
value = (Rhs::Flags&Lhs::Flags&SparseBit)==SparseBit
? SparseTimeSparseProduct
: (Lhs::Flags&SparseBit)==SparseBit
? SparseTimeDenseProduct
@ -47,15 +45,6 @@ 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>
@ -106,7 +95,7 @@ struct ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >
// RhsIsRowMajor = (RhsFlags & RowMajorBit)==RowMajorBit,
EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit),
ResultIsSparse = ProductMode==SparseTimeSparseProduct || ProductMode==DiagonalProduct,
ResultIsSparse = ProductMode==SparseTimeSparseProduct,
RemovedBits = ~( (EvalToRowMajor ? 0 : RowMajorBit) | (ResultIsSparse ? 0 : SparseBit) ),

1
test/sparse_product.cpp
View File

@ -118,7 +118,6 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType&
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()