From 94a8441949a5a9ceac7c9de1c037ca5591b7a8fb Mon Sep 17 00:00:00 2001 From: "Hugo U. R. Strand" Date: Fri, 13 Dec 2024 11:09:07 +0100 Subject: [PATCH 1/3] [itops] split convmat in allocation and a convmat_inplace builder function. --- c++/cppdlr/dlr_imtime.hpp | 81 ++++++++++++++++++++++++++++++++++----- 1 file changed, 72 insertions(+), 9 deletions(-) diff --git a/c++/cppdlr/dlr_imtime.hpp b/c++/cppdlr/dlr_imtime.hpp index abf3543..19e63c8 100644 --- a/c++/cppdlr/dlr_imtime.hpp +++ b/c++/cppdlr/dlr_imtime.hpp @@ -451,6 +451,70 @@ namespace cppdlr { template > nda::matrix convmat(double beta, statistic_t statistic, T const &fc, bool time_order = false) const { + int n, m; + + if constexpr (T::rank == 1) { // Scalar-valued Green's function + n = r; + m = r; + } else if (T::rank == 3) { // Matrix-valued Green's function + n = r * fc.shape(1); + m = r * fc.shape(2); + } else { + throw std::runtime_error("Input arrays must be rank 1 (scalar-valued Green's function) or 3 (matrix-valued Green's function)."); + } + + auto fconv = nda::matrix(n, m); // Matrix of convolution by f + convmat_inplace(nda::matrix_view(fconv), beta, statistic, fc, time_order); + + return fconv; + } + + /** + * @brief Compute matrix of convolution by an imaginary time Green's function + * + * The convolution of f and g is defined as h(t) = (f * g)(t) = int_0^beta + * f(t-t') g(t') dt', where fermionic/bosonic antiperiodicity/periodicity are + * used to define the Green's functions on (-beta, 0). This method takes the + * DLR coefficients of f as input and returns the matrix of convolution by f. + * This matrix can be applied to the values of g on the DLR imaginary time + * grid, to produce the values of h on the DLR imaginary time grid. + * + * By specifying the @p time_order flag, this method can be used to compute + * the time-ordered convolution of f and g, defined as h(t) = (f * g)(t) = + * int_0^tau f(t-t') g(t') dt'. + * + * The convolution matrix is constructed using the method described in + * Appendix A of + * + * J. Kaye, H. U. R. Strand, D. Golez, "Decomposing imaginary time Feynman + * diagrams using separable basis functions: Anderson impurity model strong + * coupling expansion," arXiv:2307.08566 (2023). + * + * @param[out] fconv Convolution matrix from DLR coefficients to DLR grid + * @param[in] beta Inverse temperature + * @param[in] statistic Fermionic ("Fermion" or 0) or bosonic ("Boson" or 1) + * @param[in] fc DLR coefficients of f + * @param[in] time_order Flag for ordinary (false or ORDINARY, default) or + * time-ordered (true or TIME_ORDERED) convolution + * + * \note Whereas the method imtime_ops::convolve takes the DLR coefficients + * of f and g as input and computes their convolution h directly, this method + * returns a matrix which should be applied to the DLR imaginary time grid values + * of g, rather than its DLR coefficients, in to order to obtain the + * convolution h. The purpose of this is to make the input and output + * representations of the convolution matrix equal, which is often convenient + * in practice. + * + * \note In the case of matrix-valued Green's functions, we think of the + * matrix of convolution by f as an r*norb x r*norb matrix, or a block r x r + * matrix of norb x norb blocks. Here r is the DLR rank and norb is the + * number of orbital indices. This matrix would then be applied to a Green's + * function g, represented as an r*norb x norb matrix, or a block r x 1 + * matrix of norb x norb blocks. + * */ + template > + void convmat_inplace(nda::matrix_view fconv, double beta, statistic_t statistic, 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 @@ -466,10 +530,9 @@ namespace cppdlr { if constexpr (T::rank == 1) { // Scalar-valued Green's function - // First construct convolution matrix from DLR coefficients to DLR grid - // values - auto fconv = nda::matrix(r, r); // Matrix of convolution by f - + if (fconv.shape(0) != r || fconv.shape(1) != r) + throw std::runtime_error("Matrix shape must be equal to DLR rank (r,r)."); + // Diagonal contribution (given by diag(tau_k) * K(tau_k, om_l) * diag(fc_l)) for (int k = 0; k < r; ++k) { for (int l = 0; l < r; ++l) { fconv(k, l) = tcf2it_v(k, l) * fc(l); } @@ -495,16 +558,16 @@ namespace cppdlr { nda::lapack::getrs(transpose(it2cf.lu), fconv, it2cf.piv); } - return beta * fconv; + fconv *= beta; } else if (T::rank == 3) { // Matrix-valued Green's function int norb1 = fc.shape(1); int norb2 = fc.shape(2); - // First construct convolution matrix from DLR coefficients to DLR grid - // values - auto fconv = nda::matrix(r * norb1, r * norb2); // Matrix of convolution by f + if (fconv.shape(0) != r*norb1 || fconv.shape(1) != r*norb2) + throw std::runtime_error("Matrix shape must be equal to DLR rank times norbs (r*norb1,r*norb2)."); + auto fconv_rs = nda::reshape(fconv, r, norb1, r, norb2); // Array view to index into fconv for conevenience // Diagonal contribution (given by diag(tau_k) * K(tau_k, om_l) * diag(fc_l)) @@ -547,7 +610,7 @@ namespace cppdlr { for (int k = 0; k < r; ++k) { fconv_rs(_, _, k, i) = fconvtmp_rs(_, _, i, k); } } - return beta * fconv; + fconv *= beta; } else { throw std::runtime_error("Input arrays must be rank 1 (scalar-valued Green's function) or 3 (matrix-valued Green's function)."); From 939e7b0ea8e2be764234447081d067d12391cf40 Mon Sep 17 00:00:00 2001 From: "Hugo U. R. Strand" Date: Tue, 4 Feb 2025 12:01:18 +0100 Subject: [PATCH 2/3] [itops] convmat_inplace add doc note on rec inplace function usage --- c++/cppdlr/dlr_imtime.hpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/c++/cppdlr/dlr_imtime.hpp b/c++/cppdlr/dlr_imtime.hpp index 19e63c8..0641e69 100644 --- a/c++/cppdlr/dlr_imtime.hpp +++ b/c++/cppdlr/dlr_imtime.hpp @@ -490,13 +490,18 @@ namespace cppdlr { * diagrams using separable basis functions: Anderson impurity model strong * coupling expansion," arXiv:2307.08566 (2023). * - * @param[out] fconv Convolution matrix from DLR coefficients to DLR grid + * @param[in/out] fconv Convolution matrix from DLR coefficients to DLR grid * @param[in] beta Inverse temperature * @param[in] statistic Fermionic ("Fermion" or 0) or bosonic ("Boson" or 1) * @param[in] fc DLR coefficients of f * @param[in] time_order Flag for ordinary (false or ORDINARY, default) or * time-ordered (true or TIME_ORDERED) convolution * + * \note This function builds the matrix of convolution, in place, in the + * provided matrix `fconv`. This makes it possible to control the memory + * allocation externally. If this is not a concern, we advice using the + * `covmat(...)` function instead (of `convmat_inplace(...)`). + * * \note Whereas the method imtime_ops::convolve takes the DLR coefficients * of f and g as input and computes their convolution h directly, this method * returns a matrix which should be applied to the DLR imaginary time grid values From 4f5cce4da055f71d9d111698988ae89ca142e95e Mon Sep 17 00:00:00 2001 From: Jason Kaye Date: Tue, 4 Feb 2025 09:22:16 -0500 Subject: [PATCH 3/3] doc typo fixes --- c++/cppdlr/dlr_imtime.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/c++/cppdlr/dlr_imtime.hpp b/c++/cppdlr/dlr_imtime.hpp index 0641e69..09c2840 100644 --- a/c++/cppdlr/dlr_imtime.hpp +++ b/c++/cppdlr/dlr_imtime.hpp @@ -471,6 +471,7 @@ namespace cppdlr { /** * @brief Compute matrix of convolution by an imaginary time Green's function + * in place * * The convolution of f and g is defined as h(t) = (f * g)(t) = int_0^beta * f(t-t') g(t') dt', where fermionic/bosonic antiperiodicity/periodicity are @@ -499,8 +500,8 @@ namespace cppdlr { * * \note This function builds the matrix of convolution, in place, in the * provided matrix `fconv`. This makes it possible to control the memory - * allocation externally. If this is not a concern, we advice using the - * `covmat(...)` function instead (of `convmat_inplace(...)`). + * allocation externally. If this is not a concern, we advise using the + * `convmat(...)` function instead of `convmat_inplace(...)`. * * \note Whereas the method imtime_ops::convolve takes the DLR coefficients * of f and g as input and computes their convolution h directly, this method