GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-23 10:09:36 +08:00
Code Issues Packages Projects Releases Wiki Activity

makes trmv works with the triangular matrix on the right

Browse Source
This commit is contained in:
Gael Guennebaud 2010-06-17 10:17:22 +02:00
parent 6bff339cc5
commit 7fdf218951
2 changed files with 62 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

64
Eigen/src/Core/products/TriangularMatrixVector.h
View File

@ -25,12 +25,26 @@
#ifndef EIGEN_TRIANGULARMATRIXVECTOR_H #ifndef EIGEN_TRIANGULARMATRIXVECTOR_H
#define EIGEN_TRIANGULARMATRIXVECTOR_H #define EIGEN_TRIANGULARMATRIXVECTOR_H
template<typename MatrixType, typename Rhs, typename Result, template<bool LhsIsTriangular, typename Lhs, typename Rhs, typename Result,
int Mode, bool ConjLhs, bool ConjRhs, int StorageOrder> int Mode, bool ConjLhs, bool ConjRhs, int StorageOrder>
struct ei_product_triangular_vector_selector; struct ei_product_triangular_vector_selector;
template<typename Lhs, typename Rhs, typename Result, int Mode, bool ConjLhs, bool ConjRhs, int StorageOrder>
struct ei_product_triangular_vector_selector<false,Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs,StorageOrder>
{
static EIGEN_DONT_INLINE void run(const Lhs& lhs, const Rhs& rhs, Result& res, typename ei_traits<Lhs>::Scalar alpha)
{
typedef Transpose<Rhs> TrRhs; TrRhs trRhs(rhs);
typedef Transpose<Lhs> TrLhs; TrLhs trLhs(lhs);
typedef Transpose<Result> TrRes; TrRes trRes(res);
ei_product_triangular_vector_selector<true,TrRhs,TrLhs,TrRes,
(Mode & UnitDiag) | (Mode & Lower) ? Upper : Lower, ConjRhs, ConjLhs, StorageOrder==RowMajor ? ColMajor : RowMajor>
::run(trRhs,trLhs,trRes,alpha);
}
};
template<typename Lhs, typename Rhs, typename Result, int Mode, bool ConjLhs, bool ConjRhs> template<typename Lhs, typename Rhs, typename Result, int Mode, bool ConjLhs, bool ConjRhs>
struct ei_product_triangular_vector_selector<Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs,ColMajor> struct ei_product_triangular_vector_selector<true,Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs,ColMajor>
{ {
typedef typename Rhs::Scalar Scalar; typedef typename Rhs::Scalar Scalar;
typedef typename Rhs::Index Index; typedef typename Rhs::Index Index;
@ -74,7 +88,7 @@ struct ei_product_triangular_vector_selector<Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs
}; };
template<typename Lhs, typename Rhs, typename Result, int Mode, bool ConjLhs, bool ConjRhs> template<typename Lhs, typename Rhs, typename Result, int Mode, bool ConjLhs, bool ConjRhs>
struct ei_product_triangular_vector_selector<Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs,RowMajor> struct ei_product_triangular_vector_selector<true,Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs,RowMajor>
{ {
typedef typename Rhs::Scalar Scalar; typedef typename Rhs::Scalar Scalar;
typedef typename Rhs::Index Index; typedef typename Rhs::Index Index;
@ -119,12 +133,17 @@ struct ei_product_triangular_vector_selector<Lhs,Rhs,Result,Mode,ConjLhs,ConjRhs
* Wrapper to ei_product_triangular_vector * Wrapper to ei_product_triangular_vector
***************************************************************************/ ***************************************************************************/
template<int Mode, /*bool LhsIsTriangular, */typename Lhs, typename Rhs> template<int Mode, bool LhsIsTriangular, typename Lhs, typename Rhs>
struct ei_traits<TriangularProduct<Mode,true,Lhs,false,Rhs,true> > struct ei_traits<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,true> >
: ei_traits<ProductBase<TriangularProduct<Mode,true,Lhs,false,Rhs,true>, Lhs, Rhs> > : ei_traits<ProductBase<TriangularProduct<Mode,LhsIsTriangular,Lhs,false,Rhs,true>, Lhs, Rhs> >
{}; {};
template<int Mode, /*bool LhsIsTriangular, */typename Lhs, typename Rhs> template<int Mode, bool LhsIsTriangular, typename Lhs, typename Rhs>
struct ei_traits<TriangularProduct<Mode,LhsIsTriangular,Lhs,true,Rhs,false> >
: ei_traits<ProductBase<TriangularProduct<Mode,LhsIsTriangular,Lhs,true,Rhs,false>, Lhs, Rhs> >
{};
template<int Mode, typename Lhs, typename Rhs>
struct TriangularProduct<Mode,true,Lhs,false,Rhs,true> struct TriangularProduct<Mode,true,Lhs,false,Rhs,true>
: public ProductBase<TriangularProduct<Mode,true,Lhs,false,Rhs,true>, Lhs, Rhs > : public ProductBase<TriangularProduct<Mode,true,Lhs,false,Rhs,true>, Lhs, Rhs >
{ {
@ -143,7 +162,7 @@ struct TriangularProduct<Mode,true,Lhs,false,Rhs,true>
* RhsBlasTraits::extractScalarFactor(m_rhs); * RhsBlasTraits::extractScalarFactor(m_rhs);
ei_product_triangular_vector_selector ei_product_triangular_vector_selector
<_ActualLhsType,_ActualRhsType,Dest, <true,_ActualLhsType,_ActualRhsType,Dest,
Mode, Mode,
LhsBlasTraits::NeedToConjugate, LhsBlasTraits::NeedToConjugate,
RhsBlasTraits::NeedToConjugate, RhsBlasTraits::NeedToConjugate,
@ -152,4 +171,33 @@ struct TriangularProduct<Mode,true,Lhs,false,Rhs,true>
} }
}; };
template<int Mode, typename Lhs, typename Rhs>
struct TriangularProduct<Mode,false,Lhs,true,Rhs,false>
: public ProductBase<TriangularProduct<Mode,false,Lhs,true,Rhs,false>, Lhs, Rhs >
{
EIGEN_PRODUCT_PUBLIC_INTERFACE(TriangularProduct)
TriangularProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {}
template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const
{
ei_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols());
const ActualLhsType lhs = LhsBlasTraits::extract(m_lhs);
const ActualRhsType rhs = RhsBlasTraits::extract(m_rhs);
Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs)
* RhsBlasTraits::extractScalarFactor(m_rhs);
ei_product_triangular_vector_selector
<false,_ActualLhsType,_ActualRhsType,Dest,
Mode,
LhsBlasTraits::NeedToConjugate,
RhsBlasTraits::NeedToConjugate,
(int(ei_traits<Rhs>::Flags)&RowMajorBit) ? RowMajor : ColMajor>
::run(lhs,rhs,dst,actualAlpha);
}
};
#endif // EIGEN_TRIANGULARMATRIXVECTOR_H #endif // EIGEN_TRIANGULARMATRIXVECTOR_H

6
test/product_trmv.cpp
View File

@ -76,6 +76,12 @@ template<typename MatrixType> void trmv(const MatrixType& m)
VERIFY((m3.adjoint() * (s1*v1.conjugate())).isApprox(m1.adjoint().template triangularView<Eigen::Upper>() * (s1*v1.conjugate()), largerEps)); VERIFY((m3.adjoint() * (s1*v1.conjugate())).isApprox(m1.adjoint().template triangularView<Eigen::Upper>() * (s1*v1.conjugate()), largerEps));
m3 = m1.template triangularView<Eigen::UnitUpper>(); m3 = m1.template triangularView<Eigen::UnitUpper>();
// check transposed cases:
m3 = m1.template triangularView<Eigen::Lower>();
VERIFY((v1.transpose() * m3).isApprox(v1.transpose() * m1.template triangularView<Eigen::Lower>(), largerEps));
VERIFY((v1.adjoint() * m3).isApprox(v1.adjoint() * m1.template triangularView<Eigen::Lower>(), largerEps));
VERIFY((v1.adjoint() * m3.adjoint()).isApprox(v1.adjoint() * m1.template triangularView<Eigen::Lower>().adjoint(), largerEps));
// TODO check with sub-matrices // TODO check with sub-matrices
} }