1 /*
2 * Copyright (C) 1998-2018 ALPS Collaboration. See COPYRIGHT.TXT
3 * All rights reserved. Use is subject to license terms. See LICENSE.TXT
4 * For use in publications, see ACKNOWLEDGE.TXT
5 */
6
7 #include <gtest/gtest.h>
8 #include <complex>
9
10 #include "alps/numeric/tensors/tensor_base.hpp"
11
12 using namespace alps::numerics::detail;
13 using namespace alps::numerics;
14
15
TEST(TensorTest,StoragesComparison)16 TEST(TensorTest, StoragesComparison) {
17 size_t N = 100;
18 data_storage<double> s(N);
19 data_view<double> v(s);
20 // copy constructed
21 data_view<double> v1(v);
22 // duplicated
23 data_view<double> v2(s);
24 for(size_t i=0; i<N; ++i) {
25 s.data(i) = i;
26 }
27 // self-comparison
28 ASSERT_EQ(s, s);
29 ASSERT_EQ(v, v);
30 // left comparison;
31 ASSERT_EQ(v, s);
32 // right comparison;
33 ASSERT_EQ(s, v);
34 // comparison with views
35 ASSERT_EQ(v1, v);
36 ASSERT_EQ(v2, v);
37
38 }
39
TEST(TensorTest,TestInitialization)40 TEST(TensorTest, TestInitialization) {
41 tensor<std::complex<double>, 3> X(1,2,3);
42 ASSERT_EQ(X.size(), 1*2*3u);
43 tensor<double, 1> U( 50 );
44 ASSERT_EQ(U.size(), 50u);
45 std::vector<double> x(100);
46 tensor<double, 2> UU(x.data(), 25, 4);
47 }
48
TEST(TensorTest,TestAssignments)49 TEST(TensorTest, TestAssignments) {
50 tensor<std::complex<double>, 3> X(std::array<size_t, 3>{{1,2,3}});
51 X( 0,0,0 ) = 10.0;
52 ASSERT_DOUBLE_EQ(X(0,0,0).real(), 10.0);
53 ASSERT_DOUBLE_EQ(X(0,0,0).imag(), 0.0);
54
55 tensor<double, 1> U(50);
56 U(3) = 555.0;
57 ASSERT_DOUBLE_EQ(U(3), 555.0);
58 ASSERT_DOUBLE_EQ(U(0), 0.0);
59
60 size_t N = 10;
61 Eigen::MatrixXd M1 = Eigen::MatrixXd::Random(N, N);
62 Eigen::MatrixXcd M2 = Eigen::MatrixXcd::Random(N, N);
63 tensor<double, 2> T1({{N, N}});
64 tensor<double, 2> T3(T1);
65 tensor<std::complex<double>, 2> T2({{N, N}});
66 for(size_t i = 0; i< N; ++i) {
67 for (size_t j = 0; j < N; ++j) {
68 T1(i,j) = M1(i,j);
69 }
70 }
71 T2 = T1;
72 ASSERT_EQ(T2, T1);
73 tensor_view<std::complex<double>, 2> V = T2;
74 ASSERT_EQ(V, T1);
75 }
76
TEST(TensorTest,TestCopyAssignments)77 TEST(TensorTest, TestCopyAssignments) {
78 size_t N = 10;
79 Eigen::MatrixXd M1 = Eigen::MatrixXd::Random(N, N);
80 Eigen::MatrixXcd M2 = Eigen::MatrixXcd::Random(N, N);
81 tensor<double, 2> T1(N, N);
82 tensor<double, 2> T3(T1);
83 tensor<double, 2> T2(N, N);
84 for(size_t i = 0; i< N; ++i) {
85 for (size_t j = 0; j < N; ++j) {
86 T1(i,j) = M1(i,j);
87 }
88 }
89 T2 = T1;
90 ASSERT_EQ(T2, T1);
91 tensor_view<double, 2> V1 = T1;
92 tensor_view<double, 2> V2 = T2;
93 ASSERT_EQ(V2, T1);
94 V1(0,0) = -15.0;
95 V2 = V1;
96 ASSERT_EQ(T2(0,0), -15.0);
97 }
98
TEST(TensorTest,TestSlices)99 TEST(TensorTest, TestSlices) {
100 tensor<std::complex<double>, 3> X(1, 2, 3);
101
102 tensor_view<std::complex<double>, 1> slice1 = X(0, 1);
103 ASSERT_EQ(slice1.dimension(), 1u);
104
105 ASSERT_EQ(slice1.shape()[0], 3u);
106
107 for (size_t i = 0; i<X.shape()[0]; ++i) {
108 tensor_view<std::complex<double>, 2> slice2 = X(i);
109 ASSERT_EQ(X.index(i), 6*i);
110 ASSERT_EQ(slice2.dimension(), 2u);
111 ASSERT_EQ(slice2.shape()[0], 2u);
112 ASSERT_EQ(slice2.shape()[1], 3u);
113 for (size_t j = 0; j < X.shape()[1]; ++j) {
114 ASSERT_EQ(X.index(i, j), 6*i + j*3);
115 for (size_t k = 0; k < X.shape()[2]; ++k) {
116 X(i, j, k) = i*X.shape()[1]*X.shape()[2] + j*X.shape()[2] + k;
117 ASSERT_DOUBLE_EQ(slice2(j, k).real(), X(i,j,k).real());
118 ASSERT_DOUBLE_EQ(slice2(j, k).imag(), X(i,j,k).imag());
119 }
120 }
121 }
122 }
123
TEST(TensorTest,TestSubSlices)124 TEST(TensorTest, TestSubSlices) {
125 tensor<double, 4> X({{3,4,5,6}});
126 for (size_t i = 0; i< X.shape()[0]; ++i) {
127 tensor_view<double, 3> Y = X(i);
128 for (size_t j = 0; j < X.shape()[1]; ++j) {
129 tensor_view<double, 2> Z = Y (j);
130 ASSERT_EQ(Z.storage().offset(), (i* X.shape()[1]+ j)* X.shape()[2]* X.shape()[3]);
131 std::vector<double> XX(X.shape()[2]* X.shape()[3], 0.0);
132 for (size_t k = 0; k < X.shape()[2]; ++k) {
133 for (size_t l = 0; l < X.shape()[3]; ++l) {
134 double value = i * X.shape()[1] * X.shape()[2] * X.shape()[3] + j * X.shape()[2] * X.shape()[3] + k * X.shape()[3] + l;
135 X(i, j, k, l) = value;
136 XX[k* X.shape()[3] + l] = value;
137 ASSERT_DOUBLE_EQ(Y(j, k, l), X(i, j, k, l));
138 ASSERT_DOUBLE_EQ(Z(k, l), X(i, j, k, l));
139 }
140 }
141 for (size_t k = 0; k < X.shape()[2]; ++k) {
142 for (size_t l = 0; l < X.shape()[3]; ++l) {
143 ASSERT_DOUBLE_EQ(XX[k* X.shape()[3] + l], X(i, j, k, l));
144 ASSERT_DOUBLE_EQ(XX[k* X.shape()[3] + l], Z(k, l));
145 }
146 }
147 }
148 }
149 }
150
TEST(TensorTest,TestMultByScalar)151 TEST(TensorTest, TestMultByScalar) {
152 tensor<double, 2> XX({{3, 4}});
153 for (size_t i = 0; i < XX.shape()[0]; ++i) {
154 for (size_t j = 0; j < XX.shape()[1]; ++j) {
155 double value = i * XX.shape()[0] + j;
156 XX(i, j) = value;
157 }
158 }
159 double mult = 12;
160 tensor<double, 2> X = XX * mult;
161 for (size_t i = 0; i < XX.shape()[0]; ++i) {
162 for (size_t j = 0; j < XX.shape()[1]; ++j) {
163 double value = i * XX.shape()[0] + j;
164 ASSERT_DOUBLE_EQ(value*mult, X(i, j));
165 ASSERT_DOUBLE_EQ(value, XX(i, j));
166 }
167 }
168 XX *= 10;
169 for (size_t i = 0; i < XX.shape()[0]; ++i) {
170 for (size_t j = 0; j < XX.shape()[1]; ++j) {
171 double value = i * XX.shape()[0] + j;
172 ASSERT_DOUBLE_EQ(value*10, XX(i, j));
173 }
174 }
175 auto Z = XX * X;
176 for (size_t i = 0; i < XX.shape()[0]; ++i) {
177 for (size_t j = 0; j < XX.shape()[1]; ++j) {
178 ASSERT_DOUBLE_EQ(Z(i, j), X(i, j) * XX(i, j));
179 }
180 }
181 }
182
TEST(TensorTest,RemoteDataRef)183 TEST(TensorTest, RemoteDataRef) {
184 size_t n = 256;
185 std::vector<double> X(n, 0.0);
186 tensor_view<double, 1> Y(X.data(), {{X.size()}});
187 tensor_view<double, 2> Z(X.data(), {{16, 16}});
188 for(size_t i = 0; i< X.size(); ++i) {
189 X[i] = i*0.5;
190 }
191 for(size_t i = 0; i< Y.shape()[0]; ++i) {
192 ASSERT_DOUBLE_EQ(X[i], Y(i));
193 }
194 for(size_t i = 0; i< Z.shape()[0]; ++i) {
195 tensor_view<double, 1> W = Z(i);
196 for (size_t j = 0; j < Z.shape()[1]; ++j) {
197 W(j) += 10;
198 ASSERT_DOUBLE_EQ(X[i* Z.shape()[1] + j], Z(i, j));
199 ASSERT_DOUBLE_EQ(X[i* Z.shape()[1] + j], W(j));
200 }
201 }
202 }
203
TEST(TensorTest,Inversion)204 TEST(TensorTest, Inversion) {
205 size_t n = 40;
206 Eigen::MatrixXd M(n,n);
207 M = Eigen::MatrixXd::Identity(n, n);
208
209 M(0, n - 1) = 1.0;
210 M(n - 1, 0) = -1.0;
211 tensor<double, 2> X({{n, n}});
212 for(size_t i = 0; i< n; ++i){
213 for (size_t j = 0; j < n; ++j) {
214 X(i, j) = M(i, j);
215 }
216 }
217 M = M.inverse();
218 X = X.inverse();
219 for(size_t i = 0; i< n; ++i){
220 for (size_t j = 0; j < n; ++j) {
221 ASSERT_DOUBLE_EQ(X(i, j), M(i, j));
222 }
223 }
224 tensor<double, 2> Z = X.dot(X.inverse());
225 for(size_t i = 0; i< n; ++i){
226 for (size_t j = 0; j < n; ++j) {
227 ASSERT_DOUBLE_EQ(Z(i, j), i==j? 1.0 : 0.0);
228 }
229 }
230 tensor<double, 2> Y = X.inverse();
231 ASSERT_DOUBLE_EQ(Y(n-1, 0), -1.0);
232 ASSERT_DOUBLE_EQ(Y(0, n-1), 1.0);
233 }
234
TEST(TensorTest,SliceInversion)235 TEST(TensorTest, SliceInversion) {
236 size_t N = 40;
237 size_t K = 10;
238 tensor<double, 3> X({{K, N, N}});
239
240 for(size_t k = 0; k < K; ++k) {
241 Eigen::MatrixXd M(N,N);
242 M = Eigen::MatrixXd::Identity(N, N);
243 M *= (k+1);
244 M(0, N - 1) = 1.0;
245 M(N - 1, 0) = -1.0;
246 tensor<double, 2> Z = X(k);
247 for(size_t i = 0; i< N; ++i){
248 for (size_t j = 0; j < N; ++j) {
249 Z(i, j) = M(i, j);
250 }
251 }
252 M = M.inverse();
253 Z = Z.inverse();
254 for(size_t i = 0; i< N; ++i){
255 for (size_t j = 0; j < N; ++j) {
256 ASSERT_DOUBLE_EQ(Z(i, j), M(i, j));
257 }
258 }
259 }
260 }
261
TEST(TensorTest,BasicArithmetics)262 TEST(TensorTest, BasicArithmetics) {
263 size_t N = 10;
264 Eigen::MatrixXd M1 = Eigen::MatrixXd::Random(N, N);
265 Eigen::MatrixXd M2 = Eigen::MatrixXd::Random(N, N);
266 tensor<double, 2> X({{N, N}});
267 tensor<double, 2> Z({{N, N}});
268 for(size_t i = 0; i< N; ++i) {
269 for (size_t j = 0; j < N; ++j) {
270 X(i,j) = M1(i,j);
271 Z(i,j) = M2(i,j);
272 }
273 }
274 auto M3 = M1 + M2;
275 auto Y = X + Z;
276 for(size_t i = 0; i< N; ++i){
277 for (size_t j = 0; j < N; ++j) {
278 ASSERT_DOUBLE_EQ(Y(i, j), M3(i, j));
279 }
280 }
281 Y -= Z;
282 for(size_t i = 0; i< N; ++i){
283 for (size_t j = 0; j < N; ++j) {
284 ASSERT_NEAR(Y(i, j), X(i, j), 1e-10);
285 }
286 }
287 }
288
TEST(TensorTest,BasicArithmeticsView)289 TEST(TensorTest, BasicArithmeticsView) {
290 size_t N = 10;
291 tensor<double, 3> W({{N,N,N}});
292 W *= 0.0;
293 for (size_t k = 0; k < N; ++k) {
294 Eigen::MatrixXd M1 = Eigen::MatrixXd::Random(N, N);
295 Eigen::MatrixXd M2 = Eigen::MatrixXd::Random(N, N);
296 tensor_view<double, 2> X = W(k);
297 tensor<double, 2> Z({{N, N}});
298 for (size_t i = 0; i < N; ++i) {
299 for (size_t j = 0; j < N; ++j) {
300 X(i, j) = M1(i, j);
301 Z(i, j) = M2(i, j);
302 }
303 }
304 auto M3 = M1 + M2;
305 tensor<double, 2> Y = X + Z;
306 for (size_t i = 0; i < N; ++i) {
307 for (size_t j = 0; j < N; ++j) {
308 ASSERT_DOUBLE_EQ(Y(i, j), M3(i, j));
309 ASSERT_DOUBLE_EQ(X(i, j), M1(i, j));
310 ASSERT_DOUBLE_EQ(Z(i, j), M2(i, j));
311 }
312 }
313 Y -= Z;
314 for (size_t i = 0; i < N; ++i) {
315 for (size_t j = 0; j < N; ++j) {
316 ASSERT_NEAR(Y(i, j), X(i, j), 1e-10);
317 }
318 }
319 }
320 }
321
TEST(TensorTest,DoubleScaleByComplex)322 TEST(TensorTest, DoubleScaleByComplex) {
323 size_t N = 10;
324 Eigen::MatrixXd M1 = Eigen::MatrixXd::Random(N, N);
325 std::complex<double> x = 5.0;
326 tensor<double, 2> X({{N, N}});
327 for(size_t i = 0; i< N; ++i) {
328 for (size_t j = 0; j < N; ++j) {
329 X(i,j) = M1(i,j);
330 }
331 }
332 auto M2 = M1 * x;
333 auto Y = X * x;
334 for(size_t i = 0; i< N; ++i){
335 for (size_t j = 0; j < N; ++j) {
336 ASSERT_DOUBLE_EQ(Y(i, j).real(), M2(i, j).real());
337 ASSERT_DOUBLE_EQ(Y(i, j).imag(), M2(i, j).imag());
338 }
339 }
340 auto M3 = M1 / x;
341 Y = X / x;
342 for(size_t i = 0; i< N; ++i){
343 for (size_t j = 0; j < N; ++j) {
344 ASSERT_DOUBLE_EQ(Y(i, j).real(), M3(i, j).real());
345 ASSERT_DOUBLE_EQ(Y(i, j).imag(), M3(i, j).imag());
346 }
347 }
348 }
349
TEST(TensorTest,DoubleScaleByComplexView)350 TEST(TensorTest, DoubleScaleByComplexView) {
351 size_t N = 10;
352 size_t M = 10;
353 std::complex<double> x = 5.0;
354 tensor<double, 3> Z({{M, N, N}});
355 for(size_t k = 0; k< M; ++k) {
356 Eigen::MatrixXd M1 = Eigen::MatrixXd::Random(N, N);
357 auto X = Z(k);
358 for (size_t i = 0; i < N; ++i) {
359 for (size_t j = 0; j < N; ++j) {
360 X(i, j) = M1(i, j);
361 }
362 }
363 auto M3 = M1 * x;
364 auto Y = X * x;
365 for (size_t i = 0; i < N; ++i) {
366 for (size_t j = 0; j < N; ++j) {
367 ASSERT_DOUBLE_EQ(Y(i, j).real(), M3(i, j).real());
368 ASSERT_DOUBLE_EQ(Y(i, j).imag(), M3(i, j).imag());
369 }
370 }
371 }
372 }
373
TEST(TensorTest,DoublePlusComplex)374 TEST(TensorTest, DoublePlusComplex) {
375 size_t N = 10;
376 Eigen::MatrixXd M1 = Eigen::MatrixXd::Random(N, N);
377 Eigen::MatrixXcd M2 = Eigen::MatrixXcd::Random(N, N);
378 tensor<double, 2> X({{N, N}});
379 tensor<std::complex<double>, 2> Z({{N, N}});
380 for(size_t i = 0; i< N; ++i) {
381 for (size_t j = 0; j < N; ++j) {
382 X(i,j) = M1(i,j);
383 Z(i,j) = M2(i,j);
384 }
385 }
386 auto M3 = M1 + M2;
387 auto Y = Z + X;
388 for(size_t i = 0; i< N; ++i){
389 for (size_t j = 0; j < N; ++j) {
390 ASSERT_NEAR(Y(i, j).real(), M3(i, j).real(), 1e-12);
391 ASSERT_NEAR(Y(i, j).imag(), M3(i, j).imag(), 1e-12);
392 }
393 }
394 Y -= Z;
395 for(size_t i = 0; i< N; ++i){
396 for (size_t j = 0; j < N; ++j) {
397 ASSERT_NEAR(Y(i, j).real(), X(i, j),1e-12);
398 }
399 }
400 }
401
TEST(TensorTest,DotProduct)402 TEST(TensorTest, DotProduct) {
403 size_t n = 40;
404 size_t l = 20;
405 // same size
406 Eigen::MatrixXd M(n,n);
407 M = Eigen::MatrixXd::Random(n, n);
408 tensor<double, 2> X1({{n, n}});
409 tensor<double, 2> X2({{n, n}});
410 for(size_t i = 0; i< n; ++i){
411 for (size_t j = 0; j < n; ++j) {
412 X1(i, j) = M(i, j);
413 X2(i, j) = M(j, i);
414 }
415 }
416 Eigen::MatrixXd M2 = M*(M.transpose());
417 auto X3 = X1.dot(X2);
418 for(size_t i = 0; i< n; ++i){
419 for (size_t j = 0; j < n; ++j) {
420 ASSERT_DOUBLE_EQ(X3(i, j), M2(i, j));
421 }
422 }
423 Eigen::MatrixXd N1(n,l);
424 Eigen::MatrixXd N2(l,n);
425 N1 = Eigen::MatrixXd::Random(n, l);
426 N2 = Eigen::MatrixXd::Random(l, n);
427 tensor<double, 2> Y1({{n, l}});
428 tensor<double, 2> Y2({{l, n}});
429 for(size_t i = 0; i< n; ++i){
430 for (size_t j = 0; j < l; ++j) {
431 Y1(i, j) = N1(i, j);
432 Y2(j, i) = N2(j, i);
433 }
434 }
435 ASSERT_ANY_THROW(Y1.dot(X1));
436 Eigen::MatrixXd N3 = N1*N2;
437 auto Y3 = Y1.dot(Y2);
438 ASSERT_EQ(Y3.shape()[0], n);
439 ASSERT_EQ(Y3.shape()[1], n);
440 for(size_t i = 0; i< n; ++i){
441 for (size_t j = 0; j < n; ++j) {
442 ASSERT_DOUBLE_EQ(Y3(i, j), N3(i, j));
443 }
444 }
445 }
446
TEST(TensorTest,StorageAssignments)447 TEST(TensorTest, StorageAssignments) {
448 size_t N = 10;
449 data_storage<double> storage_obj(N);
450 std::vector<double> buffer_obj(N, 0.0);
451 for(size_t i =0 ;i<buffer_obj.size(); ++i) {
452 buffer_obj[i] = i*15.0;
453 }
454 // data_view<double> view_obj(buffer_obj.data(), buffer_obj.size());
455 data_storage<std::complex<double> > storage_obj2 = data_view<double>(buffer_obj.data(), buffer_obj.size());
456 ASSERT_EQ(storage_obj2.size(), buffer_obj.size());
457 for(size_t i =0 ;i<buffer_obj.size(); ++i) {
458 ASSERT_DOUBLE_EQ(buffer_obj[i], storage_obj2.data(i).real());
459 }
460 ASSERT_NO_THROW(storage_obj2 = storage_obj);
461
462 }
463
TEST(TensorTest,ConstTensor)464 TEST(TensorTest, ConstTensor) {
465 size_t N = 10;
466 tensor <double, 3> X(N, N, N);
467 const tensor <double, 3> & Y = X;
468 Y( 1 ) ( 1 );
469 }
470
TEST(TensorTest,Reshape)471 TEST(TensorTest, Reshape) {
472 size_t N = 10;
473 tensor <double, 3> X(N, N, N);
474 std::array<size_t, 3> shape{1,100,5};
475 X.reshape(shape);
476 ASSERT_TRUE(X.shape()[0] == 1 && X.shape()[1] == 100 && X.shape()[2] == 5);
477 X.reshape(10,10,10);
478 ASSERT_TRUE(X.shape()[0] == 10 && X.shape()[1] == 10 && X.shape()[2] == 10);
479 }
480