blasutil.cpp 6.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
  1. // This file is part of Eigen, a lightweight C++ template library
  2. // for linear algebra.
  3. //
  4. // Copyright (C) 2020 Everton Constantino <everton.constantino@ibm.com>
  5. //
  6. // This Source Code Form is subject to the terms of the Mozilla
  7. // Public License v. 2.0. If a copy of the MPL was not distributed
  8. // with this file, You can obtain one at http://mozilla.org/MPL/2.0/
  9. #include "main.h"
  10. // Disable "ignoring attributes on template argument"
  11. // for packet_traits<Packet*>
  12. // => The only workaround would be to wrap _m128 and the likes
  13. // within wrappers.
  14. #if EIGEN_GNUC_AT_LEAST(6,0)
  15. #pragma GCC diagnostic ignored "-Wignored-attributes"
  16. #endif
  17. #define GET(i,j) (StorageOrder == RowMajor ? (i)*stride + (j) : (i) + (j)*stride)
  18. #define SCATTER(i,j,k) (StorageOrder == RowMajor ? ((i)+(k))*stride + (j) : (i) + ((j)+(k))*stride)
  19. template<typename Scalar, typename Packet>
  20. void compare(const Packet& a, const Packet& b)
  21. {
  22. int pktsz = internal::packet_traits<Scalar>::size;
  23. Scalar *buffA = new Scalar[pktsz];
  24. Scalar *buffB = new Scalar[pktsz];
  25. internal::pstoreu<Scalar, Packet>(buffA, a);
  26. internal::pstoreu<Scalar, Packet>(buffB, b);
  27. for(int i = 0; i < pktsz; i++)
  28. {
  29. VERIFY_IS_EQUAL(buffA[i], buffB[i]);
  30. }
  31. delete[] buffA;
  32. delete[] buffB;
  33. }
  34. template<typename Scalar, int StorageOrder, int n>
  35. struct PacketBlockSet
  36. {
  37. typedef typename internal::packet_traits<Scalar>::type Packet;
  38. void setPacketBlock(internal::PacketBlock<Packet,n>& block, Scalar value)
  39. {
  40. for(int idx = 0; idx < n; idx++)
  41. {
  42. block.packet[idx] = internal::pset1<Packet>(value);
  43. }
  44. }
  45. void comparePacketBlock(Scalar *data, int i, int j, int stride, internal::PacketBlock<Packet, n>& block)
  46. {
  47. for(int idx = 0; idx < n; idx++)
  48. {
  49. Packet line = internal::ploadu<Packet>(data + SCATTER(i,j,idx));
  50. compare<Scalar, Packet>(block.packet[idx], line);
  51. }
  52. }
  53. };
  54. template<typename Scalar, int StorageOrder, int BlockSize>
  55. void run_bdmp_spec_1()
  56. {
  57. typedef internal::blas_data_mapper<Scalar, int, StorageOrder> BlasDataMapper;
  58. int packetSize = internal::packet_traits<Scalar>::size;
  59. int minSize = std::max<int>(packetSize, BlockSize);
  60. typedef typename internal::packet_traits<Scalar>::type Packet;
  61. int szm = internal::random<int>(minSize,500), szn = internal::random<int>(minSize,500);
  62. int stride = StorageOrder == RowMajor ? szn : szm;
  63. Scalar *d = new Scalar[szn*szm];
  64. // Initializing with random entries
  65. for(int i = 0; i < szm*szn; i++)
  66. {
  67. d[i] = internal::random<Scalar>(static_cast<Scalar>(3), static_cast<Scalar>(10));
  68. }
  69. BlasDataMapper bdm(d, stride);
  70. // Testing operator()
  71. for(int i = 0; i < szm; i++)
  72. {
  73. for(int j = 0; j < szn; j++)
  74. {
  75. VERIFY_IS_EQUAL(d[GET(i,j)], bdm(i,j));
  76. }
  77. }
  78. // Testing getSubMapper and getLinearMapper
  79. int i0 = internal::random<int>(0,szm-2);
  80. int j0 = internal::random<int>(0,szn-2);
  81. for(int i = i0; i < szm; i++)
  82. {
  83. for(int j = j0; j < szn; j++)
  84. {
  85. const BlasDataMapper& bdmSM = bdm.getSubMapper(i0,j0);
  86. const internal::BlasLinearMapper<Scalar, int, 0>& bdmLM = bdm.getLinearMapper(i0,j0);
  87. Scalar v = bdmSM(i - i0, j - j0);
  88. Scalar vd = d[GET(i,j)];
  89. VERIFY_IS_EQUAL(vd, v);
  90. VERIFY_IS_EQUAL(vd, bdmLM(GET(i-i0, j-j0)));
  91. }
  92. }
  93. // Testing loadPacket
  94. for(int i = 0; i < szm - minSize; i++)
  95. {
  96. for(int j = 0; j < szn - minSize; j++)
  97. {
  98. Packet pktBDM = bdm.template loadPacket<Packet>(i,j);
  99. Packet pktD = internal::ploadu<Packet>(d + GET(i,j));
  100. compare<Scalar, Packet>(pktBDM, pktD);
  101. }
  102. }
  103. // Testing gatherPacket
  104. Scalar *buff = new Scalar[packetSize];
  105. for(int i = 0; i < szm - minSize; i++)
  106. {
  107. for(int j = 0; j < szn - minSize; j++)
  108. {
  109. Packet p = bdm.template gatherPacket<Packet>(i,j);
  110. internal::pstoreu<Scalar, Packet>(buff, p);
  111. for(int k = 0; k < packetSize; k++)
  112. {
  113. VERIFY_IS_EQUAL(d[SCATTER(i,j,k)], buff[k]);
  114. }
  115. }
  116. }
  117. delete[] buff;
  118. // Testing scatterPacket
  119. for(int i = 0; i < szm - minSize; i++)
  120. {
  121. for(int j = 0; j < szn - minSize; j++)
  122. {
  123. Packet p = internal::pset1<Packet>(static_cast<Scalar>(1));
  124. bdm.template scatterPacket<Packet>(i,j,p);
  125. for(int k = 0; k < packetSize; k++)
  126. {
  127. VERIFY_IS_EQUAL(d[SCATTER(i,j,k)], static_cast<Scalar>(1));
  128. }
  129. }
  130. }
  131. //Testing storePacketBlock
  132. internal::PacketBlock<Packet, BlockSize> block;
  133. PacketBlockSet<Scalar, StorageOrder, BlockSize> pbs;
  134. pbs.setPacketBlock(block, static_cast<Scalar>(2));
  135. for(int i = 0; i < szm - minSize; i++)
  136. {
  137. for(int j = 0; j < szn - minSize; j++)
  138. {
  139. bdm.template storePacketBlock<Packet, BlockSize>(i, j, block);
  140. pbs.comparePacketBlock(d, i, j, stride, block);
  141. }
  142. }
  143. delete[] d;
  144. }
  145. template<typename Scalar>
  146. void run_test()
  147. {
  148. run_bdmp_spec_1<Scalar, RowMajor, 1>();
  149. run_bdmp_spec_1<Scalar, ColMajor, 1>();
  150. run_bdmp_spec_1<Scalar, RowMajor, 2>();
  151. run_bdmp_spec_1<Scalar, ColMajor, 2>();
  152. run_bdmp_spec_1<Scalar, RowMajor, 4>();
  153. run_bdmp_spec_1<Scalar, ColMajor, 4>();
  154. run_bdmp_spec_1<Scalar, RowMajor, 8>();
  155. run_bdmp_spec_1<Scalar, ColMajor, 8>();
  156. run_bdmp_spec_1<Scalar, RowMajor, 16>();
  157. run_bdmp_spec_1<Scalar, ColMajor, 16>();
  158. }
  159. EIGEN_DECLARE_TEST(blasutil)
  160. {
  161. for(int i = 0; i < g_repeat; i++)
  162. {
  163. CALL_SUBTEST_1(run_test<numext::int8_t>());
  164. CALL_SUBTEST_2(run_test<numext::int16_t>());
  165. CALL_SUBTEST_3(run_test<numext::int32_t>());
  166. // TODO: Replace this by a call to numext::int64_t as soon as we have a way to
  167. // detect the typedef for int64_t on all platforms
  168. #if EIGEN_HAS_CXX11
  169. CALL_SUBTEST_4(run_test<signed long long>());
  170. #else
  171. CALL_SUBTEST_4(run_test<signed long>());
  172. #endif
  173. CALL_SUBTEST_5(run_test<float_t>());
  174. CALL_SUBTEST_6(run_test<double_t>());
  175. CALL_SUBTEST_7(run_test<std::complex<float> >());
  176. CALL_SUBTEST_8(run_test<std::complex<double> >());
  177. }
  178. }