You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
209 lines
7.2 KiB
209 lines
7.2 KiB
// 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 (xxxp://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 the mozilla.org home page |
|
|
|
#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="xxxp://www.holoborodko.com/pavel/mpfr">MPFR C++</a> |
|
* library which itself is built upon <a href="xxxp://www.mpfr.org/">MPFR</a>/<a href="xxxp://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); } |
|
#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; |
|
|
|
}; |
|
|
|
|
|
|
|
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
|
|
|