QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Fix branch selection in long-range ewald sum. 

Now everything is consistent except for sigma0


Former-commit-id: 86132a0db0604ed8c56f725eb95c5b22a103b804
This commit is contained in:
Marek Nečada 2018-09-13 01:16:13 +03:00
parent 0ea6aa8413
commit e204242817
3 changed files with 44 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
qpms/ewald.c
View File

@ -139,6 +139,7 @@ int ewald32_sigma0(complex double *result, double *err,
{ {
qpms_csf_result gam; qpms_csf_result gam;
int retval = complex_gamma_inc_e(-0.5, -sq(k/(2*eta)), &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) if (0 != retval)
abort(); abort();
*result = gam.val * c->legendre0[gsl_sf_legendre_array_index(0,0)] / 2 / M_SQRTPI; *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é 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) { for(qpms_l_t j = 0; j <= lMax/2; ++j) {
int retval = complex_gamma_inc_e(0.5-j, z, Gamma_pq+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(); if(!(retval==0 ||retval==GSL_EUNDRFLW)) abort();
} }
// R-DEPENDENT END // R-DEPENDENT END

4
qpms/ewaldsf.c
View File

@ -35,7 +35,7 @@ void IgnoreUnderflowsGSLErrorHandler (const char * reason,
// DLMF 8.7.3 (latter expression) for complex second argument // DLMF 8.7.3 (latter expression) for complex second argument
// BTW if a is large negative, it might take a while to evaluate. // 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) { if (a <= 0 && a == (int) a) {
result->val = NAN + NAN*I; result->val = NAN + NAN*I;
result->err = NAN; 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; result->val = NAN + NAN*I; result->err = NAN;
return retval; 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); ckahanadd(&sum, &sumc, summand);
double summanderr = fabs(fullgamma_ak.err * cabs(summand / fullgamma_ak.val)); double summanderr = fabs(fullgamma_ak.err * cabs(summand / fullgamma_ak.val));
// TODO add also the rounding error // TODO add also the rounding error

55
tests/ewalds.c
View File

@ -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.3, 100, 100, 1., TRIANGULAR_VERTICAL},
{ 6, {1.1, 0.23}, 2.3, 0.97, 2.9, 100, 100, 1., TRIANGULAR_VERTICAL}, { 6, {1.1, 0.23}, 2.3, 0.97, 2.9, 100, 100, 1., TRIANGULAR_VERTICAL},
#endif #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, 160, 1., TRIANGULAR_VERTICAL},
{ 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, 160, 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, 160, 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, 160, 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, 160, 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, 5.5, 20, 160, 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, 5.5, 10, 80, 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, 5.5, 5, 40, 1., TRIANGULAR_VERTICAL},
{ 2, {1.1, 0}, 2.3, 0.97, 0.5, 2, 2, 1., TRIANGULAR_VERTICAL}, { 2, {1.1, 0}, 2.3, 0.97, 5.5, 2, 16, 1., TRIANGULAR_VERTICAL},
{ 2, {0.3, 0}, 2.3, 0.97, 0.5, 2, 2, 1., TRIANGULAR_VERTICAL}, { 2, {1.1, 0}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
{ 2, {2.7, 0}, 2.3, 0.97, 0.5, 2, 2, 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: // end:
// { 0, {0, 0}, 0, 0, 0, 0, 0, 0, 0} // { 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); ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p);
int main() { int main() {
@ -110,15 +131,15 @@ int main() {
qpms_y_t y_conj = qpms_mn2y_sc(-m,n); qpms_y_t y_conj = qpms_mn2y_sc(-m,n);
// y n m sigma_total (err), regsigmas_416 regsigmas_415_recon // 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", 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], 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], 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], r->err_sigmas_short[y],
creal(r->regsigmas_416[y]), cimag(r->regsigmas_416[y]), san(creal(r->regsigmas_416[y])), san(cimag(r->regsigmas_416[y])),
creal(r->sigmas_total[y])+creal(r->sigmas_total[y_conj]), san(creal(r->sigmas_total[y]) + creal(r->sigmas_total[y_conj])),
cimag(r->sigmas_total[y])-cimag(r->sigmas_total[y_conj]) san(cimag(r->sigmas_total[y]) - cimag(r->sigmas_total[y_conj]))
); );
} }
} }