From 8bf9a1c54d57ea2b575c744653646f318d9b0898 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 14 May 2018 09:16:39 +0300 Subject: [PATCH] Fix besselbuf sizes, start implementing short-range parts. Former-commit-id: 7858c87377afe9e6484f6bd906d2fabfb9953945 --- qpms/translations.c | 155 ++++++++++++++++++++++++++++++++++++++++++-- qpms/translations.h | 37 +++++++++++ 2 files changed, 185 insertions(+), 7 deletions(-) diff --git a/qpms/translations.c b/qpms/translations.c index 4d741a5..8d1ff99 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -500,6 +500,9 @@ void qpms_trans_calculator_free(qpms_trans_calculator *c) { free(c->A_multipliers); free(c->B_multipliers[0]); free(c->B_multipliers); +#ifdef LATTICESUMS + free(c->hct); +#endif free(c); } @@ -908,6 +911,9 @@ qpms_trans_calculator } free(qmaxes); +#ifdef LATTICESUMS + c->hct = hankelcoefftable_init(2*lMax+1); +#endif return c; } @@ -1005,7 +1011,7 @@ complex double qpms_trans_calculator_get_B_buf(const qpms_trans_calculator *c, double costheta = cos(kdlj.theta); if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1, costheta,csphase,legendre_buf)) abort(); - if (qpms_sph_bessel_fill(J, n+nu+2, kdlj.r, bessel_buf)) abort(); + if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort(); return qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu, kdlj,r_ge_d,J,bessel_buf,legendre_buf); } @@ -1039,7 +1045,7 @@ int qpms_trans_calculator_get_AB_buf_p(const qpms_trans_calculator *c, double costheta = cos(kdlj.theta); if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1, costheta,-1,legendre_buf)) abort(); - if (qpms_sph_bessel_fill(J, n+nu+2, kdlj.r, bessel_buf)) abort(); + if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort(); *Adest = qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu, kdlj,r_ge_d,J,bessel_buf,legendre_buf); *Bdest = qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu, @@ -1077,7 +1083,7 @@ int qpms_trans_calculator_get_AB_arrays_buf(const qpms_trans_calculator *c, double costheta = cos(kdlj.theta); if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1, costheta,-1,legendre_buf)) abort(); - if (qpms_sph_bessel_fill(J, 2*c->lMax+2, kdlj.r, bessel_buf)) abort(); + if (qpms_sph_bessel_fill(J, 2*c->lMax+1, kdlj.r, bessel_buf)) abort(); size_t desti = 0, srci = 0; 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) { @@ -1107,7 +1113,7 @@ complex double qpms_trans_calculator_get_A(const qpms_trans_calculator *c, int m, int n, int mu, int nu, sph_t kdlj, bool r_ge_d, qpms_bessel_t J) { double leg[gsl_sf_legendre_array_n(n+nu)]; - complex double bes[n+nu+1]; + complex double bes[n+nu+1]; // maximum order is 2n for A coeffs, plus the zeroth. return qpms_trans_calculator_get_A_buf(c,m,n,mu,nu,kdlj,r_ge_d,J, bes,leg); } @@ -1116,7 +1122,7 @@ complex double qpms_trans_calculator_get_B(const qpms_trans_calculator *c, int m, int n, int mu, int nu, sph_t kdlj, bool r_ge_d, qpms_bessel_t J) { double leg[gsl_sf_legendre_array_n(n+nu+1)]; - complex double bes[n+nu+2]; + complex double bes[n+nu+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth. return qpms_trans_calculator_get_B_buf(c,m,n,mu,nu,kdlj,r_ge_d,J, bes,leg); } @@ -1126,7 +1132,7 @@ int qpms_trans_calculator_get_AB_p(const qpms_trans_calculator *c, int m, int n, int mu, int nu, sph_t kdlj, bool r_ge_d, qpms_bessel_t J) { double leg[gsl_sf_legendre_array_n(2*c->lMax+1)]; - complex double bes[2*c->lMax+3]; // TODO check lMax + complex double bes[2*c->lMax+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth. return qpms_trans_calculator_get_AB_buf_p(c,Adest, Bdest,m,n,mu,nu,kdlj,r_ge_d,J, bes,leg); } @@ -1136,7 +1142,7 @@ int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c, size_t deststride, size_t srcstride, sph_t kdlj, bool r_ge_d, qpms_bessel_t J) { double leg[gsl_sf_legendre_array_n(c->lMax+c->lMax+1)]; - complex double bes[2*c->lMax+3]; // TODO check lMax + complex double bes[2*c->lMax+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth. return qpms_trans_calculator_get_AB_arrays_buf(c, Adest, Bdest, deststride, srcstride, kdlj, r_ge_d, J, @@ -1144,6 +1150,141 @@ int qpms_trans_calculator_get_AB_arrays(const qpms_trans_calculator *c, } +#ifdef LATTICESUMS +int qpms_trans_calculator_get_shortrange_AB_arrays_buf(const qpms_trans_calculator *c, + complex double *Adest, complex double *Bdest, + size_t deststride, size_t srcstride, + sph_t kdlj, qpms_bessel_t J, + qpms_l_t lrcutoff, unsigned kappa, double c, // regularisation params + complex double *bessel_buf, double *legendre_buf, + ) { + assert(J == QPMS_HANKEL_PLUS); // support only J == 3 for now + if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) { + for (size_t i = 0; i < c->nelem; ++i) + for (size_t j = 0; j < c->nelem; ++j) { + *(Adest + i*srcstride + j*deststride) = NAN+I*NAN; + *(Bdest + i*srcstride + j*deststride) = NAN+I*NAN; + } + // TODO warn? different return value? + return 0; + } + switch(qpms_normalisation_t_normonly(c->normalisation)) { + case QPMS_NORMALISATION_TAYLOR: + case QPMS_NORMALISATION_POWER: + //case QPMS_NORMALISATION_NONE: // I am not sure the Hankel transform work the same way for unnormalised waves, so disallow for now + { + double costheta = cos(kdlj.theta); + if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1, + costheta,-1,legendre_buf)) abort(); + // if (qpms_sph_bessel_fill(J, 2*c->lMax+1, kdlj.r, bessel_buf)) abort(); // original + hankelparts_fill(NULL, bessel_buf, 2*c->lMax+1, lrcutoff, c->hct, kappa, c, kdlj.r); + size_t desti = 0, srci = 0; + 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) { + size_t assertindex = qpms_trans_calculator_index_mnmunu(c,m,n,mu,nu); + assert(assertindex == desti*c->nelem + srci); + *(Adest + deststride * desti + srcstride * srci) = + qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu, + kdlj,false,J,bessel_buf,legendre_buf); + *(Bdest + deststride * desti + srcstride * srci) = + qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu, + kdlj,false,J,bessel_buf,legendre_buf); + ++srci; + } + ++desti; + srci = 0; + } + return 0; + } + break; + default: + abort(); + } + assert(0); +} + +int qpms_trans_calculator_get_shortrange_AB_buf_p(const qpms_trans_calculator *c, + complex double *Adest, complex double *Bdest, + int m, int n, int mu, int nu, sph_t kdlj, + qpms_bessel_t J, + qpms_l_t lrcutoff, unsigned kappa, double c, // regularisation params + complex double *bessel_buf, double *legendre_buf) { + assert(J == QPMS_HANKEL_PLUS); // support only J == 3 for now + if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) { + *Adest = NAN+I*NAN; + *Bdest = NAN+I*NAN; + // TODO warn? different return value? + return 0; + } + switch(qpms_normalisation_t_normonly(c->normalisation)) { + case QPMS_NORMALISATION_TAYLOR: + case QPMS_NORMALISATION_KRISTENSSON: + // case QPMS_NORMALISATION_NONE: // Not sure if it would work, so disable for now + { + double costheta = cos(kdlj.theta); + if (gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,n+nu+1, + costheta,-1,legendre_buf)) abort(); + //if (qpms_sph_bessel_fill(J, n+nu+1, kdlj.r, bessel_buf)) abort(); // original + hankelparts_fill(NULL, bessel_buf, 2*c->lMax+1, lrcutoff, c->hct, kappa, c, kdlj.r); + + *Adest = qpms_trans_calculator_get_A_precalcbuf(c,m,n,mu,nu, + kdlj,r_ge_d,J,bessel_buf,legendre_buf); + *Bdest = qpms_trans_calculator_get_B_precalcbuf(c,m,n,mu,nu, + kdlj,r_ge_d,J,bessel_buf,legendre_buf); + return 0; + } + break; + default: + abort(); + } + assert(0); +} + +// Short-range parts of the translation coefficients +int qpms_trans_calculator_get_shortrange_AB_p(const qpms_trans_calculator *c, + complex double *Adest, complex double *Bdest, + qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, + qpms_bessel_t J /* Only J=3 valid for now */, + qpms_l_t lrcutoff, unsigned kappa, double c) { + double leg[gsl_sf_legendre_array_n(2*c->lMax+1)]; + complex double bes[2*c->lMax+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth. + return qpms_trans_calculator_get_shortrange_AB_buf_p(c,Adest, Bdest,m,n,mu,nu,kdlj,J, + lrcutoff, kappa, c, + bes, leg); +} + +int qpms_trans_calculator_get_shortrange_AB_arrays(const qpms_trans_calculator *c, + complex double *Adest, complex double *Bdest, + size_t deststride, size_t srcstride, + sph_t kdlj, qpms_bessel_t J /* Only J=3 valid for now */, + qpms_l_t lrcutoff, unsigned kappa, double c) { + double leg[gsl_sf_legendre_array_n(c->lMax+c->lMax+1)]; + complex double bes[2*c->lMax+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth. + return qpms_trans_calculator_get_AB_arrays_buf(c, + Adest, Bdest, deststride, srcstride, + kdlj, J, + lrcutoff, kappa, c, + bes, leg); +} + +// Fourier transforms of the long-range parts of the translation coefficients +int qpms_trans_calculator_get_Fourier_longrange_AB_p(const qpms_trans_calculator *c, + complex double *Adest, complex double *Bdest, + qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t k_sph, + qpms_bessel_t J /* Only J=3 valid for now */, + qpms_l_t lrcutoff, unsigned kappa, double cv, double k0) { + TODO; +} + +int qpms_trans_calculator_get_Fourier_longrange_AB_arrays(const qpms_trans_calculator *c, + complex double *Adest, complex double *Bdest, + size_t deststride, size_t srcstride, + sph_t k_sph, qpms_bessel_t J /* Only J=3 valid for now */, + qpms_l_t lrcutoff, unsigned kappa, double cv, double k0) { + TODO; +} +#endif // LATTICESUMS + complex double qpms_trans_calculator_get_A_ext(const qpms_trans_calculator *c, int m, int n, int mu, int nu, diff --git a/qpms/translations.h b/qpms/translations.h index 2b8a772..795b464 100644 --- a/qpms/translations.h +++ b/qpms/translations.h @@ -6,6 +6,10 @@ #include #include +#ifdef LATTICESUMS +#include "bessels.h" +#endif + // TODO replace the xplicit "Taylor" functions with general, // taking qpms_normalisation_t argument. complex double qpms_trans_single_A_Taylor(qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, @@ -41,8 +45,12 @@ typedef struct qpms_trans_calculator { // Spherical Bessel function coefficients: // TODO #endif +#ifdef LATTICESUMS + complex double *hct; // Hankel function coefficient table +#endif } qpms_trans_calculator; + qpms_trans_calculator *qpms_trans_calculator_init(qpms_l_t lMax, qpms_normalisation_t nt); void qpms_trans_calculator_free(qpms_trans_calculator *); @@ -82,6 +90,35 @@ int qpms_trans_calculator_get_AB_arrays_ext(const qpms_trans_calculator *c, double kdlj_r, double kdlj_theta, double kdlj_phi, int r_ge_d, int J); +#ifdef LATTICESUMS +// Short-range parts of the translation coefficients +int qpms_trans_calculator_get_shortrange_AB_p(const qpms_trans_calculator *c, + complex double *Adest, complex double *Bdest, + qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t kdlj, + qpms_bessel_t J /* Only J=3 valid for now */, + qpms_l_t longrange_order_cutoff, unsigned kappa, double c); +int qpms_trans_calculator_get_shortrange_AB_arrays(const qpms_trans_calculator *c, + complex double *Adest, complex double *Bdest, + size_t deststride, size_t srcstride, + sph_t kdlj, qpms_bessel_t J /* Only J=3 valid for now */, + qpms_l_t longrange_order_cutoff, unsigned kappa, double c); + +// Fourier transforms of the long-range parts of the translation coefficients +int qpms_trans_calculator_get_Fourier_longrange_AB_p(const qpms_trans_calculator *c, + complex double *Adest, complex double *Bdest, + qpms_m_t m, qpms_l_t n, qpms_m_t mu, qpms_l_t nu, sph_t k_sph, + qpms_bessel_t J /* Only J=3 valid for now */, + qpms_l_t longrange_order_cutoff, unsigned kappa, double cv, double k0); + +int qpms_trans_calculator_get_Fourier_longrange_AB_arrays(const qpms_trans_calculator *c, + complex double *Adest, complex double *Bdest, + size_t deststride, size_t srcstride, + sph_t k_sph, qpms_bessel_t J /* Only J=3 valid for now */, + qpms_l_t longrange_order_cutoff, unsigned kappa, double cv, double k0); +#endif // LATTICESUMS + + + #ifdef QPMS_COMPILE_PYTHON_EXTENSIONS #include #include