Fix bug #770: workaround thread safety in mpreal

This commit is contained in:
Gael Guennebaud 2014-07-17 12:00:56 +02:00
parent 40b74411e4
commit 84ad8ce7e3

View File

@ -146,8 +146,8 @@ int main()
}; };
}; };
template<typename Index, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs> template<typename Index, bool ConjugateLhs, bool ConjugateRhs>
struct gebp_kernel<mpfr::mpreal,mpfr::mpreal,Index,mr,nr,ConjugateLhs,ConjugateRhs> struct gebp_kernel<mpfr::mpreal,mpfr::mpreal,Index,1,1,ConjugateLhs,ConjugateRhs>
{ {
typedef mpfr::mpreal mpreal; typedef mpfr::mpreal mpreal;
@ -155,34 +155,34 @@ int main()
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) Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0)
{ {
mpreal acc1, tmp; mpreal acc1(0,mpfr_get_prec(blockA[0].mpfr_srcptr())),
tmp (0,mpfr_get_prec(blockA[0].mpfr_srcptr()));
if(strideA==-1) strideA = depth; if(strideA==-1) strideA = depth;
if(strideB==-1) strideB = depth; if(strideB==-1) strideB = depth;
for(Index j=0; j<cols; j+=nr) for(Index i=0; i<rows; ++i)
{ {
Index actual_nr = (std::min<Index>)(nr,cols-j); for(Index j=0; j<cols; ++j)
mpreal *C1 = res + j*resStride;
for(Index i=0; i<rows; i++)
{ {
mpreal *B = const_cast<mpreal*>(blockB) + j*strideB + offsetB*actual_nr; mpreal *C1 = res + j*resStride;
mpreal *A = const_cast<mpreal*>(blockA) + i*strideA + offsetA;
const mpreal *A = blockA + i*strideA + offsetA;
const mpreal *B = blockB + j*strideB + offsetB;
acc1 = 0; acc1 = 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_srcptr(), B[k].mpfr_srcptr(), 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());
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_srcptr(), alpha.mpfr_srcptr(), 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_srcptr(), acc1.mpfr_srcptr(), mpreal::get_default_rnd());
} }
} }
} }
}; };
} // end namespace internal } // end namespace internal
} }