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

Free unfinished PGens from Ewald sums. 

Also replace some bare `if (xxx) abort();` with
`QPMS_ENSURE_SUCCESS(xxx)`.


Former-commit-id: d4712292fe0f9fb397cf3a490131bcb37d542fa6
This commit is contained in:
Marek Nečada 2019-10-01 11:26:07 +03:00
parent 7e74cb100c
commit c7f2e32ee4
2 changed files with 51 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
qpms/ewald.h
View File

@ -178,6 +178,12 @@ int ewald3_sigma_short(
/** Ignored apart from asserts and possible optimisations, as the SR formula stays the same. */ /** Ignored apart from asserts and possible optimisations, as the SR formula stays the same. */
const LatticeDimensionality latdim, const LatticeDimensionality latdim,
/// Lattice point generator for the direct Bravais lattice. /// Lattice point generator for the direct Bravais lattice.
/** There is a possibility that the whole PGen is not consumed
* (this might happen if the summand start to be consistently smaller
* than the (partial) sums * DBL_EPSILON.
* In such case, it is the responsibility of the caller to deallocate
* the generator.
*/
PGen *pgen_R, PGen *pgen_R,
/// Indicates whether pgen_R already generates shifted points. /// Indicates whether pgen_R already generates shifted points.
/** If false, the behaviour corresponds to the old ewald32_sigma_short_points_and_shift(), /** If false, the behaviour corresponds to the old ewald32_sigma_short_points_and_shift(),
@ -204,6 +210,12 @@ int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, dep
/// Lattice dimensionality. /// Lattice dimensionality.
const LatticeDimensionality latdim, const LatticeDimensionality latdim,
/// Lattice point generator for the reciprocal lattice. /// Lattice point generator for the reciprocal lattice.
/** There is a possibility that the whole PGen is not consumed
* (this might happen if the summand start to be consistently smaller
* than the (partial) sums * DBL_EPSILON.
* In such case, it is the responsibility of the caller to deallocate
* the generator.
*/
PGen *pgen_K, PGen *pgen_K,
/// Indicates whether pgen_K already generates shifted points. /// Indicates whether pgen_K already generates shifted points.
/** If false, the behaviour corresponds to the old ewald32_sigma_long_points_and_shift(), /** If false, the behaviour corresponds to the old ewald32_sigma_long_points_and_shift(),

81
qpms/translations.c
View File

@ -10,7 +10,6 @@
#include "tiny_inlines.h" #include "tiny_inlines.h"
#include "assert_cython_workaround.h" #include "assert_cython_workaround.h"
#include "kahansum.h" #include "kahansum.h"
#include <stdlib.h> //abort()
#include <gsl/gsl_sf_coupling.h> #include <gsl/gsl_sf_coupling.h>
#include "qpms_error.h" #include "qpms_error.h"
#include "normalisation.h" #include "normalisation.h"
@ -145,13 +144,13 @@ complex double qpms_trans_single_A(qpms_normalisation_t norm,
double a1q[qmax+1]; double a1q[qmax+1];
int err; int err;
gaunt_xu(-m,n,mu,nu,qmax,a1q,&err); gaunt_xu(-m,n,mu,nu,qmax,a1q,&err);
QPMS_ENSURE_SUCCESS(err);
double a1q0 = a1q[0]; double a1q0 = a1q[0];
if (err) abort();
double leg[gsl_sf_legendre_array_n(n+nu)]; double leg[gsl_sf_legendre_array_n(n+nu)];
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg)) abort(); QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,costheta,-1,leg));
complex double bes[n+nu+1]; complex double bes[n+nu+1];
if (qpms_sph_bessel_fill(J, n+nu, kdlj.r, bes)) abort(); QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J, n+nu, kdlj.r, bes));
complex double sum = 0; complex double sum = 0;
for(int q = 0; q <= qmax; ++q) { for(int q = 0; q <= qmax; ++q) {
int p = n+nu-2*q; int p = n+nu-2*q;
@ -219,16 +218,18 @@ complex double qpms_trans_single_B(qpms_normalisation_t norm,
a2q0 = 1; a2q0 = 1;
} }
else { else {
gaunt_xu(-m-1,n+1,mu+1,nu,q2max,a2q,&err); if (err) abort(); gaunt_xu(-m-1,n+1,mu+1,nu,q2max,a2q,&err);
QPMS_ENSURE_SUCCESS(err);
a2q0 = a2q[0]; a2q0 = a2q[0];
} }
gaunt_xu(-m,n+1,mu,nu,Qmax,a3q,&err); if (err) abort(); gaunt_xu(-m,n+1,mu,nu,Qmax,a3q,&err);
QPMS_ENSURE_SUCCESS(err);
a3q0 = a3q[0]; a3q0 = a3q[0];
double leg[gsl_sf_legendre_array_n(n+nu+1)]; double leg[gsl_sf_legendre_array_n(n+nu+1)];
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg)) abort(); QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,costheta,-1,leg));
complex double bes[n+nu+2]; complex double bes[n+nu+2];
if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bes)) abort(); QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bes));
complex double sum = 0; complex double sum = 0;
for (int q = 0; q <= realQmax; ++q) { for (int q = 0; q <= realQmax; ++q) {
@ -317,7 +318,7 @@ static void qpms_trans_calculator_multipliers_A(
double a1q[qmax+1]; double a1q[qmax+1];
int err; int err;
gaunt_xu(-m,n,mu,nu,qmax,a1q,&err); gaunt_xu(-m,n,mu,nu,qmax,a1q,&err);
if (err) abort(); QPMS_ENSURE_SUCCESS(err);
double a1q0 = a1q[0]; double a1q0 = a1q[0];
double normlogfac = qpms_trans_normlogfac(norm,m,n,mu,nu); double normlogfac = qpms_trans_normlogfac(norm,m,n,mu,nu);
@ -502,8 +503,8 @@ qpms_trans_calculator
c->e3c = qpms_ewald3_constants_init(2 * lMax + 1, /*csphase*/ qpms_normalisation_t_csphase(normalisation)); c->e3c = qpms_ewald3_constants_init(2 * lMax + 1, /*csphase*/ qpms_normalisation_t_csphase(normalisation));
#endif #endif
c->legendre0 = malloc(gsl_sf_legendre_array_n(2*lMax+1) * sizeof(double)); c->legendre0 = malloc(gsl_sf_legendre_array_n(2*lMax+1) * sizeof(double));
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*lMax+1, QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*lMax+1,
0,-1,c->legendre0)) abort(); // TODO maybe use some "precise" analytical formula instead? 0,-1,c->legendre0)); // TODO maybe use some "precise" analytical formula instead?
return c; return c;
} }
@ -540,9 +541,9 @@ complex double qpms_trans_calculator_get_A_buf(const qpms_trans_calculator *c,
int csphase = qpms_normalisation_t_csphase(c->normalisation); int csphase = qpms_normalisation_t_csphase(c->normalisation);
double costheta = cos(kdlj.theta); double costheta = cos(kdlj.theta);
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu, QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu,
costheta,csphase,legendre_buf)) abort(); costheta,csphase,legendre_buf));
if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort(); QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf));
return qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu, return qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu,
kdlj,r_ge_d,J,bessel_buf,legendre_buf); kdlj,r_ge_d,J,bessel_buf,legendre_buf);
} }
@ -579,9 +580,9 @@ complex double qpms_trans_calculator_get_B_buf(const qpms_trans_calculator *c,
return NAN+I*NAN; return NAN+I*NAN;
int csphase = qpms_normalisation_t_csphase(c->normalisation); int csphase = qpms_normalisation_t_csphase(c->normalisation);
double costheta = cos(kdlj.theta); double costheta = cos(kdlj.theta);
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1, QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,
costheta,csphase,legendre_buf)) abort(); costheta,csphase,legendre_buf));
if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort(); QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf));
return qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu, return qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu,
kdlj,r_ge_d,J,bessel_buf,legendre_buf); kdlj,r_ge_d,J,bessel_buf,legendre_buf);
} }
@ -599,9 +600,9 @@ int qpms_trans_calculator_get_AB_buf_p(const qpms_trans_calculator *c,
return 0; return 0;
} }
double costheta = cos(kdlj.theta); double costheta = cos(kdlj.theta);
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1, QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1,
costheta,-1,legendre_buf)) abort(); costheta,-1,legendre_buf));
if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort(); QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf));
*Adest = qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu, *Adest = qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu,
kdlj,r_ge_d,J,bessel_buf,legendre_buf); kdlj,r_ge_d,J,bessel_buf,legendre_buf);
*Bdest = qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu, *Bdest = qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu,
@ -627,9 +628,9 @@ int qpms_trans_calculator_get_AB_arrays_buf(const qpms_trans_calculator *c,
} }
{ {
double costheta = cos(kdlj.theta); double costheta = cos(kdlj.theta);
if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1, QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1,
costheta,-1,legendre_buf)) abort(); costheta,-1,legendre_buf));
if (qpms_sph_bessel_fill(J, 2*c->lMax+1, kdlj.r, bessel_buf)) abort(); QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J, 2*c->lMax+1, kdlj.r, bessel_buf));
size_t desti = 0, srci = 0; size_t desti = 0, srci = 0;
for (int n = 1; n <= c->lMax; ++n) for (int m = -n; m <= n; ++m) { for (int n = 1; n <= c->lMax; ++n) for (int m = -n; m <= n; ++m) {
for (int nu = 1; nu <= c->lMax; ++nu) for (int mu = -nu; mu <= nu; ++mu) { for (int nu = 1; nu <= c->lMax; ++nu) for (int mu = -nu; mu <= nu; ++mu) {
@ -829,14 +830,11 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
serr_total = malloc(sizeof(double)*nelem2_sc); serr_total = malloc(sizeof(double)*nelem2_sc);
} else serr_short = serr_long = serr_total = NULL; } else serr_short = serr_long = serr_total = NULL;
int retval; QPMS_ENSURE_SUCCESS(ewald31z_sigma_long_points_and_shift(sigmas_long, serr_long, //DIFF21
retval = ewald31z_sigma_long_points_and_shift(sigmas_long, serr_long, //DIFF21 c->e3c, eta, k, unitcell_area, nKpoints, Kpoints, beta, particle_shift));
c->e3c, eta, k, unitcell_area, nKpoints, Kpoints, beta, particle_shift);
if (retval) abort();
retval = ewald31z_sigma_short_points_and_shift(sigmas_short, serr_short, //DIFF21 QPMS_ENSURE_SUCCESS(ewald31z_sigma_short_points_and_shift(sigmas_short, serr_short, //DIFF21
c->e3c, eta, k, nRpoints, Rpoints, beta, particle_shift); c->e3c, eta, k, nRpoints, Rpoints, beta, particle_shift));
if (retval) abort();
for(qpms_y_t y = 0; y < nelem2_sc; ++y) for(qpms_y_t y = 0; y < nelem2_sc; ++y)
sigmas_total[y] = sigmas_short[y] + sigmas_long[y]; sigmas_total[y] = sigmas_short[y] + sigmas_long[y];
@ -845,8 +843,7 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
complex double sigma0 = 0; double sigma0_err = 0; complex double sigma0 = 0; double sigma0_err = 0;
if (do_sigma0) { if (do_sigma0) {
retval = ewald31z_sigma0(&sigma0, &sigma0_err, c->e3c, eta, k);//DIFF21 QPMS_ENSURE_SUCCESS(ewald31z_sigma0(&sigma0, &sigma0_err, c->e3c, eta, k));
if(retval) abort();
const qpms_l_t y = qpms_mn2y_sc(0,0); const qpms_l_t y = qpms_mn2y_sc(0,0);
sigmas_total[y] += sigma0; sigmas_total[y] += sigma0;
if(doerr) serr_total[y] += sigma0_err; if(doerr) serr_total[y] += sigma0_err;
@ -946,7 +943,7 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
const double unitcell_area = l2d_unitcell_area(b1, b2); const double unitcell_area = l2d_unitcell_area(b1, b2);
cart2_t rb1, rb2; // reciprocal basis cart2_t rb1, rb2; // reciprocal basis
if (QPMS_SUCCESS != l2d_reciprocalBasis2pi(b1, b2, &rb1, &rb2)) abort(); QPMS_ENSURE_SUCCESS(l2d_reciprocalBasis2pi(b1, b2, &rb1, &rb2));
PGen Rgen = PGen_xyWeb_new(b1, b2, BASIS_RTOL, PGen Rgen = PGen_xyWeb_new(b1, b2, BASIS_RTOL,
#ifdef GEN_RSHIFTEDPOINTS #ifdef GEN_RSHIFTEDPOINTS
@ -963,26 +960,27 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
#endif #endif
0, true, maxK, false); 0, true, maxK, false);
int retval; QPMS_ENSURE_SUCCESS(ewald3_sigma_long(sigmas_long, serr_long, c->e3c, eta, k,
retval = ewald3_sigma_long(sigmas_long, serr_long, c->e3c, eta, k,
unitcell_area, LAT_2D_IN_3D_XYONLY, &Kgen, unitcell_area, LAT_2D_IN_3D_XYONLY, &Kgen,
#ifdef GEN_KSHIFTEDPOINTS #ifdef GEN_KSHIFTEDPOINTS
true, true,
#else #else
false, false,
#endif #endif
cart22cart3xy(beta), cart22cart3xy(particle_shift)); cart22cart3xy(beta), cart22cart3xy(particle_shift)));
if (retval) abort(); if(Kgen.stateData) // PGen not consumed entirely (converged earlier)
PGen_destroy(&Kgen);
retval = ewald3_sigma_short(sigmas_short, serr_short, c->e3c, eta, k, QPMS_ENSURE_SUCCESS(ewald3_sigma_short(sigmas_short, serr_short, c->e3c, eta, k,
LAT_2D_IN_3D_XYONLY, &Rgen, LAT_2D_IN_3D_XYONLY, &Rgen,
#ifdef GEN_RSHIFTEDPOINTS #ifdef GEN_RSHIFTEDPOINTS
true, true,
#else #else
false, false,
#endif #endif
cart22cart3xy(beta), cart22cart3xy(particle_shift)); cart22cart3xy(beta), cart22cart3xy(particle_shift)));
if (retval) abort(); if(Rgen.stateData) // PGen not consumed entirely (converged earlier)
PGen_destroy(&Rgen);
for(qpms_y_t y = 0; y < nelem2_sc; ++y) for(qpms_y_t y = 0; y < nelem2_sc; ++y)
sigmas_total[y] = sigmas_short[y] + sigmas_long[y]; sigmas_total[y] = sigmas_short[y] + sigmas_long[y];
@ -991,8 +989,7 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
complex double sigma0 = 0; double sigma0_err = 0; complex double sigma0 = 0; double sigma0_err = 0;
if (do_sigma0) { if (do_sigma0) {
retval = ewald3_sigma0(&sigma0, &sigma0_err, c->e3c, eta, k); QPMS_ENSURE_SUCCESS(ewald3_sigma0(&sigma0, &sigma0_err, c->e3c, eta, k));
if(retval) abort();
const qpms_l_t y = qpms_mn2y_sc(0,0); const qpms_l_t y = qpms_mn2y_sc(0,0);
sigmas_total[y] += sigma0; sigmas_total[y] += sigma0;
if(doerr) serr_total[y] += sigma0_err; if(doerr) serr_total[y] += sigma0_err;