mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-30 17:40:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

Big refactoring/cleaning in the spasre module with

Browse Source 
in particular the addition of a selfadjointView, and the
extension of triangularView. The rest is cleaning and does not
change/extend the API.
This commit is contained in:
Gael Guennebaud 2009-11-18 14:52:52 +01:00
parent 1e62e0b0d8
commit 0529ecfe1b
18 changed files with 451 additions and 274 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
Eigen/Sparse
View File

@ -98,7 +98,6 @@ struct Sparse {};
#include "src/Sparse/SparseVector.h"
#include "src/Sparse/CoreIterators.h"
#include "src/Sparse/SparseTranspose.h"
#include "src/Sparse/SparseCwise.h"
#include "src/Sparse/SparseCwiseUnaryOp.h"
#include "src/Sparse/SparseCwiseBinaryOp.h"
#include "src/Sparse/SparseDot.h"
@ -107,12 +106,12 @@ struct Sparse {};
#include "src/Sparse/SparseFuzzy.h"
#include "src/Sparse/SparseProduct.h"
#include "src/Sparse/SparseDiagonalProduct.h"
#include "src/Sparse/SparseTriangular.h"
#include "src/Sparse/SparseTriangularView.h"
#include "src/Sparse/SparseSelfAdjointView.h"
#include "src/Sparse/TriangularSolver.h"
#include "src/Sparse/SparseLLT.h"
#include "src/Sparse/SparseLDLT.h"
#include "src/Sparse/SparseLU.h"
#include "src/Sparse/SparseExpressionMaker.h"
#ifdef EIGEN_CHOLMOD_SUPPORT
# include "src/Sparse/CholmodSupport.h"

9
Eigen/src/Core/MatrixBase.h
View File

@ -710,15 +710,6 @@ template<typename Derived> class MatrixBase
typedef Homogeneous<Derived,MatrixBase<Derived>::ColsAtCompileTime==1?Vertical:Horizontal> HomogeneousReturnType;
const HomogeneousReturnType homogeneous() const;
/////////// Sparse module ///////////
// dense = sparse * dense
template<typename Derived1, typename Derived2>
Derived& lazyAssign(const SparseProduct<Derived1,Derived2,SparseTimeDenseProduct>& product);
// dense = dense * sparse
template<typename Derived1, typename Derived2>
Derived& lazyAssign(const SparseProduct<Derived1,Derived2,DenseTimeSparseProduct>& product);
////////// Householder module ///////////
void makeHouseholderInPlace(Scalar& tau, RealScalar& beta);

8
Eigen/src/Core/SelfAdjointView.h
View File

@ -201,15 +201,15 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, SelfAdjoint, Dynami
***************************************************************************/
template<typename Derived>
template<unsigned int Mode>
const SelfAdjointView<Derived, Mode> MatrixBase<Derived>::selfadjointView() const
template<unsigned int UpLo>
const SelfAdjointView<Derived, UpLo> MatrixBase<Derived>::selfadjointView() const
{
return derived();
}
template<typename Derived>
template<unsigned int Mode>
SelfAdjointView<Derived, Mode> MatrixBase<Derived>::selfadjointView()
template<unsigned int UpLo>
SelfAdjointView<Derived, UpLo> MatrixBase<Derived>::selfadjointView()
{
return derived();
}

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

