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>
test_evals()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>
test_scalar()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>
test_multidims()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>
test_holes()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>
test_full_redux()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>
test_contraction_of_contraction()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>
test_expr()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>
test_out_of_order_contraction()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>
test_consistency()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>
test_large_contraction()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>
test_matrix_vector()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>
test_tensor_vector()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>
test_small_blocking_factors()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>
test_tensor_product()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>
test_const_inputs()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 
test_cxx11_tensor_contraction()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