From 7b2260d31d65b448d177cb4086d4043eae98ad40 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sun, 4 Jun 2017 14:36:15 -0400
Subject: [PATCH 1/9] Steps towards generalizing the fmmlib Helmholtz expansion
 wrangler to 3D

---
 boxtree/pyfmmlib_integration.py | 74 ++++++++++++++++++++++++---------
 test/test_fmm.py                | 27 ++++++------
 2 files changed, 70 insertions(+), 31 deletions(-)

diff --git a/boxtree/pyfmmlib_integration.py b/boxtree/pyfmmlib_integration.py
index 34126f2..87eb0f6 100644
--- a/boxtree/pyfmmlib_integration.py
+++ b/boxtree/pyfmmlib_integration.py
@@ -33,7 +33,7 @@ __doc__ = """Integrates :mod:`boxtree` with
 """
 
 
-class Helmholtz2DExpansionWrangler:
+class HelmholtzExpansionWrangler:
     """Implements the :class:`boxtree.fmm.ExpansionWranglerInterface`
     by using pyfmmlib.
     """
@@ -43,8 +43,44 @@ class Helmholtz2DExpansionWrangler:
         self.helmholtz_k = helmholtz_k
         self.nterms = nterms
 
+        common_extra_kwargs = {}
+        if tree.dimensions == 3:
+            nquad = max(6, int(2.5*nterms))
+            import scipy.special as sps
+            weights = sps.legendre(nquad).weights
+
+            common_extra_kwargs = {
+                    "xnodes": weights[:, 0],
+                    "wts": weights[:, 2],
+                    }
+
+        self.common_extra_kwargs = common_extra_kwargs
+
+    def get_routine(self, name, vec=False):
+        suffix = ""
+        if self.tree.dimensions == 3:
+            suffix = "quadu"
+
+        if vec:
+            suffix += "_vec"
+
+        import pyfmmlib
+        return getattr(pyfmmlib, "h%s%s" % (
+            name % self.tree.dimensions,
+            suffix))
+
+    def get_vec_routine(self, name):
+        return self.get_routine(name, vec=True)
+
     def multipole_expansion_zeros(self):
-        return np.zeros((self.tree.nboxes, 2*self.nterms+1), dtype=np.complex128)
+        if self.tree.dimensions == 2:
+            return np.zeros((self.tree.nboxes, 2*self.nterms+1), dtype=np.complex128)
+        elif self.tree.dimensions == 3:
+            return np.zeros(
+                    (self.tree.nboxes, self.nterms+1, 2*self.nterms+1),
+                    dtype=np.complex128)
+        else:
+            raise ValueError("unsupported functionality")
 
     local_expansion_zeros = multipole_expansion_zeros
 
@@ -84,7 +120,7 @@ class Helmholtz2DExpansionWrangler:
     def form_multipoles(self, level_start_source_box_nrs, source_boxes, src_weights):
         rscale = 1  # FIXME
 
-        from pyfmmlib import h2dformmp
+        formmp = self.get_routine("%ddformmp")
 
         mpoles = self.multipole_expansion_zeros()
         for src_ibox in source_boxes:
@@ -93,7 +129,7 @@ class Helmholtz2DExpansionWrangler:
             if pslice.stop - pslice.start == 0:
                 continue
 