@ -193,7 +193,7 @@ enum { AsRequested=0, EnforceAlignedAccess=2 };
enum { ConditionalJumpCost = 5 };
enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight };
enum DirectionType { Vertical, Horizontal, BothDirections };
enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, SparseTimeSparseProduct, SparseTimeDenseProduct, DenseTimeSparseProduct };
enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct };
enum {
/** \internal Equivalent to a slice vectorization for fixed-size matrices having good alignment

4
Eigen/src/Core/util/ForwardDeclarations.h
View File

@ -146,8 +146,4 @@ template<typename Scalar,int Dim> class Translation;
template<typename Scalar> class UniformScaling;
template<typename MatrixType,int Direction> class Homogeneous;
// Sparse module:
template<typename Lhs, typename Rhs, int ProductMode> class SparseProduct;
#endif // EIGEN_FORWARDDECLARATIONS_H

36
Eigen/src/Sparse/SparseCwiseBinaryOp.h
View File

@ -42,35 +42,6 @@
// 4 - dense op dense product dense
// generic dense
// template<typename BinaryOp, typename Lhs, typename Rhs>
// struct ei_traits<SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> >
// {
// typedef typename ei_result_of<
// BinaryOp(
// typename Lhs::Scalar,
// typename Rhs::Scalar
// )
// >::type Scalar;
// typedef typename ei_promote_storage_type<typename ei_traits<Lhs>::StorageType,
// typename ei_traits<Rhs>::StorageType>::ret StorageType;
// typedef typename Lhs::Nested LhsNested;
// typedef typename Rhs::Nested RhsNested;
// typedef typename ei_unref<LhsNested>::type _LhsNested;
// typedef typename ei_unref<RhsNested>::type _RhsNested;
// enum {
// LhsCoeffReadCost = _LhsNested::CoeffReadCost,
// RhsCoeffReadCost = _RhsNested::CoeffReadCost,
// LhsFlags = _LhsNested::Flags,
// RhsFlags = _RhsNested::Flags,
// RowsAtCompileTime = Lhs::RowsAtCompileTime,
// ColsAtCompileTime = Lhs::ColsAtCompileTime,
// MaxRowsAtCompileTime = Lhs::MaxRowsAtCompileTime,
// MaxColsAtCompileTime = Lhs::MaxColsAtCompileTime,
// Flags = (int(LhsFlags) | int(RhsFlags)) & HereditaryBits,
// CoeffReadCost = LhsCoeffReadCost + RhsCoeffReadCost + ei_functor_traits<BinaryOp>::Cost
// };
// };
template<> struct ei_promote_storage_type<Dense,Sparse>
{ typedef Sparse ret; };
@ -82,16 +53,9 @@ class CwiseBinaryOpImpl<BinaryOp, Lhs, Rhs, Sparse>
: public SparseMatrixBase<CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
{
public:
class InnerIterator;
typedef CwiseBinaryOp<BinaryOp, Lhs, Rhs> Derived;
EIGEN_SPARSE_PUBLIC_INTERFACE(Derived)
// typedef typename ei_traits<SparseCwiseBinaryOp>::LhsNested LhsNested;
// typedef typename ei_traits<SparseCwiseBinaryOp>::RhsNested RhsNested;
// typedef typename ei_unref<LhsNested>::type _LhsNested;
// typedef typename ei_unref<RhsNested>::type _RhsNested;
};
template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived,

4
Eigen/src/Sparse/SparseLDLT.h
View File

@ -333,12 +333,12 @@ bool SparseLDLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
return false;
if (m_matrix.nonZeros()>0) // otherwise L==I
m_matrix.template triangular<UnitLowerTriangular>().solveInPlace(b);
m_matrix.template triangularView<UnitLowerTriangular>().solveInPlace(b);
b = b.cwise() / m_diag;
// FIXME should be .adjoint() but it fails to compile...
if (m_matrix.nonZeros()>0) // otherwise L==I
m_matrix.transpose().template triangular<UnitUpperTriangular>().solveInPlace(b);
m_matrix.transpose().template triangularView<UnitUpperTriangular>().solveInPlace(b);
return true;
}

6
Eigen/src/Sparse/SparseLLT.h
View File

@ -193,15 +193,15 @@ bool SparseLLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
const int size = m_matrix.rows();
ei_assert(size==b.rows());
m_matrix.template triangular<LowerTriangular>().solveInPlace(b);
m_matrix.template triangularView<LowerTriangular>().solveInPlace(b);
// FIXME should be simply .adjoint() but it fails to compile...
if (NumTraits<Scalar>::IsComplex)
{
CholMatrixType aux = m_matrix.conjugate();
aux.transpose().template triangular<UpperTriangular>().solveInPlace(b);
aux.transpose().template triangularView<UpperTriangular>().solveInPlace(b);
}
else
m_matrix.transpose().template triangular<UpperTriangular>().solveInPlace(b);
m_matrix.transpose().template triangularView<UpperTriangular>().solveInPlace(b);
return true;
}

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

@ -540,6 +540,7 @@ class SparseMatrix<Scalar,_Options>::InnerIterator
inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix.m_data.value(m_id)); }
inline int index() const { return m_matrix.m_data.index(m_id); }
inline int outer() const { return m_outer; }
inline int row() const { return IsRowMajor ? m_outer : index(); }
inline int col() const { return IsRowMajor ? index() : m_outer; }

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

@ -244,7 +244,7 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived
}
template<typename Lhs, typename Rhs>
inline Derived& operator=(const SparseProduct<Lhs,Rhs,SparseTimeSparseProduct>& product);
inline Derived& operator=(const SparseProduct<Lhs,Rhs>& product);
friend std::ostream & operator << (std::ostream & s, const SparseMatrixBase& m)
{
@ -337,12 +337,13 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived
// dense * sparse (return a dense object)
template<typename OtherDerived> friend
const typename SparseProductReturnType<OtherDerived,Derived>::Type
const DenseTimeSparseProduct<OtherDerived,Derived>
operator*(const MatrixBase<OtherDerived>& lhs, const Derived& rhs)
{ return typename SparseProductReturnType<OtherDerived,Derived>::Type(lhs.derived(),rhs); }
{ return DenseTimeSparseProduct<OtherDerived,Derived>(lhs.derived(),rhs); }
// sparse * dense (returns a dense object)
template<typename OtherDerived>
const typename SparseProductReturnType<Derived,OtherDerived>::Type
const SparseTimeDenseProduct<Derived,OtherDerived>
operator*(const MatrixBase<OtherDerived> &other) const;
template<typename OtherDerived>
@ -360,7 +361,10 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived
// void solveTriangularInPlace(SparseMatrixBase<OtherDerived>& other) const;
template<int Mode>
inline const SparseTriangular<Derived, Mode> triangular() const;
inline const SparseTriangularView<Derived, Mode> triangularView() const;
template<unsigned int UpLo> inline const SparseSelfAdjointView<Derived, UpLo> selfadjointView() const;
template<unsigned int UpLo> inline SparseSelfAdjointView<Derived, UpLo> selfadjointView();
template<typename OtherDerived> Scalar dot(const MatrixBase<OtherDerived>& other) const;
template<typename OtherDerived> Scalar dot(const SparseMatrixBase<OtherDerived>& other) const;

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

@ -25,22 +25,8 @@
#ifndef EIGEN_SPARSEPRODUCT_H
#define EIGEN_SPARSEPRODUCT_H
template<typename Lhs, typename Rhs> struct ei_sparse_product_mode<Lhs,Rhs,Sparse,Sparse> { enum { value = SparseTimeSparseProduct }; };
template<typename Lhs, typename Rhs> struct ei_sparse_product_mode<Lhs,Rhs,Dense, Sparse> { enum { value = DenseTimeSparseProduct }; };
template<typename Lhs, typename Rhs> struct ei_sparse_product_mode<Lhs,Rhs,Sparse,Dense> { enum { value = SparseTimeDenseProduct }; };
template<typename Lhs, typename Rhs, int ProductMode>
struct SparseProductReturnType
{
typedef const typename ei_nested<Lhs,Rhs::RowsAtCompileTime>::type LhsNested;
typedef const typename ei_nested<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
typedef SparseProduct<LhsNested, RhsNested, ProductMode> Type;
};
// sparse product return type specialization
template<typename Lhs, typename Rhs>
struct SparseProductReturnType<Lhs,Rhs,SparseTimeSparseProduct>
struct SparseProductReturnType
{
typedef typename ei_traits<Lhs>::Scalar Scalar;
enum {
@ -60,11 +46,11 @@ struct SparseProductReturnType<Lhs,Rhs,SparseTimeSparseProduct>
SparseMatrix<Scalar,0>,
const typename ei_nested<Rhs,Lhs::RowsAtCompileTime>::type>::ret RhsNested;
typedef SparseProduct<LhsNested, RhsNested, SparseTimeSparseProduct> Type;
typedef SparseProduct<LhsNested, RhsNested> Type;
};
template<typename LhsNested, typename RhsNested, int ProductMode>
struct ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >
template<typename LhsNested, typename RhsNested>
struct ei_traits<SparseProduct<LhsNested, RhsNested> >
{
// clean the nested types:
typedef typename ei_cleantype<LhsNested>::type _LhsNested;
@ -85,7 +71,6 @@ struct ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >
MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,
EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit),
ResultIsSparse = ProductMode==SparseTimeSparseProduct,
RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit),
@ -96,16 +81,14 @@ struct ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >
CoeffReadCost = Dynamic
};
typedef typename ei_meta_if<ResultIsSparse, Sparse, Dense>::ret StorageType;
typedef Sparse StorageType;
typedef typename ei_meta_if<ResultIsSparse,
SparseMatrixBase<SparseProduct<LhsNested, RhsNested, ProductMode> >,
MatrixBase<SparseProduct<LhsNested, RhsNested, ProductMode> > >::ret Base;
typedef SparseMatrixBase<SparseProduct<LhsNested, RhsNested> > Base;
};
template<typename LhsNested, typename RhsNested, int ProductMode>
template<typename LhsNested, typename RhsNested>
class SparseProduct : ei_no_assignment_operator,
public ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >::Base
public ei_traits<SparseProduct<LhsNested, RhsNested> >::Base
{
public:
@ -256,38 +239,14 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
}
};
// NOTE eventually let's transpose one argument even in this case since it might be expensive if
// the result is not dense.
// template<typename Lhs, typename Rhs, typename ResultType, int ResStorageOrder>
// struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,ColMajor,ResStorageOrder>
// {
// static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
// {
// // trivial product as lhs.row/rhs.col dot products
// // loop over the preferred order of the result
// }
// };
// NOTE the 2 others cases (col row *) must never occurs since they are caught
// by ProductReturnType which transform it to (col col *) by evaluating rhs.
// template<typename Derived>
// template<typename Lhs, typename Rhs>
// inline Derived& SparseMatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs>& product)
// {
// // std::cout << "sparse product to dense\n";
// ei_sparse_product_selector<
// typename ei_cleantype<Lhs>::type,
// typename ei_cleantype<Rhs>::type,
// typename ei_cleantype<Derived>::type>::run(product.lhs(),product.rhs(),derived());
// return derived();
// }
// sparse = sparse * sparse
template<typename Derived>
template<typename Lhs, typename Rhs>
inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs,SparseTimeSparseProduct>& product)
inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs>& product)
{
ei_sparse_product_selector<
typename ei_cleantype<Lhs>::type,
@ -296,78 +255,79 @@ inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs
return derived();
}
// dense = sparse * dense
// Note that here we force no inlining and separate the setZero() because GCC messes up otherwise
template<typename Lhs, typename Rhs, typename Dest>
EIGEN_DONT_INLINE void ei_sparse_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst)
{
typedef typename ei_cleantype<Lhs>::type _Lhs;
typedef typename ei_cleantype<Rhs>::type _Rhs;
typedef typename _Lhs::InnerIterator LhsInnerIterator;
enum {
LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit,
LhsIsSelfAdjoint = (_Lhs::Flags&SelfAdjointBit)==SelfAdjointBit,
ProcessFirstHalf = LhsIsSelfAdjoint
&& ( ((_Lhs::Flags&(UpperTriangularBit|LowerTriangularBit))==0)
|| ( (_Lhs::Flags&UpperTriangularBit) && !LhsIsRowMajor)
|| ( (_Lhs::Flags&LowerTriangularBit) && LhsIsRowMajor) ),
ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf)
};
for (int j=0; j<lhs.outerSize(); ++j)
{
LhsInnerIterator i(lhs,j);
if (ProcessSecondHalf && i && (i.index()==j))
{
dst.row(j) += i.value() * rhs.row(j);
++i;
}
typename Rhs::Scalar rhs_j = rhs.coeff(j,0);
Block<Dest,1,Dest::ColsAtCompileTime> dst_j(dst.row(LhsIsRowMajor ? j : 0));
for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
{
if(LhsIsSelfAdjoint)
{
int a = LhsIsRowMajor ? j : i.index();
int b = LhsIsRowMajor ? i.index() : j;
typename Lhs::Scalar v = i.value();
dst.row(a) += (v) * rhs.row(b);
dst.row(b) += ei_conj(v) * rhs.row(a);
}
else if(LhsIsRowMajor)
dst_j += i.value() * rhs.row(i.index());
else if(Rhs::ColsAtCompileTime==1)
dst.coeffRef(i.index()) += i.value() * rhs_j;
else
dst.row(i.index()) += i.value() * rhs.row(j);
}
if (ProcessFirstHalf && i && (i.index()==j))
dst.row(j) += i.value() * rhs.row(j);
}
}
template<typename Derived>
template<typename Lhs, typename Rhs>
Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,SparseTimeDenseProduct>& product)
struct ei_traits<SparseTimeDenseProduct<Lhs,Rhs> >
: ei_traits<ProductBase<SparseTimeDenseProduct<Lhs,Rhs>, Lhs, Rhs> >
{
derived().setZero();
ei_sparse_time_dense_product(product.lhs(), product.rhs(), derived());
return derived();
}
typedef Dense StorageType;
};
template<typename Lhs, typename Rhs>
class SparseTimeDenseProduct
: public ProductBase<SparseTimeDenseProduct<Lhs,Rhs>, Lhs, Rhs>
{
public:
EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseTimeDenseProduct)
SparseTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
{}
template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const
{
typedef typename ei_cleantype<Lhs>::type _Lhs;
typedef typename ei_cleantype<Rhs>::type _Rhs;
typedef typename _Lhs::InnerIterator LhsInnerIterator;
enum { LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit };
for(int j=0; j<m_lhs.outerSize(); ++j)
{
typename Rhs::Scalar rhs_j = alpha * m_rhs.coeff(j,0);
Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest.row(LhsIsRowMajor ? j : 0));
for(LhsInnerIterator it(m_lhs,j); it ;++it)
{
if(LhsIsRowMajor) dest_j += (alpha*it.value()) * m_rhs.row(it.index());
else if(Rhs::ColsAtCompileTime==1) dest.coeffRef(it.index()) += it.value() * rhs_j;
else dest.row(it.index()) += (alpha*it.value()) * m_rhs.row(j);
}
}
}
private:
SparseTimeDenseProduct& operator=(const SparseTimeDenseProduct&);
};
// dense = dense * sparse
template<typename Derived>
template<typename Lhs, typename Rhs>
Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,DenseTimeSparseProduct>& product)
struct ei_traits<DenseTimeSparseProduct<Lhs,Rhs> >
: ei_traits<ProductBase<DenseTimeSparseProduct<Lhs,Rhs>, Lhs, Rhs> >
{
typedef typename ei_cleantype<Rhs>::type _Rhs;
typedef typename _Rhs::InnerIterator RhsInnerIterator;
enum { RhsIsRowMajor = (_Rhs::Flags&RowMajorBit)==RowMajorBit };
derived().setZero();
for (int j=0; j<product.rhs().outerSize(); ++j)
for (RhsInnerIterator i(product.rhs(),j); i; ++i)
derived().col(RhsIsRowMajor ? i.index() : j) += i.value() * product.lhs().col(RhsIsRowMajor ? j : i.index());
return derived();
}
typedef Dense StorageType;
};
template<typename Lhs, typename Rhs>
class DenseTimeSparseProduct
: public ProductBase<DenseTimeSparseProduct<Lhs,Rhs>, Lhs, Rhs>
{
public:
EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseProduct)
DenseTimeSparseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
{}
template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const
{
typedef typename ei_cleantype<Rhs>::type _Rhs;
typedef typename _Rhs::InnerIterator RhsInnerIterator;
enum { RhsIsRowMajor = (_Rhs::Flags&RowMajorBit)==RowMajorBit };
for(int j=0; j<m_rhs.outerSize(); ++j)
for(RhsInnerIterator i(m_rhs,j); i; ++i)
dest.col(RhsIsRowMajor ? i.index() : j) += (alpha*i.value()) * m_lhs.col(RhsIsRowMajor ? j : i.index());
}
private:
DenseTimeSparseProduct& operator=(const DenseTimeSparseProduct&);
};
// sparse * sparse
template<typename Derived>
@ -381,10 +341,10 @@ SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other
// sparse * dense
template<typename Derived>
template<typename OtherDerived>
EIGEN_STRONG_INLINE const typename SparseProductReturnType<Derived,OtherDerived>::Type
EIGEN_STRONG_INLINE const SparseTimeDenseProduct<Derived,OtherDerived>
SparseMatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const
{
return typename SparseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
return SparseTimeDenseProduct<Derived,OtherDerived>(derived(), other.derived());
}
#endif // EIGEN_SPARSEPRODUCT_H

223
Eigen/src/Sparse/SparseSelfAdjointView.h Normal file
View File

@ -0,0 +1,223 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// 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_SELFADJOINTVIEW_H
#define EIGEN_SPARSE_SELFADJOINTVIEW_H
/** \class SparseSelfAdjointView
* \nonstableyet
*
* \brief Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix.
*
* \param MatrixType the type of the dense matrix storing the coefficients
* \param UpLo can be either \c LowerTriangular or \c UpperTriangular
*
* This class is an expression of a sefladjoint matrix from a triangular part of a matrix
* with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView()
* and most of the time this is the only way that it is used.
*
* \sa SparseMatrixBase::selfAdjointView()
*/
template<typename Lhs, typename Rhs, int UpLo>
class SparseSelfAdjointTimeDenseProduct;
template<typename Lhs, typename Rhs, int UpLo>
class DenseTimeSparseSelfAdjointProduct;
template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView
{
public:
typedef typename ei_traits<MatrixType>::Scalar Scalar;
inline SparseSelfAdjointView(const MatrixType& matrix) : m_matrix(matrix)
{
ei_assert(ei_are_flags_consistent<UpLo>::ret);
ei_assert(rows()==cols() && "SelfAdjointView is only for squared matrices");
}
inline int rows() const { return m_matrix.rows(); }
inline int cols() const { return m_matrix.cols(); }
/** \internal \returns a reference to the nested matrix */
const MatrixType& matrix() const { return m_matrix; }
/** Efficient sparse self-adjoint matrix times dense vector/matrix product */
template<typename OtherDerived>
SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>
operator*(const MatrixBase<OtherDerived>& rhs) const
{
return SparseSelfAdjointTimeDenseProduct<MatrixType,OtherDerived,UpLo>(m_matrix, rhs.derived());
}
/** Efficient dense vector/matrix times sparse self-adjoint matrix product */
template<typename OtherDerived> friend
DenseTimeSparseSelfAdjointProduct<OtherDerived,MatrixType,UpLo>
operator*(const MatrixBase<OtherDerived>& lhs, const SparseSelfAdjointView& rhs)
{
return DenseTimeSparseSelfAdjointProduct<OtherDerived,MatrixType,UpLo>(lhs.derived(), rhs.m_matrix);
}
/** Perform a symmetric rank K update of the selfadjoint matrix \c *this:
* \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix.
*
* \returns a reference to \c *this
*
* Note that it is faster to set alpha=0 than initializing the matrix to zero
* and then keep the default value alpha=1.
*
* To perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply
* call this function with u.adjoint().
*/
template<typename DerivedU>
SparseSelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
// const SparseLLT<PlainMatrixType, UpLo> llt() const;
// const SparseLDLT<PlainMatrixType, UpLo> ldlt() const;
protected:
const typename MatrixType::Nested m_matrix;
};
/***************************************************************************
* Implementation of SparseMatrixBase methods
***************************************************************************/
template<typename Derived>
template<unsigned int UpLo>
const SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView() const
{
return derived();
}
template<typename Derived>
template<unsigned int UpLo>
SparseSelfAdjointView<Derived, UpLo> SparseMatrixBase<Derived>::selfadjointView()
{
return derived();
}
/***************************************************************************
* Implementation of SparseSelfAdjointView methods
***************************************************************************/
template<typename MatrixType, unsigned int UpLo>
template<typename DerivedU>
SparseSelfAdjointView<MatrixType,UpLo>&
SparseSelfAdjointView<MatrixType,UpLo>::rankUpdate(const MatrixBase<DerivedU>& u, Scalar alpha)
{
SparseMatrix<Scalar,MatrixType::Flags&RowMajorBit?RowMajor:ColMajor> tmp = u * u.adjoint();
if(alpha==Scalar(0))
m_matrix = tmp.template triangularView<UpLo>();
else
m_matrix += alpha * tmp.template triangularView<UpLo>();
return this;
}
/***************************************************************************
* Implementation of sparse self-adjoint time dense matrix
***************************************************************************/
template<typename Lhs, typename Rhs, int UpLo>
struct ei_traits<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo> >
: ei_traits<ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
{};
template<typename Lhs, typename Rhs, int UpLo>
class SparseSelfAdjointTimeDenseProduct
: public ProductBase<SparseSelfAdjointTimeDenseProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
{
public:
EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseSelfAdjointTimeDenseProduct)
SparseSelfAdjointTimeDenseProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
{}
template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const
{
// TODO use alpha
ei_assert(alpha==Scalar(1) && "alpha != 1 is not implemented yet, sorry");
typedef typename ei_cleantype<Lhs>::type _Lhs;
typedef typename ei_cleantype<Rhs>::type _Rhs;
typedef typename _Lhs::InnerIterator LhsInnerIterator;
enum {
LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit,
ProcessFirstHalf =
((UpLo&(UpperTriangularBit|LowerTriangularBit))==(UpperTriangularBit|LowerTriangularBit))
|| ( (UpLo&UpperTriangularBit) && !LhsIsRowMajor)
|| ( (UpLo&LowerTriangularBit) && LhsIsRowMajor),
ProcessSecondHalf = !ProcessFirstHalf
};
for (int j=0; j<m_lhs.outerSize(); ++j)
{
LhsInnerIterator i(m_lhs,j);
if (ProcessSecondHalf && i && (i.index()==j))
{
dest.row(j) += i.value() * m_rhs.row(j);
++i;
}
Block<Dest,1,Dest::ColsAtCompileTime> dest_j(dest.row(LhsIsRowMajor ? j : 0));
for(; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
{
int a = LhsIsRowMajor ? j : i.index();
int b = LhsIsRowMajor ? i.index() : j;
typename Lhs::Scalar v = i.value();
dest.row(a) += (v) * m_rhs.row(b);
dest.row(b) += ei_conj(v) * m_rhs.row(a);
}
if (ProcessFirstHalf && i && (i.index()==j))
dest.row(j) += i.value() * m_rhs.row(j);
}
}
private:
SparseSelfAdjointTimeDenseProduct& operator=(const SparseSelfAdjointTimeDenseProduct&);
};
template<typename Lhs, typename Rhs, int UpLo>
struct ei_traits<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo> >
: ei_traits<ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs> >
{};
template<typename Lhs, typename Rhs, int UpLo>
class DenseTimeSparseSelfAdjointProduct
: public ProductBase<DenseTimeSparseSelfAdjointProduct<Lhs,Rhs,UpLo>, Lhs, Rhs>
{
public:
EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseSelfAdjointProduct)
DenseTimeSparseSelfAdjointProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
{}
template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const
{
// TODO
}
private:
DenseTimeSparseSelfAdjointProduct& operator=(const DenseTimeSparseSelfAdjointProduct&);
};
#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H

