From b944ab51dab08d8c73c6b4ccc83e481b7546afb2 Mon Sep 17 00:00:00 2001 From: Moritz Leucke Date: Tue, 26 Mar 2024 11:58:43 +0100 Subject: [PATCH] changed pade AC interface --- CMakeLists.txt | 2 +- GX-AnalyticContinuation/CMakeLists.txt | 4 +- .../api/{gx_ac.f90 => gx_ac.F90} | 192 +++++++++++------- .../test/test_pade_approximant.f90 | 82 ++++---- 4 files changed, 168 insertions(+), 112 deletions(-) rename GX-AnalyticContinuation/api/{gx_ac.f90 => gx_ac.F90} (62%) diff --git a/CMakeLists.txt b/CMakeLists.txt index b632cbe1..3f456191 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -68,7 +68,7 @@ endif () # Optional multiple precision arithmetic option(ENABLE_GNU_GMP "Enable the GNU GMP library for multiple precision arithmetic" ON) if (${ENABLE_GNU_GMP}) - find_package(GMPXX REQUIRED) + find_package(GMPXX) endif() # Optional unit testing lib diff --git a/GX-AnalyticContinuation/CMakeLists.txt b/GX-AnalyticContinuation/CMakeLists.txt index 45ec924c..9faa8837 100644 --- a/GX-AnalyticContinuation/CMakeLists.txt +++ b/GX-AnalyticContinuation/CMakeLists.txt @@ -12,11 +12,13 @@ set_target_properties(LibGXAC LIBRARY_OUTPUT_DIRECTORY ${CMAKE_Fortran_LIB_DIRECTORY}) target_include_directories(LibGXAC PUBLIC src/) -target_sources(LibGXAC PRIVATE src/pade_approximant.f90 src/pade_mp.cpp api/gx_ac.f90) if(GMPXX_FOUND) + add_definitions(-DGMPXX_FOUND) + target_sources(LibGXAC PRIVATE src/pade_approximant.f90 src/pade_mp.cpp api/gx_ac.F90) target_include_directories(LibGXAC PRIVATE ${GMPXX_INCLUDE_DIRS}) target_link_libraries(LibGXAC GXCommon ${GMPXX_LIBRARIES} gmp) else() + target_sources(LibGXAC PRIVATE src/pade_approximant.f90 api/gx_ac.F90) target_link_libraries(LibGXAC GXCommon) endif() diff --git a/GX-AnalyticContinuation/api/gx_ac.f90 b/GX-AnalyticContinuation/api/gx_ac.F90 similarity index 62% rename from GX-AnalyticContinuation/api/gx_ac.f90 rename to GX-AnalyticContinuation/api/gx_ac.F90 index 0f6d955f..e9117512 100644 --- a/GX-AnalyticContinuation/api/gx_ac.f90 +++ b/GX-AnalyticContinuation/api/gx_ac.F90 @@ -1,5 +1,5 @@ ! *************************************************************************************************** -! Copyright (C) 2020-2023 GreenX library +! Copyright (C) 2020-2024 GreenX library ! This file is distributed under the terms of the APACHE2 License. ! ! *************************************************************************************************** @@ -10,29 +10,30 @@ module gx_ac use pade_approximant, only: evaluate_thiele_pade, thiele_pade, c_zero, c_one implicit none + private public :: thiele_pade_api, & - params, & - create_thiele_pade, & - evaluate_thiele_pade_at, & - params_mp, & - create_thiele_pade_mp, & - evaluate_thiele_pade_mp, & - free_pade_model - - !> brief store arbitrary precision parameters - type :: params_mp - logical :: initialized = .False. - type(c_ptr) :: params_ptr - end type params_mp - - type :: params + params, & + create_thiele_pade, & + evaluate_thiele_pade_at, & + free_params + + !> brief store the parameters of the tiehle pade model + !> (potentially in abitrary precision floats using GMP) + type :: params logical :: initialized = .false. integer :: n_par + integer :: precision + logical :: multiprecision_used_internally + + ! for multiple precision arithmetic + type(c_ptr) :: params_ptr + + ! for double precision arithmetic complex(kind=dp), dimension(:), allocatable :: x_ref complex(kind=dp), dimension(:), allocatable :: a_par - end type params - + end type params +#ifdef GMPXX_FOUND interface !> auxiliary function to compute Thiele-Pade parameters using arbitrary precision numbers @@ -69,6 +70,7 @@ subroutine free_pade_model(params_ptr) bind(C, name="free_pade_model") end subroutine free_pade_model end interface +#endif contains @@ -104,25 +106,94 @@ subroutine thiele_pade_api(n_par, x_ref, y_ref, x_query, y_query, do_greedy) end subroutine thiele_pade_api - type(params) function create_thiele_pade(n_par, x, y, do_greedy) result(par) + + + !> API function to compute Thiele-Pade parameters + !! (potentially using arbitrary precision arithmetic) + !! + !! @param[in] n_par - order of the interpolant + !! @param[in] x - array of the reference points + !! @param[in] y - array of the reference function values + !! @param[in] do_greedy - whether to use the default greedy algorithm or the naive one + !! @param[in] precision - precision in bits (!! not bytes !!) + !! @return params - abstract type to store all parameters in arb. prec. representation + type(params) function create_thiele_pade(n_par, x, y, do_greedy, precision) result(par) integer, intent(in) :: n_par complex(kind=dp), dimension(:), intent(in) :: x, y logical, optional, intent(in) :: do_greedy + integer, optional, intent(in) :: precision + + ! Internal variables + integer :: local_do_greedy = 1 ! initialize type par%initialized = .true. par%n_par = n_par - allocate(par%a_par(n_par)) - allocate(par%x_ref(size(x))) - ! compute the coefficients - par%x_ref(:) = x - call thiele_pade(n_par, par%x_ref, y, par%a_par, do_greedy) + ! precision of arithmetic internally + if (present(precision)) then + if (precision .eq. 64) then + ! double precision case + par%precision = precision + par%multiprecision_used_internally = .false. + else + ! arbitrary precision case + par%precision = precision + par%multiprecision_used_internally = .true. +#ifdef GMPXX_FOUND +#else + print *, "*** ERROR: multiple precision float arithmetic requested but not linked against GMP library" + return +#endif + + end if + else +#ifdef GMPXX_FOUND + ! default is quadrupel precision if GMP is linked + par%precision = 128 + par%multiprecision_used_internally = .true. +#else + ! default is double precision if GMP is not linked + par%precision = 64 + par%multiprecision_used_internally = .false. +#endif + end if + + if (.not. par%multiprecision_used_internally) then + + allocate(par%a_par(n_par)) + allocate(par%x_ref(size(x))) + + ! compute the coefficients + par%x_ref(:) = x + call thiele_pade(n_par, par%x_ref, y, par%a_par, do_greedy) + + elseif (par%multiprecision_used_internally) then + + ! use integer bools for interoperability with C + if (present(do_greedy)) then + if (do_greedy) then + local_do_greedy = 1 + else + local_do_greedy = 0 + end if + end if +#ifdef GMPXX_FOUND + par%params_ptr = thiele_pade_mp_aux(n_par, x, y, local_do_greedy) +#endif + + end if end function create_thiele_pade + !> API function to evaluate the Thiele-Pade approximation + !! (potentially using arbitrary precision numbers) + !! + !! @param[in] x - point where the function is evaluated + !! @param[in] params - abstract type to store all parameters + !! @return y - interpolated function value at x function evaluate_thiele_pade_at(par, x) result(y) type(params), intent(in) :: par complex(kind=dp), dimension(:), intent(in) :: x @@ -134,6 +205,8 @@ function evaluate_thiele_pade_at(par, x) result(y) ! initialized? if (par%initialized .eqv. .false.) then print *, "WARNING: pade parameters not initialized" + y(:) = 0.0d0 + return end if ! Compute the number of query points @@ -141,69 +214,38 @@ function evaluate_thiele_pade_at(par, x) result(y) ! Evaluate the Thiele-Pade approximation at the query points do i = 1, num_query - call evaluate_thiele_pade(par%n_par, par%x_ref, x(i), par%a_par, y(i)) + + if (.not. par%multiprecision_used_internally) then + call evaluate_thiele_pade(par%n_par, par%x_ref, x(i), par%a_par, y(i)) + elseif (par%multiprecision_used_internally) then +#ifdef GMPXX_FOUND + y(i) = evaluate_thiele_pade_mp_aux(x(i), par%params_ptr) +#endif + end if + end do end function evaluate_thiele_pade_at + !> deallocate the pade model + !! + !! @param[inout] par - pade model struct subroutine free_params(par) type(params), intent(inout) :: par + if (.not. par%initialized) return + if (allocated(par%a_par)) deallocate(par%a_par) if (allocated(par%x_ref)) deallocate(par%x_ref) - end subroutine free_params - - - - - - !> API function to compute Thiele-Pade parameters using arbitrary precision numbers - !! @param[in] n_par - order of the interpolant - !! @param[in] x_ref - array of the reference points - !! @param[in] y_ref - array of the reference function values - !! @param[out] params - abstract type to store all parameters in arb. prec. representation - !! @param[in] do_greedy - whether to use the default greedy algorithm or the naive one - subroutine create_thiele_pade_mp(n_par, x_ref, y_ref, params, do_greedy) - integer, intent(in) :: n_par - complex(kind=dp), dimension(:), intent(in) :: x_ref, y_ref - type(params_mp), intent(out) :: params - logical, optional, intent(in) :: do_greedy - - ! Internal variables - integer :: local_do_greedy = 1 - - ! use integer bools for interoperability with C - if (present(do_greedy)) then - if (do_greedy) then - local_do_greedy = 1 - else - local_do_greedy = 0 - end if - end if - - ! compute coefficients - params%initialized = .True. - params%params_ptr = thiele_pade_mp_aux(n_par, x_ref, y_ref, local_do_greedy) - - end subroutine create_thiele_pade_mp + if (par%multiprecision_used_internally) then +#ifdef GMPXX_FOUND + call free_pade_model(par%params_ptr) +#endif + end if - !> API function to evaluate the Thiele-Pade approximation using arbitrary precision numbers - !! @param[in] x - point where the function is evaluated - !! @param[out] y - interpolated function value at x - !! @param[in] params - abstract type to store all parameters in arb. prec. representation - subroutine evaluate_thiele_pade_mp(x, y, params) - complex(kind=dp), intent(in) :: x - complex(kind=dp), intent(out) :: y - type(params_mp), intent(in) :: params - - y = cmplx(0.0d0, 0.0d0, kind=8) - if (.not. params%initialized) return - - y = evaluate_thiele_pade_mp_aux(x, params%params_ptr) - - end subroutine evaluate_thiele_pade_mp + end subroutine free_params end module gx_ac diff --git a/GX-AnalyticContinuation/test/test_pade_approximant.f90 b/GX-AnalyticContinuation/test/test_pade_approximant.f90 index fed9649f..b2bf9ff5 100644 --- a/GX-AnalyticContinuation/test/test_pade_approximant.f90 +++ b/GX-AnalyticContinuation/test/test_pade_approximant.f90 @@ -1,5 +1,5 @@ ! ************************************************************************************************** -! Copyright (C) 2020-2023 GreenX library +! Copyright (C) 2020-2024 GreenX library ! This file is distributed under the terms of the APACHE2 License. ! ! ************************************************************************************************** @@ -10,7 +10,8 @@ module test_pade_approximant use kinds, only: dp ! Module being tested - use gx_ac, only: thiele_pade_api, create_thiele_pade_mp, evaluate_thiele_pade_mp, params_mp + use gx_ac, only: thiele_pade_api, create_thiele_pade, evaluate_thiele_pade_at, & + free_params, params use pade_approximant, only: pade implicit none @@ -39,6 +40,8 @@ logical function is_close_complex_dp(a, b, tol) is_close_complex_dp = abs(a - b) <= tolerance end function is_close_complex_dp + + !> Test the Pade interpolant against the function -1 / (x - x0) subroutine test_pade(test) !> Test object @@ -76,6 +79,8 @@ subroutine test_pade(test) end subroutine test_pade + + !> Test the GMP Pade interpolant against the function -1 / (x - x0) subroutine test_pade_mp(test) !> Test object @@ -83,12 +88,12 @@ subroutine test_pade_mp(test) !> N sampling points integer, parameter :: n = 100 - type(params_mp) :: params + type(params) :: params_thiele !> Variable and function, respectively complex(dp), allocatable :: x(:), f(:) !> Pade approximant of f, and its reference value - complex(dp) :: xx - complex(dp) :: f_approx + complex(dp), dimension(:), allocatable :: xx + complex(dp), dimension(:), allocatable :: f_approx complex(dp) :: ref !> Some function center complex(dp), parameter :: x0 = cmplx(2.0_dp, 2.0_dp, kind=dp) @@ -97,28 +102,29 @@ subroutine test_pade_mp(test) integer :: i !> Test setup - allocate(x(n), f(n)) + allocate(x(n), f(n), xx(1), f_approx(1)) do i = 1, n - x(i) = cmplx(i, 0.0_dp, kind=dp) + x(i) = cmplx(i, 0.1_dp, kind=dp) f(i) = -1.0_dp / (x(i) - x0) end do - xx = cmplx(1.0_dp, 1.0_dp, kind=dp) - ref = -1.0_dp / (xx - x0) - - call create_thiele_pade_mp(n, x, f, params) - call evaluate_thiele_pade_mp(xx, f_approx, params) + xx = cmplx(1.5_dp, 0.1_dp, kind=dp) + ref = -1.0_dp / (xx(1) - x0) + params_thiele = create_thiele_pade(n, x, f) + f_approx = evaluate_thiele_pade_at(params_thiele, xx) !> Test execution - call test%assert(is_close(f_approx, ref, tol=tol), name = 'Test Pade GMP ~ -1 / (x - x0)') + call test%assert(is_close(f_approx(1), ref, tol=tol), name = 'Test Pade GMP ~ -1 / (x - x0)') !> Clean-up - deallocate(x) - deallocate(f) + deallocate(x, f, xx, f_approx) + call free_params(params_thiele) end subroutine test_pade_mp + + !> Test the Thiele-Pade interpolant against the function 1 / (-x^2 + 1) which has poles subroutine test_thiele_pade_poles(test) class(unit_test_type), intent(inout) :: test @@ -155,6 +161,8 @@ subroutine test_thiele_pade_poles(test) end subroutine test_thiele_pade_poles + + !> Test the GMP Thiele-Pade interpolant against the function 1 / (-x^2 + 1) which has poles subroutine test_thiele_pade_poles_mp(test) class(unit_test_type), intent(inout) :: test @@ -163,36 +171,38 @@ subroutine test_thiele_pade_poles_mp(test) integer, parameter :: n = 100 !> Variable, function, and parameters, respectively complex(dp), allocatable :: x(:), f(:) - type(params_mp) :: params + type(params) :: params_thiele !> Pade approximant of f, and its reference value - complex(dp):: xx - complex(dp) :: f_approx + complex(dp), dimension(:), allocatable :: xx + complex(dp), dimension(:), allocatable :: f_approx complex(dp) :: ref !> Tolerance real(dp) :: tol = 1.e-7_dp integer :: i !> Test setup - allocate(x(n), f(n)) + allocate(x(n), f(n), xx(1), f_approx(1)) do i = 1, n x(i) = cmplx(i, 0.05_dp, kind=dp) f(i) = 1.0_dp / (-x(i) * x(i) + 1.0_dp) end do xx = cmplx(1.0_dp, 3.0_dp, kind=dp) - ref = 1.0_dp / (-xx * xx + 1.0_dp) - call create_thiele_pade_mp(n, x, f, params) - call evaluate_thiele_pade_mp(xx, f_approx, params) + ref = 1.0_dp / (-xx(1) * xx(1) + 1.0_dp) + params_thiele = create_thiele_pade(n, x, f) + f_approx = evaluate_thiele_pade_at(params_thiele, xx) !> Test execution - call test%assert(is_close(f_approx, ref, tol=tol), name = 'Test GMP Thiele-Pade ~ 1 / (-x^2 + 1)') + call test%assert(is_close(f_approx(1), ref, tol=tol), name = 'Test GMP Thiele-Pade ~ 1 / (-x^2 + 1)') !> Clean-up - deallocate(x) - deallocate(f) + deallocate(x, f, xx, f_approx) + call free_params(params_thiele) end subroutine test_thiele_pade_poles_mp + + !> Test the Thiele-Pade interpolant against the function |x| which has a branch point subroutine test_thiele_pade_abs(test) class(unit_test_type), intent(inout) :: test @@ -237,6 +247,8 @@ subroutine test_thiele_pade_abs(test) end subroutine test_thiele_pade_abs + + !> Test the GMP Thiele-Pade interpolant against the function |x| which has a branch point subroutine test_thiele_pade_abs_mp(test) class(unit_test_type), intent(inout) :: test @@ -248,10 +260,10 @@ subroutine test_thiele_pade_abs_mp(test) real(dp), parameter :: delta_eta = 0.0005_dp !> Variable, function, and parameters, respectively complex(dp), allocatable :: x(:), f(:) - type(params_mp) :: params + type(params) :: params_thiele !> Pade approximant of f, and its reference value - complex(dp) :: xx - complex(dp) :: f_approx + complex(dp), dimension(:), allocatable :: xx + complex(dp), dimension(:), allocatable :: f_approx complex(dp) :: ref !> Test point !> Tolerance @@ -260,7 +272,7 @@ subroutine test_thiele_pade_abs_mp(test) !> Test setup npar = 2 * n - allocate(x(npar), f(npar)) + allocate(x(npar), f(npar), xx(1), f_approx(1)) !> Here we use a Newman grid with 2n points do i = 1, n @@ -271,18 +283,18 @@ subroutine test_thiele_pade_abs_mp(test) f(:) = abs(x(:)) xx = cmplx(0.7_dp, 0.0_dp, kind=dp) - ref = abs(xx) + ref = abs(xx(1)) - call create_thiele_pade_mp(npar, x, f, params) - call evaluate_thiele_pade_mp(xx, f_approx, params) + params_thiele = create_thiele_pade(npar, x, f) + f_approx = evaluate_thiele_pade_at(params_thiele, xx) !> Test execution - call test%assert(is_close(f_approx, ref, tol=tol), name = 'Test GMP Thiele-Pade ~ |x|') + call test%assert(is_close(f_approx(1), ref, tol=tol), name = 'Test GMP Thiele-Pade ~ |x|') !> Clean-up - deallocate(x) - deallocate(f) + deallocate(x, f, xx, f_approx) + call free_params(params_thiele) end subroutine test_thiele_pade_abs_mp