Start PetscObject

devito/core/cpu.py

@@ -182,10 +182,8 @@ def _specialize_iet(cls, graph, **kwargs):
 
         # Misc optimizations
         hoist_prodders(graph)
-
         # Symbol definitions
         cls._Target.DataManager(**kwargs).process(graph)
-
         # Linearize n-dimensional Indexeds
         linearize(graph, **kwargs)
 

devito/finite_differences/finite_difference.py

@@ -12,6 +12,10 @@
 # code generation
 _PRECISION = 9
 
+# could add:
+# import numpy as np
+# _PRECISION = _PRECISION = {np.float32: 9, np.float64: 16}
+
 
 @check_input
 @check_symbolic
@@ -222,6 +226,9 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, symbolic
     # Enforce fixed precision FD coefficients to avoid variations in results
     weights = [sympify(w).evalf(_PRECISION) for w in weights][::matvec.val]
 
+    #if using mapper for precision, could have:
+    # weights = [sympify(w).evalf(_PRECISION[expr.dtype]) for w in weights]
+
     # Transpose the FD, if necessary
     if matvec == transpose:
         indices = indices.transpose()

devito/ir/clusters/cluster.py

@@ -343,16 +343,17 @@ def dspace(self):
                     d = i.dim.parent
                 else:
                     d = i.dim
-                try:
-                    if i.lower < 0 or \
-                       i.upper > f._size_nodomain[d].left + f._size_halo[d].right:
-                        # It'd mean trying to access a point before the
-                        # left halo (test0) or after the right halo (test1)
-                        oobs.update(d._defines)
-                except (KeyError, TypeError):
-                    # Unable to detect presence of OOB accesses (e.g., `d` not in
-                    # `f._size_halo`, that is typical of indirect accesses `A[B[i]]`)
-                    pass
+                # TMP turn off for PETSc
+                # try:
+                #     if i.lower < 0 or \
+                #        i.upper > f._size_nodomain[d].left + f._size_halo[d].right:
+                #         # It'd mean trying to access a point before the
+                #         # left halo (test0) or after the right halo (test1)
+                #         oobs.update(d._defines)
+                # except (KeyError, TypeError):
+                #     # Unable to detect presence of OOB accesses (e.g., `d` not in
+                #     # `f._size_halo`, that is typical of indirect accesses `A[B[i]]`)
+                #     pass
 
         # Construct the `intervals` of the DataSpace, that is a global,
         # Dimension-centric view of the data space

devito/ir/equations/algorithms.py

@@ -120,9 +120,12 @@ def lower_exprs(expressions, **kwargs):
         for i in retrieve_indexed(expr):
             f = i.function
 
-            # Introduce shifting to align with the computational domain
-            indices = [(lower_exprs(a) + o) for a, o in
-                       zip(i.indices, f._size_nodomain.left)]
+            # # Introduce shifting to align with the computational domain
+            # indices = [(lower_exprs(a) + o) for a, o in
+            #            zip(i.indices, f._size_nodomain.left)]
+
+            # tmp turn off this shifting in order to utilise DMDA in PETSc
+            indices = [lower_exprs(a) for a in i.indices]
 
             # Substitute spacing (spacing only used in own dimension)
             indices = [i.xreplace({d.spacing: 1, -d.spacing: -1})

devito/ir/iet/nodes.py

@@ -27,7 +27,7 @@
            'AugmentedExpression', 'Increment', 'Return', 'While',
            'ParallelIteration', 'ParallelBlock', 'Dereference', 'Lambda',
            'SyncSpot', 'Pragma', 'DummyExpr', 'BlankLine', 'ParallelTree',
-           'BusyWait', 'CallableBody', 'Transfer']
+           'BusyWait', 'CallableBody', 'Transfer', 'CallBack']
 
 # First-class IET nodes
 
@@ -729,6 +729,19 @@ def defines(self):
         return self.all_parameters
 
 