60
Eigen/src/Sparse/SparseTriangular.h
View File

@ -1,60 +0,0 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// 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_TRIANGULAR_H
#define EIGEN_SPARSE_TRIANGULAR_H
template<typename ExpressionType, int Mode> class SparseTriangular
{
public:
typedef typename ei_traits<ExpressionType>::Scalar Scalar;
typedef typename ei_meta_if<ei_must_nest_by_value<ExpressionType>::ret,
ExpressionType, const ExpressionType&>::ret ExpressionTypeNested;
typedef CwiseUnaryOp<ei_scalar_add_op<Scalar>, ExpressionType> ScalarAddReturnType;
inline SparseTriangular(const ExpressionType& matrix) : m_matrix(matrix) {}
/** \internal */
inline const ExpressionType& _expression() const { return m_matrix; }
template<typename OtherDerived>
typename ei_plain_matrix_type_column_major<OtherDerived>::type
solve(const MatrixBase<OtherDerived>& other) const;
template<typename OtherDerived> void solveInPlace(MatrixBase<OtherDerived>& other) const;
template<typename OtherDerived> void solveInPlace(SparseMatrixBase<OtherDerived>& other) const;
protected:
ExpressionTypeNested m_matrix;
};
template<typename Derived>
template<int Mode>
inline const SparseTriangular<Derived, Mode>
SparseMatrixBase<Derived>::triangular() const
{
return derived();
}
#endif // EIGEN_SPARSE_TRIANGULAR_H

