devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

finally directly calling the low-level products is faster

Browse Source
This commit is contained in:
Gael Guennebaud
2009-07-10 10:41:26 +02:00
parent 8885d56928
commit 1a1b2e9f27
12 changed files with 111 additions and 146 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
Eigen/src/Core/GenericPacketMath.h
View File

@@ -245,5 +245,15 @@ inline void ei_palign(PacketType& first, const PacketType& second)
ei_palign_impl<Offset,PacketType>::run(first,second);
}
/***************************************************************************
* Fast complex products (GCC generates a function call which is very slow)
***************************************************************************/
template<> inline std::complex<float> ei_pmul(const std::complex<float>& a, const std::complex<float>& b)
{ return std::complex<float>(ei_real(a)*ei_real(b) - ei_imag(a)*ei_imag(b), ei_imag(a)*ei_real(b) + ei_real(a)*ei_imag(b)); }
template<> inline std::complex<double> ei_pmul(const std::complex<double>& a, const std::complex<double>& b)
{ return std::complex<double>(ei_real(a)*ei_real(b) - ei_imag(a)*ei_imag(b), ei_imag(a)*ei_real(b) + ei_real(a)*ei_imag(b)); }
#endif // EIGEN_GENERIC_PACKET_MATH_H

3
Eigen/src/Core/Product.h
View File

