1 // Copyright (c) 2018-2019 Cem Bassoy
2 //
3 // Distributed under the Boost Software License, Version 1.0. (See
4 // accompanying file LICENSE_1_0.txt or copy at
5 // http://www.boost.org/LICENSE_1_0.txt)
6 //
7 // The authors gratefully acknowledge the support of
8 // Fraunhofer and Google in producing this work
9 // which started as a Google Summer of Code project.
10 //
11
12
13 #include <iostream>
14 #include <algorithm>
15 #include <vector>
16
17 #include <boost/numeric/ublas/tensor/multiplication.hpp>
18 #include <boost/numeric/ublas/tensor/extents.hpp>
19 #include <boost/numeric/ublas/tensor/strides.hpp>
20 #include "utility.hpp"
21
22 #include <boost/test/unit_test.hpp>
23
24
25 BOOST_AUTO_TEST_SUITE (test_tensor_contraction)
26
27
28 using test_types = zip<int,long,float,double,std::complex<float>>::with_t<boost::numeric::ublas::first_order, boost::numeric::ublas::last_order>;
29
30 //using test_types = zip<int>::with_t<boost::numeric::ublas::first_order>;
31
32
33 struct fixture
34 {
35 using extents_type = boost::numeric::ublas::shape;
fixturefixture36 fixture()
37 : extents {
38 extents_type{1,1}, // 1
39 extents_type{1,2}, // 2
40 extents_type{2,1}, // 3
41 extents_type{2,3}, // 4
42 extents_type{5,4}, // 5
43 extents_type{2,3,1}, // 6
44 extents_type{4,1,3}, // 7
45 extents_type{1,2,3}, // 8
46 extents_type{4,2,3}, // 9
47 extents_type{4,2,3,5}} // 10
48 {
49 }
50 std::vector<extents_type> extents;
51 };
52
53
54
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_mtv,value,test_types,fixture)55 BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_mtv, value, test_types, fixture )
56 {
57 using namespace boost::numeric;
58 using value_type = typename value::first_type;
59 using layout_type = typename value::second_type;
60 using strides_type = ublas::strides<layout_type>;
61 using vector_type = std::vector<value_type>;
62 using extents_type = ublas::shape;
63 using extents_type_base = typename extents_type::base_type;
64 using size_type = typename extents_type_base::value_type;
65
66
67 for(auto const& na : extents) {
68
69 if(na.size() > 2)
70 continue;
71
72 auto a = vector_type(na.product(), value_type{2});
73 auto wa = strides_type(na);
74 for(auto m = 0u; m < na.size(); ++m){
75 auto nb = extents_type {na[m],1};
76 auto wb = strides_type (nb);
77 auto b = vector_type (nb.product(), value_type{1} );
78
79 auto nc_base = extents_type_base(std::max(na.size()-1, size_type{2}), 1);
80
81 for(auto i = 0u, j = 0u; i < na.size(); ++i)
82 if(i != m)
83 nc_base[j++] = na[i];
84
85 auto nc = extents_type (nc_base);
86 auto wc = strides_type (nc);
87 auto c = vector_type (nc.product(), value_type{0});
88
89 ublas::detail::recursive::mtv(
90 size_type(m),
91 c.data(), nc.data(), wc.data(),
92 a.data(), na.data(), wa.data(),
93 b.data());
94
95
96 for(auto i = 0u; i < c.size(); ++i)
97 BOOST_CHECK_EQUAL( c[i] , value_type(na[m]) * a[i] );
98
99 }
100 }
101 }
102
103
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_mtm,value,test_types,fixture)104 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_mtm, value, test_types, fixture )
105 {
106 using namespace boost::numeric;
107 using value_type = typename value::first_type;
108 using layout_type = typename value::second_type;
109 using strides_type = ublas::strides<layout_type>;
110 using vector_type = std::vector<value_type>;
111 using extents_type = ublas::shape;
112 // using extents_type_base = typename extents_type::base_type;
113
114
115 for(auto const& na : extents) {
116
117 if(na.size() != 2)
118 continue;
119
120 auto a = vector_type (na.product(), value_type{2});
121 auto wa = strides_type (na);
122
123 auto nb = extents_type {na[1],na[0]};
124 auto wb = strides_type (nb);
125 auto b = vector_type (nb.product(), value_type{1} );
126
127 auto nc = extents_type {na[0],nb[1]};
128 auto wc = strides_type (nc);
129 auto c = vector_type (nc.product());
130
131
132 ublas::detail::recursive::mtm(
133 c.data(), nc.data(), wc.data(),
134 a.data(), na.data(), wa.data(),
135 b.data(), nb.data(), wb.data());
136
137
138 for(auto i = 0u; i < c.size(); ++i)
139 BOOST_CHECK_EQUAL( c[i] , value_type(na[1]) * a[0] );
140
141
142 }
143 }
144
145
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_ttv,value,test_types,fixture)146 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_ttv, value, test_types, fixture )
147 {
148 using namespace boost::numeric;
149 using value_type = typename value::first_type;
150 using layout_type = typename value::second_type;
151 using strides_type = ublas::strides<layout_type>;
152 using vector_type = std::vector<value_type>;
153 using extents_type = ublas::shape;
154 using extents_type_base = typename extents_type::base_type;
155 using size_type = typename extents_type_base::value_type;
156
157
158 for(auto const& na : extents) {
159
160 auto a = vector_type(na.product(), value_type{2});
161 auto wa = strides_type(na);
162 for(auto m = 0u; m < na.size(); ++m){
163 auto b = vector_type (na[m], value_type{1} );
164 auto nb = extents_type {na[m],1};
165 auto wb = strides_type (nb);
166
167 auto nc_base = extents_type_base(std::max(na.size()-1, size_type(2)),1);
168
169 for(auto i = 0ul, j = 0ul; i < na.size(); ++i)
170 if(i != m)
171 nc_base[j++] = na[i];
172
173 auto nc = extents_type (nc_base);
174 auto wc = strides_type (nc);
175 auto c = vector_type (nc.product(), value_type{0});
176
177 ublas::ttv(size_type(m+1), na.size(),
178 c.data(), nc.data(), wc.data(),
179 a.data(), na.data(), wa.data(),
180 b.data(), nb.data(), wb.data());
181
182
183 for(auto i = 0u; i < c.size(); ++i)
184 BOOST_CHECK_EQUAL( c[i] , value_type(na[m]) * a[i] );
185
186 }
187 }
188 }
189
190
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_ttm,value,test_types,fixture)191 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_ttm, value, test_types, fixture )
192 {
193 using namespace boost::numeric;
194 using value_type = typename value::first_type;
195 using layout_type = typename value::second_type;
196 using strides_type = ublas::strides<layout_type>;
197 using vector_type = std::vector<value_type>;
198 using extents_type = ublas::shape;
199 using size_type = typename extents_type::value_type;
200
201
202 for(auto const& na : extents) {
203
204 auto a = vector_type(na.product(), value_type{2});
205 auto wa = strides_type(na);
206 for(auto m = 0u; m < na.size(); ++m){
207 auto nb = extents_type {na[m], na[m] };
208 auto b = vector_type (nb.product(), value_type{1} );
209 auto wb = strides_type (nb);
210
211
212 auto nc = na;
213 auto wc = strides_type (nc);
214 auto c = vector_type (nc.product(), value_type{0});
215
216 ublas::ttm(size_type(m+1), na.size(),
217 c.data(), nc.data(), wc.data(),
218 a.data(), na.data(), wa.data(),
219 b.data(), nb.data(), wb.data());
220
221 for(auto i = 0u; i < c.size(); ++i)
222 BOOST_CHECK_EQUAL( c[i] , value_type(na[m]) * a[i] );
223
224 }
225 }
226 }
227
228
229
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_ttt_permutation,value,test_types,fixture)230 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_ttt_permutation, value, test_types, fixture )
231 {
232 using namespace boost::numeric;
233 using value_type = typename value::first_type;
234 using layout_type = typename value::second_type;
235 using strides_type = ublas::strides<layout_type>;
236 using vector_type = std::vector<value_type>;
237 using extents_type = ublas::shape;
238 using size_type = typename strides_type::value_type;
239
240
241 auto compute_factorial = [](auto const& p){
242 auto f = 1ul;
243 for(auto i = 1u; i <= p; ++i)
244 f *= i;
245 return f;
246 };
247
248
249 auto compute_inverse_permutation = [](auto const& pi){
250 auto pi_inv = pi;
251 for(auto j = 0u; j < pi.size(); ++j)
252 pi_inv[pi[j]-1] = j+1;
253 return pi_inv;
254 };
255
256 auto permute_extents = [](auto const& pi, auto const& na){
257 auto nb = na;
258 assert(pi.size() == na.size());
259 for(auto j = 0u; j < pi.size(); ++j)
260 nb[j] = na[pi[j]-1];
261 return nb;
262 };
263
264
265 // left-hand and right-hand side have the
266 // the same number of elements
267
268 // computing the inner product with
269 // different permutation tuples for
270 // right-hand side
271
272 for(auto const& na : extents) {
273
274 auto wa = strides_type(na);
275 auto a = vector_type(na.product(), value_type{2});
276 auto pa = na.size();
277 auto pia = std::vector<size_type>(pa);
278 std::iota( pia.begin(), pia.end(), 1 );
279
280 auto pib = pia;
281 auto pib_inv = compute_inverse_permutation(pib);
282
283 auto f = compute_factorial(pa);
284
285 // for the number of possible permutations
286 // only permutation tuple pib is changed.
287 for(auto i = 0u; i < f; ++i) {
288
289 auto nb = permute_extents( pib, na );
290 auto wb = strides_type(nb);
291 auto b = vector_type(nb.product(), value_type{3});
292 auto pb = nb.size();
293
294 // the number of contractions is changed.
295 for( auto q = size_type(0); q <= pa; ++q) {
296
297 auto r = pa - q;
298 auto s = pb - q;
299
300 auto pc = r+s > 0 ? std::max(r+s,size_type(2)) : size_type(2);
301
302 auto nc_base = std::vector<size_type>( pc , 1 );
303
304 for(auto i = 0u; i < r; ++i)
305 nc_base[ i ] = na[ pia[i]-1 ];
306
307 for(auto i = 0u; i < s; ++i)
308 nc_base[ r + i ] = nb[ pib_inv[i]-1 ];
309
310 auto nc = extents_type ( nc_base );
311 auto wc = strides_type ( nc );
312 auto c = vector_type ( nc.product(), value_type(0) );
313
314 ublas::ttt(pa,pb,q,
315 pia.data(), pib_inv.data(),
316 c.data(), nc.data(), wc.data(),
317 a.data(), na.data(), wa.data(),
318 b.data(), nb.data(), wb.data());
319
320
321 auto acc = value_type(1);
322 for(auto i = r; i < pa; ++i)
323 acc *= value_type(na[pia[i]-1]);
324
325 for(auto i = 0ul; i < c.size(); ++i)
326 BOOST_CHECK_EQUAL( c[i] , acc * a[0] * b[0] );
327
328 }
329
330 std::next_permutation(pib.begin(), pib.end());
331 pib_inv = compute_inverse_permutation(pib);
332 }
333 }
334 }
335
336
337
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_ttt,value,test_types,fixture)338 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_ttt, value, test_types, fixture )
339 {
340 using namespace boost::numeric;
341 using value_type = typename value::first_type;
342 using layout_type = typename value::second_type;
343 using strides_type = ublas::strides<layout_type>;
344 using vector_type = std::vector<value_type>;
345 using extents_type = ublas::shape;
346 using size_type = typename strides_type::value_type;
347
348 // left-hand and right-hand side have the
349 // the same number of elements
350
351 // computing the inner product with
352 // different permutation tuples for
353 // right-hand side
354
355 for(auto const& na : extents) {
356
357 auto wa = strides_type(na);
358 auto a = vector_type(na.product(), value_type{2});
359 auto pa = na.size();
360
361 auto nb = na;
362 auto wb = strides_type(nb);
363 auto b = vector_type(nb.product(), value_type{3});
364 auto pb = nb.size();
365
366 // std::cout << "na = ";
367 // std::copy(na.begin(), na.end(), std::ostream_iterator<size_type>(std::cout, " "));
368 // std::cout << std::endl;
369
370 // std::cout << "nb = ";
371 // std::copy(nb.begin(), nb.end(), std::ostream_iterator<size_type>(std::cout, " "));
372 // std::cout << std::endl;
373
374
375 // the number of contractions is changed.
376 for( auto q = size_type(0); q <= pa; ++q) { // pa
377
378 auto r = pa - q;
379 auto s = pb - q;
380
381 auto pc = r+s > 0 ? std::max(r+s, size_type(2)) : size_type(2);
382
383 auto nc_base = std::vector<size_type>( pc , 1 );
384
385 for(auto i = 0u; i < r; ++i)
386 nc_base[ i ] = na[ i ];
387
388 for(auto i = 0u; i < s; ++i)
389 nc_base[ r + i ] = nb[ i ];
390
391 auto nc = extents_type ( nc_base );
392 auto wc = strides_type ( nc );
393 auto c = vector_type ( nc.product(), value_type{0} );
394
395 // std::cout << "nc = ";
396 // std::copy(nc.begin(), nc.end(), std::ostream_iterator<size_type>(std::cout, " "));
397 // std::cout << std::endl;
398
399 ublas::ttt(pa,pb,q,
400 c.data(), nc.data(), wc.data(),
401 a.data(), na.data(), wa.data(),
402 b.data(), nb.data(), wb.data());
403
404
405 auto acc = value_type(1);
406 for(auto i = r; i < pa; ++i)
407 acc *= value_type(na[i]);
408
409 for(auto i = 0u; i < c.size(); ++i)
410 BOOST_CHECK_EQUAL( c[i] , acc * a[0] * b[0] );
411
412 }
413
414 }
415 }
416
417
418
419
420
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_inner,value,test_types,fixture)421 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_inner, value, test_types, fixture )
422 {
423 using namespace boost::numeric;
424 using value_type = typename value::first_type;
425 using layout_type = typename value::second_type;
426 using strides_type = ublas::strides<layout_type>;
427 using vector_type = std::vector<value_type>;
428
429
430 for(auto const& n : extents) {
431
432 auto a = vector_type(n.product(), value_type{2});
433 auto b = vector_type(n.product(), value_type{3});
434 auto w = strides_type(n);
435
436 auto c = ublas::inner(n.size(), n.data(), a.data(), w.data(), b.data(), w.data(), value_type(0));
437 auto cref = std::inner_product(a.begin(), a.end(), b.begin(), value_type(0));
438
439
440 BOOST_CHECK_EQUAL( c , cref );
441
442 }
443
444 }
445
446
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_outer,value,test_types,fixture)447 BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_outer, value, test_types, fixture )
448 {
449 using namespace boost::numeric;
450 using value_type = typename value::first_type;
451 using layout_type = typename value::second_type;
452 using extents_type = ublas::shape;
453 using strides_type = ublas::strides<layout_type>;
454 using vector_type = std::vector<value_type>;
455
456
457 for(auto const& na : extents) {
458
459 auto a = vector_type(na.product(), value_type{2});
460 auto wa = strides_type(na);
461
462 for(auto const& nb : extents) {
463
464 auto b = vector_type(nb.product(), value_type{3});
465 auto wb = strides_type(nb);
466
467 auto c = vector_type(nb.product()*na.product());
468 auto nc = typename extents_type::base_type(na.size()+nb.size());
469
470 for(auto i = 0u; i < na.size(); ++i)
471 nc[i] = na[i];
472 for(auto i = 0u; i < nb.size(); ++i)
473 nc[i+na.size()] = nb[i];
474
475 auto wc = strides_type(extents_type(nc));
476
477 ublas::outer(c.data(), nc.size(), nc.data(), wc.data(),
478 a.data(), na.size(), na.data(), wa.data(),
479 b.data(), nb.size(), nb.data(), wb.data());
480
481 for(auto const& cc : c)
482 BOOST_CHECK_EQUAL( cc , a[0]*b[0] );
483 }
484
485 }
486
487 }
488
489
490 BOOST_AUTO_TEST_SUITE_END()
491
492