95
Eigen/src/Sparse/SparseTriangularView.h Normal file
View File

@ -0,0 +1,95 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// 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_TRIANGULARVIEW_H
#define EIGEN_SPARSE_TRIANGULARVIEW_H
template<typename MatrixType, int Mode>
struct ei_traits<SparseTriangularView<MatrixType,Mode> >
: public ei_traits<MatrixType>
{};
template<typename MatrixType, int Mode> class SparseTriangularView
: public SparseMatrixBase<SparseTriangularView<MatrixType,Mode> >
{
enum { SkipFirst = (Mode==LowerTriangular && (!MatrixType::Flags&RowMajorBit))
|| (Mode==UpperTriangular && ( MatrixType::Flags&RowMajorBit)) };
public:
class InnerIterator;
inline int rows() { return m_matrix.rows(); }
inline int cols() { return m_matrix.cols(); }
typedef typename ei_traits<MatrixType>::Scalar Scalar;
typedef typename ei_meta_if<ei_must_nest_by_value<MatrixType>::ret,
MatrixType, const MatrixType&>::ret MatrixTypeNested;
inline SparseTriangularView(const MatrixType& matrix) : m_matrix(matrix) {}
/** \internal */
inline const MatrixType& nestedExpression() const { return m_matrix; }
template<typename OtherDerived>
typename ei_plain_matrix_type_column_major<OtherDerived>::type
solve(const MatrixBase<OtherDerived>& other) const;
template<typename OtherDerived> void solveInPlace(MatrixBase<OtherDerived>& other) const;
template<typename OtherDerived> void solveInPlace(SparseMatrixBase<OtherDerived>& other) const;
protected:
MatrixTypeNested m_matrix;
};
template<typename MatrixType, int Mode>
class SparseTriangularView<MatrixType,Mode>::InnerIterator : public MatrixType::InnerIterator
{
typedef typename MatrixType::InnerIterator Base;
public:
EIGEN_STRONG_INLINE InnerIterator(const SparseTriangularView& view, int outer)
: Base(view.nestedExpression(), outer)
{
if(SkipFirst)
while((*this) && this->index()<outer)
++(*this);
}
inline int row() const { return Base::row(); }
inline int col() const { return Base::col(); }
EIGEN_STRONG_INLINE operator bool() const
{
return SkipFirst ? Base::operator bool() : (Base::operator bool() && this->index() < this->outer());
}
};
template<typename Derived>
template<int Mode>
inline const SparseTriangularView<Derived, Mode>
SparseMatrixBase<Derived>::triangularView() const
{
return derived();
}
#endif // EIGEN_SPARSE_TRIANGULARVIEW_H

