change semantics of imaginary check for function imaginary->real (c13a5107) · Commits · NDIP / Tool Sources / Direct-Geometry Spectroscopy / DCA / DCA Main

include/dca/function/util/real_complex_conversion.hpp

+40 −10

Original line number	Diff line number	Diff line
		@@ -19,6 +19,7 @@
		#include "dca/function/function.hpp"

		namespace dca {
		enum class ImagCheck { IGNORE, WARN, FAIL };
		namespace func {
		namespace util {
		// dca::func::util::
		@@ -34,23 +35,52 @@ auto complex(const function<Scalartype, Dmn>& f) {
		return f_complex;
		}

		// Returns a real valued function that is equal to the real part of f.
		// If check_imaginary = true, the method checks whether the imaginary part of f is zero and throws a
		// std::logic_error if this is not the case.
		template <typename Scalartype, typename Dmn>
		auto real(const function<std::complex<Scalartype>, Dmn>& f, const bool check_imaginary = false) {
		namespace detail {
		template <typename Scalartype, typename Dmn, ImagCheck IC>
		auto real(const function<std::complex<Scalartype>, Dmn>& f) {
		function<Scalartype, Dmn> f_real;

		auto checkForImaginary = [](const std::complex<Scalartype> s) -> bool {
		return (std::abs(s.imag()) > (1000 * std::numeric_limits<Scalartype>::epsilon()));
		};
		auto writeCheckFail = [](const int i, const std::complex<Scalartype> s) -> std::string {
		std::ostringstream fail_message;
		fail_message << "Element " << i << " of function conversion complex -> real "
		<< " has a non negligible imaginary component of " << std::abs(s.imag()) << '\n';
		return fail_message.str();
		};
		bool have_not_warned = true;
		for (int i = 0; i < f_real.size(); ++i) {
		if (check_imaginary && std::abs(f(i).imag()) > 1000 * std::numeric_limits<Scalartype>::epsilon()) {
		std::ostringstream err_msg;
		throw(std::logic_error(err_msg.str()));
		if constexpr (IC == ImagCheck::WARN) {
		if (have_not_warned) {
		if (checkForImaginary(f(i)))
		std::cerr << writeCheckFail(i, f(i));
		}
		}
		else if constexpr (IC == ImagCheck::FAIL) {
		if (checkForImaginary(f(i)))
		throw(std::out_of_range(writeCheckFail(i, f(i))));
		}
		f_real(i) = f(i).real();
		}

		return f_real;
		}
		} // namespace detail
		// Returns a real valued function that is equal to the real part of f.
		// If IC ==
		// ImagCheck::FAIL -> throw exception if imaginary part of f is zero
		// ImagCheck::WARN -> write one warning per conversion to std::err continue without further
		// checks ImagCheck::IGNORE -> do not check for imaginary elements in f at all
		template <typename Scalartype, typename Dmn>
		auto real(const function<std::complex<Scalartype>, Dmn>& f, ImagCheck ic = ImagCheck::IGNORE) {
		switch (ic) {
		case ImagCheck::IGNORE:
		return detail::real<Scalartype, Dmn, ImagCheck::IGNORE>(f);
		case ImagCheck::WARN:
		return detail::real<Scalartype, Dmn, ImagCheck::WARN>(f);
		case ImagCheck::FAIL:
		return detail::real<Scalartype, Dmn, ImagCheck::FAIL>(f);
		}
		}

		} // namespace util
		} // namespace func

include/dca/math/function_transform/special_transforms/momentum_to_space.hpp

+1 −1

Original line number	Diff line number	Diff line
		@@ -68,7 +68,7 @@ public:
		func::function<Real, func::dmn_variadic<NuDmn, NuDmn, RDmn, LastDmn>>& f_output) {
		func::function<Complex, func::dmn_variadic<NuDmn, NuDmn, RDmn, LastDmn>> f_out_cmplx;
		execute(f_input, f_out_cmplx);
		f_output = std::move(func::util::real(f_out_cmplx, true));
		f_output = std::move(func::util::real(f_out_cmplx, ImagCheck::WARN));
		}

		// Default to old implementation when no band is present.

include/dca/math/function_transform/special_transforms/space_to_momentum.hpp

+1 −1

Original line number	Diff line number	Diff line
		@@ -72,7 +72,7 @@ public:
		func::function<Real, func::dmn_variadic<NuDmn, NuDmn, KDmn, LastDmn>>& f_output) {
		func::function<Complex, func::dmn_variadic<NuDmn, NuDmn, KDmn, LastDmn>> f_out_cmplx;
		execute(f_input, f_out_cmplx);
		f_output = std::move(func::util::real(f_out_cmplx, false));
		f_output = std::move(func::util::real(f_out_cmplx, ImagCheck::WARN));
		}

		// Default to old implementation when no band is present.

include/dca/phys/dca_algorithms/compute_free_greens_function.hpp

+1 −1

Original line number	Diff line number	Diff line
		@@ -90,7 +90,7 @@ void compute_G0_k_t(

		compute_G0_k_t(H0_k, mu, beta, f_cmplx);

		G0_k_t = std::move(func::util::real(f_cmplx, true));
		G0_k_t = std::move(func::util::real(f_cmplx, ImagCheck::FAIL));
		}

		} // namespace phys

test/unit/function/util/real_complex_conversion_test.cpp

+4 −4

Original line number	Diff line number	Diff line
		@@ -37,16 +37,16 @@ TEST(RealComplexConversionTest, ComplexToReal) {
		f(0) = std::complex<double>(0.5, 0.);
		f(1) = std::complex<double>(1., 1.);

		auto f_real = dca::func::util::real(f, false);
		auto f_real = dca::func::util::real(f);

		EXPECT_DOUBLE_EQ(0.5, f_real(0));
		EXPECT_DOUBLE_EQ(1., f_real(1));

		// If the second argument is true, the function's imaginary part must be zero.
		EXPECT_THROW(dca::func::util::real(f, true), std::logic_error);
		// Template argument controls behavior of check, the function's imaginary part must be zero.
		EXPECT_THROW(dca::func::util::real(f, dca::ImagCheck::FAIL), std::out_of_range);

		f(1).imag(0.);
		f_real = dca::func::util::real(f, true);
		f_real = dca::func::util::real(f, dca::ImagCheck::FAIL);

		EXPECT_DOUBLE_EQ(0.5, f_real(0));
		EXPECT_DOUBLE_EQ(1., f_real(1));