| // This file is part of a joint effort between Eigen, a lightweight C++ template library |
| // for linear algebra, and MPFR C++, a C++ interface to MPFR library (http://www.holoborodko.com/pavel/) |
| // |
| // Copyright (C) 2010-2012 Pavel Holoborodko <pavel@holoborodko.com> |
| // Copyright (C) 2010 Konstantin Holoborodko <konstantin@holoborodko.com> |
| // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr> |
| // |
| // This Source Code Form is subject to the terms of the Mozilla |
| // Public License v. 2.0. If a copy of the MPL was not distributed |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| |
| #ifndef EIGEN_MPREALSUPPORT_MODULE_H |
| #define EIGEN_MPREALSUPPORT_MODULE_H |
| |
| #include <Eigen/Core> |
| #include <mpreal.h> |
| |
| namespace Eigen { |
| |
| /** |
| * \defgroup MPRealSupport_Module MPFRC++ Support module |
| * \code |
| * #include <Eigen/MPRealSupport> |
| * \endcode |
| * |
| * This module provides support for multi precision floating point numbers |
| * via the <a href="http://www.holoborodko.com/pavel/mpfr">MPFR C++</a> |
| * library which itself is built upon <a href="http://www.mpfr.org/">MPFR</a>/<a href="http://gmplib.org/">GMP</a>. |
| * |
| * You can find a copy of MPFR C++ that is known to be compatible in the unsupported/test/mpreal folder. |
| * |
| * Here is an example: |
| * |
| \code |
| #include <iostream> |
| #include <Eigen/MPRealSupport> |
| #include <Eigen/LU> |
| using namespace mpfr; |
| using namespace Eigen; |
| int main() |
| { |
| // set precision to 256 bits (double has only 53 bits) |
| mpreal::set_default_prec(256); |
| // Declare matrix and vector types with multi-precision scalar type |
| typedef Matrix<mpreal,Dynamic,Dynamic> MatrixXmp; |
| typedef Matrix<mpreal,Dynamic,1> VectorXmp; |
| |
| MatrixXmp A = MatrixXmp::Random(100,100); |
| VectorXmp b = VectorXmp::Random(100); |
| |
| // Solve Ax=b using LU |
| VectorXmp x = A.lu().solve(b); |
| std::cout << "relative error: " << (A*x - b).norm() / b.norm() << std::endl; |
| return 0; |
| } |
| \endcode |
| * |
| */ |
| |
| template<> struct NumTraits<mpfr::mpreal> |
| : GenericNumTraits<mpfr::mpreal> |
| { |
| enum { |
| IsInteger = 0, |
| IsSigned = 1, |
| IsComplex = 0, |
| RequireInitialization = 1, |
| ReadCost = 10, |
| AddCost = 10, |
| MulCost = 40 |
| }; |
| |
| typedef mpfr::mpreal Real; |
| typedef mpfr::mpreal NonInteger; |
| |
| inline static Real highest (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::maxval(Precision); } |
| inline static Real lowest (long Precision = mpfr::mpreal::get_default_prec()) { return -mpfr::maxval(Precision); } |
| |
| // Constants |
| inline static Real Pi (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_pi(Precision); } |
| inline static Real Euler (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_euler(Precision); } |
| inline static Real Log2 (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_log2(Precision); } |
| inline static Real Catalan (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_catalan(Precision); } |
| |
| inline static Real epsilon (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::machine_epsilon(Precision); } |
| inline static Real epsilon (const Real& x) { return mpfr::machine_epsilon(x); } |
| |
| inline static Real dummy_precision() |
| { |
| unsigned int weak_prec = ((mpfr::mpreal::get_default_prec()-1) * 90) / 100; |
| return mpfr::machine_epsilon(weak_prec); |
| } |
| }; |
| |
| namespace internal { |
| |
| template<> inline mpfr::mpreal random<mpfr::mpreal>() |
| { |
| return mpfr::random(); |
| } |
| |
| template<> inline mpfr::mpreal random<mpfr::mpreal>(const mpfr::mpreal& a, const mpfr::mpreal& b) |
| { |
| return a + (b-a) * random<mpfr::mpreal>(); |
| } |
| |
| inline bool isMuchSmallerThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps) |
| { |
| return mpfr::abs(a) <= mpfr::abs(b) * eps; |
| } |
| |
| inline bool isApprox(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps) |
| { |
| return mpfr::isEqualFuzzy(a,b,eps); |
| } |
| |
| inline bool isApproxOrLessThan(const mpfr::mpreal& a, const mpfr::mpreal& b, const mpfr::mpreal& eps) |
| { |
| return a <= b || mpfr::isEqualFuzzy(a,b,eps); |
| } |
| |
| template<> inline long double cast<mpfr::mpreal,long double>(const mpfr::mpreal& x) |
| { return x.toLDouble(); } |
| |
| template<> inline double cast<mpfr::mpreal,double>(const mpfr::mpreal& x) |
| { return x.toDouble(); } |
| |
| template<> inline long cast<mpfr::mpreal,long>(const mpfr::mpreal& x) |
| { return x.toLong(); } |
| |
| template<> inline int cast<mpfr::mpreal,int>(const mpfr::mpreal& x) |
| { return int(x.toLong()); } |
| |
| // Specialize GEBP kernel and traits for mpreal (no need for peeling, nor complicated stuff) |
| // This also permits to directly call mpfr's routines and avoid many temporaries produced by mpreal |
| template<> |
| class gebp_traits<mpfr::mpreal, mpfr::mpreal, false, false> |
| { |
| public: |
| typedef mpfr::mpreal ResScalar; |
| enum { |
| nr = 2, // must be 2 for proper packing... |
| mr = 1, |
| WorkSpaceFactor = nr, |
| LhsProgress = 1, |
| RhsProgress = 1 |
| }; |
| }; |
| |
| template<typename Index, int mr, int nr, bool ConjugateLhs, bool ConjugateRhs> |
| struct gebp_kernel<mpfr::mpreal,mpfr::mpreal,Index,mr,nr,ConjugateLhs,ConjugateRhs> |
| { |
| typedef mpfr::mpreal mpreal; |
| |
| EIGEN_DONT_INLINE |
| 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) |
| { |
| mpreal acc1, acc2, tmp; |
| |
| if(strideA==-1) strideA = depth; |
| if(strideB==-1) strideB = depth; |
| |
| for(Index j=0; j<cols; j+=nr) |
| { |
| Index actual_nr = (std::min<Index>)(nr,cols-j); |
| mpreal *C1 = res + j*resStride; |
| mpreal *C2 = res + (j+1)*resStride; |
| for(Index i=0; i<rows; i++) |
| { |
| mpreal *B = const_cast<mpreal*>(blockB) + j*strideB + offsetB*actual_nr; |
| mpreal *A = const_cast<mpreal*>(blockA) + i*strideA + offsetA; |
| acc1 = 0; |
| acc2 = 0; |
| 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_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; |
| } |
| |
| 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()); |
| |
| 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()); |
| } |
| } |
| } |
| } |
| }; |
| |
| } // end namespace internal |
| } |
| |
| #endif // EIGEN_MPREALSUPPORT_MODULE_H |