@@ -281,8 +281,7 @@ template<typename LhsNested, typename RhsNested, int ProductMode> class Product
*/
EIGEN_STRONG_INLINE bool _useCacheFriendlyProduct() const
{
#define EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD 16
// TODO do something more accurate here
// TODO do something more accurate here (especially for mat-vec products)
return m_lhs.cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
&& ( rows()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
|| cols()>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD);

81
Eigen/src/Core/SolveTriangular.h
View File

@@ -27,7 +27,6 @@
template<typename Lhs, typename Rhs,
int Mode, // Upper/Lower | UnitDiag
// bool ConjugateLhs, bool ConjugateRhs,
int UpLo = (Mode & LowerTriangularBit)
? LowerTriangular
: (Mode & UpperTriangularBit)
@@ -38,15 +37,20 @@ template<typename Lhs, typename Rhs,
struct ei_triangular_solver_selector;
// forward substitution, row-major
template<typename Lhs, typename Rhs, int Mode, /*bool ConjugateLhs, bool ConjugateRhs,*/ int UpLo>
struct ei_triangular_solver_selector<Lhs,Rhs,Mode,/*ConjugateLhs,ConjugateRhs,*/UpLo,RowMajor>
template<typename Lhs, typename Rhs, int Mode, int UpLo>
struct ei_triangular_solver_selector<Lhs,Rhs,Mode,UpLo,RowMajor>
{
typedef typename Rhs::Scalar Scalar;
typedef ei_product_factor_traits<Lhs> LhsProductTraits;
typedef typename LhsProductTraits::ActualXprType ActualLhsType;
enum {
IsLowerTriangular = (UpLo==LowerTriangular)
};
static void run(const Lhs& lhs, Rhs& other)
{//std::cerr << "row maj " << ConjugateLhs << " , " << ConjugateRhs
// << " " << typeid(Lhs).name() << "\n";
static const int PanelWidth = 40; // TODO make this a user definable constant
static const bool IsLowerTriangular = (UpLo==LowerTriangular);
{//std::cerr << "row maj " << LhsProductTraits::NeedToConjugate << "\n";
static const int PanelWidth = EIGEN_TUNE_TRSV_PANEL_WIDTH;
const ActualLhsType& actualLhs = LhsProductTraits::extract(lhs);
const int size = lhs.cols();
for(int c=0 ; c<other.cols() ; ++c)
{
@@ -61,15 +65,12 @@ struct ei_triangular_solver_selector<Lhs,Rhs,Mode,/*ConjugateLhs,ConjugateRhs,*/
{
int startRow = IsLowerTriangular ? pi : pi-actualPanelWidth;
int startCol = IsLowerTriangular ? 0 : pi;
// Block<Rhs,Dynamic,1> target(other,startRow,c,actualPanelWidth,1);
// ei_cache_friendly_product_rowmajor_times_vector<ConjugateLhs,ConjugateRhs>(
// &(lhs.const_cast_derived().coeffRef(startRow,startCol)), lhs.stride(),
// &(other.coeffRef(startCol, c)), r,
// target, Scalar(-1));
other.col(c).segment(startRow,actualPanelWidth) -=
lhs.block(startRow,startCol,actualPanelWidth,r)
* other.col(c).segment(startCol,r);
Block<Rhs,Dynamic,1> target(other,startRow,c,actualPanelWidth,1);
ei_cache_friendly_product_rowmajor_times_vector<LhsProductTraits::NeedToConjugate,false>(
&(actualLhs.const_cast_derived().coeffRef(startRow,startCol)), actualLhs.stride(),
&(other.coeffRef(startCol, c)), r,
target, Scalar(-1));
}
for(int k=0; k<actualPanelWidth; ++k)
@@ -83,7 +84,6 @@ struct ei_triangular_solver_selector<Lhs,Rhs,Mode,/*ConjugateLhs,ConjugateRhs,*/
if(!(Mode & UnitDiagBit))
other.coeffRef(i,c) /= lhs.coeff(i,i);
}
}
}
}
@@ -94,17 +94,23 @@ struct ei_triangular_solver_selector<Lhs,Rhs,Mode,/*ConjugateLhs,ConjugateRhs,*/
// - inv(LowerTriangular,UnitDiag,ColMajor) * Column vector
// - inv(UpperTriangular, ColMajor) * Column vector
// - inv(UpperTriangular,UnitDiag,ColMajor) * Column vector
template<typename Lhs, typename Rhs, int Mode, /*bool ConjugateLhs, bool ConjugateRhs,*/ int UpLo>
struct ei_triangular_solver_selector<Lhs,Rhs,Mode,/*ConjugateLhs,ConjugateRhs,*/UpLo,ColMajor>
template<typename Lhs, typename Rhs, int Mode, int UpLo>
struct ei_triangular_solver_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
{
typedef typename Rhs::Scalar Scalar;
typedef typename ei_packet_traits<Scalar>::type Packet;
enum { PacketSize = ei_packet_traits<Scalar>::size };
typedef ei_product_factor_traits<Lhs> LhsProductTraits;
typedef typename LhsProductTraits::ActualXprType ActualLhsType;
enum {
PacketSize = ei_packet_traits<Scalar>::size,
IsLowerTriangular = (UpLo==LowerTriangular)
};
static void run(const Lhs& lhs, Rhs& other)
{//std::cerr << "col maj " << ConjugateLhs << " , " << ConjugateRhs << "\n";
static const int PanelWidth = 4; // TODO make this a user definable constant
static const bool IsLowerTriangular = (UpLo==LowerTriangular);
{//std::cerr << "col maj " << LhsProductTraits::NeedToConjugate << "\n";
static const int PanelWidth = EIGEN_TUNE_TRSV_PANEL_WIDTH;
const ActualLhsType& actualLhs = LhsProductTraits::extract(lhs);
const int size = lhs.cols();
for(int c=0 ; c<other.cols() ; ++c)
{
@@ -133,16 +139,15 @@ struct ei_triangular_solver_selector<Lhs,Rhs,Mode,/*ConjugateLhs,ConjugateRhs,*/
int r = IsLowerTriangular ? size - endBlock : startBlock; // remaining size
if (r > 0)
{
// ei_cache_friendly_product_colmajor_times_vector<ConjugateLhs,ConjugateRhs>(
// r,
// &(lhs.const_cast_derived().coeffRef(endBlock,startBlock)), lhs.stride(),
// other.col(c).segment(startBlock, actualPanelWidth),
// &(other.coeffRef(endBlock, c)),
// Scalar(-1));
other.col(c).segment(endBlock,r) -=
lhs.block(endBlock,startBlock,r,actualPanelWidth)
* other.col(c).segment(startBlock,actualPanelWidth);
// let's directly call this function because:
// 1 - it is faster to compile
// 2 - it is slighlty faster at runtime
ei_cache_friendly_product_colmajor_times_vector<LhsProductTraits::NeedToConjugate,false>(
r,
&(actualLhs.const_cast_derived().coeffRef(endBlock,startBlock)), actualLhs.stride(),
other.col(c).segment(startBlock, actualPanelWidth),
&(other.coeffRef(endBlock, c)),
Scalar(-1));
}
}
}
@@ -168,21 +173,13 @@ void TriangularView<MatrixType,Mode>::solveInPlace(const MatrixBase<RhsDerived>&
ei_assert(!(Mode & ZeroDiagBit));
ei_assert(Mode & (UpperTriangularBit|LowerTriangularBit));
// typedef ei_product_factor_traits<MatrixType> LhsProductTraits;
// typedef ei_product_factor_traits<RhsDerived> RhsProductTraits;
// typedef typename LhsProductTraits::ActualXprType ActualLhsType;
// typedef typename RhsProductTraits::ActualXprType ActualRhsType;
// const ActualLhsType& actualLhs = LhsProductTraits::extract(_expression());
// ActualRhsType& actualRhs = const_cast<ActualRhsType&>(RhsProductTraits::extract(rhs));
enum { copy = ei_traits<RhsDerived>::Flags & RowMajorBit };
// std::cerr << typeid(MatrixType).name() << "\n";
typedef typename ei_meta_if<copy,
typename ei_plain_matrix_type_column_major<RhsDerived>::type, RhsDerived&>::ret RhsCopy;
RhsCopy rhsCopy(rhs);
ei_triangular_solver_selector<MatrixType, typename ei_unref<RhsCopy>::type,
Mode/*, LhsProductTraits::NeedToConjugate,RhsProductTraits::NeedToConjugate*/>::run(_expression(), rhsCopy);
Mode>::run(_expression(), rhsCopy);
if (copy)
rhs = rhsCopy;

2
Eigen/src/Core/arch/SSE/PacketMath.h
View File

@@ -26,7 +26,7 @@
#define EIGEN_PACKET_MATH_SSE_H
#ifndef EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD
#define EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD 16
#define EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD 8
#endif
typedef __m128 Packet4f;

6
Eigen/src/Core/products/GeneralMatrixMatrix.h
View File

@@ -47,8 +47,7 @@ template<> struct ei_conj_helper<false,true>
{ return c + pmul(x,y); }
template<typename T> std::complex<T> pmul(const std::complex<T>& x, const std::complex<T>& y) const
//{ return std::complex<T>(ei_real(x)*ei_real(y) + ei_imag(x)*ei_imag(y), ei_imag(x)*ei_real(y) - ei_real(x)*ei_imag(y)); }
{ return x * ei_conj(y); }
{ return std::complex<T>(ei_real(x)*ei_real(y) + ei_imag(x)*ei_imag(y), ei_imag(x)*ei_real(y) - ei_real(x)*ei_imag(y)); }
};
template<> struct ei_conj_helper<true,false>
@@ -68,8 +67,7 @@ template<> struct ei_conj_helper<true,true>
{ return c + pmul(x,y); }
template<typename T> std::complex<T> pmul(const std::complex<T>& x, const std::complex<T>& y) const
// { return std::complex<T>(ei_real(x)*ei_real(y) - ei_imag(x)*ei_imag(y), - ei_real(x)*ei_imag(y) - ei_imag(x)*ei_real(y)); }
{ return ei_conj(x) * ei_conj(y); }
{ return std::complex<T>(ei_real(x)*ei_real(y) - ei_imag(x)*ei_imag(y), - ei_real(x)*ei_imag(y) - ei_imag(x)*ei_real(y)); }
};
#ifndef EIGEN_EXTERN_INSTANTIATIONS

8
Eigen/src/Core/products/SelfadjointMatrixVector.h
View File

@@ -25,14 +25,6 @@
#ifndef EIGEN_SELFADJOINT_MATRIX_VECTOR_H
#define EIGEN_SELFADJOINT_MATRIX_VECTOR_H
template<bool Conjugate> struct ei_conj_if {
template<typename Scalar> Scalar operator() (const Scalar& x) const { return ei_conj(x); }
};
template<> struct ei_conj_if<false> {
template<typename Scalar> Scalar& operator() (Scalar& x) const { return x; }
};
/* Optimized col-major selfadjoint matrix * vector product:
* This algorithm processes 2 columns at onces that allows to both reduce
* the number of load/stores of the result by a factor 2 and to reduce

7
Eigen/src/Core/util/Macros.h
View File

@@ -94,6 +94,13 @@
#define EIGEN_TUNE_FOR_CPU_CACHE_SIZE (sizeof(float)*256*256)
#endif
/** Defines the maximal width of the blocks used in the triangular solver
* for vectors (level 2 blas xTRSV). The default is 8.
*/
#ifndef EIGEN_TUNE_TRSV_PANEL_WIDTH
#define EIGEN_TUNE_TRSV_PANEL_WIDTH 8
#endif
/** Allows to disable some optimizations which might affect the accuracy of the result.
* Such optimization are enabled by default, and set EIGEN_FAST_MATH to 0 to disable them.
* They currently include:

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

@@ -198,4 +198,16 @@ template<typename T> struct ei_is_diagonal<DiagonalWrapper<T> >
template<typename T, int S> struct ei_is_diagonal<DiagonalMatrix<T,S> >
{ enum { ret = true }; };
template<bool Conjugate> struct ei_conj_if;
template<> struct ei_conj_if<true> {
template<typename T>
inline T operator()(const T& x) { return ei_conj(x); }
};
template<> struct ei_conj_if<false> {
template<typename T>
inline const T& operator()(const T& x) { return x; }
};
#endif // EIGEN_META_H