Skip to content

Commit fae598a

Browse files
jasonkayejasonkaye
committed
add test for convolution of non-square matrix-valued functions
1 parent d69258e commit fae598a

File tree

1 file changed

+143
-140
lines changed

1 file changed

+143
-140
lines changed

‎test/c++/imtime_ops.cpp‎

Lines changed: 143 additions & 140 deletions
Original file line numberDiff line numberDiff line change
@@ -632,84 +632,84 @@ TEST(imtime_ops, convolve_scalar_cmplx) {
632632
*/
633633
TEST(imtime_ops, convolve_matrix_real) {
634634

635-
double lambda = 1000; // DLR cutoff
635+
double lambda = 10; // DLR cutoff
636636
double eps = 1e-12; // DLR tolerance
637637

638-
double beta = 1000; // Inverse temperature
639-
int ntst = 10000; // # imag time test points
638+
double beta = 10; // Inverse temperature
639+
int ntst = 100; // # imag time test points
640640

641641
std::cout << fmt::format("eps = {:e}, Lambda = {:e}\n", eps, lambda);
642642

643-
int norb = 2; // Orbital dimensions
644-
645-
// Specify frequency of single exponentials to be used for f and g
646-
double omf = 0.1234, omg = -0.5678;
647-
648-
// Get DLR frequencies
649-
auto dlr_rf = build_dlr_rf(lambda, eps);
650-
651-
// Get DLR imaginary time object
652-
auto itops = imtime_ops(lambda, dlr_rf);
653-
654-
// Sample Green's function G at DLR imaginary time nodes
655-
int r = itops.rank();
656-
// [Q] Is this correct or just auto?
657-
auto const &dlr_it = itops.get_itnodes();
658-
auto f = nda::array<double, 3>(r, norb, norb);
659-
auto g = nda::array<double, 3>(r, norb, norb);
660-
for (int i = 0; i < r; ++i) { f(i, _, _) = k_it(dlr_it(i), omf, beta) / sqrt(1.0 * norb); };
661-
for (int i = 0; i < r; ++i) { g(i, _, _) = k_it(dlr_it(i), omg, beta) / sqrt(1.0 * norb); };
662-
663-
// Get DLR coefficients of f and g
664-
auto fc = itops.vals2coefs(f);
665-
auto gc = itops.vals2coefs(g);
666-
667-
// Get convolution and time-ordered convolution of f and g directly
668-
auto h = itops.convolve(beta, Fermion, fc, gc);
669-
auto ht = itops.convolve(beta, Fermion, fc, gc, TIME_ORDERED);
670-
671-
// Get convolution and time-ordered convolution of f and g by first forming
672-
// matrix of convolution by f and then applying it to g
673-
auto h2 = itops.convolve(itops.convmat(beta, Fermion, fc), g);
674-
auto ht2 = itops.convolve(itops.convmat(beta, Fermion, fc, TIME_ORDERED), g);
675-
676-
// Check that the two methods give the same result
677-
EXPECT_LT(max_element(abs(h - h2)), 1e-14);
678-
EXPECT_LT(max_element(abs(ht - ht2)), 1e-14);
679-
680-
// Check error of convolution and time-ordered convolution
681-
682-
auto hc = itops.vals2coefs(h); // DLR coefficients of h
683-
auto htc = itops.vals2coefs(ht); // DLR coefficients of ht
684-
685-
// Get test points in relative format
686-
auto ttst = eqptsrel(ntst);
687-
688-
// Compute error in imaginary time
689-
auto gtru = nda::array<double, 2>(norb, norb);
690-
auto gtst = nda::array<double, 2>(norb, norb);
691-
auto gttru = nda::array<double, 2>(norb, norb);
692-
auto gttst = nda::array<double, 2>(norb, norb);
693-
double errlinf = 0, errl2 = 0, errtlinf = 0, errtl2 = 0;
694-
for (int i = 0; i < ntst; ++i) {
695-
gtru = (k_it(ttst(i), omg, beta) - k_it(ttst(i), omf, beta)) / (omg - omf); // Exact result
696-
gttru = (k_it(0.0, omf, beta) * k_it(ttst(i), omg, beta) - k_it(ttst(i), omf, beta) * k_it(0.0, omg, beta)) / (omf - omg); // Exact result
697-
gtst = itops.coefs2eval(hc, ttst(i));
698-
gttst = itops.coefs2eval(htc, ttst(i));
699-
errlinf = std::max(errlinf, max_element(abs(gtru - gtst)));
700-
errtlinf = std::max(errtlinf, max_element(abs(gttru - gttst)));
701-
errl2 += pow(frobenius_norm(gtru - gtst), 2);
702-
errtl2 += pow(frobenius_norm(gttru - gttst), 2);
643+
int norb2 = 2; // Orbital dimensions
644+
for (int norb1 = 2; norb1 <= 3; ++norb1) {
645+
std::cout << fmt::format("norb1 = {}, norb2 = {}\n", norb1, norb2);
646+
647+
// Specify frequency of single exponentials to be used for f and g
648+
double omf = 0.1234, omg = -0.5678;
649+
650+
// Get DLR frequencies
651+
auto dlr_rf = build_dlr_rf(lambda, eps);
652+
653+
// Get DLR imaginary time object
654+
auto itops = imtime_ops(lambda, dlr_rf);
655+
656+
// Sample Green's function G at DLR imaginary time nodes
657+
int r = itops.rank();
658+
// [Q] Is this correct or just auto?
659+
auto const &dlr_it = itops.get_itnodes();
660+
auto f = nda::array<double, 3>(r, norb1, norb2);
661+
auto g = nda::array<double, 3>(r, norb2, norb2);
662+
for (int i = 0; i < r; ++i) { f(i, _, _) = k_it(dlr_it(i), omf, beta) / sqrt(1.0 * norb2); };
663+
for (int i = 0; i < r; ++i) { g(i, _, _) = k_it(dlr_it(i), omg, beta) / sqrt(1.0 * norb2); };
664+
665+
// Get DLR coefficients of f and g
666+
auto fc = itops.vals2coefs(f);
667+
auto gc = itops.vals2coefs(g);
668+
669+
// Get convolution and time-ordered convolution of f and g directly
670+
auto h = itops.convolve(beta, Fermion, fc, gc);
671+
auto ht = itops.convolve(beta, Fermion, fc, gc, TIME_ORDERED);
672+
673+
auto h2 = itops.convolve(itops.convmat(beta, Fermion, fc), g);
674+
auto ht2 = itops.convolve(itops.convmat(beta, Fermion, fc, TIME_ORDERED), g);
675+
676+
// Check that the two methods give the same result
677+
EXPECT_LT(max_element(abs(h - h2)), 1e-14);
678+
EXPECT_LT(max_element(abs(ht - ht2)), 1e-14);
679+
680+
// Check error of convolution and time-ordered convolution
681+
auto hc = itops.vals2coefs(h); // DLR coefficients of h
682+
auto htc = itops.vals2coefs(ht); // DLR coefficients of ht
683+
684+
// Get test points in relative format
685+
auto ttst = eqptsrel(ntst);
686+
687+
// Compute error in imaginary time
688+
auto htru = nda::array<double, 2>(norb1, norb2);
689+
auto htst = nda::array<double, 2>(norb1, norb2);
690+
auto httru = nda::array<double, 2>(norb1, norb2);
691+
auto httst = nda::array<double, 2>(norb1, norb2);
692+
double errlinf = 0, errl2 = 0, errtlinf = 0, errtl2 = 0;
693+
for (int i = 0; i < ntst; ++i) {
694+
htru = (k_it(ttst(i), omg, beta) - k_it(ttst(i), omf, beta)) / (omg - omf); // Exact result
695+
httru = (k_it(0.0, omf, beta) * k_it(ttst(i), omg, beta) - k_it(ttst(i), omf, beta) * k_it(0.0, omg, beta)) / (omf - omg); // Exact result
696+
htst = itops.coefs2eval(hc, ttst(i));
697+
httst = itops.coefs2eval(htc, ttst(i));
698+
errlinf = std::max(errlinf, max_element(abs(htru - htst)));
699+
errtlinf = std::max(errtlinf, max_element(abs(httru - httst)));
700+
errl2 += pow(frobenius_norm(htru - htst), 2);
701+
errtl2 += pow(frobenius_norm(httru - httst), 2);
702+
}
703+
errl2 = sqrt(errl2 / ntst);
704+
errtl2 = sqrt(errtl2 / ntst);
705+
706+
EXPECT_LT(errlinf, 23 * eps);
707+
EXPECT_LT(errl2, 7 * eps);
708+
EXPECT_LT(errtlinf, 23 * eps);
709+
EXPECT_LT(errtl2, 6 * eps);
710+
std::cout << fmt::format("Ordinary convolution: L^inf err = {:e}, L^2 err = {:e}\n", errlinf, errl2);
711+
std::cout << fmt::format("Time-ordered convolution: L^inf err = {:e}, L^2 err = {:e}\n", errtlinf, errtl2);
703712
}
704-
errl2 = sqrt(errl2 / ntst);
705-
errtl2 = sqrt(errtl2 / ntst);
706-
707-
EXPECT_LT(errlinf, 23 * eps);
708-
EXPECT_LT(errl2, 7 * eps);
709-
EXPECT_LT(errtlinf, 23 * eps);
710-
EXPECT_LT(errtl2, 6 * eps);
711-
std::cout << fmt::format("Ordinary convolution: L^inf err = {:e}, L^2 err = {:e}\n", errlinf, errl2);
712-
std::cout << fmt::format("Time-ordered convolution: L^inf err = {:e}, L^2 err = {:e}\n", errtlinf, errtl2);
713713
}
714714

715715
/**
@@ -730,76 +730,79 @@ TEST(imtime_ops, convolve_matrix_cmplx) {
730730

731731
std::cout << fmt::format("eps = {:e}, Lambda = {:e}\n", eps, lambda);
732732

733-
int norb = 2; // Orbital dimensions
734-
735-
// Specify frequency of single exponentials to be used for f and g
736-
double omf = 0.1234, omg = -0.5678;
737-
738-
// Get DLR frequencies
739-
auto dlr_rf = build_dlr_rf(lambda, eps);
740-
741-
// Get DLR imaginary time object
742-
auto itops = imtime_ops(lambda, dlr_rf);
743-
744-
// Sample Green's function G at DLR imaginary time nodes
745-
int r = itops.rank();
746-
// [Q] Is this correct or just auto?
747-
auto const &dlr_it = itops.get_itnodes();
748-
auto f = nda::array<dcomplex, 3>(r, norb, norb);
749-
auto g = nda::array<dcomplex, 3>(r, norb, norb);
750-
std::complex<double> c1 = (1.0 + 2.0i) / 3.0, c2 = (2.0 + 1.0i) / 3.0;
751-
for (int i = 0; i < r; ++i) { f(i, _, _) = c1 * k_it(dlr_it(i), omf, beta) / sqrt(1.0 * norb); };
752-
for (int i = 0; i < r; ++i) { g(i, _, _) = c2 * k_it(dlr_it(i), omg, beta) / sqrt(1.0 * norb); };
753-
754-
// Get DLR coefficients of f and g
755-
auto fc = itops.vals2coefs(f);
756-
auto gc = itops.vals2coefs(g);
757-
758-
// Get convolution and time-ordered convolution of f and g directly
759-
auto h = itops.convolve(beta, Fermion, fc, gc);
760-
auto ht = itops.convolve(beta, Fermion, fc, gc, TIME_ORDERED);
761-
762-
// Get convolution and time-ordered convolution of f and g by first forming
763-
// matrix of convolution by f and then applying it to g
764-
auto h2 = itops.convolve(itops.convmat(beta, Fermion, fc), g);
765-
auto ht2 = itops.convolve(itops.convmat(beta, Fermion, fc, TIME_ORDERED), g);
766-
767-
// Check that the two methods give the same result
768-
EXPECT_LT(max_element(abs(h - h2)), 1e-14);
769-
EXPECT_LT(max_element(abs(ht - ht2)), 1e-14);
770-
771-
// Check error of convolution and time-ordered convolution
772-
773-
auto hc = itops.vals2coefs(h); // DLR coefficients of h
774-
auto htc = itops.vals2coefs(ht); // DLR coefficients of ht
775-
776-
// Get test points in relative format
777-
auto ttst = eqptsrel(ntst);
733+
int norb2 = 2; // Orbital dimensions
734+
for (int norb1 = 2; norb1 <= 3; ++norb1) {
735+
std::cout << fmt::format("norb1 = {}, norb2 = {}\n", norb1, norb2);
736+
737+
// Specify frequency of single exponentials to be used for f and g
738+
double omf = 0.1234, omg = -0.5678;
739+
740+
// Get DLR frequencies
741+
auto dlr_rf = build_dlr_rf(lambda, eps);
742+
743+
// Get DLR imaginary time object
744+
auto itops = imtime_ops(lambda, dlr_rf);
745+
746+
// Sample Green's function G at DLR imaginary time nodes
747+
int r = itops.rank();
748+
// [Q] Is this correct or just auto?
749+
auto const &dlr_it = itops.get_itnodes();
750+
auto f = nda::array<dcomplex, 3>(r, norb1, norb2);
751+
auto g = nda::array<dcomplex, 3>(r, norb2, norb2);
752+
std::complex<double> c1 = (1.0 + 2.0i) / 3.0, c2 = (2.0 + 1.0i) / 3.0;
753+
for (int i = 0; i < r; ++i) { f(i, _, _) = c1 * k_it(dlr_it(i), omf, beta) / sqrt(1.0 * norb2); };
754+
for (int i = 0; i < r; ++i) { g(i, _, _) = c2 * k_it(dlr_it(i), omg, beta) / sqrt(1.0 * norb2); };
755+
756+
// Get DLR coefficients of f and g
757+
auto fc = itops.vals2coefs(f);
758+
auto gc = itops.vals2coefs(g);
759+
760+
// Get convolution and time-ordered convolution of f and g directly
761+
auto h = itops.convolve(beta, Fermion, fc, gc);
762+
auto ht = itops.convolve(beta, Fermion, fc, gc, TIME_ORDERED);
763+
764+
// Get convolution and time-ordered convolution of f and g by first forming
765+
// matrix of convolution by f and then applying it to g
766+
auto h2 = itops.convolve(itops.convmat(beta, Fermion, fc), g);
767+
auto ht2 = itops.convolve(itops.convmat(beta, Fermion, fc, TIME_ORDERED), g);
768+
769+
// Check that the two methods give the same result
770+
EXPECT_LT(max_element(abs(h - h2)), 1e-14);
771+
EXPECT_LT(max_element(abs(ht - ht2)), 1e-14);
772+
773+
// Check error of convolution and time-ordered convolution
774+
775+
auto hc = itops.vals2coefs(h); // DLR coefficients of h
776+
auto htc = itops.vals2coefs(ht); // DLR coefficients of ht
777+
778+
// Get test points in relative format
779+
auto ttst = eqptsrel(ntst);
780+
781+
// Compute error in imaginary time
782+
auto htru = nda::array<dcomplex, 2>(norb1, norb2);
783+
auto htst = nda::array<dcomplex, 2>(norb1, norb2);
784+
auto httru = nda::array<dcomplex, 2>(norb1, norb2);
785+
auto httst = nda::array<dcomplex, 2>(norb1, norb2);
786+
double errlinf = 0, errl2 = 0, errtlinf = 0, errtl2 = 0;
787+
for (int i = 0; i < ntst; ++i) {
788+
htru = c1 * c2 * (k_it(ttst(i), omg, beta) - k_it(ttst(i), omf, beta)) / (omg - omf); // Exact result
789+
httru =
790+
c1 * c2 * (k_it(0.0, omf, beta) * k_it(ttst(i), omg, beta) - k_it(ttst(i), omf, beta) * k_it(0.0, omg, beta)) / (omf - omg); // Exact result
791+
htst = itops.coefs2eval(hc, ttst(i));
792+
httst = itops.coefs2eval(htc, ttst(i));
793+
errlinf = std::max(errlinf, max_element(abs(htru - htst)));
794+
errl2 += pow(frobenius_norm(htru - htst), 2);
795+
errtlinf = std::max(errtlinf, max_element(abs(httru - httst)));
796+
errtl2 += pow(frobenius_norm(httru - httst), 2);
797+
}
778798

779-
// Compute error in imaginary time
780-
auto gtru = nda::array<dcomplex, 2>(norb, norb);
781-
auto gtst = nda::array<dcomplex, 2>(norb, norb);
782-
auto gttru = nda::array<dcomplex, 2>(norb, norb);
783-
auto gttst = nda::array<dcomplex, 2>(norb, norb);
784-
double errlinf = 0, errl2 = 0, errtlinf = 0, errtl2 = 0;
785-
for (int i = 0; i < ntst; ++i) {
786-
gtru = c1 * c2 * (k_it(ttst(i), omg, beta) - k_it(ttst(i), omf, beta)) / (omg - omf); // Exact result
787-
gttru =
788-
c1 * c2 * (k_it(0.0, omf, beta) * k_it(ttst(i), omg, beta) - k_it(ttst(i), omf, beta) * k_it(0.0, omg, beta)) / (omf - omg); // Exact result
789-
gtst = itops.coefs2eval(hc, ttst(i));
790-
gttst = itops.coefs2eval(htc, ttst(i));
791-
errlinf = std::max(errlinf, max_element(abs(gtru - gtst)));
792-
errl2 += pow(frobenius_norm(gtru - gtst), 2);
793-
errtlinf = std::max(errtlinf, max_element(abs(gttru - gttst)));
794-
errtl2 += pow(frobenius_norm(gttru - gttst), 2);
799+
EXPECT_LT(errlinf, 13 * eps);
800+
EXPECT_LT(errl2, 2 * eps);
801+
EXPECT_LT(errtlinf, 13 * eps);
802+
EXPECT_LT(errtl2, 2 * eps);
803+
std::cout << fmt::format("Ordinary convolution: L^inf err = {:e}, L^2 err = {:e}\n", errlinf, errl2);
804+
std::cout << fmt::format("Time-ordered convolution: L^inf err = {:e}, L^2 err = {:e}\n", errtlinf, errtl2);
795805
}
796-
797-
EXPECT_LT(errlinf, 13 * eps);
798-
EXPECT_LT(errl2, 2 * eps);
799-
EXPECT_LT(errtlinf, 13 * eps);
800-
EXPECT_LT(errtl2, 2 * eps);
801-
std::cout << fmt::format("Ordinary convolution: L^inf err = {:e}, L^2 err = {:e}\n", errlinf, errl2);
802-
std::cout << fmt::format("Time-ordered convolution: L^inf err = {:e}, L^2 err = {:e}\n", errtlinf, errtl2);
803806
}
804807

805808
/**

0 commit comments

Comments
 (0)