From 5b5e77f98521568f7830f71e61a3ef57f5327875 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 13 Mar 2018 10:57:08 +0200 Subject: [PATCH] Kahan sums 3 Former-commit-id: e34f5bbb16a95c4d5517121454a4d5229bcb50c8 --- qpms/vswf.c | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/qpms/vswf.c b/qpms/vswf.c index aa03ab9..75e9d6a 100644 --- a/qpms/vswf.c +++ b/qpms/vswf.c @@ -296,21 +296,27 @@ csphvec_t qpms_eval_vswf(sph_t kr, qpms_l_t lMax, qpms_bessel_t btyp, qpms_normalisation_t norm) { qpms_y_t nelem = qpms_lMax2nelem(lMax); - csphvec_t lsum = {0, 0, 0}, msum = {0, 0, 0}, esum = {0, 0, 0}; + csphvec_t lsum, msum, esum, lcomp, mcomp, ecomp; + csphvec_kahaninit(&lsum, &lcomp); + csphvec_kahaninit(&msum, &mcomp); + csphvec_kahaninit(&esum, &ecomp); csphvec_t *lset = NULL, *mset = NULL, *eset = NULL; if(lc) lset = malloc(nelem * sizeof(csphvec_t)); if(mc) mset = malloc(nelem * sizeof(csphvec_t)); if(ec) eset = malloc(nelem * sizeof(csphvec_t)); qpms_vswf_fill(lset, mset, eset, lMax, kr, btyp, norm); if(lc) for(qpms_y_t y = 0; y < nelem; ++y) - lsum = csphvec_add(lsum, csphvec_scale(lc[y], lset[y])); + csphvec_kahanadd(&lsum, &lcomp, csphvec_scale(lc[y], lset[y])); if(mc) for(qpms_y_t y = 0; y < nelem; ++y) - msum = csphvec_add(msum, csphvec_scale(mc[y], mset[y])); + csphvec_kahanadd(&msum, &mcomp, csphvec_scale(mc[y], mset[y])); if(ec) for(qpms_y_t y = 0; y < nelem; ++y) - esum = csphvec_add(esum, csphvec_scale(ec[y], eset[y])); + csphvec_kahanadd(&esum, &ecomp, csphvec_scale(ec[y], eset[y])); if(lc) free(lset); if(mc) free(mset); if(ec) free(eset); - return csphvec_add(esum, csphvec_add(msum, lsum)); + //return csphvec_add(esum, csphvec_add(msum, lsum)); + csphvec_kahanadd(&esum, &ecomp, msum); + csphvec_kahanadd(&esum, &ecomp, lsum); + return esum; }