packetmath.cpp (55436B)
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 "packetmath_test_shared.h" 12 #include "random_without_cast_overflow.h" 13 14 template <typename T> 15 inline T REF_ADD(const T& a, const T& b) { 16 return a + b; 17 } 18 template <typename T> 19 inline T REF_SUB(const T& a, const T& b) { 20 return a - b; 21 } 22 template <typename T> 23 inline T REF_MUL(const T& a, const T& b) { 24 return a * b; 25 } 26 template <typename T> 27 inline T REF_DIV(const T& a, const T& b) { 28 return a / b; 29 } 30 template <typename T> 31 inline T REF_ABS_DIFF(const T& a, const T& b) { 32 return a > b ? a - b : b - a; 33 } 34 35 // Specializations for bool. 36 template <> 37 inline bool REF_ADD(const bool& a, const bool& b) { 38 return a || b; 39 } 40 template <> 41 inline bool REF_SUB(const bool& a, const bool& b) { 42 return a ^ b; 43 } 44 template <> 45 inline bool REF_MUL(const bool& a, const bool& b) { 46 return a && b; 47 } 48 49 template <typename T> 50 inline T REF_FREXP(const T& x, T& exp) { 51 int iexp; 52 EIGEN_USING_STD(frexp) 53 const T out = static_cast<T>(frexp(x, &iexp)); 54 exp = static_cast<T>(iexp); 55 return out; 56 } 57 58 template <typename T> 59 inline T REF_LDEXP(const T& x, const T& exp) { 60 EIGEN_USING_STD(ldexp) 61 return static_cast<T>(ldexp(x, static_cast<int>(exp))); 62 } 63 64 // Uses pcast to cast from one array to another. 65 template <typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio> 66 struct pcast_array; 67 68 template <typename SrcPacket, typename TgtPacket, int TgtCoeffRatio> 69 struct pcast_array<SrcPacket, TgtPacket, 1, TgtCoeffRatio> { 70 typedef typename internal::unpacket_traits<SrcPacket>::type SrcScalar; 71 typedef typename internal::unpacket_traits<TgtPacket>::type TgtScalar; 72 static void cast(const SrcScalar* src, size_t size, TgtScalar* dst) { 73 static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size; 74 static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size; 75 size_t i; 76 for (i = 0; i < size && i + SrcPacketSize <= size; i += TgtPacketSize) { 77 internal::pstoreu(dst + i, internal::pcast<SrcPacket, TgtPacket>(internal::ploadu<SrcPacket>(src + i))); 78 } 79 // Leftovers that cannot be loaded into a packet. 80 for (; i < size; ++i) { 81 dst[i] = static_cast<TgtScalar>(src[i]); 82 } 83 } 84 }; 85 86 template <typename SrcPacket, typename TgtPacket> 87 struct pcast_array<SrcPacket, TgtPacket, 2, 1> { 88 static void cast(const typename internal::unpacket_traits<SrcPacket>::type* src, size_t size, 89 typename internal::unpacket_traits<TgtPacket>::type* dst) { 90 static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size; 91 static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size; 92 for (size_t i = 0; i < size; i += TgtPacketSize) { 93 SrcPacket a = internal::ploadu<SrcPacket>(src + i); 94 SrcPacket b = internal::ploadu<SrcPacket>(src + i + SrcPacketSize); 95 internal::pstoreu(dst + i, internal::pcast<SrcPacket, TgtPacket>(a, b)); 96 } 97 } 98 }; 99 100 template <typename SrcPacket, typename TgtPacket> 101 struct pcast_array<SrcPacket, TgtPacket, 4, 1> { 102 static void cast(const typename internal::unpacket_traits<SrcPacket>::type* src, size_t size, 103 typename internal::unpacket_traits<TgtPacket>::type* dst) { 104 static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size; 105 static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size; 106 for (size_t i = 0; i < size; i += TgtPacketSize) { 107 SrcPacket a = internal::ploadu<SrcPacket>(src + i); 108 SrcPacket b = internal::ploadu<SrcPacket>(src + i + SrcPacketSize); 109 SrcPacket c = internal::ploadu<SrcPacket>(src + i + 2 * SrcPacketSize); 110 SrcPacket d = internal::ploadu<SrcPacket>(src + i + 3 * SrcPacketSize); 111 internal::pstoreu(dst + i, internal::pcast<SrcPacket, TgtPacket>(a, b, c, d)); 112 } 113 } 114 }; 115 116 template <typename SrcPacket, typename TgtPacket> 117 struct pcast_array<SrcPacket, TgtPacket, 8, 1> { 118 static void cast(const typename internal::unpacket_traits<SrcPacket>::type* src, size_t size, 119 typename internal::unpacket_traits<TgtPacket>::type* dst) { 120 static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size; 121 static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size; 122 for (size_t i = 0; i < size; i += TgtPacketSize) { 123 SrcPacket a = internal::ploadu<SrcPacket>(src + i); 124 SrcPacket b = internal::ploadu<SrcPacket>(src + i + SrcPacketSize); 125 SrcPacket c = internal::ploadu<SrcPacket>(src + i + 2 * SrcPacketSize); 126 SrcPacket d = internal::ploadu<SrcPacket>(src + i + 3 * SrcPacketSize); 127 SrcPacket e = internal::ploadu<SrcPacket>(src + i + 4 * SrcPacketSize); 128 SrcPacket f = internal::ploadu<SrcPacket>(src + i + 5 * SrcPacketSize); 129 SrcPacket g = internal::ploadu<SrcPacket>(src + i + 6 * SrcPacketSize); 130 SrcPacket h = internal::ploadu<SrcPacket>(src + i + 7 * SrcPacketSize); 131 internal::pstoreu(dst + i, internal::pcast<SrcPacket, TgtPacket>(a, b, c, d, e, f, g, h)); 132 } 133 } 134 }; 135 136 template <typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio, bool CanCast = false> 137 struct test_cast_helper; 138 139 template <typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio> 140 struct test_cast_helper<SrcPacket, TgtPacket, SrcCoeffRatio, TgtCoeffRatio, false> { 141 static void run() {} 142 }; 143 144 template <typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio> 145 struct test_cast_helper<SrcPacket, TgtPacket, SrcCoeffRatio, TgtCoeffRatio, true> { 146 static void run() { 147 typedef typename internal::unpacket_traits<SrcPacket>::type SrcScalar; 148 typedef typename internal::unpacket_traits<TgtPacket>::type TgtScalar; 149 static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size; 150 static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size; 151 static const int BlockSize = SrcPacketSize * SrcCoeffRatio; 152 eigen_assert(BlockSize == TgtPacketSize * TgtCoeffRatio && "Packet sizes and cast ratios are mismatched."); 153 154 static const int DataSize = 10 * BlockSize; 155 EIGEN_ALIGN_MAX SrcScalar data1[DataSize]; 156 EIGEN_ALIGN_MAX TgtScalar data2[DataSize]; 157 EIGEN_ALIGN_MAX TgtScalar ref[DataSize]; 158 159 // Construct a packet of scalars that will not overflow when casting 160 for (int i = 0; i < DataSize; ++i) { 161 data1[i] = internal::random_without_cast_overflow<SrcScalar, TgtScalar>::value(); 162 } 163 164 for (int i = 0; i < DataSize; ++i) { 165 ref[i] = static_cast<const TgtScalar>(data1[i]); 166 } 167 168 pcast_array<SrcPacket, TgtPacket, SrcCoeffRatio, TgtCoeffRatio>::cast(data1, DataSize, data2); 169 170 VERIFY(test::areApprox(ref, data2, DataSize) && "internal::pcast<>"); 171 } 172 }; 173 174 template <typename SrcPacket, typename TgtPacket> 175 struct test_cast { 176 static void run() { 177 typedef typename internal::unpacket_traits<SrcPacket>::type SrcScalar; 178 typedef typename internal::unpacket_traits<TgtPacket>::type TgtScalar; 179 typedef typename internal::type_casting_traits<SrcScalar, TgtScalar> TypeCastingTraits; 180 static const int SrcCoeffRatio = TypeCastingTraits::SrcCoeffRatio; 181 static const int TgtCoeffRatio = TypeCastingTraits::TgtCoeffRatio; 182 static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size; 183 static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size; 184 static const bool HasCast = 185 internal::unpacket_traits<SrcPacket>::vectorizable && internal::unpacket_traits<TgtPacket>::vectorizable && 186 TypeCastingTraits::VectorizedCast && (SrcPacketSize * SrcCoeffRatio == TgtPacketSize * TgtCoeffRatio); 187 test_cast_helper<SrcPacket, TgtPacket, SrcCoeffRatio, TgtCoeffRatio, HasCast>::run(); 188 } 189 }; 190 191 template <typename SrcPacket, typename TgtScalar, 192 typename TgtPacket = typename internal::packet_traits<TgtScalar>::type, 193 bool Vectorized = internal::packet_traits<TgtScalar>::Vectorizable, 194 bool HasHalf = !internal::is_same<typename internal::unpacket_traits<TgtPacket>::half, TgtPacket>::value> 195 struct test_cast_runner; 196 197 template <typename SrcPacket, typename TgtScalar, typename TgtPacket> 198 struct test_cast_runner<SrcPacket, TgtScalar, TgtPacket, true, false> { 199 static void run() { test_cast<SrcPacket, TgtPacket>::run(); } 200 }; 201 202 template <typename SrcPacket, typename TgtScalar, typename TgtPacket> 203 struct test_cast_runner<SrcPacket, TgtScalar, TgtPacket, true, true> { 204 static void run() { 205 test_cast<SrcPacket, TgtPacket>::run(); 206 test_cast_runner<SrcPacket, TgtScalar, typename internal::unpacket_traits<TgtPacket>::half>::run(); 207 } 208 }; 209 210 template <typename SrcPacket, typename TgtScalar, typename TgtPacket> 211 struct test_cast_runner<SrcPacket, TgtScalar, TgtPacket, false, false> { 212 static void run() {} 213 }; 214 215 template <typename Scalar, typename Packet, typename EnableIf = void> 216 struct packetmath_pcast_ops_runner { 217 static void run() { 218 test_cast_runner<Packet, float>::run(); 219 test_cast_runner<Packet, double>::run(); 220 test_cast_runner<Packet, int8_t>::run(); 221 test_cast_runner<Packet, uint8_t>::run(); 222 test_cast_runner<Packet, int16_t>::run(); 223 test_cast_runner<Packet, uint16_t>::run(); 224 test_cast_runner<Packet, int32_t>::run(); 225 test_cast_runner<Packet, uint32_t>::run(); 226 test_cast_runner<Packet, int64_t>::run(); 227 test_cast_runner<Packet, uint64_t>::run(); 228 test_cast_runner<Packet, bool>::run(); 229 test_cast_runner<Packet, std::complex<float> >::run(); 230 test_cast_runner<Packet, std::complex<double> >::run(); 231 test_cast_runner<Packet, half>::run(); 232 test_cast_runner<Packet, bfloat16>::run(); 233 } 234 }; 235 236 // Only some types support cast from std::complex<>. 237 template <typename Scalar, typename Packet> 238 struct packetmath_pcast_ops_runner<Scalar, Packet, typename internal::enable_if<NumTraits<Scalar>::IsComplex>::type> { 239 static void run() { 240 test_cast_runner<Packet, std::complex<float> >::run(); 241 test_cast_runner<Packet, std::complex<double> >::run(); 242 test_cast_runner<Packet, half>::run(); 243 test_cast_runner<Packet, bfloat16>::run(); 244 } 245 }; 246 247 template <typename Scalar, typename Packet> 248 void packetmath_boolean_mask_ops() { 249 const int PacketSize = internal::unpacket_traits<Packet>::size; 250 const int size = 2 * PacketSize; 251 EIGEN_ALIGN_MAX Scalar data1[size]; 252 EIGEN_ALIGN_MAX Scalar data2[size]; 253 EIGEN_ALIGN_MAX Scalar ref[size]; 254 255 for (int i = 0; i < size; ++i) { 256 data1[i] = internal::random<Scalar>(); 257 } 258 CHECK_CWISE1(internal::ptrue, internal::ptrue); 259 CHECK_CWISE2_IF(true, internal::pandnot, internal::pandnot); 260 for (int i = 0; i < PacketSize; ++i) { 261 data1[i] = Scalar(i); 262 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0); 263 } 264 265 CHECK_CWISE2_IF(true, internal::pcmp_eq, internal::pcmp_eq); 266 267 //Test (-0) == (0) for signed operations 268 for (int i = 0; i < PacketSize; ++i) { 269 data1[i] = Scalar(-0.0); 270 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0); 271 } 272 CHECK_CWISE2_IF(true, internal::pcmp_eq, internal::pcmp_eq); 273 274 //Test NaN 275 for (int i = 0; i < PacketSize; ++i) { 276 data1[i] = NumTraits<Scalar>::quiet_NaN(); 277 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0); 278 } 279 CHECK_CWISE2_IF(true, internal::pcmp_eq, internal::pcmp_eq); 280 } 281 282 template <typename Scalar, typename Packet> 283 void packetmath_boolean_mask_ops_real() { 284 const int PacketSize = internal::unpacket_traits<Packet>::size; 285 const int size = 2 * PacketSize; 286 EIGEN_ALIGN_MAX Scalar data1[size]; 287 EIGEN_ALIGN_MAX Scalar data2[size]; 288 EIGEN_ALIGN_MAX Scalar ref[size]; 289 290 for (int i = 0; i < PacketSize; ++i) { 291 data1[i] = internal::random<Scalar>(); 292 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0); 293 } 294 295 CHECK_CWISE2_IF(true, internal::pcmp_lt_or_nan, internal::pcmp_lt_or_nan); 296 297 //Test (-0) <=/< (0) for signed operations 298 for (int i = 0; i < PacketSize; ++i) { 299 data1[i] = Scalar(-0.0); 300 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0); 301 } 302 CHECK_CWISE2_IF(true, internal::pcmp_lt_or_nan, internal::pcmp_lt_or_nan); 303 304 //Test NaN 305 for (int i = 0; i < PacketSize; ++i) { 306 data1[i] = NumTraits<Scalar>::quiet_NaN(); 307 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0); 308 } 309 CHECK_CWISE2_IF(true, internal::pcmp_lt_or_nan, internal::pcmp_lt_or_nan); 310 } 311 312 template <typename Scalar, typename Packet> 313 void packetmath_boolean_mask_ops_notcomplex() { 314 const int PacketSize = internal::unpacket_traits<Packet>::size; 315 const int size = 2 * PacketSize; 316 EIGEN_ALIGN_MAX Scalar data1[size]; 317 EIGEN_ALIGN_MAX Scalar data2[size]; 318 EIGEN_ALIGN_MAX Scalar ref[size]; 319 320 for (int i = 0; i < PacketSize; ++i) { 321 data1[i] = internal::random<Scalar>(); 322 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0); 323 } 324 325 CHECK_CWISE2_IF(true, internal::pcmp_le, internal::pcmp_le); 326 CHECK_CWISE2_IF(true, internal::pcmp_lt, internal::pcmp_lt); 327 328 //Test (-0) <=/< (0) for signed operations 329 for (int i = 0; i < PacketSize; ++i) { 330 data1[i] = Scalar(-0.0); 331 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0); 332 } 333 CHECK_CWISE2_IF(true, internal::pcmp_le, internal::pcmp_le); 334 CHECK_CWISE2_IF(true, internal::pcmp_lt, internal::pcmp_lt); 335 336 //Test NaN 337 for (int i = 0; i < PacketSize; ++i) { 338 data1[i] = NumTraits<Scalar>::quiet_NaN(); 339 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0); 340 } 341 CHECK_CWISE2_IF(true, internal::pcmp_le, internal::pcmp_le); 342 CHECK_CWISE2_IF(true, internal::pcmp_lt, internal::pcmp_lt); 343 } 344 345 // Packet16b representing bool does not support ptrue, pandnot or pcmp_eq, since the scalar path 346 // (for some compilers) compute the bitwise and with 0x1 of the results to keep the value in [0,1]. 347 template<> 348 void packetmath_boolean_mask_ops<bool, internal::packet_traits<bool>::type>() {} 349 template<> 350 void packetmath_boolean_mask_ops_notcomplex<bool, internal::packet_traits<bool>::type>() {} 351 352 template <typename Scalar, typename Packet> 353 void packetmath_minus_zero_add() { 354 const int PacketSize = internal::unpacket_traits<Packet>::size; 355 const int size = 2 * PacketSize; 356 EIGEN_ALIGN_MAX Scalar data1[size]; 357 EIGEN_ALIGN_MAX Scalar data2[size]; 358 EIGEN_ALIGN_MAX Scalar ref[size]; 359 360 for (int i = 0; i < PacketSize; ++i) { 361 data1[i] = Scalar(-0.0); 362 data1[i + PacketSize] = Scalar(-0.0); 363 } 364 CHECK_CWISE2_IF(internal::packet_traits<Scalar>::HasAdd, REF_ADD, internal::padd); 365 } 366 367 // Ensure optimization barrier compiles and doesn't modify contents. 368 // Only applies to raw types, so will not work for std::complex, Eigen::half 369 // or Eigen::bfloat16. For those you would need to refer to an underlying 370 // storage element. 371 template<typename Packet, typename EnableIf = void> 372 struct eigen_optimization_barrier_test { 373 static void run() {} 374 }; 375 376 template<typename Packet> 377 struct eigen_optimization_barrier_test<Packet, typename internal::enable_if< 378 !NumTraits<Packet>::IsComplex && 379 !internal::is_same<Packet, Eigen::half>::value && 380 !internal::is_same<Packet, Eigen::bfloat16>::value 381 >::type> { 382 static void run() { 383 typedef typename internal::unpacket_traits<Packet>::type Scalar; 384 Scalar s = internal::random<Scalar>(); 385 Packet barrier = internal::pset1<Packet>(s); 386 EIGEN_OPTIMIZATION_BARRIER(barrier); 387 eigen_assert(s == internal::pfirst(barrier) && "EIGEN_OPTIMIZATION_BARRIER"); 388 } 389 }; 390 391 template <typename Scalar, typename Packet> 392 void packetmath() { 393 typedef internal::packet_traits<Scalar> PacketTraits; 394 const int PacketSize = internal::unpacket_traits<Packet>::size; 395 typedef typename NumTraits<Scalar>::Real RealScalar; 396 397 if (g_first_pass) 398 std::cerr << "=== Testing packet of type '" << typeid(Packet).name() << "' and scalar type '" 399 << typeid(Scalar).name() << "' and size '" << PacketSize << "' ===\n"; 400 401 const int max_size = PacketSize > 4 ? PacketSize : 4; 402 const int size = PacketSize * max_size; 403 EIGEN_ALIGN_MAX Scalar data1[size]; 404 EIGEN_ALIGN_MAX Scalar data2[size]; 405 EIGEN_ALIGN_MAX Scalar data3[size]; 406 EIGEN_ALIGN_MAX Scalar ref[size]; 407 RealScalar refvalue = RealScalar(0); 408 409 eigen_optimization_barrier_test<Packet>::run(); 410 eigen_optimization_barrier_test<Scalar>::run(); 411 412 for (int i = 0; i < size; ++i) { 413 data1[i] = internal::random<Scalar>() / RealScalar(PacketSize); 414 data2[i] = internal::random<Scalar>() / RealScalar(PacketSize); 415 refvalue = (std::max)(refvalue, numext::abs(data1[i])); 416 } 417 418 internal::pstore(data2, internal::pload<Packet>(data1)); 419 VERIFY(test::areApprox(data1, data2, PacketSize) && "aligned load/store"); 420 421 for (int offset = 0; offset < PacketSize; ++offset) { 422 internal::pstore(data2, internal::ploadu<Packet>(data1 + offset)); 423 VERIFY(test::areApprox(data1 + offset, data2, PacketSize) && "internal::ploadu"); 424 } 425 426 for (int offset = 0; offset < PacketSize; ++offset) { 427 internal::pstoreu(data2 + offset, internal::pload<Packet>(data1)); 428 VERIFY(test::areApprox(data1, data2 + offset, PacketSize) && "internal::pstoreu"); 429 } 430 431 if (internal::unpacket_traits<Packet>::masked_load_available) { 432 test::packet_helper<internal::unpacket_traits<Packet>::masked_load_available, Packet> h; 433 unsigned long long max_umask = (0x1ull << PacketSize); 434 435 for (int offset = 0; offset < PacketSize; ++offset) { 436 for (unsigned long long umask = 0; umask < max_umask; ++umask) { 437 h.store(data2, h.load(data1 + offset, umask)); 438 for (int k = 0; k < PacketSize; ++k) data3[k] = ((umask & (0x1ull << k)) >> k) ? data1[k + offset] : Scalar(0); 439 VERIFY(test::areApprox(data3, data2, PacketSize) && "internal::ploadu masked"); 440 } 441 } 442 } 443 444 if (internal::unpacket_traits<Packet>::masked_store_available) { 445 test::packet_helper<internal::unpacket_traits<Packet>::masked_store_available, Packet> h; 446 unsigned long long max_umask = (0x1ull << PacketSize); 447 448 for (int offset = 0; offset < PacketSize; ++offset) { 449 for (unsigned long long umask = 0; umask < max_umask; ++umask) { 450 internal::pstore(data2, internal::pset1<Packet>(Scalar(0))); 451 h.store(data2, h.loadu(data1 + offset), umask); 452 for (int k = 0; k < PacketSize; ++k) data3[k] = ((umask & (0x1ull << k)) >> k) ? data1[k + offset] : Scalar(0); 453 VERIFY(test::areApprox(data3, data2, PacketSize) && "internal::pstoreu masked"); 454 } 455 } 456 } 457 458 VERIFY((!PacketTraits::Vectorizable) || PacketTraits::HasAdd); 459 VERIFY((!PacketTraits::Vectorizable) || PacketTraits::HasSub); 460 VERIFY((!PacketTraits::Vectorizable) || PacketTraits::HasMul); 461 462 CHECK_CWISE2_IF(PacketTraits::HasAdd, REF_ADD, internal::padd); 463 CHECK_CWISE2_IF(PacketTraits::HasSub, REF_SUB, internal::psub); 464 CHECK_CWISE2_IF(PacketTraits::HasMul, REF_MUL, internal::pmul); 465 CHECK_CWISE2_IF(PacketTraits::HasDiv, REF_DIV, internal::pdiv); 466 467 if (PacketTraits::HasNegate) CHECK_CWISE1(internal::negate, internal::pnegate); 468 CHECK_CWISE1(numext::conj, internal::pconj); 469 470 for (int offset = 0; offset < 3; ++offset) { 471 for (int i = 0; i < PacketSize; ++i) ref[i] = data1[offset]; 472 internal::pstore(data2, internal::pset1<Packet>(data1[offset])); 473 VERIFY(test::areApprox(ref, data2, PacketSize) && "internal::pset1"); 474 } 475 476 { 477 for (int i = 0; i < PacketSize * 4; ++i) ref[i] = data1[i / PacketSize]; 478 Packet A0, A1, A2, A3; 479 internal::pbroadcast4<Packet>(data1, A0, A1, A2, A3); 480 internal::pstore(data2 + 0 * PacketSize, A0); 481 internal::pstore(data2 + 1 * PacketSize, A1); 482 internal::pstore(data2 + 2 * PacketSize, A2); 483 internal::pstore(data2 + 3 * PacketSize, A3); 484 VERIFY(test::areApprox(ref, data2, 4 * PacketSize) && "internal::pbroadcast4"); 485 } 486 487 { 488 for (int i = 0; i < PacketSize * 2; ++i) ref[i] = data1[i / PacketSize]; 489 Packet A0, A1; 490 internal::pbroadcast2<Packet>(data1, A0, A1); 491 internal::pstore(data2 + 0 * PacketSize, A0); 492 internal::pstore(data2 + 1 * PacketSize, A1); 493 VERIFY(test::areApprox(ref, data2, 2 * PacketSize) && "internal::pbroadcast2"); 494 } 495 496 VERIFY(internal::isApprox(data1[0], internal::pfirst(internal::pload<Packet>(data1))) && "internal::pfirst"); 497 498 if (PacketSize > 1) { 499 // apply different offsets to check that ploaddup is robust to unaligned inputs 500 for (int offset = 0; offset < 4; ++offset) { 501 for (int i = 0; i < PacketSize / 2; ++i) ref[2 * i + 0] = ref[2 * i + 1] = data1[offset + i]; 502 internal::pstore(data2, internal::ploaddup<Packet>(data1 + offset)); 503 VERIFY(test::areApprox(ref, data2, PacketSize) && "ploaddup"); 504 } 505 } 506 507 if (PacketSize > 2) { 508 // apply different offsets to check that ploadquad is robust to unaligned inputs 509 for (int offset = 0; offset < 4; ++offset) { 510 for (int i = 0; i < PacketSize / 4; ++i) 511 ref[4 * i + 0] = ref[4 * i + 1] = ref[4 * i + 2] = ref[4 * i + 3] = data1[offset + i]; 512 internal::pstore(data2, internal::ploadquad<Packet>(data1 + offset)); 513 VERIFY(test::areApprox(ref, data2, PacketSize) && "ploadquad"); 514 } 515 } 516 517 ref[0] = Scalar(0); 518 for (int i = 0; i < PacketSize; ++i) ref[0] += data1[i]; 519 VERIFY(test::isApproxAbs(ref[0], internal::predux(internal::pload<Packet>(data1)), refvalue) && "internal::predux"); 520 521 if (!internal::is_same<Packet, typename internal::unpacket_traits<Packet>::half>::value) { 522 int HalfPacketSize = PacketSize > 4 ? PacketSize / 2 : PacketSize; 523 for (int i = 0; i < HalfPacketSize; ++i) ref[i] = Scalar(0); 524 for (int i = 0; i < PacketSize; ++i) ref[i % HalfPacketSize] += data1[i]; 525 internal::pstore(data2, internal::predux_half_dowto4(internal::pload<Packet>(data1))); 526 VERIFY(test::areApprox(ref, data2, HalfPacketSize) && "internal::predux_half_dowto4"); 527 } 528 529 ref[0] = Scalar(1); 530 for (int i = 0; i < PacketSize; ++i) ref[0] = REF_MUL(ref[0], data1[i]); 531 VERIFY(internal::isApprox(ref[0], internal::predux_mul(internal::pload<Packet>(data1))) && "internal::predux_mul"); 532 533 for (int i = 0; i < PacketSize; ++i) ref[i] = data1[PacketSize - i - 1]; 534 internal::pstore(data2, internal::preverse(internal::pload<Packet>(data1))); 535 VERIFY(test::areApprox(ref, data2, PacketSize) && "internal::preverse"); 536 537 internal::PacketBlock<Packet> kernel; 538 for (int i = 0; i < PacketSize; ++i) { 539 kernel.packet[i] = internal::pload<Packet>(data1 + i * PacketSize); 540 } 541 ptranspose(kernel); 542 for (int i = 0; i < PacketSize; ++i) { 543 internal::pstore(data2, kernel.packet[i]); 544 for (int j = 0; j < PacketSize; ++j) { 545 VERIFY(test::isApproxAbs(data2[j], data1[i + j * PacketSize], refvalue) && "ptranspose"); 546 } 547 } 548 549 // GeneralBlockPanelKernel also checks PacketBlock<Packet,(PacketSize%4)==0?4:PacketSize>; 550 if (PacketSize > 4 && PacketSize % 4 == 0) { 551 internal::PacketBlock<Packet, PacketSize%4==0?4:PacketSize> kernel2; 552 for (int i = 0; i < 4; ++i) { 553 kernel2.packet[i] = internal::pload<Packet>(data1 + i * PacketSize); 554 } 555 ptranspose(kernel2); 556 int data_counter = 0; 557 for (int i = 0; i < PacketSize; ++i) { 558 for (int j = 0; j < 4; ++j) { 559 data2[data_counter++] = data1[j*PacketSize + i]; 560 } 561 } 562 for (int i = 0; i < 4; ++i) { 563 internal::pstore(data3, kernel2.packet[i]); 564 for (int j = 0; j < PacketSize; ++j) { 565 VERIFY(test::isApproxAbs(data3[j], data2[i*PacketSize + j], refvalue) && "ptranspose"); 566 } 567 } 568 } 569 570 if (PacketTraits::HasBlend) { 571 Packet thenPacket = internal::pload<Packet>(data1); 572 Packet elsePacket = internal::pload<Packet>(data2); 573 EIGEN_ALIGN_MAX internal::Selector<PacketSize> selector; 574 for (int i = 0; i < PacketSize; ++i) { 575 selector.select[i] = i; 576 } 577 578 Packet blend = internal::pblend(selector, thenPacket, elsePacket); 579 EIGEN_ALIGN_MAX Scalar result[size]; 580 internal::pstore(result, blend); 581 for (int i = 0; i < PacketSize; ++i) { 582 VERIFY(test::isApproxAbs(result[i], (selector.select[i] ? data1[i] : data2[i]), refvalue)); 583 } 584 } 585 586 { 587 for (int i = 0; i < PacketSize; ++i) { 588 // "if" mask 589 unsigned char v = internal::random<bool>() ? 0xff : 0; 590 char* bytes = (char*)(data1 + i); 591 for (int k = 0; k < int(sizeof(Scalar)); ++k) { 592 bytes[k] = v; 593 } 594 // "then" packet 595 data1[i + PacketSize] = internal::random<Scalar>(); 596 // "else" packet 597 data1[i + 2 * PacketSize] = internal::random<Scalar>(); 598 } 599 CHECK_CWISE3_IF(true, internal::pselect, internal::pselect); 600 } 601 602 for (int i = 0; i < size; ++i) { 603 data1[i] = internal::random<Scalar>(); 604 } 605 CHECK_CWISE1(internal::pzero, internal::pzero); 606 CHECK_CWISE2_IF(true, internal::por, internal::por); 607 CHECK_CWISE2_IF(true, internal::pxor, internal::pxor); 608 CHECK_CWISE2_IF(true, internal::pand, internal::pand); 609 610 packetmath_boolean_mask_ops<Scalar, Packet>(); 611 packetmath_pcast_ops_runner<Scalar, Packet>::run(); 612 packetmath_minus_zero_add<Scalar, Packet>(); 613 614 for (int i = 0; i < size; ++i) { 615 data1[i] = numext::abs(internal::random<Scalar>()); 616 } 617 CHECK_CWISE1_IF(PacketTraits::HasSqrt, numext::sqrt, internal::psqrt); 618 CHECK_CWISE1_IF(PacketTraits::HasRsqrt, numext::rsqrt, internal::prsqrt); 619 } 620 621 // Notice that this definition works for complex types as well. 622 // c++11 has std::log2 for real, but not for complex types. 623 template <typename Scalar> 624 Scalar log2(Scalar x) { 625 return Scalar(EIGEN_LOG2E) * std::log(x); 626 } 627 628 template <typename Scalar, typename Packet> 629 void packetmath_real() { 630 typedef internal::packet_traits<Scalar> PacketTraits; 631 const int PacketSize = internal::unpacket_traits<Packet>::size; 632 633 const int size = PacketSize * 4; 634 EIGEN_ALIGN_MAX Scalar data1[PacketSize * 4]; 635 EIGEN_ALIGN_MAX Scalar data2[PacketSize * 4]; 636 EIGEN_ALIGN_MAX Scalar ref[PacketSize * 4]; 637 638 for (int i = 0; i < size; ++i) { 639 data1[i] = Scalar(internal::random<double>(0, 1) * std::pow(10., internal::random<double>(-6, 6))); 640 data2[i] = Scalar(internal::random<double>(0, 1) * std::pow(10., internal::random<double>(-6, 6))); 641 } 642 643 if (internal::random<float>(0, 1) < 0.1f) data1[internal::random<int>(0, PacketSize)] = Scalar(0); 644 645 CHECK_CWISE1_IF(PacketTraits::HasLog, std::log, internal::plog); 646 CHECK_CWISE1_IF(PacketTraits::HasLog, log2, internal::plog2); 647 CHECK_CWISE1_IF(PacketTraits::HasRsqrt, numext::rsqrt, internal::prsqrt); 648 649 for (int i = 0; i < size; ++i) { 650 data1[i] = Scalar(internal::random<double>(-1, 1) * std::pow(10., internal::random<double>(-3, 3))); 651 data2[i] = Scalar(internal::random<double>(-1, 1) * std::pow(10., internal::random<double>(-3, 3))); 652 } 653 CHECK_CWISE1_IF(PacketTraits::HasSin, std::sin, internal::psin); 654 CHECK_CWISE1_IF(PacketTraits::HasCos, std::cos, internal::pcos); 655 CHECK_CWISE1_IF(PacketTraits::HasTan, std::tan, internal::ptan); 656 657 CHECK_CWISE1_EXACT_IF(PacketTraits::HasRound, numext::round, internal::pround); 658 CHECK_CWISE1_EXACT_IF(PacketTraits::HasCeil, numext::ceil, internal::pceil); 659 CHECK_CWISE1_EXACT_IF(PacketTraits::HasFloor, numext::floor, internal::pfloor); 660 CHECK_CWISE1_EXACT_IF(PacketTraits::HasRint, numext::rint, internal::print); 661 662 packetmath_boolean_mask_ops_real<Scalar,Packet>(); 663 664 // Rounding edge cases. 665 if (PacketTraits::HasRound || PacketTraits::HasCeil || PacketTraits::HasFloor || PacketTraits::HasRint) { 666 typedef typename internal::make_integer<Scalar>::type IntType; 667 // Start with values that cannot fit inside an integer, work down to less than one. 668 Scalar val = numext::mini( 669 Scalar(2) * static_cast<Scalar>(NumTraits<IntType>::highest()), 670 NumTraits<Scalar>::highest()); 671 std::vector<Scalar> values; 672 while (val > Scalar(0.25)) { 673 // Cover both even and odd, positive and negative cases. 674 values.push_back(val); 675 values.push_back(val + Scalar(0.3)); 676 values.push_back(val + Scalar(0.5)); 677 values.push_back(val + Scalar(0.8)); 678 values.push_back(val + Scalar(1)); 679 values.push_back(val + Scalar(1.3)); 680 values.push_back(val + Scalar(1.5)); 681 values.push_back(val + Scalar(1.8)); 682 values.push_back(-val); 683 values.push_back(-val - Scalar(0.3)); 684 values.push_back(-val - Scalar(0.5)); 685 values.push_back(-val - Scalar(0.8)); 686 values.push_back(-val - Scalar(1)); 687 values.push_back(-val - Scalar(1.3)); 688 values.push_back(-val - Scalar(1.5)); 689 values.push_back(-val - Scalar(1.8)); 690 values.push_back(Scalar(-1.5) + val); // Bug 1785. 691 val = val / Scalar(2); 692 } 693 values.push_back(NumTraits<Scalar>::infinity()); 694 values.push_back(-NumTraits<Scalar>::infinity()); 695 values.push_back(NumTraits<Scalar>::quiet_NaN()); 696 697 for (size_t k=0; k<values.size(); ++k) { 698 data1[0] = values[k]; 699 CHECK_CWISE1_EXACT_IF(PacketTraits::HasRound, numext::round, internal::pround); 700 CHECK_CWISE1_EXACT_IF(PacketTraits::HasCeil, numext::ceil, internal::pceil); 701 CHECK_CWISE1_EXACT_IF(PacketTraits::HasFloor, numext::floor, internal::pfloor); 702 CHECK_CWISE1_EXACT_IF(PacketTraits::HasRint, numext::rint, internal::print); 703 } 704 } 705 706 for (int i = 0; i < size; ++i) { 707 data1[i] = Scalar(internal::random<double>(-1, 1)); 708 data2[i] = Scalar(internal::random<double>(-1, 1)); 709 } 710 CHECK_CWISE1_IF(PacketTraits::HasASin, std::asin, internal::pasin); 711 CHECK_CWISE1_IF(PacketTraits::HasACos, std::acos, internal::pacos); 712 713 for (int i = 0; i < size; ++i) { 714 data1[i] = Scalar(internal::random<double>(-87, 88)); 715 data2[i] = Scalar(internal::random<double>(-87, 88)); 716 } 717 CHECK_CWISE1_IF(PacketTraits::HasExp, std::exp, internal::pexp); 718 719 CHECK_CWISE1_BYREF1_IF(PacketTraits::HasExp, REF_FREXP, internal::pfrexp); 720 if (PacketTraits::HasExp) { 721 // Check denormals: 722 for (int j=0; j<3; ++j) { 723 data1[0] = Scalar(std::ldexp(1, NumTraits<Scalar>::min_exponent()-j)); 724 CHECK_CWISE1_BYREF1_IF(PacketTraits::HasExp, REF_FREXP, internal::pfrexp); 725 data1[0] = -data1[0]; 726 CHECK_CWISE1_BYREF1_IF(PacketTraits::HasExp, REF_FREXP, internal::pfrexp); 727 } 728 729 // zero 730 data1[0] = Scalar(0); 731 CHECK_CWISE1_BYREF1_IF(PacketTraits::HasExp, REF_FREXP, internal::pfrexp); 732 733 // inf and NaN only compare output fraction, not exponent. 734 test::packet_helper<PacketTraits::HasExp,Packet> h; 735 Packet pout; 736 Scalar sout; 737 Scalar special[] = { NumTraits<Scalar>::infinity(), 738 -NumTraits<Scalar>::infinity(), 739 NumTraits<Scalar>::quiet_NaN()}; 740 for (int i=0; i<3; ++i) { 741 data1[0] = special[i]; 742 ref[0] = Scalar(REF_FREXP(data1[0], ref[PacketSize])); 743 h.store(data2, internal::pfrexp(h.load(data1), h.forward_reference(pout, sout))); 744 VERIFY(test::areApprox(ref, data2, 1) && "internal::pfrexp"); 745 } 746 } 747 748 for (int i = 0; i < PacketSize; ++i) { 749 data1[i] = Scalar(internal::random<double>(-1, 1)); 750 data2[i] = Scalar(internal::random<double>(-1, 1)); 751 } 752 for (int i = 0; i < PacketSize; ++i) { 753 data1[i+PacketSize] = Scalar(internal::random<int>(-4, 4)); 754 data2[i+PacketSize] = Scalar(internal::random<double>(-4, 4)); 755 } 756 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp); 757 if (PacketTraits::HasExp) { 758 data1[0] = Scalar(-1); 759 // underflow to zero 760 data1[PacketSize] = Scalar(NumTraits<Scalar>::min_exponent()-55); 761 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp); 762 // overflow to inf 763 data1[PacketSize] = Scalar(NumTraits<Scalar>::max_exponent()+10); 764 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp); 765 // NaN stays NaN 766 data1[0] = NumTraits<Scalar>::quiet_NaN(); 767 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp); 768 VERIFY((numext::isnan)(data2[0])); 769 // inf stays inf 770 data1[0] = NumTraits<Scalar>::infinity(); 771 data1[PacketSize] = Scalar(NumTraits<Scalar>::min_exponent()-10); 772 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp); 773 // zero stays zero 774 data1[0] = Scalar(0); 775 data1[PacketSize] = Scalar(NumTraits<Scalar>::max_exponent()+10); 776 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp); 777 // Small number big exponent. 778 data1[0] = Scalar(std::ldexp(Scalar(1.0), NumTraits<Scalar>::min_exponent()-1)); 779 data1[PacketSize] = Scalar(-NumTraits<Scalar>::min_exponent() 780 +NumTraits<Scalar>::max_exponent()); 781 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp); 782 // Big number small exponent. 783 data1[0] = Scalar(std::ldexp(Scalar(1.0), NumTraits<Scalar>::max_exponent()-1)); 784 data1[PacketSize] = Scalar(+NumTraits<Scalar>::min_exponent() 785 -NumTraits<Scalar>::max_exponent()); 786 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp); 787 } 788 789 for (int i = 0; i < size; ++i) { 790 data1[i] = Scalar(internal::random<double>(-1, 1) * std::pow(10., internal::random<double>(-6, 6))); 791 data2[i] = Scalar(internal::random<double>(-1, 1) * std::pow(10., internal::random<double>(-6, 6))); 792 } 793 data1[0] = Scalar(1e-20); 794 CHECK_CWISE1_IF(PacketTraits::HasTanh, std::tanh, internal::ptanh); 795 if (PacketTraits::HasExp && PacketSize >= 2) { 796 const Scalar small = NumTraits<Scalar>::epsilon(); 797 data1[0] = NumTraits<Scalar>::quiet_NaN(); 798 data1[1] = small; 799 test::packet_helper<PacketTraits::HasExp, Packet> h; 800 h.store(data2, internal::pexp(h.load(data1))); 801 VERIFY((numext::isnan)(data2[0])); 802 // TODO(rmlarsen): Re-enable for bfloat16. 803 if (!internal::is_same<Scalar, bfloat16>::value) { 804 VERIFY_IS_APPROX(std::exp(small), data2[1]); 805 } 806 807 data1[0] = -small; 808 data1[1] = Scalar(0); 809 h.store(data2, internal::pexp(h.load(data1))); 810 // TODO(rmlarsen): Re-enable for bfloat16. 811 if (!internal::is_same<Scalar, bfloat16>::value) { 812 VERIFY_IS_APPROX(std::exp(-small), data2[0]); 813 } 814 VERIFY_IS_EQUAL(std::exp(Scalar(0)), data2[1]); 815 816 data1[0] = (std::numeric_limits<Scalar>::min)(); 817 data1[1] = -(std::numeric_limits<Scalar>::min)(); 818 h.store(data2, internal::pexp(h.load(data1))); 819 VERIFY_IS_APPROX(std::exp((std::numeric_limits<Scalar>::min)()), data2[0]); 820 VERIFY_IS_APPROX(std::exp(-(std::numeric_limits<Scalar>::min)()), data2[1]); 821 822 data1[0] = std::numeric_limits<Scalar>::denorm_min(); 823 data1[1] = -std::numeric_limits<Scalar>::denorm_min(); 824 h.store(data2, internal::pexp(h.load(data1))); 825 VERIFY_IS_APPROX(std::exp(std::numeric_limits<Scalar>::denorm_min()), data2[0]); 826 VERIFY_IS_APPROX(std::exp(-std::numeric_limits<Scalar>::denorm_min()), data2[1]); 827 } 828 829 if (PacketTraits::HasTanh) { 830 // NOTE this test migh fail with GCC prior to 6.3, see MathFunctionsImpl.h for details. 831 data1[0] = NumTraits<Scalar>::quiet_NaN(); 832 test::packet_helper<internal::packet_traits<Scalar>::HasTanh, Packet> h; 833 h.store(data2, internal::ptanh(h.load(data1))); 834 VERIFY((numext::isnan)(data2[0])); 835 } 836 837 if (PacketTraits::HasExp) { 838 internal::scalar_logistic_op<Scalar> logistic; 839 for (int i = 0; i < size; ++i) { 840 data1[i] = Scalar(internal::random<double>(-20, 20)); 841 } 842 843 test::packet_helper<PacketTraits::HasExp, Packet> h; 844 h.store(data2, logistic.packetOp(h.load(data1))); 845 for (int i = 0; i < PacketSize; ++i) { 846 VERIFY_IS_APPROX(data2[i], logistic(data1[i])); 847 } 848 } 849 850 #if EIGEN_HAS_C99_MATH && (EIGEN_COMP_CXXVER >= 11) 851 data1[0] = NumTraits<Scalar>::infinity(); 852 data1[1] = Scalar(-1); 853 CHECK_CWISE1_IF(PacketTraits::HasLog1p, std::log1p, internal::plog1p); 854 data1[0] = NumTraits<Scalar>::infinity(); 855 data1[1] = -NumTraits<Scalar>::infinity(); 856 CHECK_CWISE1_IF(PacketTraits::HasExpm1, std::expm1, internal::pexpm1); 857 #endif 858 859 if (PacketSize >= 2) { 860 data1[0] = NumTraits<Scalar>::quiet_NaN(); 861 data1[1] = NumTraits<Scalar>::epsilon(); 862 if (PacketTraits::HasLog) { 863 test::packet_helper<PacketTraits::HasLog, Packet> h; 864 h.store(data2, internal::plog(h.load(data1))); 865 VERIFY((numext::isnan)(data2[0])); 866 // TODO(cantonios): Re-enable for bfloat16. 867 if (!internal::is_same<Scalar, bfloat16>::value) { 868 VERIFY_IS_APPROX(std::log(data1[1]), data2[1]); 869 } 870 871 data1[0] = -NumTraits<Scalar>::epsilon(); 872 data1[1] = Scalar(0); 873 h.store(data2, internal::plog(h.load(data1))); 874 VERIFY((numext::isnan)(data2[0])); 875 VERIFY_IS_EQUAL(std::log(Scalar(0)), data2[1]); 876 877 data1[0] = (std::numeric_limits<Scalar>::min)(); 878 data1[1] = -(std::numeric_limits<Scalar>::min)(); 879 h.store(data2, internal::plog(h.load(data1))); 880 // TODO(cantonios): Re-enable for bfloat16. 881 if (!internal::is_same<Scalar, bfloat16>::value) { 882 VERIFY_IS_APPROX(std::log((std::numeric_limits<Scalar>::min)()), data2[0]); 883 } 884 VERIFY((numext::isnan)(data2[1])); 885 886 // Note: 32-bit arm always flushes denorms to zero. 887 #if !EIGEN_ARCH_ARM 888 if (std::numeric_limits<Scalar>::has_denorm == std::denorm_present) { 889 data1[0] = std::numeric_limits<Scalar>::denorm_min(); 890 data1[1] = -std::numeric_limits<Scalar>::denorm_min(); 891 h.store(data2, internal::plog(h.load(data1))); 892 // TODO(rmlarsen): Reenable. 893 // VERIFY_IS_EQUAL(std::log(std::numeric_limits<Scalar>::denorm_min()), data2[0]); 894 VERIFY((numext::isnan)(data2[1])); 895 } 896 #endif 897 898 data1[0] = Scalar(-1.0f); 899 h.store(data2, internal::plog(h.load(data1))); 900 VERIFY((numext::isnan)(data2[0])); 901 902 data1[0] = NumTraits<Scalar>::infinity(); 903 h.store(data2, internal::plog(h.load(data1))); 904 VERIFY((numext::isinf)(data2[0])); 905 } 906 if (PacketTraits::HasLog1p) { 907 test::packet_helper<PacketTraits::HasLog1p, Packet> h; 908 data1[0] = Scalar(-2); 909 data1[1] = -NumTraits<Scalar>::infinity(); 910 h.store(data2, internal::plog1p(h.load(data1))); 911 VERIFY((numext::isnan)(data2[0])); 912 VERIFY((numext::isnan)(data2[1])); 913 } 914 if (PacketTraits::HasSqrt) { 915 test::packet_helper<PacketTraits::HasSqrt, Packet> h; 916 data1[0] = Scalar(-1.0f); 917 if (std::numeric_limits<Scalar>::has_denorm == std::denorm_present) { 918 data1[1] = -std::numeric_limits<Scalar>::denorm_min(); 919 } else { 920 data1[1] = -NumTraits<Scalar>::epsilon(); 921 } 922 h.store(data2, internal::psqrt(h.load(data1))); 923 VERIFY((numext::isnan)(data2[0])); 924 VERIFY((numext::isnan)(data2[1])); 925 } 926 // TODO(rmlarsen): Re-enable for half and bfloat16. 927 if (PacketTraits::HasCos 928 && !internal::is_same<Scalar, half>::value 929 && !internal::is_same<Scalar, bfloat16>::value) { 930 test::packet_helper<PacketTraits::HasCos, Packet> h; 931 for (Scalar k = Scalar(1); k < Scalar(10000) / NumTraits<Scalar>::epsilon(); k *= Scalar(2)) { 932 for (int k1 = 0; k1 <= 1; ++k1) { 933 data1[0] = Scalar((2 * double(k) + k1) * double(EIGEN_PI) / 2 * internal::random<double>(0.8, 1.2)); 934 data1[1] = Scalar((2 * double(k) + 2 + k1) * double(EIGEN_PI) / 2 * internal::random<double>(0.8, 1.2)); 935 h.store(data2, internal::pcos(h.load(data1))); 936 h.store(data2 + PacketSize, internal::psin(h.load(data1))); 937 VERIFY(data2[0] <= Scalar(1.) && data2[0] >= Scalar(-1.)); 938 VERIFY(data2[1] <= Scalar(1.) && data2[1] >= Scalar(-1.)); 939 VERIFY(data2[PacketSize + 0] <= Scalar(1.) && data2[PacketSize + 0] >= Scalar(-1.)); 940 VERIFY(data2[PacketSize + 1] <= Scalar(1.) && data2[PacketSize + 1] >= Scalar(-1.)); 941 942 VERIFY_IS_APPROX(data2[0], std::cos(data1[0])); 943 VERIFY_IS_APPROX(data2[1], std::cos(data1[1])); 944 VERIFY_IS_APPROX(data2[PacketSize + 0], std::sin(data1[0])); 945 VERIFY_IS_APPROX(data2[PacketSize + 1], std::sin(data1[1])); 946 947 VERIFY_IS_APPROX(numext::abs2(data2[0]) + numext::abs2(data2[PacketSize + 0]), Scalar(1)); 948 VERIFY_IS_APPROX(numext::abs2(data2[1]) + numext::abs2(data2[PacketSize + 1]), Scalar(1)); 949 } 950 } 951 952 data1[0] = NumTraits<Scalar>::infinity(); 953 data1[1] = -NumTraits<Scalar>::infinity(); 954 h.store(data2, internal::psin(h.load(data1))); 955 VERIFY((numext::isnan)(data2[0])); 956 VERIFY((numext::isnan)(data2[1])); 957 958 h.store(data2, internal::pcos(h.load(data1))); 959 VERIFY((numext::isnan)(data2[0])); 960 VERIFY((numext::isnan)(data2[1])); 961 962 data1[0] = NumTraits<Scalar>::quiet_NaN(); 963 h.store(data2, internal::psin(h.load(data1))); 964 VERIFY((numext::isnan)(data2[0])); 965 h.store(data2, internal::pcos(h.load(data1))); 966 VERIFY((numext::isnan)(data2[0])); 967 968 data1[0] = -Scalar(0.); 969 h.store(data2, internal::psin(h.load(data1))); 970 VERIFY(internal::biteq(data2[0], data1[0])); 971 h.store(data2, internal::pcos(h.load(data1))); 972 VERIFY_IS_EQUAL(data2[0], Scalar(1)); 973 } 974 } 975 } 976 977 #define CAST_CHECK_CWISE1_IF(COND, REFOP, POP, SCALAR, REFTYPE) if(COND) { \ 978 test::packet_helper<COND,Packet> h; \ 979 for (int i=0; i<PacketSize; ++i) \ 980 ref[i] = SCALAR(REFOP(static_cast<REFTYPE>(data1[i]))); \ 981 h.store(data2, POP(h.load(data1))); \ 982 VERIFY(test::areApprox(ref, data2, PacketSize) && #POP); \ 983 } 984 985 template <typename Scalar> 986 Scalar propagate_nan_max(const Scalar& a, const Scalar& b) { 987 if ((numext::isnan)(a)) return a; 988 if ((numext::isnan)(b)) return b; 989 return (numext::maxi)(a,b); 990 } 991 992 template <typename Scalar> 993 Scalar propagate_nan_min(const Scalar& a, const Scalar& b) { 994 if ((numext::isnan)(a)) return a; 995 if ((numext::isnan)(b)) return b; 996 return (numext::mini)(a,b); 997 } 998 999 template <typename Scalar> 1000 Scalar propagate_number_max(const Scalar& a, const Scalar& b) { 1001 if ((numext::isnan)(a)) return b; 1002 if ((numext::isnan)(b)) return a; 1003 return (numext::maxi)(a,b); 1004 } 1005 1006 template <typename Scalar> 1007 Scalar propagate_number_min(const Scalar& a, const Scalar& b) { 1008 if ((numext::isnan)(a)) return b; 1009 if ((numext::isnan)(b)) return a; 1010 return (numext::mini)(a,b); 1011 } 1012 1013 template <typename Scalar, typename Packet> 1014 void packetmath_notcomplex() { 1015 typedef internal::packet_traits<Scalar> PacketTraits; 1016 const int PacketSize = internal::unpacket_traits<Packet>::size; 1017 1018 EIGEN_ALIGN_MAX Scalar data1[PacketSize * 4]; 1019 EIGEN_ALIGN_MAX Scalar data2[PacketSize * 4]; 1020 EIGEN_ALIGN_MAX Scalar ref[PacketSize * 4]; 1021 1022 Array<Scalar, Dynamic, 1>::Map(data1, PacketSize * 4).setRandom(); 1023 1024 VERIFY((!PacketTraits::Vectorizable) || PacketTraits::HasMin); 1025 VERIFY((!PacketTraits::Vectorizable) || PacketTraits::HasMax); 1026 1027 CHECK_CWISE2_IF(PacketTraits::HasMin, (std::min), internal::pmin); 1028 CHECK_CWISE2_IF(PacketTraits::HasMax, (std::max), internal::pmax); 1029 1030 CHECK_CWISE2_IF(PacketTraits::HasMin, propagate_number_min, internal::pmin<PropagateNumbers>); 1031 CHECK_CWISE2_IF(PacketTraits::HasMax, propagate_number_max, internal::pmax<PropagateNumbers>); 1032 CHECK_CWISE1(numext::abs, internal::pabs); 1033 CHECK_CWISE2_IF(PacketTraits::HasAbsDiff, REF_ABS_DIFF, internal::pabsdiff); 1034 1035 ref[0] = data1[0]; 1036 for (int i = 0; i < PacketSize; ++i) ref[0] = internal::pmin(ref[0], data1[i]); 1037 VERIFY(internal::isApprox(ref[0], internal::predux_min(internal::pload<Packet>(data1))) && "internal::predux_min"); 1038 ref[0] = data1[0]; 1039 for (int i = 0; i < PacketSize; ++i) ref[0] = internal::pmax(ref[0], data1[i]); 1040 VERIFY(internal::isApprox(ref[0], internal::predux_max(internal::pload<Packet>(data1))) && "internal::predux_max"); 1041 1042 for (int i = 0; i < PacketSize; ++i) ref[i] = data1[0] + Scalar(i); 1043 internal::pstore(data2, internal::plset<Packet>(data1[0])); 1044 VERIFY(test::areApprox(ref, data2, PacketSize) && "internal::plset"); 1045 1046 { 1047 unsigned char* data1_bits = reinterpret_cast<unsigned char*>(data1); 1048 // predux_all - not needed yet 1049 // for (unsigned int i=0; i<PacketSize*sizeof(Scalar); ++i) data1_bits[i] = 0xff; 1050 // VERIFY(internal::predux_all(internal::pload<Packet>(data1)) && "internal::predux_all(1111)"); 1051 // for(int k=0; k<PacketSize; ++k) 1052 // { 1053 // for (unsigned int i=0; i<sizeof(Scalar); ++i) data1_bits[k*sizeof(Scalar)+i] = 0x0; 1054 // VERIFY( (!internal::predux_all(internal::pload<Packet>(data1))) && "internal::predux_all(0101)"); 1055 // for (unsigned int i=0; i<sizeof(Scalar); ++i) data1_bits[k*sizeof(Scalar)+i] = 0xff; 1056 // } 1057 1058 // predux_any 1059 for (unsigned int i = 0; i < PacketSize * sizeof(Scalar); ++i) data1_bits[i] = 0x0; 1060 VERIFY((!internal::predux_any(internal::pload<Packet>(data1))) && "internal::predux_any(0000)"); 1061 for (int k = 0; k < PacketSize; ++k) { 1062 for (unsigned int i = 0; i < sizeof(Scalar); ++i) data1_bits[k * sizeof(Scalar) + i] = 0xff; 1063 VERIFY(internal::predux_any(internal::pload<Packet>(data1)) && "internal::predux_any(0101)"); 1064 for (unsigned int i = 0; i < sizeof(Scalar); ++i) data1_bits[k * sizeof(Scalar) + i] = 0x00; 1065 } 1066 } 1067 1068 1069 // Test NaN propagation. 1070 if (!NumTraits<Scalar>::IsInteger) { 1071 // Test reductions with no NaNs. 1072 ref[0] = data1[0]; 1073 for (int i = 0; i < PacketSize; ++i) ref[0] = internal::pmin<PropagateNumbers>(ref[0], data1[i]); 1074 VERIFY(internal::isApprox(ref[0], internal::predux_min<PropagateNumbers>(internal::pload<Packet>(data1))) && "internal::predux_min<PropagateNumbers>"); 1075 ref[0] = data1[0]; 1076 for (int i = 0; i < PacketSize; ++i) ref[0] = internal::pmin<PropagateNaN>(ref[0], data1[i]); 1077 VERIFY(internal::isApprox(ref[0], internal::predux_min<PropagateNaN>(internal::pload<Packet>(data1))) && "internal::predux_min<PropagateNaN>"); 1078 ref[0] = data1[0]; 1079 for (int i = 0; i < PacketSize; ++i) ref[0] = internal::pmax<PropagateNumbers>(ref[0], data1[i]); 1080 VERIFY(internal::isApprox(ref[0], internal::predux_max<PropagateNumbers>(internal::pload<Packet>(data1))) && "internal::predux_max<PropagateNumbers>"); 1081 ref[0] = data1[0]; 1082 for (int i = 0; i < PacketSize; ++i) ref[0] = internal::pmax<PropagateNaN>(ref[0], data1[i]); 1083 VERIFY(internal::isApprox(ref[0], internal::predux_max<PropagateNaN>(internal::pload<Packet>(data1))) && "internal::predux_max<PropagateNumbers>"); 1084 // A single NaN. 1085 const size_t index = std::numeric_limits<size_t>::quiet_NaN() % PacketSize; 1086 data1[index] = NumTraits<Scalar>::quiet_NaN(); 1087 VERIFY(PacketSize==1 || !(numext::isnan)(internal::predux_min<PropagateNumbers>(internal::pload<Packet>(data1)))); 1088 VERIFY((numext::isnan)(internal::predux_min<PropagateNaN>(internal::pload<Packet>(data1)))); 1089 VERIFY(PacketSize==1 || !(numext::isnan)(internal::predux_max<PropagateNumbers>(internal::pload<Packet>(data1)))); 1090 VERIFY((numext::isnan)(internal::predux_max<PropagateNaN>(internal::pload<Packet>(data1)))); 1091 // All NaNs. 1092 for (int i = 0; i < 4 * PacketSize; ++i) data1[i] = NumTraits<Scalar>::quiet_NaN(); 1093 VERIFY((numext::isnan)(internal::predux_min<PropagateNumbers>(internal::pload<Packet>(data1)))); 1094 VERIFY((numext::isnan)(internal::predux_min<PropagateNaN>(internal::pload<Packet>(data1)))); 1095 VERIFY((numext::isnan)(internal::predux_max<PropagateNumbers>(internal::pload<Packet>(data1)))); 1096 VERIFY((numext::isnan)(internal::predux_max<PropagateNaN>(internal::pload<Packet>(data1)))); 1097 1098 // Test NaN propagation for coefficient-wise min and max. 1099 for (int i = 0; i < PacketSize; ++i) { 1100 data1[i] = internal::random<bool>() ? NumTraits<Scalar>::quiet_NaN() : Scalar(0); 1101 data1[i + PacketSize] = internal::random<bool>() ? NumTraits<Scalar>::quiet_NaN() : Scalar(0); 1102 } 1103 // Note: NaN propagation is implementation defined for pmin/pmax, so we do not test it here. 1104 CHECK_CWISE2_IF(PacketTraits::HasMin, propagate_number_min, (internal::pmin<PropagateNumbers>)); 1105 CHECK_CWISE2_IF(PacketTraits::HasMax, propagate_number_max, internal::pmax<PropagateNumbers>); 1106 CHECK_CWISE2_IF(PacketTraits::HasMin, propagate_nan_min, (internal::pmin<PropagateNaN>)); 1107 CHECK_CWISE2_IF(PacketTraits::HasMax, propagate_nan_max, internal::pmax<PropagateNaN>); 1108 } 1109 1110 packetmath_boolean_mask_ops_notcomplex<Scalar, Packet>(); 1111 } 1112 1113 template <typename Scalar, typename Packet, bool ConjLhs, bool ConjRhs> 1114 void test_conj_helper(Scalar* data1, Scalar* data2, Scalar* ref, Scalar* pval) { 1115 const int PacketSize = internal::unpacket_traits<Packet>::size; 1116 1117 internal::conj_if<ConjLhs> cj0; 1118 internal::conj_if<ConjRhs> cj1; 1119 internal::conj_helper<Scalar, Scalar, ConjLhs, ConjRhs> cj; 1120 internal::conj_helper<Packet, Packet, ConjLhs, ConjRhs> pcj; 1121 1122 for (int i = 0; i < PacketSize; ++i) { 1123 ref[i] = cj0(data1[i]) * cj1(data2[i]); 1124 VERIFY(internal::isApprox(ref[i], cj.pmul(data1[i], data2[i])) && "conj_helper pmul"); 1125 } 1126 internal::pstore(pval, pcj.pmul(internal::pload<Packet>(data1), internal::pload<Packet>(data2))); 1127 VERIFY(test::areApprox(ref, pval, PacketSize) && "conj_helper pmul"); 1128 1129 for (int i = 0; i < PacketSize; ++i) { 1130 Scalar tmp = ref[i]; 1131 ref[i] += cj0(data1[i]) * cj1(data2[i]); 1132 VERIFY(internal::isApprox(ref[i], cj.pmadd(data1[i], data2[i], tmp)) && "conj_helper pmadd"); 1133 } 1134 internal::pstore( 1135 pval, pcj.pmadd(internal::pload<Packet>(data1), internal::pload<Packet>(data2), internal::pload<Packet>(pval))); 1136 VERIFY(test::areApprox(ref, pval, PacketSize) && "conj_helper pmadd"); 1137 } 1138 1139 template <typename Scalar, typename Packet> 1140 void packetmath_complex() { 1141 typedef internal::packet_traits<Scalar> PacketTraits; 1142 typedef typename Scalar::value_type RealScalar; 1143 const int PacketSize = internal::unpacket_traits<Packet>::size; 1144 1145 const int size = PacketSize * 4; 1146 EIGEN_ALIGN_MAX Scalar data1[PacketSize * 4]; 1147 EIGEN_ALIGN_MAX Scalar data2[PacketSize * 4]; 1148 EIGEN_ALIGN_MAX Scalar ref[PacketSize * 4]; 1149 EIGEN_ALIGN_MAX Scalar pval[PacketSize * 4]; 1150 1151 for (int i = 0; i < size; ++i) { 1152 data1[i] = internal::random<Scalar>() * Scalar(1e2); 1153 data2[i] = internal::random<Scalar>() * Scalar(1e2); 1154 } 1155 1156 test_conj_helper<Scalar, Packet, false, false>(data1, data2, ref, pval); 1157 test_conj_helper<Scalar, Packet, false, true>(data1, data2, ref, pval); 1158 test_conj_helper<Scalar, Packet, true, false>(data1, data2, ref, pval); 1159 test_conj_helper<Scalar, Packet, true, true>(data1, data2, ref, pval); 1160 1161 // Test pcplxflip. 1162 { 1163 for (int i = 0; i < PacketSize; ++i) ref[i] = Scalar(std::imag(data1[i]), std::real(data1[i])); 1164 internal::pstore(pval, internal::pcplxflip(internal::pload<Packet>(data1))); 1165 VERIFY(test::areApprox(ref, pval, PacketSize) && "pcplxflip"); 1166 } 1167 1168 if (PacketTraits::HasSqrt) { 1169 for (int i = 0; i < size; ++i) { 1170 data1[i] = Scalar(internal::random<RealScalar>(), internal::random<RealScalar>()); 1171 } 1172 CHECK_CWISE1_N(numext::sqrt, internal::psqrt, size); 1173 1174 // Test misc. corner cases. 1175 const RealScalar zero = RealScalar(0); 1176 const RealScalar one = RealScalar(1); 1177 const RealScalar inf = std::numeric_limits<RealScalar>::infinity(); 1178 const RealScalar nan = std::numeric_limits<RealScalar>::quiet_NaN(); 1179 data1[0] = Scalar(zero, zero); 1180 data1[1] = Scalar(-zero, zero); 1181 data1[2] = Scalar(one, zero); 1182 data1[3] = Scalar(zero, one); 1183 CHECK_CWISE1_N(numext::sqrt, internal::psqrt, 4); 1184 data1[0] = Scalar(-one, zero); 1185 data1[1] = Scalar(zero, -one); 1186 data1[2] = Scalar(one, one); 1187 data1[3] = Scalar(-one, -one); 1188 CHECK_CWISE1_N(numext::sqrt, internal::psqrt, 4); 1189 data1[0] = Scalar(inf, zero); 1190 data1[1] = Scalar(zero, inf); 1191 data1[2] = Scalar(-inf, zero); 1192 data1[3] = Scalar(zero, -inf); 1193 CHECK_CWISE1_N(numext::sqrt, internal::psqrt, 4); 1194 data1[0] = Scalar(inf, inf); 1195 data1[1] = Scalar(-inf, inf); 1196 data1[2] = Scalar(inf, -inf); 1197 data1[3] = Scalar(-inf, -inf); 1198 CHECK_CWISE1_N(numext::sqrt, internal::psqrt, 4); 1199 data1[0] = Scalar(nan, zero); 1200 data1[1] = Scalar(zero, nan); 1201 data1[2] = Scalar(nan, one); 1202 data1[3] = Scalar(one, nan); 1203 CHECK_CWISE1_N(numext::sqrt, internal::psqrt, 4); 1204 data1[0] = Scalar(nan, nan); 1205 data1[1] = Scalar(inf, nan); 1206 data1[2] = Scalar(nan, inf); 1207 data1[3] = Scalar(-inf, nan); 1208 CHECK_CWISE1_N(numext::sqrt, internal::psqrt, 4); 1209 } 1210 } 1211 1212 template <typename Scalar, typename Packet> 1213 void packetmath_scatter_gather() { 1214 typedef typename NumTraits<Scalar>::Real RealScalar; 1215 const int PacketSize = internal::unpacket_traits<Packet>::size; 1216 EIGEN_ALIGN_MAX Scalar data1[PacketSize]; 1217 RealScalar refvalue = RealScalar(0); 1218 for (int i = 0; i < PacketSize; ++i) { 1219 data1[i] = internal::random<Scalar>() / RealScalar(PacketSize); 1220 } 1221 1222 int stride = internal::random<int>(1, 20); 1223 1224 // Buffer of zeros. 1225 EIGEN_ALIGN_MAX Scalar buffer[PacketSize * 20] = {}; 1226 1227 Packet packet = internal::pload<Packet>(data1); 1228 internal::pscatter<Scalar, Packet>(buffer, packet, stride); 1229 1230 for (int i = 0; i < PacketSize * 20; ++i) { 1231 if ((i % stride) == 0 && i < stride * PacketSize) { 1232 VERIFY(test::isApproxAbs(buffer[i], data1[i / stride], refvalue) && "pscatter"); 1233 } else { 1234 VERIFY(test::isApproxAbs(buffer[i], Scalar(0), refvalue) && "pscatter"); 1235 } 1236 } 1237 1238 for (int i = 0; i < PacketSize * 7; ++i) { 1239 buffer[i] = internal::random<Scalar>() / RealScalar(PacketSize); 1240 } 1241 packet = internal::pgather<Scalar, Packet>(buffer, 7); 1242 internal::pstore(data1, packet); 1243 for (int i = 0; i < PacketSize; ++i) { 1244 VERIFY(test::isApproxAbs(data1[i], buffer[i * 7], refvalue) && "pgather"); 1245 } 1246 } 1247 1248 namespace Eigen { 1249 namespace test { 1250 1251 template <typename Scalar, typename PacketType> 1252 struct runall<Scalar, PacketType, false, false> { // i.e. float or double 1253 static void run() { 1254 packetmath<Scalar, PacketType>(); 1255 packetmath_scatter_gather<Scalar, PacketType>(); 1256 packetmath_notcomplex<Scalar, PacketType>(); 1257 packetmath_real<Scalar, PacketType>(); 1258 } 1259 }; 1260 1261 template <typename Scalar, typename PacketType> 1262 struct runall<Scalar, PacketType, false, true> { // i.e. int 1263 static void run() { 1264 packetmath<Scalar, PacketType>(); 1265 packetmath_scatter_gather<Scalar, PacketType>(); 1266 packetmath_notcomplex<Scalar, PacketType>(); 1267 } 1268 }; 1269 1270 template <typename Scalar, typename PacketType> 1271 struct runall<Scalar, PacketType, true, false> { // i.e. complex 1272 static void run() { 1273 packetmath<Scalar, PacketType>(); 1274 packetmath_scatter_gather<Scalar, PacketType>(); 1275 packetmath_complex<Scalar, PacketType>(); 1276 } 1277 }; 1278 1279 } // namespace test 1280 } // namespace Eigen 1281 1282 EIGEN_DECLARE_TEST(packetmath) { 1283 g_first_pass = true; 1284 for (int i = 0; i < g_repeat; i++) { 1285 CALL_SUBTEST_1(test::runner<float>::run()); 1286 CALL_SUBTEST_2(test::runner<double>::run()); 1287 CALL_SUBTEST_3(test::runner<int8_t>::run()); 1288 CALL_SUBTEST_4(test::runner<uint8_t>::run()); 1289 CALL_SUBTEST_5(test::runner<int16_t>::run()); 1290 CALL_SUBTEST_6(test::runner<uint16_t>::run()); 1291 CALL_SUBTEST_7(test::runner<int32_t>::run()); 1292 CALL_SUBTEST_8(test::runner<uint32_t>::run()); 1293 CALL_SUBTEST_9(test::runner<int64_t>::run()); 1294 CALL_SUBTEST_10(test::runner<uint64_t>::run()); 1295 CALL_SUBTEST_11(test::runner<std::complex<float> >::run()); 1296 CALL_SUBTEST_12(test::runner<std::complex<double> >::run()); 1297 CALL_SUBTEST_13(test::runner<half>::run()); 1298 CALL_SUBTEST_14((packetmath<bool, internal::packet_traits<bool>::type>())); 1299 CALL_SUBTEST_15(test::runner<bfloat16>::run()); 1300 g_first_pass = false; 1301 } 1302 }