28
Eigen/src/Sparse/SparseUtil.h
View File

@ -118,25 +118,29 @@ enum {
};
template<typename Derived> class SparseMatrixBase;
template<typename _Scalar, int _Flags = 0> class SparseMatrix;
template<typename _Scalar, int _Flags = 0> class DynamicSparseMatrix;
template<typename _Scalar, int _Flags = 0> class SparseVector;
template<typename _Scalar, int _Flags = 0> class MappedSparseMatrix;
template<typename _Scalar, int _Flags = 0> class SparseMatrix;
template<typename _Scalar, int _Flags = 0> class DynamicSparseMatrix;
template<typename _Scalar, int _Flags = 0> class SparseVector;
template<typename _Scalar, int _Flags = 0> class MappedSparseMatrix;
template<typename MatrixType> class SparseNestByValue;
template<typename MatrixType, int Size> class SparseInnerVectorSet;
template<typename ExpressionType,
unsigned int Added, unsigned int Removed> class SparseFlagged;
template<typename ExpressionType, int Mode> class SparseTriangular;
template<typename Lhs, typename Rhs> class SparseDiagonalProduct;
template<typename MatrixType> class SparseNestByValue;
template<typename MatrixType, int Size> class SparseInnerVectorSet;
template<typename MatrixType, int Mode> class SparseTriangularView;
template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView;
template<typename Lhs, typename Rhs> class SparseDiagonalProduct;
template<typename Lhs, typename Rhs> class SparseProduct;
template<typename Lhs, typename Rhs> class SparseTimeDenseProduct;
template<typename Lhs, typename Rhs> class DenseTimeSparseProduct;
template<typename Lhs, typename Rhs,
typename LhsStorage = typename ei_traits<Lhs>::StorageType,
typename RhsStorage = typename ei_traits<Rhs>::StorageType> struct ei_sparse_product_mode;
template<typename Lhs, typename Rhs, int ProductMode = ei_sparse_product_mode<Lhs,Rhs>::value> struct SparseProductReturnType;
template<typename Lhs, typename Rhs> struct SparseProductReturnType;
const int CoherentAccessPattern = 0x1;
const int CoherentAccessPattern = 0x1;
const int InnerRandomAccessPattern = 0x2 | CoherentAccessPattern;
const int OuterRandomAccessPattern = 0x4 | CoherentAccessPattern;
const int RandomAccessPattern = 0x8 | OuterRandomAccessPattern | InnerRandomAccessPattern;

