From e2de1dbce333244fb46fbbcd6546006c86fc6ac6 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Tue, 16 May 2017 12:04:22 +0200 Subject: [PATCH 01/18] yask: Change output format --- devito/dle/backends/yask.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/devito/dle/backends/yask.py b/devito/dle/backends/yask.py index 6f0e460f98..0cf627a8c4 100644 --- a/devito/dle/backends/yask.py +++ b/devito/dle/backends/yask.py @@ -60,15 +60,15 @@ def _yaskize(self, state): try: ast = transformer(i.stencil) # Scalar - print(ast.format_simple()) + # print(ast.format_simple()) # AVX2 intrinsics # print soln.format('avx2') # AVX2 intrinsics to file - # import os - # path = os.path.join(os.environ.get('YASK_HOME', '.'), 'src') - # soln.write(os.path.join(path, 'stencil_code.hpp'), 'avx2') + import os + path = os.path.join(os.environ.get('YASK_HOME', '.'), 'src') + soln.write(os.path.join(path, 'stencil_code.hpp'), 'avx2') except: pass From 5792500b6779e4d1c719aa1ec6af8844d23e30b6 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Wed, 17 May 2017 12:18:04 +0200 Subject: [PATCH 02/18] yask: Add yaskarray, subclass of numpy.ndarray --- devito/dle/__init__.py | 1 + devito/dle/backends/__init__.py | 2 +- devito/dle/backends/yask.py | 22 ++++++++++++++++++++++ devito/memory.py | 13 +++++++++++++ tests/test_yask.py | 19 +++++++++++++++++++ 5 files changed, 56 insertions(+), 1 deletion(-) create mode 100644 tests/test_yask.py diff --git a/devito/dle/__init__.py b/devito/dle/__init__.py index 16585c7fad..a5d3d53514 100644 --- a/devito/dle/__init__.py +++ b/devito/dle/__init__.py @@ -1,3 +1,4 @@ from devito.dle.inspection import * # noqa from devito.dle.manipulation import * # noqa from devito.dle.transformer import * # noqa +from devito.dle.backends import yaskarray diff --git a/devito/dle/backends/__init__.py b/devito/dle/backends/__init__.py index e73c80d7b7..73f1615b2d 100644 --- a/devito/dle/backends/__init__.py +++ b/devito/dle/backends/__init__.py @@ -3,4 +3,4 @@ from devito.dle.backends.basic import BasicRewriter # noqa from devito.dle.backends.advanced import (DevitoRewriter, DevitoSpeculativeRewriter, # noqa DevitoCustomRewriter) # noqa -from devito.dle.backends.yask import YaskRewriter # noqa +from devito.dle.backends.yask import * # noqa diff --git a/devito/dle/backends/yask.py b/devito/dle/backends/yask.py index 0cf627a8c4..b3cb70798d 100644 --- a/devito/dle/backends/yask.py +++ b/devito/dle/backends/yask.py @@ -1,3 +1,4 @@ +import numpy as np from sympy import Indexed from devito.dimension import LoweredDimension @@ -14,6 +15,8 @@ except ImportError: yask = None +__all__ = ['YaskRewriter', 'yaskarray'] + class YaskRewriter(BasicRewriter): @@ -117,3 +120,22 @@ def run(expr): raise NotImplementedError return run(expr) + + +class yaskarray(np.ndarray): + + """ + An implementation of a ``numpy.ndarray`` suitable for the YASK storage layout. + + This subclass follows the ``numpy`` rules for subclasses detailed at: :: + + https://docs.scipy.org/doc/numpy/user/basics.subclassing.html + """ + + def __new__(cls, array): + # Input array is an already formed ndarray instance + # We first cast to be our class type + return np.asarray(array).view(cls) + + def __array_finalize__(self, obj): + if obj is None: return diff --git a/devito/memory.py b/devito/memory.py index 15727c6eb3..247d92bc32 100644 --- a/devito/memory.py +++ b/devito/memory.py @@ -10,13 +10,26 @@ from devito.dimension import t from devito.logger import error +from devito.parameters import parameters from devito.tools import convert_dtype_to_ctype class CMemory(object): + def __init__(self, shape, dtype=np.float32, alignment=None): self.ndpointer, self.data_pointer = malloc_aligned(shape, alignment, dtype) + @classmethod + def cast(cls, ndpointer): + """ + Pick a different ndarray type depending on how the backend will access the data. + """ + if parameters['dle']['mode'] == 'yask': + from devito.dle import yaskarray + return yaskarray(ndpointer) + else: + return ndpointer + def __del__(self): free(self.data_pointer) self.data_pointer = None diff --git a/tests/test_yask.py b/tests/test_yask.py new file mode 100644 index 0000000000..dd3fc70833 --- /dev/null +++ b/tests/test_yask.py @@ -0,0 +1,19 @@ +from sympy import Eq + +import pytest + +from devito import DLE_DEFAULT, Operator, DenseData, parameters, x, y +from devito.dle.backends import yaskarray + + +def setup_module(module): + parameters['dle']['mode'] = 'yask' + + +def teardown_module(module): + parameters['dle']['mode'] = DLE_DEFAULT + + +def test_data_type(): + u = DenseData(name='yu', shape=(10, 10), dimensions=(x, y)) + assert type(u.data) == yaskarray From 61090269e33f8eeccd498c45f03b1360c4163e65 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Wed, 17 May 2017 12:47:53 +0200 Subject: [PATCH 03/18] yask: Add basic get/setitem into yaskarray --- devito/dle/backends/yask.py | 10 ++++++++++ tests/test_yask.py | 22 +++++++++++++++++++++- 2 files changed, 31 insertions(+), 1 deletion(-) diff --git a/devito/dle/backends/yask.py b/devito/dle/backends/yask.py index b3cb70798d..5433e92092 100644 --- a/devito/dle/backends/yask.py +++ b/devito/dle/backends/yask.py @@ -127,6 +127,9 @@ class yaskarray(np.ndarray): """ An implementation of a ``numpy.ndarray`` suitable for the YASK storage layout. + WIP: Currently, the YASK storage layout is assumed transposed w.r.t. the + usual row-major format. + This subclass follows the ``numpy`` rules for subclasses detailed at: :: https://docs.scipy.org/doc/numpy/user/basics.subclassing.html @@ -139,3 +142,10 @@ def __new__(cls, array): def __array_finalize__(self, obj): if obj is None: return + + def __getitem__(self, index): + expected_layout = self.transpose() + return super(yaskarray, expected_layout).__getitem__(index) + + def __setitem__(self, index, val): + super(yaskarray, self).__setitem__(index, val) diff --git a/tests/test_yask.py b/tests/test_yask.py index dd3fc70833..9495293242 100644 --- a/tests/test_yask.py +++ b/tests/test_yask.py @@ -1,8 +1,9 @@ +import numpy as np from sympy import Eq import pytest -from devito import DLE_DEFAULT, Operator, DenseData, parameters, x, y +from devito import DLE_DEFAULT, Operator, DenseData, parameters, x, y, z from devito.dle.backends import yaskarray @@ -17,3 +18,22 @@ def teardown_module(module): def test_data_type(): u = DenseData(name='yu', shape=(10, 10), dimensions=(x, y)) assert type(u.data) == yaskarray + + +def test_data_swap(): + u = DenseData(name='yu1D', shape=(10,), dimensions=(x,)) + u.data[1] = 1. + assert u.data[1] == 1. + u = DenseData(name='yu2D', shape=(10, 10), dimensions=(x, y)) + u.data[0, 1] = 1. + assert u.data[1, 0] == 1. + u = DenseData(name='yu3D', shape=(10, 10, 10), dimensions=(x, y, z)) + u.data[0, 1, 1] = 1. + assert u.data[1, 1, 0] == 1. + + +def test_storage_layout(): + u = DenseData(name='yu', shape=(10, 10), dimensions=(x, y)) + op = Operator(Eq(u, 1.), dse='noop', dle='noop') + op.apply(u) + assert np.allclose(u.data, 1) From 9bfca1d3c7096c5e838478e179168fdc88c3afb5 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Mon, 29 May 2017 10:40:06 +0200 Subject: [PATCH 04/18] Leftover fixes after non-trivial rebase --- devito/memory.py | 7 ++++--- tests/test_yask.py | 7 ++++--- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/devito/memory.py b/devito/memory.py index 247d92bc32..f022b40592 100644 --- a/devito/memory.py +++ b/devito/memory.py @@ -10,7 +10,6 @@ from devito.dimension import t from devito.logger import error -from devito.parameters import parameters from devito.tools import convert_dtype_to_ctype @@ -24,8 +23,10 @@ def cast(cls, ndpointer): """ Pick a different ndarray type depending on how the backend will access the data. """ - if parameters['dle']['mode'] == 'yask': - from devito.dle import yaskarray + from devito.parameters import configuration + from devito.dle import yaskarray + + if configuration['dle'] == 'yask': return yaskarray(ndpointer) else: return ndpointer diff --git a/tests/test_yask.py b/tests/test_yask.py index 9495293242..ea92d80e50 100644 --- a/tests/test_yask.py +++ b/tests/test_yask.py @@ -3,16 +3,17 @@ import pytest -from devito import DLE_DEFAULT, Operator, DenseData, parameters, x, y, z +from devito import Operator, DenseData, x, y, z +from devito.parameters import configuration, defaults from devito.dle.backends import yaskarray def setup_module(module): - parameters['dle']['mode'] = 'yask' + configuration['dle'] = 'yask' def teardown_module(module): - parameters['dle']['mode'] = DLE_DEFAULT + configuration['dle'] = defaults['dle'] def test_data_type(): From 047a8a8f69ebaaceffa895160a47b6eea577834e Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Mon, 29 May 2017 10:40:20 +0200 Subject: [PATCH 05/18] yask: Update to new YASK compiler API --- devito/dle/backends/yask.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/devito/dle/backends/yask.py b/devito/dle/backends/yask.py index 5433e92092..f7ce91b0b8 100644 --- a/devito/dle/backends/yask.py +++ b/devito/dle/backends/yask.py @@ -10,8 +10,8 @@ try: import yask_compiler as yask # Factories to interact with YASK - cfac = yask.yask_compiler_factory() - fac = yask.node_factory() + cfac = yask.yc_factory() + fac = yask.yc_node_factory() except ImportError: yask = None @@ -43,7 +43,7 @@ def _yaskize(self, state): candidate = tree[-1] # Set up the YASK solution - soln = cfac.new_stencil_solution("solution") + soln = cfac.new_solution("solution") # Set up the YASK grids grids = FindSymbols(mode='symbolics').visit(candidate) @@ -55,13 +55,17 @@ def _yaskize(self, state): soln.set_fold_len('y', 1) soln.set_fold_len('z', 8) + # Set necessary run-time parameters + soln.set_step_dim("t") + soln.set_elem_bytes(4) + # Perform the translation on an expression basis transformer = sympy2yask(grids) expressions = [e for e in candidate.nodes if e.is_Expression] # yaskASTs = [transformer(e.stencil) for e in expressions] for i in expressions: try: - ast = transformer(i.stencil) + ast = transformer(i.expr) # Scalar # print(ast.format_simple()) From 7d8ae2544a6009c7234e1797e5a24f621e2cb3dd Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Wed, 31 May 2017 16:59:14 +0200 Subject: [PATCH 06/18] yask: Initial DSE support; update to new API --- devito/dle/backends/yask.py | 66 +++++++++++++++++++++++-------------- 1 file changed, 41 insertions(+), 25 deletions(-) diff --git a/devito/dle/backends/yask.py b/devito/dle/backends/yask.py index f7ce91b0b8..5d3b228e6a 100644 --- a/devito/dle/backends/yask.py +++ b/devito/dle/backends/yask.py @@ -1,10 +1,12 @@ +import os + import numpy as np from sympy import Indexed from devito.dimension import LoweredDimension from devito.dle import retrieve_iteration_tree from devito.dle.backends import BasicRewriter, dle_pass -from devito.logger import dle_warning +from devito.logger import dle, dle_warning from devito.visitors import FindSymbols try: @@ -43,41 +45,47 @@ def _yaskize(self, state): candidate = tree[-1] # Set up the YASK solution - soln = cfac.new_solution("solution") + soln = cfac.new_solution("devito-test-solution") # Set up the YASK grids grids = FindSymbols(mode='symbolics').visit(candidate) grids = {g.name: soln.new_grid(g.name, *[str(i) for i in g.indices]) for g in grids} + # Perform the translation on an expression basis + transform = sympy2yask(grids) + expressions = [e for e in candidate.nodes if e.is_Expression] + try: + for i in expressions: + transform(i.expr) + dle("Converted %s into YASK format", str(i.expr)) + except: + dle_warning("Cannot convert %s into YASK format", str(i.expr)) + continue + + # Print some useful information to screen about the YASK conversion + dle("Solution '" + soln.get_name() + "' contains " + + str(soln.get_num_grids()) + " grid(s), and " + + str(soln.get_num_equations()) + " equation(s).") + + # Provide stuff to YASK-land + # ========================== + # Scalar: print(ast.format_simple()) + # AVX2 intrinsics: print soln.format('avx2') + # AVX2 intrinsics to file (active now) + path = os.path.join(os.environ.get('YASK_HOME', '.'), + 'src', 'kernel', 'gen') + soln.write(os.path.join(path, 'yask_stencil_code.hpp'), 'avx2', True) + + # Set kernel parameters + # ===================== # Vector folding API usage example soln.set_fold_len('x', 1) soln.set_fold_len('y', 1) soln.set_fold_len('z', 8) - # Set necessary run-time parameters soln.set_step_dim("t") - soln.set_elem_bytes(4) - - # Perform the translation on an expression basis - transformer = sympy2yask(grids) - expressions = [e for e in candidate.nodes if e.is_Expression] - # yaskASTs = [transformer(e.stencil) for e in expressions] - for i in expressions: - try: - ast = transformer(i.expr) - # Scalar - # print(ast.format_simple()) - - # AVX2 intrinsics - # print soln.format('avx2') - - # AVX2 intrinsics to file - import os - path = os.path.join(os.environ.get('YASK_HOME', '.'), 'src') - soln.write(os.path.join(path, 'stencil_code.hpp'), 'avx2') - except: - pass + soln.set_element_bytes(4) dle_warning("Falling back to basic DLE optimizations...") @@ -91,6 +99,7 @@ class sympy2yask(object): def __init__(self, grids): self.grids = grids + self.mapper = {} def __call__(self, expr): @@ -103,6 +112,9 @@ def run(expr): return fac.new_const_number_node(int(expr)) elif expr.is_Float: return fac.new_const_number_node(float(expr)) + elif expr.is_Symbol: + assert expr in self.mapper + return self.mapper[expr] elif isinstance(expr, Indexed): function = expr.base.function assert function.name in self.grids @@ -118,7 +130,11 @@ def run(expr): if num == 1: return fac.new_divide_node(run(num), run(den)) elif expr.is_Equality: - return fac.new_equation_node(*[run(i) for i in expr.args]) + if expr.lhs.is_Symbol: + assert expr.lhs not in self.mapper + self.mapper[expr.lhs] = run(expr.rhs) + else: + return fac.new_equation_node(*[run(i) for i in expr.args]) else: dle_warning("Missing handler in Devito-YASK translation") raise NotImplementedError From c362900893971fdd02eb5ad2fa4976facb728d96 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Tue, 6 Jun 2017 09:35:40 +0200 Subject: [PATCH 07/18] yask: Strengthen input checks --- devito/dle/__init__.py | 2 +- devito/dle/backends/yask.py | 30 ++++++++++++++++++++---------- tests/test_yask.py | 2 +- 3 files changed, 22 insertions(+), 12 deletions(-) diff --git a/devito/dle/__init__.py b/devito/dle/__init__.py index a5d3d53514..1002b7ec53 100644 --- a/devito/dle/__init__.py +++ b/devito/dle/__init__.py @@ -1,4 +1,4 @@ from devito.dle.inspection import * # noqa from devito.dle.manipulation import * # noqa from devito.dle.transformer import * # noqa -from devito.dle.backends import yaskarray +from devito.dle.backends import yaskarray # noqa diff --git a/devito/dle/backends/yask.py b/devito/dle/backends/yask.py index 5d3b228e6a..e630e7d1a9 100644 --- a/devito/dle/backends/yask.py +++ b/devito/dle/backends/yask.py @@ -9,22 +9,33 @@ from devito.logger import dle, dle_warning from devito.visitors import FindSymbols +__all__ = ['YaskRewriter', 'yaskarray'] + + +# Init + try: - import yask_compiler as yask + import yask_compiler as yc # Factories to interact with YASK - cfac = yask.yc_factory() - fac = yask.yc_node_factory() + cfac = yc.yc_factory() + fac = yc.yc_node_factory() except ImportError: - yask = None + yc, emsg = None, "[Python bindings]" -__all__ = ['YaskRewriter', 'yaskarray'] +try: + # Set directory for generated code + path = os.path.join(os.environ['YASK_HOME'], 'src', 'kernel', 'gen') + if not os.path.exists(path): + os.makedirs(path) +except KeyError: + yc, emsg = None, "[Missing YASK_HOME]" class YaskRewriter(BasicRewriter): def _pipeline(self, state): - if yask is None: - dle_warning("Cannot find YASK. Skipping DLE optimizations...") + if yc is None: + dle_warning("Cannot find YASK %s. Skipping DLE optimizations..." % emsg) super(YaskRewriter, self)._pipeline(state) return self._avoid_denormals(state) @@ -73,8 +84,6 @@ def _yaskize(self, state): # Scalar: print(ast.format_simple()) # AVX2 intrinsics: print soln.format('avx2') # AVX2 intrinsics to file (active now) - path = os.path.join(os.environ.get('YASK_HOME', '.'), - 'src', 'kernel', 'gen') soln.write(os.path.join(path, 'yask_stencil_code.hpp'), 'avx2', True) # Set kernel parameters @@ -161,7 +170,8 @@ def __new__(cls, array): return np.asarray(array).view(cls) def __array_finalize__(self, obj): - if obj is None: return + if obj is None: + return def __getitem__(self, index): expected_layout = self.transpose() diff --git a/tests/test_yask.py b/tests/test_yask.py index ea92d80e50..502ac00f43 100644 --- a/tests/test_yask.py +++ b/tests/test_yask.py @@ -1,7 +1,7 @@ import numpy as np from sympy import Eq -import pytest +import pytest # noqa from devito import Operator, DenseData, x, y, z from devito.parameters import configuration, defaults From 4a5753d32ab15f9188fcf02c9b5ef7e88f5681f4 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Tue, 6 Jun 2017 12:32:53 +0200 Subject: [PATCH 08/18] compiler: Add utility to invoke Makefiles --- devito/compiler.py | 29 +++++++++++++++++++++++++++++ devito/exceptions.py | 4 ++++ devito/tools.py | 20 ++++++++++++++++++++ 3 files changed, 53 insertions(+) diff --git a/devito/compiler.py b/devito/compiler.py index 3b2cc634d3..8fb61210be 100644 --- a/devito/compiler.py +++ b/devito/compiler.py @@ -4,12 +4,15 @@ from tempfile import gettempdir from time import time from sys import platform +import subprocess import numpy.ctypeslib as npct from codepy.jit import extension_file_from_string from codepy.toolchain import GCCToolchain +from devito.exceptions import CompilationError from devito.logger import log +from devito.tools import change_directory __all__ = ['get_tmp_dir', 'set_compiler', 'jit_compile', 'load', 'GNUCompiler'] @@ -275,3 +278,29 @@ def jit_compile(ccode, compiler=GNUCompiler): log("%s: compiled %s [%.2f s]" % (compiler, src_file, toc-tic)) return basename + + +def make(loc, args): + """ + Invoke ``make`` command from within ``loc`` with arguments ``args``. + """ + hash_key = sha1(loc + str(args).encode()).hexdigest() + logfile = path.join(get_tmp_dir(), "%s.log" % hash_key) + errfile = path.join(get_tmp_dir(), "%s.err" % hash_key) + + with change_directory(loc): + with open(logfile, "w") as log: + with open(errfile, "w") as err: + + command = ['make'] + args + log.write("Compilation command:\n") + log.write(" ".join(command)) + log.write("\n\n") + try: + subprocess.check_call(command, stderr=err, stdout=log) + except subprocess.CalledProcessError as e: + raise CompilationError('Command "%s" return error status %d. ' + 'Unable to compile code.\n' + 'Compile log in %s\n' + 'Compile errors in %s\n' % + (e.cmd, e.returncode, logfile, errfile)) diff --git a/devito/exceptions.py b/devito/exceptions.py index 6c3f814c55..abbe5e9ad0 100644 --- a/devito/exceptions.py +++ b/devito/exceptions.py @@ -1,3 +1,7 @@ +class CompilationError(Exception): + pass + + class InvalidArgument(Exception): pass diff --git a/devito/tools.py b/devito/tools.py index 91e40d57b4..9295af0e59 100644 --- a/devito/tools.py +++ b/devito/tools.py @@ -1,3 +1,4 @@ +import os import ctypes from collections import Callable, Iterable, OrderedDict try: @@ -195,3 +196,22 @@ def copy(self): def __copy__(self): return type(self)(self.default_factory, self) + + +class change_directory(object): + """ + Context manager for changing the current working directory. + + Adapted from: :: + + https://stackoverflow.com/questions/431684/how-do-i-cd-in-python/ + """ + def __init__(self, newPath): + self.newPath = os.path.expanduser(newPath) + + def __enter__(self): + self.savedPath = os.getcwd() + os.chdir(self.newPath) + + def __exit__(self, etype, value, traceback): + os.chdir(self.savedPath) From 8c8c2c211be0cc3825b3ec6fff90cb119715ec52 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Thu, 8 Jun 2017 12:23:24 +0200 Subject: [PATCH 09/18] yask: Introduce YaskState concept This involves a proper init routine Also update to newest YASK API --- devito/dle/__init__.py | 2 +- devito/dle/backends/yask.py | 241 +++++++++++++++++++++++++++--------- devito/interfaces.py | 25 +++- devito/memory.py | 13 -- tests/test_yask.py | 4 +- 5 files changed, 204 insertions(+), 81 deletions(-) diff --git a/devito/dle/__init__.py b/devito/dle/__init__.py index 1002b7ec53..f7f213c87f 100644 --- a/devito/dle/__init__.py +++ b/devito/dle/__init__.py @@ -1,4 +1,4 @@ from devito.dle.inspection import * # noqa from devito.dle.manipulation import * # noqa from devito.dle.transformer import * # noqa -from devito.dle.backends import yaskarray # noqa +from devito.dle.backends import init, make_grid # noqa diff --git a/devito/dle/backends/yask.py b/devito/dle/backends/yask.py index e630e7d1a9..f8846f16dc 100644 --- a/devito/dle/backends/yask.py +++ b/devito/dle/backends/yask.py @@ -1,43 +1,101 @@ import os +import sys import numpy as np from sympy import Indexed +from devito.compiler import make from devito.dimension import LoweredDimension from devito.dle import retrieve_iteration_tree from devito.dle.backends import BasicRewriter, dle_pass -from devito.logger import dle, dle_warning +from devito.exceptions import CompilationError +from devito.logger import dle, dle_warning, error from devito.visitors import FindSymbols -__all__ = ['YaskRewriter', 'yaskarray'] +__all__ = ['YaskRewriter', 'init', 'make_grid'] -# Init +YASK = None +"""Global state for generation of YASK kernels.""" -try: - import yask_compiler as yc - # Factories to interact with YASK - cfac = yc.yc_factory() - fac = yc.yc_node_factory() -except ImportError: - yc, emsg = None, "[Python bindings]" -try: - # Set directory for generated code - path = os.path.join(os.environ['YASK_HOME'], 'src', 'kernel', 'gen') - if not os.path.exists(path): - os.makedirs(path) -except KeyError: - yc, emsg = None, "[Missing YASK_HOME]" +class YaskState(object): + + def __init__(self, cfac, nfac, path, env, settings, hook_soln): + """ + Global state to interact with YASK. + """ + self.cfac = cfac # YASK compiler factory, to create Solutions + self.nfac = nfac # YASK node factory, to create ASTs + self.path = path # Generated code dump directory + self.env = env # Global environment (e.g., MPI) + self.settings = settings # Dimensions, grid sizes, etc. + self.hook_soln = hook_soln # "Fake" solution to track YASK grids + + @property + def dimensions(self): + return self.hook_soln.get_domain_dim_names() + + @property + def grids(self): + mapper = {} + for i in range(self.hook_soln.get_num_grids()): + grid = self.hook_soln.get_grid(i) + mapper[grid.get_name()] = grid + return mapper + + def setdefault(self, name, vals=0.0): + """ + Add and return a new grid ``name``. If a grid ``name`` already exists, + then return it without performing any other actions. + """ + grids = self.grids + if name in grids: + return grids[name] + else: + # new_grid() also modifies the /hook_soln/ state + grid = self.hook_soln.new_grid(name, *self.dimensions) + # Allocate memory + self.hook_soln.prepare_solution() + # Initialization + grid.set_all_elements(vals) + # TODO : return YaskGrid (subclass of NumPy array) + return grid + + +class YaskGrid(np.ndarray): + + """ + An implementation of a ``numpy.ndarray`` suitable for the YASK storage layout. + + WIP: Currently, the YASK storage layout is assumed transposed w.r.t. the + usual row-major format. + + This subclass follows the ``numpy`` rules for subclasses detailed at: :: + + https://docs.scipy.org/doc/numpy/user/basics.subclassing.html + """ + + def __new__(cls, array): + # Input array is an already formed ndarray instance + # We first cast to be our class type + return np.asarray(array).view(cls) + + def __array_finalize__(self, obj): + if obj is None: + return + + def __getitem__(self, index): + expected_layout = self.transpose() + return super(YaskGrid, expected_layout).__getitem__(index) + + def __setitem__(self, index, val): + super(YaskGrid, self).__setitem__(index, val) class YaskRewriter(BasicRewriter): def _pipeline(self, state): - if yc is None: - dle_warning("Cannot find YASK %s. Skipping DLE optimizations..." % emsg) - super(YaskRewriter, self)._pipeline(state) - return self._avoid_denormals(state) self._yaskize(state) self._create_elemental_functions(state) @@ -56,12 +114,11 @@ def _yaskize(self, state): candidate = tree[-1] # Set up the YASK solution - soln = cfac.new_solution("devito-test-solution") + soln = YASK.cfac.new_solution("devito-test-solution") # Set up the YASK grids grids = FindSymbols(mode='symbolics').visit(candidate) - grids = {g.name: soln.new_grid(g.name, *[str(i) for i in g.indices]) - for g in grids} + grids = [YASK.setdefault(i.name) for i in grids] # Perform the translation on an expression basis transform = sympy2yask(grids) @@ -81,17 +138,11 @@ def _yaskize(self, state): # Provide stuff to YASK-land # ========================== - # Scalar: print(ast.format_simple()) + # Scalar: print(ast.format_simple()) # AVX2 intrinsics: print soln.format('avx2') # AVX2 intrinsics to file (active now) - soln.write(os.path.join(path, 'yask_stencil_code.hpp'), 'avx2', True) - - # Set kernel parameters - # ===================== - # Vector folding API usage example - soln.set_fold_len('x', 1) - soln.set_fold_len('y', 1) - soln.set_fold_len('z', 8) + soln.write(os.path.join(YASK.path, 'yask_stencil_code.hpp'), 'avx2', True) + # Set necessary run-time parameters soln.set_step_dim("t") soln.set_element_bytes(4) @@ -118,9 +169,9 @@ def nary2binary(args, op): def run(expr): if expr.is_Integer: - return fac.new_const_number_node(int(expr)) + return YASK.nfac.new_const_number_node(int(expr)) elif expr.is_Float: - return fac.new_const_number_node(float(expr)) + return YASK.nfac.new_const_number_node(float(expr)) elif expr.is_Symbol: assert expr in self.mapper return self.mapper[expr] @@ -131,19 +182,19 @@ def run(expr): for i, j in zip(expr.indices, function.indices)] return self.grids[function.name].new_relative_grid_point(*indices) elif expr.is_Add: - return nary2binary(expr.args, fac.new_add_node) + return nary2binary(expr.args, YASK.nfac.new_add_node) elif expr.is_Mul: - return nary2binary(expr.args, fac.new_multiply_node) + return nary2binary(expr.args, YASK.nfac.new_multiply_node) elif expr.is_Pow: num, den = expr.as_numer_denom() if num == 1: - return fac.new_divide_node(run(num), run(den)) + return YASK.nfac.new_divide_node(run(num), run(den)) elif expr.is_Equality: if expr.lhs.is_Symbol: assert expr.lhs not in self.mapper self.mapper[expr.lhs] = run(expr.rhs) else: - return fac.new_equation_node(*[run(i) for i in expr.args]) + return YASK.nfac.new_equation_node(*[run(i) for i in expr.args]) else: dle_warning("Missing handler in Devito-YASK translation") raise NotImplementedError @@ -151,31 +202,101 @@ def run(expr): return run(expr) -class yaskarray(np.ndarray): +# YASK interface +def init(dimensions, shape, dtype, architecture='hsw', isa='avx2'): """ - An implementation of a ``numpy.ndarray`` suitable for the YASK storage layout. - - WIP: Currently, the YASK storage layout is assumed transposed w.r.t. the - usual row-major format. + To be called prior to any YASK-related operation. - This subclass follows the ``numpy`` rules for subclasses detailed at: :: + A new bootstrap is required wheneven any of the following change: :: - https://docs.scipy.org/doc/numpy/user/basics.subclassing.html + * YASK version + * Target architecture (``architecture`` param) + * Floating-point precision (``dtype`` param) + * Domain dimensions (``dimensions`` param) + * Folding + * Grid memory layout scheme """ + global YASK + + if YASK is not None: + return + + dle("Initializing YASK...") + + try: + import yask_compiler as yc + # YASK compiler factories + cfac = yc.yc_factory() + nfac = yc.yc_node_factory() + except ImportError: + _force_exit("Python YASK compiler bindings") + + try: + # Set directory for generated code + path = os.path.join(os.environ['YASK_HOME'], 'src', 'kernel', 'gen') + if not os.path.exists(path): + os.makedirs(path) + except KeyError: + _force_exit("Missing YASK_HOME") + + # Create a new stencil solution + soln = cfac.new_solution("Hook") + soln.set_step_dim("t") + soln.set_domain_dims(*[str(i) for i in dimensions]) # TODO: YASK only accepts x,y,z + + # Number of bytes in each FP value + soln.set_element_bytes(dtype().itemsize) + + # Generate YASK output + soln.write(os.path.join(path, 'yask_stencil_code.hpp'), isa, True) + + # Build YASK output, and load the corresponding YASK kernel + try: + make(os.environ['YASK_HOME'], + ['-j', 'stencil=Hook', 'arch=%s' % architecture, 'yk-api']) + except CompilationError: + _force_exit("Hook solution compilation") + try: + import yask_kernel as yk + except ImportError: + _force_exit("Python YASK kernel bindings") + + # YASK Hook kernel factory + kfac = yk.yk_factory() + + # Initalize MPI, etc + env = kfac.new_env() + + # Set global settings and create hook solution + settings = kfac.new_settings() + hook_soln = kfac.new_solution(env, settings) + for dm, ds in zip(hook_soln.get_domain_dim_names(), shape): + # Set domain size in each dim. + settings.set_rank_domain_size(dm, ds) + # Set block size to 64 in z dim and 32 in other dims. + settings.set_block_size(dm, min(64 if dm == "z" else 32, ds)) + + # Simple rank configuration in 1st dim only. # TODO Improve me + settings.set_num_ranks(hook_soln.get_domain_dim_name(0), env.get_num_ranks()) + + # Finish off by initializing YASK + YASK = YaskState(cfac, nfac, path, env, settings, hook_soln) + + dle("YASK backend successfully initialized!") + + +def make_grid(name, shape, dimensions, dtype): + """ + Create a new YASK Grid and attach it to a "fake" solution. + """ + init(dimensions, shape, dtype) + return YASK.setdefault(name) - def __new__(cls, array): - # Input array is an already formed ndarray instance - # We first cast to be our class type - return np.asarray(array).view(cls) - - def __array_finalize__(self, obj): - if obj is None: - return - - def __getitem__(self, index): - expected_layout = self.transpose() - return super(yaskarray, expected_layout).__getitem__(index) - def __setitem__(self, index, val): - super(yaskarray, self).__setitem__(index, val) +def _force_exit(emsg): + """ + Handle fatal errors. + """ + error("Couldn't startup YASK [%s]. Exiting..." % emsg) + sys.exit(0) diff --git a/devito/interfaces.py b/devito/interfaces.py index 7454b247e0..1d3a6e3da2 100644 --- a/devito/interfaces.py +++ b/devito/interfaces.py @@ -302,12 +302,27 @@ def _indices(cls, **kwargs): def _allocate_memory(self): """Allocate memory in terms of numpy ndarrays.""" - debug("Allocating memory for %s (%s)" % (self.name, str(self.shape))) - self._data_object = CMemory(self.shape, dtype=self.dtype) - if self.numa: - first_touch(self) + from devito.parameters import configuration + + if configuration['dle'] == 'yask': + from devito.dle import make_grid + self._data = make_grid(self.name, self.shape, self.indices, self.dtype) + + # self._data = self._data_object.ndpointer + # TODO: DELETE THE BELOW + debug("Allocating memory for %s (%s)" % (self.name, str(self.shape))) + self._data_object = CMemory(self.shape, dtype=self.dtype) + if self.numa: + first_touch(self) + else: + self.data.fill(0) else: - self.data.fill(0) + debug("Allocating memory for %s (%s)" % (self.name, str(self.shape))) + self._data_object = CMemory(self.shape, dtype=self.dtype) + if self.numa: + first_touch(self) + else: + self.data.fill(0) @property def data(self): diff --git a/devito/memory.py b/devito/memory.py index f022b40592..209ac10bb2 100644 --- a/devito/memory.py +++ b/devito/memory.py @@ -18,19 +18,6 @@ class CMemory(object): def __init__(self, shape, dtype=np.float32, alignment=None): self.ndpointer, self.data_pointer = malloc_aligned(shape, alignment, dtype) - @classmethod - def cast(cls, ndpointer): - """ - Pick a different ndarray type depending on how the backend will access the data. - """ - from devito.parameters import configuration - from devito.dle import yaskarray - - if configuration['dle'] == 'yask': - return yaskarray(ndpointer) - else: - return ndpointer - def __del__(self): free(self.data_pointer) self.data_pointer = None diff --git a/tests/test_yask.py b/tests/test_yask.py index 502ac00f43..a8d48577f1 100644 --- a/tests/test_yask.py +++ b/tests/test_yask.py @@ -5,7 +5,7 @@ from devito import Operator, DenseData, x, y, z from devito.parameters import configuration, defaults -from devito.dle.backends import yaskarray +from devito.dle.backends import YaskGrid def setup_module(module): @@ -18,7 +18,7 @@ def teardown_module(module): def test_data_type(): u = DenseData(name='yu', shape=(10, 10), dimensions=(x, y)) - assert type(u.data) == yaskarray + assert type(u.data) == YaskGrid def test_data_swap(): From a5158093a7e9536feff72ade344dcb76ae907856 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Thu, 8 Jun 2017 14:51:56 +0200 Subject: [PATCH 10/18] yask: Initial support for YaskGrid [Unfinished] --- devito/dle/backends/yask.py | 145 +++++++++++++++++++++++++----------- devito/interfaces.py | 9 --- tests/test_yask.py | 24 ++++-- 3 files changed, 118 insertions(+), 60 deletions(-) diff --git a/devito/dle/backends/yask.py b/devito/dle/backends/yask.py index f8846f16dc..b06c54af5f 100644 --- a/devito/dle/backends/yask.py +++ b/devito/dle/backends/yask.py @@ -1,5 +1,4 @@ import os -import sys import numpy as np from sympy import Indexed @@ -8,9 +7,10 @@ from devito.dimension import LoweredDimension from devito.dle import retrieve_iteration_tree from devito.dle.backends import BasicRewriter, dle_pass -from devito.exceptions import CompilationError -from devito.logger import dle, dle_warning, error +from devito.exceptions import CompilationError, DLEException +from devito.logger import debug, dle, dle_warning, error from devito.visitors import FindSymbols +from devito.tools import as_tuple __all__ = ['YaskRewriter', 'init', 'make_grid'] @@ -21,16 +21,25 @@ class YaskState(object): - def __init__(self, cfac, nfac, path, env, settings, hook_soln): + def __init__(self, cfac, nfac, path, env, shape, dtype, hook_soln): """ Global state to interact with YASK. + + :param cfac: YASK compiler factory, to create Solutions. + :param nfac: YASK node factory, to create ASTs. + :param path: Generated code dump directory. + :param env: Global environment (e.g., MPI). + :param shape: Domain size along each dimension. + :param dtype: The data type used in kernels, as a NumPy dtype. + :param hook_soln: "Fake" solution to track YASK grids. """ - self.cfac = cfac # YASK compiler factory, to create Solutions - self.nfac = nfac # YASK node factory, to create ASTs - self.path = path # Generated code dump directory - self.env = env # Global environment (e.g., MPI) - self.settings = settings # Dimensions, grid sizes, etc. - self.hook_soln = hook_soln # "Fake" solution to track YASK grids + self.cfac = cfac + self.nfac = nfac + self.path = path + self.env = env + self.shape = shape + self.dtype = dtype + self.hook_soln = hook_soln @property def dimensions(self): @@ -44,7 +53,7 @@ def grids(self): mapper[grid.get_name()] = grid return mapper - def setdefault(self, name, vals=0.0): + def setdefault(self, name, vals=None): """ Add and return a new grid ``name``. If a grid ``name`` already exists, then return it without performing any other actions. @@ -57,40 +66,85 @@ def setdefault(self, name, vals=0.0): grid = self.hook_soln.new_grid(name, *self.dimensions) # Allocate memory self.hook_soln.prepare_solution() - # Initialization - grid.set_all_elements(vals) - # TODO : return YaskGrid (subclass of NumPy array) - return grid - + # TODO : return YaskGrid "wrapping" the allocated data + return YaskGrid(name, self.shape, self.dtype, grid, vals) -class YaskGrid(np.ndarray): - """ - An implementation of a ``numpy.ndarray`` suitable for the YASK storage layout. +class YaskGrid(object): - WIP: Currently, the YASK storage layout is assumed transposed w.r.t. the - usual row-major format. + def __init__(self, name, shape, dtype, grid, val=None): + self.name = name + self.shape = shape + self.dtype = dtype + self.grid = grid - This subclass follows the ``numpy`` rules for subclasses detailed at: :: - - https://docs.scipy.org/doc/numpy/user/basics.subclassing.html - """ + # Always init the grid, at least with 0.0 + self[:] = 0.0 if val is None else val - def __new__(cls, array): - # Input array is an already formed ndarray instance - # We first cast to be our class type - return np.asarray(array).view(cls) + def _convert_multislice(self, multislice): + start = [] + stop = [] + for i, idx in enumerate(multislice): + if isinstance(idx, slice): + start.append(idx.start or 0) + stop.append((idx.stop or self.shape[i]) - 1) + else: + start.append(idx) + stop.append(idx) + shape = [j - i + 1 for i, j in zip(start, stop)] + return start, stop, shape - def __array_finalize__(self, obj): - if obj is None: - return + def __repr__(self): + # TODO: need to return a VIEW + return repr(self[:]) def __getitem__(self, index): - expected_layout = self.transpose() - return super(YaskGrid, expected_layout).__getitem__(index) + if isinstance(index, int) and len(self.shape) > 1 and self.shape != val.shape: + _force_exit("Retrieval from a YaskGrid, non-matching shapes.") + if isinstance(index, slice): + debug("YaskGrid: Getting whole array or block via single slice") + if index.step is not None: + _force_exit("Retrieval from a YaskGrid, unsupported stepping != 1.") + start = [index.start or 0 for j in self.shape] + stop = [(index.stop or j) - 1 for j in self.shape] + shape = [j - i + 1 for i, j in zip(start, stop)] + # TODO: using empty ndarray for now, will need to return views... + output = np.empty(shape, self.dtype, 'C'); + self.grid.get_elements_in_slice(output.data, start, stop) + return output + elif all(isinstance(i, int) for i in as_tuple(index)): + debug("YaskGrid: Getting single entry") + return self.grid.get_element(*as_tuple(index)) + else: + debug("YaskGrid: Getting whole array or block via multiple slices/indices") + start, stop, shape = self._convert_multislice(index) + output = np.empty(shape, self.dtype, 'C'); + self.grid.get_elements_in_slice(output.data, start, stop) + return output def __setitem__(self, index, val): - super(YaskGrid, self).__setitem__(index, val) + # TODO: ATM, no MPI support. + if isinstance(index, int) and len(self.shape) > 1 and self.shape != val.shape: + _force_exit("Insertion into a YaskGrid, non-matching shapes.") + if isinstance(index, slice): + debug("YaskGrid: Setting whole array or block via single slice") + if isinstance(val, np.ndarray): + start = [index.start or 0 for j in self.shape] + stop = [(index.stop or j) - 1 for j in self.shape] + self.grid.set_elements_in_slice(val, start, stop) + else: + self.grid.set_all_elements_same(val) + elif all(isinstance(i, int) for i in as_tuple(index)): + debug("YaskGrid: Setting single entry") + self.grid.set_element(val, *as_tuple(index)) + else: + debug("YaskGrid: Setting whole array or block via multiple slices/indices") + start, stop, _ = self._convert_multislice(index) + if isinstance(val, np.ndarray): + self.grid.set_elements_in_slice(val, start, stop) + else: + # TODO: NEED set_elements_in_slice acceptign single value + assert False, "STIIL WIP" class YaskRewriter(BasicRewriter): @@ -268,20 +322,22 @@ def init(dimensions, shape, dtype, architecture='hsw', isa='avx2'): # Initalize MPI, etc env = kfac.new_env() - # Set global settings and create hook solution - settings = kfac.new_settings() - hook_soln = kfac.new_solution(env, settings) + # Create hook solution + hook_soln = kfac.new_solution(env) for dm, ds in zip(hook_soln.get_domain_dim_names(), shape): # Set domain size in each dim. - settings.set_rank_domain_size(dm, ds) + hook_soln.set_rank_domain_size(dm, ds) + # TODO: Add something like: hook_soln.set_min_pad_size(dm, 16) # Set block size to 64 in z dim and 32 in other dims. - settings.set_block_size(dm, min(64 if dm == "z" else 32, ds)) + hook_soln.set_block_size(dm, min(64 if dm == "z" else 32, ds)) - # Simple rank configuration in 1st dim only. # TODO Improve me - settings.set_num_ranks(hook_soln.get_domain_dim_name(0), env.get_num_ranks()) + # Simple rank configuration in 1st dim only. + # In production runs, the ranks would be distributed along all domain dimensions. + # TODO Improve me + hook_soln.set_num_ranks(hook_soln.get_domain_dim_name(0), env.get_num_ranks()) # Finish off by initializing YASK - YASK = YaskState(cfac, nfac, path, env, settings, hook_soln) + YASK = YaskState(cfac, nfac, path, env, shape, dtype, hook_soln) dle("YASK backend successfully initialized!") @@ -298,5 +354,4 @@ def _force_exit(emsg): """ Handle fatal errors. """ - error("Couldn't startup YASK [%s]. Exiting..." % emsg) - sys.exit(0) + raise DLEException("Couldn't startup YASK [%s]. Exiting..." % emsg) diff --git a/devito/interfaces.py b/devito/interfaces.py index 1d3a6e3da2..f2063e32c1 100644 --- a/devito/interfaces.py +++ b/devito/interfaces.py @@ -307,15 +307,6 @@ def _allocate_memory(self): if configuration['dle'] == 'yask': from devito.dle import make_grid self._data = make_grid(self.name, self.shape, self.indices, self.dtype) - - # self._data = self._data_object.ndpointer - # TODO: DELETE THE BELOW - debug("Allocating memory for %s (%s)" % (self.name, str(self.shape))) - self._data_object = CMemory(self.shape, dtype=self.dtype) - if self.numa: - first_touch(self) - else: - self.data.fill(0) else: debug("Allocating memory for %s (%s)" % (self.name, str(self.shape))) self._data_object = CMemory(self.shape, dtype=self.dtype) diff --git a/tests/test_yask.py b/tests/test_yask.py index a8d48577f1..86c51bc4e8 100644 --- a/tests/test_yask.py +++ b/tests/test_yask.py @@ -5,7 +5,7 @@ from devito import Operator, DenseData, x, y, z from devito.parameters import configuration, defaults -from devito.dle.backends import YaskGrid +from devito.dle.backends.yask import YaskGrid def setup_module(module): @@ -21,16 +21,28 @@ def test_data_type(): assert type(u.data) == YaskGrid -def test_data_swap(): +@pytest.mark.xfail(reason="YASK always seems to use 3D grids") +def test_data_movement_1D(): u = DenseData(name='yu1D', shape=(10,), dimensions=(x,)) u.data[1] = 1. + assert u.data[0] == 0. assert u.data[1] == 1. - u = DenseData(name='yu2D', shape=(10, 10), dimensions=(x, y)) - u.data[0, 1] = 1. - assert u.data[1, 0] == 1. + assert all(i == 0 for i in u.data[2:]) + + +@pytest.mark.xfail(reason="FAIL on block insertion") +def test_data_movement_nD(): u = DenseData(name='yu3D', shape=(10, 10, 10), dimensions=(x, y, z)) u.data[0, 1, 1] = 1. - assert u.data[1, 1, 0] == 1. + assert u.data[0, 0, 0] == 0. + assert u.data[0, 1, 1] == 1. + assert np.all(u.data[:] == u.data[:,:,:]) + + # Test block insertion + block = np.ndarray(shape=(1, 10, 1)) + block.fill(5.) + u.data[5,:,5] = block + assert np.all(u.data[5,:,5] == block) def test_storage_layout(): From 9375196fc9eea890835a09176d02c4fdf9b23caf Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Mon, 19 Jun 2017 17:18:27 +0200 Subject: [PATCH 11/18] yask: Initial support for YaskGrid [Unfinished] Support for negative indices --- devito/dle/backends/yask.py | 72 +++++++++++++++++++++---------------- tests/test_yask.py | 26 ++++++++++---- 2 files changed, 61 insertions(+), 37 deletions(-) diff --git a/devito/dle/backends/yask.py b/devito/dle/backends/yask.py index b06c54af5f..57dcb805d8 100644 --- a/devito/dle/backends/yask.py +++ b/devito/dle/backends/yask.py @@ -66,7 +66,6 @@ def setdefault(self, name, vals=None): grid = self.hook_soln.new_grid(name, *self.dimensions) # Allocate memory self.hook_soln.prepare_solution() - # TODO : return YaskGrid "wrapping" the allocated data return YaskGrid(name, self.shape, self.dtype, grid, vals) @@ -81,46 +80,36 @@ def __init__(self, name, shape, dtype, grid, val=None): # Always init the grid, at least with 0.0 self[:] = 0.0 if val is None else val - def _convert_multislice(self, multislice): - start = [] - stop = [] - for i, idx in enumerate(multislice): - if isinstance(idx, slice): - start.append(idx.start or 0) - stop.append((idx.stop or self.shape[i]) - 1) - else: - start.append(idx) - stop.append(idx) - shape = [j - i + 1 for i, j in zip(start, stop)] - return start, stop, shape + @property + def ndim(self): + return len(self.shape) def __repr__(self): # TODO: need to return a VIEW return repr(self[:]) def __getitem__(self, index): + # TODO: ATM, no MPI support. if isinstance(index, int) and len(self.shape) > 1 and self.shape != val.shape: _force_exit("Retrieval from a YaskGrid, non-matching shapes.") if isinstance(index, slice): debug("YaskGrid: Getting whole array or block via single slice") if index.step is not None: _force_exit("Retrieval from a YaskGrid, unsupported stepping != 1.") - start = [index.start or 0 for j in self.shape] - stop = [(index.stop or j) - 1 for j in self.shape] - shape = [j - i + 1 for i, j in zip(start, stop)] + start, stop, shape = convert_multislice(as_tuple(index), self.shape) # TODO: using empty ndarray for now, will need to return views... - output = np.empty(shape, self.dtype, 'C'); - self.grid.get_elements_in_slice(output.data, start, stop) - return output + out = np.empty(shape*self.ndim, self.dtype, 'C'); + self.grid.get_elements_in_slice(out.data, start*self.ndim, stop*self.ndim) + return out elif all(isinstance(i, int) for i in as_tuple(index)): debug("YaskGrid: Getting single entry") - return self.grid.get_element(*as_tuple(index)) + return self.grid.get_element(*as_tuple(normalize_index(index, self.shape))) else: debug("YaskGrid: Getting whole array or block via multiple slices/indices") - start, stop, shape = self._convert_multislice(index) - output = np.empty(shape, self.dtype, 'C'); - self.grid.get_elements_in_slice(output.data, start, stop) - return output + start, stop, shape = convert_multislice(index, self.shape) + out = np.empty(shape, self.dtype, 'C'); + self.grid.get_elements_in_slice(out.data, start, stop) + return out def __setitem__(self, index, val): # TODO: ATM, no MPI support. @@ -129,22 +118,22 @@ def __setitem__(self, index, val): if isinstance(index, slice): debug("YaskGrid: Setting whole array or block via single slice") if isinstance(val, np.ndarray): - start = [index.start or 0 for j in self.shape] - stop = [(index.stop or j) - 1 for j in self.shape] - self.grid.set_elements_in_slice(val, start, stop) + start, stop, _ = convert_multislice(as_tuple(index), self.shape) + self.grid.set_elements_in_slice(val, start*self.ndim, stop*self.ndim) else: self.grid.set_all_elements_same(val) elif all(isinstance(i, int) for i in as_tuple(index)): debug("YaskGrid: Setting single entry") - self.grid.set_element(val, *as_tuple(index)) + self.grid.set_element(val, *as_tuple(normalize_index(index, self.shape))) else: debug("YaskGrid: Setting whole array or block via multiple slices/indices") - start, stop, _ = self._convert_multislice(index) + start, stop, _ = convert_multislice(index, self.shape) if isinstance(val, np.ndarray): self.grid.set_elements_in_slice(val, start, stop) else: - # TODO: NEED set_elements_in_slice acceptign single value - assert False, "STIIL WIP" + self.grid.set_elements_in_slice_same(val, start, stop) + + # TODO: override __incr__ etc class YaskRewriter(BasicRewriter): @@ -355,3 +344,24 @@ def _force_exit(emsg): Handle fatal errors. """ raise DLEException("Couldn't startup YASK [%s]. Exiting..." % emsg) + + +# Generic utility functions + +def convert_multislice(multislice, shape): + start = [] + stop = [] + for i, idx in enumerate(multislice): + if isinstance(idx, slice): + start.append(idx.start or 0) + stop.append((idx.stop or shape[i]) - 1) + else: + start.append(normalize_index(idx or 0, shape)) + stop.append(normalize_index(idx or (shape[i] - 1), shape)) + shape = [j - i + 1 for i, j in zip(start, stop)] + return start, stop, shape + + +def normalize_index(index, shape): + normalized = [i if i >= 0 else j + i for i, j in zip(as_tuple(index), shape)] + return normalized[0] if len(normalized) == 1 else normalized diff --git a/tests/test_yask.py b/tests/test_yask.py index 86c51bc4e8..b81ffc0543 100644 --- a/tests/test_yask.py +++ b/tests/test_yask.py @@ -30,19 +30,33 @@ def test_data_movement_1D(): assert all(i == 0 for i in u.data[2:]) -@pytest.mark.xfail(reason="FAIL on block insertion") def test_data_movement_nD(): u = DenseData(name='yu3D', shape=(10, 10, 10), dimensions=(x, y, z)) + + # Test simple insertion and extraction u.data[0, 1, 1] = 1. assert u.data[0, 0, 0] == 0. assert u.data[0, 1, 1] == 1. assert np.all(u.data[:] == u.data[:,:,:]) - # Test block insertion - block = np.ndarray(shape=(1, 10, 1)) - block.fill(5.) - u.data[5,:,5] = block - assert np.all(u.data[5,:,5] == block) + # Test negative indices + assert u.data[0, -9, -9] == 1. + u.data[6,0,0] = 1. + assert u.data[-4,:,:].sum() == 1. + + # Test setting whole array to given value + u.data[:] = 3. + assert np.all(u.data[:] == 3.) + + # Test insertion of single value into block + u.data[5,:,5] = 5. + assert np.all(u.data[5,:,5] == 5.) + + # Test insertion of block into block + block = np.ndarray(shape=(1, 10, 1), dtype=np.float32) + block.fill(4.) + u.data[4,:,4] = block + assert np.all(u.data[4,:,4] == block) def test_storage_layout(): From fbbb278831a337077873ab339fdc2ef3ae6197a7 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Tue, 20 Jun 2017 18:21:17 +0200 Subject: [PATCH 12/18] yask: Clean up and enhance YaskGrid [Unfinished] --- devito/dle/backends/yask.py | 127 ++++++++++++++++++++---------------- tests/test_yask.py | 12 +++- 2 files changed, 82 insertions(+), 57 deletions(-) diff --git a/devito/dle/backends/yask.py b/devito/dle/backends/yask.py index 57dcb805d8..cad227dd90 100644 --- a/devito/dle/backends/yask.py +++ b/devito/dle/backends/yask.py @@ -1,4 +1,6 @@ +import numbers import os +import sys import numpy as np from sympy import Indexed @@ -53,7 +55,7 @@ def grids(self): mapper[grid.get_name()] = grid return mapper - def setdefault(self, name, vals=None): + def setdefault(self, name, buffer=None): """ Add and return a new grid ``name``. If a grid ``name`` already exists, then return it without performing any other actions. @@ -66,74 +68,73 @@ def setdefault(self, name, vals=None): grid = self.hook_soln.new_grid(name, *self.dimensions) # Allocate memory self.hook_soln.prepare_solution() - return YaskGrid(name, self.shape, self.dtype, grid, vals) + return YaskGrid(name, grid, self.shape, self.dtype, buffer) class YaskGrid(object): - def __init__(self, name, shape, dtype, grid, val=None): + """ + An implementation of an array that behaves similarly to a ``numpy.ndarray``, + suitable for the YASK storage layout. + + Subclassing ``numpy.ndarray`` would have led to shadow data copies, because + of the different storage layout. + """ + + def __init__(self, name, grid, shape, dtype, buffer=None): self.name = name self.shape = shape self.dtype = dtype self.grid = grid # Always init the grid, at least with 0.0 - self[:] = 0.0 if val is None else val - - @property - def ndim(self): - return len(self.shape) - - def __repr__(self): - # TODO: need to return a VIEW - return repr(self[:]) + self[:] = 0.0 if buffer is None else val def __getitem__(self, index): # TODO: ATM, no MPI support. - if isinstance(index, int) and len(self.shape) > 1 and self.shape != val.shape: - _force_exit("Retrieval from a YaskGrid, non-matching shapes.") - if isinstance(index, slice): - debug("YaskGrid: Getting whole array or block via single slice") - if index.step is not None: - _force_exit("Retrieval from a YaskGrid, unsupported stepping != 1.") - start, stop, shape = convert_multislice(as_tuple(index), self.shape) - # TODO: using empty ndarray for now, will need to return views... - out = np.empty(shape*self.ndim, self.dtype, 'C'); - self.grid.get_elements_in_slice(out.data, start*self.ndim, stop*self.ndim) - return out - elif all(isinstance(i, int) for i in as_tuple(index)): + start, stop, shape = convert_multislice(as_tuple(index), self.shape) + if not shape: debug("YaskGrid: Getting single entry") - return self.grid.get_element(*as_tuple(normalize_index(index, self.shape))) + assert start == stop + out = self.grid.get_element(*start) else: - debug("YaskGrid: Getting whole array or block via multiple slices/indices") - start, stop, shape = convert_multislice(index, self.shape) - out = np.empty(shape, self.dtype, 'C'); + debug("YaskGrid: Getting full-array/block via slices/indices") + out = np.empty(shape, self.dtype, 'C') self.grid.get_elements_in_slice(out.data, start, stop) - return out + return out def __setitem__(self, index, val): # TODO: ATM, no MPI support. - if isinstance(index, int) and len(self.shape) > 1 and self.shape != val.shape: - _force_exit("Insertion into a YaskGrid, non-matching shapes.") - if isinstance(index, slice): - debug("YaskGrid: Setting whole array or block via single slice") - if isinstance(val, np.ndarray): - start, stop, _ = convert_multislice(as_tuple(index), self.shape) - self.grid.set_elements_in_slice(val, start*self.ndim, stop*self.ndim) - else: - self.grid.set_all_elements_same(val) - elif all(isinstance(i, int) for i in as_tuple(index)): + start, stop, shape = convert_multislice(as_tuple(index), self.shape, 'set') + if all(i == 1 for i in shape): debug("YaskGrid: Setting single entry") - self.grid.set_element(val, *as_tuple(normalize_index(index, self.shape))) + assert start == stop + self.grid.set_element(val, *start) else: - debug("YaskGrid: Setting whole array or block via multiple slices/indices") - start, stop, _ = convert_multislice(index, self.shape) + debug("YaskGrid: Setting full-array/block via multiple slices/indices") if isinstance(val, np.ndarray): self.grid.set_elements_in_slice(val, start, stop) else: self.grid.set_elements_in_slice_same(val, start, stop) - # TODO: override __incr__ etc + def __getslice__(self, start, stop): + if stop == sys.maxint: + stop = None + return self.__getitem__(slice(start, stop)) + + def __setslice__(self, start, stop, val): + if stop == sys.maxint: + stop = None + self.__setitem__(slice(start, stop), val) + + def __repr__(self): + return repr(self[:]) + + def __eq__(self, other): + if isinstance(other, (np.ndarray, numbers.Number)): + return self[:] == other + else: + raise NotImplementedError class YaskRewriter(BasicRewriter): @@ -343,23 +344,39 @@ def _force_exit(emsg): """ Handle fatal errors. """ - raise DLEException("Couldn't startup YASK [%s]. Exiting..." % emsg) + raise DLEException("YASK Error [%s]. Exiting..." % emsg) # Generic utility functions -def convert_multislice(multislice, shape): - start = [] - stop = [] - for i, idx in enumerate(multislice): - if isinstance(idx, slice): - start.append(idx.start or 0) - stop.append((idx.stop or shape[i]) - 1) +def convert_multislice(multislice, shape, mode='get'): + assert mode in ['get', 'set'] + multislice = as_tuple(multislice) + + # Convert dimensions + cstart = [] + cstop = [] + cshape = [] + for i, v in enumerate(multislice): + if isinstance(v, slice): + if v.step is not None: + _force_exit("Unsupported stepping != 1.") + cstart.append(v.start or 0) + cstop.append((v.stop or shape[i]) - 1) + cshape.append(cstop[-1] - cstart[-1] + 1) else: - start.append(normalize_index(idx or 0, shape)) - stop.append(normalize_index(idx or (shape[i] - 1), shape)) - shape = [j - i + 1 for i, j in zip(start, stop)] - return start, stop, shape + cstart.append(normalize_index(v if v is not None else 0, shape)) + cstop.append(normalize_index(v if v is not None else (shape[i]-1), shape)) + if mode == 'set': + cshape.append(1) + + # Remainder (e.g., requesting A[1] and A has shape (3,3)) + nremainder = len(shape) - len(multislice) + cstart.extend([0]*nremainder) + cstop.extend([shape[i + j] - 1 for j in range(nremainder)]) + cshape.extend([shape[i + j] for j in range(nremainder)]) + + return cstart, cstop, cshape def normalize_index(index, shape): diff --git a/tests/test_yask.py b/tests/test_yask.py index b81ffc0543..d7053523a1 100644 --- a/tests/test_yask.py +++ b/tests/test_yask.py @@ -37,7 +37,10 @@ def test_data_movement_nD(): u.data[0, 1, 1] = 1. assert u.data[0, 0, 0] == 0. assert u.data[0, 1, 1] == 1. - assert np.all(u.data[:] == u.data[:,:,:]) + print u.data + assert np.all(u.data == u.data[:,:,:]) + assert 1. in u.data[0] + assert 1. in u.data[0, 1] # Test negative indices assert u.data[0, -9, -9] == 1. @@ -46,7 +49,7 @@ def test_data_movement_nD(): # Test setting whole array to given value u.data[:] = 3. - assert np.all(u.data[:] == 3.) + assert np.all(u.data == 3.) # Test insertion of single value into block u.data[5,:,5] = 5. @@ -59,6 +62,11 @@ def test_data_movement_nD(): assert np.all(u.data[4,:,4] == block) +def test_data_arithmetic_nD(): + u = DenseData(name='yu3D', shape=(10, 10, 10), dimensions=(x, y, z)) + assert np.all(u.data == 0) + + def test_storage_layout(): u = DenseData(name='yu', shape=(10, 10), dimensions=(x, y)) op = Operator(Eq(u, 1.), dse='noop', dle='noop') From cea4e139b5d5891c228b129073aeba653be6d7e8 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Wed, 21 Jun 2017 10:26:13 +0200 Subject: [PATCH 13/18] yask: Support binary ops in YaskGrid [Unfinished] --- devito/dle/backends/yask.py | 32 +++++++++++++++++++++++++++----- tests/test_yask.py | 23 ++++++++++++++++++++++- 2 files changed, 49 insertions(+), 6 deletions(-) diff --git a/devito/dle/backends/yask.py b/devito/dle/backends/yask.py index cad227dd90..fc1291084c 100644 --- a/devito/dle/backends/yask.py +++ b/devito/dle/backends/yask.py @@ -81,6 +81,9 @@ class YaskGrid(object): of the different storage layout. """ + # Force __rOP__ methods (OP={add,mul,...) to get arrays, not scalars, for efficiency + __array_priority__ = 1000 + def __init__(self, name, grid, shape, dtype, buffer=None): self.name = name self.shape = shape @@ -130,11 +133,30 @@ def __setslice__(self, start, stop, val): def __repr__(self): return repr(self[:]) - def __eq__(self, other): - if isinstance(other, (np.ndarray, numbers.Number)): - return self[:] == other - else: - raise NotImplementedError + def __meta_binop(op): + # Used to build all binary operations such as __eq__, __add__, etc. + # These all boil down to calling the numpy equivalents + def f(self, other): + return getattr(self[:], op)(other) + return f + __eq__ = __meta_binop('__eq__') + __ne__ = __meta_binop('__ne__') + __le__ = __meta_binop('__le__') + __lt__ = __meta_binop('__lt__') + __ge__ = __meta_binop('__ge__') + __gt__ = __meta_binop('__gt__') + __add__ = __meta_binop('__add__') + __radd__ = __meta_binop('__add__') + __sub__ = __meta_binop('__sub__') + __rsub__ = __meta_binop('__sub__') + __mul__ = __meta_binop('__mul__') + __rmul__ = __meta_binop('__mul__') + __div__ = __meta_binop('__div__') + __rdiv__ = __meta_binop('__div__') + __truediv__ = __meta_binop('__truediv__') + __rtruediv__ = __meta_binop('__truediv__') + __mod__ = __meta_binop('__mod__') + __rmod__ = __meta_binop('__mod__') class YaskRewriter(BasicRewriter): diff --git a/tests/test_yask.py b/tests/test_yask.py index d7053523a1..0d3b8d6370 100644 --- a/tests/test_yask.py +++ b/tests/test_yask.py @@ -64,7 +64,28 @@ def test_data_movement_nD(): def test_data_arithmetic_nD(): u = DenseData(name='yu3D', shape=(10, 10, 10), dimensions=(x, y, z)) - assert np.all(u.data == 0) + + # Simple arithmetic + u.data[:] = 1 + assert np.all(u.data == 1) + assert np.all(u.data + 2. == 3.) + assert np.all(u.data - 2. == -1.) + assert np.all(u.data * 2. == 2.) + assert np.all(u.data / 2. == 0.5) + assert np.all(u.data % 2 == 1.) + + # Increments and parital increments + u.data[:] += 2. + assert np.all(u.data == 3.) + u.data[9,:,:] += 1. + assert all(np.all(u.data[i,:,:] == 3.) for i in range(9)) + assert np.all(u.data[9,:,:] == 4.) + + # Right operations __rOP__ + u.data[:] = 1. + arr = np.ndarray(shape=(10, 10, 10), dtype=np.float32) + arr.fill(2.) + assert np.all(arr - u.data == -1.) def test_storage_layout(): From 40f23d8f84a4d86b3db5d95fb11ee65c1dacdcc9 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Wed, 21 Jun 2017 10:45:02 +0200 Subject: [PATCH 14/18] yask: Specialize __setitem__ in YaskGrid [Unfinished] --- devito/dle/backends/yask.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/devito/dle/backends/yask.py b/devito/dle/backends/yask.py index fc1291084c..2de3d88367 100644 --- a/devito/dle/backends/yask.py +++ b/devito/dle/backends/yask.py @@ -101,7 +101,7 @@ def __getitem__(self, index): assert start == stop out = self.grid.get_element(*start) else: - debug("YaskGrid: Getting full-array/block via slices/indices") + debug("YaskGrid: Getting full-array/block via index [%s]" % str(index)) out = np.empty(shape, self.dtype, 'C') self.grid.get_elements_in_slice(out.data, start, stop) return out @@ -113,12 +113,15 @@ def __setitem__(self, index, val): debug("YaskGrid: Setting single entry") assert start == stop self.grid.set_element(val, *start) + elif isinstance(val, np.ndarray): + debug("YaskGrid: Setting full-array/block via index [%s]" % str(index)) + self.grid.set_elements_in_slice(val, start, stop) + elif all(i == j-1 for i, j in zip(shape, self.shape)): + debug("YaskGrid: Setting full-array to given scalar via single grid sweep") + self.grid.set_all_elements_same(val) else: - debug("YaskGrid: Setting full-array/block via multiple slices/indices") - if isinstance(val, np.ndarray): - self.grid.set_elements_in_slice(val, start, stop) - else: - self.grid.set_elements_in_slice_same(val, start, stop) + debug("YaskGrid: Setting block to given scalar via index [%s]" % str(index)) + self.grid.set_elements_in_slice_same(val, start, stop) def __getslice__(self, start, stop): if stop == sys.maxint: From 3cf50bd8155431c369c880044c4e8b2177ff9bcd Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Thu, 22 Jun 2017 09:30:46 +0200 Subject: [PATCH 15/18] yask: Refactor YaskGrid instantiation --- devito/dle/__init__.py | 2 +- devito/dle/backends/yask.py | 39 ++++++++++++++++++++++++------------- devito/interfaces.py | 26 +++++++++++++++++-------- 3 files changed, 45 insertions(+), 22 deletions(-) diff --git a/devito/dle/__init__.py b/devito/dle/__init__.py index f7f213c87f..e1e708e079 100644 --- a/devito/dle/__init__.py +++ b/devito/dle/__init__.py @@ -1,4 +1,4 @@ from devito.dle.inspection import * # noqa from devito.dle.manipulation import * # noqa from devito.dle.transformer import * # noqa -from devito.dle.backends import init, make_grid # noqa +from devito.dle.backends import YaskGrid, init # noqa diff --git a/devito/dle/backends/yask.py b/devito/dle/backends/yask.py index 2de3d88367..1a690cda4a 100644 --- a/devito/dle/backends/yask.py +++ b/devito/dle/backends/yask.py @@ -14,7 +14,7 @@ from devito.visitors import FindSymbols from devito.tools import as_tuple -__all__ = ['YaskRewriter', 'init', 'make_grid'] +__all__ = ['YaskRewriter', 'init', 'YaskGrid'] YASK = None @@ -55,7 +55,7 @@ def grids(self): mapper[grid.get_name()] = grid return mapper - def setdefault(self, name, buffer=None): + def setdefault(self, name): """ Add and return a new grid ``name``. If a grid ``name`` already exists, then return it without performing any other actions. @@ -68,7 +68,7 @@ def setdefault(self, name, buffer=None): grid = self.hook_soln.new_grid(name, *self.dimensions) # Allocate memory self.hook_soln.prepare_solution() - return YaskGrid(name, grid, self.shape, self.dtype, buffer) + return grid class YaskGrid(object): @@ -84,11 +84,27 @@ class YaskGrid(object): # Force __rOP__ methods (OP={add,mul,...) to get arrays, not scalars, for efficiency __array_priority__ = 1000 - def __init__(self, name, grid, shape, dtype, buffer=None): + def __new__(cls, name, shape, dimensions, dtype, buffer=None): + """ + Create a new YASK Grid and attach it to a "fake" solution. + """ + # Init YASK if not initialized already + init(dimensions, shape, dtype) + # Only create a YaskGrid if the requested grid is a dense one + if tuple(i.name for i in dimensions) == YASK.dimensions: + obj = super(YaskGrid, cls).__new__(cls) + obj.__init__(name, shape, dimensions, dtype, buffer) + return obj + else: + return None + + def __init__(self, name, shape, dimensions, dtype, buffer=None): self.name = name self.shape = shape + self.dimensions = dimensions self.dtype = dtype - self.grid = grid + + self.grid = YASK.setdefault(name) # Always init the grid, at least with 0.0 self[:] = 0.0 if buffer is None else val @@ -161,6 +177,11 @@ def f(self, other): __mod__ = __meta_binop('__mod__') __rmod__ = __meta_binop('__mod__') + @property + def ndpointer(self): + # TODO: see corresponding comment in interfaces.py about CMemory + return self + class YaskRewriter(BasicRewriter): @@ -357,14 +378,6 @@ def init(dimensions, shape, dtype, architecture='hsw', isa='avx2'): dle("YASK backend successfully initialized!") -def make_grid(name, shape, dimensions, dtype): - """ - Create a new YASK Grid and attach it to a "fake" solution. - """ - init(dimensions, shape, dtype) - return YASK.setdefault(name) - - def _force_exit(emsg): """ Handle fatal errors. diff --git a/devito/interfaces.py b/devito/interfaces.py index f2063e32c1..892c69be99 100644 --- a/devito/interfaces.py +++ b/devito/interfaces.py @@ -305,15 +305,25 @@ def _allocate_memory(self): from devito.parameters import configuration if configuration['dle'] == 'yask': - from devito.dle import make_grid - self._data = make_grid(self.name, self.shape, self.indices, self.dtype) + # TODO: Use inheritance + # TODO: Refactor CMemory to be our _data_object, while _data will + # be the YaskGrid itself. + from devito.dle import YaskGrid + debug("Allocating YaskGrid for %s (%s)" % (self.name, str(self.shape))) + + self._data_object = YaskGrid(self.name, self.shape, self.indices, self.dtype) + if self._data_object is not None: + return + + debug("Failed. Reverting to plain allocation...") + + debug("Allocating memory for %s (%s)" % (self.name, str(self.shape))) + + self._data_object = CMemory(self.shape, dtype=self.dtype) + if self.numa: + first_touch(self) else: - debug("Allocating memory for %s (%s)" % (self.name, str(self.shape))) - self._data_object = CMemory(self.shape, dtype=self.dtype) - if self.numa: - first_touch(self) - else: - self.data.fill(0) + self.data.fill(0) @property def data(self): From 80a662014ac5a3ec445f45b1b8b82e76d9f0d915 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Thu, 22 Jun 2017 09:31:00 +0200 Subject: [PATCH 16/18] yask: Run tests only if YASK is available --- tests/test_yask.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/test_yask.py b/tests/test_yask.py index 0d3b8d6370..46c920ecaa 100644 --- a/tests/test_yask.py +++ b/tests/test_yask.py @@ -4,8 +4,10 @@ import pytest # noqa from devito import Operator, DenseData, x, y, z +from devito.dle import YaskGrid from devito.parameters import configuration, defaults -from devito.dle.backends.yask import YaskGrid + +pexpect = pytest.importorskip('yask_compiler') # Run only if YASK is available def setup_module(module): @@ -37,7 +39,6 @@ def test_data_movement_nD(): u.data[0, 1, 1] = 1. assert u.data[0, 0, 0] == 0. assert u.data[0, 1, 1] == 1. - print u.data assert np.all(u.data == u.data[:,:,:]) assert 1. in u.data[0] assert 1. in u.data[0, 1] From e49f0071dfc3fd5d2325182a48108cacf08a5247 Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Thu, 22 Jun 2017 09:35:03 +0200 Subject: [PATCH 17/18] flake8 fixes --- devito/dle/backends/yask.py | 5 ++--- tests/test_yask.py | 20 ++++++++++---------- 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/devito/dle/backends/yask.py b/devito/dle/backends/yask.py index 1a690cda4a..b11025395c 100644 --- a/devito/dle/backends/yask.py +++ b/devito/dle/backends/yask.py @@ -1,4 +1,3 @@ -import numbers import os import sys @@ -10,7 +9,7 @@ from devito.dle import retrieve_iteration_tree from devito.dle.backends import BasicRewriter, dle_pass from devito.exceptions import CompilationError, DLEException -from devito.logger import debug, dle, dle_warning, error +from devito.logger import debug, dle, dle_warning from devito.visitors import FindSymbols from devito.tools import as_tuple @@ -107,7 +106,7 @@ def __init__(self, name, shape, dimensions, dtype, buffer=None): self.grid = YASK.setdefault(name) # Always init the grid, at least with 0.0 - self[:] = 0.0 if buffer is None else val + self[:] = 0.0 if buffer is None else buffer def __getitem__(self, index): # TODO: ATM, no MPI support. diff --git a/tests/test_yask.py b/tests/test_yask.py index 46c920ecaa..0f7f98715f 100644 --- a/tests/test_yask.py +++ b/tests/test_yask.py @@ -39,28 +39,28 @@ def test_data_movement_nD(): u.data[0, 1, 1] = 1. assert u.data[0, 0, 0] == 0. assert u.data[0, 1, 1] == 1. - assert np.all(u.data == u.data[:,:,:]) + assert np.all(u.data == u.data[:, :, :]) assert 1. in u.data[0] assert 1. in u.data[0, 1] # Test negative indices assert u.data[0, -9, -9] == 1. - u.data[6,0,0] = 1. - assert u.data[-4,:,:].sum() == 1. + u.data[6, 0, 0] = 1. + assert u.data[-4, :, :].sum() == 1. # Test setting whole array to given value u.data[:] = 3. assert np.all(u.data == 3.) # Test insertion of single value into block - u.data[5,:,5] = 5. - assert np.all(u.data[5,:,5] == 5.) + u.data[5, :, 5] = 5. + assert np.all(u.data[5, :, 5] == 5.) # Test insertion of block into block block = np.ndarray(shape=(1, 10, 1), dtype=np.float32) block.fill(4.) - u.data[4,:,4] = block - assert np.all(u.data[4,:,4] == block) + u.data[4, :, 4] = block + assert np.all(u.data[4, :, 4] == block) def test_data_arithmetic_nD(): @@ -78,9 +78,9 @@ def test_data_arithmetic_nD(): # Increments and parital increments u.data[:] += 2. assert np.all(u.data == 3.) - u.data[9,:,:] += 1. - assert all(np.all(u.data[i,:,:] == 3.) for i in range(9)) - assert np.all(u.data[9,:,:] == 4.) + u.data[9, :, :] += 1. + assert all(np.all(u.data[i, :, :] == 3.) for i in range(9)) + assert np.all(u.data[9, :, :] == 4.) # Right operations __rOP__ u.data[:] = 1. From 888d6e641bb6f48d9374e6a2b9b69edb2e75342f Mon Sep 17 00:00:00 2001 From: Fabio Luporini Date: Tue, 27 Jun 2017 12:31:20 +0200 Subject: [PATCH 18/18] yask: Improve code documentation --- devito/dle/backends/yask.py | 42 ++++++++++++++++++++++++++++++++++++- 1 file changed, 41 insertions(+), 1 deletion(-) diff --git a/devito/dle/backends/yask.py b/devito/dle/backends/yask.py index b11025395c..f56d918a88 100644 --- a/devito/dle/backends/yask.py +++ b/devito/dle/backends/yask.py @@ -140,11 +140,13 @@ def __setitem__(self, index, val): def __getslice__(self, start, stop): if stop == sys.maxint: + # Emulate default NumPy behaviour stop = None return self.__getitem__(slice(start, stop)) def __setslice__(self, start, stop, val): if stop == sys.maxint: + # Emulate default NumPy behaviour stop = None self.__setitem__(slice(start, stop), val) @@ -332,7 +334,11 @@ def init(dimensions, shape, dtype, architecture='hsw', isa='avx2'): # Create a new stencil solution soln = cfac.new_solution("Hook") soln.set_step_dim("t") - soln.set_domain_dims(*[str(i) for i in dimensions]) # TODO: YASK only accepts x,y,z + + dimensions = [str(i) for i in dimensions] + if any(i not in ['x', 'y', 'z'] for i in dimensions): + _force_exit("Need a DenseData[x,y,z] for initialization") + soln.set_domain_dims(*dimensions) # TODO: YASK only accepts x,y,z # Number of bytes in each FP value soln.set_element_bytes(dtype().itemsize) @@ -387,6 +393,40 @@ def _force_exit(emsg): # Generic utility functions def convert_multislice(multislice, shape, mode='get'): + """ + Convert a multislice into a format suitable to YASK's get_elements_{...} + and set_elements_{...} grid routines. + + A multislice is the typical object received by NumPy ndarray's __getitem__ + and __setitem__ methods; this function, therefore, converts NumPy indices + into YASK indices. + + In particular, a multislice is either a single element or an iterable of + elements. An element can be a slice object, an integer index, or a tuple + of integer indices. + + In the general case in which ``multislice`` is an iterable, each element in + the iterable corresponds to a dimension in ``shape``. In this case, an element + can be either a slice or an integer index, but not a tuple of integers. + + If ``multislice`` is a single element, then it is interpreted as follows: :: + + * slice object: the slice spans the whole shape; + * single (integer) index: shape is one-dimensional, and the index + represents a specific entry; + * a tuple of (integer) indices: it must be ``len(multislice) == len(shape)``, + and each entry in ``multislice`` corresponds to a specific entry in a + dimension in ``shape``. + + The returned value is a 3-tuple ``(starts, ends, shapes)``, where ``starts, + ends, shapes`` are lists of length ``len(shape)``. By taking ``starts[i]`` and + `` ends[i]``, one gets the start and end points of the section of elements to + be accessed along dimension ``i``; ``shapes[i]`` gives the size of the section. + """ + + # Note: the '-1' below are because YASK uses '<=', rather than '<', to check + # bounds when iterating over grid dimensions + assert mode in ['get', 'set'] multislice = as_tuple(multislice)