Skip to content

Remove unnecessary statistic parameter from convolution-related functions #17

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Remove unnecessary statistic parameter from convolution-related functions #17

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions c++/cppdlr/dlr_dyson.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ namespace cppdlr {
auto g0c = itops_ptr->vals2coefs(g0); // DLR coefficients of free Green's function

// Get matrix of convolution by free Green's function
g0mat = itops_ptr->convmat(beta, Fermion, g0c, time_order);
g0mat = itops_ptr->convmat(beta, g0c, time_order);

// Get right hand side of Dyson equation
if constexpr (std::floating_point<Ht>) { // If h is real scalar, rhs is a vector
Expand All @@ -86,7 +86,7 @@ namespace cppdlr {
* \note Hamiltonian must either be a symmetric matrix, a Hermitian matrix,
* or a real scalar.
*/
dyson_it(double beta, imtime_ops itops, Ht const &h, bool time_order) : dyson_it(beta, itops, h, 0, time_order) {};
dyson_it(double beta, imtime_ops itops, Ht const &h, bool time_order) : dyson_it(beta, itops, h, 0, time_order){};

/**
* @brief Solve Dyson equation for given self-energy
Expand All @@ -111,7 +111,7 @@ namespace cppdlr {
// Obtain Dyson equation system matrix I - G0 * Sig, where G0 and Sig are the
// matrices of convolution by the free Green's function and self-energy,
// respectively.
auto sysmat = make_regular(nda::eye<double>(r * norb) - g0mat * itops_ptr->convmat(beta, Fermion, sigc, time_order));
auto sysmat = make_regular(nda::eye<double>(r * norb) - g0mat * itops_ptr->convmat(beta, sigc, time_order));

// Factorize system matrix
auto ipiv = nda::vector<int>(r * norb);
Expand Down
62 changes: 52 additions & 10 deletions c++/cppdlr/dlr_imtime.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -304,13 +304,10 @@ namespace cppdlr {
* @return Values of h = f * g on DLR imaginary time grid
* */
template <nda::MemoryArray T, nda::Scalar S = nda::get_value_t<T>>
typename T::regular_type convolve(double beta, statistic_t statistic, T const &fc, T const &gc, bool time_order = false) const {
typename T::regular_type convolve(double beta, T const &fc, T const &gc, bool time_order = false) const {

if (r != fc.shape(0) || r != gc.shape(0)) throw std::runtime_error("First dim of input arrays must be equal to DLR rank r.");

// TODO: implement bosonic case and remove
if (statistic == 0) throw std::runtime_error("imtime_ops::convolve not yet implemented for bosonic Green's functions.");

// Initialize convolution, if it hasn't been done already
if (!time_order & hilb.empty()) { convolve_init(); }
if (time_order & thilb.empty()) { tconvolve_init(); }
Expand Down Expand Up @@ -353,6 +350,22 @@ namespace cppdlr {
}
}

/**
* @brief Deprecated overload of convolve method
*
* This version includes an unused extra parameter for backward
* compatibility.
*
* @deprecated Use convolve(beta, fc, gc, time_order) instead; this works for
* both fermionic and bosonic functions.
*/
template <nda::MemoryArray T, nda::Scalar S = nda::get_value_t<T>>
[[deprecated("Use convolve(beta, fc, gc, time_order) instead.")]] typename T::regular_type
convolve(double beta, statistic_t statistic, T const &fc, T const &gc, bool time_order = false) const {
(void)statistic; // Unused parameter, kept for backward compatibility
return convolve(beta, fc, gc, time_order);
}

/**
* @brief Compute convolution of two imaginary time Green's functions,
* given matrix of convolution by one of them
Expand Down Expand Up @@ -449,7 +462,7 @@ namespace cppdlr {
* r*norb2 x r*norb3 matrix, or a block r x 1 matrix of norb2 x norb3 blocks.
* */
template <nda::MemoryArray T, nda::Scalar S = nda::get_value_t<T>>
nda::matrix<S> convmat(double beta, statistic_t statistic, T const &fc, bool time_order = false) const {
nda::matrix<S> convmat(double beta, T const &fc, bool time_order = false) const {

int n, m;

Expand All @@ -464,11 +477,27 @@ namespace cppdlr {
}

auto fconv = nda::matrix<S, nda::C_layout>(n, m); // Matrix of convolution by f
convmat_inplace(nda::matrix_view<S, nda::C_layout>(fconv), beta, statistic, fc, time_order);
convmat_inplace(nda::matrix_view<S, nda::C_layout>(fconv), beta, fc, time_order);

return fconv;
}

/**
* @brief Deprecated overload of convmat method
*
* This version includes an unused extra parameter for backward
* compatibility.
*
* @deprecated Use convmat(beta, fc, time_order) instead; this works for both
* fermionic and bosonic functions.
*/
template <nda::MemoryArray T, nda::Scalar S = nda::get_value_t<T>>
[[deprecated("Use convmat(beta, fc, time_order) instead.")]] nda::matrix<S> convmat(double beta, statistic_t statistic, T const &fc,
bool time_order = false) const {
(void)statistic; // Unused parameter, kept for backward compatibility
return convmat(beta, fc, time_order);
}

/**
* @brief Compute matrix of convolution by an imaginary time Green's function
* in place
Expand Down Expand Up @@ -519,13 +548,10 @@ namespace cppdlr {
* r*norb2 x r*norb3 matrix, or a block r x 1 matrix of norb2 x norb3 blocks.
* */
template <nda::MemoryArray T, nda::Scalar S = nda::get_value_t<T>>
void convmat_inplace(nda::matrix_view<S, nda::C_layout> fconv, double beta, statistic_t statistic, T const &fc, bool time_order = false) const {
void convmat_inplace(nda::matrix_view<S, nda::C_layout> fconv, double beta, T const &fc, bool time_order = false) const {

if (r != fc.shape(0)) throw std::runtime_error("First dim of input array must be equal to DLR rank r.");

// TODO: implement bosonic case and remove
if (statistic == 0) throw std::runtime_error("imtime_ops::convmat not yet implemented for bosonic Green's functions.");

// Initialize convolution, if it hasn't been done already
if (!time_order & hilb.empty()) { convolve_init(); }
if (time_order & thilb.empty()) { tconvolve_init(); }
Expand Down Expand Up @@ -622,6 +648,22 @@ namespace cppdlr {
}
}

/**
* @brief Deprecated overload of convmat_inplace method
*
* This version includes an unused extra parameter for backward
* compatibility.
*
* @deprecated Use convmat_inplace(fconv, beta, fc, time_order) instead; this
* works for both fermionic and bosonic functions.
*/
template <nda::MemoryArray T, nda::Scalar S = nda::get_value_t<T>>
[[deprecated("Use convmat_inplace(fconv, beta, fc, time_order) instead.")]] void
convmat_inplace(nda::matrix_view<S, nda::C_layout> fconv, double beta, statistic_t statistic, T const &fc, bool time_order = false) const {
(void)statistic; // Unused parameter, kept for backward compatibility
return convmat_inplace(fconv, beta, fc, time_order);
}

/**
* @brief Compute inner product of two imaginary time Green's functions
*
Expand Down
32 changes: 16 additions & 16 deletions test/c++/imtime_ops.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -492,13 +492,13 @@ TEST(imtime_ops, convolve_scalar_real) {
auto gc = itops.vals2coefs(g);

// Get convolution and time-ordered convolution of f and g directly
auto h = itops.convolve(beta, Fermion, fc, gc);
auto ht = itops.convolve(beta, Fermion, fc, gc, TIME_ORDERED);
auto h = itops.convolve(beta, fc, gc);
auto ht = itops.convolve(beta, fc, gc, TIME_ORDERED);

// Get convolution and time-ordered convolution of f and g by first forming
// matrix of convolution by f and then applying it to g
auto h2 = itops.convolve(itops.convmat(beta, Fermion, fc), g);
auto ht2 = itops.convolve(itops.convmat(beta, Fermion, fc, TIME_ORDERED), g);
auto h2 = itops.convolve(itops.convmat(beta, fc), g);
auto ht2 = itops.convolve(itops.convmat(beta, fc, TIME_ORDERED), g);

// Check that the two methods give the same result
EXPECT_LT(max_element(abs(h - h2)), 1e-14);
Expand Down Expand Up @@ -577,13 +577,13 @@ TEST(imtime_ops, convolve_scalar_cmplx) {
auto gc = itops.vals2coefs(g);

// Get convolution and time-ordered convolution of f and g directly
auto h = itops.convolve(beta, Fermion, fc, gc);
auto ht = itops.convolve(beta, Fermion, fc, gc, TIME_ORDERED);
auto h = itops.convolve(beta, fc, gc);
auto ht = itops.convolve(beta, fc, gc, TIME_ORDERED);

// Get convolution and time-ordered convolution of f and g by first forming
// matrix of convolution by f and then applying it to g
auto h2 = itops.convolve(itops.convmat(beta, Fermion, fc), g);
auto ht2 = itops.convolve(itops.convmat(beta, Fermion, fc, TIME_ORDERED), g);
auto h2 = itops.convolve(itops.convmat(beta, fc), g);
auto ht2 = itops.convolve(itops.convmat(beta, fc, TIME_ORDERED), g);

// Check that the two methods give the same result
EXPECT_LT(max_element(abs(h - h2)), 1e-14);
Expand Down Expand Up @@ -667,11 +667,11 @@ TEST(imtime_ops, convolve_matrix_real) {
auto gc = itops.vals2coefs(g);

// Get convolution and time-ordered convolution of f and g directly
auto h = itops.convolve(beta, Fermion, fc, gc);
auto ht = itops.convolve(beta, Fermion, fc, gc, TIME_ORDERED);
auto h = itops.convolve(beta, fc, gc);
auto ht = itops.convolve(beta, fc, gc, TIME_ORDERED);

auto h2 = itops.convolve(itops.convmat(beta, Fermion, fc), g);
auto ht2 = itops.convolve(itops.convmat(beta, Fermion, fc, TIME_ORDERED), g);
auto h2 = itops.convolve(itops.convmat(beta, fc), g);
auto ht2 = itops.convolve(itops.convmat(beta, fc, TIME_ORDERED), g);

// Check that the two methods give the same result
EXPECT_LT(max_element(abs(h - h2)), 1e-14);
Expand Down Expand Up @@ -758,13 +758,13 @@ TEST(imtime_ops, convolve_matrix_cmplx) {
auto gc = itops.vals2coefs(g);

// Get convolution and time-ordered convolution of f and g directly
auto h = itops.convolve(beta, Fermion, fc, gc);
auto ht = itops.convolve(beta, Fermion, fc, gc, TIME_ORDERED);
auto h = itops.convolve(beta, fc, gc);
auto ht = itops.convolve(beta, fc, gc, TIME_ORDERED);

// Get convolution and time-ordered convolution of f and g by first forming
// matrix of convolution by f and then applying it to g
auto h2 = itops.convolve(itops.convmat(beta, Fermion, fc), g);
auto ht2 = itops.convolve(itops.convmat(beta, Fermion, fc, TIME_ORDERED), g);
auto h2 = itops.convolve(itops.convmat(beta, fc), g);
auto ht2 = itops.convolve(itops.convmat(beta, fc, TIME_ORDERED), g);

// Check that the two methods give the same result
EXPECT_LT(max_element(abs(h - h2)), 1e-14);
Expand Down