cxx11_tensor_contraction.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) 2014 Benoit Steiner <benoit.steiner.goog@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  Eigen::array<DimPair, 0> dims;
475  Tensor<float, 4, DataLayout> result = mat1.contract(mat2, dims);
476 
477  VERIFY_IS_EQUAL(result.dimension(0), 2);
478  VERIFY_IS_EQUAL(result.dimension(1), 3);
479  VERIFY_IS_EQUAL(result.dimension(2), 4);
480  VERIFY_IS_EQUAL(result.dimension(3), 1);
481  for (int i = 0; i < result.dimension(0); ++i) {
482  for (int j = 0; j < result.dimension(1); ++j) {
483  for (int k = 0; k < result.dimension(2); ++k) {
484  for (int l = 0; l < result.dimension(3); ++l) {
485  VERIFY_IS_APPROX(result(i, j, k, l), mat1(i, j) * mat2(k, l) );
486  }
487  }
488  }
489  }
490 }
491 
492 
493 template<int DataLayout>
494 static void test_const_inputs()
495 {
496  Tensor<float, 2, DataLayout> in1(2, 3);
497  Tensor<float, 2, DataLayout> in2(3, 2);
498  in1.setRandom();
499  in2.setRandom();
500 
501  TensorMap<Tensor<const float, 2, DataLayout> > mat1(in1.data(), 2, 3);
502  TensorMap<Tensor<const float, 2, DataLayout> > mat2(in2.data(), 3, 2);
503  Tensor<float, 2, DataLayout> mat3(2,2);
504 
505  Eigen::array<DimPair, 1> dims = {{DimPair(1, 0)}};
506  mat3 = mat1.contract(mat2, dims);
507 
508  VERIFY_IS_APPROX(mat3(0,0), mat1(0,0)*mat2(0,0) + mat1(0,1)*mat2(1,0) + mat1(0,2)*mat2(2,0));
509  VERIFY_IS_APPROX(mat3(0,1), mat1(0,0)*mat2(0,1) + mat1(0,1)*mat2(1,1) + mat1(0,2)*mat2(2,1));
510  VERIFY_IS_APPROX(mat3(1,0), mat1(1,0)*mat2(0,0) + mat1(1,1)*mat2(1,0) + mat1(1,2)*mat2(2,0));
511  VERIFY_IS_APPROX(mat3(1,1), mat1(1,0)*mat2(0,1) + mat1(1,1)*mat2(1,1) + mat1(1,2)*mat2(2,1));
512 }
513 
514 // Apply Sqrt to all output elements.
515 struct SqrtOutputKernel {
516  template <typename Index, typename Scalar>
517  EIGEN_ALWAYS_INLINE void operator()(
518  const internal::blas_data_mapper<Scalar, Index, ColMajor>& output_mapper,
519  const TensorContractionParams&, Index, Index, Index num_rows,
520  Index num_cols) const {
521  for (int i = 0; i < num_rows; ++i) {
522  for (int j = 0; j < num_cols; ++j) {
523  output_mapper(i, j) = std::sqrt(output_mapper(i, j));
524  }
525  }
526  }
527 };
528 
529 template <int DataLayout>
530 static void test_large_contraction_with_output_kernel() {
531  Tensor<float, 4, DataLayout> t_left(30, 50, 8, 31);
532  Tensor<float, 5, DataLayout> t_right(8, 31, 7, 20, 10);
533  Tensor<float, 5, DataLayout> t_result(30, 50, 7, 20, 10);
534 
535  t_left.setRandom();
536  t_right.setRandom();
537  // Put trash in mat4 to verify contraction clears output memory.
538  t_result.setRandom();
539 
540  // Add a little offset so that the results won't be close to zero.
541  t_left += t_left.constant(1.0f);
542  t_right += t_right.constant(1.0f);
543 
544  typedef Map<Eigen::Matrix<float, Dynamic, Dynamic, DataLayout>> MapXf;
545  MapXf m_left(t_left.data(), 1500, 248);
546  MapXf m_right(t_right.data(), 248, 1400);
547  Eigen::Matrix<float, Dynamic, Dynamic, DataLayout> m_result(1500, 1400);
548 
549  // this contraction should be equivalent to a single matrix multiplication
550  Eigen::array<DimPair, 2> dims({{DimPair(2, 0), DimPair(3, 1)}});
551 
552  // compute results by separate methods
553  t_result = t_left.contract(t_right, dims, SqrtOutputKernel());
554 
555  m_result = m_left * m_right;
556 
557  for (std::ptrdiff_t i = 0; i < t_result.dimensions().TotalSize(); i++) {
558  VERIFY(&t_result.data()[i] != &m_result.data()[i]);
559  VERIFY_IS_APPROX(t_result.data()[i], std::sqrt(m_result.data()[i]));
560  }
561 }
562 
563 EIGEN_DECLARE_TEST(cxx11_tensor_contraction)
564 {
565  CALL_SUBTEST_1(test_evals<ColMajor>());
566  CALL_SUBTEST_1(test_evals<RowMajor>());
567  CALL_SUBTEST_1(test_scalar<ColMajor>());
568  CALL_SUBTEST_1(test_scalar<RowMajor>());
569  CALL_SUBTEST_2(test_multidims<ColMajor>());
570  CALL_SUBTEST_2(test_multidims<RowMajor>());
571  CALL_SUBTEST_2(test_holes<ColMajor>());
572  CALL_SUBTEST_2(test_holes<RowMajor>());
573  CALL_SUBTEST_3(test_full_redux<ColMajor>());
574  CALL_SUBTEST_3(test_full_redux<RowMajor>());
575  CALL_SUBTEST_3(test_contraction_of_contraction<ColMajor>());
576  CALL_SUBTEST_3(test_contraction_of_contraction<RowMajor>());
577  CALL_SUBTEST_4(test_expr<ColMajor>());
578  CALL_SUBTEST_4(test_expr<RowMajor>());
579  CALL_SUBTEST_4(test_out_of_order_contraction<ColMajor>());
580  CALL_SUBTEST_4(test_out_of_order_contraction<RowMajor>());
581  CALL_SUBTEST_5(test_consistency<ColMajor>());
582  CALL_SUBTEST_5(test_consistency<RowMajor>());
583  CALL_SUBTEST_5(test_large_contraction<ColMajor>());
584  CALL_SUBTEST_5(test_large_contraction<RowMajor>());
585  CALL_SUBTEST_6(test_matrix_vector<ColMajor>());
586  CALL_SUBTEST_6(test_matrix_vector<RowMajor>());
587  CALL_SUBTEST_6(test_tensor_vector<ColMajor>());
588  CALL_SUBTEST_6(test_tensor_vector<RowMajor>());
589  CALL_SUBTEST_7(test_small_blocking_factors<ColMajor>());
590  CALL_SUBTEST_7(test_small_blocking_factors<RowMajor>());
591  CALL_SUBTEST_7(test_tensor_product<ColMajor>());
592  CALL_SUBTEST_7(test_tensor_product<RowMajor>());
593  CALL_SUBTEST_8(test_const_inputs<ColMajor>());
594  CALL_SUBTEST_8(test_const_inputs<RowMajor>());
595  CALL_SUBTEST_8(test_large_contraction_with_output_kernel<ColMajor>());
596  CALL_SUBTEST_8(test_large_contraction_with_output_kernel<RowMajor>());
597 
598  // Force CMake to split this test.
599  // EIGEN_SUFFIXES;1;2;3;4;5;6;7;8
600 
601 }
SqrtOutputKernel::operator()
EIGEN_ALWAYS_INLINE void operator()(const internal::blas_data_mapper< Scalar, Index, ColMajor > &output_mapper, const TensorContractionParams &, Index, Index, Index num_rows, Index num_cols) const
Definition: cxx11_tensor_contraction.cpp:517
Eigen::Tensor
The tensor class.
Definition: Tensor.h:63
test_holes
static void test_holes()
Definition: cxx11_tensor_contraction.cpp:159
test_out_of_order_contraction
static void test_out_of_order_contraction()
Definition: cxx11_tensor_contraction.cpp:271
Eigen::array
Definition: EmulateArray.h:21
test_large_contraction
static void test_large_contraction()
Definition: cxx11_tensor_contraction.cpp:353
mat1
MatrixXd mat1(size, size)
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
DataLayout
static const int DataLayout
Definition: cxx11_tensor_image_patch_sycl.cpp:24
test_scalar
static void test_scalar()
Definition: cxx11_tensor_contraction.cpp:82
m1
Matrix3d m1
Definition: IOFormat.cpp:2
Eigen::internal::isApprox
EIGEN_DEVICE_FUNC bool isApprox(const Scalar &x, const Scalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
Definition: Eigen/src/Core/MathFunctions.h:1947
Eigen::setCpuCacheSizes
void setCpuCacheSizes(std::ptrdiff_t l1, std::ptrdiff_t l2, std::ptrdiff_t l3)
Definition: products/GeneralBlockPanelKernel.h:2638
test_large_contraction_with_output_kernel
static void test_large_contraction_with_output_kernel()
Definition: cxx11_tensor_contraction.cpp:530
test_expr
static void test_expr()
Definition: cxx11_tensor_contraction.cpp:252
result
Values result
Definition: OdometryOptimize.cpp:8
Eigen::TensorContractionParams
Definition: TensorContraction.h:281
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition: split_test_helper.h:22
Eigen::DefaultDevice
Definition: TensorDeviceDefault.h:17
test_small_blocking_factors
static void test_small_blocking_factors()
Definition: cxx11_tensor_contraction.cpp:437
test_consistency
static void test_consistency()
Definition: cxx11_tensor_contraction.cpp:316
m2
MatrixType m2(n_dims)
CALL_SUBTEST_3
#define CALL_SUBTEST_3(FUNC)
Definition: split_test_helper.h:16
SqrtOutputKernel
Definition: cxx11_tensor_contraction.cpp:515
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
Eigen::Tensor::dimensions
EIGEN_DEVICE_FUNC const EIGEN_STRONG_INLINE Dimensions & dimensions() const
Definition: Tensor.h:102
test_tensor_product
static void test_tensor_product()
Definition: cxx11_tensor_contraction.cpp:467
scalar
mxArray * scalar(mxClassID classid)
Definition: matlab.h:82
test_full_redux
static void test_full_redux()
Definition: cxx11_tensor_contraction.cpp:193
l
static const Line3 l(Rot3(), 1, 1)
test_evals
static void test_evals()
Definition: cxx11_tensor_contraction.cpp:20
cholesky::expected
Matrix expected
Definition: testMatrix.cpp:971
test_tensor_vector
static void test_tensor_vector()
Definition: cxx11_tensor_contraction.cpp:413
Eigen::TensorMap
A tensor expression mapping an existing array of data.
Definition: TensorForwardDeclarations.h:52
CALL_SUBTEST_5
#define CALL_SUBTEST_5(FUNC)
Definition: split_test_helper.h:28
test_multidims
static void test_multidims()
Definition: cxx11_tensor_contraction.cpp:101
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
EIGEN_ALWAYS_INLINE
#define EIGEN_ALWAYS_INLINE
Definition: Macros.h:932
CALL_SUBTEST_6
#define CALL_SUBTEST_6(FUNC)
Definition: split_test_helper.h:34
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
test_contraction_of_contraction
static void test_contraction_of_contraction()
Definition: cxx11_tensor_contraction.cpp:219
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(cxx11_tensor_contraction)
Definition: cxx11_tensor_contraction.cpp:563
tree::f
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Definition: testExpression.cpp:218
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
Eigen::TensorBase< Tensor< Scalar_, NumIndices_, Options_, IndexType_ > >::setRandom
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Tensor< Scalar_, NumIndices_, Options_, IndexType_ > & setRandom()
Definition: TensorBase.h:996
m3
static const DiscreteKey m3(M(3), 2)
EIGEN_STATIC_ASSERT
#define EIGEN_STATIC_ASSERT(CONDITION, MSG)
Definition: StaticAssert.h:127
main.h
test_matrix_vector
static void test_matrix_vector()
Definition: cxx11_tensor_contraction.cpp:385
Eigen::Tensor::data
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar * data()
Definition: Tensor.h:104
test_const_inputs
static void test_const_inputs()
Definition: cxx11_tensor_contraction.cpp:494
Eigen::PlainObjectBase< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::data
EIGEN_DEVICE_FUNC const EIGEN_STRONG_INLINE Scalar * data() const
Definition: PlainObjectBase.h:247
Eigen::TensorEvaluator
A cost model used to limit the number of threads used for evaluating tensor expression.
Definition: TensorEvaluator.h:28
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
Eigen::ColMajor
@ ColMajor
Definition: Constants.h:319
CALL_SUBTEST_7
#define CALL_SUBTEST_7(FUNC)
Definition: split_test_helper.h:40
CALL_SUBTEST_8
#define CALL_SUBTEST_8(FUNC)
Definition: split_test_helper.h:46
DimPair
Tensor< float, 1 >::DimensionPair DimPair
Definition: cxx11_tensor_contraction.cpp:17
Eigen::TensorBase< Tensor< Scalar_, NumIndices_, Options_, IndexType_ > >::setZero
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Tensor< Scalar_, NumIndices_, Options_, IndexType_ > & setZero()
Definition: TensorBase.h:988
vec1
RowVectorXd vec1(3)
eval
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Definition: sparse_permutations.cpp:38
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
VERIFY
#define VERIFY(a)
Definition: main.h:380
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74

gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:20