mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-15 07:10:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug #842: fix specialized product for mpreal

Browse Source
This commit is contained in:
Gael Guennebaud 2014-07-17 09:41:33 +02:00
parent a53f2b0e43
commit 424c3ad266
2 changed files with 5 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
unsupported/Eigen/MPRealSupport
View File

@ -139,9 +139,8 @@ int main()
public: public:
typedef mpfr::mpreal ResScalar; typedef mpfr::mpreal ResScalar;
enum { enum {
nr = 2, // must be 2 for proper packing... nr = 1,
mr = 1, mr = 1,
WorkSpaceFactor = nr,
LhsProgress = 1, LhsProgress = 1,
RhsProgress = 1 RhsProgress = 1
}; };
@ -154,9 +153,9 @@ int main()
EIGEN_DONT_INLINE EIGEN_DONT_INLINE
void operator()(mpreal* res, Index resStride, const mpreal* blockA, const mpreal* blockB, Index rows, Index depth, Index cols, mpreal alpha, void operator()(mpreal* res, Index resStride, const mpreal* blockA, const mpreal* blockB, Index rows, Index depth, Index cols, mpreal alpha,
Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0, mpreal* /*unpackedB*/ = 0) Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0)
{ {
mpreal acc1, acc2, tmp; mpreal acc1, tmp;
if(strideA==-1) strideA = depth; if(strideA==-1) strideA = depth;
if(strideB==-1) strideB = depth; if(strideB==-1) strideB = depth;
@ -165,33 +164,20 @@ int main()
{ {
Index actual_nr = (std::min<Index>)(nr,cols-j); Index actual_nr = (std::min<Index>)(nr,cols-j);
mpreal *C1 = res + j*resStride; mpreal *C1 = res + j*resStride;
mpreal *C2 = res + (j+1)*resStride;
for(Index i=0; i<rows; i++) for(Index i=0; i<rows; i++)
{ {
mpreal *B = const_cast<mpreal*>(blockB) + j*strideB + offsetB*actual_nr; mpreal *B = const_cast<mpreal*>(blockB) + j*strideB + offsetB*actual_nr;
mpreal *A = const_cast<mpreal*>(blockA) + i*strideA + offsetA; mpreal *A = const_cast<mpreal*>(blockA) + i*strideA + offsetA;
acc1 = 0; acc1 = 0;
acc2 = 0;
for(Index k=0; k<depth; k++) for(Index k=0; k<depth; k++)
{ {
mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_ptr(), B[0].mpfr_ptr(), mpreal::get_default_rnd()); mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_ptr(), B[0].mpfr_ptr(), mpreal::get_default_rnd());
mpfr_add(acc1.mpfr_ptr(), acc1.mpfr_ptr(), tmp.mpfr_ptr(), mpreal::get_default_rnd()); mpfr_add(acc1.mpfr_ptr(), acc1.mpfr_ptr(), tmp.mpfr_ptr(), mpreal::get_default_rnd());
if(actual_nr==2) {
mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_ptr(), B[1].mpfr_ptr(), mpreal::get_default_rnd());
mpfr_add(acc2.mpfr_ptr(), acc2.mpfr_ptr(), tmp.mpfr_ptr(), mpreal::get_default_rnd());
}
B+=actual_nr; B+=actual_nr;
} }
mpfr_mul(acc1.mpfr_ptr(), acc1.mpfr_ptr(), alpha.mpfr_ptr(), mpreal::get_default_rnd()); mpfr_mul(acc1.mpfr_ptr(), acc1.mpfr_ptr(), alpha.mpfr_ptr(), mpreal::get_default_rnd());
mpfr_add(C1[i].mpfr_ptr(), C1[i].mpfr_ptr(), acc1.mpfr_ptr(), mpreal::get_default_rnd()); mpfr_add(C1[i].mpfr_ptr(), C1[i].mpfr_ptr(), acc1.mpfr_ptr(), mpreal::get_default_rnd());
if(actual_nr==2) {
mpfr_mul(acc2.mpfr_ptr(), acc2.mpfr_ptr(), alpha.mpfr_ptr(), mpreal::get_default_rnd());
mpfr_add(C2[i].mpfr_ptr(), C2[i].mpfr_ptr(), acc2.mpfr_ptr(), mpreal::get_default_rnd());
}
} }
} }
} }

4
unsupported/test/mpreal_support.cpp
View File

@ -1,3 +1,4 @@
#include <bench/BenchTimer.h>
#include "main.h" #include "main.h"
#include <Eigen/MPRealSupport> #include <Eigen/MPRealSupport>
#include <Eigen/LU> #include <Eigen/LU>
@ -32,12 +33,11 @@ void test_mpreal_support()
VERIFY_IS_APPROX(A.array().abs2().sqrt(), A.array().abs()); VERIFY_IS_APPROX(A.array().abs2().sqrt(), A.array().abs());
VERIFY_IS_APPROX(A.array().sin(), sin(A.array())); VERIFY_IS_APPROX(A.array().sin(), sin(A.array()));
VERIFY_IS_APPROX(A.array().cos(), cos(A.array())); VERIFY_IS_APPROX(A.array().cos(), cos(A.array()));
// Cholesky // Cholesky
X = S.selfadjointView<Lower>().llt().solve(B); X = S.selfadjointView<Lower>().llt().solve(B);
VERIFY_IS_APPROX((S.selfadjointView<Lower>()*X).eval(),B); VERIFY_IS_APPROX((S.selfadjointView<Lower>()*X).eval(),B);
// partial LU // partial LU
X = A.lu().solve(B); X = A.lu().solve(B);
VERIFY_IS_APPROX((A*X).eval(),B); VERIFY_IS_APPROX((A*X).eval(),B);