1 //
2 // Copyright (c) 2018-2019, Cem Bassoy, cem.bassoy@gmail.com
3 //
4 // Distributed under the Boost Software License, Version 1.0. (See
5 // accompanying file LICENSE_1_0.txt or copy at
6 // http://www.boost.org/LICENSE_1_0.txt)
7 //
8 // The authors gratefully acknowledge the support of
9 // Fraunhofer IOSB, Ettlingen, Germany
10 //
11
12
13 #include <boost/numeric/ublas/tensor.hpp>
14 #include <boost/numeric/ublas/matrix.hpp>
15 #include <boost/numeric/ublas/vector.hpp>
16 #include <iostream>
17
main()18 int main()
19 {
20 using namespace boost::numeric::ublas;
21
22 using format_t = column_major;
23 using value_t = float; // std::complex<double>;
24 using tensor_t = tensor<value_t,format_t>;
25 using matrix_t = matrix<value_t,format_t>;
26 using vector_t = vector<value_t>;
27
28 // Tensor-Vector-Multiplications - Including Transposition
29 {
30
31 auto n = shape{3,4,2};
32 auto A = tensor_t(n,2);
33 auto q = 0u; // contraction mode
34
35 // C1(j,k) = T2(j,k) + A(i,j,k)*T1(i);
36 q = 1u;
37 tensor_t C1 = matrix_t(n[1],n[2],2) + prod(A,vector_t(n[q-1],1),q);
38
39 // C2(i,k) = A(i,j,k)*T1(j) + 4;
40 q = 2u;
41 tensor_t C2 = prod(A,vector_t(n[q-1],1),q) + 4;
42
43 // C3() = A(i,j,k)*T1(i)*T2(j)*T2(k);
44 tensor_t C3 = prod(prod(prod(A,vector_t(n[0],1),1),vector_t(n[1],1),1),vector_t(n[2],1),1);
45
46 // C4(i,j) = A(k,i,j)*T1(k) + 4;
47 q = 1u;
48 tensor_t C4 = prod(trans(A,{2,3,1}),vector_t(n[2],1),q) + 4;
49
50
51 // formatted output
52 std::cout << "% --------------------------- " << std::endl;
53 std::cout << "% --------------------------- " << std::endl << std::endl;
54 std::cout << "% C1(j,k) = T2(j,k) + A(i,j,k)*T1(i);" << std::endl << std::endl;
55 std::cout << "C1=" << C1 << ";" << std::endl << std::endl;
56
57 // formatted output
58 std::cout << "% --------------------------- " << std::endl;
59 std::cout << "% --------------------------- " << std::endl << std::endl;
60 std::cout << "% C2(i,k) = A(i,j,k)*T1(j) + 4;" << std::endl << std::endl;
61 std::cout << "C2=" << C2 << ";" << std::endl << std::endl;
62
63 // formatted output
64 std::cout << "% --------------------------- " << std::endl;
65 std::cout << "% --------------------------- " << std::endl << std::endl;
66 std::cout << "% C3() = A(i,j,k)*T1(i)*T2(j)*T2(k);" << std::endl << std::endl;
67 std::cout << "C3()=" << C3(0) << ";" << std::endl << std::endl;
68
69 // formatted output
70 std::cout << "% --------------------------- " << std::endl;
71 std::cout << "% --------------------------- " << std::endl << std::endl;
72 std::cout << "% C4(i,j) = A(k,i,j)*T1(k) + 4;" << std::endl << std::endl;
73 std::cout << "C4=" << C4 << ";" << std::endl << std::endl;
74
75 }
76
77
78 // Tensor-Matrix-Multiplications - Including Transposition
79 {
80
81 auto n = shape{3,4,2};
82 auto A = tensor_t(n,2);
83 auto m = 5u;
84 auto q = 0u; // contraction mode
85
86 // C1(l,j,k) = T2(l,j,k) + A(i,j,k)*T1(l,i);
87 q = 1u;
88 tensor_t C1 = tensor_t(shape{m,n[1],n[2]},2) + prod(A,matrix_t(m,n[q-1],1),q);
89
90 // C2(i,l,k) = A(i,j,k)*T1(l,j) + 4;
91 q = 2u;
92 tensor_t C2 = prod(A,matrix_t(m,n[q-1],1),q) + 4;
93
94 // C3(i,l1,l2) = A(i,j,k)*T1(l1,j)*T2(l2,k);
95 q = 3u;
96 tensor_t C3 = prod(prod(A,matrix_t(m+1,n[q-2],1),q-1),matrix_t(m+2,n[q-1],1),q);
97
98 // C4(i,l1,l2) = A(i,j,k)*T2(l2,k)*T1(l1,j);
99 tensor_t C4 = prod(prod(A,matrix_t(m+2,n[q-1],1),q),matrix_t(m+1,n[q-2],1),q-1);
100
101 // C5(i,k,l) = A(i,k,j)*T1(l,j) + 4;
102 q = 3u;
103 tensor_t C5 = prod(trans(A,{1,3,2}),matrix_t(m,n[1],1),q) + 4;
104
105 // formatted output
106 std::cout << "% --------------------------- " << std::endl;
107 std::cout << "% --------------------------- " << std::endl << std::endl;
108 std::cout << "% C1(l,j,k) = T2(l,j,k) + A(i,j,k)*T1(l,i);" << std::endl << std::endl;
109 std::cout << "C1=" << C1 << ";" << std::endl << std::endl;
110
111 // formatted output
112 std::cout << "% --------------------------- " << std::endl;
113 std::cout << "% --------------------------- " << std::endl << std::endl;
114 std::cout << "% C2(i,l,k) = A(i,j,k)*T1(l,j) + 4;" << std::endl << std::endl;
115 std::cout << "C2=" << C2 << ";" << std::endl << std::endl;
116
117 // formatted output
118 std::cout << "% --------------------------- " << std::endl;
119 std::cout << "% --------------------------- " << std::endl << std::endl;
120 std::cout << "% C3(i,l1,l2) = A(i,j,k)*T1(l1,j)*T2(l2,k);" << std::endl << std::endl;
121 std::cout << "C3=" << C3 << ";" << std::endl << std::endl;
122
123 // formatted output
124 std::cout << "% --------------------------- " << std::endl;
125 std::cout << "% --------------------------- " << std::endl << std::endl;
126 std::cout << "% C4(i,l1,l2) = A(i,j,k)*T2(l2,k)*T1(l1,j);" << std::endl << std::endl;
127 std::cout << "C4=" << C4 << ";" << std::endl << std::endl;
128 std::cout << "% C3 and C4 should have the same values, true? " << std::boolalpha << (C3 == C4) << "!" << std::endl;
129
130
131 // formatted output
132 std::cout << "% --------------------------- " << std::endl;
133 std::cout << "% --------------------------- " << std::endl << std::endl;
134 std::cout << "% C5(i,k,l) = A(i,k,j)*T1(l,j) + 4;" << std::endl << std::endl;
135 std::cout << "C5=" << C5 << ";" << std::endl << std::endl;
136 }
137
138
139
140
141
142 // Tensor-Tensor-Multiplications Including Transposition
143 {
144
145 using perm_t = std::vector<std::size_t>;
146
147 auto na = shape{3,4,5};
148 auto nb = shape{4,6,3,2};
149 auto A = tensor_t(na,2);
150 auto B = tensor_t(nb,3);
151
152
153 // C1(j,l) = T(j,l) + A(i,j,k)*A(i,j,l) + 5;
154 tensor_t C1 = tensor_t(shape{na[2],na[2]},2) + prod(A,A,perm_t{1,2}) + 5;
155
156 // formatted output
157 std::cout << "% --------------------------- " << std::endl;
158 std::cout << "% --------------------------- " << std::endl << std::endl;
159 std::cout << "% C1(k,l) = T(k,l) + A(i,j,k)*A(i,j,l) + 5;" << std::endl << std::endl;
160 std::cout << "C1=" << C1 << ";" << std::endl << std::endl;
161
162
163 // C2(k,l,m) = T(k,l,m) + A(i,j,k)*B(j,l,i,m) + 5;
164 tensor_t C2 = tensor_t(shape{na[2],nb[1],nb[3]},2) + prod(A,B,perm_t{1,2},perm_t{3,1}) + 5;
165
166 // formatted output
167 std::cout << "% --------------------------- " << std::endl;
168 std::cout << "% --------------------------- " << std::endl << std::endl;
169 std::cout << "% C2(k,l,m) = T(k,l,m) + A(i,j,k)*B(j,l,i,m) + 5;" << std::endl << std::endl;
170 std::cout << "C2=" << C2 << ";" << std::endl << std::endl;
171
172
173 // C3(k,l,m) = T(k,l,m) + A(i,j,k)*trans(B(j,l,i,m),{2,3,1,4})+ 5;
174 tensor_t C3 = tensor_t(shape{na[2],nb[1],nb[3]},2) + prod(A,trans(B,{2,3,1,4}),perm_t{1,2}) + 5;
175
176 // formatted output
177 std::cout << "% --------------------------- " << std::endl;
178 std::cout << "% --------------------------- " << std::endl << std::endl;
179 std::cout << "% C3(k,l,m) = T(k,l,m) + A(i,j,k)*trans(B(j,l,i,m),{2,3,1,4})+ 5;" << std::endl << std::endl;
180 std::cout << "C3=" << C3 << ";" << std::endl << std::endl;
181
182 }
183 }
184