186 auto core_log = logging::core::get();
188 LogManager::createSink(LogManager::getStrmSelf(),
"ATOM_TEST"));
189 LogManager::setLog(
"ATOM_TEST");
190 BOOST_LOG_SCOPED_THREAD_ATTR(
"Timeline", attrs::timer());
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
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);
219 for (
int ii = 0; ii != 3; ++ii)
220 for (
int jj = 0; jj != 3; ++jj)
222 << ii + 1 <<
" " << jj + 1 <<
" : " << t(ii, jj);
225 enum swap { swap12, swap01 };
226 auto run_lapack = [](
auto &a, swap s = swap12) {
241 info =
lapack_dsyev(
'V',
'U', 3, a.data(), 3,
w, &wkopt, lwork);
243 THROW_MESSAGE(
"The algorithm failed to compute eigenvalues.");
245 std::vector<double> work(lwork);
247 info =
lapack_dsyev(
'V',
'U', 3, a.data(), 3,
w, &*work.begin(), lwork);
249 THROW_MESSAGE(
"The algorithm failed to compute eigenvalues.");
254 a[0 * 3 + 0], a[0 * 3 + 1], a[0 * 3 + 2],
256 a[2 * 3 + 0], a[2 * 3 + 1], a[2 * 3 + 2],
258 a[1 * 3 + 0], a[1 * 3 + 1], a[1 * 3 + 2]};
261 return std::make_tuple(t_a, t_eig_vec, t_eig_vals);
266 a[1 * 3 + 0], a[1 * 3 + 1], a[1 * 3 + 2],
268 a[0 * 3 + 0], a[0 * 3 + 1], a[0 * 3 + 2],
270 a[2 * 3 + 0], a[2 * 3 + 1], a[2 * 3 + 2],
275 return std::make_tuple(t_a, t_eig_vec, t_eig_vals);
291 0.969837, -0.0860972, 0.228042,
293 0.0790574, 0.996073, 0.0398449,
295 -0.230577, -0.0206147, 0.972836};
307 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"Diff A";
312 auto f = [](
double v) {
return exp(v); };
313 auto d_f = [](
double v) {
return exp(v); };
314 auto dd_f = [](
double v) {
return exp(v); };
317 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"Reconstruct mat";
322 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"Diff";
323 print_ddg_direction(t_d, 0, 2);
327 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"Diff Diff";
328 print_ddg_direction(t_dd, 0, 2);
330 auto norm2_t_b = t_b(
i,
j) * t_b(
i,
j);
331 MOFEM_LOG(
"ATOM_TEST", Sev::inform) <<
"norm2_t_b " << norm2_t_b;
334 MOFEM_LOG(
"ATOM_TEST", Sev::inform) <<
"norm2_t_d " << norm2_t_d;
337 MOFEM_LOG(
"ATOM_TEST", Sev::inform) <<
"norm2_t_dd " << norm2_t_dd;
339 constexpr
double regression_t_b = 572.543;
340 constexpr
double regression_t_d = 859.939;
341 constexpr
double regression_t_dd = 859.939;
343 constexpr
double eps = 1e-2;
344 if (std::abs(norm2_t_b - regression_t_b) >
eps)
346 if (std::abs(norm2_t_d - regression_t_d) >
eps)
348 if (std::abs(norm2_t_dd - regression_t_dd) >
eps)
357 std::array<double, 9> a{1., 0.1, -0.5,
363 auto tuple = run_lapack(a);
364 auto &t_a = std::get<0>(tuple);
365 auto &t_eig_vec = std::get<1>(tuple);
366 auto &t_eig_vals = std::get<2>(tuple);
368 auto t_eig_val_diff =
369 (t_eig_vals(
i) - t_L(
i)) * (t_eig_vals(
i) - t_L(
i));
371 <<
"t_eig_val_diff " << t_eig_val_diff;
373 auto f = [](
double v) {
return exp(v); };
374 auto d_f = [](
double v) {
return exp(v); };
375 auto dd_f = [](
double v) {
return exp(v); };
379 t_c(
i,
j) -= t_b(
i,
j);
382 auto norm2_t_c = t_c(
i,
j) * t_c(
i,
j);
384 <<
"Reconstruct mat difference with lapack eigen valyes and "
388 constexpr
double eps = 1e-8;
389 if (fabs(norm2_t_c) >
eps)
391 "Matrix not reeconstructed");
396 t_one(
i,
j,
k,
l) = (t_kd(
i,
k) ^ t_kd(
j,
l)) / 4.;
400 auto f = [](
double v) {
return v; };
401 auto d_f = [](
double v) {
return 1; };
402 auto dd_f = [](
double v) {
return 0; };
404 constexpr
double eps = 1e-10;
407 t_b(
i,
j) -= (t_A(
i,
j) || t_A(
j,
i)) / 2;
408 auto norm2_t_b = t_b(
i,
j) * t_b(
i,
j);
410 <<
"Result should be matrix itself " << norm2_t_b;
413 "This norm should be zero");
421 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
422 print_ddg(t_d_a,
"hand ");
423 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
424 print_ddg(t_d,
"code ");
428 <<
"Direvarive hand calculation minus code " << nrm2_t_d_a;
429 if (nrm2_t_d_a >
eps)
431 "This norm should be zero");
447 MOFEM_LOG(
"ATOM_TEST", Sev::inform) <<
"norm2_t_dd " << norm2_t_dd;
448 if (norm2_t_dd >
eps)
450 "This norm should be zero");
456 auto f = [](
double v) {
return v * v; };
457 auto d_f = [](
double v) {
return 2 * v; };
458 auto dd_f = [](
double v) {
return 2; };
460 constexpr
double eps = 1e-9;
466 t_a(
i,
j) = t_b(
i,
j) - t_A(
i,
k) * t_A(
k,
j);
468 auto norm2_t_a = t_a(
i,
j) * t_a(
i,
j);
470 <<
"Result should be matrix times matrix " << norm2_t_a;
473 "This norm should be zero");
479 print_ddg_direction(t_d, 0, 2);
481 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
482 print_ddg(t_d_a,
"hand ");
483 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
484 print_ddg(t_d,
"code ");
487 <<
"Direvarive hand calculation minus code " << nrm2_t_d_a;
488 if (nrm2_t_d_a >
eps)
490 "This norm should be zero");
497 1., 1. / 2., 1. / 3.,
499 2. / 2., 1., 2. / 3.,
501 3. / 2., 1., 3. / 3.};
507 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd_a";
508 print_ddg(t_dd_a,
"hand ");
509 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd";
510 print_ddg(t_dd,
"code ");
514 <<
"Direvarive hand calculation minus code " << nrm2_t_dd_a;
515 if (nrm2_t_dd_a >
eps)
517 "This norm should be zero");
525 std::array<double, 9> a{5., 4., 0,
531 auto tuple = run_lapack(a, swap01);
532 auto &t_a = std::get<0>(tuple);
533 auto &t_eig_vecs = std::get<1>(tuple);
534 auto &t_eig_vals = std::get<2>(tuple);
536 auto f = [](
double v) {
return v; };
537 auto d_f = [](
double v) {
return 1; };
538 auto dd_f = [](
double v) {
return 0; };
540 constexpr
double eps = 1e-10;
543 t_b(
i,
j) -= (t_a(
i,
j) || t_a(
j,
i)) / 2;
544 auto norm2_t_b = t_b(
i,
j) * t_b(
i,
j);
546 <<
"Result should be matrix itself " << norm2_t_b;
549 "This norm should be zero");
555 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
556 print_ddg(t_d_a,
"hand ");
557 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
558 print_ddg(t_d,
"code ");
561 <<
"Direvarive hand calculation minus code " << nrm2_t_d_a;
562 if (nrm2_t_d_a >
eps)
564 "This norm should be zero");
568 auto f = [](
double v) {
return v * v; };
569 auto d_f = [](
double v) {
return 2 * v; };
570 auto dd_f = [](
double v) {
return 2; };
573 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
574 print_ddg(t_d_a,
"hand ");
575 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
576 print_ddg(t_d,
"code ");
579 <<
"Direvarive hand calculation minus code " << nrm2_t_d_a;
580 if (nrm2_t_d_a >
eps)
582 "This norm should be zero");
589 std::array<double, 9> a{4., 0., 0,
595 auto f = [](
double v) {
return v; };
596 auto d_f = [](
double v) {
return 1; };
597 auto dd_f = [](
double v) {
return 0; };
599 auto tuple = run_lapack(a);
600 auto &t_a = std::get<0>(tuple);
601 auto &t_eig_vecs = std::get<1>(tuple);
602 auto &t_eig_vals = std::get<2>(tuple);
604 constexpr
double eps = 1e-10;
607 t_b(
i,
j) -= (t_a(
i,
j) || t_a(
j,
i)) / 2;
608 auto norm2_t_b = t_b(
i,
j) * t_b(
i,
j);
610 <<
"Result should be matrix itself " << norm2_t_b;
613 "This norm should be zero");
619 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
620 print_ddg(t_d_a,
"hand ");
621 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
622 print_ddg(t_d,
"code ");
625 <<
"Direvarive hand calculation minus code " << nrm2_t_d_a;
626 if (nrm2_t_d_a >
eps)
628 "This norm should be zero");
632 auto f = [](
double v) {
return v * v; };
633 auto d_f = [](
double v) {
return 2 * v; };
634 auto dd_f = [](
double v) {
return 2; };
637 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d_a";
638 print_ddg(t_d_a,
"hand ");
639 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_d";
640 print_ddg(t_d,
"code ");
643 <<
"Direvarive hand calculation minus code " << nrm2_t_d_a;
644 if (nrm2_t_d_a >
eps)
646 "This norm should be zero");
653 std::array<double, 9> a{0.1, 0., 0.,
659 auto tuple = run_lapack(a);
660 auto &t_a = std::get<0>(tuple);
661 auto &t_eig_vecs = std::get<1>(tuple);
662 auto &t_eig_vals = std::get<2>(tuple);
664 t_eig_vals(0) -= 1e-4;
665 t_eig_vals(2) += 1e-4;
667 constexpr
double eps = 1e-10;
669 auto f = [](
double v) {
return v; };
670 auto d_f = [](
double v) {
return 1; };
671 auto dd_f = [](
double v) {
return 0; };
675 1., 1. / 2., 1. / 3.,
677 2. / 2., 1., 2. / 3.,
679 3. / 2., 1., 3. / 3.};
684 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd";
685 print_ddg(t_dd,
"test ");
689 <<
"Direvarive hand calculation minus code " << nrm2_t_dd;
692 "This norm should be zero");
698 std::array<double, 9> a{2, 0., 0.,
704 auto tuple = run_lapack(a);
705 auto &t_a = std::get<0>(tuple);
706 auto &t_eig_vecs = std::get<1>(tuple);
707 auto &t_eig_vals = std::get<2>(tuple);
709 constexpr
double eps = 1e-10;
711 auto f = [](
double v) {
return v * v; };
712 auto d_f = [](
double v) {
return 2 * v; };
713 auto dd_f = [](
double v) {
return 2; };
716 1., 1. / 2., 1. / 3.,
718 2. / 1., 1., 2. / 3.,
720 3. / 1., 3. / 1., 1.};
730 <<
"Direvarive hand calculation minus code " << nrm2_t_dd_a;
731 if (nrm2_t_dd_a >
eps)
733 "This norm should be zero");
739 std::array<double, 9> a{5., 4., 0.,
745 auto tuple = run_lapack(a, swap01);
746 auto &t_a = std::get<0>(tuple);
747 auto &t_eig_vecs = std::get<1>(tuple);
748 auto &t_eig_vals = std::get<2>(tuple);
750 constexpr
double eps = 1e-10;
752 auto f = [](
double v) {
return v * v; };
753 auto d_f = [](
double v) {
return 2 * v; };
754 auto dd_f = [](
double v) {
return 2; };
758 1., 1. / 2., 1. / 3.,
760 2. / 1., 1., 2. / 3.,
762 3. / 1., 3. / 1., 1.};
766 print_ddg(t_dd,
"test ");
772 <<
"Direvarive hand calculation minus code " << nrm2_t_dd_a;
773 if (nrm2_t_dd_a >
eps)
775 "This norm should be zero");
781 std::array<double, 9> a{2, 0., 0.,
787 auto tuple = run_lapack(a);
788 auto &t_a = std::get<0>(tuple);
789 auto &t_eig_vecs = std::get<1>(tuple);
790 auto &t_eig_vals = std::get<2>(tuple);
792 t_eig_vals(0) -= 1e-5;
793 t_eig_vals(2) += 1e-5;
795 constexpr
double eps = 1e-7;
797 auto f = [](
double v) {
return exp(v); };
798 auto d_f = [](
double v) {
return exp(v); };
799 auto dd_f = [](
double v) {
return exp(v); };
802 1., 1. / 2., 1. / 3.,
804 2. / 1., 1., 2. / 3.,
806 3. / 1., 3. / 1., 1.};
815 <<
"Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
819 <<
"Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
821 print_ddg(t_dd_1,
"t_dd_1 ");
822 print_ddg(t_dd_2,
"t_dd_2 ");
825 t_dd_3(
i,
j,
k,
l) = t_dd_1(
i,
j,
k,
l) - t_dd_2(
i,
j,
k,
l);
827 for (
int ii = 0; ii != 3; ++ii)
828 for (
int jj = 0; jj != 3; ++jj)
829 for (
int kk = 0; kk != 3; ++kk)
830 for (
int ll = 0; ll != 3; ++ll) {
831 constexpr
double eps = 1e-4;
832 if (std::abs(t_dd_3(ii, jj, kk, ll)) >
eps)
834 <<
"Error " << ii <<
" " << jj <<
" " << kk <<
" " << ll
835 <<
" " << t_dd_1(ii, jj, kk, ll) <<
" "
836 << t_dd_2(ii, jj, kk, ll);
841 <<
"Direvarive approx. calculation minus code " << nrm2_t_dd_3;
842 if (nrm2_t_dd_3 >
eps)
844 "This norm should be zero");
850 std::array<double, 9> a{5., 4., 0.,
856 auto tuple = run_lapack(a, swap01);
857 auto &t_a = std::get<0>(tuple);
858 auto &t_eig_vecs = std::get<1>(tuple);
859 auto &t_eig_vals = std::get<2>(tuple);
861 t_eig_vals(0) -= 1e-4;
862 t_eig_vals(2) += 1e-4;
864 constexpr
double eps = 1e-4;
866 auto f = [](
double v) {
return v * v; };
867 auto d_f = [](
double v) {
return 2 * v; };
868 auto dd_f = [](
double v) {
return 2; };
872 1., 1. / 2., 1. / 3.,
874 2. / 1., 1., 2. / 3.,
876 3. / 1., 3. / 1., 1.};
885 <<
"Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
889 <<
"Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
891 print_ddg(t_dd_1,
"t_dd_1 ");
892 print_ddg(t_dd_2,
"t_dd_2 ");
895 t_dd_3(
i,
j,
k,
l) = t_dd_1(
i,
j,
k,
l) - t_dd_2(
i,
j,
k,
l);
897 for (
int ii = 0; ii != 3; ++ii)
898 for (
int jj = 0; jj != 3; ++jj)
899 for (
int kk = 0; kk != 3; ++kk)
900 for (
int ll = 0; ll != 3; ++ll) {
901 constexpr
double eps = 1e-3;
902 if (std::abs(t_dd_3(ii, jj, kk, ll)) >
eps)
904 <<
"Error " << ii <<
" " << jj <<
" " << kk <<
" " << ll
905 <<
" " << t_dd_1(ii, jj, kk, ll) <<
" "
906 << t_dd_2(ii, jj, kk, ll);
911 <<
"Direvarive approx. calculation minus code " << nrm2_t_dd_3;
912 if (nrm2_t_dd_3 >
eps)
914 "This norm should be zero");
920 std::array<double, 9> a{5., 4., 0.,
926 auto tuple = run_lapack(a, swap01);
927 auto &t_a = std::get<0>(tuple);
928 auto &t_eig_vecs = std::get<1>(tuple);
929 auto &t_eig_vals = std::get<2>(tuple);
931 constexpr
double eps = 1e-4;
934 auto f = [](
double v) {
return pow(v,
p); };
935 auto d_f = [](
double v) {
return p * pow(v,
p - 1); };
936 auto dd_f = [](
double v) {
937 return p * (
p - 1) * pow(v, std::max(0,
p - 2));
942 1., 1. / 2., 1. / 3.,
944 2. / 1., 1., 2. / 3.,
946 3. / 1., 3. / 1., 1.};
948 t_eig_vals(0) += 2e-5;
949 t_eig_vals(2) -= 2e-5;
957 <<
"Direvarive nor t_dd_1 " << nrm2_t_dd_t1;
961 <<
"Direvarive norm t_dd_2 " << nrm2_t_dd_t2;
963 print_ddg(t_dd_1,
"t_dd_1 ");
964 print_ddg(t_dd_2,
"t_dd_2 ");
967 t_dd_3(
i,
j,
k,
l) = t_dd_1(
i,
j,
k,
l) - t_dd_2(
i,
j,
k,
l);
969 for (
int ii = 0; ii != 3; ++ii)
970 for (
int jj = 0; jj != 3; ++jj)
971 for (
int kk = 0; kk != 3; ++kk)
972 for (
int ll = 0; ll != 3; ++ll) {
973 constexpr
double eps = 1e-3;
974 if (std::abs(t_dd_3(ii, jj, kk, ll)) >
eps)
976 <<
"Error " << ii <<
" " << jj <<
" " << kk <<
" " << ll
977 <<
" " << t_dd_1(ii, jj, kk, ll) <<
" "
978 << t_dd_2(ii, jj, kk, ll) <<
" " << t_dd_3(ii, jj, kk, ll);
983 <<
"Direvarive approx. calculation minus code " << nrm2_t_dd_3;
984 if (nrm2_t_dd_3 >
eps)
986 "This norm should be zero");
992 std::array<double, 9> a{1., 0.1, -0.5,
998 auto tuple = run_lapack(a);
999 auto &t_a = std::get<0>(tuple);
1000 auto &t_eig_vecs = std::get<1>(tuple);
1001 auto &t_eig_vals = std::get<2>(tuple);
1003 auto f = [](
double v) {
return exp(v); };
1004 auto d_f = [](
double v) {
return exp(v); };
1005 auto dd_f = [](
double v) {
return exp(v); };
1009 1., 1. / 2., 1. / 3.,
1011 2. / 1., 1., 2. / 3.,
1013 3. / 1., 3. / 1., 1.};
1015 MOFEM_LOG(
"ATOM_TEST", Sev::inform) <<
"Start";
1016 for (
int ii = 0; ii != 1000; ++ii) {
1021 MOFEM_LOG(
"ATOM_TEST", Sev::inform) <<
"End";
1026 auto run_lapack_2d = [](
auto &a) {
1039 info =
lapack_dsyev(
'V',
'U', 2, a.data(), 2,
w, &wkopt, lwork);
1041 THROW_MESSAGE(
"The algorithm failed to compute eigenvalues.");
1043 std::vector<double> work(lwork);
1045 info =
lapack_dsyev(
'V',
'U', 2, a.data(), 2,
w, &*work.begin(), lwork);
1047 THROW_MESSAGE(
"The algorithm failed to compute eigenvalues.");
1051 a[0 * 2 + 0], a[0 * 2 + 1],
1053 a[1 * 2 + 0], a[1 * 2 + 1]};
1057 return std::make_tuple(t_a, t_eig_vecs, t_eig_vals);
1063 std::array<double, 9> a{1., 0.1,
1067 auto tuple = run_lapack_2d(a);
1068 auto &t_A = std::get<0>(tuple);
1069 auto &t_eig_vecs = std::get<1>(tuple);
1070 auto &t_eig_vals = std::get<2>(tuple);
1072 auto f = [](
double v) {
return v * v; };
1073 auto d_f = [](
double v) {
return 2 * v; };
1074 auto dd_f = [](
double v) {
return 2; };
1076 constexpr
double eps = 1e-6;
1087 t_a(
i,
j) = t_b(
i,
j) - t_A(
i,
k) * t_A(
k,
j);
1089 auto norm2_t_a = t_a(
i,
j) * t_a(
i,
j);
1091 <<
"Result should be matrix times matrix " << norm2_t_a;
1092 if (norm2_t_a >
eps)
1094 "This norm should be zero");
1103 <<
"Direvarive hand calculation minus code " << nrm2_t_d_a;
1104 if (nrm2_t_d_a >
eps)
1106 "This norm should be zero");
1121 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd_a";
1122 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd";
1126 <<
"Direvarive hand calculation minus code " << nrm2_t_dd_a;
1127 if (nrm2_t_dd_a >
eps)
1129 "This norm should be zero");
1136 std::array<double, 9> a{2., 0,
1140 auto tuple = run_lapack_2d(a);
1141 auto &t_A = std::get<0>(tuple);
1142 auto &t_eig_vecs = std::get<1>(tuple);
1143 auto &t_eig_vals = std::get<2>(tuple);
1145 auto f = [](
double v) {
return v * v; };
1146 auto d_f = [](
double v) {
return 2 * v; };
1147 auto dd_f = [](
double v) {
return 2; };
1149 constexpr
double eps = 1e-6;
1160 t_a(
i,
j) = t_b(
i,
j) - t_A(
i,
k) * t_A(
k,
j);
1162 auto norm2_t_a = t_a(
i,
j) * t_a(
i,
j);
1164 <<
"Result should be matrix times matrix " << norm2_t_a;
1165 if (norm2_t_a >
eps)
1167 "This norm should be zero");
1176 <<
"Direvarive hand calculation minus code " << nrm2_t_d_a;
1177 if (nrm2_t_d_a >
eps)
1179 "This norm should be zero");
1194 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd_a";
1195 MOFEM_LOG(
"ATOM_TEST", Sev::verbose) <<
"t_dd";
1199 <<
"Direvarive hand calculation minus code " << nrm2_t_dd_a;
1200 if (nrm2_t_dd_a >
eps)
1202 "This norm should be zero");