-            ier, mpoles[src_ibox] = h2dformmp(
+            ier, mpoles[src_ibox] = formmp(
                     self.helmholtz_k, rscale, self._get_sources(pslice),
                     src_weights[pslice],
                     self.tree.box_centers[:, src_ibox], self.nterms)
@@ -107,7 +143,7 @@ class Helmholtz2DExpansionWrangler:
         tree = self.tree
         rscale = 1  # FIXME
 
-        from pyfmmlib import h2dmpmp_vec
+        mpmp = self.get_vec_routine("%ddmpmp")
 
         # 2 is the last relevant source_level.
         # 1 is the last relevant target_level.
@@ -121,7 +157,7 @@ class Helmholtz2DExpansionWrangler:
                     if child:
                         child_center = tree.box_centers[:, child]
 
-                        new_mp = h2dmpmp_vec(
+                        new_mp = mpmp(
                                 self.helmholtz_k,
                                 rscale, child_center, mpoles[child],
                                 rscale, parent_center, self.nterms)
@@ -132,7 +168,7 @@ class Helmholtz2DExpansionWrangler:
             neighbor_sources_lists, src_weights):
         pot = self.potential_zeros()
 
-        from pyfmmlib import hpotgrad2dall_vec
+        potgradall = self.get_vec_routine("potgrad%ddall")
 
         for itgt_box, tgt_ibox in enumerate(target_boxes):
             tgt_pslice = self._get_target_slice(tgt_ibox)
@@ -148,7 +184,7 @@ class Helmholtz2DExpansionWrangler:
                 if src_pslice.stop - src_pslice.start == 0:
                     continue
 
-                tmp_pot, _, _ = hpotgrad2dall_vec(
+                tmp_pot, _, _ = potgradall(
                         ifgrad=False, ifhess=False,
                         sources=self._get_sources(src_pslice),
                         charge=src_weights[src_pslice],
@@ -169,7 +205,7 @@ class Helmholtz2DExpansionWrangler:
 
         rscale = 1
 
-        from pyfmmlib import h2dmploc_vec
+        mploc = self.get_vec_routine("%ddmploc")
 
         for itgt_box, tgt_ibox in enumerate(target_or_target_parent_boxes):
             start, end = starts[itgt_box:itgt_box+2]
@@ -181,7 +217,7 @@ class Helmholtz2DExpansionWrangler:
             for src_ibox in lists[start:end]:
                 src_center = tree.box_centers[:, src_ibox]
 
-                tgt_loc = tgt_loc + h2dmploc_vec(
+                tgt_loc = tgt_loc + mploc(
                         self.helmholtz_k,
                         rscale, src_center, mpole_exps[src_ibox],
                         rscale, tgt_center, self.nterms)[:, 0]
@@ -196,7 +232,7 @@ class Helmholtz2DExpansionWrangler:
 
         rscale = 1
 
-        from pyfmmlib import h2dmpeval_vec
+        mpeval = self.get_vec_routine("%ddmpeval")
 
         for ssn in sep_smaller_nonsiblings_by_level:
             for itgt_box, tgt_ibox in enumerate(target_boxes):
@@ -209,7 +245,7 @@ class Helmholtz2DExpansionWrangler:
                 start, end = ssn.starts[itgt_box:itgt_box+2]
                 for src_ibox in ssn.lists[start:end]:
 
-                    tmp_pot, _, _ = h2dmpeval_vec(self.helmholtz_k, rscale, self.
+                    tmp_pot, _, _ = mpeval(self.helmholtz_k, rscale, self.
                             tree.box_centers[:, src_ibox], mpole_exps[src_ibox],
                             self._get_targets(tgt_pslice),
                             ifgrad=False, ifhess=False)
@@ -226,7 +262,7 @@ class Helmholtz2DExpansionWrangler:
         rscale = 1  # FIXME
         local_exps = self.local_expansion_zeros()
 
-        from pyfmmlib import h2dformta
+        formta = self.get_routine("%ddformta")
 
         for itgt_box, tgt_ibox in enumerate(target_or_target_parent_boxes):
             start, end = starts[itgt_box:itgt_box+2]
@@ -240,12 +276,12 @@ class Helmholtz2DExpansionWrangler:
                 if src_pslice.stop - src_pslice.start == 0:
                     continue
 
-                ier, mpole = h2dformta(
+                ier, mpole = formta(
                         self.helmholtz_k, rscale,
                         self._get_sources(src_pslice), src_weights[src_pslice],
                         tgt_center, self.nterms)
                 if ier:
-                    raise RuntimeError("h2dformta failed")
+                    raise RuntimeError("formta failed")
 
                 contrib = contrib + mpole
 
@@ -257,7 +293,7 @@ class Helmholtz2DExpansionWrangler:
             target_or_target_parent_boxes, local_exps):
         rscale = 1  # FIXME
 
-        from pyfmmlib import h2dlocloc_vec
+        locloc = self.get_vec_routine("%ddlocloc")
 
         for target_lev in range(1, self.tree.nlevels):
             start, stop = level_start_target_or_target_parent_box_nrs[
@@ -268,7 +304,7 @@ class Helmholtz2DExpansionWrangler:
                 src_ibox = self.tree.box_parent_ids[tgt_ibox]
                 src_center = self.tree.box_centers[:, src_ibox]
 
-                tmp_loc_exp = h2dlocloc_vec(
+                tmp_loc_exp = locloc(
                             self.helmholtz_k,
                             rscale, src_center, local_exps[src_ibox],
                             rscale, tgt_center, self.nterms)[:, 0]
@@ -281,7 +317,7 @@ class Helmholtz2DExpansionWrangler:
         pot = self.potential_zeros()
         rscale = 1  # FIXME
 
-        from pyfmmlib import h2dtaeval_vec
+        taeval = self.get_vec_routine("%ddtaeval")
 
         for tgt_ibox in target_boxes:
             tgt_pslice = self._get_target_slice(tgt_ibox)
@@ -289,7 +325,7 @@ class Helmholtz2DExpansionWrangler:
             if tgt_pslice.stop - tgt_pslice.start == 0:
                 continue
 
-            tmp_pot, _, _ = h2dtaeval_vec(self.helmholtz_k, rscale,
+            tmp_pot, _, _ = taeval(self.helmholtz_k, rscale,
                     self.tree.box_centers[:, tgt_ibox], local_exps[tgt_ibox],
                     self._get_targets(tgt_pslice), ifgrad=False, ifhess=False)
 
diff --git a/test/test_fmm.py b/test/test_fmm.py
index 4fedf29..1bf3d84 100644
--- a/test/test_fmm.py
+++ b/test/test_fmm.py
@@ -454,7 +454,8 @@ def test_fmm_completeness(ctx_getter, dims, nsources_req, ntargets_req,
 
 # {{{ test Helmholtz fmm with pyfmmlib
 
-def test_pyfmmlib_fmm(ctx_getter):
+@pytest.mark.parametrize("dims", [2])
+def test_pyfmmlib_fmm(ctx_getter, dims):
     logging.basicConfig(level=logging.INFO)
 
     from pytest import importorskip
@@ -465,7 +466,6 @@ def test_pyfmmlib_fmm(ctx_getter):
 
     nsources = 3000
     ntargets = 1000
-    dims = 2
     dtype = np.float64
 
     helmholtz_k = 2
@@ -473,7 +473,7 @@ def test_pyfmmlib_fmm(ctx_getter):
     sources = p_normal(queue, nsources, dims, dtype, seed=15)
     targets = (
             p_normal(queue, ntargets, dims, dtype, seed=18)
-            + np.array([2, 0]))
+            + np.array([2, 0, 0])[:dims])
 
     sources_host = particle_array_to_host(sources)
     targets_host = particle_array_to_host(targets)
@@ -496,19 +496,22 @@ def test_pyfmmlib_fmm(ctx_getter):
     weights = rng.uniform(queue, nsources, dtype=np.float64).get()
     #weights = np.ones(nsources)
 
-    logger.info("computing direct (reference) result")
-
-    from pyfmmlib import hpotgrad2dall_vec
-    ref_pot, _, _ = hpotgrad2dall_vec(ifgrad=False, ifhess=False,
-            sources=sources_host.T, charge=weights,
-            targets=targets_host.T, zk=helmholtz_k)
-
-    from boxtree.pyfmmlib_integration import Helmholtz2DExpansionWrangler
-    wrangler = Helmholtz2DExpansionWrangler(trav.tree, helmholtz_k, nterms=10)
+    from boxtree.pyfmmlib_integration import HelmholtzExpansionWrangler
+    wrangler = HelmholtzExpansionWrangler(trav.tree, helmholtz_k, nterms=10)
 
     from boxtree.fmm import drive_fmm
     pot = drive_fmm(trav, wrangler, weights)
 
+    logger.info("computing direct (reference) result")
+
+    if dims == 2:
+        from pyfmmlib import hpotgrad2dall_vec
+        ref_pot, _, _ = hpotgrad2dall_vec(ifgrad=False, ifhess=False,
+                sources=sources_host.T, charge=weights,
+                targets=targets_host.T, zk=helmholtz_k)
+    else:
+        raise NotImplementedError()
+
     rel_err = la.norm(pot - ref_pot) / la.norm(ref_pot)
     logger.info("relative l2 error: %g" % rel_err)
     assert rel_err < 1e-5
-- 
GitLab


From 7fcf506f108442ea606ef2d2622ec2056d52fb22 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sun, 4 Jun 2017 18:07:34 -0400
Subject: [PATCH 2/9] More steps towards Helmholtz 3D generalization

---
 boxtree/pyfmmlib_integration.py | 92 ++++++++++++++++++++++++++-------
 1 file changed, 73 insertions(+), 19 deletions(-)

diff --git a/boxtree/pyfmmlib_integration.py b/boxtree/pyfmmlib_integration.py
index 87eb0f6..6b6b1fd 100644
--- a/boxtree/pyfmmlib_integration.py
+++ b/boxtree/pyfmmlib_integration.py
@@ -26,6 +26,7 @@ THE SOFTWARE.
 
 
 import numpy as np
+from functools import partial
 
 
 __doc__ = """Integrates :mod:`boxtree` with
@@ -56,21 +57,64 @@ class HelmholtzExpansionWrangler:
 
         self.common_extra_kwargs = common_extra_kwargs
 
-    def get_routine(self, name, vec=False):
-        suffix = ""
-        if self.tree.dimensions == 3:
-            suffix = "quadu"
-
-        if vec:
-            suffix += "_vec"
-
+    def get_routine(self, name, suffix=""):
         import pyfmmlib
         return getattr(pyfmmlib, "h%s%s" % (
             name % self.tree.dimensions,
             suffix))
 
     def get_vec_routine(self, name):
-        return self.get_routine(name, vec=True)
+        return self.get_routine(name, "_vec")
+
+    def get_translation_routine(self, name):
+        suffix = ""
+        if self.tree.dimensions == 3:
+            suffix = "quadu"
+        suffix += "_vec"
+
+        rout = self.get_routine(name, suffix)
+
+        if self.tree.dimensions == 2:
+            return rout
+        else:
+
+            def wrapper(*args, **kwargs):
+                kwargs.update(self.common_extra_kwargs)
+                val, ier = rout(*args, **kwargs)
+                if (ier != 0).any():
+                    raise RuntimeError("%s failed with nonzero ier" % name)
+
+                return val
+
+            from functools import update_wrapper
+            update_wrapper(wrapper, rout)
+            return wrapper
+
+    def get_direct_eval_routine(self):
+        if self.tree.dimensions == 2:
+            rout = self.get_routine("potgrad%ddall", "_vec")
+
+            def wrapper(*args, **kwargs):
+                pot, grad, hess = rout(*args, **kwargs, ifgrad=False, ifhess=False)
+                return pot
+
+            from functools import update_wrapper
+            update_wrapper(wrapper, rout)
+            return wrapper
+
+        elif self.tree.dimensions == 3:
+            rout = self.get_routine("potfld%ddall", "_vec")
+
+            def wrapper(*args, **kwargs):
+                pot, fld = rout(*args, **kwargs, iffld=False)
+                # grad = -fld
+                return pot
+
+            from functools import update_wrapper
+            update_wrapper(wrapper, rout)
+            return wrapper
+        else:
+            raise ValueError("unsupported dimensionality")
 
     def multipole_expansion_zeros(self):
         if self.tree.dimensions == 2:
@@ -143,7 +187,7 @@ class HelmholtzExpansionWrangler:
         tree = self.tree
         rscale = 1  # FIXME
 
-        mpmp = self.get_vec_routine("%ddmpmp")
+        mpmp = self.get_translation_routine("%ddmpmp")
 
         # 2 is the last relevant source_level.
         # 1 is the last relevant target_level.
@@ -151,24 +195,31 @@ class HelmholtzExpansionWrangler:
         for source_level in range(tree.nlevels-1, 1, -1):
             start, stop = level_start_source_parent_box_nrs[
                             source_level:source_level+2]
+            target_level = source_level - 1
+
             for ibox in source_parent_boxes[start:stop]:
                 parent_center = tree.box_centers[:, ibox]
                 for child in tree.box_child_ids[:, ibox]:
                     if child:
                         child_center = tree.box_centers[:, child]
 
+                        kwargs = {}
+                        if self.tree.dimensions == 3:
+                            kwargs["radius"] = tree.root_extent * 2**(-target_level)
+
                         new_mp = mpmp(
                                 self.helmholtz_k,
                                 rscale, child_center, mpoles[child],
-                                rscale, parent_center, self.nterms)
+                                rscale, parent_center, self.nterms,
+                                **kwargs)
 
-                        mpoles[ibox] += new_mp[:, 0]
+                        mpoles[ibox] += new_mp[..., 0]
 
     def eval_direct(self, target_boxes, neighbor_sources_starts,
             neighbor_sources_lists, src_weights):
         pot = self.potential_zeros()
 
-        potgradall = self.get_vec_routine("potgrad%ddall")
+        ev = self.get_direct_eval_routine()
 
         for itgt_box, tgt_ibox in enumerate(target_boxes):
             tgt_pslice = self._get_target_slice(tgt_ibox)
@@ -184,8 +235,7 @@ class HelmholtzExpansionWrangler:
                 if src_pslice.stop - src_pslice.start == 0:
                     continue
 
-                tmp_pot, _, _ = potgradall(
-                        ifgrad=False, ifhess=False,
+                tmp_pot = ev(
                         sources=self._get_sources(src_pslice),
                         charge=src_weights[src_pslice],
                         targets=self._get_targets(tgt_pslice), zk=self.helmholtz_k)
@@ -205,7 +255,7 @@ class HelmholtzExpansionWrangler:
 
         rscale = 1
 
-        mploc = self.get_vec_routine("%ddmploc")
+        mploc = self.get_translation_routine("%ddmploc")
 
         for itgt_box, tgt_ibox in enumerate(target_or_target_parent_boxes):
             start, end = starts[itgt_box:itgt_box+2]
@@ -217,10 +267,14 @@ class HelmholtzExpansionWrangler:
             for src_ibox in lists[start:end]:
                 src_center = tree.box_centers[:, src_ibox]
 
+                kwargs = {}
+                if self.tree.dimensions == 3:
+                    kwargs["radius"] = tree.root_extent * 2**(-target_level)
+
                 tgt_loc = tgt_loc + mploc(
                         self.helmholtz_k,
                         rscale, src_center, mpole_exps[src_ibox],
-                        rscale, tgt_center, self.nterms)[:, 0]
+                        rscale, tgt_center, self.nterms, **kwargs)[..., 0]
 
             local_exps[tgt_ibox] += tgt_loc
 
@@ -293,7 +347,7 @@ class HelmholtzExpansionWrangler:
             target_or_target_parent_boxes, local_exps):
         rscale = 1  # FIXME
 
-        locloc = self.get_vec_routine("%ddlocloc")
+        locloc = self.get_translation_routine("%ddlocloc")
 
         for target_lev in range(1, self.tree.nlevels):
             start, stop = level_start_target_or_target_parent_box_nrs[
@@ -307,7 +361,7 @@ class HelmholtzExpansionWrangler:
                 tmp_loc_exp = locloc(
                             self.helmholtz_k,
                             rscale, src_center, local_exps[src_ibox],
-                            rscale, tgt_center, self.nterms)[:, 0]
+                            rscale, tgt_center, self.nterms)[..., 0]
 
                 local_exps[tgt_ibox] += tmp_loc_exp
 
-- 
GitLab


From 1c7e8e3a437c99d6354c6ea1bb5dd09ce0b1bd10 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sun, 4 Jun 2017 18:19:48 -0400
Subject: [PATCH 3/9] Helmholtz 3D pyfmmlib: Per-level M2Ls

---
 boxtree/pyfmmlib_integration.py | 36 +++++++++++++++++++--------------
 1 file changed, 21 insertions(+), 15 deletions(-)

diff --git a/boxtree/pyfmmlib_integration.py b/boxtree/pyfmmlib_integration.py
index 6b6b1fd..3adccf3 100644
--- a/boxtree/pyfmmlib_integration.py
+++ b/boxtree/pyfmmlib_integration.py
@@ -257,26 +257,32 @@ class HelmholtzExpansionWrangler:
 
         mploc = self.get_translation_routine("%ddmploc")
 
-        for itgt_box, tgt_ibox in enumerate(target_or_target_parent_boxes):
-            start, end = starts[itgt_box:itgt_box+2]
-            tgt_center = tree.box_centers[:, tgt_ibox]
+        for lev in range(self.tree.nlevels):
+            lstart, lstop = level_start_target_or_target_parent_box_nrs[lev:lev+2]
+            if lstart == lstop:
+                continue
 
-            #print tgt_ibox, "<-", lists[start:end]
-            tgt_loc = 0
+            for itgt_box, tgt_ibox in enumerate(
+                    target_or_target_parent_boxes[lstart:lstop]):
+                start, end = starts[lstart + itgt_box:lstart + itgt_box+2]
+                tgt_center = tree.box_centers[:, tgt_ibox]
 
-            for src_ibox in lists[start:end]:
-                src_center = tree.box_centers[:, src_ibox]
+                #print tgt_ibox, "<-", lists[start:end]
+                tgt_loc = 0
+
+                for src_ibox in lists[start:end]:
+                    src_center = tree.box_centers[:, src_ibox]
 
-                kwargs = {}
-                if self.tree.dimensions == 3:
-                    kwargs["radius"] = tree.root_extent * 2**(-target_level)
+                    kwargs = {}
+                    if self.tree.dimensions == 3:
+                        kwargs["radius"] = tree.root_extent * 2**(-lev)
 
-                tgt_loc = tgt_loc + mploc(
-                        self.helmholtz_k,
-                        rscale, src_center, mpole_exps[src_ibox],
-                        rscale, tgt_center, self.nterms, **kwargs)[..., 0]
+                    tgt_loc = tgt_loc + mploc(
+                            self.helmholtz_k,
+                            rscale, src_center, mpole_exps[src_ibox],
+                            rscale, tgt_center, self.nterms, **kwargs)[..., 0]
 
-            local_exps[tgt_ibox] += tgt_loc
+                local_exps[tgt_ibox] += tgt_loc
 
         return local_exps
 
-- 
GitLab


From d0347df49799b8346c3b927fc025b44393bc6f41 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sun, 4 Jun 2017 20:05:29 -0400
Subject: [PATCH 4/9] Use https (rather than the git protocol) to fetch
 requirements

---
 requirements.txt | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/requirements.txt b/requirements.txt
index e36a43a..2a67d20 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,5 @@
 numpy
-git+git://github.com/pyopencl/pyopencl
-git+git://github.com/inducer/islpy
-git+git://github.com/inducer/loopy
-git+git://github.com/inducer/pyfmmlib
+https+git://github.com/pyopencl/pyopencl
+https+git://github.com/inducer/islpy
+https+git://github.com/inducer/loopy
+https+git://github.com/inducer/pyfmmlib
-- 
GitLab


From 12166bb535339777e51c21504d95f8afa5ee21ac Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sun, 4 Jun 2017 20:06:53 -0400
Subject: [PATCH 5/9] Apparently working h3d FMM routine wrappers

---
 boxtree/pyfmmlib_integration.py | 47 ++++++++++++++++++++++++++-------
 test/test_fmm.py                |  7 +++--
 2 files changed, 43 insertions(+), 11 deletions(-)

diff --git a/boxtree/pyfmmlib_integration.py b/boxtree/pyfmmlib_integration.py
index 3adccf3..0c848ee 100644
--- a/boxtree/pyfmmlib_integration.py
+++ b/boxtree/pyfmmlib_integration.py
@@ -26,7 +26,6 @@ THE SOFTWARE.
 
 
 import numpy as np
-from functools import partial
 
 
 __doc__ = """Integrates :mod:`boxtree` with
@@ -116,6 +115,33 @@ class HelmholtzExpansionWrangler:
         else:
             raise ValueError("unsupported dimensionality")
 
+    def get_expn_eval_routine(self, expn_kind):
+        name = "%%dd%seval" % expn_kind
+        rout = self.get_routine(name, "_vec")
+
+        if self.tree.dimensions == 2:
+            def wrapper(*args, **kwargs):
+                pot, grad, hess = rout(*args, **kwargs, ifgrad=False, ifhess=False)
+                return pot
+
+            from functools import update_wrapper
+            update_wrapper(wrapper, rout)
+            return wrapper
+
+        elif self.tree.dimensions == 3:
+            def wrapper(*args, **kwargs):
+                pot, fld, ier = rout(*args, **kwargs, iffld=False)
+                # grad = -fld
+                if (ier != 0).any():
+                    raise RuntimeError("%s failed with nonzero ier" % name)
+                return pot
+
+            from functools import update_wrapper
+            update_wrapper(wrapper, rout)
+            return wrapper
+        else:
+            raise ValueError("unsupported dimensionality")
+
     def multipole_expansion_zeros(self):
         if self.tree.dimensions == 2:
             return np.zeros((self.tree.nboxes, 2*self.nterms+1), dtype=np.complex128)
@@ -292,7 +318,7 @@ class HelmholtzExpansionWrangler:
 
         rscale = 1
 
-        mpeval = self.get_vec_routine("%ddmpeval")
+        mpeval = self.get_expn_eval_routine("mp")
 
         for ssn in sep_smaller_nonsiblings_by_level:
             for itgt_box, tgt_ibox in enumerate(target_boxes):
@@ -305,10 +331,9 @@ class HelmholtzExpansionWrangler:
                 start, end = ssn.starts[itgt_box:itgt_box+2]
                 for src_ibox in ssn.lists[start:end]:
 
-                    tmp_pot, _, _ = mpeval(self.helmholtz_k, rscale, self.
+                    tmp_pot = mpeval(self.helmholtz_k, rscale, self.
                             tree.box_centers[:, src_ibox], mpole_exps[src_ibox],
-                            self._get_targets(tgt_pslice),
-                            ifgrad=False, ifhess=False)
+                            self._get_targets(tgt_pslice))
 
                     tgt_pot = tgt_pot + tmp_pot
 
@@ -364,10 +389,14 @@ class HelmholtzExpansionWrangler:
                 src_ibox = self.tree.box_parent_ids[tgt_ibox]
                 src_center = self.tree.box_centers[:, src_ibox]
 
+                kwargs = {}
+                if self.tree.dimensions == 3:
+                    kwargs["radius"] = self.tree.root_extent * 2**(-target_lev)
+
                 tmp_loc_exp = locloc(
                             self.helmholtz_k,
                             rscale, src_center, local_exps[src_ibox],
-                            rscale, tgt_center, self.nterms)[..., 0]
+                            rscale, tgt_center, self.nterms, **kwargs)[..., 0]
 
                 local_exps[tgt_ibox] += tmp_loc_exp
 
@@ -377,7 +406,7 @@ class HelmholtzExpansionWrangler:
         pot = self.potential_zeros()
         rscale = 1  # FIXME
 
-        taeval = self.get_vec_routine("%ddtaeval")
+        taeval = self.get_expn_eval_routine("ta")
 
         for tgt_ibox in target_boxes:
             tgt_pslice = self._get_target_slice(tgt_ibox)
@@ -385,9 +414,9 @@ class HelmholtzExpansionWrangler:
             if tgt_pslice.stop - tgt_pslice.start == 0:
                 continue
 
-            tmp_pot, _, _ = taeval(self.helmholtz_k, rscale,
+            tmp_pot = taeval(self.helmholtz_k, rscale,
                     self.tree.box_centers[:, tgt_ibox], local_exps[tgt_ibox],
-                    self._get_targets(tgt_pslice), ifgrad=False, ifhess=False)
+                    self._get_targets(tgt_pslice))
 
             pot[tgt_pslice] += tmp_pot
 
diff --git a/test/test_fmm.py b/test/test_fmm.py
index 1bf3d84..29ad5fc 100644
--- a/test/test_fmm.py
+++ b/test/test_fmm.py
@@ -454,7 +454,7 @@ def test_fmm_completeness(ctx_getter, dims, nsources_req, ntargets_req,
 
 # {{{ test Helmholtz fmm with pyfmmlib
 
-@pytest.mark.parametrize("dims", [2])
+@pytest.mark.parametrize("dims", [2, 3])
 def test_pyfmmlib_fmm(ctx_getter, dims):
     logging.basicConfig(level=logging.INFO)
 
@@ -510,7 +510,10 @@ def test_pyfmmlib_fmm(ctx_getter, dims):
                 sources=sources_host.T, charge=weights,
                 targets=targets_host.T, zk=helmholtz_k)
     else:
-        raise NotImplementedError()
+        from pyfmmlib import hpotfld3dall_vec
+        ref_pot, _ = hpotfld3dall_vec(iffld=False,
+                sources=sources_host.T, charge=weights,
+                targets=targets_host.T, zk=helmholtz_k)
 
     rel_err = la.norm(pot - ref_pot) / la.norm(ref_pot)
     logger.info("relative l2 error: %g" % rel_err)
-- 
GitLab


From 17de894d982249dacd71f84317732d453cc8fdc5 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sun, 4 Jun 2017 20:09:45 -0400
Subject: [PATCH 6/9] Fix: Use https (rather than the git protocol) to fetch
 requirements

---
 requirements.txt | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/requirements.txt b/requirements.txt
index 2a67d20..b6fe300 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,5 @@
 numpy
-https+git://github.com/pyopencl/pyopencl
-https+git://github.com/inducer/islpy
-https+git://github.com/inducer/loopy
-https+git://github.com/inducer/pyfmmlib
+git+https://github.com/pyopencl/pyopencl
+git+https://github.com/inducer/islpy
+git+https://github.com/inducer/loopy
+git+https://github.com/inducer/pyfmmlib
-- 
GitLab


From b382b9a2326832555c35605f21dc17f3a1307615 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sun, 4 Jun 2017 23:06:22 -0400
Subject: [PATCH 7/9] Use legewhts instead of scipy to get Legendre nodes

---
 boxtree/pyfmmlib_integration.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/boxtree/pyfmmlib_integration.py b/boxtree/pyfmmlib_integration.py
index 0c848ee..875ac9e 100644
--- a/boxtree/pyfmmlib_integration.py
+++ b/boxtree/pyfmmlib_integration.py
@@ -46,12 +46,12 @@ class HelmholtzExpansionWrangler:
         common_extra_kwargs = {}
         if tree.dimensions == 3:
             nquad = max(6, int(2.5*nterms))
-            import scipy.special as sps
-            weights = sps.legendre(nquad).weights
+            from pyfmmlib import legewhts
+            xnodes, weights = legewhts(nquad, ifwhts=1)
 
             common_extra_kwargs = {
-                    "xnodes": weights[:, 0],
-                    "wts": weights[:, 2],
+                    "xnodes": xnodes,
+                    "wts": weights,
                     }
 
         self.common_extra_kwargs = common_extra_kwargs
-- 
GitLab


From a2930215ea26eebefe02dfbd93a641be59553ee6 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sun, 4 Jun 2017 23:26:19 -0400
Subject: [PATCH 8/9] Fix Python syntax for 2.7

---
 boxtree/pyfmmlib_integration.py | 14 ++++++++++----
 1 file changed, 10 insertions(+), 4 deletions(-)

diff --git a/boxtree/pyfmmlib_integration.py b/boxtree/pyfmmlib_integration.py
index 875ac9e..08b6d81 100644
--- a/boxtree/pyfmmlib_integration.py
+++ b/boxtree/pyfmmlib_integration.py
@@ -94,7 +94,9 @@ class HelmholtzExpansionWrangler:
             rout = self.get_routine("potgrad%ddall", "_vec")
 
             def wrapper(*args, **kwargs):
-                pot, grad, hess = rout(*args, **kwargs, ifgrad=False, ifhess=False)
+                kwargs["ifgrad"] = False
+                kwargs["ifhess"] = False
+                pot, grad, hess = rout(*args, **kwargs)
                 return pot
 
             from functools import update_wrapper
@@ -105,7 +107,8 @@ class HelmholtzExpansionWrangler:
             rout = self.get_routine("potfld%ddall", "_vec")
 
             def wrapper(*args, **kwargs):
-                pot, fld = rout(*args, **kwargs, iffld=False)
+                kwargs["iffld"] = False
+                pot, fld = rout(*args, **kwargs)
                 # grad = -fld
                 return pot
 
@@ -121,7 +124,9 @@ class HelmholtzExpansionWrangler:
 
         if self.tree.dimensions == 2:
             def wrapper(*args, **kwargs):
-                pot, grad, hess = rout(*args, **kwargs, ifgrad=False, ifhess=False)
+                kwargs["ifgrad"] = False
+                kwargs["ifhess"] = False
+                pot, grad, hess = rout(*args, **kwargs)
                 return pot
 
             from functools import update_wrapper
@@ -130,7 +135,8 @@ class HelmholtzExpansionWrangler:
 
         elif self.tree.dimensions == 3:
             def wrapper(*args, **kwargs):
-                pot, fld, ier = rout(*args, **kwargs, iffld=False)
+                kwargs["iffld"] = False
+                pot, fld, ier = rout(*args, **kwargs)
                 # grad = -fld
                 if (ier != 0).any():
                     raise RuntimeError("%s failed with nonzero ier" % name)
-- 
GitLab


From 0b5d48308174bcc97b667466b4c8096f4df9156b Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Mon, 5 Jun 2017 00:48:10 -0400
Subject: [PATCH 9/9] Don't use update_wrapper in h3d FMM wrapper: confuses Py2

---
 boxtree/pyfmmlib_integration.py | 25 +++++++++++++++----------
 1 file changed, 15 insertions(+), 10 deletions(-)

diff --git a/boxtree/pyfmmlib_integration.py b/boxtree/pyfmmlib_integration.py
index 08b6d81..6330b70 100644
--- a/boxtree/pyfmmlib_integration.py
+++ b/boxtree/pyfmmlib_integration.py
@@ -85,8 +85,9 @@ class HelmholtzExpansionWrangler:
 
                 return val
 
-            from functools import update_wrapper
-            update_wrapper(wrapper, rout)
+            # Doesn't work in in Py2
+            # from functools import update_wrapper
+            # update_wrapper(wrapper, rout)
             return wrapper
 
     def get_direct_eval_routine(self):
@@ -99,8 +100,9 @@ class HelmholtzExpansionWrangler:
                 pot, grad, hess = rout(*args, **kwargs)
                 return pot
 
-            from functools import update_wrapper
-            update_wrapper(wrapper, rout)
+            # Doesn't work in in Py2
+            # from functools import update_wrapper
+            # update_wrapper(wrapper, rout)
             return wrapper
 
         elif self.tree.dimensions == 3:
@@ -112,8 +114,9 @@ class HelmholtzExpansionWrangler:
                 # grad = -fld
                 return pot
 
-            from functools import update_wrapper
-            update_wrapper(wrapper, rout)
+            # Doesn't work in in Py2
+            # from functools import update_wrapper
+            # update_wrapper(wrapper, rout)
             return wrapper
         else:
             raise ValueError("unsupported dimensionality")
@@ -129,8 +132,9 @@ class HelmholtzExpansionWrangler:
                 pot, grad, hess = rout(*args, **kwargs)
                 return pot
 
-            from functools import update_wrapper
-            update_wrapper(wrapper, rout)
+            # Doesn't work in in Py2
+            # from functools import update_wrapper
+            # update_wrapper(wrapper, rout)
             return wrapper
 
         elif self.tree.dimensions == 3:
@@ -142,8 +146,9 @@ class HelmholtzExpansionWrangler:
                     raise RuntimeError("%s failed with nonzero ier" % name)
                 return pot
 
-            from functools import update_wrapper
-            update_wrapper(wrapper, rout)
+            # Doesn't work in in Py2
+            # from functools import update_wrapper
+            # update_wrapper(wrapper, rout)
             return wrapper
         else:
             raise ValueError("unsupported dimensionality")
-- 
GitLab