18
Eigen/src/Sparse/TriangularSolver.h
View File

@ -160,7 +160,7 @@ struct ei_sparse_solve_triangular_selector<Lhs,Rhs,Mode,UpperTriangular,ColMajor
template<typename ExpressionType,int Mode>
template<typename OtherDerived>
void SparseTriangular<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDerived>& other) const
void SparseTriangularView<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDerived>& other) const
{
ei_assert(m_matrix.cols() == m_matrix.rows());
ei_assert(m_matrix.cols() == other.rows());
@ -182,7 +182,7 @@ void SparseTriangular<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDerived
template<typename ExpressionType,int Mode>
template<typename OtherDerived>
typename ei_plain_matrix_type_column_major<OtherDerived>::type
SparseTriangular<ExpressionType,Mode>::solve(const MatrixBase<OtherDerived>& other) const
SparseTriangularView<ExpressionType,Mode>::solve(const MatrixBase<OtherDerived>& other) const
{
typename ei_plain_matrix_type_column_major<OtherDerived>::type res(other);
solveInPlace(res);
@ -210,10 +210,10 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
const bool IsLowerTriangular = (UpLo==LowerTriangular);
AmbiVector<Scalar> tempVector(other.rows()*2);
tempVector.setBounds(0,other.rows());
Rhs res(other.rows(), other.cols());
res.reserve(other.nonZeros());
for(int col=0 ; col<other.cols() ; ++col)
{
// FIXME estimate number of non zeros
@ -224,7 +224,7 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
{
tempVector.coeffRef(rhsIt.index()) = rhsIt.value();
}
for(int i=IsLowerTriangular?0:lhs.cols()-1;
IsLowerTriangular?i<lhs.cols():i>=0;
i+=IsLowerTriangular?1:-1)
@ -233,7 +233,7 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
Scalar& ci = tempVector.coeffRef(i);
if (ci!=Scalar(0))
{
// find
// find
typename Lhs::InnerIterator it(lhs, i);
if(!(Mode & UnitDiagBit))
{
@ -260,8 +260,8 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
}
}
}
int count = 0;
// FIXME compute a reference value to filter zeros
for (typename AmbiVector<Scalar>::Iterator it(tempVector/*,1e-12*/); it; ++it)
@ -281,7 +281,7 @@ struct ei_sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
template<typename ExpressionType,int Mode>
template<typename OtherDerived>
void SparseTriangular<ExpressionType,Mode>::solveInPlace(SparseMatrixBase<OtherDerived>& other) const
void SparseTriangularView<ExpressionType,Mode>::solveInPlace(SparseMatrixBase<OtherDerived>& other) const
{
ei_assert(m_matrix.cols() == m_matrix.rows());
ei_assert(m_matrix.cols() == other.rows());

8
test/sparse_product.cpp
View File

@ -114,10 +114,10 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType&
VERIFY_IS_APPROX(mS.transpose().conjugate(), mS);
VERIFY_IS_APPROX(mS, refS);
VERIFY_IS_APPROX(x=mS*b, refX=refS*b);
// TODO properly implement triangular/selfadjoint views
// 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);
VERIFY_IS_APPROX(x=mUp.template selfadjointView<UpperTriangular>()*b, refX=refS*b);
VERIFY_IS_APPROX(x=mLo.template selfadjointView<LowerTriangular>()*b, refX=refS*b);
VERIFY_IS_APPROX(x=mS.template selfadjointView<UpperTriangular|LowerTriangular>()*b, refX=refS*b);
}
}

10
test/sparse_solvers.cpp
View File

@ -66,12 +66,12 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
// lower - dense
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
VERIFY_IS_APPROX(refMat2.template triangularView<LowerTriangular>().solve(vec2),
m2.template triangular<LowerTriangular>().solve(vec3));
m2.template triangularView<LowerTriangular>().solve(vec3));
// upper - dense
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
VERIFY_IS_APPROX(refMat2.template triangularView<UpperTriangular>().solve(vec2),
m2.template triangular<UpperTriangular>().solve(vec3));
m2.template triangularView<UpperTriangular>().solve(vec3));
// TODO test row major
@ -82,20 +82,20 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular);
initSparse<Scalar>(density, refMatB, matB);
refMat2.template triangularView<LowerTriangular>().solveInPlace(refMatB);
m2.template triangular<LowerTriangular>().solveInPlace(matB);
m2.template triangularView<LowerTriangular>().solveInPlace(matB);
VERIFY_IS_APPROX(matB.toDense(), refMatB);
// upper - sparse
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular);
initSparse<Scalar>(density, refMatB, matB);
refMat2.template triangularView<UpperTriangular>().solveInPlace(refMatB);
m2.template triangular<UpperTriangular>().solveInPlace(matB);
m2.template triangularView<UpperTriangular>().solveInPlace(matB);
VERIFY_IS_APPROX(matB, refMatB);
// test deprecated API
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
VERIFY_IS_APPROX(refMat2.template triangularView<LowerTriangular>().solve(vec2),
m2.template triangular<LowerTriangular>().solve(vec3));
m2.template triangularView<LowerTriangular>().solve(vec3));
}
// test LLT