Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
36 changes: 35 additions & 1 deletion pyop2/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -694,12 +694,17 @@ def data(self, value):

@property
def soa(self):
"""Are the data in SoA format? This is always false for :class:`Global`
objects."""
return False

#FIXME: Part of kernel API, but must be declared before Map for the validation.

class IterationIndex(object):
"""OP2 iteration space index"""
"""OP2 iteration space index

Users should not directly instantiate :class:`IterationIndex` objects. Use
``op2.i`` instead."""

def __init__(self, index=None):
assert index is None or isinstance(index, int), "i must be an int"
Expand All @@ -713,6 +718,7 @@ def __repr__(self):

@property
def index(self):
"""Return the integer value of this index."""
return self._index

def __getitem__(self, idx):
Expand Down Expand Up @@ -1032,6 +1038,7 @@ def code(self):

@property
def md5(self):
"""MD5 digest of kernel code and name."""
if not hasattr(self, '_md5'):
import md5
self._md5 = md5.new(self._code + self._name).hexdigest()
Expand All @@ -1052,6 +1059,11 @@ def _parloop_cache_size():
return len(_parloop_cache)

class ParLoop(object):
"""Represents the kernel, iteration space and arguments of a parallel loop
invocation.

Users should not directly construct :class:`ParLoop` objects, but use
``op2.par_loop()`` instead."""
def __init__(self, kernel, itspace, *args):
self._kernel = kernel
if isinstance(itspace, IterationSpace):
Expand Down Expand Up @@ -1086,16 +1098,27 @@ def reduction_begin(self):
arg.reduction_begin()

def reduction_end(self):
"""End reductions"""
for arg in self.args:
if arg._is_global_reduction:
arg.reduction_end()

def maybe_set_halo_update_needed(self):
"""Set halo update needed for :class:`Dat` arguments that are written to
in this parallel loop."""
for arg in self.args:
if arg._is_dat and arg.access in [INC, WRITE, RW]:
arg.data.needs_halo_update = True

def check_args(self):
"""Checks the following:

1. That the iteration set of the :class:`ParLoop` matches the iteration
set of all its arguments.
2. For each argument, check that the dataset of the map used to access
it matches the dataset it is defined on.

A :class:`MapValueError` is raised if these conditions are not met."""
iterset = self._it_space._iterset
for i, arg in enumerate(self._actual_args):
if arg._is_global or arg._map == IdentityMap:
Expand All @@ -1116,27 +1139,34 @@ def generate_code(self):

@property
def it_space(self):
"""Iteration space of the parallel loop."""
return self._it_space

@property
def is_direct(self):
"""Is this parallel loop direct? I.e. are all the arguments either
:class:Dats accessed through the identity map, or :class:Global?"""
return all(a.map in [None, IdentityMap] for a in self.args)

@property
def is_indirect(self):
"""Is the parallel loop indirect?"""
return not self.is_direct

@property
def needs_exec_halo(self):
"""Does the parallel loop need an exec halo?"""
return any(arg._is_indirect_and_not_read or arg._is_mat
for arg in self.args)

@property
def kernel(self):
"""Kernel executed by this parallel loop."""
return self._kernel

@property
def args(self):
"""Arguments to this parallel loop."""
return self._actual_args

@property
Expand Down Expand Up @@ -1184,6 +1214,9 @@ def _cache_key(self):
'plot_prefix': '',
'error_on_nonconvergence': True,
'gmres_restart': 30}
"""The default parameters for the solver are the same as those used in PETSc
3.3. Note that the parameters accepted by :class:`op2.Solver` are only a subset
of all PETSc parameters."""

class Solver(object):
"""OP2 Solver object. The :class:`Solver` holds a set of parameters that are
Expand All @@ -1208,6 +1241,7 @@ class Solver(object):
:arg plot_convergence: plot a graph of the convergence history after the
solve has finished and save it to file (False, implies monitor_convergence)
:arg plot_prefix: filename prefix for plot files ('')
:arg gmres_restart: restart period when using GMRES

"""

Expand Down
49 changes: 48 additions & 1 deletion pyop2/runtime_base.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,13 @@
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.

""" Base classes for OP2 objects. The versions here extend those from the :mod:`base` module to include runtime data information which is backend independent. Individual runtime backends should subclass these as required to implement backend-specific features."""
"""Base classes for OP2 objects. The versions here extend those from the
:mod:`base` module to include runtime data information which is backend
independent. Individual runtime backends should subclass these as
required to implement backend-specific features.

.. _MatMPIAIJSetPreallocation: http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatMPIAIJSetPreallocation.html
"""

import numpy as np
import operator
Expand All @@ -50,6 +56,7 @@
PYOP2_COMM = None

def get_mpi_communicator():
"""The MPI Communicator used by PyOP2."""
global PYOP2_COMM
return PYOP2_COMM

