mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-11-27 06:30:28 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix 4x4 inverse via SSE for submatrices

Browse Source
This commit is contained in:
Gael Guennebaud 2014-07-31 16:24:29 +02:00
parent db183ca7b3
commit 26d2cdefd4
3 changed files with 19 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
Eigen/src/LU/arch/Inverse_SSE.h
View File

@ -39,9 +39,11 @@ struct compute_inverse_size4<Architecture::SSE, float, MatrixType, ResultType>
ResultAlignment = bool(ResultType::Flags&AlignedBit),
StorageOrdersMatch = (MatrixType::Flags&RowMajorBit) == (ResultType::Flags&RowMajorBit)
};
typedef typename conditional<(MatrixType::Flags&LinearAccessBit),MatrixType const &,typename MatrixType::PlainObject>::type ActualMatrixType;
static void run(const MatrixType& matrix, ResultType& result)
static void run(const MatrixType& mat, ResultType& result)
{
ActualMatrixType matrix(mat);
EIGEN_ALIGN16 const unsigned int _Sign_PNNP[4] = { 0x00000000, 0x80000000, 0x80000000, 0x00000000 };
// Load the full matrix into registers
@ -167,14 +169,17 @@ struct compute_inverse_size4<Architecture::SSE, double, MatrixType, ResultType>
ResultAlignment = bool(ResultType::Flags&AlignedBit),
StorageOrdersMatch = (MatrixType::Flags&RowMajorBit) == (ResultType::Flags&RowMajorBit)
};
static void run(const MatrixType& matrix, ResultType& result)
typedef typename conditional<(MatrixType::Flags&LinearAccessBit),MatrixType const &,typename MatrixType::PlainObject>::type ActualMatrixType;
static void run(const MatrixType& mat, ResultType& result)
{
ActualMatrixType matrix(mat);
const __m128d _Sign_NP = _mm_castsi128_pd(_mm_set_epi32(0x0,0x0,0x80000000,0x0));
const __m128d _Sign_PN = _mm_castsi128_pd(_mm_set_epi32(0x80000000,0x0,0x0,0x0));
// The inverse is calculated using "Divide and Conquer" technique. The
// original matrix is divide into four 2x2 sub-matrices. Since each
// register of the matrix holds two element, the smaller matrices are
// register of the matrix holds two elements, the smaller matrices are
// consisted of two registers. Hence we get a better locality of the
// calculations.

8
test/CMakeLists.txt
View File

@ -220,11 +220,9 @@ ei_add_test(geo_quaternion)
ei_add_test(geo_eulerangles)
ei_add_test(geo_parametrizedline)
ei_add_test(geo_alignedbox)
if(NOT EIGEN_TEST_EVALUATORS)
ei_add_test(geo_hyperplane)
ei_add_test(geo_transformations)
ei_add_test(geo_homogeneous)
endif(NOT EIGEN_TEST_EVALUATORS)
ei_add_test(geo_hyperplane)
ei_add_test(geo_transformations)
ei_add_test(geo_homogeneous)
ei_add_test(stdvector)
ei_add_test(stdvector_overload)
ei_add_test(stdlist)

8
test/inverse.cpp
View File

@ -68,6 +68,14 @@ template<typename MatrixType> void inverse(const MatrixType& m)
VERIFY_IS_MUCH_SMALLER_THAN(abs(det-m3.determinant()), RealScalar(1));
m3.computeInverseWithCheck(m4, invertible);
VERIFY( rows==1 ? invertible : !invertible );
// check with submatrices
{
Matrix<Scalar, MatrixType::RowsAtCompileTime+1, MatrixType::RowsAtCompileTime+1, MatrixType::Options> m3;
m3.setRandom();
m2 = m3.template topLeftCorner<MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime>().inverse();
VERIFY_IS_APPROX( (m3.template topLeftCorner<MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime>()), m2.inverse() );
}
#endif
// check in-place inversion