+class CallBack(Call):
+
+    """
+    Implement a function pointer argument for a callback function.
+    """
+    def __init__(self, name, return_type, parameter_type, **kwargs):
+
+        super().__init__(name=name)
+
+        self.return_type = return_type
+        self.parameter_type = parameter_type
+
+
 class CallableBody(Node):
 
     """

devito/ir/iet/utils.py

@@ -121,7 +121,6 @@ def derive_parameters(iet, drop_locals=False):
     # Maybe filter out all other compiler-generated objects
     if drop_locals:
         parameters = [p for p in parameters if not (p.is_ArrayBasic or p.is_LocalObject)]
-
     return parameters
 
 

devito/ir/iet/visitors.py

@@ -10,12 +10,12 @@
 import ctypes
 
 import cgen as c
-from sympy import IndexedBase
+from sympy import IndexedBase, sympify
 from sympy.core.function import Application
 
 from devito.exceptions import VisitorException
 from devito.ir.iet.nodes import (Node, Iteration, Expression, ExpressionBundle,
-                                 Call, Lambda, BlankLine, Section)
+                                 Call, CallBack, Lambda, BlankLine, Section)
 from devito.ir.support.space import Backward
 from devito.symbolics import ListInitializer, ccode, uxreplace
 from devito.tools import (GenericVisitor, as_tuple, ctypes_to_cstr, filter_ordered,
@@ -196,7 +196,6 @@ def _gen_struct_decl(self, obj, masked=()):
 
         if not issubclass(ctype, ctypes.Structure):
             return None
-
         try:
             return obj._C_typedecl
         except AttributeError:
@@ -299,6 +298,7 @@ def _gen_signature(self, o):
             tparams = ', '.join([i.inline() for i in self._args_decl(o.templates)])
             signature = c.Template(tparams, signature)
         return signature
+
 
     def _blankline_logic(self, children):
         """
@@ -483,6 +483,12 @@ def visit_Call(self, o, nested_call=False):
                 rettype = self._gen_rettype(retobj)
                 return c.Initializer(c.Value(rettype, retobj._C_name), call)
 
+    def visit_CallBack(self, o, nested_call=False):
+        name = o.name
+        return_type = o.return_type
+        parameter_type = o.parameter_type
+        return FunctionPointerArg(name, return_type, parameter_type)
+
     def visit_Conditional(self, o):
         try:
             then_body, else_body = self._blankline_logic(o.children)
@@ -646,7 +652,6 @@ def visit_Operator(self, o, mode='all'):
         globs = self._operator_globals(o, mode)
         if globs:
             globs.append(blankline)
-
         return c.Module(headers + includes + typedecls + globs +
                         esigns + [blankline, kernel] + efuncs)
 
@@ -1253,3 +1258,16 @@ def generate(self):
         if self.cast:
             tip = '(%s)%s' % (self.cast, tip)
         yield tip
+
+
+class FunctionPointerArg(c.Generable):
+
+    def __init__(self, name, return_type, parameter_type):
+        self.name = name
+        self.return_type = return_type
+        self.parameter_type = parameter_type
+
+    def generate(self):
+        yield "(%s (*)(%s))%s" % (self.return_type, self.parameter_type, self.name)
+
+

devito/operator/operator.py

