mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-07 18:27:40 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix vectorized Jacobi Rotation

Browse Source
This commit is contained in:
Arthur 2022-08-08 19:29:56 +00:00 committed by Rasmus Munk Larsen
parent 7a87ed1b6a
commit be20207d10
2 changed files with 18 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
Eigen/src/Jacobi/Jacobi.h
View File

@ -344,15 +344,18 @@ struct apply_rotation_in_the_plane_selector<Scalar,OtherScalar,SizeAtCompileTime
{
static inline void run(Scalar *x, Index incrx, Scalar *y, Index incry, Index size, OtherScalar c, OtherScalar s)
{
enum {
PacketSize = packet_traits<Scalar>::size,
OtherPacketSize = packet_traits<OtherScalar>::size
};
typedef typename packet_traits<Scalar>::type Packet;
typedef typename packet_traits<OtherScalar>::type OtherPacket;
enum {
RequiredAlignment = plain_enum_max(unpacket_traits<Packet>::alignment,
unpacket_traits<OtherPacket>::alignment),
PacketSize = packet_traits<Scalar>::size,
OtherPacketSize = packet_traits<OtherScalar>::size
};
/*** dynamic-size vectorized paths ***/
if(SizeAtCompileTime == Dynamic && ((incrx==1 && incry==1) || PacketSize == 1))
if(size >= 2 * PacketSize && SizeAtCompileTime == Dynamic && ((incrx == 1 && incry == 1) || PacketSize == 1))
{
// both vectors are sequentially stored in memory => vectorization
enum { Peeling = 2 };
@ -423,11 +426,11 @@ struct apply_rotation_in_the_plane_selector<Scalar,OtherScalar,SizeAtCompileTime
}
/*** fixed-size vectorized path ***/
else if(SizeAtCompileTime != Dynamic && MinAlignment>0) // FIXME should be compared to the required alignment
else if(SizeAtCompileTime != Dynamic && MinAlignment >= RequiredAlignment)
{
const OtherPacket pc = pset1<OtherPacket>(c);
const OtherPacket ps = pset1<OtherPacket>(s);
conj_helper<OtherPacket,Packet,NumTraits<OtherPacket>::IsComplex,false> pcj;
conj_helper<OtherPacket,Packet,NumTraits<OtherScalar>::IsComplex,false> pcj;
conj_helper<OtherPacket,Packet,false,false> pm;
Scalar* EIGEN_RESTRICT px = x;
Scalar* EIGEN_RESTRICT py = y;
@ -452,11 +455,11 @@ struct apply_rotation_in_the_plane_selector<Scalar,OtherScalar,SizeAtCompileTime
template<typename VectorX, typename VectorY, typename OtherScalar>
EIGEN_DEVICE_FUNC
void /*EIGEN_DONT_INLINE*/ apply_rotation_in_the_plane(DenseBase<VectorX>& xpr_x, DenseBase<VectorY>& xpr_y, const JacobiRotation<OtherScalar>& j)
void inline apply_rotation_in_the_plane(DenseBase<VectorX>& xpr_x, DenseBase<VectorY>& xpr_y, const JacobiRotation<OtherScalar>& j)
{
typedef typename VectorX::Scalar Scalar;
const bool Vectorizable = (int(VectorX::Flags) & int(VectorY::Flags) & PacketAccessBit)
&& (int(packet_traits<Scalar>::size) == int(packet_traits<OtherScalar>::size));
constexpr bool Vectorizable = (int(evaluator<VectorX>::Flags) & int(evaluator<VectorY>::Flags) & PacketAccessBit) &&
(int(packet_traits<Scalar>::size) == int(packet_traits<OtherScalar>::size));
eigen_assert(xpr_x.size() == xpr_y.size());
Index size = xpr_x.size();

5
test/jacobi.cpp
View File

@ -65,6 +65,11 @@ EIGEN_DECLARE_TEST(jacobi)
CALL_SUBTEST_3(( jacobi<Matrix4cf, float>() ));
CALL_SUBTEST_3(( jacobi<Matrix4cf, std::complex<float> >() ));
CALL_SUBTEST_1(( jacobi<Matrix<float, 3, 3, RowMajor>, float>() ));
CALL_SUBTEST_2(( jacobi<Matrix<double, 4, 4, RowMajor>, double>() ));
CALL_SUBTEST_3(( jacobi<Matrix<std::complex<float>, 4, 4, RowMajor>, float>() ));
CALL_SUBTEST_3(( jacobi<Matrix<std::complex<float>, 4, 4, RowMajor>, std::complex<float> >() ));
int r = internal::random<int>(2, internal::random<int>(1,EIGEN_TEST_MAX_SIZE)/2),
c = internal::random<int>(2, internal::random<int>(1,EIGEN_TEST_MAX_SIZE)/2);
CALL_SUBTEST_4(( jacobi<MatrixXf, float>(MatrixXf(r,c)) ));