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 auto &t_a = std::get<0>(tuple);
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 auto d_f = [](
double v) {
return exp(
v); };
378 auto dd_f = [](
double v) {
return exp(
v); };
379
382 t_c(
i,
j) -= t_b(
i,
j);
384
385 auto norm2_t_c = t_c(
i,
j) * t_c(
i,
j);
387 << "Reconstruct mat difference with lapack eigen valyes and "
388 "vectors "
389 << norm2_t_c;
390
391 constexpr double eps = 1e-8;
392 if (fabs(norm2_t_c) >
eps)
394 "Matrix not reeconstructed");
395 }
396
400
401
402 {
403 auto f = [](
double v) {
return v; };
404 auto d_f = [](
double v) {
return 1; };
405 auto dd_f = [](
double v) {
return 0; };
406
407 constexpr double eps = 1e-10;
408 {
410 t_b(
i,
j) -= (t_A(
i,
j) || t_A(
j,
i)) / 2;
411 auto norm2_t_b = t_b(
i,
j) * t_b(
i,
j);
413 << "Result should be matrix itself " << norm2_t_b;
416 "This norm should be zero");
417 }
418
419 {
420
423
424 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
425 print_ddg(t_d_a, "hand ");
426 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
427 print_ddg(t_d, "code ");
428
431 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
432 if (nrm2_t_d_a >
eps)
434 "This norm should be zero");
435 }
436
437 {
439
440 1., 0., 0.,
441
442 0., 1., 0.,
443
444 0., 0., 1.};
445
447 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
448
449 auto t_dd =
451
453 MOFEM_LOG(
"ATOM_TEST", Sev::inform) <<
"norm2_t_dd " << norm2_t_dd;
454 if (norm2_t_dd >
eps)
456 "This norm should be zero");
457 }
458 }
459
460
461 {
462 auto f = [](
double v) {
return v *
v; };
463 auto d_f = [](
double v) {
return 2 *
v; };
464 auto dd_f = [](
double v) {
return 2; };
465
466 constexpr double eps = 1e-9;
467
468
469 {
472 t_a(
i,
j) = t_b(
i,
j) - t_A(
i,
k) * t_A(
k,
j);
474 auto norm2_t_a = t_a(
i,
j) * t_a(
i,
j);
476 << "Result should be matrix times matrix " << norm2_t_a;
479 "This norm should be zero");
480 }
481
482
483 {
485 print_ddg_direction(t_d, 0, 2);
487 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
488 print_ddg(t_d_a, "hand ");
489 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
490 print_ddg(t_d, "code ");
493 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
494 if (nrm2_t_d_a >
eps)
496 "This norm should be zero");
497 }
498
499
500 {
502
503 1., 1. / 2., 1. / 3.,
504
505 2. / 2., 1., 2. / 3.,
506
507 3. / 2., 1., 3. / 3.};
508
510 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
511
512 auto t_dd =
515
516 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd_a";
517 print_ddg(t_dd_a, "hand ");
518 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd";
519 print_ddg(t_dd, "code ");
520
523 << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
524 if (nrm2_t_dd_a >
eps)
526 "This norm should be zero");
527 }
528 }
529 }
530
531
532 {
533
534 std::array<double, 9>
a{5., 4., 0,
535
536 4., 5, 0.,
537
538 0.0, 0., 9};
539
540 auto tuple = run_lapack(
a, swap01);
541 auto &t_a = std::get<0>(tuple);
542 auto &t_eig_vecs = std::get<1>(tuple);
543 auto &t_eig_vals = std::get<2>(tuple);
544
545 auto f = [](
double v) {
return v; };
546 auto d_f = [](
double v) {
return 1; };
547 auto dd_f = [](
double v) {
return 0; };
548
549 constexpr double eps = 1e-10;
550 {
552 t_b(
i,
j) -= (t_a(
i,
j) || t_a(
j,
i)) / 2;
553 auto norm2_t_b = t_b(
i,
j) * t_b(
i,
j);
555 << "Result should be matrix itself " << norm2_t_b;
558 "This norm should be zero");
559 }
560
561 {
564 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
565 print_ddg(t_d_a, "hand ");
566 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
567 print_ddg(t_d, "code ");
570 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
571 if (nrm2_t_d_a >
eps)
573 "This norm should be zero");
574 }
575
576 {
577 auto f = [](
double v) {
return v *
v; };
578 auto d_f = [](
double v) {
return 2 *
v; };
579 auto dd_f = [](
double v) {
return 2; };
582 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
583 print_ddg(t_d_a, "hand ");
584 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
585 print_ddg(t_d, "code ");
588 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
589 if (nrm2_t_d_a >
eps)
591 "This norm should be zero");
592 }
593 }
594
595
596 {
597
598 std::array<double, 9>
a{4., 0., 0,
599
600 0., 4., 0.,
601
602 0.0, 0., 4.};
603
604 auto f = [](
double v) {
return v; };
605 auto d_f = [](
double v) {
return 1; };
606 auto dd_f = [](
double v) {
return 0; };
607
608 auto tuple = run_lapack(
a);
609 auto &t_a = std::get<0>(tuple);
610 auto &t_eig_vecs = std::get<1>(tuple);
611 auto &t_eig_vals = std::get<2>(tuple);
612
613 constexpr double eps = 1e-10;
614 {
616 t_b(
i,
j) -= (t_a(
i,
j) || t_a(
j,
i)) / 2;
617 auto norm2_t_b = t_b(
i,
j) * t_b(
i,
j);
619 << "Result should be matrix itself " << norm2_t_b;
622 "This norm should be zero");
623 }
624
625 {
628 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
629 print_ddg(t_d_a, "hand ");
630 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
631 print_ddg(t_d, "code ");
634 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
635 if (nrm2_t_d_a >
eps)
637 "This norm should be zero");
638 }
639
640 {
641 auto f = [](
double v) {
return v *
v; };
642 auto d_f = [](
double v) {
return 2 *
v; };
643 auto dd_f = [](
double v) {
return 2; };
646 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
647 print_ddg(t_d_a, "hand ");
648 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
649 print_ddg(t_d, "code ");
652 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
653 if (nrm2_t_d_a >
eps)
655 "This norm should be zero");
656 }
657 }
658
659
660 {
661
662 std::array<double, 9>
a{0.1, 0., 0.,
663
664 0., 0.1, 0.,
665
666 0., 0., 0.1};
667
668 auto tuple = run_lapack(
a);
669 auto &t_a = std::get<0>(tuple);
670 auto &t_eig_vecs = std::get<1>(tuple);
671 auto &t_eig_vals = std::get<2>(tuple);
672
673 t_eig_vals(0) -= 1e-4;
674 t_eig_vals(2) += 1e-4;
675
676 constexpr double eps = 1e-10;
677
678 auto f = [](
double v) {
return v; };
679 auto d_f = [](
double v) {
return 1; };
680 auto dd_f = [](
double v) {
return 0; };
681
683
684 1., 1. / 2., 1. / 3.,
685
686 2. / 2., 1., 2. / 3.,
687
688 3. / 2., 1., 3. / 3.};
689
691 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
692
694 dd_f, t_S_sym, 1);
695
696 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd";
697 print_ddg(t_dd, "test ");
698
701 << "Direvarive hand calculation minus code " << nrm2_t_dd;
704 "This norm should be zero");
705 }
706
707
708 {
709
710 std::array<double, 9>
a{2, 0., 0.,
711
712 0., 2, 0.,
713
714 0., 0., 2};
715
716 auto tuple = run_lapack(
a);
717 auto &t_a = std::get<0>(tuple);
718 auto &t_eig_vecs = std::get<1>(tuple);
719 auto &t_eig_vals = std::get<2>(tuple);
720
721 constexpr double eps = 1e-10;
722
723 auto f = [](
double v) {
return v *
v; };
724 auto d_f = [](
double v) {
return 2 *
v; };
725 auto dd_f = [](
double v) {
return 2; };
727
728 1., 1. / 2., 1. / 3.,
729
730 2. / 1., 1., 2. / 3.,
731
732 3. / 1., 3. / 1., 1.};
733
735 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
736
738 dd_f, t_S_sym, 1);
739
740
742
745 << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
746 if (nrm2_t_dd_a >
eps)
748 "This norm should be zero");
749 }
750
751
752 {
753
754 std::array<double, 9>
a{5., 4., 0.,
755
756 4., 5., 0.,
757
758 0., 0., 9};
759
760 auto tuple = run_lapack(
a, swap01);
761 auto &t_a = std::get<0>(tuple);
762 auto &t_eig_vecs = std::get<1>(tuple);
763 auto &t_eig_vals = std::get<2>(tuple);
764
765 constexpr double eps = 1e-10;
766
767 auto f = [](
double v) {
return v *
v; };
768 auto d_f = [](
double v) {
return 2 *
v; };
769 auto dd_f = [](
double v) {
return 2; };
770
772
773 1., 1. / 2., 1. / 3.,
774
775 2. / 1., 1., 2. / 3.,
776
777 3. / 1., 3. / 1., 1.};
778
780 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
781
783 dd_f, t_S_sym, 2);
784 print_ddg(t_dd, "test ");
785
787
790 << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
791 if (nrm2_t_dd_a >
eps)
793 "This norm should be zero");
794 }
795
796
797 {
798
799 std::array<double, 9>
a{2, 0., 0.,
800
801 0., 2, 0.,
802
803 0., 0., 2};
804
805 auto tuple = run_lapack(
a);
806 auto &t_a = std::get<0>(tuple);
807 auto &t_eig_vecs = std::get<1>(tuple);
808 auto &t_eig_vals = std::get<2>(tuple);
809
810 t_eig_vals(0) -= 1e-5;
811 t_eig_vals(2) += 1e-5;
812
813 constexpr double eps = 1e-7;
814
815 auto f = [](
double v) {
return exp(
v); };
816 auto d_f = [](
double v) {
return exp(
v); };
817 auto dd_f = [](
double v) {
return exp(
v); };
819
820 1., 1. / 2., 1. / 3.,
821
822 2. / 1., 1., 2. / 3.,
823
824 3. / 1., 3. / 1., 1.};
825
827 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
828
830 dd_f, t_S_sym, 3);
832 dd_f, t_S_sym, 1);
833
836 << "Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
837
840 << "Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
841
842 print_ddg(t_dd_1, "t_dd_1 ");
843 print_ddg(t_dd_2, "t_dd_2 ");
844
846 t_dd_3(
i,
j,
k,
l) = t_dd_1(
i,
j,
k,
l) - t_dd_2(
i,
j,
k,
l);
847
848 for (int ii = 0; ii != 3; ++ii)
849 for (int jj = 0; jj != 3; ++jj)
850 for (int kk = 0; kk != 3; ++kk)
851 for (int ll = 0; ll != 3; ++ll) {
852 constexpr double eps = 1e-4;
853 if (std::abs(t_dd_3(ii, jj, kk, ll)) >
eps)
855 << "Error " << ii << " " << jj << " " << kk << " " << ll
856 << " " << t_dd_1(ii, jj, kk, ll) << " "
857 << t_dd_2(ii, jj, kk, ll);
858 }
859
862 << "Direvarive approx. calculation minus code " << nrm2_t_dd_3;
863 if (nrm2_t_dd_3 >
eps)
865 "This norm should be zero");
866 }
867
868
869 {
870
871 std::array<double, 9>
a{5., 4., 0.,
872
873 4., 5., 0.,
874
875 0., 0., 9};
876
877 auto tuple = run_lapack(
a, swap01);
878 auto &t_a = std::get<0>(tuple);
879 auto &t_eig_vecs = std::get<1>(tuple);
880 auto &t_eig_vals = std::get<2>(tuple);
881
882 t_eig_vals(0) -= 1e-4;
883 t_eig_vals(2) += 1e-4;
884
885 constexpr double eps = 1e-4;
886
887 auto f = [](
double v) {
return v *
v; };
888 auto d_f = [](
double v) {
return 2 *
v; };
889 auto dd_f = [](
double v) {
return 2; };
890
892
893 1., 1. / 2., 1. / 3.,
894
895 2. / 1., 1., 2. / 3.,
896
897 3. / 1., 3. / 1., 1.};
898
900 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
901
903 dd_f, t_S_sym, 3);
905 dd_f, t_S_sym, 2);
906
909 << "Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
910
913 << "Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
914
915 print_ddg(t_dd_1, "t_dd_1 ");
916 print_ddg(t_dd_2, "t_dd_2 ");
917
919 t_dd_3(
i,
j,
k,
l) = t_dd_1(
i,
j,
k,
l) - t_dd_2(
i,
j,
k,
l);
920
921 for (int ii = 0; ii != 3; ++ii)
922 for (int jj = 0; jj != 3; ++jj)
923 for (int kk = 0; kk != 3; ++kk)
924 for (int ll = 0; ll != 3; ++ll) {
925 constexpr double eps = 1e-3;
926 if (std::abs(t_dd_3(ii, jj, kk, ll)) >
eps)
928 << "Error " << ii << " " << jj << " " << kk << " " << ll
929 << " " << t_dd_1(ii, jj, kk, ll) << " "
930 << t_dd_2(ii, jj, kk, ll);
931 }
932
935 << "Direvarive approx. calculation minus code " << nrm2_t_dd_3;
936 if (nrm2_t_dd_3 >
eps)
938 "This norm should be zero");
939 }
940
941
942 {
943
944 std::array<double, 9>
a{5., 4., 0.,
945
946 4., 5., 0.,
947
948 0., 0., 9};
949
950 auto tuple = run_lapack(
a, swap01);
951 auto &t_a = std::get<0>(tuple);
952 auto &t_eig_vecs = std::get<1>(tuple);
953 auto &t_eig_vals = std::get<2>(tuple);
954
955 constexpr double eps = 1e-4;
957
958 auto f = [](
double v) {
return pow(
v,
p); };
959 auto d_f = [](
double v) {
return p * pow(
v,
p - 1); };
960 auto dd_f = [](
double v) {
961 return p * (
p - 1) * pow(
v, std::max(0,
p - 2));
962 };
963
965
966 1., 1. / 2., 1. / 3.,
967
968 2. / 1., 1., 2. / 3.,
969
970 3. / 1., 3. / 1., 1.};
971
972
973
974
975 t_eig_vals(0) += 2e-5;
976 t_eig_vals(2) -= 2e-5;
978 dd_f, t_S, 3);
980 dd_f, t_S, 2);
981
984 << "Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
985
988 << "Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
989
990 print_ddg(t_dd_1, "t_dd_1 ");
991 print_ddg(t_dd_2, "t_dd_2 ");
992
994 t_dd_3(
i,
j,
k,
l) = t_dd_1(
i,
j,
k,
l) - t_dd_2(
i,
j,
k,
l);
995
996 for (int ii = 0; ii != 3; ++ii)
997 for (int jj = 0; jj != 3; ++jj)
998 for (int kk = 0; kk != 3; ++kk)
999 for (int ll = 0; ll != 3; ++ll) {
1000 constexpr double eps = 1e-3;
1001 if (std::abs(t_dd_3(ii, jj, kk, ll)) >
eps)
1003 << "Error " << ii << " " << jj << " " << kk << " " << ll
1004 << " " << t_dd_1(ii, jj, kk, ll) << " "
1005 << t_dd_2(ii, jj, kk, ll) << " " << t_dd_3(ii, jj, kk, ll);
1006 }
1007
1010 << "Direvarive approx. calculation minus code " << nrm2_t_dd_3;
1011 if (nrm2_t_dd_3 >
eps)
1013 "This norm should be zero");
1014 }
1015
1016
1017 {
1018
1019 std::array<double, 9>
a{1., 0.1, -0.5,
1020
1021 0.1, 2., 0.,
1022
1023 -0.5, 0., 3.};
1024
1025 auto tuple = run_lapack(
a);
1026 auto &t_a = std::get<0>(tuple);
1027 auto &t_eig_vecs = std::get<1>(tuple);
1028 auto &t_eig_vals = std::get<2>(tuple);
1029
1030 auto f = [](
double v) {
return exp(
v); };
1031 auto d_f = [](
double v) {
return exp(
v); };
1032 auto dd_f = [](
double v) {
return exp(
v); };
1033
1035
1036 1., 1. / 2., 1. / 3.,
1037
1038 2. / 1., 1., 2. / 3.,
1039
1040 3. / 1., 3. / 1., 1.};
1041
1043 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
1044
1045 MOFEM_LOG(
"ATOM_TEST", Sev::inform) <<
"Start";
1046 for (int ii = 0; ii != 1000; ++ii) {
1049 dd_f, t_S_sym, 3);
1050 }
1051 MOFEM_LOG(
"ATOM_TEST", Sev::inform) <<
"End";
1052 }
1053
1054
1055
1056 auto run_lapack_2d = [](
auto &
a) {
1057 int info;
1058 double wkopt;
1060
1062
1064
1066
1067
1068 int lwork = -1;
1069 info =
lapack_dsyev(
'V',
'U', 2,
a.data(), 2, w, &wkopt, lwork);
1070 if (info > 0)
1071 THROW_MESSAGE(
"The algorithm failed to compute eigenvalues.");
1073 std::vector<double> work(lwork);
1074
1075 info =
lapack_dsyev(
'V',
'U', 2,
a.data(), 2, w, &*work.begin(), lwork);
1076 if (info > 0)
1077 THROW_MESSAGE(
"The algorithm failed to compute eigenvalues.");
1078
1080
1081 a[0 * 2 + 0],
a[0 * 2 + 1],
1082
1083 a[1 * 2 + 0],
a[1 * 2 + 1]};
1084
1086
1087 return std::make_tuple(t_a, t_eig_vecs, t_eig_vals);
1088 };
1089
1090
1091 {
1092
1093 std::array<double, 9>
a{1., 0.1,
1094
1095 0.1, 2.};
1096
1097 auto tuple = run_lapack_2d(
a);
1098 auto &t_A = std::get<0>(tuple);
1099 auto &t_eig_vecs = std::get<1>(tuple);
1100 auto &t_eig_vals = std::get<2>(tuple);
1101
1102 auto f = [](
double v) {
return v *
v; };
1103 auto d_f = [](
double v) {
return 2 *
v; };
1104 auto dd_f = [](
double v) {
return 2; };
1105
1106 constexpr double eps = 1e-6;
1107
1112
1113
1114 {
1117 t_a(
i,
j) = t_b(
i,
j) - t_A(
i,
k) * t_A(
k,
j);
1119 auto norm2_t_a = t_a(
i,
j) * t_a(
i,
j);
1121 << "Result should be matrix times matrix " << norm2_t_a;
1122 if (norm2_t_a >
eps)
1124 "This norm should be zero");
1125 }
1126
1127
1128 {
1133 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
1134 if (nrm2_t_d_a >
eps)
1136 "This norm should be zero");
1137 }
1138
1139
1140 {
1142
1143 1., 1. / 2,
1144
1145 2. / 2., 1.};
1146
1147
1148
1149
1150
1151
1153 dd_f, t_S, 2);
1155
1156 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd_a";
1157 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd";
1158
1161 << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
1162 if (nrm2_t_dd_a >
eps)
1164 "This norm should be zero");
1165 }
1166 }
1167
1168
1169 {
1170
1171 std::array<double, 9>
a{2., 0,
1172
1173 0, 2.};
1174
1175 auto tuple = run_lapack_2d(
a);
1176 auto &t_A = std::get<0>(tuple);
1177 auto &t_eig_vecs = std::get<1>(tuple);
1178 auto &t_eig_vals = std::get<2>(tuple);
1179
1180 auto f = [](
double v) {
return v *
v; };
1181 auto d_f = [](
double v) {
return 2 *
v; };
1182 auto dd_f = [](
double v) {
return 2; };
1183
1184 constexpr double eps = 1e-6;
1185
1190
1191
1192 {
1195 t_a(
i,
j) = t_b(
i,
j) - t_A(
i,
k) * t_A(
k,
j);
1197 auto norm2_t_a = t_a(
i,
j) * t_a(
i,
j);
1199 << "Result should be matrix times matrix " << norm2_t_a;
1200 if (norm2_t_a >
eps)
1202 "This norm should be zero");
1203 }
1204
1205
1206 {
1211 << "Direvarive hand calculation minus code " << nrm2_t_d_a;
1212 if (nrm2_t_d_a >
eps)
1214 "This norm should be zero");
1215 }
1216
1217
1218 {
1220
1221 1., 1. / 2,
1222
1223 2. / 2., 1.};
1224
1228 t_S_sym(
i,
j) = t_S(
i,
j) || t_S(
j,
i);
1229
1231 dd_f, t_S_sym, 1);
1233
1234 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd_a";
1235 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd";
1236
1239 << "Direvarive hand calculation minus code " << nrm2_t_dd_a;
1240 if (nrm2_t_dd_a >
eps)
1242 "This norm should be zero");
1243 }
1244 }
1245 }
1247
1249}
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.