From e204242817ca1fcc5359a294664d48635ec9cdc1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Thu, 13 Sep 2018 01:16:13 +0300 Subject: [PATCH] Fix branch selection in long-range ewald sum. Now everything is consistent except for sigma0 Former-commit-id: 86132a0db0604ed8c56f725eb95c5b22a103b804 --- qpms/ewald.c | 4 ++++ qpms/ewaldsf.c | 4 ++-- tests/ewalds.c | 55 ++++++++++++++++++++++++++++++++++---------------- 3 files changed, 44 insertions(+), 19 deletions(-) diff --git a/qpms/ewald.c b/qpms/ewald.c index d8915e6..bcac6fd 100644 --- a/qpms/ewald.c +++ b/qpms/ewald.c @@ -139,6 +139,7 @@ int ewald32_sigma0(complex double *result, double *err, { qpms_csf_result gam; int retval = complex_gamma_inc_e(-0.5, -sq(k/(2*eta)), &gam); + gam.val = conj(gam.val); // We take the other branch, cf. [Linton, p. 642 in the middle] if (0 != retval) abort(); *result = gam.val * c->legendre0[gsl_sf_legendre_array_index(0,0)] / 2 / M_SQRTPI; @@ -191,6 +192,9 @@ int ewald32_sigma_long_shiftedpoints ( complex double z = csq(gamma_pq*k/(2*eta)); // Když o tom tak přemýšlím, tak tohle je vlastně vždy reálné for(qpms_l_t j = 0; j <= lMax/2; ++j) { int retval = complex_gamma_inc_e(0.5-j, z, Gamma_pq+j); + // we take the other branch, cf. [Linton, p. 642 in the middle]: FIXME instead use the C11 CMPLX macros and fill in -O*I part to z in the line above + if(creal(z) < 0) + Gamma_pq[j].val = conj(Gamma_pq[j].val); //FIXME as noted above if(!(retval==0 ||retval==GSL_EUNDRFLW)) abort(); } // R-DEPENDENT END diff --git a/qpms/ewaldsf.c b/qpms/ewaldsf.c index 0f71725..20879e5 100644 --- a/qpms/ewaldsf.c +++ b/qpms/ewaldsf.c @@ -35,7 +35,7 @@ void IgnoreUnderflowsGSLErrorHandler (const char * reason, // DLMF 8.7.3 (latter expression) for complex second argument // BTW if a is large negative, it might take a while to evaluate. -int cx_gamma_inc_series_e(double a, complex z, qpms_csf_result * result) { +int cx_gamma_inc_series_e(const double a, const complex double z, qpms_csf_result * result) { if (a <= 0 && a == (int) a) { result->val = NAN + NAN*I; result->err = NAN; @@ -62,7 +62,7 @@ int cx_gamma_inc_series_e(double a, complex z, qpms_csf_result * result) { result->val = NAN + NAN*I; result->err = NAN; return retval; } - complex double summand = - cpow(z, k) / fullgamma_ak.val; + complex double summand = - cpow(z, k) / fullgamma_ak.val; // TODO test branch selection here with cimag(z) = -0.0 ckahanadd(&sum, &sumc, summand); double summanderr = fabs(fullgamma_ak.err * cabs(summand / fullgamma_ak.val)); // TODO add also the rounding error diff --git a/tests/ewalds.c b/tests/ewalds.c index 9d7908f..bfd3e2e 100644 --- a/tests/ewalds.c +++ b/tests/ewalds.c @@ -52,17 +52,34 @@ ewaldtest_triang_params paramslist[] = { { 6, {1.1, 0.23}, 2.3, 0.97, 2.3, 100, 100, 1., TRIANGULAR_VERTICAL}, { 6, {1.1, 0.23}, 2.3, 0.97, 2.9, 100, 100, 1., TRIANGULAR_VERTICAL}, #endif - { 2, {1.1, 0.23}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL}, - { 2, {-1.1, -0.23}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL}, - { 2, {0, 1.1}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL}, - { 2, {0, -1.1}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL}, - { 2, {-1.1, 0}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL}, - { 2, {1.1, 0}, 2.3, 0.97, 0.5, 20, 20, 1., TRIANGULAR_VERTICAL}, - { 2, {1.1, 0}, 2.3, 0.97, 0.5, 10, 10, 1., TRIANGULAR_VERTICAL}, - { 2, {1.1, 0}, 2.3, 0.97, 0.5, 5, 5, 1., TRIANGULAR_VERTICAL}, - { 2, {1.1, 0}, 2.3, 0.97, 0.5, 2, 2, 1., TRIANGULAR_VERTICAL}, - { 2, {0.3, 0}, 2.3, 0.97, 0.5, 2, 2, 1., TRIANGULAR_VERTICAL}, - { 2, {2.7, 0}, 2.3, 0.97, 0.5, 2, 2, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0.23}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + { 2, {-1.1, -0.23}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + { 2, {0, 1.1}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + { 2, {0, -1.1}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + { 2, {-1.1, 0}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 5.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 5.5, 10, 80, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 5.5, 5, 40, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 5.5, 2, 16, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 2.5, 10, 80, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 2.5, 5, 40, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 2.5, 2, 16, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 0.5, 10, 80, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 0.5, 5, 40, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 0}, 2.3, 0.97, 0.5, 2, 16, 1., TRIANGULAR_VERTICAL}, + { 2, {0.3, 0}, 2.3, 0.97, 0.5, 2, 16, 1., TRIANGULAR_VERTICAL}, + { 2, {2.7, 1}, 2.3, 0.97, 0.5, 5, 40, 1., TRIANGULAR_VERTICAL}, + { 2, {2.7, 1}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + { 2, {2.7, 1}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + { 2, {2.7, 1}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + { 2, {2.7, 1}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 1}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 1}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 1}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + { 2, {1.1, 1}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL}, + // end: // { 0, {0, 0}, 0, 0, 0, 0, 0, 0, 0} @@ -91,6 +108,10 @@ void dump_points2d_rordered(const points2d_rordered_t *ps, char *filename) { } +static inline double san(double x) { + return fabs(x) < 1e-13 ? 0 : x; +} + ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p); int main() { @@ -110,15 +131,15 @@ int main() { qpms_y_t y_conj = qpms_mn2y_sc(-m,n); // y n m sigma_total (err), regsigmas_416 regsigmas_415_recon printf("%zd %d %d: T:%.16g%+.16gj(%.3g) L:%.16g%+.16gj(%.3g) S:%.16g%+.16gj(%.3g) \n| predict %.16g%+.16gj \n| actual %.16g%+.16gj\n", - y, n, m, creal(r->sigmas_total[y]), cimag(r->sigmas_total[y]), + y, n, m, creal(san(r->sigmas_total[y])), san(cimag(r->sigmas_total[y])), r->err_sigmas_total[y], - creal(r->sigmas_long[y]), cimag(r->sigmas_long[y]), + san(creal(r->sigmas_long[y])), san(cimag(r->sigmas_long[y])), r->err_sigmas_long[y], - creal(r->sigmas_short[y]), cimag(r->sigmas_short[y]), + san(creal(r->sigmas_short[y])), san(cimag(r->sigmas_short[y])), r->err_sigmas_short[y], - creal(r->regsigmas_416[y]), cimag(r->regsigmas_416[y]), - creal(r->sigmas_total[y])+creal(r->sigmas_total[y_conj]), - cimag(r->sigmas_total[y])-cimag(r->sigmas_total[y_conj]) + san(creal(r->regsigmas_416[y])), san(cimag(r->regsigmas_416[y])), + san(creal(r->sigmas_total[y]) + creal(r->sigmas_total[y_conj])), + san(cimag(r->sigmas_total[y]) - cimag(r->sigmas_total[y_conj])) ); } }