From f53c8875cb0723888c8e1cbc822466f359753066 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 27 Aug 2017 22:28:24 -0500 Subject: [PATCH 1/5] Make formmp in fmmlib execute level-by-level --- boxtree/pyfmmlib_integration.py | 40 ++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 16 deletions(-) diff --git a/boxtree/pyfmmlib_integration.py b/boxtree/pyfmmlib_integration.py index 833918f..34636f6 100644 --- a/boxtree/pyfmmlib_integration.py +++ b/boxtree/pyfmmlib_integration.py @@ -365,27 +365,35 @@ class FMMLibExpansionWrangler(object): formmp = self.get_routine("%ddformmp" + self.dp_suffix) mpoles = self.multipole_expansion_zeros() - for src_ibox in source_boxes: - pslice = self._get_source_slice(src_ibox) - - if pslice.stop - pslice.start == 0: + for lev in range(self.tree.nlevels): + start, stop = level_start_source_box_nrs[lev:lev+2] + if start == stop: continue - kwargs = {} - kwargs.update(self.kernel_kwargs) - kwargs.update(self.get_source_kwargs(src_weights, pslice)) + level_start_ibox, mpoles_view = self.multipole_expansions_view( + mpoles, lev) - ier, mpole = formmp( - rscale=rscale, - source=self._get_sources(pslice), - center=self.tree.box_centers[:, src_ibox], - nterms=self.nterms, - **kwargs) + for src_ibox in source_boxes[start:stop]: + pslice = self._get_source_slice(src_ibox) + + if pslice.stop - pslice.start == 0: + continue - if ier: - raise RuntimeError("formmp failed") + kwargs = {} + kwargs.update(self.kernel_kwargs) + kwargs.update(self.get_source_kwargs(src_weights, pslice)) + + ier, mpole = formmp( + rscale=rscale, + source=self._get_sources(pslice), + center=self.tree.box_centers[:, src_ibox], + nterms=self.nterms, + **kwargs) + + if ier: + raise RuntimeError("formmp failed") - mpoles[src_ibox] = mpole.T + mpoles_view[src_ibox-level_start_ibox] = mpole.T return mpoles -- GitLab From 8c060ac996801c1571748c959733a9c8a9f3e0ea Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 27 Aug 2017 22:50:17 -0500 Subject: [PATCH 2/5] Make formta in fmmlib execute level-by-level --- boxtree/pyfmmlib_integration.py | 54 +++++++++++++++++++-------------- 1 file changed, 32 insertions(+), 22 deletions(-) diff --git a/boxtree/pyfmmlib_integration.py b/boxtree/pyfmmlib_integration.py index 34636f6..c173ff1 100644 --- a/boxtree/pyfmmlib_integration.py +++ b/boxtree/pyfmmlib_integration.py @@ -564,7 +564,7 @@ class FMMLibExpansionWrangler(object): mpeval = self.get_expn_eval_routine("mp") - for ssn in sep_smaller_nonsiblings_by_level: + for isrc_level, ssn in enumerate(sep_smaller_nonsiblings_by_level): for itgt_box, tgt_ibox in enumerate(target_boxes): tgt_pslice = self._get_target_slice(tgt_ibox) @@ -599,34 +599,44 @@ class FMMLibExpansionWrangler(object): formta = self.get_routine("%ddformta" + self.dp_suffix) - for itgt_box, tgt_ibox in enumerate(target_or_target_parent_boxes): - start, end = starts[itgt_box:itgt_box+2] + for lev in range(self.tree.nlevels): + lev_start, lev_stop = \ + level_start_target_or_target_parent_box_nrs[lev:lev+2] + if lev_start == lev_stop: + continue - contrib = 0 + target_level_start_ibox, target_local_exps_view = \ + self.local_expansions_view(local_exps, lev) - for src_ibox in lists[start:end]: - src_pslice = self._get_source_slice(src_ibox) - tgt_center = self.tree.box_centers[:, tgt_ibox] + for itgt_box, tgt_ibox in enumerate( + target_or_target_parent_boxes[lev_start:lev_stop]): + start, end = starts[lev_start+itgt_box:lev_start+itgt_box+2] - if src_pslice.stop - src_pslice.start == 0: - continue + contrib = 0 - kwargs = {} - kwargs.update(self.kernel_kwargs) - kwargs.update(self.get_source_kwargs(src_weights, src_pslice)) + for src_ibox in lists[start:end]: + src_pslice = self._get_source_slice(src_ibox) + tgt_center = self.tree.box_centers[:, tgt_ibox] - ier, mpole = formta( - rscale=rscale, - source=self._get_sources(src_pslice), - center=tgt_center, - nterms=self.nterms, - **kwargs) - if ier: - raise RuntimeError("formta failed") + if src_pslice.stop - src_pslice.start == 0: + continue + + kwargs = {} + kwargs.update(self.kernel_kwargs) + kwargs.update(self.get_source_kwargs(src_weights, src_pslice)) + + ier, mpole = formta( + rscale=rscale, + source=self._get_sources(src_pslice), + center=tgt_center, + nterms=self.nterms, + **kwargs) + if ier: + raise RuntimeError("formta failed") - contrib = contrib + mpole.T + contrib = contrib + mpole.T - local_exps[tgt_ibox] = contrib + target_local_exps_view[tgt_ibox-target_level_start_ibox] = contrib return local_exps -- GitLab From b99ddad19df5942b3e9fcafe94ca5e0aa83e3d5f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 27 Aug 2017 23:12:49 -0500 Subject: [PATCH 3/5] Make upward propagation and M2L in fmmlib use expansion views --- boxtree/pyfmmlib_integration.py | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/boxtree/pyfmmlib_integration.py b/boxtree/pyfmmlib_integration.py index c173ff1..d02beb4 100644 --- a/boxtree/pyfmmlib_integration.py +++ b/boxtree/pyfmmlib_integration.py @@ -415,6 +415,11 @@ class FMMLibExpansionWrangler(object): start, stop = level_start_source_parent_box_nrs[ target_level:target_level+2] + source_level_start_ibox, source_mpoles_view = \ + self.multipole_expansions_view(mpoles, source_level) + target_level_start_ibox, target_mpoles_view = \ + self.multipole_expansions_view(mpoles, target_level) + for ibox in source_parent_boxes[start:stop]: parent_center = tree.box_centers[:, ibox] for child in tree.box_child_ids[:, ibox]: @@ -430,7 +435,8 @@ class FMMLibExpansionWrangler(object): new_mp = mpmp( rscale1=rscale, center1=child_center, - expn1=mpoles[child].T, + expn1=source_mpoles_view[ + child - source_level_start_ibox].T, rscale2=rscale, center2=parent_center, @@ -438,7 +444,8 @@ class FMMLibExpansionWrangler(object): **kwargs) - mpoles[ibox] += new_mp[..., 0].T + target_mpoles_view[ + ibox - target_level_start_ibox] += new_mp[..., 0].T def eval_direct(self, target_boxes, neighbor_sources_starts, neighbor_sources_lists, src_weights): @@ -496,6 +503,11 @@ class FMMLibExpansionWrangler(object): starts_on_lvl = starts[lstart:lstop+1] + source_level_start_ibox, source_mpoles_view = \ + self.multipole_expansions_view(mpole_exps, lev) + target_level_start_ibox, target_local_exps_view = \ + self.local_expansions_view(local_exps, lev) + ntgt_boxes = lstop-lstart itgt_box_vec = np.arange(ntgt_boxes) tgt_ibox_vec = target_or_target_parent_boxes[lstart:lstop] @@ -542,8 +554,8 @@ class FMMLibExpansionWrangler(object): center1_offsets=lists, center1_starts=starts_on_lvl, - expn1=mpole_exps.T, - expn1_offsets=lists, + expn1=source_mpoles_view.T, + expn1_offsets=lists - source_level_start_ibox, expn1_starts=starts_on_lvl, rscale2=rscale2, @@ -552,7 +564,7 @@ class FMMLibExpansionWrangler(object): expn2=expn2.T, **kwargs).T - local_exps[tgt_ibox_vec] += expn2 + target_local_exps_view[tgt_ibox_vec - target_level_start_ibox] += expn2 return local_exps -- GitLab From 57ed422c1ac59f210d397a4d57f8be400b30465e Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 27 Aug 2017 23:38:52 -0500 Subject: [PATCH 4/5] Make pyfmmlib fmm completely level-by-level and have it use expansion views throughout --- boxtree/pyfmmlib_integration.py | 52 +++++++++++++++++++++++---------- 1 file changed, 37 insertions(+), 15 deletions(-) diff --git a/boxtree/pyfmmlib_integration.py b/boxtree/pyfmmlib_integration.py index d02beb4..89bde84 100644 --- a/boxtree/pyfmmlib_integration.py +++ b/boxtree/pyfmmlib_integration.py @@ -577,6 +577,9 @@ class FMMLibExpansionWrangler(object): mpeval = self.get_expn_eval_routine("mp") for isrc_level, ssn in enumerate(sep_smaller_nonsiblings_by_level): + source_level_start_ibox, source_mpoles_view = \ + self.multipole_expansions_view(mpole_exps, isrc_level) + for itgt_box, tgt_ibox in enumerate(target_boxes): tgt_pslice = self._get_target_slice(tgt_ibox) @@ -591,7 +594,8 @@ class FMMLibExpansionWrangler(object): tmp_pot, tmp_grad = mpeval( rscale=rscale, center=self.tree.box_centers[:, src_ibox], - expn=mpole_exps[src_ibox].T, + expn=source_mpoles_view[ + src_ibox - source_level_start_ibox].T, ztarg=self._get_targets(tgt_pslice), **self.kernel_kwargs) @@ -662,6 +666,13 @@ class FMMLibExpansionWrangler(object): start, stop = level_start_target_or_target_parent_box_nrs[ target_lev:target_lev+2] + source_lev = target_lev - 1 + + source_level_start_ibox, source_local_exps_view = \ + self.local_expansions_view(local_exps, source_lev) + target_level_start_ibox, target_local_exps_view = \ + self.local_expansions_view(local_exps, target_lev) + for tgt_ibox in target_or_target_parent_boxes[start:stop]: tgt_center = self.tree.box_centers[:, tgt_ibox] src_ibox = self.tree.box_parent_ids[tgt_ibox] @@ -675,7 +686,8 @@ class FMMLibExpansionWrangler(object): tmp_loc_exp = locloc( rscale1=rscale, center1=src_center, - expn1=local_exps[src_ibox].T, + expn1=source_local_exps_view[ + src_ibox - source_level_start_ibox].T, rscale2=rscale, center2=tgt_center, @@ -683,7 +695,8 @@ class FMMLibExpansionWrangler(object): **kwargs)[..., 0] - local_exps[tgt_ibox] += tmp_loc_exp.T + target_local_exps_view[ + tgt_ibox - target_level_start_ibox] += tmp_loc_exp.T return local_exps @@ -693,22 +706,31 @@ class FMMLibExpansionWrangler(object): taeval = self.get_expn_eval_routine("ta") - for tgt_ibox in target_boxes: - tgt_pslice = self._get_target_slice(tgt_ibox) - - if tgt_pslice.stop - tgt_pslice.start == 0: + for lev in range(self.tree.nlevels): + start, stop = level_start_target_box_nrs[lev:lev+2] + if start == stop: continue - tmp_pot, tmp_grad = taeval( - rscale=rscale, - center=self.tree.box_centers[:, tgt_ibox], - expn=local_exps[tgt_ibox].T, - ztarg=self._get_targets(tgt_pslice), + source_level_start_ibox, source_local_exps_view = \ + self.local_expansions_view(local_exps, lev) - **self.kernel_kwargs) + for tgt_ibox in target_boxes[start:stop]: + tgt_pslice = self._get_target_slice(tgt_ibox) - self.add_potgrad_onto_output( - output, tgt_pslice, tmp_pot, tmp_grad) + if tgt_pslice.stop - tgt_pslice.start == 0: + continue + + tmp_pot, tmp_grad = taeval( + rscale=rscale, + center=self.tree.box_centers[:, tgt_ibox], + expn=source_local_exps_view[ + tgt_ibox - source_level_start_ibox].T, + ztarg=self._get_targets(tgt_pslice), + + **self.kernel_kwargs) + + self.add_potgrad_onto_output( + output, tgt_pslice, tmp_pot, tmp_grad) return output -- GitLab From 256f77042288692752c87eb01aa7477163679499 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 27 Aug 2017 23:49:33 -0500 Subject: [PATCH 5/5] Enable use of rscale in fmmlib FMM --- boxtree/pyfmmlib_integration.py | 42 ++++++++++++++++++++------------- 1 file changed, 25 insertions(+), 17 deletions(-) diff --git a/boxtree/pyfmmlib_integration.py b/boxtree/pyfmmlib_integration.py index 89bde84..9afebaf 100644 --- a/boxtree/pyfmmlib_integration.py +++ b/boxtree/pyfmmlib_integration.py @@ -34,6 +34,10 @@ __doc__ = """Integrates :mod:`boxtree` with """ +def level_to_rscale(tree, level): + return tree.root_extent * 2 ** -level + + class FMMLibExpansionWrangler(object): """Implements the :class:`boxtree.fmm.ExpansionWranglerInterface` by using pyfmmlib. @@ -360,8 +364,6 @@ class FMMLibExpansionWrangler(object): } def form_multipoles(self, level_start_source_box_nrs, source_boxes, src_weights): - rscale = 1 # FIXME - formmp = self.get_routine("%ddformmp" + self.dp_suffix) mpoles = self.multipole_expansion_zeros() @@ -373,6 +375,8 @@ class FMMLibExpansionWrangler(object): level_start_ibox, mpoles_view = self.multipole_expansions_view( mpoles, lev) + rscale = level_to_rscale(self.tree, lev) + for src_ibox in source_boxes[start:stop]: pslice = self._get_source_slice(src_ibox) @@ -400,7 +404,6 @@ class FMMLibExpansionWrangler(object): def coarsen_multipoles(self, level_start_source_parent_box_nrs, source_parent_boxes, mpoles): tree = self.tree - rscale = 1 # FIXME mpmp = self.get_translation_routine("%ddmpmp") @@ -420,6 +423,9 @@ class FMMLibExpansionWrangler(object): target_level_start_ibox, target_mpoles_view = \ self.multipole_expansions_view(mpoles, target_level) + source_rscale = level_to_rscale(tree, source_level) + target_rscale = level_to_rscale(tree, target_level) + for ibox in source_parent_boxes[start:stop]: parent_center = tree.box_centers[:, ibox] for child in tree.box_child_ids[:, ibox]: @@ -433,12 +439,12 @@ class FMMLibExpansionWrangler(object): kwargs.update(self.kernel_kwargs) new_mp = mpmp( - rscale1=rscale, + rscale1=source_rscale, center1=child_center, expn1=source_mpoles_view[ child - source_level_start_ibox].T, - rscale2=rscale, + rscale2=target_rscale, center2=parent_center, nterms2=self.nterms, @@ -521,8 +527,9 @@ class FMMLibExpansionWrangler(object): src_boxes_starts[0] = 0 src_boxes_starts[1:] = np.cumsum(nsrc_boxes_per_tgt_box) - # FIXME - rscale1 = np.ones(nsrc_boxes) + rscale = level_to_rscale(tree, lev) + + rscale1 = np.ones(nsrc_boxes) * rscale rscale1_offsets = np.arange(nsrc_boxes) kwargs = {} @@ -531,8 +538,7 @@ class FMMLibExpansionWrangler(object): tree.root_extent * 2**(-lev) * np.ones(ntgt_boxes)) - # FIXME - rscale2 = np.ones(ntgt_boxes, np.float64) + rscale2 = np.ones(ntgt_boxes, np.float64) * rscale # These get max'd/added onto: pass initialized versions. if self.dim == 3: @@ -572,14 +578,14 @@ class FMMLibExpansionWrangler(object): sep_smaller_nonsiblings_by_level, mpole_exps): output = self.output_zeros() - rscale = 1 - mpeval = self.get_expn_eval_routine("mp") for isrc_level, ssn in enumerate(sep_smaller_nonsiblings_by_level): source_level_start_ibox, source_mpoles_view = \ self.multipole_expansions_view(mpole_exps, isrc_level) + rscale = level_to_rscale(self.tree, isrc_level) + for itgt_box, tgt_ibox in enumerate(target_boxes): tgt_pslice = self._get_target_slice(tgt_ibox) @@ -610,7 +616,6 @@ class FMMLibExpansionWrangler(object): def form_locals(self, level_start_target_or_target_parent_box_nrs, target_or_target_parent_boxes, starts, lists, src_weights): - rscale = 1 # FIXME local_exps = self.local_expansion_zeros() formta = self.get_routine("%ddformta" + self.dp_suffix) @@ -624,6 +629,8 @@ class FMMLibExpansionWrangler(object): target_level_start_ibox, target_local_exps_view = \ self.local_expansions_view(local_exps, lev) + rscale = level_to_rscale(self.tree, lev) + for itgt_box, tgt_ibox in enumerate( target_or_target_parent_boxes[lev_start:lev_stop]): start, end = starts[lev_start+itgt_box:lev_start+itgt_box+2] @@ -658,7 +665,6 @@ class FMMLibExpansionWrangler(object): def refine_locals(self, level_start_target_or_target_parent_box_nrs, target_or_target_parent_boxes, local_exps): - rscale = 1 # FIXME locloc = self.get_translation_routine("%ddlocloc") @@ -672,6 +678,8 @@ class FMMLibExpansionWrangler(object): self.local_expansions_view(local_exps, source_lev) target_level_start_ibox, target_local_exps_view = \ self.local_expansions_view(local_exps, target_lev) + source_rscale = level_to_rscale(self.tree, source_lev) + target_rscale = level_to_rscale(self.tree, target_lev) for tgt_ibox in target_or_target_parent_boxes[start:stop]: tgt_center = self.tree.box_centers[:, tgt_ibox] @@ -684,12 +692,12 @@ class FMMLibExpansionWrangler(object): kwargs.update(self.kernel_kwargs) tmp_loc_exp = locloc( - rscale1=rscale, + rscale1=source_rscale, center1=src_center, expn1=source_local_exps_view[ src_ibox - source_level_start_ibox].T, - rscale2=rscale, + rscale2=target_rscale, center2=tgt_center, nterms2=self.nterms, @@ -702,8 +710,6 @@ class FMMLibExpansionWrangler(object): def eval_locals(self, level_start_target_box_nrs, target_boxes, local_exps): output = self.output_zeros() - rscale = 1 # FIXME - taeval = self.get_expn_eval_routine("ta") for lev in range(self.tree.nlevels): @@ -714,6 +720,8 @@ class FMMLibExpansionWrangler(object): source_level_start_ibox, source_local_exps_view = \ self.local_expansions_view(local_exps, lev) + rscale = level_to_rscale(self.tree, lev) + for tgt_ibox in target_boxes[start:stop]: tgt_pslice = self._get_target_slice(tgt_ibox) -- GitLab