From 947daeb85321ed804bc3142623844b2617cb1b3e Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Wed, 30 Apr 2025 11:23:38 -0300 Subject: Support long double in the integrator --- .gitignore | 1 + Makefile | 11 +++++++++- log-fourier.cpp | 55 +++++++++++++++++++++++----------------------- log-fourier.hpp | 4 ++-- log-fourier_integrator.cpp | 10 ++++----- types.hpp | 39 +++++++++++++++++++++++++++++++- 6 files changed, 84 insertions(+), 36 deletions(-) diff --git a/.gitignore b/.gitignore index fe5f9c3..35be3ad 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,7 @@ correlations fourier_integrator log-fourier_integrator +log-fourier_integrator_long vgcore.* walk get_energy diff --git a/Makefile b/Makefile index 6cf8a52..c5cac94 100644 --- a/Makefile +++ b/Makefile @@ -1,4 +1,4 @@ -all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator log_get_energy +all: walk correlations integrator fourier_integrator get_energy log-fourier_integrator log-fourier_integrator_long log_get_energy CC := clang++ -std=c++17 -flto -Wno-mathematical-notation-identifier-extension -O3 -march=native -mtune=native -fopenmp -DFFTW_THREADS=1 @@ -20,12 +20,21 @@ p-spin.o: p-spin.cpp p-spin.hpp types.hpp log-fourier.o: log-fourier.cpp log-fourier.hpp types.hpp $(CC) log-fourier.cpp -c -o log-fourier.o +p-spin_long.o: p-spin.cpp p-spin.hpp types.hpp + $(CC) p-spin.cpp -D QUAD_PRECISION -c -o p-spin_long.o + +log-fourier_long.o: log-fourier.cpp log-fourier.hpp types.hpp + $(CC) log-fourier.cpp -D QUAD_PRECISION -c -o log-fourier_long.o + fourier_integrator: fourier_integrator.cpp fourier.hpp fourier.o p-spin.o $(CC) fourier_integrator.cpp fourier.o p-spin.o -lfftw3 -lgsl -o fourier_integrator log-fourier_integrator: log-fourier_integrator.cpp log-fourier.hpp log-fourier.o p-spin.o $(CC) log-fourier_integrator.cpp log-fourier.o p-spin.o -lfftw3 -lgsl -o log-fourier_integrator +log-fourier_integrator_long: log-fourier_integrator.cpp log-fourier.hpp log-fourier_long.o p-spin_long.o + $(CC) log-fourier_integrator.cpp log-fourier_long.o p-spin_long.o -D QUAD_PRECISION -lfftw3l -lgsl -o log-fourier_integrator_long + get_energy: get_energy.cpp fourier.hpp fourier.o p-spin.o $(CC) get_energy.cpp fourier.o p-spin.o -lfftw3 -lgsl -o get_energy diff --git a/log-fourier.cpp b/log-fourier.cpp index 47d5c5c..9ae70be 100644 --- a/log-fourier.cpp +++ b/log-fourier.cpp @@ -2,25 +2,26 @@ #include "p-spin.hpp" #include #include +#include LogarithmicFourierTransform::LogarithmicFourierTransform(unsigned N, Real k, Real Δτ, unsigned pad) : N(N), pad(pad), k(k), Δτ(Δτ) { τₛ = -0.5 * N; ωₛ = -0.5 * N; sₛ = -0.5 * pad * N; - a = reinterpret_cast(fftw_alloc_complex(pad*N)); - â = reinterpret_cast(fftw_alloc_complex(pad*N)); - fftw_import_wisdom_from_filename("fftw.wisdom"); - a_to_â = fftw_plan_dft_1d(pad*N, reinterpret_cast(a), reinterpret_cast(â), FFTW_BACKWARD, 0); - â_to_a = fftw_plan_dft_1d(pad*N, reinterpret_cast(â), reinterpret_cast(a), FFTW_BACKWARD, 0); - fftw_export_wisdom_to_filename("fftw.wisdom"); + a = reinterpret_cast(FFTW_ALLOC_COMPLEX(pad*N)); + â = reinterpret_cast(FFTW_ALLOC_COMPLEX(pad*N)); + FFTW_IMPORT_WISDOM("fftw.wisdom"); + a_to_â = FFTW_PLAN_DFT_1D(pad*N, reinterpret_cast(a), reinterpret_cast(â), FFTW_BACKWARD, 0); + â_to_a = FFTW_PLAN_DFT_1D(pad*N, reinterpret_cast(â), reinterpret_cast(a), FFTW_BACKWARD, 0); + FFTW_EXPORT_WISDOM("fftw.wisdom"); } LogarithmicFourierTransform::~LogarithmicFourierTransform() { - fftw_destroy_plan(a_to_â); - fftw_destroy_plan(â_to_a); - fftw_free(a); - fftw_free(â); - fftw_cleanup(); + FFTW_DESTROY_PLAN(a_to_â); + FFTW_DESTROY_PLAN(â_to_a); + FFTW_FREE(a); + FFTW_FREE(â); + FFTW_CLEANUP(); } Real LogarithmicFourierTransform::τ(unsigned n) const { @@ -47,9 +48,9 @@ Complex Γ(Complex z) { gsl_sf_result logΓ; gsl_sf_result argΓ; - gsl_sf_lngamma_complex_e(z.real(), z.imag(), &logΓ, &argΓ); + gsl_sf_lngamma_complex_e((double)z.real(), (double)z.imag(), &logΓ, &argΓ); - return exp(logΓ.val + 1i * argΓ.val); + return exp((Real)logΓ.val + II * (Real)argΓ.val); } std::vector LogarithmicFourierTransform::fourier(const std::vector& c, bool symmetric) { @@ -67,13 +68,13 @@ std::vector LogarithmicFourierTransform::fourier(const std::vector LogarithmicFourierTransform::inverse(const std::vector& C, const std::vector& R, const std::vector& Ct, const std::vector& Rt, unsigned p, unsigned s, Real λ, Real τ₀, Real β, unsigned log2n, Real Δτ, Real k) { - unsigned N = pow(2, log2n); + unsigned N = std::pow(2, log2n); std::ofstream outfile(logFourierFile("C", p, s, λ, τ₀, β, log2n, Δτ, k), std::ios::out | std::ios::binary); outfile.write((const char*)(C.data()), N * sizeof(Real)); outfile.close(); @@ -141,7 +142,7 @@ bool logFourierLoad(std::vector& C, std::vector& R, std::vector& C, con auto [RddfCt, dfCt] = RddfCtdfCt(fft, C, R, p, s, λ); Real Γ₀ = 1.0; - return ((2 * Γ₀ * std::conj(Rt[0]) + pow(β, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real(); + return ((2 * Γ₀ * std::conj(Rt[0]) + std::pow(β, 2) * (RddfCt[0] * Ct[0] + dfCt[0] * std::conj(Rt[0]))) / Ct[0]).real(); } Real energy(const LogarithmicFourierTransform& fft, std::vector& C, const std::vector& R, unsigned p, unsigned s, Real λ, Real β) { @@ -188,7 +189,7 @@ Real energy(const LogarithmicFourierTransform& fft, std::vector& C, const Real f₂ₙ₊₂ = R[2*n+2] * df(λ, p, s, C[2*n+2]); E += (h₂ₙ + h₂ₙ₊₁) / 6 * ( (2 - h₂ₙ₊₁ / h₂ₙ) * f₂ₙ - + pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁ + + std::pow(h₂ₙ + h₂ₙ₊₁, 2) / (h₂ₙ * h₂ₙ₊₁) * f₂ₙ₊₁ + (2 - h₂ₙ / h₂ₙ₊₁) * f₂ₙ₊₂ ); } diff --git a/log-fourier.hpp b/log-fourier.hpp index e232bd6..97c0e15 100644 --- a/log-fourier.hpp +++ b/log-fourier.hpp @@ -11,8 +11,8 @@ class LogarithmicFourierTransform { private: Complex* a; Complex* â; - fftw_plan a_to_â; - fftw_plan â_to_a; + FFTW_PLAN a_to_â; + FFTW_PLAN â_to_a; unsigned N; unsigned pad; Real k; diff --git a/log-fourier_integrator.cpp b/log-fourier_integrator.cpp index 9bbf554..7ac76b5 100644 --- a/log-fourier_integrator.cpp +++ b/log-fourier_integrator.cpp @@ -79,7 +79,7 @@ int main(int argc, char* argv[]) { LogarithmicFourierTransform fft(N, k, Δτ, 4); - Real Γ₀ = 1.0; + Real Γ₀ = 1; Real μₜ₋₁ = Γ₀; if (τ₀ > 0) { μₜ₋₁ = (sqrt(1+4*Γ₀*τ₀)-1)/(2*τ₀); @@ -105,7 +105,7 @@ int main(int argc, char* argv[]) { Rₜ₋₁[n] = exp(-μₜ₋₁ * fft.t(n)); Ĉₜ₋₁[n] = 2 * Γ₀ / (pow(μₜ₋₁, 2) + pow(fft.ν(n), 2)) / (1 + pow(τ₀ * fft.ν(n), 2)); - Ȓₜ₋₁[n] = 1.0 / (μₜ₋₁ + 1i * fft.ν(n)); + Ȓₜ₋₁[n] = (Real)1.0 / (μₜ₋₁ + II * fft.ν(n)); } } else { logFourierLoad(Cₜ₋₁, Rₜ₋₁, Ĉₜ₋₁, Ȓₜ₋₁, p, s, λ, τ₀, β₀, log2n, Δτ, k); @@ -130,8 +130,8 @@ int main(int argc, char* argv[]) { std::vector Ĉₜ₊₁(N); std::vector Ȓₜ₊₁(N); for (unsigned n = 0; n < N; n++) { - Ȓₜ₊₁[n] = (1.0 + pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μₜ + 1i * fft.ν(n)); - Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + pow(τ₀ * fft.ν(n), 2)) + pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + 1i * fft.ν(n)); + Ȓₜ₊₁[n] = ((Real)1.0 + std::pow(β, 2) * RddfCt[n] * Ȓₜ[n]) / (μₜ + II * fft.ν(n)); + Ĉₜ₊₁[n] = (2 * Γ₀ * std::conj(Ȓₜ[n]) / (1 + std::pow(τ₀ * fft.ν(n), 2)) + std::pow(β, 2) * (RddfCt[n] * Ĉₜ[n] + dfCt[n] * std::conj(Ȓₜ[n]))) / (μₜ + II * fft.ν(n)); } std::vector Rₜ₊₁ = fft.inverse(Ȓₜ₊₁); std::vector Cₜ₊₁ = fft.inverse(Ĉₜ₊₁); @@ -153,7 +153,7 @@ int main(int argc, char* argv[]) { } if (stepsUp > stepsToRespond) { - γ = std::max(γ/2, 1e-4); + γ = std::max(γ/2, (Real)1e-4); stepsUp = 0; ΔCmin = ΔCₜ; } diff --git a/types.hpp b/types.hpp index 05eceb4..71556fd 100644 --- a/types.hpp +++ b/types.hpp @@ -1,6 +1,43 @@ #pragma once #include +#include + +using namespace std::complex_literals; + +#ifndef QUAD_PRECISION + +#define FFTW_PLAN fftw_plan +#define FFTW_ALLOC_COMPLEX fftw_alloc_complex +#define FFTW_IMPORT_WISDOM fftw_import_wisdom_from_filename +#define FFTW_EXPORT_WISDOM fftw_export_wisdom_to_filename +#define FFTW_PLAN_DFT_1D fftw_plan_dft_1d +#define FFTW_DESTROY_PLAN fftw_destroy_plan +#define FFTW_FREE fftw_free +#define FFTW_CLEANUP fftw_cleanup +#define FFTW_EXECUTE fftw_execute +#define FFTW_COMPLEX fftw_complex using Real = double; + +#define II 1i + +#else + +#define FFTW_PLAN fftwl_plan +#define FFTW_ALLOC_COMPLEX fftwl_alloc_complex +#define FFTW_IMPORT_WISDOM fftwl_import_wisdom_from_filename +#define FFTW_EXPORT_WISDOM fftwl_export_wisdom_to_filename +#define FFTW_PLAN_DFT_1D fftwl_plan_dft_1d +#define FFTW_DESTROY_PLAN fftwl_destroy_plan +#define FFTW_FREE fftwl_free +#define FFTW_CLEANUP fftwl_cleanup +#define FFTW_EXECUTE fftwl_execute +#define FFTW_COMPLEX fftwl_complex + +using Real = long double; + +#define II 1il + +#endif + using Complex = std::complex; -using namespace std::complex_literals; -- cgit v1.2.3-70-g09d2