Expand All @@ -76,6 +83,8 @@ class Arg(base.Arg):
"""

def halo_exchange_begin(self):
"""Begin halo exchange for the argument if a halo update is required.
Doing halo exchanges only makes sense for :class:`Dat` objects."""
assert self._is_dat, "Doing halo exchanges only makes sense for Dats"
assert not self._in_flight, \
"Halo exchange already in flight for Arg %s" % self
Expand All @@ -85,12 +94,16 @@ def halo_exchange_begin(self):
self.data.halo_exchange_begin()

def halo_exchange_end(self):
"""End halo exchange if it is in flight.
Doing halo exchanges only makes sense for :class:`Dat` objects."""
assert self._is_dat, "Doing halo exchanges only makes sense for Dats"
if self.access in [READ, RW] and self._in_flight:
self._in_flight = False
self.data.halo_exchange_end()

def reduction_begin(self):
"""Begin reduction for the argument if its access is INC, MIN, or MAX.
Doing a reduction only makes sense for :class:`Global` objects."""
assert self._is_global, \
"Doing global reduction only makes sense for Globals"
assert not self._in_flight, \
Expand All @@ -112,6 +125,8 @@ def reduction_begin(self):
PYOP2_COMM.Allreduce(self.data._data, self.data._buf, op=op)

def reduction_end(self):
"""End reduction for the argument if it is in flight.
Doing a reduction only makes sense for :class:`Global` objects."""
assert self._is_global, \
"Doing global reduction only makes sense for Globals"
if self.access is not READ and self._in_flight:
Expand All @@ -135,6 +150,7 @@ def __init__(self, size, name=None, halo=None):

@classmethod
def fromhdf5(cls, f, name):
"""Construct a :class:`Set` from set named ``name`` in HDF5 data ``f``"""
slot = f[name]
size = slot.value.astype(np.int)
shape = slot.shape
Expand Down Expand Up @@ -269,6 +285,7 @@ def __idiv__(self, other):
return self._iop(other, operator.idiv)

def halo_exchange_begin(self):
"""Begin halo exchange."""
halo = self.dataset.halo
if halo is None:
return
Expand All @@ -292,6 +309,7 @@ def halo_exchange_begin(self):
source=source, tag=self._id)

def halo_exchange_end(self):
"""End halo exchange. Waits on MPI recv."""
halo = self.dataset.halo
if halo is None:
return
Expand All @@ -310,6 +328,7 @@ def norm(self):

@classmethod
def fromhdf5(cls, dataset, f, name):
"""Construct a :class:`Dat` from a Dat named ``name`` in HDF5 data ``f``"""
slot = f[name]
data = slot.value
dim = slot.shape[1:]
Expand All @@ -321,6 +340,7 @@ def fromhdf5(cls, dataset, f, name):

@property
def vec(self):
"""PETSc Vec appropriate for this Dat."""
if not hasattr(self, '_vec'):
size = (self.dataset.size * self.cdim, None)
self._vec = PETSc.Vec().createWithArray(self._data, size=size)
Expand All @@ -337,6 +357,7 @@ class Const(base.Const):

@classmethod
def fromhdf5(cls, f, name):
"""Construct a :class:`Const` from const named ``name`` in HDF5 data ``f``"""
slot = f[name]
dim = slot.shape
data = slot.value
Expand All @@ -355,6 +376,7 @@ def _c_handle(self):

@classmethod
def fromhdf5(cls, iterset, dataset, f, name):
"""Construct a :class:`Map` from set named ``name`` in HDF5 data ``f``"""
slot = f[name]
values = slot.value
dim = slot.shape[1:]
Expand Down Expand Up @@ -398,26 +420,48 @@ def __del__(self):

@property
def rowptr(self):
"""Row pointer array of CSR data structure."""
return self._rowptr

@property
def colidx(self):
"""Column indices array of CSR data structure."""
return self._colidx

@property
def nnz(self):
"""Array containing the number of non-zeroes in the various rows of the
diagonal portion of the local submatrix.

This is the same as the parameter `d_nnz` used for preallocation in
PETSc's MatMPIAIJSetPreallocation_."""
return self._d_nnz

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe mention that these things are all the same as PETSc's description and forward to the relevant page in the PETSc docs?


@property
def onnz(self):
"""Array containing the number of non-zeroes in the various rows of the
off-diagonal portion of the local submatrix.

This is the same as the parameter `o_nnz` used for preallocation in
PETSc's MatMPIAIJSetPreallocation_."""
return self._o_nnz

@property
def nz(self):
"""Number of non-zeroes per row in diagonal portion of the local
submatrix.

This is the same as the parameter `d_nz` used for preallocation in
PETSc's MatMPIAIJSetPreallocation_."""
return int(self._d_nz)

@property
def onz(self):
"""Number of non-zeroes per row in off-diagonal portion of the local
submatrix.

This is the same as the parameter o_nz used for preallocation in
PETSc's MatMPIAIJSetPreallocation_."""
return int(self._o_nz)

class Mat(base.Mat):
Expand Down Expand Up @@ -479,6 +523,7 @@ def _assemble(self):

@property
def array(self):
"""Array of non-zero values."""
if not hasattr(self, '_array'):
self._init()
return self._array
Expand All @@ -489,12 +534,14 @@ def values(self):

@property
def handle(self):
"""Petsc4py Mat holding matrix data."""
if self._handle is None:
self._init()
return self._handle

class ParLoop(base.ParLoop):
def compute(self):
"""Executes the kernel over all members of the iteration space."""
raise RuntimeError('Must select a backend')

# FIXME: Eventually (when we have a proper OpenCL solver) this wants to go in
Expand Down