GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-10 02:39:03 +08:00
Code Issues Packages Projects Releases Wiki Activity

merge

Browse Source
This commit is contained in:
Gael Guennebaud 2018-10-08 17:35:18 +02:00
commit 649d4758a6
23 changed files with 473 additions and 345 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Eigen/src/Core/DenseStorage.h
View File

@ -420,7 +420,7 @@ template<typename T, int _Options> class DenseStorage<T, Dynamic, Dynamic, Dynam
if(size != m_rows*m_cols)
{
internal::conditional_aligned_delete_auto<T,(_Options&DontAlign)==0>(m_data, m_rows*m_cols);
if (size)
if (size>0) // >0 and not simply !=0 to let the compiler knows that size cannot be negative
m_data = internal::conditional_aligned_new_auto<T,(_Options&DontAlign)==0>(size);
else
m_data = 0;
@ -497,7 +497,7 @@ template<typename T, int _Rows, int _Options> class DenseStorage<T, Dynamic, _Ro
if(size != _Rows*m_cols)
{
internal::conditional_aligned_delete_auto<T,(_Options&DontAlign)==0>(m_data, _Rows*m_cols);
if (size)
if (size>0) // >0 and not simply !=0 to let the compiler knows that size cannot be negative
m_data = internal::conditional_aligned_new_auto<T,(_Options&DontAlign)==0>(size);
else
m_data = 0;
@ -573,7 +573,7 @@ template<typename T, int _Cols, int _Options> class DenseStorage<T, Dynamic, Dyn
if(size != m_rows*_Cols)
{
internal::conditional_aligned_delete_auto<T,(_Options&DontAlign)==0>(m_data, _Cols*m_rows);
if (size)
if (size>0) // >0 and not simply !=0 to let the compiler knows that size cannot be negative
m_data = internal::conditional_aligned_new_auto<T,(_Options&DontAlign)==0>(size);
else
m_data = 0;

3
Eigen/src/Core/MathFunctions.h
View File

@ -1217,7 +1217,8 @@ inline int log2(int x)
/** \returns the square root of \a x.
*
* It is essentially equivalent to \code using std::sqrt; return sqrt(x); \endcode,
* It is essentially equivalent to
* \code using std::sqrt; return sqrt(x); \endcode
* but slightly faster for float/double and some compilers (e.g., gcc), thanks to
* specializations when SSE is enabled.
*

2
Eigen/src/Core/arch/SSE/PacketMath.h
View File

@ -28,7 +28,7 @@ namespace internal {
#endif
#endif
#if (defined EIGEN_VECTORIZE_AVX) && (EIGEN_COMP_GNUC_STRICT || EIGEN_COMP_MINGW) && (__GXX_ABI_VERSION < 1004)
#if ((defined EIGEN_VECTORIZE_AVX) && (EIGEN_COMP_GNUC_STRICT || EIGEN_COMP_MINGW) && (__GXX_ABI_VERSION < 1004)) || EIGEN_OS_QNX
// With GCC's default ABI version, a __m128 or __m256 are the same types and therefore we cannot
// have overloads for both types without linking error.
// One solution is to increase ABI version using -fabi-version=4 (or greater).

6
Eigen/src/Core/util/ConfigureVectorization.h
View File

@ -379,10 +379,12 @@
#include <cuda_fp16.h>
#endif
#if defined(EIGEN_HIP_DEVICE_COMPILE)
#if defined(EIGEN_HIPCC)
#define EIGEN_VECTORIZE_GPU
#include <hip/hip_vector_types.h>
#endif
#if defined(EIGEN_HIP_DEVICE_COMPILE)
#define EIGEN_HAS_HIP_FP16
#include <hip/hip_fp16.h>

4
Eigen/src/Core/util/IntegralConstant.h
View File

@ -55,7 +55,9 @@ public:
operator int() const { return value; }
FixedInt() {}
FixedInt( VariableAndFixedInt<N> other) {
EIGEN_ONLY_USED_FOR_DEBUG(other);
#ifndef EIGEN_INTERNAL_DEBUGGING
EIGEN_UNUSED_VARIABLE(other);
#endif
eigen_internal_assert(int(other)==N);
}

27
Eigen/src/Core/util/Memory.h
View File

@ -96,10 +96,16 @@ inline void throw_std_bad_alloc()
/** \internal Like malloc, but the returned pointer is guaranteed to be 16-byte aligned.
* Fast, but wastes 16 additional bytes of memory. Does not throw any exception.
*/
inline void* handmade_aligned_malloc(std::size_t size, std::size_t alignment = EIGEN_DEFAULT_ALIGN_BYTES)
EIGEN_DEVICE_FUNC inline void* handmade_aligned_malloc(std::size_t size, std::size_t alignment = EIGEN_DEFAULT_ALIGN_BYTES)
{
eigen_assert(alignment >= sizeof(void*) && (alignment & -alignment) == alignment && "Alignment must be at least sizeof(void*) and a power of 2");
#if defined(EIGEN_HIP_DEVICE_COMPILE)
void *original = ::malloc(size+alignment);
#else
void *original = std::malloc(size+alignment);
#endif
if (original == 0) return 0;
void *aligned = reinterpret_cast<void*>((reinterpret_cast<std::size_t>(original) & ~(std::size_t(alignment-1))) + alignment);
*(reinterpret_cast<void**>(aligned) - 1) = original;
@ -107,9 +113,15 @@ inline void* handmade_aligned_malloc(std::size_t size, std::size_t alignment = E
}
/** \internal Frees memory allocated with handmade_aligned_malloc */
inline void handmade_aligned_free(void *ptr)
EIGEN_DEVICE_FUNC inline void handmade_aligned_free(void *ptr)
{
if (ptr) std::free(*(reinterpret_cast<void**>(ptr) - 1));
if (ptr) {
#if defined(EIGEN_HIP_DEVICE_COMPILE)
::free(*(reinterpret_cast<void**>(ptr) - 1));
#else
std::free(*(reinterpret_cast<void**>(ptr) - 1));
#endif
}
}
/** \internal
@ -872,6 +884,15 @@ public:
~aligned_allocator() {}
#if EIGEN_COMP_GNUC_STRICT && EIGEN_GNUC_AT_LEAST(7,0)
// In gcc std::allocator::max_size() is bugged making gcc triggers a warning:
// eigen/Eigen/src/Core/util/Memory.h:189:12: warning: argument 1 value '18446744073709551612' exceeds maximum object size 9223372036854775807
// See https://gcc.gnu.org/bugzilla/show_bug.cgi?id=87544
size_type max_size() const {
return (std::numeric_limits<std::ptrdiff_t>::max)()/sizeof(T);
}
#endif
pointer allocate(size_type num, const void* /*hint*/ = 0)
{
internal::check_size_for_overflow<T>(num);

4
Eigen/src/plugins/ReshapedMethods.h
View File

@ -105,7 +105,7 @@ EIGEN_DEVICE_FUNC
inline Reshaped<EIGEN_RESHAPED_METHOD_CONST Derived,
internal::get_compiletime_reshape_size<NRowsType,NColsType,SizeAtCompileTime>::value,
internal::get_compiletime_reshape_size<NColsType,NRowsType,SizeAtCompileTime>::value,
Order==AutoOrder?Flags&RowMajorBit:Order>
(Order==AutoOrder?Flags&RowMajorBit:Order)>
reshaped(NRowsType nRows, NColsType nCols) EIGEN_RESHAPED_METHOD_CONST
{
return Reshaped<EIGEN_RESHAPED_METHOD_CONST Derived,
@ -128,7 +128,7 @@ reshaped() EIGEN_RESHAPED_METHOD_CONST
template<int Order>
EIGEN_DEVICE_FUNC
inline Reshaped<EIGEN_RESHAPED_METHOD_CONST Derived, SizeAtCompileTime, 1, Order==AutoOrder?Flags&RowMajorBit:Order>
inline Reshaped<EIGEN_RESHAPED_METHOD_CONST Derived, SizeAtCompileTime, 1, (Order==AutoOrder?Flags&RowMajorBit:Order)>
reshaped() EIGEN_RESHAPED_METHOD_CONST
{
EIGEN_STATIC_ASSERT(Order==RowMajor || Order==ColMajor || Order==AutoOrder, INVALID_TEMPLATE_PARAMETER);

8
test/indexed_view.cpp
View File

@ -15,6 +15,14 @@
#ifdef EIGEN_TEST_PART_3
// Make sure we also check c++98 max implementation
#define EIGEN_MAX_CPP_VER 03
// We need to disable this warning when compiling with c++11 while limiting Eigen to c++98
// Ideally we would rather configure the compiler to build in c++98 mode but this needs
// to be done at the CMakeLists.txt level.
#if defined(__GNUC__) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 8))
#pragma GCC diagnostic ignored "-Wdeprecated"
#endif
#endif
#include <valarray>

4
test/ref.cpp
View File

@ -255,8 +255,8 @@ void test_ref_overloads()
void test_ref_fixed_size_assert()
{
Vector4f v4;
VectorXf vx(10);
Vector4f v4 = Vector4f::Random();
VectorXf vx = VectorXf::Random(10);
VERIFY_RAISES_STATIC_ASSERT( Ref<Vector3f> y = v4; (void)y; );
VERIFY_RAISES_STATIC_ASSERT( Ref<Vector3f> y = vx.head<4>(); (void)y; );
VERIFY_RAISES_STATIC_ASSERT( Ref<const Vector3f> y = v4; (void)y; );

16
test/stddeque.cpp
View File

@ -18,7 +18,7 @@ void check_stddeque_matrix(const MatrixType& m)
Index rows = m.rows();
Index cols = m.cols();
MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols);
std::deque<MatrixType,Eigen::aligned_allocator<MatrixType> > v(10, MatrixType(rows,cols)), w(20, y);
std::deque<MatrixType,Eigen::aligned_allocator<MatrixType> > v(10, MatrixType::Zero(rows,cols)), w(20, y);
v.front() = x;
w.front() = w.back();
VERIFY_IS_APPROX(w.front(), w.back());
@ -33,7 +33,7 @@ void check_stddeque_matrix(const MatrixType& m)
++wi;
}
v.resize(21);
v.resize(21,MatrixType::Zero(rows,cols));
v.back() = x;
VERIFY_IS_APPROX(v.back(), x);
v.resize(22,y);
@ -46,8 +46,8 @@ template<typename TransformType>
void check_stddeque_transform(const TransformType&)
{
typedef typename TransformType::MatrixType MatrixType;
TransformType x(MatrixType::Random()), y(MatrixType::Random());
std::deque<TransformType,Eigen::aligned_allocator<TransformType> > v(10), w(20, y);
TransformType x(MatrixType::Random()), y(MatrixType::Random()), ti=TransformType::Identity();
std::deque<TransformType,Eigen::aligned_allocator<TransformType> > v(10,ti), w(20, y);
v.front() = x;
w.front() = w.back();
VERIFY_IS_APPROX(w.front(), w.back());
@ -62,7 +62,7 @@ void check_stddeque_transform(const TransformType&)
++wi;
}
v.resize(21);
v.resize(21,ti);
v.back() = x;
VERIFY_IS_APPROX(v.back(), x);
v.resize(22,y);
@ -75,8 +75,8 @@ template<typename QuaternionType>
void check_stddeque_quaternion(const QuaternionType&)
{
typedef typename QuaternionType::Coefficients Coefficients;
QuaternionType x(Coefficients::Random()), y(Coefficients::Random());
std::deque<QuaternionType,Eigen::aligned_allocator<QuaternionType> > v(10), w(20, y);
QuaternionType x(Coefficients::Random()), y(Coefficients::Random()), qi=QuaternionType::Identity();
std::deque<QuaternionType,Eigen::aligned_allocator<QuaternionType> > v(10,qi), w(20, y);
v.front() = x;
w.front() = w.back();
VERIFY_IS_APPROX(w.front(), w.back());
@ -91,7 +91,7 @@ void check_stddeque_quaternion(const QuaternionType&)
++wi;
}
v.resize(21);
v.resize(21,qi);
v.back() = x;
VERIFY_IS_APPROX(v.back(), x);
v.resize(22,y);

14
test/stddeque_overload.cpp
View File

@ -31,7 +31,7 @@ void check_stddeque_matrix(const MatrixType& m)
Index rows = m.rows();
Index cols = m.cols();
MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols);
std::deque<MatrixType> v(10, MatrixType(rows,cols)), w(20, y);
std::deque<MatrixType> v(10, MatrixType::Zero(rows,cols)), w(20, y);
v[5] = x;
w[6] = v[5];
VERIFY_IS_APPROX(w[6], v[5]);
@ -64,8 +64,8 @@ template<typename TransformType>
void check_stddeque_transform(const TransformType&)
{
typedef typename TransformType::MatrixType MatrixType;
TransformType x(MatrixType::Random()), y(MatrixType::Random());
std::deque<TransformType> v(10), w(20, y);
TransformType x(MatrixType::Random()), y(MatrixType::Random()), ti=TransformType::Identity();
std::deque<TransformType> v(10,ti), w(20, y);
v[5] = x;
w[6] = v[5];
VERIFY_IS_APPROX(w[6], v[5]);
@ -75,7 +75,7 @@ void check_stddeque_transform(const TransformType&)
VERIFY_IS_APPROX(w[i], v[i]);
}
v.resize(21);
v.resize(21,ti);
v[20] = x;
VERIFY_IS_APPROX(v[20], x);
v.resize(22,y);
@ -98,8 +98,8 @@ template<typename QuaternionType>
void check_stddeque_quaternion(const QuaternionType&)
{
typedef typename QuaternionType::Coefficients Coefficients;
QuaternionType x(Coefficients::Random()), y(Coefficients::Random());
std::deque<QuaternionType> v(10), w(20, y);
QuaternionType x(Coefficients::Random()), y(Coefficients::Random()), qi=QuaternionType::Identity();
std::deque<QuaternionType> v(10,qi), w(20, y);
v[5] = x;
w[6] = v[5];
VERIFY_IS_APPROX(w[6], v[5]);
@ -109,7 +109,7 @@ void check_stddeque_quaternion(const QuaternionType&)
VERIFY_IS_APPROX(w[i], v[i]);
}
v.resize(21);
v.resize(21,qi);
v[20] = x;
VERIFY_IS_APPROX(v[20], x);
v.resize(22,y);

16
test/stdlist.cpp
View File

@ -18,7 +18,7 @@ void check_stdlist_matrix(const MatrixType& m)
Index rows = m.rows();
Index cols = m.cols();
MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols);
std::list<MatrixType,Eigen::aligned_allocator<MatrixType> > v(10, MatrixType(rows,cols)), w(20, y);
std::list<MatrixType,Eigen::aligned_allocator<MatrixType> > v(10, MatrixType::Zero(rows,cols)), w(20, y);
v.front() = x;
w.front() = w.back();
VERIFY_IS_APPROX(w.front(), w.back());
@ -33,7 +33,7 @@ void check_stdlist_matrix(const MatrixType& m)
++wi;
}
v.resize(21);
v.resize(21, MatrixType::Zero(rows,cols));
v.back() = x;
VERIFY_IS_APPROX(v.back(), x);
v.resize(22,y);
@ -46,8 +46,8 @@ template<typename TransformType>
void check_stdlist_transform(const TransformType&)
{
typedef typename TransformType::MatrixType MatrixType;
TransformType x(MatrixType::Random()), y(MatrixType::Random());
std::list<TransformType,Eigen::aligned_allocator<TransformType> > v(10), w(20, y);
TransformType x(MatrixType::Random()), y(MatrixType::Random()), ti=TransformType::Identity();
std::list<TransformType,Eigen::aligned_allocator<TransformType> > v(10,ti), w(20, y);
v.front() = x;
w.front() = w.back();
VERIFY_IS_APPROX(w.front(), w.back());
@ -62,7 +62,7 @@ void check_stdlist_transform(const TransformType&)
++wi;
}
v.resize(21);
v.resize(21, ti);
v.back() = x;
VERIFY_IS_APPROX(v.back(), x);
v.resize(22,y);
@ -75,8 +75,8 @@ template<typename QuaternionType>
void check_stdlist_quaternion(const QuaternionType&)
{
typedef typename QuaternionType::Coefficients Coefficients;
QuaternionType x(Coefficients::Random()), y(Coefficients::Random());
std::list<QuaternionType,Eigen::aligned_allocator<QuaternionType> > v(10), w(20, y);
QuaternionType x(Coefficients::Random()), y(Coefficients::Random()), qi=QuaternionType::Identity();
std::list<QuaternionType,Eigen::aligned_allocator<QuaternionType> > v(10,qi), w(20, y);
v.front() = x;
w.front() = w.back();
VERIFY_IS_APPROX(w.front(), w.back());
@ -91,7 +91,7 @@ void check_stdlist_quaternion(const QuaternionType&)
++wi;
}
v.resize(21);
v.resize(21,qi);
v.back() = x;
VERIFY_IS_APPROX(v.back(), x);
v.resize(22,y);

14
test/stdlist_overload.cpp
View File

@ -47,7 +47,7 @@ void check_stdlist_matrix(const MatrixType& m)
Index rows = m.rows();
Index cols = m.cols();
MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols);
std::list<MatrixType> v(10, MatrixType(rows,cols)), w(20, y);
std::list<MatrixType> v(10, MatrixType::Zero(rows,cols)), w(20, y);
typename std::list<MatrixType>::iterator itv = get(v, 5);
typename std::list<MatrixType>::iterator itw = get(w, 6);
*itv = x;
@ -86,8 +86,8 @@ template<typename TransformType>
void check_stdlist_transform(const TransformType&)
{
typedef typename TransformType::MatrixType MatrixType;
TransformType x(MatrixType::Random()), y(MatrixType::Random());
std::list<TransformType> v(10), w(20, y);
TransformType x(MatrixType::Random()), y(MatrixType::Random()), ti=TransformType::Identity();
std::list<TransformType> v(10,ti), w(20, y);
typename std::list<TransformType>::iterator itv = get(v, 5);
typename std::list<TransformType>::iterator itw = get(w, 6);
*itv = x;
@ -103,7 +103,7 @@ void check_stdlist_transform(const TransformType&)
++itw;
}
v.resize(21);
v.resize(21, ti);
set(v, 20, x);
VERIFY_IS_APPROX(*get(v, 20), x);
v.resize(22,y);
@ -126,8 +126,8 @@ template<typename QuaternionType>
void check_stdlist_quaternion(const QuaternionType&)
{
typedef typename QuaternionType::Coefficients Coefficients;
QuaternionType x(Coefficients::Random()), y(Coefficients::Random());
std::list<QuaternionType> v(10), w(20, y);
QuaternionType x(Coefficients::Random()), y(Coefficients::Random()), qi=QuaternionType::Identity();
std::list<QuaternionType> v(10,qi), w(20, y);
typename std::list<QuaternionType>::iterator itv = get(v, 5);
typename std::list<QuaternionType>::iterator itw = get(w, 6);
*itv = x;
@ -143,7 +143,7 @@ void check_stdlist_quaternion(const QuaternionType&)
++itw;
}
v.resize(21);
v.resize(21,qi);
set(v, 20, x);
VERIFY_IS_APPROX(*get(v, 20), x);
v.resize(22,y);

16
test/stdvector.cpp
View File

@ -17,7 +17,7 @@ void check_stdvector_matrix(const MatrixType& m)
Index rows = m.rows();
Index cols = m.cols();
MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols);
std::vector<MatrixType,Eigen::aligned_allocator<MatrixType> > v(10, MatrixType(rows,cols)), w(20, y);
std::vector<MatrixType,Eigen::aligned_allocator<MatrixType> > v(10, MatrixType::Zero(rows,cols)), w(20, y);
v[5] = x;
w[6] = v[5];
VERIFY_IS_APPROX(w[6], v[5]);
@ -86,8 +86,8 @@ template<typename QuaternionType>
void check_stdvector_quaternion(const QuaternionType&)
{
typedef typename QuaternionType::Coefficients Coefficients;
QuaternionType x(Coefficients::Random()), y(Coefficients::Random());
std::vector<QuaternionType,Eigen::aligned_allocator<QuaternionType> > v(10), w(20, y);
QuaternionType x(Coefficients::Random()), y(Coefficients::Random()), qi=QuaternionType::Identity();
std::vector<QuaternionType,Eigen::aligned_allocator<QuaternionType> > v(10,qi), w(20, y);
v[5] = x;
w[6] = v[5];
VERIFY_IS_APPROX(w[6], v[5]);
@ -117,6 +117,16 @@ void check_stdvector_quaternion(const QuaternionType&)
}
}
// the code below triggered an invalid warning with gcc >= 7
// eigen/Eigen/src/Core/util/Memory.h:189:12: warning: argument 1 value '18446744073709551612' exceeds maximum object size 9223372036854775807
// This has been reported to gcc there: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=87544
void std_vector_gcc_warning()
{
typedef Eigen::Vector3f T;
std::vector<T, Eigen::aligned_allocator<T> > v;
v.push_back(T());
}
EIGEN_DECLARE_TEST(stdvector)
{
// some non vectorizable fixed sizes

6
test/stdvector_overload.cpp
View File

@ -31,7 +31,7 @@ void check_stdvector_matrix(const MatrixType& m)
Index rows = m.rows();
Index cols = m.cols();
MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols);
std::vector<MatrixType> v(10, MatrixType(rows,cols)), w(20, y);
std::vector<MatrixType> v(10, MatrixType::Zero(rows,cols)), w(20, y);
v[5] = x;
w[6] = v[5];
VERIFY_IS_APPROX(w[6], v[5]);
@ -100,8 +100,8 @@ template<typename QuaternionType>
void check_stdvector_quaternion(const QuaternionType&)
{
typedef typename QuaternionType::Coefficients Coefficients;
QuaternionType x(Coefficients::Random()), y(Coefficients::Random());
std::vector<QuaternionType> v(10), w(20, y);
QuaternionType x(Coefficients::Random()), y(Coefficients::Random()), qi=QuaternionType::Identity();
std::vector<QuaternionType> v(10,qi), w(20, y);
v[5] = x;
w[6] = v[5];
VERIFY_IS_APPROX(w[6], v[5]);

12
unsupported/Eigen/CXX11/src/Tensor/TensorContraction.h
View File

@ -186,21 +186,21 @@ struct TensorContractionKernel {
/*ConjugateLhs*/ false, /*ConjugateRhs*/ false>
GebpKernel;
EIGEN_DONT_INLINE
EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE
static void packLhs(LhsScalar* lhsBlock,
const typename LhsMapper::SubMapper& data_mapper,
const StorageIndex depth, const StorageIndex rows) {
LhsPacker()(lhsBlock, data_mapper, depth, rows, /*stride*/ 0, /*offset*/ 0);
}
EIGEN_DONT_INLINE
EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE
static void packRhs(RhsScalar* rhsBlock,
const typename RhsMapper::SubMapper& data_mapper,
const StorageIndex depth, const StorageIndex cols) {
RhsPacker()(rhsBlock, data_mapper, depth, cols);
}
EIGEN_DONT_INLINE
EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE
static void invoke(const OutputMapper& output_mapper,
const LhsScalar* lhsBlock, const RhsScalar* rhsBlock,
const StorageIndex rows, const StorageIndex depth,
@ -667,8 +667,8 @@ struct TensorContractionEvaluatorBase
this->m_device.memset(buffer, 0, m * n * sizeof(Scalar));
this->template evalGemmPartial<lhs_inner_dim_contiguous,
rhs_inner_dim_contiguous,
rhs_inner_dim_reordered, Alignment>(buffer,
0, k, 1);
rhs_inner_dim_reordered,
Alignment, true>(buffer, 0, k, 1);
}
template <bool lhs_inner_dim_contiguous, bool rhs_inner_dim_contiguous,
@ -681,7 +681,7 @@ struct TensorContractionEvaluatorBase
num_threads);
}
template <bool lhs_inner_dim_contiguous, bool rhs_inner_dim_contiguous, bool rhs_inner_dim_reordered, int Alignment, bool use_output_kernel = true>
template <bool lhs_inner_dim_contiguous, bool rhs_inner_dim_contiguous, bool rhs_inner_dim_reordered, int Alignment, bool use_output_kernel>
EIGEN_DEVICE_FUNC void evalGemmPartial(Scalar* buffer, Index k_start, Index k_end, int num_threads) const {
eigen_assert(k_end >= k_start && k_start >= 0 && k_end <= this->m_k_size);
// columns in slice on left side, rows on right side

2
unsupported/Eigen/CXX11/src/Tensor/TensorContractionThreadPool.h
View File

@ -794,7 +794,7 @@ struct TensorEvaluator<const TensorContractionOp<Indices, LeftArgType, RightArgT
Index num_blocks = divup<Index>(k, block_size);
// we use 'result' for the first block's partial result.
MaxSizeVector<Scalar*> block_buffers(num_blocks - 1);
Barrier barrier(num_blocks);
Barrier barrier(internal::convert_index<int>(num_blocks));
auto process_block = [=, &barrier](Scalar* buf, Index begin, Index end) {
::memset(buf, 0, m * n * sizeof(Scalar));
TENSOR_CONTRACTION_DISPATCH(

14
unsupported/Eigen/CXX11/src/Tensor/TensorMorphing.h
View File

@ -195,6 +195,14 @@ struct TensorEvaluator<const TensorReshapingOp<NewDimensions, ArgType>, Device>
m_impl.getResourceRequirements(resources);
}
// required in block(OutputTensorBlock* output_block) const
// For C++03 compatibility this must be defined outside the method
struct BlockIteratorState {
Index stride;
Index span;
Index size;
Index count;
};
// TODO(andydavis) Reduce the overhead of this function.
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void block(
OutputTensorBlock* output_block) const {
@ -219,12 +227,6 @@ struct TensorEvaluator<const TensorReshapingOp<NewDimensions, ArgType>, Device>
}
// Initialize output block iterator state.
struct BlockIteratorState {
Index stride;
Index span;
Index size;
Index count;
};
array<BlockIteratorState, NumOutputDims> block_iter_state;
for (Index i = 0; i < NumOutputDims; ++i) {

2
unsupported/Eigen/CXX11/src/Tensor/TensorReduction.h
View File

@ -218,6 +218,7 @@ struct InnerMostDimReducer<Self, Op, false, true> {
}
};
#if !defined(EIGEN_HIPCC)
template <typename Self, typename Op>
struct InnerMostDimReducer<Self, Op, true, true> {
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE typename Self::CoeffReturnType
@ -257,6 +258,7 @@ struct InnerMostDimReducer<Self, Op, true, true> {
}
}
};
#endif
template <int DimIndex, typename Self, typename Op, bool vectorizable = (Self::InputPacketAccess && Self::ReducerTraits::PacketAccess)>
struct InnerMostDimPreserver {

2
unsupported/Eigen/CXX11/src/Tensor/TensorReductionGpu.h
View File

@ -292,7 +292,7 @@ __global__ void FullReductionKernelHalfFloat(Reducer reducer, const Self input,
}
template <typename Op>
__global__ void ReductionCleanupKernelHalfFloat(Op& reducer, half* output, half2* scratch) {
__global__ void ReductionCleanupKernelHalfFloat(Op reducer, half* output, half2* scratch) {
eigen_assert(threadIdx.x == 1);
half tmp = __low2half(*scratch);
reducer.reduce(__high2half(*scratch), &tmp);

6
unsupported/Eigen/CXX11/src/Tensor/TensorScan.h
View File

@ -124,7 +124,11 @@ struct TensorEvaluator<const TensorScanOp<Op, ArgType>, Device> {
m_stride = m_stride * dims[i];
}
} else {
for (int i = NumDims - 1; i > op.axis(); --i) {
// dims can only be indexed through unsigned integers,
// so let's use an unsigned type to let the compiler knows.
// This prevents stupid warnings: ""'*((void*)(& evaluator)+64)[18446744073709551615]' may be used uninitialized in this function"
unsigned int axis = internal::convert_index<unsigned int>(op.axis());
for (unsigned int i = NumDims - 1; i > axis; --i) {
m_stride = m_stride * dims[i];
}
}

10
unsupported/test/cxx11_tensor_reduction.cpp
View File

@ -225,11 +225,11 @@ static void test_simple_reductions() {
Tensor<int, 1> ints(10);
std::iota(ints.data(), ints.data() + ints.dimension(0), 0);
TensorFixedSize<bool, Sizes<> > all;
all = ints.all();
VERIFY(!all());
all = (ints >= ints.constant(0)).all();
VERIFY(all());
TensorFixedSize<bool, Sizes<> > all_;
all_ = ints.all();
VERIFY(!all_());
all_ = (ints >= ints.constant(0)).all();
VERIFY(all_());
TensorFixedSize<bool, Sizes<> > any;
any = (ints > ints.constant(10)).any();

214
unsupported/test/mpreal/mpreal.h
View File

@ -5,7 +5,7 @@
Project homepage: http://www.holoborodko.com/pavel/mpfr
Contact e-mail: pavel@holoborodko.com
Copyright (c) 2008-2015 Pavel Holoborodko
Copyright (c) 2008-2016 Pavel Holoborodko
Contributors:
Dmitriy Gubanov, Konstantin Holoborodko, Brian Gladman,
@ -13,7 +13,7 @@
Pere Constans, Peter van Hoof, Gael Guennebaud, Tsai Chia Cheng,
Alexei Zubanov, Jauhien Piatlicki, Victor Berger, John Westwood,
Petr Aleksandrov, Orion Poplawski, Charles Karney, Arash Partow,
Rodney James, Jorge Leitao.
Rodney James, Jorge Leitao, Jerome Benoit.
Licensing:
(A) MPFR C++ is under GNU General Public License ("GPL").
@ -58,6 +58,7 @@
#include <limits>
#include <complex>
#include <algorithm>
#include <stdint.h>
// Options
#define MPREAL_HAVE_MSVC_DEBUGVIEW // Enable Debugger Visualizer for "Debug" builds in MSVC.
@ -68,16 +69,18 @@
// Library version
#define MPREAL_VERSION_MAJOR 3
#define MPREAL_VERSION_MINOR 6
#define MPREAL_VERSION_PATCHLEVEL 2
#define MPREAL_VERSION_STRING "3.6.2"
#define MPREAL_VERSION_PATCHLEVEL 5
#define MPREAL_VERSION_STRING "3.6.5"
// Detect compiler using signatures from http://predef.sourceforge.net/
#if defined(__GNUC__)
#define IsInf(x) (isinf)(x) // GNU C++/Intel ICC compiler on Linux
#if defined(__GNUC__) && defined(__INTEL_COMPILER)
#define IsInf(x) isinf EIGEN_NOT_A_MACRO (x) // Intel ICC compiler on Linux
#elif defined(_MSC_VER) // Microsoft Visual C++
#define IsInf(x) (!_finite(x))
#else
#define IsInf(x) (std::isinf)(x) // GNU C/C++ (and/or other compilers), just hope for C99 conformance
#define IsInf(x) std::isinf EIGEN_NOT_A_MACRO (x) // GNU C/C++ (and/or other compilers), just hope for C99 conformance
#endif
// A Clang feature extension to determine compiler features.
@ -85,10 +88,13 @@
#define __has_feature(x) 0
#endif
// Detect support for r-value references (move semantic). Borrowed from Eigen.
// Detect support for r-value references (move semantic).
// Move semantic should be enabled with great care in multi-threading environments,
// especially if MPFR uses custom memory allocators.
// Everything should be thread-safe and support passing ownership over thread boundary.
#if (__has_feature(cxx_rvalue_references) || \
defined(__GXX_EXPERIMENTAL_CXX0X__) || __cplusplus >= 201103L || \
(defined(_MSC_VER) && _MSC_VER >= 1600))
(defined(_MSC_VER) && _MSC_VER >= 1600) && !defined(MPREAL_DISABLE_MOVE_SEMANTIC))
#define MPREAL_HAVE_MOVE_SUPPORT
@ -99,8 +105,9 @@
// Detect support for explicit converters.
#if (__has_feature(cxx_explicit_conversions) || \
(defined(__GXX_EXPERIMENTAL_CXX0X__) && __GNUC_MINOR__ >= 5) || __cplusplus >= 201103L || \
(defined(_MSC_VER) && _MSC_VER >= 1800))
(defined(__GXX_EXPERIMENTAL_CXX0X__) && __GNUC_MINOR >= 5) || __cplusplus >= 201103L || \
(defined(_MSC_VER) && _MSC_VER >= 1800) || \
(defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1300))
#define MPREAL_HAVE_EXPLICIT_CONVERTERS
#endif
@ -127,7 +134,7 @@
// = -1 disables overflow checks (default)
// Fast replacement for mpfr_set_zero(x, +1):
// (a) uses low-level data members, might not be compatible with new versions of MPFR
// (a) uses low-level data members, might not be forward compatible
// (b) sign is not set, add (x)->_mpfr_sign = 1;
#define mpfr_set_zero_fast(x) ((x)->_mpfr_exp = __MPFR_EXP_ZERO)
@ -147,7 +154,7 @@ public:
// Get default rounding mode & precision
inline static mp_rnd_t get_default_rnd() { return (mp_rnd_t)(mpfr_get_default_rounding_mode()); }
inline static mp_prec_t get_default_prec() { return mpfr_get_default_prec(); }
inline static mp_prec_t get_default_prec() { return (mpfr_get_default_prec)(); }
// Constructors && type conversions
mpreal();
@ -354,6 +361,8 @@ public:
friend const mpreal log1p(const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal expm1(const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal nextpow2(const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode);
friend const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode);
@ -405,7 +414,7 @@ public:
friend const mpreal fms (const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode);
friend const mpreal agm (const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode);
friend const mpreal sum (const mpreal tab[], const unsigned long int n, int& status, mp_rnd_t rnd_mode);
friend int sgn(const mpreal& v); // returns -1 or +1
friend int sgn (const mpreal& v);
// MPFR 2.4.0 Specifics
#if (MPFR_VERSION >= MPFR_VERSION_NUM(2,4,0))
@ -482,7 +491,7 @@ public:
#endif
// Instance Checkers
friend bool (isnan) (const mpreal& v);
friend bool isnan EIGEN_NOT_A_MACRO (const mpreal& v);
friend bool (isinf) (const mpreal& v);
friend bool (isfinite) (const mpreal& v);
@ -509,7 +518,7 @@ public:
mpreal& setSign (int Sign, mp_rnd_t RoundingMode = get_default_rnd());
//Exponent
mp_exp_t get_exp();
mp_exp_t get_exp() const;
int set_exp(mp_exp_t e);
int check_range (int t, mp_rnd_t rnd_mode = get_default_rnd());
int subnormalize (int t, mp_rnd_t rnd_mode = get_default_rnd());
@ -580,7 +589,7 @@ inline mpreal::mpreal(const mpreal& u)
#ifdef MPREAL_HAVE_MOVE_SUPPORT
inline mpreal::mpreal(mpreal&& other)
{
mpfr_set_uninitialized(mpfr_ptr()); // make sure "other" holds no pointer to actual data
mpfr_set_uninitialized(mpfr_ptr()); // make sure "other" holds null-pointer (in uninitialized state)
mpfr_swap(mpfr_ptr(), other.mpfr_ptr());
MPREAL_MSVC_DEBUGVIEW_CODE;
@ -588,9 +597,11 @@ inline mpreal::mpreal(mpreal&& other)
inline mpreal& mpreal::operator=(mpreal&& other)
{
mpfr_swap(mpfr_ptr(), other.mpfr_ptr());
MPREAL_MSVC_DEBUGVIEW_CODE;
if (this != &other)
{
mpfr_swap(mpfr_ptr(), other.mpfr_ptr()); // destructor for "other" will be called just afterwards
MPREAL_MSVC_DEBUGVIEW_CODE;
}
return *this;
}
#endif
@ -639,13 +650,13 @@ inline mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode)
mpfr_init2(mpfr_ptr(), prec);
#if (MPREAL_DOUBLE_BITS_OVERFLOW > -1)
if(fits_in_bits(u, MPREAL_DOUBLE_BITS_OVERFLOW))
{
mpfr_set_d(mpfr_ptr(), u, mode);
}else
throw conversion_overflow();
if(fits_in_bits(u, MPREAL_DOUBLE_BITS_OVERFLOW))
{
mpfr_set_d(mpfr_ptr(), u, mode);
}else
throw conversion_overflow();
#else
mpfr_set_d(mpfr_ptr(), u, mode);
mpfr_set_d(mpfr_ptr(), u, mode);
#endif
MPREAL_MSVC_DEBUGVIEW_CODE;
@ -908,13 +919,13 @@ inline mpreal& mpreal::operator=(const mpreal& v)
{
if (this != &v)
{
mp_prec_t tp = mpfr_get_prec( mpfr_srcptr());
mp_prec_t vp = mpfr_get_prec(v.mpfr_srcptr());
mp_prec_t tp = mpfr_get_prec( mpfr_srcptr());
mp_prec_t vp = mpfr_get_prec(v.mpfr_srcptr());
if(tp != vp){
clear(mpfr_ptr());
mpfr_init2(mpfr_ptr(), vp);
}
if(tp != vp){
clear(mpfr_ptr());
mpfr_init2(mpfr_ptr(), vp);
}
mpfr_set(mpfr_ptr(), v.mpfr_srcptr(), mpreal::get_default_rnd());
@ -958,16 +969,16 @@ inline mpreal& mpreal::operator=(const long double v)
inline mpreal& mpreal::operator=(const double v)
{
#if (MPREAL_DOUBLE_BITS_OVERFLOW > -1)
if(fits_in_bits(v, MPREAL_DOUBLE_BITS_OVERFLOW))
{
mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd());
}else
throw conversion_overflow();
if(fits_in_bits(v, MPREAL_DOUBLE_BITS_OVERFLOW))
{
mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd());
}else
throw conversion_overflow();
#else
mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd());
mpfr_set_d(mpfr_ptr(),v,mpreal::get_default_rnd());
#endif
MPREAL_MSVC_DEBUGVIEW_CODE;
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@ -1161,9 +1172,9 @@ inline const mpreal mpreal::operator+()const { return mpreal(*this); }
inline const mpreal operator+(const mpreal& a, const mpreal& b)
{
mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
mpfr_add(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
return c;
mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
mpfr_add(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
return c;
}
inline mpreal& mpreal::operator++()
@ -1269,9 +1280,9 @@ inline const mpreal mpreal::operator-()const
inline const mpreal operator-(const mpreal& a, const mpreal& b)
{
mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
mpfr_sub(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
return c;
mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
mpfr_sub(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
return c;
}
inline const mpreal operator-(const double b, const mpreal& a)
@ -1386,9 +1397,9 @@ inline mpreal& mpreal::operator*=(const int v)
inline const mpreal operator*(const mpreal& a, const mpreal& b)
{
mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
mpfr_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
return c;
mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr())));
mpfr_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
return c;
}
//////////////////////////////////////////////////////////////////////////
@ -1462,9 +1473,9 @@ inline mpreal& mpreal::operator/=(const int v)
inline const mpreal operator/(const mpreal& a, const mpreal& b)
{
mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_srcptr()), mpfr_get_prec(b.mpfr_srcptr())));
mpfr_div(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
return c;
mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_srcptr()), mpfr_get_prec(b.mpfr_srcptr())));
mpfr_div(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd());
return c;
}
inline const mpreal operator/(const unsigned long int b, const mpreal& a)
@ -1659,7 +1670,7 @@ inline bool operator > (const mpreal& a, const double b ){ return !
inline bool operator >= (const mpreal& a, const mpreal& b ){ return (mpfr_greaterequal_p(a.mpfr_srcptr(),b.mpfr_srcptr()) != 0 ); }
inline bool operator >= (const mpreal& a, const unsigned long int b ){ return !isnan EIGEN_NOT_A_MACRO (a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) >= 0 ); }
// inline bool operator >= (const mpreal& a, const unsigned int b ){ return !isnan EIGEN_NOT_A_MACRO (isnan()a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) >= 0 ); }
inline bool operator >= (const mpreal& a, const unsigned int b ){ return !isnan EIGEN_NOT_A_MACRO (a) && (mpfr_cmp_ui(a.mpfr_srcptr(),b) >= 0 ); }
inline bool operator >= (const mpreal& a, const long int b ){ return !isnan EIGEN_NOT_A_MACRO (a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) >= 0 ); }
inline bool operator >= (const mpreal& a, const int b ){ return !isnan EIGEN_NOT_A_MACRO (a) && (mpfr_cmp_si(a.mpfr_srcptr(),b) >= 0 ); }
inline bool operator >= (const mpreal& a, const long double b ){ return !isnan EIGEN_NOT_A_MACRO (a) && (b == b) && (mpfr_cmp_ld(a.mpfr_srcptr(),b) >= 0 ); }
@ -1697,7 +1708,7 @@ inline bool operator != (const mpreal& a, const int b ){ return !
inline bool operator != (const mpreal& a, const long double b ){ return !(a == b); }
inline bool operator != (const mpreal& a, const double b ){ return !(a == b); }
inline bool (isnan) (const mpreal& op){ return (mpfr_nan_p (op.mpfr_srcptr()) != 0 ); }
inline bool isnan EIGEN_NOT_A_MACRO (const mpreal& op){ return (mpfr_nan_p (op.mpfr_srcptr()) != 0 ); }
inline bool (isinf) (const mpreal& op){ return (mpfr_inf_p (op.mpfr_srcptr()) != 0 ); }
inline bool (isfinite) (const mpreal& op){ return (mpfr_number_p (op.mpfr_srcptr()) != 0 ); }
inline bool iszero (const mpreal& op){ return (mpfr_zero_p (op.mpfr_srcptr()) != 0 ); }
@ -1762,7 +1773,7 @@ inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const
std::ostringstream format;
int digits = (n >= 0) ? n : 1 + bits2digits(mpfr_get_prec(mpfr_srcptr()));
int digits = (n >= 0) ? n : 2 + bits2digits(mpfr_get_prec(mpfr_srcptr()));
format << "%." << digits << "RNg";
@ -1931,14 +1942,9 @@ inline int bits2digits(mp_prec_t b)
//////////////////////////////////////////////////////////////////////////
// Set/Get number properties
inline int sgn(const mpreal& op)
{
return mpfr_sgn(op.mpfr_srcptr());
}
inline mpreal& mpreal::setSign(int sign, mp_rnd_t RoundingMode)
{
mpfr_setsign(mpfr_ptr(), mpfr_srcptr(), (sign < 0 ? 1 : 0), RoundingMode);
mpfr_setsign(mpfr_ptr(), mpfr_srcptr(), sign < 0, RoundingMode);
MPREAL_MSVC_DEBUGVIEW_CODE;
return *this;
}
@ -1993,7 +1999,7 @@ inline void mpreal::set_prec(mp_prec_t prec, mp_rnd_t rnd_mode)
MPREAL_MSVC_DEBUGVIEW_CODE;
}
inline mp_exp_t mpreal::get_exp ()
inline mp_exp_t mpreal::get_exp () const
{
return mpfr_get_exp(mpfr_srcptr());
}
@ -2091,6 +2097,12 @@ inline bool signbit(const mpreal& x)
return mpfr_signbit(x.mpfr_srcptr());
}
inline mpreal& setsignbit(mpreal& x, bool minus, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpfr_setsign(x.mpfr_ptr(), x.mpfr_srcptr(), minus, rnd_mode);
return x;
}
inline const mpreal modf(const mpreal& v, mpreal& n)
{
mpreal f(v);
@ -2194,7 +2206,7 @@ inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode)
inline const mpreal root(const mpreal& x, unsigned long int k, mp_rnd_t r = mpreal::get_default_rnd())
{
mpreal y(0, mpfr_get_prec(x.mpfr_srcptr()));
mpfr_root(y.mpfr_ptr(), x.mpfr_srcptr(), k, r);
mpfr_rootn_ui(y.mpfr_ptr(), x.mpfr_srcptr(), k, r);
return y;
}
@ -2270,6 +2282,16 @@ inline const mpreal besselj1(const mpreal& x, mp_rnd_t r = mpreal::get_default_r
inline const mpreal bessely0(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(y0 ); }
inline const mpreal bessely1(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { MPREAL_UNARY_MATH_FUNCTION_BODY(y1 ); }
inline const mpreal nextpow2(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd())
{
mpreal y(0, x.getPrecision());
if(!iszero(x))
y = ceil(log2(abs(x,r),r));
return y;
}
inline const mpreal atan2 (const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));
@ -2284,6 +2306,44 @@ inline const mpreal hypot (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode =
return a;
}
inline const mpreal hypot(const mpreal& a, const mpreal& b, const mpreal& c)
{
if(isnan EIGEN_NOT_A_MACRO (a) || isnan EIGEN_NOT_A_MACRO (b) || isnan EIGEN_NOT_A_MACRO(c)) return mpreal().setNan();
else
{
mpreal absa = abs(a), absb = abs(b), absc = abs(c);
mpreal w = (std::max)(absa, (std::max)(absb, absc));
mpreal r;
if (!iszero(w))
{
mpreal iw = 1/w;
r = w * sqrt(sqr(absa*iw) + sqr(absb*iw) + sqr(absc*iw));
}
return r;
}
}
inline const mpreal hypot(const mpreal& a, const mpreal& b, const mpreal& c, const mpreal& d)
{
if(isnan EIGEN_NOT_A_MACRO (a) || isnan EIGEN_NOT_A_MACRO (b) || isnan EIGEN_NOT_A_MACRO (c) || isnan EIGEN_NOT_A_MACRO (d)) return mpreal().setNan();
else
{
mpreal absa = abs(a), absb = abs(b), absc = abs(c), absd = abs(d);
mpreal w = (std::max)(absa, (std::max)(absb, (std::max)(absc, absd)));
mpreal r;
if (!iszero(w))
{
mpreal iw = 1/w;
r = w * sqrt(sqr(absa*iw) + sqr(absb*iw) + sqr(absc*iw) + sqr(absd*iw));
}
return r;
}
}
inline const mpreal remainder (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal a(0,(std::max)(y.getPrecision(), x.getPrecision()));
@ -2299,7 +2359,7 @@ inline const mpreal remquo (long* q, const mpreal& x, const mpreal& y, mp_rnd_t
}
inline const mpreal fac_ui (unsigned long int v, mp_prec_t prec = mpreal::get_default_prec(),
mp_rnd_t rnd_mode = mpreal::get_default_rnd())
mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x(0, prec);
mpfr_fac_ui(x.mpfr_ptr(),v,rnd_mode);
@ -2432,9 +2492,7 @@ inline const mpreal mod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = m
mpreal m = x - floor(x / y) * y;
m.setSign(sgn(y)); // make sure result has the same sign as Y
return m;
return copysign(abs(m),y); // make sure result has the same sign as Y
}
inline const mpreal fmod (const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
@ -2543,7 +2601,16 @@ inline const mpreal frac (const mpreal& x, mp_rnd_t r = mpreal::get_defaul
//////////////////////////////////////////////////////////////////////////
// Miscellaneous Functions
inline void swap (mpreal& a, mpreal& b) { mpfr_swap(a.mp,b.mp); }
inline int sgn(const mpreal& op)
{
// Please note, this is classic signum function which ignores sign of zero.
// Use signbit if you need sign of zero.
return mpfr_sgn(op.mpfr_srcptr());
}
//////////////////////////////////////////////////////////////////////////
// Miscellaneous Functions
inline void swap (mpreal& a, mpreal& b) { mpfr_swap(a.mpfr_ptr(),b.mpfr_ptr()); }
inline const mpreal (max)(const mpreal& x, const mpreal& y){ return (x>y?x:y); }
inline const mpreal (min)(const mpreal& x, const mpreal& y){ return (x<y?x:y); }
@ -2634,14 +2701,23 @@ inline const mpreal random(unsigned int seed = 0)
}
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,1,0))
#if (MPFR_VERSION >= MPFR_VERSION_NUM(3,1,0) && MPFR_VERSION < MPFR_VERSION_NUM(4,0,0))
// TODO:
// Use mpfr_nrandom since mpfr_grandom is deprecated
#if defined(_MSC_VER)
#pragma warning( push )
#pragma warning( disable : 1478)
#endif
inline const mpreal grandom (gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd())
{
mpreal x;
mpfr_grandom(x.mpfr_ptr(), NULL, state, rnd_mode);
return x;
}
#if defined(_MSC_VER)
#pragma warning( pop )
#endif
inline const mpreal grandom(unsigned int seed = 0)
{
@ -2981,7 +3057,7 @@ inline const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode)
// Non-throwing swap C++ idiom: http://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Non-throwing_swap
namespace std
{
// we are allowed to extend namespace std with specializations only
// we are allowed to extend namespace std with specializations only
template <>
inline void swap(mpfr::mpreal& x, mpfr::mpreal& y)
{