182 {
183
185
186 auto core_log = logging::core::get();
187 core_log->add_sink(
190 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
192
193 try {
194
195 auto print_ddg = [](
auto &
t,
auto str =
"") {
196 constexpr double eps = 1e-6;
197 for (int ii = 0; ii != 3; ++ii)
198 for (int jj = 0; jj != 3; ++jj)
199 for (int kk = 0; kk != 3; ++kk)
200 for (int ll = 0; ll != 3; ++ll) {
201 double v =
t(ii, jj, kk, ll);
202 double w = std::abs(
v) <
eps ? 0 :
v;
204 << str << std::fixed << std::setprecision(3) << std::showpos
205 << ii + 1 << " " << jj + 1 << " " << kk + 1 << " " << ll + 1
207 }
208 };
209
210 auto print_ddg_direction = [](
auto &
t,
auto kk,
int ll) {
211 for (int ii = 0; ii != 3; ++ii)
212 for (int jj = 0; jj <= ii; ++jj)
214 << ii + 1 << " " << jj + 1 << " " << kk + 1 << " " << ll + 1
215 <<
" : " <<
t(ii, jj, kk, ll);
216 };
217
219 for (int ii = 0; ii != 3; ++ii)
220 for (int jj = 0; jj != 3; ++jj)
222 << ii + 1 <<
" " << jj + 1 <<
" : " <<
t(ii, jj);
223 };
224
225 enum swap { swap12, swap01 };
226 auto run_lapack = [](
auto &
a, swap s = swap12) {
227 int info;
228 double wkopt;
230
232
234
236
238
239
240 int lwork = -1;
241 info =
lapack_dsyev(
'V',
'U', 3,
a.data(), 3, w, &wkopt, lwork);
242 if (info > 0)
243 THROW_MESSAGE(
"The algorithm failed to compute eigenvalues.");
245 std::vector<double> work(lwork);
246
247 info =
lapack_dsyev(
'V',
'U', 3,
a.data(), 3, w, &*work.begin(), lwork);
248 if (info > 0)
249 THROW_MESSAGE(
"The algorithm failed to compute eigenvalues.");
250
251 if (s == swap12) {
253
254 a[0 * 3 + 0],
a[0 * 3 + 1],
a[0 * 3 + 2],
255
256 a[2 * 3 + 0],
a[2 * 3 + 1],
a[2 * 3 + 2],
257
258 a[1 * 3 + 0],
a[1 * 3 + 1],
a[1 * 3 + 2]};
259
261 return std::make_tuple(t_a, t_eig_vec, t_eig_vals);
262 }
263
265
266 a[1 * 3 + 0],
a[1 * 3 + 1],
a[1 * 3 + 2],
267
268 a[0 * 3 + 0],
a[0 * 3 + 1],
a[0 * 3 + 2],
269
270 a[2 * 3 + 0],
a[2 * 3 + 1],
a[2 * 3 + 2],
271
272 };
273
275 return std::make_tuple(t_a, t_eig_vec, t_eig_vals);
276 };
277
278
279 {
281
282 1., 0.1, -0.5,
283
284 0.1, 2., 0.,
285
286 -0.5, 0., 3.};
287
289
290 0.969837, -0.0860972, 0.228042,
291
292 0.0790574, 0.996073, 0.0398449,
293
294 -0.230577, -0.0206147, 0.972836};
295
297
299
300 1., 0., 0.,
301
302 0., 1., 0.,
303
304 0., 0., 1.};
305
307 t_S_sym(
i,
j) = (t_S(
i,
j) || t_S(
j,
i)) / 2.;
308
309 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"Diff A";
311
312
313 {
314 auto f = [](
double v) {
return exp(
v); };
315 auto d_f = [](
double v) {
return exp(
v); };
316 auto dd_f = [](
double v) {
return exp(
v); };
317
319 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"Reconstruct mat";
321
323
324 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"Diff";
325 print_ddg_direction(t_d, 0, 2);
326
327 auto t_dd =
329
330 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"Diff Diff";
331 print_ddg_direction(t_dd, 0, 2);
332
333 auto norm2_t_b = t_b(
i,
j) * t_b(
i,
j);
334 MOFEM_LOG(
"ATOM_TEST", Sev::inform) <<
"norm2_t_b " << norm2_t_b;
335
337 MOFEM_LOG(
"ATOM_TEST", Sev::inform) <<
"norm2_t_d " << norm2_t_d;
338
340 MOFEM_LOG(
"ATOM_TEST", Sev::inform) <<
"norm2_t_dd " << norm2_t_dd;
341
342 constexpr double regression_t_b = 572.543;
343 constexpr double regression_t_d = 859.939;
344 constexpr double regression_t_dd = 859.939;
345
346 constexpr double eps = 1e-2;
347 if (std::abs(norm2_t_b - regression_t_b) >
eps)
349 if (std::abs(norm2_t_d - regression_t_d) >
eps)
351 if (std::abs(norm2_t_dd - regression_t_dd) >
eps)
353 }
354
355
356
357
358 {
359
360 std::array<double, 9>
a{1., 0.1, -0.5,
361
362 0.1, 2., 0.,
363
364 -0.5, 0., 3.};
365
366 auto tuple = run_lapack(a);
367
368 auto &t_eig_vec = std::get<1>(tuple);
369 auto &t_eig_vals = std::get<2>(tuple);
370
371 auto t_eig_val_diff =
372 (t_eig_vals(
i) - t_L(
i)) * (t_eig_vals(
i) - t_L(
i));
374 << "t_eig_val_diff " << t_eig_val_diff;
375
376 auto f = [](
double v) {
return exp(
v); };
377
380 t_c(
i,
j) -= t_b(
i,
j);
382
383 auto norm2_t_c = t_c(
i,
j) * t_c(
i,
j);
385 << "Reconstruct mat difference with lapack eigen valyes and "
386 "vectors "
387 << norm2_t_c;
388
389 constexpr double eps = 1e-8;
390 if (fabs(norm2_t_c) >
eps)
392 "Matrix not reeconstructed");
393 }
394
398
399
400 {
401 auto f = [](
double v) {
return v; };
402 auto d_f = [](
double v) {
return 1; };
403 auto dd_f = [](
double v) {
return 0; };
404
405 constexpr double eps = 1e-10;
406 {
408 t_b(
i,
j) -= (t_A(
i,
j) || t_A(
j,
i)) / 2;
409 auto norm2_t_b = t_b(
i,
j) * t_b(
i,
j);
411 << "Result should be matrix itself " << norm2_t_b;
414 "This norm should be zero");
415 }
416
417 {
418
421
422 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
423 print_ddg(t_d_a, "hand ");
424 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
425 print_ddg(t_d, "code ");
426
429 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
430 if (nrm2_t_d_a >
eps)
432 "This norm should be zero");
433 }
434
435 {
437
438 1., 0., 0.,
439
440 0., 1., 0.,
441
442 0., 0., 1.};
443
445 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
446
447 auto t_dd =
449
451 MOFEM_LOG(
"ATOM_TEST", Sev::inform) <<
"norm2_t_dd " << norm2_t_dd;
452 if (norm2_t_dd >
eps)
454 "This norm should be zero");
455 }
456 }
457
458
459 {
460 auto f = [](
double v) {
return v *
v; };
461 auto d_f = [](
double v) {
return 2 *
v; };
462 auto dd_f = [](
double v) {
return 2; };
463
464 constexpr double eps = 1e-9;
465
466
467 {
470 t_a(
i,
j) = t_b(
i,
j) - t_A(
i,
k) * t_A(
k,
j);
472 auto norm2_t_a = t_a(
i,
j) * t_a(
i,
j);
474 << "Result should be matrix times matrix " << norm2_t_a;
477 "This norm should be zero");
478 }
479
480
481 {
483 print_ddg_direction(t_d, 0, 2);
485 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
486 print_ddg(t_d_a, "hand ");
487 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
488 print_ddg(t_d, "code ");
491 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
492 if (nrm2_t_d_a >
eps)
494 "This norm should be zero");
495 }
496
497
498 {
500
501 1., 1. / 2., 1. / 3.,
502
503 2. / 2., 1., 2. / 3.,
504
505 3. / 2., 1., 3. / 3.};
506
508 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
509
510 auto t_dd =
513
514 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd_a";
515 print_ddg(t_dd_a, "hand ");
516 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd";
517 print_ddg(t_dd, "code ");
518
521 << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
522 if (nrm2_t_dd_a >
eps)
524 "This norm should be zero");
525 }
526 }
527 }
528
529
530 {
531
532 std::array<double, 9>
a{5., 4., 0,
533
534 4., 5, 0.,
535
536 0.0, 0., 9};
537
538 auto tuple = run_lapack(a, swap01);
539 auto &t_a = std::get<0>(tuple);
540 auto &t_eig_vecs = std::get<1>(tuple);
541 auto &t_eig_vals = std::get<2>(tuple);
542
543 auto f = [](
double v) {
return v; };
544 auto d_f = [](
double v) {
return 1; };
545
546 constexpr double eps = 1e-10;
547 {
549 t_b(
i,
j) -= (t_a(
i,
j) || t_a(
j,
i)) / 2;
550 auto norm2_t_b = t_b(
i,
j) * t_b(
i,
j);
552 << "Result should be matrix itself " << norm2_t_b;
555 "This norm should be zero");
556 }
557
558 {
561 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
562 print_ddg(t_d_a, "hand ");
563 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
564 print_ddg(t_d, "code ");
567 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
568 if (nrm2_t_d_a >
eps)
570 "This norm should be zero");
571 }
572
573 {
574 auto f = [](
double v) {
return v *
v; };
575 auto d_f = [](
double v) {
return 2 *
v; };
578 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
579 print_ddg(t_d_a, "hand ");
580 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
581 print_ddg(t_d, "code ");
584 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
585 if (nrm2_t_d_a >
eps)
587 "This norm should be zero");
588 }
589 }
590
591
592 {
593
594 std::array<double, 9>
a{4., 0., 0,
595
596 0., 4., 0.,
597
598 0.0, 0., 4.};
599
600 auto f = [](
double v) {
return v; };
601 auto d_f = [](
double v) {
return 1; };
602
603 auto tuple = run_lapack(a);
604 auto &t_a = std::get<0>(tuple);
605 auto &t_eig_vecs = std::get<1>(tuple);
606 auto &t_eig_vals = std::get<2>(tuple);
607
608 constexpr double eps = 1e-10;
609 {
611 t_b(
i,
j) -= (t_a(
i,
j) || t_a(
j,
i)) / 2;
612 auto norm2_t_b = t_b(
i,
j) * t_b(
i,
j);
614 << "Result should be matrix itself " << norm2_t_b;
617 "This norm should be zero");
618 }
619
620 {
623 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
624 print_ddg(t_d_a, "hand ");
625 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
626 print_ddg(t_d, "code ");
629 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
630 if (nrm2_t_d_a >
eps)
632 "This norm should be zero");
633 }
634
635 {
636 auto f = [](
double v) {
return v *
v; };
637 auto d_f = [](
double v) {
return 2 *
v; };
640 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
641 print_ddg(t_d_a, "hand ");
642 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
643 print_ddg(t_d, "code ");
646 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
647 if (nrm2_t_d_a >
eps)
649 "This norm should be zero");
650 }
651 }
652
653
654 {
655
656 std::array<double, 9>
a{0.1, 0., 0.,
657
658 0., 0.1, 0.,
659
660 0., 0., 0.1};
661
662 auto tuple = run_lapack(a);
663
664 auto &t_eig_vecs = std::get<1>(tuple);
665 auto &t_eig_vals = std::get<2>(tuple);
666
667 t_eig_vals(0) -= 1e-4;
668 t_eig_vals(2) += 1e-4;
669
670 constexpr double eps = 1e-10;
671
672 auto f = [](
double v) {
return v; };
673 auto d_f = [](
double v) {
return 1; };
674 auto dd_f = [](
double v) {
return 0; };
675
677
678 1., 1. / 2., 1. / 3.,
679
680 2. / 2., 1., 2. / 3.,
681
682 3. / 2., 1., 3. / 3.};
683
685 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
686
688 dd_f, t_S_sym, 1);
689
690 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd";
691 print_ddg(t_dd, "test ");
692
695 << "Direvarive hand calculation minus code " << nrm2_t_dd;
698 "This norm should be zero");
699 }
700
701
702 {
703
704 std::array<double, 9>
a{2, 0., 0.,
705
706 0., 2, 0.,
707
708 0., 0., 2};
709
710 auto tuple = run_lapack(a);
711
712 auto &t_eig_vecs = std::get<1>(tuple);
713 auto &t_eig_vals = std::get<2>(tuple);
714
715 constexpr double eps = 1e-10;
716
717 auto f = [](
double v) {
return v *
v; };
718 auto d_f = [](
double v) {
return 2 *
v; };
719 auto dd_f = [](
double v) {
return 2; };
721
722 1., 1. / 2., 1. / 3.,
723
724 2. / 1., 1., 2. / 3.,
725
726 3. / 1., 3. / 1., 1.};
727
729 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
730
732 dd_f, t_S_sym, 1);
733
734
736
739 << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
740 if (nrm2_t_dd_a >
eps)
742 "This norm should be zero");
743 }
744
745
746 {
747
748 std::array<double, 9>
a{5., 4., 0.,
749
750 4., 5., 0.,
751
752 0., 0., 9};
753
754 auto tuple = run_lapack(a, swap01);
755
756 auto &t_eig_vecs = std::get<1>(tuple);
757 auto &t_eig_vals = std::get<2>(tuple);
758
759 constexpr double eps = 1e-10;
760
761 auto f = [](
double v) {
return v *
v; };
762 auto d_f = [](
double v) {
return 2 *
v; };
763 auto dd_f = [](
double v) {
return 2; };
764
766
767 1., 1. / 2., 1. / 3.,
768
769 2. / 1., 1., 2. / 3.,
770
771 3. / 1., 3. / 1., 1.};
772
774 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
775
777 dd_f, t_S_sym, 2);
778 print_ddg(t_dd, "test ");
779
781
784 << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
785 if (nrm2_t_dd_a >
eps)
787 "This norm should be zero");
788 }
789
790
791 {
792
793 std::array<double, 9>
a{2, 0., 0.,
794
795 0., 2, 0.,
796
797 0., 0., 2};
798
799 auto tuple = run_lapack(a);
800
801 auto &t_eig_vecs = std::get<1>(tuple);
802 auto &t_eig_vals = std::get<2>(tuple);
803
804 t_eig_vals(0) -= 1e-5;
805 t_eig_vals(2) += 1e-5;
806
807 constexpr double eps = 1e-7;
808
809 auto f = [](
double v) {
return exp(
v); };
810 auto d_f = [](
double v) {
return exp(
v); };
811 auto dd_f = [](
double v) {
return exp(
v); };
813
814 1., 1. / 2., 1. / 3.,
815
816 2. / 1., 1., 2. / 3.,
817
818 3. / 1., 3. / 1., 1.};
819
821 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
822
824 dd_f, t_S_sym, 3);
826 dd_f, t_S_sym, 1);
827
830 << "Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
831
834 << "Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
835
836 print_ddg(t_dd_1, "t_dd_1 ");
837 print_ddg(t_dd_2, "t_dd_2 ");
838
840 t_dd_3(
i,
j,
k,
l) = t_dd_1(
i,
j,
k,
l) - t_dd_2(
i,
j,
k,
l);
841
842 for (int ii = 0; ii != 3; ++ii)
843 for (int jj = 0; jj != 3; ++jj)
844 for (int kk = 0; kk != 3; ++kk)
845 for (int ll = 0; ll != 3; ++ll) {
846 constexpr double eps = 1e-4;
847 if (std::abs(t_dd_3(ii, jj, kk, ll)) >
eps)
849 << "Error " << ii << " " << jj << " " << kk << " " << ll
850 << " " << t_dd_1(ii, jj, kk, ll) << " "
851 << t_dd_2(ii, jj, kk, ll);
852 }
853
856 << "Direvarive approx. calculation minus code " << nrm2_t_dd_3;
857 if (nrm2_t_dd_3 >
eps)
859 "This norm should be zero");
860 }
861
862
863 {
864
865 std::array<double, 9>
a{5., 4., 0.,
866
867 4., 5., 0.,
868
869 0., 0., 9};
870
871 auto tuple = run_lapack(a, swap01);
872
873 auto &t_eig_vecs = std::get<1>(tuple);
874 auto &t_eig_vals = std::get<2>(tuple);
875
876 t_eig_vals(0) -= 1e-4;
877 t_eig_vals(2) += 1e-4;
878
879 constexpr double eps = 1e-4;
880
881 auto f = [](
double v) {
return v *
v; };
882 auto d_f = [](
double v) {
return 2 *
v; };
883 auto dd_f = [](
double v) {
return 2; };
884
886
887 1., 1. / 2., 1. / 3.,
888
889 2. / 1., 1., 2. / 3.,
890
891 3. / 1., 3. / 1., 1.};
892
894 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
895
897 dd_f, t_S_sym, 3);
899 dd_f, t_S_sym, 2);
900
903 << "Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
904
907 << "Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
908
909 print_ddg(t_dd_1, "t_dd_1 ");
910 print_ddg(t_dd_2, "t_dd_2 ");
911
913 t_dd_3(
i,
j,
k,
l) = t_dd_1(
i,
j,
k,
l) - t_dd_2(
i,
j,
k,
l);
914
915 for (int ii = 0; ii != 3; ++ii)
916 for (int jj = 0; jj != 3; ++jj)
917 for (int kk = 0; kk != 3; ++kk)
918 for (int ll = 0; ll != 3; ++ll) {
919 constexpr double eps = 1e-3;
920 if (std::abs(t_dd_3(ii, jj, kk, ll)) >
eps)
922 << "Error " << ii << " " << jj << " " << kk << " " << ll
923 << " " << t_dd_1(ii, jj, kk, ll) << " "
924 << t_dd_2(ii, jj, kk, ll);
925 }
926
929 << "Direvarive approx. calculation minus code " << nrm2_t_dd_3;
930 if (nrm2_t_dd_3 >
eps)
932 "This norm should be zero");
933 }
934
935
936 {
937
938 std::array<double, 9>
a{5., 4., 0.,
939
940 4., 5., 0.,
941
942 0., 0., 9};
943
944 auto tuple = run_lapack(a, swap01);
945
946 auto &t_eig_vecs = std::get<1>(tuple);
947 auto &t_eig_vals = std::get<2>(tuple);
948
949 constexpr double eps = 1e-4;
951
952 auto f = [](
double v) {
return pow(
v,
p); };
953 auto d_f = [](
double v) {
return p * pow(
v,
p - 1); };
954 auto dd_f = [](
double v) {
955 return p * (
p - 1) * pow(
v, std::max(0,
p - 2));
956 };
957
959
960 1., 1. / 2., 1. / 3.,
961
962 2. / 1., 1., 2. / 3.,
963
964 3. / 1., 3. / 1., 1.};
965
966
967
968
969 t_eig_vals(0) += 2e-5;
970 t_eig_vals(2) -= 2e-5;
972 dd_f, t_S, 3);
974 dd_f, t_S, 2);
975
978 << "Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
979
982 << "Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
983
984 print_ddg(t_dd_1, "t_dd_1 ");
985 print_ddg(t_dd_2, "t_dd_2 ");
986
988 t_dd_3(
i,
j,
k,
l) = t_dd_1(
i,
j,
k,
l) - t_dd_2(
i,
j,
k,
l);
989
990 for (int ii = 0; ii != 3; ++ii)
991 for (int jj = 0; jj != 3; ++jj)
992 for (int kk = 0; kk != 3; ++kk)
993 for (int ll = 0; ll != 3; ++ll) {
994 constexpr double eps = 1e-3;
995 if (std::abs(t_dd_3(ii, jj, kk, ll)) >
eps)
997 << "Error " << ii << " " << jj << " " << kk << " " << ll
998 << " " << t_dd_1(ii, jj, kk, ll) << " "
999 << t_dd_2(ii, jj, kk, ll) << " " << t_dd_3(ii, jj, kk, ll);
1000 }
1001
1004 << "Direvarive approx. calculation minus code " << nrm2_t_dd_3;
1005 if (nrm2_t_dd_3 >
eps)
1007 "This norm should be zero");
1008 }
1009
1010
1011 {
1012
1013 std::array<double, 9>
a{1., 0.1, -0.5,
1014
1015 0.1, 2., 0.,
1016
1017 -0.5, 0., 3.};
1018
1019 auto tuple = run_lapack(a);
1020
1021 auto &t_eig_vecs = std::get<1>(tuple);
1022 auto &t_eig_vals = std::get<2>(tuple);
1023
1024 auto f = [](
double v) {
return exp(
v); };
1025 auto d_f = [](
double v) {
return exp(
v); };
1026 auto dd_f = [](
double v) {
return exp(
v); };
1027
1029
1030 1., 1. / 2., 1. / 3.,
1031
1032 2. / 1., 1., 2. / 3.,
1033
1034 3. / 1., 3. / 1., 1.};
1035
1037 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
1038
1039 MOFEM_LOG(
"ATOM_TEST", Sev::inform) <<
"Start";
1040 for (int ii = 0; ii != 1000; ++ii) {
1043 dd_f, t_S_sym, 3);
1044 std::ignore = t_d;
1045 std::ignore = t_dd;
1046 }
1047 MOFEM_LOG(
"ATOM_TEST", Sev::inform) <<
"End";
1048 }
1049
1050
1051
1052 auto run_lapack_2d = [](
auto &
a) {
1053 int info;
1054 double wkopt;
1056
1058
1060
1062
1063
1064 int lwork = -1;
1065 info =
lapack_dsyev(
'V',
'U', 2,
a.data(), 2, w, &wkopt, lwork);
1066 if (info > 0)
1067 THROW_MESSAGE(
"The algorithm failed to compute eigenvalues.");
1069 std::vector<double> work(lwork);
1070
1071 info =
lapack_dsyev(
'V',
'U', 2,
a.data(), 2, w, &*work.begin(), lwork);
1072 if (info > 0)
1073 THROW_MESSAGE(
"The algorithm failed to compute eigenvalues.");
1074
1076
1077 a[0 * 2 + 0],
a[0 * 2 + 1],
1078
1079 a[1 * 2 + 0],
a[1 * 2 + 1]};
1080
1082
1083 return std::make_tuple(t_a, t_eig_vecs, t_eig_vals);
1084 };
1085
1086
1087 {
1088
1089 std::array<double, 9>
a{1., 0.1,
1090
1091 0.1, 2.};
1092
1093 auto tuple = run_lapack_2d(a);
1094 auto &t_A = std::get<0>(tuple);
1095 auto &t_eig_vecs = std::get<1>(tuple);
1096 auto &t_eig_vals = std::get<2>(tuple);
1097
1098 auto f = [](
double v) {
return v *
v; };
1099 auto d_f = [](
double v) {
return 2 *
v; };
1100 auto dd_f = [](
double v) {
return 2; };
1101
1102 constexpr double eps = 1e-6;
1103
1108
1109
1110 {
1113 t_a(
i,
j) = t_b(
i,
j) - t_A(
i,
k) * t_A(
k,
j);
1115 auto norm2_t_a = t_a(
i,
j) * t_a(
i,
j);
1117 << "Result should be matrix times matrix " << norm2_t_a;
1118 if (norm2_t_a >
eps)
1120 "This norm should be zero");
1121 }
1122
1123
1124 {
1129 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
1130 if (nrm2_t_d_a >
eps)
1132 "This norm should be zero");
1133 }
1134
1135
1136 {
1138
1139 1., 1. / 2,
1140
1141 2. / 2., 1.};
1142
1143
1144
1145
1146
1147
1149 dd_f, t_S, 2);
1151
1152 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd_a";
1153 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd";
1154
1157 << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
1158 if (nrm2_t_dd_a >
eps)
1160 "This norm should be zero");
1161 }
1162 }
1163
1164
1165 {
1166
1167 std::array<double, 9>
a{2., 0,
1168
1169 0, 2.};
1170
1171 auto tuple = run_lapack_2d(a);
1172 auto &t_A = std::get<0>(tuple);
1173 auto &t_eig_vecs = std::get<1>(tuple);
1174 auto &t_eig_vals = std::get<2>(tuple);
1175
1176 auto f = [](
double v) {
return v *
v; };
1177 auto d_f = [](
double v) {
return 2 *
v; };
1178 auto dd_f = [](
double v) {
return 2; };
1179
1180 constexpr double eps = 1e-6;
1181
1186
1187
1188 {
1191 t_a(
i,
j) = t_b(
i,
j) - t_A(
i,
k) * t_A(
k,
j);
1193 auto norm2_t_a = t_a(
i,
j) * t_a(
i,
j);
1195 << "Result should be matrix times matrix " << norm2_t_a;
1196 if (norm2_t_a >
eps)
1198 "This norm should be zero");
1199 }
1200
1201
1202 {
1207 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
1208 if (nrm2_t_d_a >
eps)
1210 "This norm should be zero");
1211 }
1212
1213
1214 {
1216
1217 1., 1. / 2,
1218
1219 2. / 2., 1.};
1220
1224 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
1225
1227 dd_f, t_S_sym, 1);
1229
1230 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd_a";
1231 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd";
1232
1235 << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
1236 if (nrm2_t_dd_a >
eps)
1238 "This norm should be zero");
1239 }
1240 }
1241 }
1243
1245}
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
@ MOFEM_ATOM_TEST_INVALID
#define CHKERR
Inline error check.
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG_ATTRIBUTES(channel, bit)
Add attributes to channel.
const double v
phase velocity of light in medium (cm/ns)
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
auto get_diff_matrix(T1 &t_d, const FTensor::Number< DIM > &)
auto get_diff2_matrix2(T1 &t_s, T2 &t_dd, const FTensor::Number< DIM > &)
auto get_norm_t4(T &t, const FTensor::Number< DIM > &)
auto get_diff_matrix2(T1 &t_a, T2 &t_d, const FTensor::Number< DIM > &)
FTensor::Ddg< double, 3, 3 > getDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
FTensor::Tensor2_symmetric< double, 3 > getMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f)
Get the Mat object.
FTensor::Ddg< double, 3, 3 > getDiffDiffMat(Val< double, 3 > &t_val, Vec< double, 3 > &t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, FTensor::Tensor2< double, 3, 3 > &t_S, const int nb)
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmSelf()
Get the strm self object.