special_packetmath.cpp
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #include <limits>
12 #include "packetmath_test_shared.h"
13 #include "../Eigen/SpecialFunctions"
14 
15 template<typename Scalar,typename Packet> void packetmath_real()
16 {
17  using std::abs;
18  typedef internal::packet_traits<Scalar> PacketTraits;
19  const int PacketSize = internal::unpacket_traits<Packet>::size;
20 
21  const int size = PacketSize*4;
22  EIGEN_ALIGN_MAX Scalar data1[PacketSize*4];
23  EIGEN_ALIGN_MAX Scalar data2[PacketSize*4];
24  EIGEN_ALIGN_MAX Scalar ref[PacketSize*4];
25 
26 #if EIGEN_HAS_C99_MATH
27  {
28  data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
29  test::packet_helper<internal::packet_traits<Scalar>::HasLGamma,Packet> h;
30  h.store(data2, internal::plgamma(h.load(data1)));
31  VERIFY((numext::isnan)(data2[0]));
32  }
33  if (internal::packet_traits<Scalar>::HasErf) {
34  data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
35  test::packet_helper<internal::packet_traits<Scalar>::HasErf,Packet> h;
36  h.store(data2, internal::perf(h.load(data1)));
37  VERIFY((numext::isnan)(data2[0]));
38  }
39  {
40  data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
41  test::packet_helper<internal::packet_traits<Scalar>::HasErfc,Packet> h;
42  h.store(data2, internal::perfc(h.load(data1)));
43  VERIFY((numext::isnan)(data2[0]));
44  }
45  {
46  for (int i=0; i<size; ++i) {
47  data1[i] = internal::random<Scalar>(Scalar(0),Scalar(1));
48  }
49  CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasNdtri, numext::ndtri, internal::pndtri);
50  }
51 #endif // EIGEN_HAS_C99_MATH
52 
53  // For bessel_i*e and bessel_j*, the valid range is negative reals.
54  {
55  const int max_exponent = numext::mini(std::numeric_limits<Scalar>::max_exponent10-1, 6);
56  for (int i=0; i<size; ++i)
57  {
58  data1[i] = internal::random<Scalar>(Scalar(-1),Scalar(1)) * Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-max_exponent),Scalar(max_exponent))));
59  data2[i] = internal::random<Scalar>(Scalar(-1),Scalar(1)) * Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-max_exponent),Scalar(max_exponent))));
60  }
61 
62  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0e, internal::pbessel_i0e);
63  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1e, internal::pbessel_i1e);
64  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_j0, internal::pbessel_j0);
65  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_j1, internal::pbessel_j1);
66  }
67 
68  // Use a smaller data range for the bessel_i* as these can become very large.
69  // Following #1693, we also restrict this range further to avoid inf's due to
70  // differences in pexp and exp.
71  for (int i=0; i<size; ++i) {
72  data1[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) *
73  Scalar(std::pow(Scalar(9), internal::random<Scalar>(Scalar(-1),Scalar(2))));
74  data2[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) *
75  Scalar(std::pow(Scalar(9), internal::random<Scalar>(Scalar(-1),Scalar(2))));
76  }
77  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0, internal::pbessel_i0);
78  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1, internal::pbessel_i1);
79 
80 
81  // y_i, and k_i are valid for x > 0.
82  {
83  const int max_exponent = numext::mini(std::numeric_limits<Scalar>::max_exponent10-1, 5);
84  for (int i=0; i<size; ++i)
85  {
86  data1[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) * Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-2),Scalar(max_exponent))));
87  data2[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) * Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-2),Scalar(max_exponent))));
88  }
89  }
90 
91  // TODO(srvasude): Re-enable this test once properly investigated why the
92  // scalar and vector paths differ.
93  // CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_y0, internal::pbessel_y0);
94  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_y1, internal::pbessel_y1);
95  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k0e, internal::pbessel_k0e);
96  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k1e, internal::pbessel_k1e);
97 
98  // Following #1693, we restrict the range for exp to avoid zeroing out too
99  // fast.
100  for (int i=0; i<size; ++i) {
101  data1[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) *
102  Scalar(std::pow(Scalar(9), internal::random<Scalar>(Scalar(-1),Scalar(2))));
103  data2[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) *
104  Scalar(std::pow(Scalar(9), internal::random<Scalar>(Scalar(-1),Scalar(2))));
105  }
106  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k0, internal::pbessel_k0);
107  CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k1, internal::pbessel_k1);
108 
109 
110  for (int i=0; i<size; ++i) {
111  data1[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) *
112  Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-1),Scalar(2))));
113  data2[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) *
114  Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-1),Scalar(2))));
115  }
116 
117 #if EIGEN_HAS_C99_MATH && (EIGEN_COMP_CXXVER >= 11)
118  CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLGamma, std::lgamma, internal::plgamma);
119  CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErf, std::erf, internal::perf);
120  CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErfc, std::erfc, internal::perfc);
121 #endif
122 
123 }
124 
125 namespace Eigen {
126 namespace test {
127 
128 template<typename Scalar,typename PacketType, bool IsComplex, bool IsInteger>
129 struct runall {
130  static void run() {
131  packetmath_real<Scalar,PacketType>();
132  }
133 };
134 
135 }
136 }
137 
138 EIGEN_DECLARE_TEST(special_packetmath)
139 {
140  g_first_pass = true;
141  for(int i = 0; i < g_repeat; i++) {
142 
147  g_first_pass = false;
148  }
149 }
EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp< Eigen::internal::scalar_bessel_k1e_op< typename Derived::Scalar >, const Derived > bessel_k1e(const Eigen::ArrayBase< Derived > &x)
EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp< Eigen::internal::scalar_bessel_i1e_op< typename Derived::Scalar >, const Derived > bessel_i1e(const Eigen::ArrayBase< Derived > &x)
SCALAR Scalar
Definition: bench_gemm.cpp:46
EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pndtri(const Packet &a)
#define CALL_SUBTEST_4(FUNC)
EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_j1(const Packet &x)
EIGEN_DEVICE_FUNC const NdtriReturnType ndtri() const
EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet perf(const Packet &a)
EIGEN_DEVICE_FUNC const ErfReturnType erf() const
#define CALL_SUBTEST_3(FUNC)
EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp< Eigen::internal::scalar_bessel_j1_op< typename Derived::Scalar >, const Derived > bessel_j1(const Eigen::ArrayBase< Derived > &x)
Definition: test.py:1
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
EIGEN_DECLARE_TEST(special_packetmath)
#define EIGEN_ALIGN_MAX
EIGEN_DEVICE_FUNC const ErfcReturnType erfc() const
EIGEN_DEVICE_FUNC const LgammaReturnType lgamma() const
EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_i0e(const Packet &x)
bool g_first_pass
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
#define CHECK_CWISE1_IF(COND, REFOP, POP)
EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plgamma(const Packet &a)
#define CALL_SUBTEST_1(FUNC)
EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_i0(const Packet &x)
static int g_repeat
Definition: main.h:169
EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_y1(const Packet &x)
EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_k0(const Packet &x)
EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp< Eigen::internal::scalar_bessel_i0e_op< typename Derived::Scalar >, const Derived > bessel_i0e(const Eigen::ArrayBase< Derived > &x)
EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet perfc(const Packet &a)
Reference counting helper.
Definition: object.h:67
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T mini(const T &x, const T &y)
EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_k1e(const Packet &x)
const double h
EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp< Eigen::internal::scalar_bessel_i0_op< typename Derived::Scalar >, const Derived > bessel_i0(const Eigen::ArrayBase< Derived > &x)
EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp< Eigen::internal::scalar_bessel_k1_op< typename Derived::Scalar >, const Derived > bessel_k1(const Eigen::ArrayBase< Derived > &x)
EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp< Eigen::internal::scalar_bessel_k0e_op< typename Derived::Scalar >, const Derived > bessel_k0e(const Eigen::ArrayBase< Derived > &x)
#define VERIFY(a)
Definition: main.h:380
EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_i1e(const Packet &x)
EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_j0(const Packet &x)
EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_k1(const Packet &x)
void packetmath_real()
#define CALL_SUBTEST_2(FUNC)
EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp< Eigen::internal::scalar_bessel_y1_op< typename Derived::Scalar >, const Derived > bessel_y1(const Eigen::ArrayBase< Derived > &x)
Jet< T, N > pow(const Jet< T, N > &f, double g)
Definition: jet.h:570
EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_i1(const Packet &x)
#define abs(x)
Definition: datatypes.h:17
EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp< Eigen::internal::scalar_bessel_k0_op< typename Derived::Scalar >, const Derived > bessel_k0(const Eigen::ArrayBase< Derived > &x)
EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp< Eigen::internal::scalar_bessel_j0_op< typename Derived::Scalar >, const Derived > bessel_j0(const Eigen::ArrayBase< Derived > &x)
#define isnan(X)
Definition: main.h:93
EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp< Eigen::internal::scalar_bessel_i1_op< typename Derived::Scalar >, const Derived > bessel_i1(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet pbessel_k0e(const Packet &x)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:36:18