diff options
Diffstat (limited to 'src/EigenUnsupported/MPRealSupport')
-rw-r--r-- | src/EigenUnsupported/MPRealSupport | 213 |
1 files changed, 0 insertions, 213 deletions
diff --git a/src/EigenUnsupported/MPRealSupport b/src/EigenUnsupported/MPRealSupport deleted file mode 100644 index c4ea4ec..0000000 --- a/src/EigenUnsupported/MPRealSupport +++ /dev/null @@ -1,213 +0,0 @@ -// 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>. - * - * \warning MPFR C++ is licensed under the GPL. - * - * 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 = HugeCost, - AddCost = HugeCost, - MulCost = HugeCost - }; - - typedef mpfr::mpreal Real; - typedef mpfr::mpreal NonInteger; - - static inline Real highest (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::maxval(Precision); } - static inline Real lowest (long Precision = mpfr::mpreal::get_default_prec()) { return -mpfr::maxval(Precision); } - - // Constants - static inline Real Pi (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_pi(Precision); } - static inline Real Euler (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_euler(Precision); } - static inline Real Log2 (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_log2(Precision); } - static inline Real Catalan (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::const_catalan(Precision); } - - static inline Real epsilon (long Precision = mpfr::mpreal::get_default_prec()) { return mpfr::machine_epsilon(Precision); } - static inline Real epsilon (const Real& x) { return mpfr::machine_epsilon(x); } - -#ifdef MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS - static inline int digits10 (long Precision = mpfr::mpreal::get_default_prec()) { return std::numeric_limits<Real>::digits10(Precision); } - static inline int digits10 (const Real& x) { return std::numeric_limits<Real>::digits10(x); } - - static inline int digits () { return std::numeric_limits<Real>::digits(); } - static inline int digits (const Real& x) { return std::numeric_limits<Real>::digits(x); } -#endif - - static inline Real dummy_precision() - { - mpfr_prec_t 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 { - Vectorizable = false, - LhsPacketSize = 1, - RhsPacketSize = 1, - ResPacketSize = 1, - NumberOfRegisters = 1, - nr = 1, - mr = 1, - LhsProgress = 1, - RhsProgress = 1 - }; - typedef ResScalar LhsPacket; - typedef ResScalar RhsPacket; - typedef ResScalar ResPacket; - typedef LhsPacket LhsPacket4Packing; - - }; - - - - template<typename Index, typename DataMapper, bool ConjugateLhs, bool ConjugateRhs> - struct gebp_kernel<mpfr::mpreal,mpfr::mpreal,Index,DataMapper,1,1,ConjugateLhs,ConjugateRhs> - { - typedef mpfr::mpreal mpreal; - - EIGEN_DONT_INLINE - void operator()(const DataMapper& res, const mpreal* blockA, const mpreal* blockB, - Index rows, Index depth, Index cols, const mpreal& alpha, - Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0) - { - if(rows==0 || cols==0 || depth==0) - return; - - 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(strideB==-1) strideB = depth; - - for(Index i=0; i<rows; ++i) - { - for(Index j=0; j<cols; ++j) - { - const mpreal *A = blockA + i*strideA + offsetA; - const mpreal *B = blockB + j*strideB + offsetB; - - acc1 = 0; - for(Index k=0; k<depth; k++) - { - 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_mul(acc1.mpfr_ptr(), acc1.mpfr_srcptr(), alpha.mpfr_srcptr(), mpreal::get_default_rnd()); - mpfr_add(res(i,j).mpfr_ptr(), res(i,j).mpfr_srcptr(), acc1.mpfr_srcptr(), mpreal::get_default_rnd()); - } - } - } - }; - } // end namespace internal -} - -#endif // EIGEN_MPREALSUPPORT_MODULE_H |