Skip to content

Commit

Permalink
add final symmetry support in multi-precision pade routines
Browse files Browse the repository at this point in the history
  • Loading branch information
@moritzleucke
moritzleucke committed Sep 11, 2024
1 parent e1f33fd commit 0ba594b
Show file tree
Hide file tree
Showing 4 changed files with 79 additions and 30 deletions.
57 changes: 43 additions & 14 deletions GX-AnalyticContinuation/api/gx_ac.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,13 +85,14 @@ module gx_ac
!! @param[in] y_ref - array of the reference function values
!! @param[in] do_greedy - whether to use the default greedy algorithm or the naive one
!! @return - pointer to abstract type to store all parameters
function thiele_pade_mp_aux(n_par, x_ref, y_ref, do_greedy, precision) bind(C, name="thiele_pade_mp")
function thiele_pade_mp_aux(n_par, x_ref, y_ref, do_greedy, precision, symmetry) bind(C, name="thiele_pade_mp")
import :: c_double_complex, c_int, c_ptr
integer(c_int), value :: n_par
complex(c_double_complex), dimension(*) :: x_ref
complex(c_double_complex), dimension(*) :: y_ref
integer(c_int), value :: do_greedy
integer(c_int), value :: precision
integer(c_int), value :: symmetry
type(c_ptr) :: thiele_pade_mp_aux
end function thiele_pade_mp_aux

Expand Down Expand Up @@ -178,13 +179,14 @@ type(params) function create_thiele_pade(n_par, x, y, do_greedy, precision, enfo
character(*), optional, intent(in) :: enforce_symmetry

! Internal variables
integer :: local_do_greedy
integer :: c_do_greedy
integer :: c_symmetry_label

! initialize type
par%initialized = .true.
par%n_par = n_par

! precision of arithmetic internally
! precision of internal arithmetic
if (present(precision)) then
if (precision .eq. 64) then
! double precision case
Expand Down Expand Up @@ -213,7 +215,7 @@ type(params) function create_thiele_pade(n_par, x, y, do_greedy, precision, enfo
#endif
end if

! Symmetry
! Symmetry consistency check
if (present(enforce_symmetry)) then
if ((enforce_symmetry.eq."x") &
.or. (enforce_symmetry.eq."y") &
Expand All @@ -222,40 +224,67 @@ type(params) function create_thiele_pade(n_par, x, y, do_greedy, precision, enfo
.or. (enforce_symmetry.eq."odd") &
.or. (enforce_symmetry.eq."conjugate") &
.or. (enforce_symmetry.eq."anti-conjugate") &
.or. (enforce_symmetry.eq."none") &
) then
.or. (enforce_symmetry.eq."none")) then
! character for fortran
par%enforced_symmetry = enforce_symmetry
! integer for c
select case (enforce_symmetry)
case ("y")
c_symmetry_label = 1
case ("x")
c_symmetry_label = 2
case ("xy")
c_symmetry_label = 3
case ("even")
c_symmetry_label = 4
case ("odd")
c_symmetry_label = 5
case ("conjugate")
c_symmetry_label = 6
case ("anti-conjugate")
c_symmetry_label = 7
case ("none")
c_symmetry_label = 0
case default
print *, "symmetry not known!"
stop
end select
else
print *, "*** create_thiele_pade: enorce_symmetry=", enforce_symmetry, &
" not known or not supported! Aborting..."
stop
end if
else
par%enforced_symmetry = "none"
c_symmetry_label = 0
end if

! greedy algorithm
if (present(do_greedy)) then
! actual bool for fortran
par%use_greedy = do_greedy
! use integer bools for interoperability with C
if (do_greedy) then
c_do_greedy = 1
else
c_do_greedy = 0
end if
else
par%use_greedy = .true.
c_do_greedy = 0
end if

! create the pade model
if (.not. par%multiprecision_used_internally) then
allocate(par%a_par(n_par))
allocate(par%x_ref(size(x)))
par%x_ref(:) = x
call thiele_pade(n_par, par%x_ref, y, par%a_par, par%use_greedy, par%enforced_symmetry)
call thiele_pade(n_par, par%x_ref, y, par%a_par, &
par%use_greedy, par%enforced_symmetry)
elseif (par%multiprecision_used_internally) then
! use integer bools for interoperability with C
if (par%use_greedy) then
local_do_greedy = 1
else
local_do_greedy = 0
end if
#ifdef GMPXX_FOUND
par%params_ptr = thiele_pade_mp_aux(n_par, x, y, local_do_greedy, par%precision)
par%params_ptr = thiele_pade_mp_aux(n_par, x, y, c_do_greedy, &
par%precision, c_symmetry_label)
#endif
end if

Expand Down
14 changes: 7 additions & 7 deletions GX-AnalyticContinuation/src/Symmetry_pade.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,35 +38,35 @@ std::complex<double> Symmetry_pade::apply_sym_x(std::complex<double> x) {
break;
}
case sym_y: {
x_symm = std::complex<double>(abs(x.real()), x.imag());
x_symm = std::complex<double>(std::abs(x.real()), x.imag());
break;
}
case sym_x: {
x_symm = std::complex<double>(x.real(), abs(x.imag()));
x_symm = std::complex<double>(x.real(), std::abs(x.imag()));
break;
}
case sym_xy: {
x_symm = std::complex<double>(abs(x.real()), abs(x.imag()));
x_symm = std::complex<double>(std::abs(x.real()), std::abs(x.imag()));
break;
}
case sym_even: {
double new_im = std::copysign(1.0, x.real()) * x.imag();
x_symm = std::complex<double>(abs(x.real()), new_im);
x_symm = std::complex<double>(std::abs(x.real()), new_im);
break;
}
case sym_odd: {
double new_im = std::copysign(1.0, x.real()) * x.imag();
x_symm = std::complex<double>(abs(x.real()), new_im);
x_symm = std::complex<double>(std::abs(x.real()), new_im);
break;
}
case sym_conjugate: {
double new_im = std::copysign(1.0, x.real()) * x.imag();
x_symm = std::complex<double>(abs(x.real()), new_im);
x_symm = std::complex<double>(std::abs(x.real()), new_im);
break;
}
case sym_anti_conjugate: {
double new_im = std::copysign(1.0, x.real()) * x.imag();
x_symm = std::complex<double>(abs(x.real()), new_im);
x_symm = std::complex<double>(std::abs(x.real()), new_im);
break;
}
}
Expand Down
31 changes: 23 additions & 8 deletions GX-AnalyticContinuation/src/pade_mp.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,8 +69,11 @@ std::complex<double> evaluate_thiele_pade_mp(const std::complex<double> x,
pade_model *params) {
mpf_set_default_prec(params->precision);

// projection of input x to enforce symmetry in pade
std::complex<double> x_projected = params->symmetry.apply_sym_x(x);

// conversion from double to gmp
ComplexGMP x_mp(x);
ComplexGMP x_mp(x_projected);

// Define constants
const ComplexGMP c_one(std::complex<double>(1.0, 0.0));
Expand All @@ -93,8 +96,12 @@ std::complex<double> evaluate_thiele_pade_mp(const std::complex<double> x,
}
}

return std::complex<double>(acoef[params->n_par - 1] /
bcoef[params->n_par - 1]);
// evalueate projected y
std::complex<double> y = std::complex<double>(acoef[params->n_par - 1] /
bcoef[params->n_par - 1]);
// project y back to get symmetric pade interpolant
std::complex<double> y_projected_back = params->symmetry.apply_sym_y(x, y);
return y_projected_back;
}

