1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2014 Benoit Steiner <benoit.steiner.goog (at) gmail.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 10 #include "main.h" 11 12 #include <Eigen/CXX11/Tensor> 13 14 using Eigen::DefaultDevice; 15 using Eigen::Tensor; 16 17 typedef Tensor<float, 1>::DimensionPair DimPair; 18 19 template<int DataLayout> 20 static void test_evals() 21 { 22 Tensor<float, 2, DataLayout> mat1(2, 3); 23 Tensor<float, 2, DataLayout> mat2(2, 3); 24 Tensor<float, 2, DataLayout> mat3(3, 2); 25 26 mat1.setRandom(); 27 mat2.setRandom(); 28 mat3.setRandom(); 29 30 Tensor<float, 2, DataLayout> mat4(3,3); 31 mat4.setZero(); 32 Eigen::array<DimPair, 1> dims3 = {{DimPair(0, 0)}}; 33 typedef TensorEvaluator<decltype(mat1.contract(mat2, dims3)), DefaultDevice> Evaluator; 34 Evaluator eval(mat1.contract(mat2, dims3), DefaultDevice()); 35 eval.evalTo(mat4.data()); 36 EIGEN_STATIC_ASSERT(Evaluator::NumDims==2ul, YOU_MADE_A_PROGRAMMING_MISTAKE); 37 VERIFY_IS_EQUAL(eval.dimensions()[0], 3); 38 VERIFY_IS_EQUAL(eval.dimensions()[1], 3); 39 40 VERIFY_IS_APPROX(mat4(0,0), mat1(0,0)*mat2(0,0) + mat1(1,0)*mat2(1,0)); 41 VERIFY_IS_APPROX(mat4(0,1), mat1(0,0)*mat2(0,1) + mat1(1,0)*mat2(1,1)); 42 VERIFY_IS_APPROX(mat4(0,2), mat1(0,0)*mat2(0,2) + mat1(1,0)*mat2(1,2)); 43 VERIFY_IS_APPROX(mat4(1,0), mat1(0,1)*mat2(0,0) + mat1(1,1)*mat2(1,0)); 44 VERIFY_IS_APPROX(mat4(1,1), mat1(0,1)*mat2(0,1) + mat1(1,1)*mat2(1,1)); 45 VERIFY_IS_APPROX(mat4(1,2), mat1(0,1)*mat2(0,2) + mat1(1,1)*mat2(1,2)); 46 VERIFY_IS_APPROX(mat4(2,0), mat1(0,2)*mat2(0,0) + mat1(1,2)*mat2(1,0)); 47 VERIFY_IS_APPROX(mat4(2,1), mat1(0,2)*mat2(0,1) + mat1(1,2)*mat2(1,1)); 48 VERIFY_IS_APPROX(mat4(2,2), mat1(0,2)*mat2(0,2) + mat1(1,2)*mat2(1,2)); 49 50 Tensor<float, 2, DataLayout> mat5(2,2); 51 mat5.setZero(); 52 Eigen::array<DimPair, 1> dims4 = {{DimPair(1, 1)}}; 53 typedef TensorEvaluator<decltype(mat1.contract(mat2, dims4)), DefaultDevice> Evaluator2; 54 Evaluator2 eval2(mat1.contract(mat2, dims4), DefaultDevice()); 55 eval2.evalTo(mat5.data()); 56 EIGEN_STATIC_ASSERT(Evaluator2::NumDims==2ul, YOU_MADE_A_PROGRAMMING_MISTAKE); 57 VERIFY_IS_EQUAL(eval2.dimensions()[0], 2); 58 VERIFY_IS_EQUAL(eval2.dimensions()[1], 2); 59 60 VERIFY_IS_APPROX(mat5(0,0), mat1(0,0)*mat2(0,0) + mat1(0,1)*mat2(0,1) + mat1(0,2)*mat2(0,2)); 61 VERIFY_IS_APPROX(mat5(0,1), mat1(0,0)*mat2(1,0) + mat1(0,1)*mat2(1,1) + mat1(0,2)*mat2(1,2)); 62 VERIFY_IS_APPROX(mat5(1,0), mat1(1,0)*mat2(0,0) + mat1(1,1)*mat2(0,1) + mat1(1,2)*mat2(0,2)); 63 VERIFY_IS_APPROX(mat5(1,1), mat1(1,0)*mat2(1,0) + mat1(1,1)*mat2(1,1) + mat1(1,2)*mat2(1,2)); 64 65 Tensor<float, 2, DataLayout> mat6(2,2); 66 mat6.setZero(); 67 Eigen::array<DimPair, 1> dims6 = {{DimPair(1, 0)}}; 68 typedef TensorEvaluator<decltype(mat1.contract(mat3, dims6)), DefaultDevice> Evaluator3; 69 Evaluator3 eval3(mat1.contract(mat3, dims6), DefaultDevice()); 70 eval3.evalTo(mat6.data()); 71 EIGEN_STATIC_ASSERT(Evaluator3::NumDims==2ul, YOU_MADE_A_PROGRAMMING_MISTAKE); 72 VERIFY_IS_EQUAL(eval3.dimensions()[0], 2); 73 VERIFY_IS_EQUAL(eval3.dimensions()[1], 2); 74 75 VERIFY_IS_APPROX(mat6(0,0), mat1(0,0)*mat3(0,0) + mat1(0,1)*mat3(1,0) + mat1(0,2)*mat3(2,0)); 76 VERIFY_IS_APPROX(mat6(0,1), mat1(0,0)*mat3(0,1) + mat1(0,1)*mat3(1,1) + mat1(0,2)*mat3(2,1)); 77 VERIFY_IS_APPROX(mat6(1,0), mat1(1,0)*mat3(0,0) + mat1(1,1)*mat3(1,0) + mat1(1,2)*mat3(2,0)); 78 VERIFY_IS_APPROX(mat6(1,1), mat1(1,0)*mat3(0,1) + mat1(1,1)*mat3(1,1) + mat1(1,2)*mat3(2,1)); 79 } 80 81 template<int DataLayout> 82 static void test_scalar() 83 { 84 Tensor<float, 1, DataLayout> vec1({6}); 85 Tensor<float, 1, DataLayout> vec2({6}); 86 87 vec1.setRandom(); 88 vec2.setRandom(); 89 90 Eigen::array<DimPair, 1> dims = {{DimPair(0, 0)}}; 91 Tensor<float, 0, DataLayout> scalar = vec1.contract(vec2, dims); 92 93 float expected = 0.0f; 94 for (int i = 0; i < 6; ++i) { 95 expected += vec1(i) * vec2(i); 96 } 97 VERIFY_IS_APPROX(scalar(), expected); 98 } 99 100 template<int DataLayout> 101 static void test_multidims() 102 { 103 Tensor<float, 3, DataLayout> mat1(2, 2, 2); 104 Tensor<float, 4, DataLayout> mat2(2, 2, 2, 2); 105 106 mat1.setRandom(); 107 mat2.setRandom(); 108 109 Tensor<float, 3, DataLayout> mat3(2, 2, 2); 110 mat3.setZero(); 111 Eigen::array<DimPair, 2> dims = {{DimPair(1, 2), DimPair(2, 3)}}; 112 typedef TensorEvaluator<decltype(mat1.contract(mat2, dims)), DefaultDevice> Evaluator; 113 Evaluator eval(mat1.contract(mat2, dims), DefaultDevice()); 114 eval.evalTo(mat3.data()); 115 EIGEN_STATIC_ASSERT(Evaluator::NumDims==3ul, YOU_MADE_A_PROGRAMMING_MISTAKE); 116 VERIFY_IS_EQUAL(eval.dimensions()[0], 2); 117 VERIFY_IS_EQUAL(eval.dimensions()[1], 2); 118 VERIFY_IS_EQUAL(eval.dimensions()[2], 2); 119 120 VERIFY_IS_APPROX(mat3(0,0,0), mat1(0,0,0)*mat2(0,0,0,0) + mat1(0,1,0)*mat2(0,0,1,0) + 121 mat1(0,0,1)*mat2(0,0,0,1) + mat1(0,1,1)*mat2(0,0,1,1)); 122 VERIFY_IS_APPROX(mat3(0,0,1), mat1(0,0,0)*mat2(0,1,0,0) + mat1(0,1,0)*mat2(0,1,1,0) + 123 mat1(0,0,1)*mat2(0,1,0,1) + mat1(0,1,1)*mat2(0,1,1,1)); 124 VERIFY_IS_APPROX(mat3(0,1,0), mat1(0,0,0)*mat2(1,0,0,0) + mat1(0,1,0)*mat2(1,0,1,0) + 125 mat1(0,0,1)*mat2(1,0,0,1) + mat1(0,1,1)*mat2(1,0,1,1)); 126 VERIFY_IS_APPROX(mat3(0,1,1), mat1(0,0,0)*mat2(1,1,0,0) + mat1(0,1,0)*mat2(1,1,1,0) + 127 mat1(0,0,1)*mat2(1,1,0,1) + mat1(0,1,1)*mat2(1,1,1,1)); 128 VERIFY_IS_APPROX(mat3(1,0,0), mat1(1,0,0)*mat2(0,0,0,0) + mat1(1,1,0)*mat2(0,0,1,0) + 129 mat1(1,0,1)*mat2(0,0,0,1) + mat1(1,1,1)*mat2(0,0,1,1)); 130 VERIFY_IS_APPROX(mat3(1,0,1), mat1(1,0,0)*mat2(0,1,0,0) + mat1(1,1,0)*mat2(0,1,1,0) + 131 mat1(1,0,1)*mat2(0,1,0,1) + mat1(1,1,1)*mat2(0,1,1,1)); 132 VERIFY_IS_APPROX(mat3(1,1,0), mat1(1,0,0)*mat2(1,0,0,0) + mat1(1,1,0)*mat2(1,0,1,0) + 133 mat1(1,0,1)*mat2(1,0,0,1) + mat1(1,1,1)*mat2(1,0,1,1)); 134 VERIFY_IS_APPROX(mat3(1,1,1), mat1(1,0,0)*mat2(1,1,0,0) + mat1(1,1,0)*mat2(1,1,1,0) + 135 mat1(1,0,1)*mat2(1,1,0,1) + mat1(1,1,1)*mat2(1,1,1,1)); 136 137 Tensor<float, 2, DataLayout> mat4(2, 2); 138 Tensor<float, 3, DataLayout> mat5(2, 2, 2); 139 140 mat4.setRandom(); 141 mat5.setRandom(); 142 143 Tensor<float, 1, DataLayout> mat6(2); 144 mat6.setZero(); 145 Eigen::array<DimPair, 2> dims2({{DimPair(0, 1), DimPair(1, 0)}}); 146 typedef TensorEvaluator<decltype(mat4.contract(mat5, dims2)), DefaultDevice> Evaluator2; 147 Evaluator2 eval2(mat4.contract(mat5, dims2), DefaultDevice()); 148 eval2.evalTo(mat6.data()); 149 EIGEN_STATIC_ASSERT(Evaluator2::NumDims==1ul, YOU_MADE_A_PROGRAMMING_MISTAKE); 150 VERIFY_IS_EQUAL(eval2.dimensions()[0], 2); 151 152 VERIFY_IS_APPROX(mat6(0), mat4(0,0)*mat5(0,0,0) + mat4(1,0)*mat5(0,1,0) + 153 mat4(0,1)*mat5(1,0,0) + mat4(1,1)*mat5(1,1,0)); 154 VERIFY_IS_APPROX(mat6(1), mat4(0,0)*mat5(0,0,1) + mat4(1,0)*mat5(0,1,1) + 155 mat4(0,1)*mat5(1,0,1) + mat4(1,1)*mat5(1,1,1)); 156 } 157 158 template<int DataLayout> 159 static void test_holes() { 160 Tensor<float, 4, DataLayout> t1(2, 5, 7, 3); 161 Tensor<float, 5, DataLayout> t2(2, 7, 11, 13, 3); 162 t1.setRandom(); 163 t2.setRandom(); 164 165 Eigen::array<DimPair, 2> dims = {{DimPair(0, 0), DimPair(3, 4)}}; 166 Tensor<float, 5, DataLayout> result = t1.contract(t2, dims); 167 VERIFY_IS_EQUAL(result.dimension(0), 5); 168 VERIFY_IS_EQUAL(result.dimension(1), 7); 169 VERIFY_IS_EQUAL(result.dimension(2), 7); 170 VERIFY_IS_EQUAL(result.dimension(3), 11); 171 VERIFY_IS_EQUAL(result.dimension(4), 13); 172 173 for (int i = 0; i < 5; ++i) { 174 for (int j = 0; j < 5; ++j) { 175 for (int k = 0; k < 5; ++k) { 176 for (int l = 0; l < 5; ++l) { 177 for (int m = 0; m < 5; ++m) { 178 VERIFY_IS_APPROX(result(i, j, k, l, m), 179 t1(0, i, j, 0) * t2(0, k, l, m, 0) + 180 t1(1, i, j, 0) * t2(1, k, l, m, 0) + 181 t1(0, i, j, 1) * t2(0, k, l, m, 1) + 182 t1(1, i, j, 1) * t2(1, k, l, m, 1) + 183 t1(0, i, j, 2) * t2(0, k, l, m, 2) + 184 t1(1, i, j, 2) * t2(1, k, l, m, 2)); 185 } 186 } 187 } 188 } 189 } 190 } 191 192 template<int DataLayout> 193 static void test_full_redux() 194 { 195 Tensor<float, 2, DataLayout> t1(2, 2); 196 Tensor<float, 3, DataLayout> t2(2, 2, 2); 197 t1.setRandom(); 198 t2.setRandom(); 199 200 Eigen::array<DimPair, 2> dims = {{DimPair(0, 0), DimPair(1, 1)}}; 201 Tensor<float, 1, DataLayout> result = t1.contract(t2, dims); 202 VERIFY_IS_EQUAL(result.dimension(0), 2); 203 VERIFY_IS_APPROX(result(0), t1(0, 0) * t2(0, 0, 0) + t1(1, 0) * t2(1, 0, 0) 204 + t1(0, 1) * t2(0, 1, 0) + t1(1, 1) * t2(1, 1, 0)); 205 VERIFY_IS_APPROX(result(1), t1(0, 0) * t2(0, 0, 1) + t1(1, 0) * t2(1, 0, 1) 206 + t1(0, 1) * t2(0, 1, 1) + t1(1, 1) * t2(1, 1, 1)); 207 208 dims[0] = DimPair(1, 0); 209 dims[1] = DimPair(2, 1); 210 result = t2.contract(t1, dims); 211 VERIFY_IS_EQUAL(result.dimension(0), 2); 212 VERIFY_IS_APPROX(result(0), t1(0, 0) * t2(0, 0, 0) + t1(1, 0) * t2(0, 1, 0) 213 + t1(0, 1) * t2(0, 0, 1) + t1(1, 1) * t2(0, 1, 1)); 214 VERIFY_IS_APPROX(result(1), t1(0, 0) * t2(1, 0, 0) + t1(1, 0) * t2(1, 1, 0) 215 + t1(0, 1) * t2(1, 0, 1) + t1(1, 1) * t2(1, 1, 1)); 216 } 217 218 template<int DataLayout> 219 static void test_contraction_of_contraction() 220 { 221 Tensor<float, 2, DataLayout> t1(2, 2); 222 Tensor<float, 2, DataLayout> t2(2, 2); 223 Tensor<float, 2, DataLayout> t3(2, 2); 224 Tensor<float, 2, DataLayout> t4(2, 2); 225 t1.setRandom(); 226 t2.setRandom(); 227 t3.setRandom(); 228 t4.setRandom(); 229 230 Eigen::array<DimPair, 1> dims = {{DimPair(1, 0)}}; 231 auto contract1 = t1.contract(t2, dims); 232 auto diff = t3 - contract1; 233 auto contract2 = t1.contract(t4, dims); 234 Tensor<float, 2, DataLayout> result = contract2.contract(diff, dims); 235 236 VERIFY_IS_EQUAL(result.dimension(0), 2); 237 VERIFY_IS_EQUAL(result.dimension(1), 2); 238 239 Eigen::Map<Eigen::Matrix<float, Dynamic, Dynamic, DataLayout>> 240 m1(t1.data(), 2, 2), m2(t2.data(), 2, 2), m3(t3.data(), 2, 2), 241 m4(t4.data(), 2, 2); 242 Eigen::Matrix<float, Dynamic, Dynamic, DataLayout> 243 expected = (m1 * m4) * (m3 - m1 * m2); 244 245 VERIFY_IS_APPROX(result(0, 0), expected(0, 0)); 246 VERIFY_IS_APPROX(result(0, 1), expected(0, 1)); 247 VERIFY_IS_APPROX(result(1, 0), expected(1, 0)); 248 VERIFY_IS_APPROX(result(1, 1), expected(1, 1)); 249 } 250 251 template<int DataLayout> 252 static void test_expr() 253 { 254 Tensor<float, 2, DataLayout> mat1(2, 3); 255 Tensor<float, 2, DataLayout> mat2(3, 2); 256 mat1.setRandom(); 257 mat2.setRandom(); 258 259 Tensor<float, 2, DataLayout> mat3(2,2); 260 261 Eigen::array<DimPair, 1> dims = {{DimPair(1, 0)}}; 262 mat3 = mat1.contract(mat2, dims); 263 264 VERIFY_IS_APPROX(mat3(0,0), mat1(0,0)*mat2(0,0) + mat1(0,1)*mat2(1,0) + mat1(0,2)*mat2(2,0)); 265 VERIFY_IS_APPROX(mat3(0,1), mat1(0,0)*mat2(0,1) + mat1(0,1)*mat2(1,1) + mat1(0,2)*mat2(2,1)); 266 VERIFY_IS_APPROX(mat3(1,0), mat1(1,0)*mat2(0,0) + mat1(1,1)*mat2(1,0) + mat1(1,2)*mat2(2,0)); 267 VERIFY_IS_APPROX(mat3(1,1), mat1(1,0)*mat2(0,1) + mat1(1,1)*mat2(1,1) + mat1(1,2)*mat2(2,1)); 268 } 269 270 template<int DataLayout> 271 static void test_out_of_order_contraction() 272 { 273 Tensor<float, 3, DataLayout> mat1(2, 2, 2); 274 Tensor<float, 3, DataLayout> mat2(2, 2, 2); 275 276 mat1.setRandom(); 277 mat2.setRandom(); 278 279 Tensor<float, 2, DataLayout> mat3(2, 2); 280 281 Eigen::array<DimPair, 2> dims = {{DimPair(2, 0), DimPair(0, 2)}}; 282 mat3 = mat1.contract(mat2, dims); 283 284 VERIFY_IS_APPROX(mat3(0, 0), 285 mat1(0,0,0)*mat2(0,0,0) + mat1(1,0,0)*mat2(0,0,1) + 286 mat1(0,0,1)*mat2(1,0,0) + mat1(1,0,1)*mat2(1,0,1)); 287 VERIFY_IS_APPROX(mat3(1, 0), 288 mat1(0,1,0)*mat2(0,0,0) + mat1(1,1,0)*mat2(0,0,1) + 289 mat1(0,1,1)*mat2(1,0,0) + mat1(1,1,1)*mat2(1,0,1)); 290 VERIFY_IS_APPROX(mat3(0, 1), 291 mat1(0,0,0)*mat2(0,1,0) + mat1(1,0,0)*mat2(0,1,1) + 292 mat1(0,0,1)*mat2(1,1,0) + mat1(1,0,1)*mat2(1,1,1)); 293 VERIFY_IS_APPROX(mat3(1, 1), 294 mat1(0,1,0)*mat2(0,1,0) + mat1(1,1,0)*mat2(0,1,1) + 295 mat1(0,1,1)*mat2(1,1,0) + mat1(1,1,1)*mat2(1,1,1)); 296 297 Eigen::array<DimPair, 2> dims2 = {{DimPair(0, 2), DimPair(2, 0)}}; 298 mat3 = mat1.contract(mat2, dims2); 299 300 VERIFY_IS_APPROX(mat3(0, 0), 301 mat1(0,0,0)*mat2(0,0,0) + mat1(1,0,0)*mat2(0,0,1) + 302 mat1(0,0,1)*mat2(1,0,0) + mat1(1,0,1)*mat2(1,0,1)); 303 VERIFY_IS_APPROX(mat3(1, 0), 304 mat1(0,1,0)*mat2(0,0,0) + mat1(1,1,0)*mat2(0,0,1) + 305 mat1(0,1,1)*mat2(1,0,0) + mat1(1,1,1)*mat2(1,0,1)); 306 VERIFY_IS_APPROX(mat3(0, 1), 307 mat1(0,0,0)*mat2(0,1,0) + mat1(1,0,0)*mat2(0,1,1) + 308 mat1(0,0,1)*mat2(1,1,0) + mat1(1,0,1)*mat2(1,1,1)); 309 VERIFY_IS_APPROX(mat3(1, 1), 310 mat1(0,1,0)*mat2(0,1,0) + mat1(1,1,0)*mat2(0,1,1) + 311 mat1(0,1,1)*mat2(1,1,0) + mat1(1,1,1)*mat2(1,1,1)); 312 313 } 314 315 template<int DataLayout> 316 static void test_consistency() 317 { 318 // this does something like testing (A*B)^T = (B^T * A^T) 319 320 Tensor<float, 3, DataLayout> mat1(4, 3, 5); 321 Tensor<float, 5, DataLayout> mat2(3, 2, 1, 5, 4); 322 mat1.setRandom(); 323 mat2.setRandom(); 324 325 Tensor<float, 4, DataLayout> mat3(5, 2, 1, 5); 326 Tensor<float, 4, DataLayout> mat4(2, 1, 5, 5); 327 328 // contract on dimensions of size 4 and 3 329 Eigen::array<DimPair, 2> dims1 = {{DimPair(0, 4), DimPair(1, 0)}}; 330 Eigen::array<DimPair, 2> dims2 = {{DimPair(4, 0), DimPair(0, 1)}}; 331 332 mat3 = mat1.contract(mat2, dims1); 333 mat4 = mat2.contract(mat1, dims2); 334 335 // check that these are equal except for ordering of dimensions 336 if (DataLayout == ColMajor) { 337 for (size_t i = 0; i < 5; i++) { 338 for (size_t j = 0; j < 10; j++) { 339 VERIFY_IS_APPROX(mat3.data()[i + 5 * j], mat4.data()[j + 10 * i]); 340 } 341 } 342 } else { 343 // Row major 344 for (size_t i = 0; i < 5; i++) { 345 for (size_t j = 0; j < 10; j++) { 346 VERIFY_IS_APPROX(mat3.data()[10 * i + j], mat4.data()[i + 5 * j]); 347 } 348 } 349 } 350 } 351 352 template<int DataLayout> 353 static void test_large_contraction() 354 { 355 Tensor<float, 4, DataLayout> t_left(30, 50, 8, 31); 356 Tensor<float, 5, DataLayout> t_right(8, 31, 7, 20, 10); 357 Tensor<float, 5, DataLayout> t_result(30, 50, 7, 20, 10); 358 359 t_left.setRandom(); 360 t_right.setRandom(); 361 362 // Add a little offset so that the results won't be close to zero. 363 t_left += t_left.constant(1.0f); 364 t_right += t_right.constant(1.0f); 365 366 typedef Map<Eigen::Matrix<float, Dynamic, Dynamic, DataLayout>> MapXf; 367 MapXf m_left(t_left.data(), 1500, 248); 368 MapXf m_right(t_right.data(), 248, 1400); 369 Eigen::Matrix<float, Dynamic, Dynamic, DataLayout> m_result(1500, 1400); 370 371 // this contraction should be equivalent to a single matrix multiplication 372 Eigen::array<DimPair, 2> dims = {{DimPair(2, 0), DimPair(3, 1)}}; 373 374 // compute results by separate methods 375 t_result = t_left.contract(t_right, dims); 376 m_result = m_left * m_right; 377 378 for (int i = 0; i < t_result.dimensions().TotalSize(); i++) { 379 VERIFY(&t_result.data()[i] != &m_result.data()[i]); 380 VERIFY_IS_APPROX(t_result.data()[i], m_result.data()[i]); 381 } 382 } 383 384 template<int DataLayout> 385 static void test_matrix_vector() 386 { 387 Tensor<float, 2, DataLayout> t_left(30, 50); 388 Tensor<float, 1, DataLayout> t_right(50); 389 Tensor<float, 1, DataLayout> t_result(30); 390 391 t_left.setRandom(); 392 t_right.setRandom(); 393 394 typedef Map<Eigen::Matrix<float, Dynamic, Dynamic, DataLayout>> MapXf; 395 MapXf m_left(t_left.data(), 30, 50); 396 MapXf m_right(t_right.data(), 50, 1); 397 Eigen::Matrix<float, Dynamic, Dynamic, DataLayout> m_result(30, 1); 398 399 // this contraction should be equivalent to a single matrix multiplication 400 Eigen::array<DimPair, 1> dims{{DimPair(1, 0)}}; 401 402 // compute results by separate methods 403 t_result = t_left.contract(t_right, dims); 404 m_result = m_left * m_right; 405 406 for (int i = 0; i < t_result.dimensions().TotalSize(); i++) { 407 VERIFY(internal::isApprox(t_result(i), m_result(i, 0), 1)); 408 } 409 } 410 411 412 template<int DataLayout> 413 static void test_tensor_vector() 414 { 415 Tensor<float, 3, DataLayout> t_left(7, 13, 17); 416 Tensor<float, 2, DataLayout> t_right(1, 7); 417 418 t_left.setRandom(); 419 t_right.setRandom(); 420 421 typedef typename Tensor<float, 1, DataLayout>::DimensionPair DimensionPair; 422 Eigen::array<DimensionPair, 1> dim_pair01{{{0, 1}}}; 423 Tensor<float, 3, DataLayout> t_result = t_left.contract(t_right, dim_pair01); 424 425 typedef Map<Eigen::Matrix<float, Dynamic, Dynamic, DataLayout>> MapXf; 426 MapXf m_left(t_left.data(), 7, 13*17); 427 MapXf m_right(t_right.data(), 1, 7); 428 Eigen::Matrix<float, Dynamic, Dynamic, DataLayout> m_result = m_left.transpose() * m_right.transpose(); 429 430 for (int i = 0; i < t_result.dimensions().TotalSize(); i++) { 431 VERIFY(internal::isApprox(t_result(i), m_result(i, 0), 1)); 432 } 433 } 434 435 436 template<int DataLayout> 437 static void test_small_blocking_factors() 438 { 439 Tensor<float, 4, DataLayout> t_left(30, 5, 3, 31); 440 Tensor<float, 5, DataLayout> t_right(3, 31, 7, 20, 1); 441 t_left.setRandom(); 442 t_right.setRandom(); 443 444 // Add a little offset so that the results won't be close to zero. 445 t_left += t_left.constant(1.0f); 446 t_right += t_right.constant(1.0f); 447 448 // Force the cache sizes, which results in smaller blocking factors. 449 Eigen::setCpuCacheSizes(896, 1920, 2944); 450 451 // this contraction should be equivalent to a single matrix multiplication 452 Eigen::array<DimPair, 2> dims = {{DimPair(2, 0), DimPair(3, 1)}}; 453 Tensor<float, 5, DataLayout> t_result; 454 t_result = t_left.contract(t_right, dims); 455 456 // compute result using a simple eigen matrix product 457 Map<Eigen::Matrix<float, Dynamic, Dynamic, DataLayout>> m_left(t_left.data(), 150, 93); 458 Map<Eigen::Matrix<float, Dynamic, Dynamic, DataLayout>> m_right(t_right.data(), 93, 140); 459 Eigen::Matrix<float, Dynamic, Dynamic, DataLayout> m_result = m_left * m_right; 460 461 for (int i = 0; i < t_result.dimensions().TotalSize(); i++) { 462 VERIFY_IS_APPROX(t_result.data()[i], m_result.data()[i]); 463 } 464 } 465 466 template<int DataLayout> 467 static void test_tensor_product() 468 { 469 Tensor<float, 2, DataLayout> mat1(2, 3); 470 Tensor<float, 2, DataLayout> mat2(4, 1); 471 mat1.setRandom(); 472 mat2.setRandom(); 473 474 Tensor<float, 4, DataLayout> result = mat1.contract(mat2, Eigen::array<DimPair, 0>{{}}); 475 476 VERIFY_IS_EQUAL(result.dimension(0), 2); 477 VERIFY_IS_EQUAL(result.dimension(1), 3); 478 VERIFY_IS_EQUAL(result.dimension(2), 4); 479 VERIFY_IS_EQUAL(result.dimension(3), 1); 480 for (int i = 0; i < result.dimension(0); ++i) { 481 for (int j = 0; j < result.dimension(1); ++j) { 482 for (int k = 0; k < result.dimension(2); ++k) { 483 for (int l = 0; l < result.dimension(3); ++l) { 484 VERIFY_IS_APPROX(result(i, j, k, l), mat1(i, j) * mat2(k, l) ); 485 } 486 } 487 } 488 } 489 } 490 491 492 template<int DataLayout> 493 static void test_const_inputs() 494 { 495 Tensor<float, 2, DataLayout> in1(2, 3); 496 Tensor<float, 2, DataLayout> in2(3, 2); 497 in1.setRandom(); 498 in2.setRandom(); 499 500 TensorMap<Tensor<const float, 2, DataLayout> > mat1(in1.data(), 2, 3); 501 TensorMap<Tensor<const float, 2, DataLayout> > mat2(in2.data(), 3, 2); 502 Tensor<float, 2, DataLayout> mat3(2,2); 503 504 Eigen::array<DimPair, 1> dims = {{DimPair(1, 0)}}; 505 mat3 = mat1.contract(mat2, dims); 506 507 VERIFY_IS_APPROX(mat3(0,0), mat1(0,0)*mat2(0,0) + mat1(0,1)*mat2(1,0) + mat1(0,2)*mat2(2,0)); 508 VERIFY_IS_APPROX(mat3(0,1), mat1(0,0)*mat2(0,1) + mat1(0,1)*mat2(1,1) + mat1(0,2)*mat2(2,1)); 509 VERIFY_IS_APPROX(mat3(1,0), mat1(1,0)*mat2(0,0) + mat1(1,1)*mat2(1,0) + mat1(1,2)*mat2(2,0)); 510 VERIFY_IS_APPROX(mat3(1,1), mat1(1,0)*mat2(0,1) + mat1(1,1)*mat2(1,1) + mat1(1,2)*mat2(2,1)); 511 } 512 513 void test_cxx11_tensor_contraction() 514 { 515 CALL_SUBTEST(test_evals<ColMajor>()); 516 CALL_SUBTEST(test_evals<RowMajor>()); 517 CALL_SUBTEST(test_scalar<ColMajor>()); 518 CALL_SUBTEST(test_scalar<RowMajor>()); 519 CALL_SUBTEST(test_multidims<ColMajor>()); 520 CALL_SUBTEST(test_multidims<RowMajor>()); 521 CALL_SUBTEST(test_holes<ColMajor>()); 522 CALL_SUBTEST(test_holes<RowMajor>()); 523 CALL_SUBTEST(test_full_redux<ColMajor>()); 524 CALL_SUBTEST(test_full_redux<RowMajor>()); 525 CALL_SUBTEST(test_contraction_of_contraction<ColMajor>()); 526 CALL_SUBTEST(test_contraction_of_contraction<RowMajor>()); 527 CALL_SUBTEST(test_expr<ColMajor>()); 528 CALL_SUBTEST(test_expr<RowMajor>()); 529 CALL_SUBTEST(test_out_of_order_contraction<ColMajor>()); 530 CALL_SUBTEST(test_out_of_order_contraction<RowMajor>()); 531 CALL_SUBTEST(test_consistency<ColMajor>()); 532 CALL_SUBTEST(test_consistency<RowMajor>()); 533 CALL_SUBTEST(test_large_contraction<ColMajor>()); 534 CALL_SUBTEST(test_large_contraction<RowMajor>()); 535 CALL_SUBTEST(test_matrix_vector<ColMajor>()); 536 CALL_SUBTEST(test_matrix_vector<RowMajor>()); 537 CALL_SUBTEST(test_tensor_vector<ColMajor>()); 538 CALL_SUBTEST(test_tensor_vector<RowMajor>()); 539 CALL_SUBTEST(test_small_blocking_factors<ColMajor>()); 540 CALL_SUBTEST(test_small_blocking_factors<RowMajor>()); 541 CALL_SUBTEST(test_tensor_product<ColMajor>()); 542 CALL_SUBTEST(test_tensor_product<RowMajor>()); 543 CALL_SUBTEST(test_const_inputs<ColMajor>()); 544 CALL_SUBTEST(test_const_inputs<RowMajor>()); 545 } 546
