1 /*
2 * @BEGIN LICENSE
3 *
4 * ambit: C++ library for the implementation of tensor product calculations
5 * through a clean, concise user interface.
6 *
7 * Copyright (c) 2014-2017 Ambit developers.
8 *
9 * The copyrights for code used from other parties are included in
10 * the corresponding files.
11 *
12 * This file is part of ambit.
13 *
14 * Ambit is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU Lesser General Public License as published by
16 * the Free Software Foundation, version 3.
17 *
18 * Ambit is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU Lesser General Public License for more details.
22 *
23 * You should have received a copy of the GNU Lesser General Public License along
24 * with ambit; if not, write to the Free Software Foundation, Inc.,
25 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 *
27 * @END LICENSE
28 */
29
30 #include <ambit/blocked_tensor.h>
31 #include <cmath>
32 #include <cstdlib>
33 #include <stdexcept>
34
35 #define ANSI_COLOR_RED "\x1b[31m"
36 #define ANSI_COLOR_GREEN "\x1b[32m"
37 #define ANSI_COLOR_YELLOW "\x1b[33m"
38 #define ANSI_COLOR_BLUE "\x1b[34m"
39 #define ANSI_COLOR_MAGENTA "\x1b[35m"
40 #define ANSI_COLOR_CYAN "\x1b[36m"
41 #define ANSI_COLOR_RESET "\x1b[0m"
42
43 #define MAXTWO 10
44 #define MAXFOUR 10
45
46 double a1[MAXTWO];
47 double a2[MAXTWO][MAXTWO];
48 double b2[MAXTWO][MAXTWO];
49 double c2[MAXTWO][MAXTWO];
50 double d2[MAXTWO][MAXTWO];
51 double e2[MAXTWO][MAXTWO];
52 double f2[MAXTWO][MAXTWO];
53 double a4[MAXFOUR][MAXFOUR][MAXFOUR][MAXFOUR];
54 double b4[MAXFOUR][MAXFOUR][MAXFOUR][MAXFOUR];
55 double c4[MAXFOUR][MAXFOUR][MAXFOUR][MAXFOUR];
56 double d4[MAXFOUR][MAXFOUR][MAXFOUR][MAXFOUR];
57
58 using namespace ambit;
59
60 /// Expected relative accuracy for numerical exactness
61 const double epsilon = 1.0E-14;
62 const double zero = 1.0E-14;
63
64 /// Scheme to categorize expected vs. actual op behavior
65 enum TestResult
66 {
67 kPass,
68 kFail,
69 kException
70 };
71
72 /// Global alpha parameter
73 double alpha = 1.0;
74 /// Global beta parameter
75 double beta = 0.0;
76 /// 0 - explicit call, 1 - = OO, 2 - += OO, 3 - -= OO
77 int mode = 0;
78
79 TensorType tensor_type = CoreTensor;
80
initialize_random(Tensor & tensor,double matrix[MAXTWO])81 void initialize_random(Tensor &tensor, double matrix[MAXTWO])
82 {
83 size_t n0 = tensor.dims()[0];
84 std::vector<double> &vec = tensor.data();
85 for (size_t i = 0; i < n0; ++i)
86 {
87 double randnum = double(std::rand()) / double(RAND_MAX);
88 matrix[i] = randnum;
89 vec[i] = randnum;
90 }
91 }
92
initialize_random(Tensor & tensor,double matrix[MAXTWO][MAXTWO])93 void initialize_random(Tensor &tensor, double matrix[MAXTWO][MAXTWO])
94 {
95 size_t n0 = tensor.dims()[0];
96 size_t n1 = tensor.dims()[1];
97 std::vector<double> &vec = tensor.data();
98 for (size_t i = 0, ij = 0; i < n0; ++i)
99 {
100 for (size_t j = 0; j < n1; ++j, ++ij)
101 {
102 double randnum = double(std::rand()) / double(RAND_MAX);
103 matrix[i][j] = randnum;
104 vec[ij] = randnum;
105 }
106 }
107 }
108
initialize_random(Tensor & tensor,double matrix[MAXFOUR][MAXFOUR][MAXFOUR][MAXFOUR])109 void initialize_random(Tensor &tensor,
110 double matrix[MAXFOUR][MAXFOUR][MAXFOUR][MAXFOUR])
111 {
112 size_t n0 = tensor.dims()[0];
113 size_t n1 = tensor.dims()[1];
114 size_t n2 = tensor.dims()[2];
115 size_t n3 = tensor.dims()[3];
116
117 std::vector<double> &vec = tensor.data();
118 for (size_t i = 0, ijkl = 0; i < n0; ++i)
119 {
120 for (size_t j = 0; j < n1; ++j)
121 {
122 for (size_t k = 0; k < n2; ++k)
123 {
124 for (size_t l = 0; l < n3; ++l, ++ijkl)
125 {
126 double randnum = double(std::rand()) / double(RAND_MAX);
127 matrix[i][j][k][l] = randnum;
128 vec[ijkl] = randnum;
129 }
130 }
131 }
132 }
133 }
134
difference(Tensor & tensor,double matrix[MAXTWO])135 std::pair<double, double> difference(Tensor &tensor, double matrix[MAXTWO])
136 {
137 size_t n0 = tensor.dims()[0];
138
139 const std::vector<double> &result = tensor.data();
140
141 double sum_diff = 0.0;
142 double max_diff = 0.0;
143 for (size_t i = 0; i < n0; ++i)
144 {
145 double diff = std::fabs(matrix[i] - result[i]);
146 sum_diff += diff;
147 max_diff = std::max(diff, max_diff);
148 }
149 return std::make_pair(sum_diff, max_diff);
150 }
151
difference(Tensor & tensor,double matrix[MAXTWO][MAXTWO])152 std::pair<double, double> difference(Tensor &tensor,
153 double matrix[MAXTWO][MAXTWO])
154 {
155 size_t n0 = tensor.dims()[0];
156 size_t n1 = tensor.dims()[1];
157
158 const std::vector<double> &result = tensor.data();
159
160 double sum_diff = 0.0;
161 double max_diff = 0.0;
162 for (size_t i = 0, ij = 0; i < n0; ++i)
163 {
164 for (size_t j = 0; j < n1; ++j, ++ij)
165 {
166 double diff = std::fabs(matrix[i][j] - result[ij]);
167 sum_diff += diff;
168 max_diff = std::max(diff, max_diff);
169 }
170 }
171 return std::make_pair(sum_diff, max_diff);
172 }
173
174 std::pair<double, double>
difference(Tensor & tensor,double matrix[MAXFOUR][MAXFOUR][MAXFOUR][MAXFOUR])175 difference(Tensor &tensor, double matrix[MAXFOUR][MAXFOUR][MAXFOUR][MAXFOUR])
176 {
177 size_t n0 = tensor.dims()[0];
178 size_t n1 = tensor.dims()[1];
179 size_t n2 = tensor.dims()[2];
180 size_t n3 = tensor.dims()[3];
181
182 const std::vector<double> &result = tensor.data();
183
184 double sum_diff = 0.0;
185 double max_diff = 0.0;
186
187 for (size_t i = 0, ijkl = 0; i < n0; ++i)
188 {
189 for (size_t j = 0; j < n1; ++j)
190 {
191 for (size_t k = 0; k < n2; ++k)
192 {
193 for (size_t l = 0; l < n3; ++l, ++ijkl)
194 {
195 double diff = std::fabs(matrix[i][j][k][l] - result[ijkl]);
196 sum_diff += diff;
197 max_diff = std::max(diff, max_diff);
198 }
199 }
200 }
201 }
202 return std::make_pair(sum_diff, max_diff);
203 }
204
build_and_fill(const std::string & name,const Dimension & dims,double matrix[MAXTWO])205 Tensor build_and_fill(const std::string &name, const Dimension &dims,
206 double matrix[MAXTWO])
207 {
208 Tensor T = Tensor::build(tensor_type, name, dims);
209 initialize_random(T, matrix);
210 std::pair<double, double> a_diff = difference(T, matrix);
211 if (std::fabs(a_diff.second) > zero)
212 throw std::runtime_error("Tensor and standard matrix don't match.");
213 return T;
214 }
215
build_and_fill(const std::string & name,const Dimension & dims,double matrix[MAXTWO][MAXTWO])216 Tensor build_and_fill(const std::string &name, const Dimension &dims,
217 double matrix[MAXTWO][MAXTWO])
218 {
219 Tensor T = Tensor::build(tensor_type, name, dims);
220 initialize_random(T, matrix);
221 std::pair<double, double> a_diff = difference(T, matrix);
222 if (std::fabs(a_diff.second) > zero)
223 throw std::runtime_error("Tensor and standard matrix don't match.");
224 return T;
225 }
226
build_and_fill(const std::string & name,const Dimension & dims,double matrix[MAXFOUR][MAXFOUR][MAXFOUR][MAXFOUR])227 Tensor build_and_fill(const std::string &name, const Dimension &dims,
228 double matrix[MAXFOUR][MAXFOUR][MAXFOUR][MAXFOUR])
229 {
230 Tensor T = Tensor::build(tensor_type, name, dims);
231 initialize_random(T, matrix);
232 std::pair<double, double> a_diff = difference(T, matrix);
233 if (std::fabs(a_diff.second) > zero)
234 throw std::runtime_error("Tensor and standard matrix don't match.");
235 return T;
236 }
237
test_mo_space()238 double test_mo_space()
239 {
240 MOSpace alpha_occ("o", "i,j,k,l", {0, 1, 2, 3, 4}, AlphaSpin);
241 MOSpace alpha_vir("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
242 return 0.0;
243 }
244
test_add_mo_space()245 double test_add_mo_space()
246 {
247 BlockedTensor::reset_mo_spaces();
248 BlockedTensor::add_mo_space("o", "i,j,k,l", {0, 1, 2, 3, 4}, AlphaSpin);
249 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
250 BlockedTensor::add_composite_mo_space("g", "p,q,r,s,t", {"o", "v"});
251 return 0.0;
252 }
253
test_add_mo_space_nonexisting_space()254 double test_add_mo_space_nonexisting_space()
255 {
256 BlockedTensor::reset_mo_spaces();
257 BlockedTensor::add_mo_space("o", "i,j,k,l", {0, 1, 2, 3, 4}, AlphaSpin);
258 BlockedTensor::add_composite_mo_space("g", "p,q,r,s,t", {"o", "v"});
259 return 0.0;
260 }
261
test_add_mo_space_repeated_index1()262 double test_add_mo_space_repeated_index1()
263 {
264 BlockedTensor::reset_mo_spaces();
265 BlockedTensor::add_mo_space("o", "i,j,k,l", {0, 1, 2, 3, 4}, AlphaSpin);
266 BlockedTensor::add_mo_space("v", "a,b,c,i", {5, 6, 7, 8, 9}, AlphaSpin);
267 return 0.0;
268 }
269
test_add_mo_space_repeated_index2()270 double test_add_mo_space_repeated_index2()
271 {
272 BlockedTensor::reset_mo_spaces();
273 BlockedTensor::add_mo_space("o", "i,j,k,l", {0, 1, 2, 3, 4}, AlphaSpin);
274 BlockedTensor::add_mo_space("v", "a,b,c,a", {5, 6, 7, 8, 9}, AlphaSpin);
275 return 0.0;
276 }
277
test_add_mo_space_repeated_index3()278 double test_add_mo_space_repeated_index3()
279 {
280 BlockedTensor::reset_mo_spaces();
281 BlockedTensor::add_mo_space("o", "i,j,k,l", {0, 1, 2, 3, 4}, AlphaSpin);
282 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
283 BlockedTensor::add_composite_mo_space("g", "p,q,r,s,c", {"o", "v"});
284 return 0.0;
285 }
286
test_add_mo_space_no_name1()287 double test_add_mo_space_no_name1()
288 {
289 BlockedTensor::reset_mo_spaces();
290 BlockedTensor::add_mo_space("", "i,j,k", {0, 1, 2, 3, 4}, AlphaSpin);
291 return 0.0;
292 }
293
test_add_mo_space_no_name2()294 double test_add_mo_space_no_name2()
295 {
296 BlockedTensor::reset_mo_spaces();
297 BlockedTensor::add_mo_space("o", "i,j,k,l", {0, 1, 2, 3, 4}, AlphaSpin);
298 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
299 BlockedTensor::add_composite_mo_space("", "p,q,r,s", {"o", "v"});
300 return 0.0;
301 }
302
test_add_mo_space_no_index1()303 double test_add_mo_space_no_index1()
304 {
305 BlockedTensor::reset_mo_spaces();
306 BlockedTensor::add_mo_space("o", "", {0, 1, 2, 3, 4}, AlphaSpin);
307 return 0.0;
308 }
309
test_add_mo_space_no_index2()310 double test_add_mo_space_no_index2()
311 {
312 BlockedTensor::reset_mo_spaces();
313 BlockedTensor::add_mo_space("o", "i,j,k,l", {0, 1, 2, 3, 4}, AlphaSpin);
314 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
315 BlockedTensor::add_composite_mo_space("g", "", {"o", "v"});
316 return 0.0;
317 }
318
test_add_mo_space_no_mos()319 double test_add_mo_space_no_mos()
320 {
321 BlockedTensor::reset_mo_spaces();
322 BlockedTensor::add_mo_space("o", "i,j,k,l", {}, AlphaSpin);
323 return 0.0;
324 }
325
test_add_mo_space_repeated_space1()326 double test_add_mo_space_repeated_space1()
327 {
328 BlockedTensor::reset_mo_spaces();
329 BlockedTensor::add_mo_space("o", "i,j,k,l", {0, 1, 2, 3, 4}, AlphaSpin);
330 BlockedTensor::add_mo_space("o", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
331 return 0.0;
332 }
333
test_add_mo_space_repeated_space2()334 double test_add_mo_space_repeated_space2()
335 {
336 BlockedTensor::reset_mo_spaces();
337 BlockedTensor::add_mo_space("o", "i,j,k,l", {0, 1, 2, 3, 4}, AlphaSpin);
338 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
339 BlockedTensor::add_composite_mo_space("o", "p,q,r,s", {"o", "v"});
340 return 0.0;
341 }
342
test_block_creation1()343 double test_block_creation1()
344 {
345 BlockedTensor::reset_mo_spaces();
346 BlockedTensor::add_mo_space("o", "i,j", {0, 1, 2}, AlphaSpin);
347 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
348 BlockedTensor::build(CoreTensor, "T", {"oo", "vv"});
349 return 0.0;
350 }
351
test_block_creation2()352 double test_block_creation2()
353 {
354 BlockedTensor::reset_mo_spaces();
355 BlockedTensor::add_mo_space("o", "i,j", {0, 1, 2}, AlphaSpin);
356 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
357 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
358 BlockedTensor::build(CoreTensor, "F", {"gg"});
359 BlockedTensor::build(CoreTensor, "V", {"gggg"});
360 return 0.0;
361 }
362
test_block_creation3()363 double test_block_creation3()
364 {
365 BlockedTensor::reset_mo_spaces();
366 BlockedTensor::add_mo_space("c", "m,n", {0, 1, 2}, AlphaSpin);
367 BlockedTensor::add_mo_space("a", "u,v", {3, 4}, AlphaSpin);
368 BlockedTensor::add_mo_space("v", "e,f", {5, 6, 7, 8, 9}, AlphaSpin);
369 BlockedTensor::add_composite_mo_space("h", "i,j,k,l", {"c", "a"});
370 BlockedTensor::add_composite_mo_space("p", "a,b,c,d", {"a", "v"});
371 BlockedTensor::build(CoreTensor, "T1", {"hp"});
372 BlockedTensor::build(CoreTensor, "T2", {"hhpp"});
373 return 0.0;
374 }
375
test_block_creation_bad_rank()376 double test_block_creation_bad_rank()
377 {
378 BlockedTensor::reset_mo_spaces();
379 BlockedTensor::add_mo_space("o", "i,j", {0, 1, 2}, AlphaSpin);
380 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
381 BlockedTensor::build(CoreTensor, "T", {"oo", "ovv"});
382 return 0.0;
383 }
384
test_block_norm_1()385 double test_block_norm_1()
386 {
387 BlockedTensor::reset_mo_spaces();
388 BlockedTensor::add_mo_space("a", "u,v", {2, 3, 4}, NoSpin);
389 BlockedTensor::add_mo_space("v", "e,f", {5, 6, 7, 8, 9}, NoSpin);
390 BlockedTensor T2 = BlockedTensor::build(CoreTensor, "T2", {"aavv"});
391 T2.set(0.5);
392 double diff = T2.norm(1) - 112.5;
393 return diff;
394 }
395
test_block_norm_2()396 double test_block_norm_2()
397 {
398 BlockedTensor::reset_mo_spaces();
399 BlockedTensor::add_mo_space("a", "u,v", {2, 3, 4}, AlphaSpin);
400 BlockedTensor::add_mo_space("v", "e,f", {5, 6, 7, 8, 9}, AlphaSpin);
401 BlockedTensor T2 = BlockedTensor::build(CoreTensor, "T2", {"aavv"});
402 T2.set(0.5);
403 double diff = T2.norm(2) - 7.5;
404 return diff;
405 }
406
test_block_norm_3()407 double test_block_norm_3()
408 {
409 BlockedTensor::reset_mo_spaces();
410 BlockedTensor::add_mo_space("a", "u,v", {2, 3, 4}, AlphaSpin);
411 BlockedTensor::add_mo_space("v", "e,f", {5, 6, 7, 8, 9}, AlphaSpin);
412 BlockedTensor T2 = BlockedTensor::build(CoreTensor, "T2", {"aavv"});
413 T2.set(0.5);
414 double diff = T2.norm(0) - 0.5;
415 return diff;
416 }
417
test_block_zero()418 double test_block_zero()
419 {
420 BlockedTensor::reset_mo_spaces();
421 BlockedTensor::add_mo_space("a", "u,v", {2, 3, 4}, AlphaSpin);
422 BlockedTensor::add_mo_space("v", "e,f", {5, 6, 7, 8, 9}, AlphaSpin);
423 BlockedTensor T2 = BlockedTensor::build(CoreTensor, "T2", {"aavv"});
424 T2.set(0.5);
425 T2.zero();
426 return T2.norm(2);
427 }
428
test_block_scale()429 double test_block_scale()
430 {
431 BlockedTensor::reset_mo_spaces();
432 BlockedTensor::add_mo_space("a", "u,v", {2, 3, 4}, AlphaSpin);
433 BlockedTensor::add_mo_space("v", "e,f", {5, 6, 7, 8, 9}, AlphaSpin);
434 BlockedTensor T2 = BlockedTensor::build(CoreTensor, "T2", {"aavv"});
435 T2.set(2.0);
436 T2.scale(0.25);
437 double diff = T2.norm(2) - 7.5;
438 return diff;
439 }
440
test_block_labels1()441 double test_block_labels1()
442 {
443 BlockedTensor::reset_mo_spaces();
444 BlockedTensor::add_mo_space("o", "i,j", {0, 1, 2}, AlphaSpin);
445 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
446 BlockedTensor T = BlockedTensor::build(CoreTensor, "T", {"oo", "vv"});
447 T("ij");
448 return 0.0;
449 }
450
test_block_retrive_block1()451 double test_block_retrive_block1()
452 {
453 BlockedTensor::reset_mo_spaces();
454 BlockedTensor::add_mo_space("o", "i,j", {0, 1, 2}, AlphaSpin);
455 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
456 BlockedTensor T = BlockedTensor::build(CoreTensor, "T", {"oo", "vv"});
457 T.block("oo");
458 return 0.0;
459 }
460
test_block_retrive_block2()461 double test_block_retrive_block2()
462 {
463 BlockedTensor::reset_mo_spaces();
464 BlockedTensor::add_mo_space("o", "i,j", {0, 1, 2}, AlphaSpin);
465 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
466 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
467 BlockedTensor T = BlockedTensor::build(CoreTensor, "T", {"oo", "vv"});
468 T.block("og");
469 return 0.0;
470 }
471
test_block_retrive_block3()472 double test_block_retrive_block3()
473 {
474 BlockedTensor::reset_mo_spaces();
475 BlockedTensor::add_mo_space("o", "i,j", {0, 1, 2}, AlphaSpin);
476 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
477 BlockedTensor T = BlockedTensor::build(CoreTensor, "T", {"oo", "vv"});
478 T.block("ov");
479 return 0.0;
480 }
481
test_block_retrive_block4()482 double test_block_retrive_block4()
483 {
484 BlockedTensor::reset_mo_spaces();
485 BlockedTensor::add_mo_space("o", "i,j", {0, 1, 2}, AlphaSpin);
486 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
487 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
488 BlockedTensor T = BlockedTensor::build(CoreTensor, "T", {"oo", "vv"});
489 T.block("");
490 return 0.0;
491 }
492
493 // double test_copy()
494 //{
495 // BlockedTensor::reset_mo_spaces();
496 // BlockedTensor::add_mo_space("o","i,j",{0,1,2},AlphaSpin);
497 // BlockedTensor::add_mo_space("v","a,b,c,d",{5,6,7,8,9},AlphaSpin);
498
499 // BlockedTensor A =
500 // BlockedTensor::build(CoreTensor,"A",{"oo","vv","ov","vo"});
501 // BlockedTensor C =
502 // BlockedTensor::build(CoreTensor,"C",{"oo","vv","ov","vo"});
503
504 // size_t no = 3;
505 // size_t nv = 5;
506
507 // Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
508 // Tensor Aov_t = build_and_fill("Aov", {no, nv}, b2);
509 // Tensor Avo_t = build_and_fill("Avo", {nv, no}, c2);
510 // Tensor Avv_t = build_and_fill("Avv", {nv, nv}, d2);
511
512 // A.block("oo")("pq") = Aoo_t("pq");
513 // A.block("ov")("pq") = Aov_t("pq");
514 // A.block("vo")("pq") = Avo_t("pq");
515 // A.block("vv")("pq") = Avv_t("pq");
516
517 // C.copy(A);
518
519 // A.scale(2.0);
520
521 // Tensor Coo = C.block("oo");
522 // Tensor Cov = C.block("ov");
523 // Tensor Cvo = C.block("vo");
524 // Tensor Cvv = C.block("vv");
525
526 // double diff_oo = difference(Coo, a2).second;
527 // double diff_vo = difference(Cvo, b2).second;
528 // double diff_ov = difference(Cov, c2).second;
529 // double diff_vv = difference(Cvv, d2).second;
530
531 // return diff_oo + diff_vo + diff_ov + diff_vv;
532 //}
533
test_block_iterator_1()534 double test_block_iterator_1()
535 {
536 BlockedTensor::reset_mo_spaces();
537 BlockedTensor::add_mo_space("a", "u,v", {2, 3, 4}, NoSpin);
538 BlockedTensor::add_mo_space("v", "e,f", {5, 6, 7, 8, 9}, NoSpin);
539 BlockedTensor T2 = BlockedTensor::build(CoreTensor, "T2", {"aaaa", "vvvv"});
540 T2.iterate([](const std::vector<size_t> &indices,
541 const std::vector<SpinType> &spin, double &value)
542 {
543 bool add = true;
544 for (size_t k : indices)
545 {
546 if (k > 4)
547 add = false;
548 }
549 if (add)
550 {
551 value = 1.0;
552 }
553 });
554 double diff = T2.norm(1) - 81.0;
555 return diff;
556 }
557
test_Cij_equal_Aji()558 double test_Cij_equal_Aji()
559 {
560 BlockedTensor::reset_mo_spaces();
561 BlockedTensor::add_mo_space("o", "i,j", {0, 1, 2}, AlphaSpin);
562 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
563
564 BlockedTensor A =
565 BlockedTensor::build(CoreTensor, "A", {"oo", "vv", "ov", "vo"});
566 BlockedTensor C =
567 BlockedTensor::build(CoreTensor, "C", {"oo", "vv", "ov", "vo"});
568
569 Tensor Aoo = A.block("oo");
570 Tensor Coo = C.block("oo");
571
572 size_t no = 3;
573 size_t nv = 5;
574
575 Tensor Aoo_t = build_and_fill("A", {no, no}, a2);
576 Tensor Coo_t = build_and_fill("C", {no, no}, c2);
577
578 Aoo("ij") = Aoo_t("ij");
579 Coo("ij") = Coo_t("ij");
580
581 C("ij") = A("ji");
582
583 for (size_t i = 0; i < no; ++i)
584 {
585 for (size_t j = 0; j < no; ++j)
586 {
587 c2[i][j] = a2[j][i];
588 }
589 }
590
591 double diff_oo = difference(Coo, c2).second;
592
593 Tensor Aov = A.block("ov");
594 Tensor Cvo = C.block("vo");
595
596 Tensor Aov_t = build_and_fill("A", {no, nv}, a2);
597 Tensor Cvo_t = build_and_fill("C", {nv, no}, c2);
598
599 Aov("ij") = Aov_t("ij");
600 Cvo("ij") = Cvo_t("ij");
601
602 C("ai") = A("ia");
603
604 for (size_t i = 0; i < no; ++i)
605 {
606 for (size_t a = 0; a < nv; ++a)
607 {
608 c2[a][i] = a2[i][a];
609 }
610 }
611
612 double diff_vo = difference(Cvo, c2).second;
613
614 Tensor Avv = A.block("vv");
615 Tensor Cvv = C.block("vv");
616
617 Tensor Avv_t = build_and_fill("A", {nv, nv}, a2);
618 Tensor Cvv_t = build_and_fill("C", {nv, nv}, c2);
619
620 Avv("ij") = Avv_t("ij");
621 Cvv("ij") = Cvv_t("ij");
622
623 C("ab") = A("ab");
624
625 for (size_t a = 0; a < nv; ++a)
626 {
627 for (size_t b = 0; b < nv; ++b)
628 {
629 c2[a][b] = a2[a][b];
630 }
631 }
632
633 double diff_vv = difference(Cvv, c2).second;
634
635 return diff_oo + diff_vo + diff_vv;
636 }
637
test_Aij_equal_Aji()638 double test_Aij_equal_Aji()
639 {
640 BlockedTensor::reset_mo_spaces();
641 BlockedTensor::add_mo_space("o", "i,j", {0, 1, 2}, AlphaSpin);
642 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
643
644 BlockedTensor A =
645 BlockedTensor::build(CoreTensor, "A", {"oo", "vv", "ov", "vo"});
646
647 A("ij") = A("ji");
648
649 return 0.0;
650 }
651
test_Cijab_plus_equal_Aaibj()652 double test_Cijab_plus_equal_Aaibj()
653 {
654 BlockedTensor::reset_mo_spaces();
655 BlockedTensor::add_mo_space("o", "i,j", {0, 1, 2}, AlphaSpin);
656 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
657
658 BlockedTensor A =
659 BlockedTensor::build(CoreTensor, "A", {"vovo", "ovvo", "voov"});
660 BlockedTensor C =
661 BlockedTensor::build(CoreTensor, "C", {"oovv", "ovvo", "voov"});
662
663 size_t no = 3;
664 size_t nv = 5;
665
666 Tensor Avovo = A.block("vovo");
667 Tensor Coovv = C.block("oovv");
668
669 Tensor Avovo_t = build_and_fill("A", {nv, no, nv, no}, a4);
670 Tensor Coovv_t = build_and_fill("C", {no, no, nv, nv}, c4);
671
672 Avovo("pqrs") = Avovo_t("pqrs");
673 Coovv("pqrs") = Coovv_t("pqrs");
674
675 C("ijab") += A("aibj");
676
677 for (size_t i = 0; i < no; ++i)
678 {
679 for (size_t j = 0; j < no; ++j)
680 {
681 for (size_t a = 0; a < nv; ++a)
682 {
683 for (size_t b = 0; b < nv; ++b)
684 {
685 c4[i][j][a][b] += a4[a][i][b][j];
686 }
687 }
688 }
689 }
690
691 return difference(Coovv, c4).second;
692 }
693
test_Cbija_minus_equal_Ajabi()694 double test_Cbija_minus_equal_Ajabi()
695 {
696 BlockedTensor::reset_mo_spaces();
697 BlockedTensor::add_mo_space("o", "i,j", {0, 1, 2}, AlphaSpin);
698 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
699
700 BlockedTensor A =
701 BlockedTensor::build(CoreTensor, "A", {"vovo", "ovvo", "voov"});
702 BlockedTensor C =
703 BlockedTensor::build(CoreTensor, "C", {"oovv", "ovvo", "voov"});
704
705 size_t no = 3;
706 size_t nv = 5;
707
708 Tensor Aovvo = A.block("ovvo");
709 Tensor Cvoov = C.block("voov");
710
711 Tensor Aovvo_t = build_and_fill("A", {no, nv, nv, no}, a4);
712 Tensor Cvoov_t = build_and_fill("C", {nv, no, no, nv}, c4);
713
714 Aovvo("pqrs") = Aovvo_t("pqrs");
715 Cvoov("pqrs") = Cvoov_t("pqrs");
716
717 C("bija") -= A("jabi");
718
719 for (size_t i = 0; i < no; ++i)
720 {
721 for (size_t j = 0; j < no; ++j)
722 {
723 for (size_t a = 0; a < nv; ++a)
724 {
725 for (size_t b = 0; b < nv; ++b)
726 {
727 c4[b][i][j][a] -= a4[j][a][b][i];
728 }
729 }
730 }
731 }
732
733 return difference(Cvoov, c4).second;
734 }
735
test_Cij_times_equal_double()736 double test_Cij_times_equal_double()
737 {
738 BlockedTensor::reset_mo_spaces();
739 BlockedTensor::add_mo_space("o", "i,j", {0, 1, 2}, AlphaSpin);
740 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
741
742 BlockedTensor A =
743 BlockedTensor::build(CoreTensor, "A", {"oo", "vv", "ov", "vo"});
744
745 size_t no = 3;
746 size_t nv = 5;
747
748 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
749 Tensor Aov_t = build_and_fill("Aov", {no, nv}, b2);
750
751 A.block("oo")("pq") = Aoo_t("pq");
752 A.block("ov")("pq") = Aov_t("pq");
753
754 A("ij") *= std::exp(1.0);
755
756 for (size_t i = 0; i < no; ++i)
757 {
758 for (size_t j = 0; j < no; ++j)
759 {
760 c2[i][j] = std::exp(1.0) * a2[i][j];
761 }
762 }
763
764 Tensor Aoo = A.block("oo");
765 double diff_oo = difference(Aoo, c2).second;
766
767 // A("ia") /= std::exp(1.0);
768
769 for (size_t i = 0; i < no; ++i)
770 {
771 for (size_t a = 0; a < nv; ++a)
772 {
773 c2[i][a] = b2[i][a];
774 }
775 }
776
777 Tensor Aov = A.block("ov");
778 double diff_ov = difference(Aov, c2).second;
779
780 return diff_oo + diff_ov;
781 }
782
test_Cip_times_equal_double()783 double test_Cip_times_equal_double()
784 {
785 BlockedTensor::reset_mo_spaces();
786 BlockedTensor::add_mo_space("o", "i,j", {0, 1, 2}, AlphaSpin);
787 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
788 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
789
790 BlockedTensor A =
791 BlockedTensor::build(CoreTensor, "A", {"oo", "vv", "ov", "vo"});
792
793 size_t no = 3;
794 size_t nv = 5;
795
796 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
797 Tensor Aov_t = build_and_fill("Aov", {no, nv}, b2);
798
799 A.block("oo")("pq") = Aoo_t("pq");
800 A.block("ov")("pq") = Aov_t("pq");
801
802 A("ip") *= std::exp(1.0);
803
804 for (size_t i = 0; i < no; ++i)
805 {
806 for (size_t j = 0; j < no; ++j)
807 {
808 c2[i][j] = std::exp(1.0) * a2[i][j];
809 }
810 }
811
812 Tensor Aoo = A.block("oo");
813 double diff_oo = difference(Aoo, c2).second;
814
815 for (size_t i = 0; i < no; ++i)
816 {
817 for (size_t a = 0; a < nv; ++a)
818 {
819 c2[i][a] = std::exp(1.0) * b2[i][a];
820 }
821 }
822
823 Tensor Aov = A.block("ov");
824 double diff_ov = difference(Aov, c2).second;
825
826 return diff_oo + diff_ov;
827 }
828
test_Cij_equal_Aik_B_jk()829 double test_Cij_equal_Aik_B_jk()
830 {
831 BlockedTensor::reset_mo_spaces();
832 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
833 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
834 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
835
836 BlockedTensor A =
837 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
838 BlockedTensor B =
839 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
840 BlockedTensor C =
841 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
842
843 size_t no = 3;
844
845 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
846 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
847 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
848
849 A.block("oo")("pq") = Aoo_t("pq");
850 B.block("oo")("pq") = Boo_t("pq");
851 C.block("oo")("pq") = Coo_t("pq");
852
853 for (size_t i = 0; i < no; ++i)
854 {
855 for (size_t j = 0; j < no; ++j)
856 {
857 c2[i][j] = 0.0;
858 for (size_t k = 0; k < no; ++k)
859 {
860 c2[i][j] += a2[i][k] * b2[j][k];
861 }
862 }
863 }
864
865 C("ij") = A("ik") * B("jk");
866
867 Tensor Coo = C.block("oo");
868 double diff_oo = difference(Coo, c2).second;
869
870 return diff_oo;
871 }
872
test_Cij_equal_Aip_B_jp()873 double test_Cij_equal_Aip_B_jp()
874 {
875 BlockedTensor::reset_mo_spaces();
876 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 10, 12}, AlphaSpin);
877 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9, 3, 4},
878 AlphaSpin);
879 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
880
881 BlockedTensor A =
882 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
883 BlockedTensor B =
884 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
885 BlockedTensor C =
886 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
887
888 size_t no = 5;
889 size_t nv = 7;
890
891 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
892 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
893 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
894 Tensor Aov_t = build_and_fill("Aov", {no, nv}, d2);
895 Tensor Bov_t = build_and_fill("Bov", {no, nv}, e2);
896 Tensor Cov_t = build_and_fill("Cov", {no, nv}, f2);
897
898 A.block("oo")("pq") = Aoo_t("pq");
899 B.block("oo")("pq") = Boo_t("pq");
900 C.block("oo")("pq") = Coo_t("pq");
901 A.block("ov")("pq") = Aov_t("pq");
902 B.block("ov")("pq") = Bov_t("pq");
903 C.block("ov")("pq") = Cov_t("pq");
904
905 for (size_t i = 0; i < no; ++i)
906 {
907 for (size_t j = 0; j < no; ++j)
908 {
909 c2[i][j] = 0.0;
910 for (size_t k = 0; k < no; ++k)
911 {
912 c2[i][j] += a2[i][k] * b2[j][k];
913 }
914 for (size_t a = 0; a < nv; ++a)
915 {
916 c2[i][j] += d2[i][a] * e2[j][a];
917 }
918 }
919 }
920
921 C("ij") = A("ip") * B("jp");
922 C("ab") = A("ap") * B("bp");
923
924 Tensor Coo = C.block("oo");
925 double diff_oo = difference(Coo, c2).second;
926
927 return diff_oo;
928 }
929
test_Cij_equal_Aip_B_jp_fail()930 double test_Cij_equal_Aip_B_jp_fail()
931 {
932 BlockedTensor::reset_mo_spaces();
933 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 10, 12}, AlphaSpin);
934 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9, 3, 4},
935 AlphaSpin);
936 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
937
938 BlockedTensor A = BlockedTensor::build(CoreTensor, "A", {"oo", "vo", "vv"});
939 BlockedTensor B = BlockedTensor::build(CoreTensor, "B", {"oo", "vo", "vv"});
940 BlockedTensor C =
941 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
942
943 size_t no = 5;
944 size_t nv = 7;
945
946 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
947 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
948 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
949 Tensor Cov_t = build_and_fill("Cov", {no, nv}, f2);
950
951 A.block("oo")("pq") = Aoo_t("pq");
952 B.block("oo")("pq") = Boo_t("pq");
953 C.block("oo")("pq") = Coo_t("pq");
954 C.block("ov")("pq") = Cov_t("pq");
955
956 for (size_t i = 0; i < no; ++i)
957 {
958 for (size_t j = 0; j < no; ++j)
959 {
960 c2[i][j] = 0.0;
961 for (size_t k = 0; k < no; ++k)
962 {
963 c2[i][j] += a2[i][k] * b2[j][k];
964 }
965 }
966 }
967
968 C("ij") = A("ip") * B("jp");
969
970 Tensor Coo = C.block("oo");
971 double diff_oo = difference(Coo, c2).second;
972
973 return diff_oo;
974 }
975
test_Cij_equal_Aip_B_jp_expert()976 double test_Cij_equal_Aip_B_jp_expert()
977 {
978 BlockedTensor::set_expert_mode(true);
979 BlockedTensor::reset_mo_spaces();
980 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 10, 12}, AlphaSpin);
981 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9, 3, 4},
982 AlphaSpin);
983 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
984
985 BlockedTensor A = BlockedTensor::build(CoreTensor, "A", {"oo", "vo", "vv"});
986 BlockedTensor B = BlockedTensor::build(CoreTensor, "B", {"oo", "vo", "vv"});
987 BlockedTensor C =
988 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
989
990 size_t no = 5;
991 size_t nv = 7;
992
993 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
994 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
995 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
996 Tensor Cov_t = build_and_fill("Cov", {no, nv}, f2);
997
998 A.block("oo")("pq") = Aoo_t("pq");
999 B.block("oo")("pq") = Boo_t("pq");
1000 C.block("oo")("pq") = Coo_t("pq");
1001 C.block("ov")("pq") = Cov_t("pq");
1002
1003 for (size_t i = 0; i < no; ++i)
1004 {
1005 for (size_t j = 0; j < no; ++j)
1006 {
1007 c2[i][j] = 0.0;
1008 for (size_t k = 0; k < no; ++k)
1009 {
1010 c2[i][j] += a2[i][k] * b2[j][k];
1011 }
1012 }
1013 }
1014
1015 C("ij") = A("ip") * B("jp");
1016
1017 Tensor Coo = C.block("oo");
1018 double diff_oo = difference(Coo, c2).second;
1019
1020 BlockedTensor::set_expert_mode(false);
1021
1022 return diff_oo;
1023 }
1024
test_Cpq_equal_Apq_B_pq_expert()1025 double test_Cpq_equal_Apq_B_pq_expert()
1026 {
1027 BlockedTensor::set_expert_mode(true);
1028 BlockedTensor::reset_mo_spaces();
1029 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 10, 12}, AlphaSpin);
1030 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9, 3, 4},
1031 AlphaSpin);
1032 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
1033
1034 BlockedTensor A = BlockedTensor::build(CoreTensor, "A", {"oo", "vo", "ov"});
1035 BlockedTensor B = BlockedTensor::build(CoreTensor, "B", {"oo", "vo", "ov"});
1036 BlockedTensor C = BlockedTensor::build(CoreTensor, "C", {"oo"});
1037
1038 size_t no = 5;
1039 size_t nv = 7;
1040
1041 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
1042 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
1043 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1044 Tensor Cov_t = build_and_fill("Cov", {no, nv}, f2);
1045
1046 A.block("oo")("pq") = Aoo_t("pq");
1047 B.block("oo")("pq") = Boo_t("pq");
1048 C.block("oo")("pq") = Coo_t("pq");
1049
1050 for (size_t i = 0; i < no; ++i)
1051 {
1052 for (size_t j = 0; j < no; ++j)
1053 {
1054 c2[i][j] = a2[i][j] * b2[i][j];
1055 }
1056 }
1057
1058 C("pq") = A("pq") * B("pq");
1059
1060 Tensor Coo = C.block("oo");
1061 double diff_oo = difference(Coo, c2).second;
1062
1063 BlockedTensor::set_expert_mode(false);
1064
1065 return diff_oo;
1066 }
1067
test_Cij_equal_half_Aia_B_aj()1068 double test_Cij_equal_half_Aia_B_aj()
1069 {
1070 BlockedTensor::reset_mo_spaces();
1071 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
1072 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
1073 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
1074
1075 BlockedTensor A =
1076 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1077 BlockedTensor B =
1078 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1079 BlockedTensor C =
1080 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1081
1082 size_t no = 3;
1083 size_t nv = 5;
1084
1085 Tensor Aov_t = build_and_fill("Aov", {no, nv}, a2);
1086 Tensor Bvo_t = build_and_fill("Bvo", {nv, no}, b2);
1087 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1088
1089 A.block("ov")("pq") = Aov_t("pq");
1090 B.block("vo")("pq") = Bvo_t("pq");
1091 C.block("oo")("pq") = Coo_t("pq");
1092
1093 for (size_t i = 0; i < no; ++i)
1094 {
1095 for (size_t j = 0; j < no; ++j)
1096 {
1097 c2[i][j] = 0.0;
1098 for (size_t a = 0; a < nv; ++a)
1099 {
1100 c2[i][j] += 0.5 * a2[i][a] * b2[a][j];
1101 }
1102 }
1103 }
1104
1105 C("ij") = 0.5 * A("ia") * B("aj");
1106
1107 Tensor Coo = C.block("oo");
1108 double diff_oo = difference(Coo, c2).second;
1109
1110 return diff_oo;
1111 }
1112
test_Cij_plus_equal_half_Aai_B_ja()1113 double test_Cij_plus_equal_half_Aai_B_ja()
1114 {
1115 BlockedTensor::reset_mo_spaces();
1116 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
1117 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
1118 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
1119
1120 BlockedTensor A =
1121 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1122 BlockedTensor B =
1123 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1124 BlockedTensor C =
1125 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1126
1127 size_t no = 3;
1128 size_t nv = 5;
1129
1130 Tensor Avo_t = build_and_fill("Avo", {nv, no}, a2);
1131 Tensor Bov_t = build_and_fill("Bov", {no, nv}, b2);
1132 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1133
1134 A.block("vo")("pq") = Avo_t("pq");
1135 B.block("ov")("pq") = Bov_t("pq");
1136 C.block("oo")("pq") = Coo_t("pq");
1137
1138 for (size_t i = 0; i < no; ++i)
1139 {
1140 for (size_t j = 0; j < no; ++j)
1141 {
1142 for (size_t a = 0; a < nv; ++a)
1143 {
1144 c2[i][j] += 0.5 * a2[a][i] * b2[j][a];
1145 }
1146 }
1147 }
1148
1149 C("ij") += A("ai") * (0.5 * B("ja"));
1150
1151 Tensor Coo = C.block("oo");
1152 double diff_oo = difference(Coo, c2).second;
1153
1154 return diff_oo;
1155 }
1156
test_Cij_minus_equal_Aik_B_jk()1157 double test_Cij_minus_equal_Aik_B_jk()
1158 {
1159 BlockedTensor::reset_mo_spaces();
1160 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
1161 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
1162 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
1163
1164 BlockedTensor A =
1165 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1166 BlockedTensor B =
1167 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1168 BlockedTensor C =
1169 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1170
1171 size_t no = 3;
1172
1173 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
1174 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
1175 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1176
1177 A.block("oo")("pq") = Aoo_t("pq");
1178 B.block("oo")("pq") = Boo_t("pq");
1179 C.block("oo")("pq") = Coo_t("pq");
1180
1181 for (size_t i = 0; i < no; ++i)
1182 {
1183 for (size_t j = 0; j < no; ++j)
1184 {
1185 for (size_t k = 0; k < no; ++k)
1186 {
1187 c2[i][j] -= a2[i][k] * b2[j][k];
1188 }
1189 }
1190 }
1191
1192 C("ij") -= A("ik") * B("jk");
1193
1194 Tensor Coo = C.block("oo");
1195 double diff_oo = difference(Coo, c2).second;
1196
1197 return diff_oo;
1198 }
1199
test_greek_Cij_equal_Aik_B_jk()1200 double test_greek_Cij_equal_Aik_B_jk()
1201 {
1202 BlockedTensor::reset_mo_spaces();
1203 BlockedTensor::add_mo_space("o", "ι,κ,λ,μ,i,j,k,l", {0, 1, 2}, AlphaSpin);
1204 BlockedTensor::add_mo_space("v", "α,β,γ,δ", {5, 6, 7, 8, 9}, AlphaSpin);
1205 BlockedTensor::add_composite_mo_space("g", "ρ,σ", {"o", "v"});
1206
1207 BlockedTensor A =
1208 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1209 BlockedTensor B =
1210 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1211 BlockedTensor C =
1212 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1213
1214 size_t no = 3;
1215
1216 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
1217 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
1218 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1219
1220 A.block("oo")("pq") = Aoo_t("pq");
1221 B.block("oo")("pq") = Boo_t("pq");
1222 C.block("oo")("pq") = Coo_t("pq");
1223
1224 for (size_t i = 0; i < no; ++i)
1225 {
1226 for (size_t j = 0; j < no; ++j)
1227 {
1228 c2[i][j] = 0.0;
1229 for (size_t k = 0; k < no; ++k)
1230 {
1231 c2[i][j] += a2[i][k] * b2[j][k];
1232 }
1233 }
1234 }
1235
1236 C("ι,κ") = A("ι,λ") * B("κ,λ");
1237
1238 Tensor Coo = C.block("oo");
1239 double diff_oo = difference(Coo, c2).second;
1240
1241 return diff_oo;
1242 }
1243
test_Aij_equal_Aik_B_jk()1244 double test_Aij_equal_Aik_B_jk()
1245 {
1246 BlockedTensor::reset_mo_spaces();
1247 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
1248 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
1249
1250 BlockedTensor A =
1251 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1252 BlockedTensor B =
1253 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1254
1255 A("ij") = A("ik") * B("jk");
1256
1257 return 0.0;
1258 }
1259
test_chain_multiply()1260 double test_chain_multiply()
1261 {
1262 BlockedTensor::reset_mo_spaces();
1263 BlockedTensor::add_mo_space("o", "i,j,k,l", {0, 1, 2}, AlphaSpin);
1264 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
1265 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
1266
1267 BlockedTensor A =
1268 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1269 BlockedTensor B =
1270 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1271 BlockedTensor C =
1272 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1273 BlockedTensor D =
1274 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1275
1276 size_t no = 3;
1277
1278 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
1279 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
1280 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1281 Tensor Doo_t = build_and_fill("Doo", {no, no}, d2);
1282
1283 A.block("oo")("pq") = Aoo_t("pq");
1284 B.block("oo")("pq") = Boo_t("pq");
1285 C.block("oo")("pq") = Coo_t("pq");
1286 D.block("oo")("pq") = Doo_t("pq");
1287
1288 for (size_t i = 0; i < no; ++i)
1289 {
1290 for (size_t j = 0; j < no; ++j)
1291 {
1292 d2[i][j] = 0.0;
1293 for (size_t k = 0; k < no; ++k)
1294 {
1295 for (size_t l = 0; l < no; ++l)
1296 {
1297 d2[i][j] += a2[l][j] * b2[i][k] * c2[k][l];
1298 }
1299 }
1300 }
1301 }
1302
1303 D("ij") = B("ik") * C("kl") * A("lj");
1304
1305 Tensor Doo = D.block("oo");
1306 return difference(Doo, d2).second;
1307 }
1308
test_chain_multiply2()1309 double test_chain_multiply2()
1310 {
1311 BlockedTensor::reset_mo_spaces();
1312 BlockedTensor::add_mo_space("o", "i,j,k,l", {0, 1, 2, 3, 4}, AlphaSpin);
1313 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9, 10, 15, 20},
1314 AlphaSpin);
1315 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
1316
1317 BlockedTensor A = BlockedTensor::build(CoreTensor, "A", {"vvoo"});
1318 BlockedTensor B =
1319 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1320 BlockedTensor C =
1321 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1322 BlockedTensor D = BlockedTensor::build(CoreTensor, "D", {"oovv", "ovvo"});
1323
1324 size_t no = 5;
1325 size_t nv = 8;
1326
1327 Tensor Avvoo_t = build_and_fill("Aoo", {nv, nv, no, no}, a4);
1328 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
1329 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1330 Tensor Doovv_t = build_and_fill("Doo", {no, no, nv, nv}, d4);
1331
1332 A.block("vvoo")("pqrs") = Avvoo_t("pqrs");
1333 B.block("oo")("pq") = Boo_t("pq");
1334 C.block("oo")("pq") = Coo_t("pq");
1335 D.block("oovv")("pqrs") = Doovv_t("pqrs");
1336
1337 for (size_t i = 0; i < no; ++i)
1338 {
1339 for (size_t j = 0; j < no; ++j)
1340 {
1341 for (size_t a = 0; a < nv; ++a)
1342 {
1343 for (size_t b = 0; b < nv; ++b)
1344 {
1345 d4[i][j][a][b] = 0.0;
1346 for (size_t k = 0; k < no; ++k)
1347 {
1348 for (size_t l = 0; l < no; ++l)
1349 {
1350 d4[i][j][a][b] +=
1351 a4[a][b][l][j] * b2[i][k] * c2[k][l];
1352 }
1353 }
1354 }
1355 }
1356 }
1357 }
1358
1359 D("ijab") = B("ik") * C("kl") * A("ablj");
1360
1361 Tensor Doo = D.block("oovv");
1362 return difference(Doo, d4).second;
1363 }
1364
test_Cij_equal_Aij_plus_Bij()1365 double test_Cij_equal_Aij_plus_Bij()
1366 {
1367 BlockedTensor::reset_mo_spaces();
1368 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
1369 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
1370
1371 BlockedTensor A =
1372 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1373 BlockedTensor B =
1374 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1375 BlockedTensor C =
1376 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1377
1378 size_t no = 3;
1379
1380 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
1381 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
1382 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1383
1384 A.block("oo")("pq") = Aoo_t("pq");
1385 B.block("oo")("pq") = Boo_t("pq");
1386 C.block("oo")("pq") = Coo_t("pq");
1387
1388 for (size_t i = 0; i < no; ++i)
1389 {
1390 for (size_t j = 0; j < no; ++j)
1391 {
1392 c2[i][j] = a2[i][j] + b2[i][j];
1393 }
1394 }
1395
1396 C("ij") = A("ij") + B("ij");
1397
1398 Tensor Coo = C.block("oo");
1399 return difference(Coo, c2).second;
1400 }
1401
test_Cia_plus_equal_Aia_minus_three_Bai()1402 double test_Cia_plus_equal_Aia_minus_three_Bai()
1403 {
1404 BlockedTensor::reset_mo_spaces();
1405 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
1406 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
1407 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
1408
1409 BlockedTensor A =
1410 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1411 BlockedTensor B =
1412 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1413 BlockedTensor C =
1414 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1415
1416 size_t no = 3;
1417 size_t nv = 5;
1418
1419 Tensor Aov_t = build_and_fill("Aoo", {no, nv}, a2);
1420 Tensor Bvo_t = build_and_fill("Boo", {nv, no}, b2);
1421 Tensor Cov_t = build_and_fill("Coo", {no, nv}, c2);
1422
1423 A.block("ov")("pq") = Aov_t("pq");
1424 B.block("vo")("pq") = Bvo_t("pq");
1425 C.block("ov")("pq") = Cov_t("pq");
1426
1427 for (size_t i = 0; i < no; ++i)
1428 {
1429 for (size_t a = 0; a < nv; ++a)
1430 {
1431 c2[i][a] += a2[i][a] - 3.0 * b2[a][i];
1432 }
1433 }
1434
1435 C("ia") += A("ia") - 3.0 * B("ai");
1436
1437 Tensor Cov = C.block("ov");
1438 return difference(Cov, c2).second;
1439 }
1440
test_Dij_equal_Aij_times_Bij_plus_Cij()1441 double test_Dij_equal_Aij_times_Bij_plus_Cij()
1442 {
1443 BlockedTensor::reset_mo_spaces();
1444 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 4, 5}, AlphaSpin);
1445 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
1446
1447 BlockedTensor A =
1448 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1449 BlockedTensor B =
1450 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1451 BlockedTensor C =
1452 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1453 BlockedTensor D =
1454 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1455
1456 size_t no = 5;
1457
1458 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
1459 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
1460 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1461 Tensor Doo_t = build_and_fill("Doo", {no, no}, d2);
1462
1463 A.block("oo")("pq") = Aoo_t("pq");
1464 B.block("oo")("pq") = Boo_t("pq");
1465 C.block("oo")("pq") = Coo_t("pq");
1466 D.block("oo")("pq") = Doo_t("pq");
1467
1468 for (size_t i = 0; i < no; ++i)
1469 {
1470 for (size_t j = 0; j < no; ++j)
1471 {
1472 d2[i][j] = a2[i][j] * (2.0 * b2[i][j] - c2[i][j]);
1473 }
1474 }
1475
1476 D("ij") = A("ij") * (2.0 * B("ij") - C("ij"));
1477
1478 Tensor Doo = D.block("oo");
1479 return difference(Doo, d2).second;
1480 }
1481
test_Dij_equal_Bij_plus_Cij_times_Aij()1482 double test_Dij_equal_Bij_plus_Cij_times_Aij()
1483 {
1484 BlockedTensor::reset_mo_spaces();
1485 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 4, 5}, AlphaSpin);
1486 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
1487
1488 BlockedTensor A =
1489 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1490 BlockedTensor B =
1491 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1492 BlockedTensor C =
1493 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1494 BlockedTensor D =
1495 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1496
1497 size_t no = 5;
1498
1499 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
1500 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
1501 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1502 Tensor Doo_t = build_and_fill("Doo", {no, no}, d2);
1503
1504 A.block("oo")("pq") = Aoo_t("pq");
1505 B.block("oo")("pq") = Boo_t("pq");
1506 C.block("oo")("pq") = Coo_t("pq");
1507 D.block("oo")("pq") = Doo_t("pq");
1508
1509 for (size_t i = 0; i < no; ++i)
1510 {
1511 for (size_t j = 0; j < no; ++j)
1512 {
1513 d2[i][j] = a2[i][j] * (2.0 * b2[i][j] - c2[i][j]);
1514 }
1515 }
1516
1517 D("ij") = (2.0 * B("ij") - C("ij")) * A("ij");
1518
1519 Tensor Doo = D.block("oo");
1520 return difference(Doo, d2).second;
1521 }
1522
test_Dij_plus_equal_Bij_plus_Cij_times_Aij()1523 double test_Dij_plus_equal_Bij_plus_Cij_times_Aij()
1524 {
1525 BlockedTensor::reset_mo_spaces();
1526 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 4, 5}, AlphaSpin);
1527 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
1528
1529 BlockedTensor A =
1530 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1531 BlockedTensor B =
1532 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1533 BlockedTensor C =
1534 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1535 BlockedTensor D =
1536 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1537
1538 size_t no = 5;
1539
1540 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
1541 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
1542 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1543 Tensor Doo_t = build_and_fill("Doo", {no, no}, d2);
1544
1545 A.block("oo")("pq") = Aoo_t("pq");
1546 B.block("oo")("pq") = Boo_t("pq");
1547 C.block("oo")("pq") = Coo_t("pq");
1548 D.block("oo")("pq") = Doo_t("pq");
1549
1550 for (size_t i = 0; i < no; ++i)
1551 {
1552 for (size_t j = 0; j < no; ++j)
1553 {
1554 d2[i][j] += a2[i][j] * (2.0 * b2[i][j] - c2[i][j]);
1555 }
1556 }
1557
1558 D("ij") += (2.0 * B("ij") - C("ij")) * A("ij");
1559
1560 Tensor Doo = D.block("oo");
1561 return difference(Doo, d2).second;
1562 }
1563
test_F_equal_D_times_2g_minus_g()1564 double test_F_equal_D_times_2g_minus_g()
1565 {
1566 BlockedTensor::reset_mo_spaces();
1567 BlockedTensor::add_mo_space("o", "i,j,k,l", {0, 1, 2, 4, 5}, AlphaSpin);
1568 BlockedTensor::add_mo_space("v", "a,b,c,d", {6, 7, 8, 9}, AlphaSpin);
1569
1570 BlockedTensor F =
1571 BlockedTensor::build(CoreTensor, "F", {"oo", "ov", "vo", "vv"});
1572 BlockedTensor D =
1573 BlockedTensor::build(CoreTensor, "D", {"oo", "ov", "vo", "vv"});
1574 BlockedTensor g = BlockedTensor::build(CoreTensor, "g", {"oooo", "vvvv"});
1575
1576 size_t no = 5;
1577
1578 Tensor Foo_t = build_and_fill("Foo", {no, no}, a2);
1579 Tensor Doo_t = build_and_fill("Doo", {no, no}, b2);
1580 Tensor goo_t = build_and_fill("goo", {no, no, no, no}, c4);
1581
1582 F.block("oo")("pq") = Foo_t("pq");
1583 D.block("oo")("pq") = Doo_t("pq");
1584 g.block("oooo")("pqrs") = goo_t("pqrs");
1585
1586 for (size_t i = 0; i < no; ++i)
1587 {
1588 for (size_t j = 0; j < no; ++j)
1589 {
1590 a2[i][j] = 0.0;
1591 for (size_t k = 0; k < no; ++k)
1592 {
1593 for (size_t l = 0; l < no; ++l)
1594 {
1595 a2[i][j] +=
1596 b2[k][l] * (2.0 * c4[i][j][k][l] - c4[i][k][j][l]);
1597 }
1598 }
1599 }
1600 }
1601
1602 F("i,j") = D("k,l") * (2.0 * g("i,j,k,l") - g("i,k,j,l"));
1603 F("c,d") = D("a,b") * (2.0 * g("a,b,c,d") - g("a,c,b,d"));
1604
1605 Tensor Foo = F.block("oo");
1606 return difference(Foo, a2).second;
1607 }
1608
test_Dij_equal_2_times_Aij_plus_Bij()1609 double test_Dij_equal_2_times_Aij_plus_Bij()
1610 {
1611 BlockedTensor::reset_mo_spaces();
1612 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
1613 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
1614
1615 BlockedTensor A =
1616 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1617 BlockedTensor B =
1618 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1619 BlockedTensor C =
1620 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1621
1622 size_t no = 3;
1623
1624 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
1625 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
1626 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1627
1628 A.block("oo")("pq") = Aoo_t("pq");
1629 B.block("oo")("pq") = Boo_t("pq");
1630 C.block("oo")("pq") = Coo_t("pq");
1631
1632 for (size_t i = 0; i < no; ++i)
1633 {
1634 for (size_t j = 0; j < no; ++j)
1635 {
1636 c2[i][j] = 2.0 * (a2[i][j] - b2[i][j]);
1637 }
1638 }
1639
1640 C("ij") = 2.0 * (A("ij") - B("ij"));
1641
1642 Tensor Coo = C.block("oo");
1643 return difference(Coo, c2).second;
1644 }
1645
test_Cij_minus_equal_3_times_Aij_minus_2_Bji()1646 double test_Cij_minus_equal_3_times_Aij_minus_2_Bji()
1647 {
1648 BlockedTensor::reset_mo_spaces();
1649 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
1650 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
1651
1652 BlockedTensor A =
1653 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1654 BlockedTensor B =
1655 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1656 BlockedTensor C =
1657 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1658
1659 size_t no = 3;
1660
1661 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
1662 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
1663 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1664
1665 A.block("oo")("pq") = Aoo_t("pq");
1666 B.block("oo")("pq") = Boo_t("pq");
1667 C.block("oo")("pq") = Coo_t("pq");
1668
1669 for (size_t i = 0; i < no; ++i)
1670 {
1671 for (size_t j = 0; j < no; ++j)
1672 {
1673 c2[i][j] -= 3.0 * (a2[i][j] - 2.0 * b2[j][i]);
1674 }
1675 }
1676
1677 C("ij") -= 3.0 * (A("ij") - 2.0 * B("ji"));
1678
1679 Tensor Coo = C.block("oo");
1680 return difference(Coo, c2).second;
1681 }
1682
test_Cij_equal_negate_Aij_plus_Bij()1683 double test_Cij_equal_negate_Aij_plus_Bij()
1684 {
1685 BlockedTensor::reset_mo_spaces();
1686 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
1687 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
1688
1689 BlockedTensor A =
1690 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1691 BlockedTensor B =
1692 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1693 BlockedTensor C =
1694 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1695
1696 size_t no = 3;
1697
1698 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
1699 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
1700 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1701
1702 A.block("oo")("pq") = Aoo_t("pq");
1703 B.block("oo")("pq") = Boo_t("pq");
1704 C.block("oo")("pq") = Coo_t("pq");
1705
1706 for (size_t i = 0; i < no; ++i)
1707 {
1708 for (size_t j = 0; j < no; ++j)
1709 {
1710 c2[i][j] = -(a2[i][j] + b2[i][j]);
1711 }
1712 }
1713
1714 C("ij") = -(A("ij") + B("ij"));
1715
1716 Tensor Coo = C.block("oo");
1717 return difference(Coo, c2).second;
1718 }
1719
test_dot_product()1720 double test_dot_product()
1721 {
1722 BlockedTensor::reset_mo_spaces();
1723 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 3, 4}, AlphaSpin);
1724 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9, 10, 11},
1725 AlphaSpin);
1726 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
1727
1728 BlockedTensor A =
1729 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1730 BlockedTensor B =
1731 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1732
1733 size_t no = 5;
1734 size_t nv = 7;
1735
1736 Tensor Aov_t = build_and_fill("Aov", {no, nv}, a2);
1737 Tensor Bvo_t = build_and_fill("Bvo", {nv, no}, b2);
1738
1739 A.block("ov")("pq") = Aov_t("pq");
1740 B.block("vo")("pq") = Bvo_t("pq");
1741
1742 double c = 0.0;
1743 for (size_t i = 0; i < no; ++i)
1744 {
1745 for (size_t a = 0; a < nv; ++a)
1746 {
1747 c += a2[i][a] * b2[a][i];
1748 }
1749 }
1750
1751 double C = A("ia") * B("ai");
1752
1753 return std::fabs(C - c);
1754 }
1755
test_dot_product_fail1()1756 double test_dot_product_fail1()
1757 {
1758 BlockedTensor::reset_mo_spaces();
1759 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 3, 4}, AlphaSpin);
1760 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9, 10, 11},
1761 AlphaSpin);
1762 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
1763
1764 BlockedTensor A =
1765 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1766 BlockedTensor B =
1767 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1768
1769 size_t no = 5;
1770 size_t nv = 7;
1771
1772 Tensor Aov_t = build_and_fill("Aov", {no, nv}, a2);
1773 Tensor Bvo_t = build_and_fill("Bvo", {nv, no}, b2);
1774
1775 A.block("ov")("pq") = Aov_t("pq");
1776 B.block("vo")("pq") = Bvo_t("pq");
1777
1778 double c = 0.0;
1779 for (size_t i = 0; i < no; ++i)
1780 {
1781 for (size_t a = 0; a < nv; ++a)
1782 {
1783 c += a2[i][a] * b2[a][i];
1784 }
1785 }
1786
1787 double C = A("ia") * B("bi");
1788
1789 return std::fabs(C - c);
1790 }
1791
test_dot_product_fail2()1792 double test_dot_product_fail2()
1793 {
1794 BlockedTensor::reset_mo_spaces();
1795 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 3, 4}, AlphaSpin);
1796 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9, 10, 11},
1797 AlphaSpin);
1798 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
1799
1800 BlockedTensor A =
1801 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1802 BlockedTensor B =
1803 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1804
1805 size_t no = 5;
1806 size_t nv = 7;
1807
1808 Tensor Aov_t = build_and_fill("Aov", {no, nv}, a2);
1809 Tensor Bvo_t = build_and_fill("Bvo", {nv, no}, b2);
1810
1811 A.block("ov")("pq") = Aov_t("pq");
1812 B.block("vo")("pq") = Bvo_t("pq");
1813
1814 double c = 0.0;
1815 for (size_t i = 0; i < no; ++i)
1816 {
1817 for (size_t a = 0; a < nv; ++a)
1818 {
1819 c += a2[i][a] * b2[a][i];
1820 }
1821 }
1822
1823 double C = A("ia") * B("aij");
1824
1825 return std::fabs(C - c);
1826 }
1827
test_contraction_exception1()1828 double test_contraction_exception1()
1829 {
1830 BlockedTensor::reset_mo_spaces();
1831 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
1832 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
1833 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
1834
1835 BlockedTensor A =
1836 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1837 BlockedTensor B =
1838 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1839 BlockedTensor C = BlockedTensor::build(CoreTensor, "C", {"ov", "vo", "vv"});
1840
1841 C("ij") = A("ia") * B("aj");
1842
1843 return 0.0;
1844 }
1845
test_Cij_equal_Aik_B_jk_batched()1846 double test_Cij_equal_Aik_B_jk_batched()
1847 {
1848 BlockedTensor::reset_mo_spaces();
1849 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
1850 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
1851 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
1852
1853 BlockedTensor A =
1854 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1855 BlockedTensor B =
1856 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1857 BlockedTensor C =
1858 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1859
1860 size_t no = 3;
1861
1862 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
1863 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
1864 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1865
1866 A.block("oo")("pq") = Aoo_t("pq");
1867 B.block("oo")("pq") = Boo_t("pq");
1868 C.block("oo")("pq") = Coo_t("pq");
1869
1870 for (size_t i = 0; i < no; ++i)
1871 {
1872 for (size_t j = 0; j < no; ++j)
1873 {
1874 c2[i][j] = 0.0;
1875 for (size_t k = 0; k < no; ++k)
1876 {
1877 c2[i][j] += a2[i][k] * b2[j][k];
1878 }
1879 }
1880 }
1881
1882 C("ij") = batched("i", A("ik") * B("jk"));
1883
1884 Tensor Coo = C.block("oo");
1885 double diff_oo = difference(Coo, c2).second;
1886
1887 return diff_oo;
1888 }
1889
test_Cij_equal_Aip_B_jp_batched()1890 double test_Cij_equal_Aip_B_jp_batched()
1891 {
1892 BlockedTensor::reset_mo_spaces();
1893 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 10, 12}, AlphaSpin);
1894 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9, 3, 4},
1895 AlphaSpin);
1896 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
1897
1898 BlockedTensor A =
1899 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
1900 BlockedTensor B =
1901 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
1902 BlockedTensor C =
1903 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1904
1905 size_t no = 5;
1906 size_t nv = 7;
1907
1908 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
1909 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
1910 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1911 Tensor Aov_t = build_and_fill("Aov", {no, nv}, d2);
1912 Tensor Bov_t = build_and_fill("Bov", {no, nv}, e2);
1913 Tensor Cov_t = build_and_fill("Cov", {no, nv}, f2);
1914
1915 A.block("oo")("pq") = Aoo_t("pq");
1916 B.block("oo")("pq") = Boo_t("pq");
1917 C.block("oo")("pq") = Coo_t("pq");
1918 A.block("ov")("pq") = Aov_t("pq");
1919 B.block("ov")("pq") = Bov_t("pq");
1920 C.block("ov")("pq") = Cov_t("pq");
1921
1922 for (size_t i = 0; i < no; ++i)
1923 {
1924 for (size_t j = 0; j < no; ++j)
1925 {
1926 c2[i][j] = 0.0;
1927 for (size_t k = 0; k < no; ++k)
1928 {
1929 c2[i][j] += a2[i][k] * b2[j][k];
1930 }
1931 for (size_t a = 0; a < nv; ++a)
1932 {
1933 c2[i][j] += d2[i][a] * e2[j][a];
1934 }
1935 }
1936 }
1937
1938 C("ij") = batched("j", A("ip") * B("jp"));
1939 C("ab") = batched("b", A("ap") * B("bp"));
1940
1941 Tensor Coo = C.block("oo");
1942 double diff_oo = difference(Coo, c2).second;
1943
1944 return diff_oo;
1945 }
1946
test_Cij_equal_Aip_B_jp_fail_batched()1947 double test_Cij_equal_Aip_B_jp_fail_batched()
1948 {
1949 BlockedTensor::reset_mo_spaces();
1950 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 10, 12}, AlphaSpin);
1951 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9, 3, 4},
1952 AlphaSpin);
1953 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
1954
1955 BlockedTensor A = BlockedTensor::build(CoreTensor, "A", {"oo", "vo", "vv"});
1956 BlockedTensor B = BlockedTensor::build(CoreTensor, "B", {"oo", "vo", "vv"});
1957 BlockedTensor C =
1958 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
1959
1960 size_t no = 5;
1961 size_t nv = 7;
1962
1963 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
1964 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
1965 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
1966 Tensor Cov_t = build_and_fill("Cov", {no, nv}, f2);
1967
1968 A.block("oo")("pq") = Aoo_t("pq");
1969 B.block("oo")("pq") = Boo_t("pq");
1970 C.block("oo")("pq") = Coo_t("pq");
1971 C.block("ov")("pq") = Cov_t("pq");
1972
1973 for (size_t i = 0; i < no; ++i)
1974 {
1975 for (size_t j = 0; j < no; ++j)
1976 {
1977 c2[i][j] = 0.0;
1978 for (size_t k = 0; k < no; ++k)
1979 {
1980 c2[i][j] += a2[i][k] * b2[j][k];
1981 }
1982 }
1983 }
1984
1985 C("ij") = batched("i", A("ip") * B("jp"));
1986
1987 Tensor Coo = C.block("oo");
1988 double diff_oo = difference(Coo, c2).second;
1989
1990 return diff_oo;
1991 }
1992
test_Cij_equal_Aip_B_jp_expert_batched()1993 double test_Cij_equal_Aip_B_jp_expert_batched()
1994 {
1995 BlockedTensor::set_expert_mode(true);
1996 BlockedTensor::reset_mo_spaces();
1997 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 10, 12}, AlphaSpin);
1998 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9, 3, 4},
1999 AlphaSpin);
2000 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
2001
2002 BlockedTensor A = BlockedTensor::build(CoreTensor, "A", {"oo", "vo", "vv"});
2003 BlockedTensor B = BlockedTensor::build(CoreTensor, "B", {"oo", "vo", "vv"});
2004 BlockedTensor C =
2005 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
2006
2007 size_t no = 5;
2008 size_t nv = 7;
2009
2010 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
2011 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
2012 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
2013 Tensor Cov_t = build_and_fill("Cov", {no, nv}, f2);
2014
2015 A.block("oo")("pq") = Aoo_t("pq");
2016 B.block("oo")("pq") = Boo_t("pq");
2017 C.block("oo")("pq") = Coo_t("pq");
2018 C.block("ov")("pq") = Cov_t("pq");
2019
2020 for (size_t i = 0; i < no; ++i)
2021 {
2022 for (size_t j = 0; j < no; ++j)
2023 {
2024 c2[i][j] = 0.0;
2025 for (size_t k = 0; k < no; ++k)
2026 {
2027 c2[i][j] += a2[i][k] * b2[j][k];
2028 }
2029 }
2030 }
2031
2032 C("ij") = batched("i", A("ip") * B("jp"));
2033
2034 Tensor Coo = C.block("oo");
2035 double diff_oo = difference(Coo, c2).second;
2036
2037 BlockedTensor::set_expert_mode(false);
2038
2039 return diff_oo;
2040 }
2041
test_Cpq_equal_Apq_B_pq_expert_batched()2042 double test_Cpq_equal_Apq_B_pq_expert_batched()
2043 {
2044 BlockedTensor::set_expert_mode(true);
2045 BlockedTensor::reset_mo_spaces();
2046 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 10, 12}, AlphaSpin);
2047 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9, 3, 4},
2048 AlphaSpin);
2049 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
2050
2051 BlockedTensor A = BlockedTensor::build(CoreTensor, "A", {"oo", "vo", "ov"});
2052 BlockedTensor B = BlockedTensor::build(CoreTensor, "B", {"oo", "vo", "ov"});
2053 BlockedTensor C = BlockedTensor::build(CoreTensor, "C", {"oo"});
2054
2055 size_t no = 5;
2056 size_t nv = 7;
2057
2058 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
2059 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
2060 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
2061 Tensor Cov_t = build_and_fill("Cov", {no, nv}, f2);
2062
2063 A.block("oo")("pq") = Aoo_t("pq");
2064 B.block("oo")("pq") = Boo_t("pq");
2065 C.block("oo")("pq") = Coo_t("pq");
2066
2067 for (size_t i = 0; i < no; ++i)
2068 {
2069 for (size_t j = 0; j < no; ++j)
2070 {
2071 c2[i][j] = a2[i][j] * b2[i][j];
2072 }
2073 }
2074
2075 C("pq") = batched("p", A("pq") * B("pq"));
2076
2077 Tensor Coo = C.block("oo");
2078 double diff_oo = difference(Coo, c2).second;
2079
2080 BlockedTensor::set_expert_mode(false);
2081
2082 return diff_oo;
2083 }
2084
test_Cij_equal_half_Aia_B_aj_batched()2085 double test_Cij_equal_half_Aia_B_aj_batched()
2086 {
2087 BlockedTensor::reset_mo_spaces();
2088 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
2089 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
2090 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
2091
2092 BlockedTensor A =
2093 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
2094 BlockedTensor B =
2095 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
2096 BlockedTensor C =
2097 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
2098
2099 size_t no = 3;
2100 size_t nv = 5;
2101
2102 Tensor Aov_t = build_and_fill("Aov", {no, nv}, a2);
2103 Tensor Bvo_t = build_and_fill("Bvo", {nv, no}, b2);
2104 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
2105
2106 A.block("ov")("pq") = Aov_t("pq");
2107 B.block("vo")("pq") = Bvo_t("pq");
2108 C.block("oo")("pq") = Coo_t("pq");
2109
2110 for (size_t i = 0; i < no; ++i)
2111 {
2112 for (size_t j = 0; j < no; ++j)
2113 {
2114 c2[i][j] = 0.0;
2115 for (size_t a = 0; a < nv; ++a)
2116 {
2117 c2[i][j] += 0.5 * a2[i][a] * b2[a][j];
2118 }
2119 }
2120 }
2121
2122 C("ij") = batched("j", 0.5 * A("ia") * B("aj"));
2123
2124 Tensor Coo = C.block("oo");
2125 double diff_oo = difference(Coo, c2).second;
2126
2127 return diff_oo;
2128 }
2129
test_Cij_plus_equal_half_Aai_B_ja_batched()2130 double test_Cij_plus_equal_half_Aai_B_ja_batched()
2131 {
2132 BlockedTensor::reset_mo_spaces();
2133 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
2134 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
2135 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
2136
2137 BlockedTensor A =
2138 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
2139 BlockedTensor B =
2140 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
2141 BlockedTensor C =
2142 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
2143
2144 size_t no = 3;
2145 size_t nv = 5;
2146
2147 Tensor Avo_t = build_and_fill("Avo", {nv, no}, a2);
2148 Tensor Bov_t = build_and_fill("Bov", {no, nv}, b2);
2149 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
2150
2151 A.block("vo")("pq") = Avo_t("pq");
2152 B.block("ov")("pq") = Bov_t("pq");
2153 C.block("oo")("pq") = Coo_t("pq");
2154
2155 for (size_t i = 0; i < no; ++i)
2156 {
2157 for (size_t j = 0; j < no; ++j)
2158 {
2159 for (size_t a = 0; a < nv; ++a)
2160 {
2161 c2[i][j] += 0.5 * a2[a][i] * b2[j][a];
2162 }
2163 }
2164 }
2165
2166 C("ij") += batched("i", A("ai") * (0.5 * B("ja")));
2167
2168 Tensor Coo = C.block("oo");
2169 double diff_oo = difference(Coo, c2).second;
2170
2171 return diff_oo;
2172 }
2173
test_Cij_minus_equal_Aik_B_jk_batched()2174 double test_Cij_minus_equal_Aik_B_jk_batched()
2175 {
2176 BlockedTensor::reset_mo_spaces();
2177 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
2178 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
2179 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
2180
2181 BlockedTensor A =
2182 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
2183 BlockedTensor B =
2184 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
2185 BlockedTensor C =
2186 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
2187
2188 size_t no = 3;
2189
2190 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
2191 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
2192 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
2193
2194 A.block("oo")("pq") = Aoo_t("pq");
2195 B.block("oo")("pq") = Boo_t("pq");
2196 C.block("oo")("pq") = Coo_t("pq");
2197
2198 for (size_t i = 0; i < no; ++i)
2199 {
2200 for (size_t j = 0; j < no; ++j)
2201 {
2202 for (size_t k = 0; k < no; ++k)
2203 {
2204 c2[i][j] -= a2[i][k] * b2[j][k];
2205 }
2206 }
2207 }
2208
2209 C("ij") -= batched("i", A("ik") * B("jk"));
2210
2211 Tensor Coo = C.block("oo");
2212 double diff_oo = difference(Coo, c2).second;
2213
2214 return diff_oo;
2215 }
2216
test_greek_Cij_equal_Aik_B_jk_batched()2217 double test_greek_Cij_equal_Aik_B_jk_batched()
2218 {
2219 BlockedTensor::reset_mo_spaces();
2220 BlockedTensor::add_mo_space("o", "ι,κ,λ,μ,i,j,k,l", {0, 1, 2}, AlphaSpin);
2221 BlockedTensor::add_mo_space("v", "α,β,γ,δ", {5, 6, 7, 8, 9}, AlphaSpin);
2222 BlockedTensor::add_composite_mo_space("g", "ρ,σ", {"o", "v"});
2223
2224 BlockedTensor A =
2225 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
2226 BlockedTensor B =
2227 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
2228 BlockedTensor C =
2229 BlockedTensor::build(CoreTensor, "C", {"oo", "ov", "vo", "vv"});
2230
2231 size_t no = 3;
2232
2233 Tensor Aoo_t = build_and_fill("Aoo", {no, no}, a2);
2234 Tensor Boo_t = build_and_fill("Boo", {no, no}, b2);
2235 Tensor Coo_t = build_and_fill("Coo", {no, no}, c2);
2236
2237 A.block("oo")("pq") = Aoo_t("pq");
2238 B.block("oo")("pq") = Boo_t("pq");
2239 C.block("oo")("pq") = Coo_t("pq");
2240
2241 for (size_t i = 0; i < no; ++i)
2242 {
2243 for (size_t j = 0; j < no; ++j)
2244 {
2245 c2[i][j] = 0.0;
2246 for (size_t k = 0; k < no; ++k)
2247 {
2248 c2[i][j] += a2[i][k] * b2[j][k];
2249 }
2250 }
2251 }
2252
2253 C("ι,κ") = batched("κ,", A("ι,λ") * B("κ,λ"));
2254
2255 Tensor Coo = C.block("oo");
2256 double diff_oo = difference(Coo, c2).second;
2257
2258 return diff_oo;
2259 }
2260
test_Aij_equal_Aik_B_jk_batched()2261 double test_Aij_equal_Aik_B_jk_batched()
2262 {
2263 BlockedTensor::reset_mo_spaces();
2264 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2}, AlphaSpin);
2265 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9}, AlphaSpin);
2266
2267 BlockedTensor A =
2268 BlockedTensor::build(CoreTensor, "A", {"oo", "ov", "vo", "vv"});
2269 BlockedTensor B =
2270 BlockedTensor::build(CoreTensor, "B", {"oo", "ov", "vo", "vv"});
2271
2272 A("ij") = batched("j", A("ik") * B("jk"));
2273
2274 return 0.0;
2275 }
2276
test_batched_with_factor()2277 double test_batched_with_factor()
2278 {
2279 BlockedTensor::reset_mo_spaces();
2280 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 10, 12}, AlphaSpin);
2281 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9, 3, 4},
2282 AlphaSpin);
2283 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
2284
2285 BlockedTensor A =
2286 BlockedTensor::build(CoreTensor, "A", {"gggg"});
2287 BlockedTensor B =
2288 BlockedTensor::build(CoreTensor, "B", {"gggg"});
2289 BlockedTensor C =
2290 BlockedTensor::build(CoreTensor, "C", {"gggg"});
2291 BlockedTensor C2 =
2292 BlockedTensor::build(CoreTensor, "C2", {"gggg"});
2293
2294 size_t no = 5;
2295 size_t nv = 7;
2296
2297 for (const std::string & bl : A.block_labels()) {
2298 A.block(bl)("pqrs") = build_and_fill("A"+bl, A.block(bl).dims(), a4)("pqrs");
2299 }
2300
2301 for (const std::string & bl : B.block_labels()) {
2302 B.block(bl)("pqrs") = build_and_fill("B"+bl, B.block(bl).dims(), b4)("pqrs");
2303 }
2304
2305 for (const std::string & bl : C.block_labels()) {
2306 C.block(bl)("pqrs") = build_and_fill("C"+bl, C.block(bl).dims(), c4)("pqrs");
2307 }
2308
2309 C2["pqrs"] = C["pqrs"];
2310
2311 C["ijrs"] = 0.5 * A["abrs"] * B["ijab"];
2312 C2["ijrs"] = batched("r", 0.5 * A["abrs"] * B["ijab"]);
2313
2314 C2["pqrs"] -= C["pqrs"];
2315
2316 return C2.norm(0);
2317 }
2318
test_batch_all_indices()2319 double test_batch_all_indices()
2320 {
2321 BlockedTensor::reset_mo_spaces();
2322 BlockedTensor::add_mo_space("o", "i,j,k", {0, 1, 2, 10, 12}, AlphaSpin);
2323 BlockedTensor::add_mo_space("v", "a,b,c,d", {5, 6, 7, 8, 9, 3, 4},
2324 AlphaSpin);
2325 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"o", "v"});
2326
2327 BlockedTensor A =
2328 BlockedTensor::build(CoreTensor, "A", {"gggg"});
2329 BlockedTensor B =
2330 BlockedTensor::build(CoreTensor, "B", {"gggg"});
2331 BlockedTensor C =
2332 BlockedTensor::build(CoreTensor, "C", {"gggg"});
2333 BlockedTensor C2 =
2334 BlockedTensor::build(CoreTensor, "C2", {"gggg"});
2335
2336 size_t no = 5;
2337 size_t nv = 7;
2338
2339 for (const std::string & bl : A.block_labels()) {
2340 A.block(bl)("pqrs") = build_and_fill("A"+bl, A.block(bl).dims(), a4)("pqrs");
2341 }
2342
2343 for (const std::string & bl : B.block_labels()) {
2344 B.block(bl)("pqrs") = build_and_fill("B"+bl, B.block(bl).dims(), b4)("pqrs");
2345 }
2346
2347 for (const std::string & bl : C.block_labels()) {
2348 C.block(bl)("pqrs") = build_and_fill("C"+bl, C.block(bl).dims(), c4)("pqrs");
2349 }
2350
2351 C2["pqrs"] = C["pqrs"];
2352
2353 C["ijrs"] = 0.5 * A["abrs"] * B["ijab"];
2354 C2["ijrs"] = batched("srji", 0.5 * A["abrs"] * B["ijab"]);
2355
2356 C2["pqrs"] -= C["pqrs"];
2357
2358 return C2.norm(0);
2359 }
2360
test_Oia_equal_Cbu_Guv_Tivab_expert()2361 double test_Oia_equal_Cbu_Guv_Tivab_expert()
2362 {
2363 BlockedTensor::set_expert_mode(true);
2364 BlockedTensor::reset_mo_spaces();
2365
2366 const size_t nc = 3, na = 2, nv = 5;
2367
2368 // define space labels
2369 std::string acore_label_ = "c";
2370 std::string aactv_label_ = "a";
2371 std::string avirt_label_ = "v";
2372
2373 std::vector<size_t> core_mos_(nc);
2374 for (size_t i = 0; i < nc; ++i) {
2375 core_mos_[i] = i;
2376 }
2377 std::vector<size_t> actv_mos_(na);
2378 for (size_t i = 0; i < na; ++i) {
2379 actv_mos_[i] = i;
2380 }
2381 std::vector<size_t> virt_mos_(nv);
2382 for (size_t i = 0; i < nv; ++i) {
2383 virt_mos_[i] = i;
2384 }
2385 // add Ambit index labels
2386 BlockedTensor::add_mo_space(acore_label_, "mn", core_mos_, AlphaSpin);
2387 BlockedTensor::add_mo_space(aactv_label_, "uv", actv_mos_, AlphaSpin);
2388 BlockedTensor::add_mo_space(avirt_label_, "ef", virt_mos_, AlphaSpin);
2389
2390 // define composite spaces
2391 BlockedTensor::add_composite_mo_space("h", "ijkl", {acore_label_, aactv_label_});
2392 BlockedTensor::add_composite_mo_space("p", "abcd", {aactv_label_, avirt_label_});
2393
2394 BlockedTensor G = BlockedTensor::build(CoreTensor, "Gamma1", {"aa"});
2395 BlockedTensor O = BlockedTensor::build(CoreTensor, "O1 PT3 1/3", {"pc", "va"});
2396 BlockedTensor C = BlockedTensor::build(CoreTensor, "C1", {"cp", "av", "pc", "va"});
2397 BlockedTensor T = BlockedTensor::build(CoreTensor, "T2 Amplitudes", {"hhpp"});
2398
2399 Tensor Oac_t = build_and_fill("Oac", {na, nc}, a2);
2400 Tensor Ovc_t = build_and_fill("Ovc", {nv, nc}, b2);
2401 Tensor Ova_t = build_and_fill("Ova", {nv, na}, c2);
2402
2403 O.block("ac")("ia") = Oac_t("ia");
2404 O.block("vc")("ia") = Ovc_t("ia");
2405 O.block("va")("ia") = Ova_t("ia");
2406
2407 O["ia"] = C["bu"] * G["uv"] * T["ivab"];
2408
2409 Tensor Oac = O.block("ac");
2410 double diff_oo = difference(Oac, a2).second;
2411
2412 BlockedTensor::set_expert_mode(false);
2413
2414 return diff_oo;
2415 }
2416
test_batched_incomplete_result_blocks_expert()2417 double test_batched_incomplete_result_blocks_expert() {
2418 BlockedTensor::set_expert_mode(true);
2419 BlockedTensor::reset_mo_spaces();
2420
2421 std::vector<size_t> core_mos{0, 1, 2};
2422 std::vector<size_t> actv_mos{3, 4};
2423 std::vector<size_t> virt_mos{5, 6, 7, 8, 9};
2424
2425 std::vector<size_t> aux_mos{0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
2426
2427 BlockedTensor::add_mo_space("c", "m,n", core_mos, NoSpin);
2428 BlockedTensor::add_mo_space("a", "u,v", actv_mos, NoSpin);
2429 BlockedTensor::add_mo_space("v", "e,f", virt_mos, NoSpin);
2430
2431 BlockedTensor::add_composite_mo_space("h", "i,j,k,l", {"c", "a"});
2432 BlockedTensor::add_composite_mo_space("p", "a,b,c,d", {"a", "v"});
2433 BlockedTensor::add_composite_mo_space("g", "p,q,r,s", {"c", "a", "v"});
2434
2435 BlockedTensor::add_mo_space("L", "g", aux_mos, NoSpin);
2436
2437 std::vector<std::string> od_hhpp_blocks{"hhpv", "hhva", "hcaa", "caaa"};
2438
2439 auto B = BlockedTensor::build(tensor_type, "B", {"Lgg"});
2440 auto T2 = BlockedTensor::build(tensor_type, "T2", {"hhpp"});
2441 auto C2 = BlockedTensor::build(tensor_type, "C2", od_hhpp_blocks);
2442 auto X2 = BlockedTensor::build(tensor_type, "X2", od_hhpp_blocks);
2443
2444 auto init_random = [](BlockedTensor A) {
2445 for (const std::string& block : A.block_labels()) {
2446 std::vector<double>& data = A.block(block).data();
2447 for (size_t i = 0, size = data.size(); i < size; ++i) {
2448 data[i] = double(std::rand()) / double(RAND_MAX);
2449 }
2450 }
2451 };
2452
2453 init_random(B);
2454 init_random(T2);
2455 init_random(C2);
2456 X2["pqrs"] = C2["pqrs"];
2457
2458 C2["ijes"] += batched("e", B["gae"] * B["gbs"] * T2["ijab"]);
2459 C2["ijus"] += batched("u", B["gau"] * B["gbs"] * T2["ijab"]);
2460 C2["ijms"] += batched("m", B["gam"] * B["gbs"] * T2["ijab"]);
2461
2462 X2["ijrs"] += batched("r,", B["gar"] * B["gbs"] * T2["ijab"]);
2463 X2["pqrs"] -= C2["pqrs"];
2464
2465 BlockedTensor::set_expert_mode(false);
2466
2467 return X2.norm(0);
2468 }
2469
main(int argc,char * argv[])2470 int main(int argc, char *argv[])
2471 {
2472 printf(ANSI_COLOR_RESET);
2473 srand(time(nullptr));
2474 ambit::initialize(argc, argv);
2475
2476 printf("==> Simple Operations <==\n\n");
2477
2478 auto test_functions = {
2479 // Expectation, test function, User friendly description
2480 std::make_tuple(kPass, test_mo_space, "Test"),
2481 std::make_tuple(kPass, test_add_mo_space, "Testing composite spaces"),
2482 std::make_tuple(kException, test_add_mo_space_nonexisting_space,
2483 "Testing adding nonexisting space"),
2484 std::make_tuple(kException, test_add_mo_space_repeated_index1,
2485 "Testing adding repeated orbital indices (1)"),
2486 std::make_tuple(kException, test_add_mo_space_repeated_index2,
2487 "Testing adding repeated orbital indices (2)"),
2488 std::make_tuple(kException, test_add_mo_space_repeated_index3,
2489 "Testing adding repeated orbital indices (3)"),
2490 std::make_tuple(kException, test_add_mo_space_repeated_space1,
2491 "Testing adding repeated orbital spaces (1)"),
2492 std::make_tuple(kException, test_add_mo_space_repeated_space2,
2493 "Testing adding repeated orbital spaces (2)"),
2494 std::make_tuple(kException, test_add_mo_space_no_name1,
2495 "Testing adding orbital space with no name (1)"),
2496 std::make_tuple(kException, test_add_mo_space_no_name2,
2497 "Testing adding orbital space with no name (2)"),
2498 std::make_tuple(kException, test_add_mo_space_no_index1,
2499 "Testing adding orbital space with no indices (1)"),
2500 std::make_tuple(kException, test_add_mo_space_no_index2,
2501 "Testing adding orbital space with no indices (2)"),
2502 std::make_tuple(kPass, test_add_mo_space_no_mos,
2503 "Testing adding orbital space with no orbital list"),
2504 std::make_tuple(kPass, test_block_creation1,
2505 "Testing blocked tensor creation (1)"),
2506 std::make_tuple(kPass, test_block_creation2,
2507 "Testing blocked tensor creation (2)"),
2508 std::make_tuple(kPass, test_block_creation3,
2509 "Testing blocked tensor creation (3)"),
2510 std::make_tuple(kException, test_block_creation_bad_rank,
2511 "Testing blocked tensor creation with variable rank"),
2512 std::make_tuple(kPass, test_block_norm_1,
2513 "Testing blocked tensor 1-norm"),
2514 std::make_tuple(kPass, test_block_norm_2,
2515 "Testing blocked tensor 2-norm"),
2516 std::make_tuple(kPass, test_block_norm_3,
2517 "Testing blocked tensor inf-norm"),
2518 std::make_tuple(kPass, test_block_zero, "Testing blocked tensor zero"),
2519 std::make_tuple(kPass, test_block_scale,
2520 "Testing blocked tensor scale"),
2521 std::make_tuple(kPass, test_block_labels1,
2522 "Testing blocked tensor labeling (1)"),
2523 std::make_tuple(kPass, test_block_retrive_block1,
2524 "Testing blocked tensor retrieve existing block"),
2525 std::make_tuple(kException, test_block_retrive_block2,
2526 "Testing blocked tensor retrieve ambiguous block"),
2527 std::make_tuple(kException, test_block_retrive_block3,
2528 "Testing blocked tensor retrieve null block (1)"),
2529 std::make_tuple(kException, test_block_retrive_block4,
2530 "Testing blocked tensor retrieve null block (2)"),
2531 std::make_tuple(kPass, test_block_iterator_1,
2532 "Testing blocked tensor iterator (1)"),
2533 // std::make_tuple(kException, test_copy,
2534 // "Testing blocked tensor copy"),
2535 std::make_tuple(kPass, test_Cij_equal_Aji,
2536 "Testing blocked tensor C(\"ij\") = A(\"ji\")"),
2537 std::make_tuple(kException, test_Aij_equal_Aji,
2538 "Testing blocked tensor A(\"ij\") = A(\"ji\")"),
2539 std::make_tuple(kPass, test_Cijab_plus_equal_Aaibj,
2540 "Testing blocked tensor C(\"ijab\") += A(\"aibj\")"),
2541 std::make_tuple(kPass, test_Cbija_minus_equal_Ajabi,
2542 "Testing blocked tensor C(\"bija\") -= A(\"jabi\")"),
2543 std::make_tuple(kPass, test_Cij_times_equal_double,
2544 "Testing blocked tensor A(\"ij\") *= double"),
2545 std::make_tuple(kPass, test_Cip_times_equal_double,
2546 "Testing blocked tensor A(\"ip\") *= double"),
2547 std::make_tuple(
2548 kPass, test_Cij_equal_Aik_B_jk,
2549 "Testing blocked tensor C(\"ij\") = A(\"ik\") * B(\"jk\")"),
2550 std::make_tuple(
2551 kPass, test_Cij_equal_Aip_B_jp,
2552 "Testing blocked tensor C(\"ij\") = A(\"ip\") * B(\"jp\") (1)"),
2553 std::make_tuple(
2554 kException, test_Cij_equal_Aip_B_jp_fail,
2555 "Testing blocked tensor C(\"ij\") = A(\"ip\") * B(\"jp\") (2)"),
2556 std::make_tuple(
2557 kPass, test_Cij_equal_Aip_B_jp_expert,
2558 "Testing blocked tensor C(\"ij\") = A(\"ip\") * B(\"jp\") (3)"),
2559 std::make_tuple(
2560 kPass, test_Cpq_equal_Apq_B_pq_expert,
2561 "Testing blocked tensor C(\"pq\") = A(\"pq\") * B(\"pq\")"),
2562 std::make_tuple(
2563 kPass, test_Cij_equal_half_Aia_B_aj,
2564 "Testing blocked tensor C(\"ij\") = 0.5 * A(\"ia\") * B(\"aj\")"),
2565 std::make_tuple(kPass, test_Cij_plus_equal_half_Aai_B_ja,
2566 "Testing blocked tensor C(\"ij\") += A(\"ai\") * (0.5 "
2567 "* B(\"ja\"))"),
2568 std::make_tuple(
2569 kPass, test_Cij_minus_equal_Aik_B_jk,
2570 "Testing blocked tensor C(\"ij\") -= A(\"ik\") * B(\"jk\")"),
2571 std::make_tuple(
2572 kPass, test_greek_Cij_equal_Aik_B_jk,
2573 "Testing blocked tensor C(\"ij\") = A(\"ik\") * B(\"jk\") [Greek]"),
2574 std::make_tuple(
2575 kException, test_Aij_equal_Aik_B_jk,
2576 "Testing blocked tensor A(\"ij\") = A(\"ik\") * B(\"jk\")"),
2577 std::make_tuple(kPass, test_chain_multiply,
2578 "Testing blocked tensor chain multiply (1)"),
2579 std::make_tuple(kPass, test_chain_multiply2,
2580 "Testing blocked tensor chain multiply (2)"),
2581 std::make_tuple(
2582 kPass, test_Cij_equal_Aij_plus_Bij,
2583 "Testing blocked tensor C(\"ij\") = A(\"ij\") + B(\"ij\")"),
2584 std::make_tuple(
2585 kPass, test_Dij_equal_2_times_Aij_plus_Bij,
2586 "Testing blocked tensor C(\"ij\") = 2 * (A(\"ij\") - B(\"ij\"))"),
2587 std::make_tuple(kPass, test_Cij_minus_equal_3_times_Aij_minus_2_Bji,
2588 "Testing blocked tensor C(\"ij\") -= 3 * (A(\"ij\") - "
2589 "2 B(\"ji\"))"),
2590 std::make_tuple(
2591 kPass, test_Cij_equal_negate_Aij_plus_Bij,
2592 "Testing blocked tensor C(\"ij\") = - (A(\"ij\") + B(\"ij\"))"),
2593 std::make_tuple(
2594 kPass, test_Cia_plus_equal_Aia_minus_three_Bai,
2595 "Testing blocked tensor C(\"ia\") += A(\"ia\") - 3 * B(\"ai\")"),
2596 std::make_tuple(kPass, test_Dij_equal_Aij_times_Bij_plus_Cij,
2597 "Testing blocked tensor distributive (1)"),
2598 std::make_tuple(kPass, test_Dij_plus_equal_Bij_plus_Cij_times_Aij,
2599 "Testing blocked tensor distributive (2)"),
2600 std::make_tuple(kPass, test_Dij_equal_Bij_plus_Cij_times_Aij,
2601 "Testing blocked tensor distributive (3)"),
2602 std::make_tuple(kPass, test_F_equal_D_times_2g_minus_g,
2603 "Testing blocked tensor distributive (4)"),
2604 std::make_tuple(
2605 kPass, test_dot_product,
2606 "Testing blocked tensor dot product C = A(\"ia\") * B(\"ai\")"),
2607 std::make_tuple(
2608 kException, test_dot_product_fail1,
2609 "Testing blocked tensor dot product index mismatch (1)"),
2610 std::make_tuple(
2611 kException, test_dot_product_fail2,
2612 "Testing blocked tensor dot product index mismatch (2)"),
2613 std::make_tuple(kException, test_contraction_exception1,
2614 "Testing blocked tensor contraction exception (1)"),
2615 std::make_tuple(
2616 kPass, test_Cij_equal_Aik_B_jk_batched,
2617 "C(\"ij\") = batched(\"i\", A(\"ik\") * B(\"jk\"))"),
2618 std::make_tuple(
2619 kPass, test_Cij_equal_Aip_B_jp_batched,
2620 "C(\"ij\") = batched(\"j\", A(\"ip\") * B(\"jp\")) (1)"),
2621 std::make_tuple(
2622 kException, test_Cij_equal_Aip_B_jp_fail_batched,
2623 "C(\"ij\") = batched(\"i\", A(\"ip\") * B(\"jp\")) (2)"),
2624 std::make_tuple(
2625 kPass, test_Cij_equal_Aip_B_jp_expert_batched,
2626 "C(\"ij\") = batched(\"i\", A(\"ip\") * B(\"jp\")) (3)"),
2627 std::make_tuple(
2628 kPass, test_Cpq_equal_Apq_B_pq_expert_batched,
2629 "C(\"pq\") = batched(\"p\", A(\"pq\") * B(\"pq\"))"),
2630 std::make_tuple(
2631 kPass, test_Cij_equal_half_Aia_B_aj_batched,
2632 "C(\"ij\") = batched(\"j\", 0.5 * A(\"ia\") * B(\"aj\"))"),
2633 std::make_tuple(kPass, test_Cij_plus_equal_half_Aai_B_ja_batched,
2634 "C(\"ij\") += batched(\"i\", A(\"ai\") * (0.5 "
2635 "* B(\"ja\")))"),
2636 std::make_tuple(
2637 kPass, test_Cij_minus_equal_Aik_B_jk_batched,
2638 "C(\"ij\") -= batched(\"i\", A(\"ik\") * B(\"jk\"))"),
2639 std::make_tuple(
2640 kPass, test_greek_Cij_equal_Aik_B_jk_batched,
2641 "C(\"ij\") = batched(\"j\", A(\"ik\") * B(\"jk\")) [Greek]"),
2642 std::make_tuple(
2643 kException, test_Aij_equal_Aik_B_jk_batched,
2644 "A(\"ij\") = batched(\"j\", A(\"ik\") * B(\"jk\"))"),
2645 std::make_tuple(
2646 kPass, test_batched_with_factor,
2647 "C2[\"ijrs\"] = batched(\"r\", 0.5 * A[\"abrs\"] * B[\"ijab\"])"),
2648 std::make_tuple(
2649 kPass, test_batch_all_indices,
2650 "C2[\"ijrs\"] = batched(\"srji\", 0.5 * A[\"abrs\"] * B[\"ijab\"])"),
2651 std::make_tuple(
2652 kPass, test_Oia_equal_Cbu_Guv_Tivab_expert,
2653 "O[\"ia\"] = C[\"bu\"] * G[\"uv\"] * T[\"ivab\"]"),
2654 std::make_tuple(kPass, test_batched_incomplete_result_blocks_expert,
2655 "C[\"ijrs\"] += batched(\"r,\", B[\"gar\"] * B[\"gbs\"] * T[\"ijab\"])")
2656 };
2657
2658 std::vector<std::tuple<std::string, TestResult, double>> results;
2659
2660 printf(ANSI_COLOR_RESET);
2661
2662 printf("\n %-60s %12s %s", "Description", "Max. error", "Result");
2663 printf("\n %s", std::string(83, '-').c_str());
2664
2665 bool success = true;
2666 for (auto test_function : test_functions)
2667 {
2668 printf("\n %-60s", std::get<2>(test_function));
2669 double result = 0.0;
2670 TestResult tresult = kPass, report_result = kPass;
2671 std::string exception;
2672 try
2673 {
2674 result = std::get<1>(test_function)();
2675
2676 // Did the test pass based on returned value?
2677 tresult = std::fabs(result) < epsilon ? kPass : kFail;
2678 // Was the tresult the expected result? If so color green else red.
2679 report_result =
2680 tresult == std::get<0>(test_function) ? kPass : kFail;
2681 }
2682 catch (std::exception &e)
2683 {
2684 // was an exception expected?
2685 tresult = kException;
2686 report_result =
2687 tresult == std::get<0>(test_function) ? kPass : kException;
2688
2689 // printf("\n %s",e.what());
2690 if (report_result == kException)
2691 {
2692 exception = e.what();
2693 }
2694 }
2695 printf(" %7e", result);
2696 switch (report_result)
2697 {
2698 case kPass:
2699 printf(ANSI_COLOR_GREEN);
2700 break;
2701 case kFail:
2702 printf(ANSI_COLOR_RED);
2703 break;
2704 default:
2705 printf(ANSI_COLOR_YELLOW);
2706 }
2707 switch (tresult)
2708 {
2709 case kPass:
2710 printf(" Passed" ANSI_COLOR_RESET);
2711 break;
2712 case kFail:
2713 printf(" Failed" ANSI_COLOR_RESET);
2714 break;
2715 default:
2716 printf(" Exception" ANSI_COLOR_RESET);
2717 }
2718
2719 if (report_result == kException)
2720 printf("\n Unexpected: %s", exception.c_str());
2721 if (report_result != kPass)
2722 success = false;
2723 }
2724 printf("\n %s", std::string(83, '-').c_str());
2725 printf("\n Tests: %s\n", success ? "All passed" : "Some failed");
2726
2727 // printf("%s\n",std::string(82,'-').c_str());
2728 // printf("%-50s %-9s %-9s %11s\n", "Description", "Expected",
2729 // "Observed", "Delta");
2730 // printf("%s\n",std::string(82,'-').c_str());
2731 // success = true;
2732 // success &= test_function(try_relative_difference, "Relative
2733 // Difference", kExact);
2734 // success &= test_function(try_1_norm , "1-Norm"
2735 // , kEpsilon);
2736 // success &= test_function(try_2_norm , "2-Norm"
2737 // , kEpsilon);
2738 // success &= test_function(try_inf_norm , "Inf-Norm"
2739 // , kEpsilon);
2740 // success &= test_function(try_zero , "Zero" , kExact);
2741 // success &= test_function(try_copy , "Copy" , kExact);
2742 // success &= test_function(try_scale , "Scale" , kExact);
2743 // printf("%s\n",std::string(82,'-').c_str());
2744 // printf("Tests: %s\n\n",success ? "All passed" : "Some failed");
2745
2746 // printf("==> Slice Operations <==\n\n");
2747 // printf("%s\n",std::string(82,'-').c_str());
2748 // printf("%-50s %-9s %-9s %11s\n", "Description", "Expected",
2749 // "Observed", "Delta");
2750 // printf("%s\n",std::string(82,'-').c_str());
2751 // success = true;
2752 // success &= test_function(try_slice_rank0 , "Full Slice Rank-0"
2753 // , kExact);
2754 // success &= test_function(try_slice_rank3 , "Full Slice Rank-3"
2755 // , kExact);
2756 // printf("%s\n",std::string(82,'-').c_str());
2757 // printf("Tests: %s\n\n",success ? "All passed" : "Some failed");
2758
2759 // printf("==> Permute Operations <==\n\n");
2760 // success = true;
2761 // printf("%s\n",std::string(82,'-').c_str());
2762 // printf("%-50s %-9s %-9s %11s\n", "Description", "Expected",
2763 // "Observed", "Delta");
2764 // mode = 0; alpha = 1.0; beta = 0.0;
2765 // printf("%s\n",std::string(82,'-').c_str());
2766 // printf("Explicit: alpha = %11.3E, beta = %11.3E\n", alpha, beta);
2767 // printf("%s\n",std::string(82,'-').c_str());
2768 // success &= test_function(try_permute_rank0 , "Permute Rank-0"
2769 // , kExact);
2770 // success &= test_function(try_permute_rank1_i , "Permute Rank-1 i"
2771 // , kExact);
2772 // success &= test_function(try_permute_rank2_ij , "Permute Rank-2 ij"
2773 // , kExact);
2774 // success &= test_function(try_permute_rank2_ji , "Permute Rank-2 ji"
2775 // , kExact);
2776 // success &= test_function(try_permute_rank3_ijk , "Permute Rank-3 ijk"
2777 // , kExact);
2778 // success &= test_function(try_permute_rank3_ikj , "Permute Rank-3 ikj"
2779 // , kExact);
2780 // success &= test_function(try_permute_rank3_jik , "Permute Rank-3 jik"
2781 // , kExact);
2782 // success &= test_function(try_permute_rank3_jki , "Permute Rank-3 jki"
2783 // , kExact);
2784 // success &= test_function(try_permute_rank3_kij , "Permute Rank-3 kij"
2785 // , kExact);
2786 // success &= test_function(try_permute_rank3_kji , "Permute Rank-3 kji"
2787 // , kExact);
2788 // success &= test_function(try_permute_rank4_ijkl , "Permute Rank-4
2789 // ijkl" , kExact);
2790 // success &= test_function(try_permute_rank4_lkji , "Permute Rank-4
2791 // lkji" , kExact);
2792 // success &= test_function(try_permute_rank4_ijlk , "Permute Rank-4
2793 // ijlk" , kExact);
2794 // success &= test_function(try_permute_rank4_jikl , "Permute Rank-4
2795 // jikl" , kExact);
2796 // success &= test_function(try_permute_rank4_ikjl , "Permute Rank-4
2797 // ikjl" , kExact);
2798 // success &= test_function(try_permute_rank4_lkji , "Permute Rank-4
2799 // lkji" , kExact);
2800 // mode = 0; alpha = random_double(); beta = random_double();
2801 // printf("%s\n",std::string(82,'-').c_str());
2802 // printf("Explicit: alpha = %11.3E, beta = %11.3E\n", alpha, beta);
2803 // printf("%s\n",std::string(82,'-').c_str());
2804 // success &= test_function(try_permute_rank0 , "Permute Rank-0"
2805 // , kExact);
2806 // success &= test_function(try_permute_rank1_i , "Permute Rank-1 i"
2807 // , kExact);
2808 // success &= test_function(try_permute_rank2_ij , "Permute Rank-2 ij"
2809 // , kExact);
2810 // success &= test_function(try_permute_rank2_ji , "Permute Rank-2 ji"
2811 // , kExact);
2812 // success &= test_function(try_permute_rank3_ijk , "Permute Rank-3 ijk"
2813 // , kExact);
2814 // success &= test_function(try_permute_rank3_ikj , "Permute Rank-3 ikj"
2815 // , kExact);
2816 // success &= test_function(try_permute_rank3_jik , "Permute Rank-3 jik"
2817 // , kExact);
2818 // success &= test_function(try_permute_rank3_jki , "Permute Rank-3 jki"
2819 // , kExact);
2820 // success &= test_function(try_permute_rank3_kij , "Permute Rank-3 kij"
2821 // , kExact);
2822 // success &= test_function(try_permute_rank3_kji , "Permute Rank-3 kji"
2823 // , kExact);
2824 // success &= test_function(try_permute_rank4_ijkl , "Permute Rank-4
2825 // ijkl" , kExact);
2826 // success &= test_function(try_permute_rank4_ijlk , "Permute Rank-4
2827 // ijlk" , kExact);
2828 // success &= test_function(try_permute_rank4_jikl , "Permute Rank-4
2829 // jikl" , kExact);
2830 // success &= test_function(try_permute_rank4_ikjl , "Permute Rank-4
2831 // ikjl" , kExact);
2832 // success &= test_function(try_permute_rank4_lkji , "Permute Rank-4
2833 // lkji" , kExact);
2834 // mode = 1; alpha = 1.0; beta = 0.0;
2835 // printf("%s\n",std::string(82,'-').c_str());
2836 // printf("Operator Overloading: =\n");
2837 // printf("%s\n",std::string(82,'-').c_str());
2838 // success &= test_function(try_permute_rank0 , "Permute Rank-0"
2839 // , kExact);
2840 // success &= test_function(try_permute_rank1_i , "Permute Rank-1 i"
2841 // , kExact);
2842 // success &= test_function(try_permute_rank2_ij , "Permute Rank-2 ij"
2843 // , kExact);
2844 // success &= test_function(try_permute_rank2_ji , "Permute Rank-2 ji"
2845 // , kExact);
2846 // success &= test_function(try_permute_rank3_ijk , "Permute Rank-3 ijk"
2847 // , kExact);
2848 // success &= test_function(try_permute_rank3_ikj , "Permute Rank-3 ikj"
2849 // , kExact);
2850 // success &= test_function(try_permute_rank3_jik , "Permute Rank-3 jik"
2851 // , kExact);
2852 // success &= test_function(try_permute_rank3_jki , "Permute Rank-3 jki"
2853 // , kExact);
2854 // success &= test_function(try_permute_rank3_kij , "Permute Rank-3 kij"
2855 // , kExact);
2856 // success &= test_function(try_permute_rank3_kji , "Permute Rank-3 kji"
2857 // , kExact);
2858 // success &= test_function(try_permute_rank4_ijkl , "Permute Rank-4
2859 // ijkl" , kExact);
2860 // success &= test_function(try_permute_rank4_ijlk , "Permute Rank-4
2861 // ijlk" , kExact);
2862 // success &= test_function(try_permute_rank4_jikl , "Permute Rank-4
2863 // jikl" , kExact);
2864 // success &= test_function(try_permute_rank4_ikjl , "Permute Rank-4
2865 // ikjl" , kExact);
2866 // success &= test_function(try_permute_rank4_lkji , "Permute Rank-4
2867 // lkji" , kExact);
2868 // mode = 2; alpha = 1.0; beta = 1.0;
2869 // printf("%s\n",std::string(82,'-').c_str());
2870 // printf("Operator Overloading: +=\n");
2871 // printf("%s\n",std::string(82,'-').c_str());
2872 // success &= test_function(try_permute_rank0 , "Permute Rank-0"
2873 // , kExact);
2874 // success &= test_function(try_permute_rank1_i , "Permute Rank-1 i"
2875 // , kExact);
2876 // success &= test_function(try_permute_rank2_ij , "Permute Rank-2 ij"
2877 // , kExact);
2878 // success &= test_function(try_permute_rank2_ji , "Permute Rank-2 ji"
2879 // , kExact);
2880 // success &= test_function(try_permute_rank3_ijk , "Permute Rank-3 ijk"
2881 // , kExact);
2882 // success &= test_function(try_permute_rank3_ikj , "Permute Rank-3 ikj"
2883 // , kExact);
2884 // success &= test_function(try_permute_rank3_jik , "Permute Rank-3 jik"
2885 // , kExact);
2886 // success &= test_function(try_permute_rank3_jki , "Permute Rank-3 jki"
2887 // , kExact);
2888 // success &= test_function(try_permute_rank3_kij , "Permute Rank-3 kij"
2889 // , kExact);
2890 // success &= test_function(try_permute_rank3_kji , "Permute Rank-3 kji"
2891 // , kExact);
2892 // success &= test_function(try_permute_rank4_ijkl , "Permute Rank-4
2893 // ijkl" , kExact);
2894 // success &= test_function(try_permute_rank4_ijlk , "Permute Rank-4
2895 // ijlk" , kExact);
2896 // success &= test_function(try_permute_rank4_jikl , "Permute Rank-4
2897 // jikl" , kExact);
2898 // success &= test_function(try_permute_rank4_ikjl , "Permute Rank-4
2899 // ikjl" , kExact);
2900 // success &= test_function(try_permute_rank4_lkji , "Permute Rank-4
2901 // lkji" , kExact);
2902 // mode = 3; alpha = -1.0; beta = 1.0;
2903 // printf("%s\n",std::string(82,'-').c_str());
2904 // printf("Operator Overloading: -=\n");
2905 // printf("%s\n",std::string(82,'-').c_str());
2906 // success &= test_function(try_permute_rank0 , "Permute Rank-0"
2907 // , kExact);
2908 // success &= test_function(try_permute_rank1_i , "Permute Rank-1 i"
2909 // , kExact);
2910 // success &= test_function(try_permute_rank2_ij , "Permute Rank-2 ij"
2911 // , kExact);
2912 // success &= test_function(try_permute_rank2_ji , "Permute Rank-2 ji"
2913 // , kExact);
2914 // success &= test_function(try_permute_rank3_ijk , "Permute Rank-3 ijk"
2915 // , kExact);
2916 // success &= test_function(try_permute_rank3_ikj , "Permute Rank-3 ikj"
2917 // , kExact);
2918 // success &= test_function(try_permute_rank3_jik , "Permute Rank-3 jik"
2919 // , kExact);
2920 // success &= test_function(try_permute_rank3_jki , "Permute Rank-3 jki"
2921 // , kExact);
2922 // success &= test_function(try_permute_rank3_kij , "Permute Rank-3 kij"
2923 // , kExact);
2924 // success &= test_function(try_permute_rank3_kji , "Permute Rank-3 kji"
2925 // , kExact);
2926 // success &= test_function(try_permute_rank4_ijkl , "Permute Rank-4
2927 // ijkl" , kExact);
2928 // success &= test_function(try_permute_rank4_ijlk , "Permute Rank-4
2929 // ijlk" , kExact);
2930 // success &= test_function(try_permute_rank4_jikl , "Permute Rank-4
2931 // jikl" , kExact);
2932 // success &= test_function(try_permute_rank4_ikjl , "Permute Rank-4
2933 // ikjl" , kExact);
2934 // success &= test_function(try_permute_rank4_lkji , "Permute Rank-4
2935 // lkji" , kExact);
2936 // printf("%s\n",std::string(82,'-').c_str());
2937 // printf("Tests: %s\n\n",success ? "All Passed" : "Some Failed");
2938
2939 // printf("==> Permute Exceptions <==\n\n");
2940 // success = true;
2941 // printf("%s\n",std::string(82,'-').c_str());
2942 // printf("%-50s %-9s %-9s %11s\n", "Description", "Expected",
2943 // "Observed", "Delta");
2944 // mode = 0; alpha = 1.0; beta = 0.0;
2945 // printf("%s\n",std::string(82,'-').c_str());
2946 // printf("Explicit: alpha = %11.3E, beta = %11.3E\n", alpha, beta);
2947 // printf("%s\n",std::string(82,'-').c_str());
2948 // success &= test_function(try_permute_label_fail , "Permute Label
2949 // Fail" , kException);
2950 // success &= test_function(try_permute_rank_fail , "Permute Rank Fail"
2951 // , kException);
2952 // success &= test_function(try_permute_index_fail , "Permute Index
2953 // Fail" , kException);
2954 // success &= test_function(try_permute_size_fail , "Permute Size Fail"
2955 // , kException);
2956 // printf("%s\n",std::string(82,'-').c_str());
2957 // printf("Tests: %s\n\n",success ? "All Passed" : "Some Failed");
2958
2959 ambit::finalize();
2960
2961 return success ? EXIT_SUCCESS : EXIT_FAILURE;
2962 }
2963