@@ -20,7 +20,7 @@
 from devito.mpi import MPI
 from devito.parameters import configuration
 from devito.passes import (Graph, lower_index_derivatives, generate_implicit,
-                           generate_macros, minimize_symbols, unevaluate)
+                           generate_macros, minimize_symbols, unevaluate, lower_petsc)
 from devito.symbolics import estimate_cost
 from devito.tools import (DAG, OrderedSet, Signer, ReducerMap, as_tuple, flatten,
                           filter_sorted, frozendict, is_integer, split, timed_pass,
@@ -175,7 +175,6 @@ def _check_kwargs(cls, **kwargs):
     def _build(cls, expressions, **kwargs):
         # Python- (i.e., compile-) and C-level (i.e., run-time) performance
         profiler = create_profile('timers')
-
         # Lower the input expressions into an IET
         irs, byproduct = cls._lower(expressions, profiler=profiler, **kwargs)
 
@@ -247,10 +246,8 @@ def _lower(cls, expressions, **kwargs):
         # This may be used by a compilation pass that constructs a new
         # expression for which a partial or complete lowering is desired
         kwargs['rcompile'] = cls._rcompile_wrapper(**kwargs)
-
         # [Eq] -> [LoweredEq]
         expressions = cls._lower_exprs(expressions, **kwargs)
-
         # [LoweredEq] -> [Clusters]
         clusters = cls._lower_clusters(expressions, **kwargs)
 
@@ -327,9 +324,7 @@ def _lower_exprs(cls, expressions, **kwargs):
 
         # "True" lowering (indexification, shifting, ...)
         expressions = lower_exprs(expressions, **kwargs)
-
         processed = [LoweredEq(i) for i in expressions]
-
         return processed
 
     # Compilation -- Cluster level
@@ -448,7 +443,6 @@ def _lower_iet(cls, uiet, profiler=None, **kwargs):
         """
         name = kwargs.get("name", "Kernel")
         sregistry = kwargs['sregistry']
-
         # Wrap the IET with an EntryFunction (a special Callable representing
         # the entry point of the generated library)
         parameters = derive_parameters(uiet, True)
@@ -461,14 +455,15 @@ def _lower_iet(cls, uiet, profiler=None, **kwargs):
         # Instrument the IET for C-level profiling
         # Note: this is postponed until after _specialize_iet because during
         # specialization further Sections may be introduced
-        cls._Target.instrument(graph, profiler=profiler, **kwargs)
 
+        cls._Target.instrument(graph, profiler=profiler, **kwargs)
         # Extract the necessary macros from the symbolic objects
         generate_macros(graph)
-
         # Target-independent optimizations
         minimize_symbols(graph)
 
+        lower_petsc(graph, **kwargs)
+
         return graph.root, graph
 
     # Read-only properties exposed to the outside world
@@ -748,7 +743,6 @@ def cinterface(self, force=False):
 
         cfile = name.with_suffix(".%s" % self._compiler.src_ext)
         hfile = name.with_suffix('.h')
-
         # Generate the .c and .h code
         ccode, hcode = CInterface().visit(self)
 
@@ -832,7 +826,6 @@ def apply(self, **kwargs):
         # Build the arguments list to invoke the kernel function
         with self._profiler.timer_on('arguments'):
             args = self.arguments(**kwargs)
-
         # Invoke kernel function with args
         arg_values = [args[p.name] for p in self.parameters]
         try:

devito/passes/iet/__init__.py

@@ -7,3 +7,4 @@
 from .asynchrony import *  # noqa
 from .instrument import *  # noqa
 from .languages import *  # noqa
+from .petsc import *  # noqa

devito/passes/iet/definitions.py

@@ -322,6 +322,7 @@ def place_definitions(self, iet, globs=None, **kwargs):
                 self._alloc_scalar_on_low_lat_mem((iet,) + v, k, storage)
 
         iet, _ = self._inject_definitions(iet, storage)
+
 
         # Process all other definitions, essentially all temporary objects
         # created by the compiler up to this point (Array, LocalObject, etc.)
@@ -359,8 +360,12 @@ def place_definitions(self, iet, globs=None, **kwargs):
                 if v:
                     includes.add(v)
 
+        # it is at this point storage has been populated with the efunc foo and 
+        # a <devito.passes.iet....> so after calling _inject_definitions again
+        # the ccode inclues cudaStreamCreate etc
         iet, efuncs = self._inject_definitions(iet, storage)
 
+
         return iet, {'efuncs': efuncs,
                      'globals': as_tuple(globs),
                      'includes': as_tuple(includes)}

devito/passes/iet/engine.py

@@ -225,7 +225,6 @@ def reuse_efuncs(root, efuncs, sregistry=None):
     root = efuncs.pop(root.name)
     processed = [afunc if len(efuncs) > 1 else efuncs.pop()
                  for afunc, efuncs in mapper.values()]
-
     return root, processed
 
 
@@ -245,7 +244,6 @@ def abstract_efunc(efunc):
 
     efunc = Uxreplace(mapper).visit(efunc)
     efunc = efunc._rebuild(name='foo')
-
     return efunc
 
 

devito/passes/iet/instrument.py

@@ -19,7 +19,6 @@ def instrument(graph, **kwargs):
     if profiler is None:
         return
     timer = Timer(profiler.name, list(profiler.all_sections))
-
     instrument_sections(graph, timer=timer, **kwargs)
     sync_sections(graph, **kwargs)
 
@@ -103,14 +102,11 @@ def instrument_sections(iet, **kwargs):
     """
     profiler = kwargs['profiler']
     timer = kwargs['timer']
-
     piet = profiler.instrument(iet, timer)
-
     if piet is iet:
         return piet, {}
 
     headers = [TimedList._start_timer_header(), TimedList._stop_timer_header()]
-
     return piet, {'headers': headers}
 
 

devito/passes/iet/mpi.py

@@ -297,7 +297,6 @@ def make_mpi(iet, mpimode=None, **kwargs):
     """
     # To produce unique object names
     generators = {'msg': generator(), 'comm': generator(), 'comp': generator()}
-
     sync_heb = HaloExchangeBuilder('basic', generators, **kwargs)
     user_heb = HaloExchangeBuilder(mpimode, generators, **kwargs)
     mapper = {}
@@ -307,7 +306,6 @@ def make_mpi(iet, mpimode=None, **kwargs):
 
     efuncs = sync_heb.efuncs + user_heb.efuncs
     iet = Transformer(mapper, nested=True).visit(iet)
-
     # Must drop the PARALLEL tag from the Iterations within which halo
     # exchanges are performed
     mapper = {}