Scalar m,n indexing to allow n = 0.

Former-commit-id: 4b6f2f3611a00e2019b54bfe73c407baafdd8355
This commit is contained in:
Marek Nečada 2018-09-05 12:39:02 +00:00
parent 16ce3a6ba8
commit e11f995b52
4 changed files with 88 additions and 43 deletions

View File

@ -66,15 +66,15 @@ qpms_ewald32_constants_t *qpms_ewald32_constants_init(const qpms_l_t lMax /*, co
qpms_ewald32_constants_t *c = malloc(sizeof(qpms_ewald32_constants_t)); qpms_ewald32_constants_t *c = malloc(sizeof(qpms_ewald32_constants_t));
//if (c == NULL) return NULL; // Do I really want to do this? //if (c == NULL) return NULL; // Do I really want to do this?
c->lMax = lMax; c->lMax = lMax;
c->nelem = qpms_lMax2nelem(lMax); c->nelem_sc = qpms_lMax2nelem_sc(lMax);
c->s1_jMaxes = malloc(c->nelem * sizeof(qpms_l_t)); c->s1_jMaxes = malloc(c->nelem_sc * sizeof(qpms_l_t));
c->s1_constfacs = malloc(c->nelem * sizeof(complex double *)); c->s1_constfacs = malloc(c->nelem_sc * sizeof(complex double *));
//if (c->s1_jMaxes == NULL) return NULL; //if (c->s1_jMaxes == NULL) return NULL;
//determine sizes //determine sizes
size_t s1_constfacs_sz = 0; size_t s1_constfacs_sz = 0;
for (qpms_y_t y = 0; y < c->nelem; ++y) { for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
qpms_l_t n; qpms_m_t m; qpms_y2mn_p(y, &m, &n); qpms_l_t n; qpms_m_t m; qpms_y2mn_sc_p(y, &m, &n);
if ((m + n) % 2 == 0) if ((m + n) % 2 == 0)
s1_constfacs_sz += 1 + (c->s1_jMaxes[y] = (n-abs(m))/2); s1_constfacs_sz += 1 + (c->s1_jMaxes[y] = (n-abs(m))/2);
else else
@ -82,10 +82,10 @@ qpms_ewald32_constants_t *qpms_ewald32_constants_init(const qpms_l_t lMax /*, co
} }
c->s1_constfacs[0]; //WTF??? c->s1_constfacs[0]; //WTF???
c->s1_constfacs_base = malloc(c->nelem * sizeof(complex double)); c->s1_constfacs_base = malloc(c->nelem_sc * sizeof(complex double));
size_t s1_constfacs_sz_cumsum = 0; size_t s1_constfacs_sz_cumsum = 0;
for (qpms_y_t y = 0; y < c->nelem; ++y) { for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
qpms_l_t n; qpms_m_t m; qpms_y2mn_p(y, &m, &n); qpms_l_t n; qpms_m_t m; qpms_y2mn_sc_p(y, &m, &n);
if ((m + n) % 2 == 0) { if ((m + n) % 2 == 0) {
c->s1_constfacs[y] = c->s1_constfacs_base + s1_constfacs_sz_cumsum; c->s1_constfacs[y] = c->s1_constfacs_base + s1_constfacs_sz_cumsum;
// and here comes the actual calculation // and here comes the actual calculation
@ -151,7 +151,7 @@ int ewald32_sigma0(complex double *result, double *err,
int ewald32_sigma_long_shiftedpoints ( int ewald32_sigma_long_shiftedpoints (
complex double *target, // must be c->nelem long complex double *target, // must be c->nelem_sc long
double *err, double *err,
const qpms_ewald32_constants_t *c, const qpms_ewald32_constants_t *c,
const double eta, const double k, const double unitcell_area, const double eta, const double k, const double unitcell_area,
@ -159,16 +159,16 @@ int ewald32_sigma_long_shiftedpoints (
const point2d particle_shift // target - src const point2d particle_shift // target - src
) )
{ {
const qpms_y_t nelem = c->nelem; const qpms_y_t nelem_sc = c->nelem_sc;
const qpms_l_t lMax = c->lMax; const qpms_l_t lMax = c->lMax;
// Manual init of the ewald summation targets // Manual init of the ewald summation targets
complex double *target_c = calloc(nelem, sizeof(complex double)); complex double *target_c = calloc(nelem_sc, sizeof(complex double));
memset(target, 0, nelem * sizeof(complex double)); memset(target, 0, nelem_sc * sizeof(complex double));
double *err_c = NULL; double *err_c = NULL;
if (err) { if (err) {
err_c = calloc(nelem, sizeof(double)); err_c = calloc(nelem_sc, sizeof(double));
memset(err, 0, nelem * sizeof(double)); memset(err, 0, nelem_sc * sizeof(double));
} }
const double commonfac = 1/(k*k*unitcell_area); // used in the very end (CFC) const double commonfac = 1/(k*k*unitcell_area); // used in the very end (CFC)
@ -197,11 +197,11 @@ int ewald32_sigma_long_shiftedpoints (
// TODO optimisations: all the j-dependent powers can be done for each j only once, stored in array // TODO optimisations: all the j-dependent powers can be done for each j only once, stored in array
// and just fetched for each n, m pair // and just fetched for each n, m pair
for(qpms_l_t n = 1; n <= lMax; ++n) for(qpms_l_t n = 0; n <= lMax; ++n)
for(qpms_m_t m = -n; m <= n; ++m) { for(qpms_m_t m = -n; m <= n; ++m) {
if((m+n) % 2 != 0) // odd coefficients are zero. if((m+n) % 2 != 0) // odd coefficients are zero.
continue; continue;
qpms_y_t y = qpms_mn2y(m, n); qpms_y_t y = qpms_mn2y_sc(m, n);
complex double e_imalpha_pq = cexp(I*m*arg_pq); complex double e_imalpha_pq = cexp(I*m*arg_pq);
complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c); complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c);
double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors? double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors?
@ -225,10 +225,10 @@ int ewald32_sigma_long_shiftedpoints (
free(err_c); free(err_c);
free(target_c); free(target_c);
for(qpms_y_t y = 0; y < nelem; ++y) // CFC common factor from above for(qpms_y_t y = 0; y < nelem_sc; ++y) // CFC common factor from above
target[y] *= commonfac; target[y] *= commonfac;
if(err) if(err)
for(qpms_y_t y = 0; y < nelem; ++y) for(qpms_y_t y = 0; y < nelem_sc; ++y)
err[y] *= commonfac; err[y] *= commonfac;
return 0; return 0;
} }
@ -270,7 +270,7 @@ static int ewald32_sr_integral(double r, double k, int n, double eta,
} }
int ewald32_sigma_short_shiftedpoints( int ewald32_sigma_short_shiftedpoints(
complex double *target, // must be c->nelem long complex double *target, // must be c->nelem_sc long
double *err, double *err,
const qpms_ewald32_constants_t *c, // N.B. not too useful here const qpms_ewald32_constants_t *c, // N.B. not too useful here
const double eta, const double k, const double eta, const double k,
@ -279,18 +279,18 @@ int ewald32_sigma_short_shiftedpoints(
const point2d particle_shift // used only in the very end to multiply it by the phase const point2d particle_shift // used only in the very end to multiply it by the phase
) )
{ {
const qpms_y_t nelem = c->nelem; const qpms_y_t nelem_sc = c->nelem_sc;
const qpms_l_t lMax = c->lMax; const qpms_l_t lMax = c->lMax;
gsl_integration_workspace *workspace = gsl_integration_workspace *workspace =
gsl_integration_workspace_alloc(INTEGRATION_WORKSPACE_LIMIT); gsl_integration_workspace_alloc(INTEGRATION_WORKSPACE_LIMIT);
// Manual init of the ewald summation targets // Manual init of the ewald summation targets
complex double *target_c = calloc(nelem, sizeof(complex double)); complex double *target_c = calloc(nelem_sc, sizeof(complex double));
memset(target, 0, nelem * sizeof(complex double)); memset(target, 0, nelem_sc * sizeof(complex double));
double *err_c = NULL; double *err_c = NULL;
if (err) { if (err) {
err_c = calloc(nelem, sizeof(double)); err_c = calloc(nelem_sc, sizeof(double));
memset(err, 0, nelem * sizeof(double)); memset(err, 0, nelem_sc * sizeof(double));
} }
@ -304,7 +304,7 @@ int ewald32_sigma_short_shiftedpoints(
double Rpq_shifted_arg = atan2(Rpq_shifted.x, Rpq_shifted.y); // POINT-DEPENDENT double Rpq_shifted_arg = atan2(Rpq_shifted.x, Rpq_shifted.y); // POINT-DEPENDENT
complex double e_beta_Rpq = cexp(I*cart2_dot(beta, Rpq_shifted)); // POINT-DEPENDENT complex double e_beta_Rpq = cexp(I*cart2_dot(beta, Rpq_shifted)); // POINT-DEPENDENT
for(qpms_l_t n = 1; n <= lMax; ++n) { for(qpms_l_t n = 0; n <= lMax; ++n) {
double complex prefacn = - I * pow(2./k, n+1) * M_2_SQRTPI / 2; // TODO put outside the R-loop and multiply in the end double complex prefacn = - I * pow(2./k, n+1) * M_2_SQRTPI / 2; // TODO put outside the R-loop and multiply in the end
double R_pq_pown = pow(r_pq_shifted, n); double R_pq_pown = pow(r_pq_shifted, n);
// TODO the integral here // TODO the integral here
@ -317,7 +317,7 @@ int ewald32_sigma_short_shiftedpoints(
continue; // nothing needed, already done by memset continue; // nothing needed, already done by memset
complex double e_imf = cexp(I*m*Rpq_shifted_arg); complex double e_imf = cexp(I*m*Rpq_shifted_arg);
double leg = c->legendre0[gsl_sf_legendre_array_index(n, m)]; double leg = c->legendre0[gsl_sf_legendre_array_index(n, m)];
qpms_y_t y = qpms_mn2y(m,n); qpms_y_t y = qpms_mn2y_sc(m,n);
if(err) if(err)
kahanadd(err + y, err_c + y, cabs(leg * (prefacn / I) * R_pq_pown kahanadd(err + y, err_c + y, cabs(leg * (prefacn / I) * R_pq_pown
* interr)); // TODO include also other errors * interr)); // TODO include also other errors
@ -339,7 +339,7 @@ int ewald32_sigma_short_shiftedpoints(
int ewald32_sigma_long_points_and_shift ( int ewald32_sigma_long_points_and_shift (
complex double *target_sigmalr_y, // must be c->nelem long complex double *target_sigmalr_y, // must be c->nelem_sc long
const qpms_ewald32_constants_t *c, const qpms_ewald32_constants_t *c,
double eta, double k, double unitcell_area, double eta, double k, double unitcell_area,
size_t npoints, const point2d *Kpoints, size_t npoints, const point2d *Kpoints,
@ -347,21 +347,21 @@ int ewald32_sigma_long_points_and_shift (
point2d particle_shift point2d particle_shift
); );
int ewald32_sigma_long_shiftedpoints_rordered( int ewald32_sigma_long_shiftedpoints_rordered(
complex double *target_sigmalr_y, // must be c->nelem long complex double *target_sigmalr_y, // must be c->nelem_sc long
const qpms_ewald32_constants_t *c, const qpms_ewald32_constants_t *c,
double eta, double k, double unitcell_area, double eta, double k, double unitcell_area,
const points2d_rordered_t *Kpoints_plus_beta_rordered, const points2d_rordered_t *Kpoints_plus_beta_rordered,
point2d particle_shift point2d particle_shift
); );
int ewald32_sigma_short_points_and_shift( int ewald32_sigma_short_points_and_shift(
complex double *target_sigmasr_y, // must be c->nelem long complex double *target_sigmasr_y, // must be c->nelem_sc long
const qpms_ewald32_constants_t *c, // N.B. not too useful here const qpms_ewald32_constants_t *c, // N.B. not too useful here
double eta, double k, double eta, double k,
size_t npoints, const point2d *Rpoints, size_t npoints, const point2d *Rpoints,
point2d particle_shift point2d particle_shift
); );
int ewald32_sigma_short_points_rordered( int ewald32_sigma_short_points_rordered(
complex double *target_sigmasr_y, // must be c->nelem long complex double *target_sigmasr_y, // must be c->nelem_sc long
const qpms_ewald32_constants_t *c, // N.B. not too useful here const qpms_ewald32_constants_t *c, // N.B. not too useful here
double eta, double k, double eta, double k,
const points2d_rordered_t *Rpoints_plus_particle_shift_rordered, const points2d_rordered_t *Rpoints_plus_particle_shift_rordered,

View File

@ -27,7 +27,7 @@
/* Object holding the constant factors from [1, (4.5)] */ /* Object holding the constant factors from [1, (4.5)] */
typedef struct { typedef struct {
qpms_l_t lMax; qpms_l_t lMax;
qpms_y_t nelem; qpms_y_t nelem_sc;
qpms_l_t *s1_jMaxes; qpms_l_t *s1_jMaxes;
complex double **s1_constfacs; // indices [y][j] where j is same as in [1, (4.5)] complex double **s1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
// TODO probably normalisation and equatorial legendre polynomials should be included, too // TODO probably normalisation and equatorial legendre polynomials should be included, too
@ -88,8 +88,8 @@ int ewald32_sigma0(complex double *result, double *err,
// are not included. // are not included.
int ewald32_sigma_long_shiftedpoints_e ( int ewald32_sigma_long_shiftedpoints_e (
complex double *target_sigmalr_y, // must be c->nelem long complex double *target_sigmalr_y, // must be c->nelem_sc long
double *target_sigmalr_y_err, // must be c->nelem long or NULL double *target_sigmalr_y_err, // must be c->nelem_sc long or NULL
const qpms_ewald32_constants_t *c, const qpms_ewald32_constants_t *c,
double eta, double k, double unitcell_area, double eta, double k, double unitcell_area,
size_t npoints, const point2d *Kpoints_plus_beta, size_t npoints, const point2d *Kpoints_plus_beta,
@ -97,8 +97,8 @@ int ewald32_sigma_long_shiftedpoints_e (
point2d particle_shift point2d particle_shift
); );
int ewald32_sigma_long_points_and_shift (//NI int ewald32_sigma_long_points_and_shift (//NI
complex double *target_sigmalr_y, // must be c->nelem long complex double *target_sigmalr_y, // must be c->nelem_sc long
double *target_sigmalr_y_err, // must be c->nelem long or NULL double *target_sigmalr_y_err, // must be c->nelem_sc long or NULL
const qpms_ewald32_constants_t *c, const qpms_ewald32_constants_t *c,
double eta, double k, double unitcell_area, double eta, double k, double unitcell_area,
size_t npoints, const point2d *Kpoints, size_t npoints, const point2d *Kpoints,
@ -106,8 +106,8 @@ int ewald32_sigma_long_points_and_shift (//NI
point2d particle_shift point2d particle_shift
); );
int ewald32_sigma_long_shiftedpoints_rordered(//NI int ewald32_sigma_long_shiftedpoints_rordered(//NI
complex double *target_sigmalr_y, // must be c->nelem long complex double *target_sigmalr_y, // must be c->nelem_sc long
double *target_sigmalr_y_err, // must be c->nelem long or NULL double *target_sigmalr_y_err, // must be c->nelem_sc long or NULL
const qpms_ewald32_constants_t *c, const qpms_ewald32_constants_t *c,
double eta, double k, double unitcell_area, double eta, double k, double unitcell_area,
const points2d_rordered_t *Kpoints_plus_beta_rordered, const points2d_rordered_t *Kpoints_plus_beta_rordered,
@ -115,8 +115,8 @@ int ewald32_sigma_long_shiftedpoints_rordered(//NI
); );
int ewald32_sigma_short_shiftedpoints( int ewald32_sigma_short_shiftedpoints(
complex double *target_sigmasr_y, // must be c->nelem long complex double *target_sigmasr_y, // must be c->nelem_sc long
double *target_sigmasr_y_err, // must be c->nelem long or NULL double *target_sigmasr_y_err, // must be c->nelem_sc long or NULL
const qpms_ewald32_constants_t *c, // N.B. not too useful here const qpms_ewald32_constants_t *c, // N.B. not too useful here
double eta, double k, double eta, double k,
size_t npoints, const point2d *Rpoints_plus_particle_shift, size_t npoints, const point2d *Rpoints_plus_particle_shift,
@ -124,16 +124,16 @@ int ewald32_sigma_short_shiftedpoints(
point2d particle_shift // used only in the very end to multiply it by the phase point2d particle_shift // used only in the very end to multiply it by the phase
); );
int ewald32_sigma_short_points_and_shift(//NI int ewald32_sigma_short_points_and_shift(//NI
complex double *target_sigmasr_y, // must be c->nelem long complex double *target_sigmasr_y, // must be c->nelem_sc long
double *target_sigmasr_y_err, // must be c->nelem long or NULL double *target_sigmasr_y_err, // must be c->nelem_sc long or NULL
const qpms_ewald32_constants_t *c, // N.B. not too useful here const qpms_ewald32_constants_t *c, // N.B. not too useful here
double eta, double k, double eta, double k,
size_t npoints, const point2d *Rpoints, size_t npoints, const point2d *Rpoints,
point2d particle_shift point2d particle_shift
); );
int ewald32_sigma_short_points_rordered(//NI int ewald32_sigma_short_points_rordered(//NI
complex double *target_sigmasr_y, // must be c->nelem long complex double *target_sigmasr_y, // must be c->nelem_sc long
double *target_sigmasr_y_err, // must be c->nelem long or NULL double *target_sigmasr_y_err, // must be c->nelem_sc long or NULL
const qpms_ewald32_constants_t *c, // N.B. not too useful here const qpms_ewald32_constants_t *c, // N.B. not too useful here
double eta, double k, double eta, double k,
const points2d_rordered_t *Rpoints_plus_particle_shift_rordered, const points2d_rordered_t *Rpoints_plus_particle_shift_rordered,

View File

@ -25,4 +25,28 @@ static inline qpms_y_t qpms_lMax2nelem(qpms_l_t lmax){
return lmax * ((qpms_y_t)lmax + 2); return lmax * ((qpms_y_t)lmax + 2);
} }
// Scalar versions: they have a place for the 0, 0 term in the beginning
static inline qpms_y_t qpms_mn2y_sc(qpms_m_t m, qpms_l_t n) {
return n * (n + 1) + m;
}
static inline qpms_lm_t qpms_y2n_sc(qpms_y_t y) {
//return (sqrt(5+y)-2)/2; // the cast will truncate the fractional part, which is what we want
return sqrt(y);
}
static inline qpms_m_t qpms_yn2m_sc(qpms_y_t y, qpms_l_t n) {
return y-qpms_mn2y_sc(0,n);
}
static inline void qpms_y2mn_sc_p(qpms_y_t y, qpms_m_t *m, qpms_l_t *n){
*m=qpms_yn2m_sc(y,*n=qpms_y2n_sc(y));
}
static inline qpms_y_t qpms_lMax2nelem_sc(qpms_l_t lmax){
return lmax * ((qpms_y_t)lmax + 2) + 1;
}
#endif //QPMS_INDEXING_H #endif //QPMS_INDEXING_H

21
tests/ewalds.c Normal file
View File

@ -0,0 +1,21 @@
// implementation of the [LT(4.16)] test
#include <qpms/ewald.h>
typedef struct ewaldtest_hex_params {
qpms_l_t lMax;
poinnt2d beta;
double k;
double h;
double eta;
} ewaldtest_hex_params;
ewaldtest_hex_paraps paramslist = {
{ 3, {1.1, 0.23}, 2.3, 0.97, 0.3},
// end:
{ 0, {0, 0}, 0, 0, 0}
}