/// @brief Gets the Pade approximant of a meromorphic function F
Expand All @@ -103,18 +110,26 @@ std::complex<double> evaluate_thiele_pade_mp(const std::complex<double> x,
/// @param y_ref array of the reference function values
/// @param do_greedy whether to use the default greedy algorithm or the naive
/// @param precision floating point arythmetic precision in bits (!! not bytes!!)
/// @param symmetry symmetry label as int
/// @return a pointer to the struct holding all model parameters
pade_model *thiele_pade_mp(int n_par, const std::complex<double> *x_ref,
const std::complex<double> *y_ref, int do_greedy, int precision) {
const std::complex<double> *y_ref, int do_greedy,
int precision, int symmetry) {
// set floating point precision of GMP lib
mpf_set_default_prec(precision);

// store settings in struct
pade_model *params = new pade_model;
params->n_par = n_par;
params->precision = precision;
params->symmetry.set_symmetry(1);
params->symmetry.set_symmetry(symmetry);

// apply projection to enforce symmetry of pade model
std::vector<std::complex<double>> x_projected(n_par), y_projected(n_par);
for (int i_par = 0; i_par < n_par; i_par++) {
x_projected[i_par] = params->symmetry.apply_sym_x(x_ref[i_par]);
y_projected[i_par] = params->symmetry.apply_sym_y(x_ref[i_par], y_ref[i_par]);
}

// Define constants
const ComplexGMP c_one(std::complex<double>(1.0, 0.0));
Expand All @@ -136,9 +151,9 @@ pade_model *thiele_pade_mp(int n_par, const std::complex<double> *x_ref,

// initialize variables
for (int i_par = 0; i_par < n_par; i_par++) {
x_ref_mp[i_par] = x_ref[i_par];
x[i_par] = x_ref[i_par];
y_ref_mp.push_back(y_ref[i_par]);
x_ref_mp[i_par] = x_projected[i_par];
x[i_par] = x_projected[i_par];
y_ref_mp.push_back(y_projected[i_par]);
g_func.push_back(std::vector<ComplexGMP>());
a_par_mp[i_par] = ComplexGMP();
n_rem_idx.push_back(i_par);
Expand Down
7 changes: 6 additions & 1 deletion GX-AnalyticContinuation/src/pade_mp.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,15 +26,20 @@ struct pade_model {

// function prototype for fortran interface
extern "C" {

std::complex<double> evaluate_thiele_pade_mp(const std::complex<double> x,
pade_model *params_ptr);

pade_model *thiele_pade_mp(int n_par, const std::complex<double> *x_ref,
const std::complex<double> *y_ref, int do_greedy, int precision);
const std::complex<double> *y_ref, int do_greedy,
int precision, int symmetry);

void free_pade_model(pade_model *model) {
if (model != nullptr) {
delete[] model->a_par;
delete[] model->xref;
delete model;
}
}

}

0 comments on commit 0ba594b

Please sign in to comment.