doc/sphinx/api.rst

@@ -5,6 +5,7 @@ API Reference
    :maxdepth: 3
 
    kernel_compile_and_call.rst
+   enums.rst
    simplifications.rst
    datahandling.rst
    configuration.rst

doc/sphinx/enums.rst

+************
+Enumerations
+************
+
+.. automodule:: pystencils.enums
+   :members:

doc/sphinx/kernel_compile_and_call.rst

@@ -8,9 +8,14 @@ Creating kernels
 
 .. autofunction:: pystencils.create_kernel
 
-.. autofunction:: pystencils.create_indexed_kernel
+.. autoclass:: pystencils.CreateKernelConfig
+    :members:
 
-.. autofunction:: pystencils.create_staggered_kernel
+.. autofunction:: pystencils.kernelcreation.create_domain_kernel
+
+.. autofunction:: pystencils.kernelcreation.create_indexed_kernel
+
+.. autofunction:: pystencils.kernelcreation.create_staggered_kernel
 
 
 Code printing
@@ -22,11 +27,11 @@ Code printing
 GPU Indexing
 -------------
 
-.. autoclass:: pystencils.gpucuda.AbstractIndexing
+.. autoclass:: pystencils.gpu.AbstractIndexing
    :members:
 
-.. autoclass:: pystencils.gpucuda.BlockIndexing
+.. autoclass:: pystencils.gpu.BlockIndexing
    :members:
 
-.. autoclass:: pystencils.gpucuda.LineIndexing
+.. autoclass:: pystencils.gpu.LineIndexing
    :members:

doc/sphinx/simplifications.rst

@@ -10,13 +10,27 @@ AssignmentCollection
    :members:
 
 
+SimplificationStrategy
+======================
+
+.. autoclass:: pystencils.simp.SimplificationStrategy
+    :members:
+
 Simplifications
 ===============
 
-.. automodule:: pystencils.simp
-   :members:
+.. automodule:: pystencils.simp.simplifications
+    :members:
+
+Subexpression insertion
+=======================
 
+The subexpression insertions have the goal to insert subexpressions which will not reduce the number of FLOPs.
+For example a constant value kept as subexpression will lead to a new variable in the code which will occupy
+a register slot. On the other side a single variable could just be inserted in all assignments.
 
+.. automodule:: pystencils.simp.subexpression_insertion
+    :members:
 
 
 

pyproject.toml

+[project]
+name = "pystencils"
+description = "Speeding up stencil computations on CPUs and GPUs"
+dynamic = ["version"]
+readme = "README.md"
+authors = [
+    { name = "Martin Bauer" },
+    { name = "Jan Hönig " },
+    { name = "Markus Holzer" },
+    { name = "Frederik Hennig" },
+    { email = "cs10-codegen@fau.de" },
+]
+license = { file = "COPYING.txt" }
+requires-python = ">=3.10"
+dependencies = ["sympy>=1.9,<=1.12.1", "numpy>=1.8.0", "appdirs", "joblib", "pyyaml", "fasteners"]
+classifiers = [
+    "Development Status :: 4 - Beta",
+    "Framework :: Jupyter",
+    "Topic :: Software Development :: Code Generators",
+    "Topic :: Scientific/Engineering :: Physics",
+    "Intended Audience :: Developers",
+    "Intended Audience :: Science/Research",
+    "License :: OSI Approved :: GNU Affero General Public License v3 or later (AGPLv3+)",
+]
+
+[project.urls]
+"Bug Tracker" = "https://i10git.cs.fau.de/pycodegen/pystencils/-/issues"
+"Documentation" = "https://pycodegen.pages.i10git.cs.fau.de/pystencils/"
+"Source Code" = "https://i10git.cs.fau.de/pycodegen/pystencils"
+
+[project.optional-dependencies]
+gpu = ['cupy']
+alltrafos = ['islpy', 'py-cpuinfo']
+bench_db = ['blitzdb', 'pymongo', 'pandas']
+interactive = [
+    'matplotlib',
+    'ipy_table',
+    'imageio',
+    'jupyter',
+    'pyevtk',
+    'rich',
+    'graphviz',
+]
+use_cython = [
+    'Cython'
+]
+doc = [
+    'sphinx',
+    'sphinx_rtd_theme',
+    'nbsphinx',
+    'sphinxcontrib-bibtex',
+    'sphinx_autodoc_typehints',
+    'pandoc',
+]
+tests = [
+    'pytest',
+    'pytest-cov',
+    'pytest-html',
+    'ansi2html',
+    'pytest-xdist',
+    'flake8',
+    'nbformat',
+    'nbconvert',
+    'ipython',
+    'matplotlib',
+    'py-cpuinfo',
+    'randomgen>=1.18',
+]
+
+[build-system]
+requires = [
+    "setuptools>=61",
+    "versioneer[toml]>=0.29",
+    # 'Cython'
+]
+build-backend = "setuptools.build_meta"
+
+[tool.setuptools.package-data]
+pystencils = [
+    "include/*.h",
+    "boundaries/createindexlistcython.pyx"
+]
+
+[tool.setuptools.packages.find]
+where = ["src"]
+include = ["pystencils", "pystencils.*"]
+namespaces = false
+
+[tool.versioneer]
+# See the docstring in versioneer.py for instructions. Note that you must
+# re-run 'versioneer.py setup' after changing this section, and commit the
+# resulting files.
+VCS = "git"
+style = "pep440"
+versionfile_source = "src/pystencils/_version.py"
+versionfile_build = "pystencils/_version.py"
+tag_prefix = "release/"
+parentdir_prefix = "pystencils-"

pystencils/autodiff.py

-"""
-Provides tools for generation of auto-differentiable operations.
-
-See https://github.com/theHamsta/pystencils_autodiff
-
-Installation:
-
-.. code-block:: bash
-    pip install pystencils-autodiff
-"""
-
-raise NotImplementedError('pystencils-autodiff is not installed. Run `pip install pystencils-autodiff`')

pystencils/backends/cuda_known_functions.txt

-__prof_trigger
-printf
-
-__syncthreads
-__syncthreads_count
-__syncthreads_and
-__syncthreads_or
-__syncwarp
-__threadfence
-__threadfence_block
-__threadfence_system
-
-atomicAdd
-atomicSub
-atomicExch
-atomicMin
-atomicMax
-atomicInc
-atomicDec
-atomicAnd
-atomicOr
-atomicXor
-atomicCAS
-
-__all_sync
-__any_sync
-__ballot_sync
-__active_mask
-
-__shfl_sync
-__shfl_up_sync
-__shfl_down_sync
-__shfl_xor_sync
-
-__match_any_sync
-__match_all_sync
-
-__isGlobal
-__isShared
-__isConstant
-__isLocal
-
-tex1Dfetch
-tex1D
-tex2D
-tex3D
-
-sqrtf
-rsqrtf
-cbrtf
-rcbrtf
-hypotf
-rhypotf
-norm3df
-rnorm3df
-norm4df
-rnorm4df
-normf
-rnormf
-expf
-exp2f
-exp10f
-expm1f
-logf
-log2f
-log10f
-log1pf
-sinf
-cosf
-tanf
-sincosf
-sinpif
-cospif
-sincospif
-asinf
-acosf
-atanf
-atan2f
-sinhf
-coshf
-tanhf
-asinhf
-acoshf
-atanhf
-powf
-erff
-erfcf
-erfinvf
-erfcinvf
-erfcxf
-normcdff
-normcdfinvf
-lgammaf
-tgammaf
-fmaf
-frexpf
-ldexpf
-scalbnf
-scalblnf
-logbf
-ilogbf
-j0f
-j1f
-jnf
-y0f
-y1f
-ynf
-cyl_bessel_i0f
-cyl_bessel_i1f
-fmodf
-remainderf
-remquof
-modff
-fdimf
-truncf
-roundf
-rintf
-nearbyintf
-ceilf
-floorf
-lrintf
-lroundf
-llrintf
-llroundf
-
-sqrt
-rsqrt
-cbrt
-rcbrt
-hypot
-rhypot
-norm3d
-rnorm3d
-norm4d
-rnorm4d
-norm
-rnorm
-exp
-exp2
-exp10
-expm1
-log
-log2
-log10
-log1p
-sin
-cos
-tan
-sincos
-sinpi
-cospi
-sincospi
-asin
-acos
-atan
-atan2
-sinh
-cosh
-tanh
-asinh
-acosh
-atanh
-pow
-erf
-erfc
-erfinv
-erfcinv
-erfcx
-normcdf
-normcdfinv
-lgamma
-tgamma
-fma
-frexp
-ldexp
-scalbn
-scalbln
-logb
-ilogb
-j0
-j1
-jn
-y0
-y1
-yn
-cyl_bessel_i0
-cyl_bessel_i1
-fmod
-remainder
-remquo
-mod
-fdim
-trunc
-round
-rint
-nearbyint
-ceil
-floor
-lrint
-lround
-llrint
-llround
-
-__fdividef
-__sinf
-__cosf
-__tanf
-__sincosf
-__logf
-__log2f
-__log10f
-__expf
-__exp10f
-__powf
-
-__fadd_rn
-__fsub_rn
-__fmul_rn
-__fmaf_rn
-__frcp_rn
-__fsqrt_rn
-__frsqrt_rn
-__fdiv_rn
-
-__fadd_rz
-__fsub_rz
-__fmul_rz
-__fmaf_rz
-__frcp_rz
-__fsqrt_rz
-__frsqrt_rz
-__fdiv_rz
-
-__fadd_ru
-__fsub_ru
-__fmul_ru
-__fmaf_ru
-__frcp_ru
-__fsqrt_ru
-__frsqrt_ru
-__fdiv_ru
-
-__fadd_rd
-__fsub_rd
-__fmul_rd
-__fmaf_rd
-__frcp_rd
-__fsqrt_rd
-__frsqrt_rd
-__fdiv_rd
-
-__fdividef
-__expf
-__exp10f
-__logf
-__log2f
-__log10f
-__sinf
-__cosf
-__sincosf
-__tanf
-__powf
-
-__dadd_rn
-__dsub_rn
-__dmul_rn
-__fma_rn
-__ddiv_rn
-__drcp_rn
-__dsqrt_rn
-
-__dadd_rz
-__dsub_rz
-__dmul_rz
-__fma_rz
-__ddiv_rz
-__drcp_rz
-__dsqrt_rz
-
-__dadd_ru
-__dsub_ru
-__dmul_ru
-__fma_ru
-__ddiv_ru
-__drcp_ru
-__dsqrt_ru
-
-__dadd_rd
-__dsub_rd
-__dmul_rd
-__fma_rd
-__ddiv_rd
-__drcp_rd
-__dsqrt_rd

pystencils/backends/opencl1.1_known_functions.txt

-acos
-acosh
-acospi
-asin
-asinh
-asinpi
-atan
-atan2
-atanh
-atanpi
-atan2pi
-cbrt
-ceil
-copysign
-cos
-cosh
-cospi
-erfc
-erf
-exp
-exp2
-exp10
-expm1
-fabs
-fdim
-floor
-fma
-fmax
-fmax
-fmin45
-fmin
-fmod
-fract
-frexp
-hypot
-ilogb
-ldexp
-lgamma
-lgamma_r
-log
-log2
-log10
-log1p
-logb
-mad 
-maxmag 
-minmag 
-modf 
-nextafter
-pow 
-pown
-powr
-remquo 
-intn
-remquo
-rint
-rootn
-rootn
-round
-rsqrt
-sin
-sincos
-sinh
-sinpi
-sqrt
-tan
-tanh
-tanpi
-tgamma
-trunc
-
-
-half_cos 
-half_divide
-half_exp 
-half_exp2 
-half_exp10 
-half_log 
-half_log2 
-half_log10 
-half_powr 
-half_recip 
-half_rsqrt 
-half_sin 
-half_sqrt 
-half_tan 
-native_cos 
-native_divide
-native_exp 
-native_exp2 
-native_exp10 
-native_log 
-native_log2 
-native_log10 
-native_powr 
-native_recip 
-native_rsqrt 
-native_sin 
-native_sqrt 
-native_tan 

pystencils/backends/opencl_backend.py

-from os.path import dirname, join
-
-import pystencils.data_types
-from pystencils.astnodes import Node
-from pystencils.backends.cbackend import CustomSympyPrinter, generate_c
-from pystencils.backends.cuda_backend import CudaBackend, CudaSympyPrinter
-from pystencils.fast_approximation import fast_division, fast_inv_sqrt, fast_sqrt
-
-with open(join(dirname(__file__), 'opencl1.1_known_functions.txt')) as f:
-    lines = f.readlines()
-    OPENCL_KNOWN_FUNCTIONS = {l.strip(): l.strip() for l in lines if l}
-
-
-def generate_opencl(ast_node: Node, signature_only: bool = False, custom_backend=None, with_globals=True) -> str:
-    """Prints an abstract syntax tree node (made for target 'gpu') as OpenCL code.
-
-    Args:
-        ast_node: ast representation of kernel
-        signature_only: generate signature without function body
-        custom_backend: use own custom printer for code generation
-        with_globals: enable usage of global variables
-
-    Returns:
-        OpenCL code for the ast node and its descendants
-    """
-    return generate_c(ast_node, signature_only, dialect='opencl',
-                      custom_backend=custom_backend, with_globals=with_globals)
-
-
-class OpenClBackend(CudaBackend):
-
-    def __init__(self,
-                 sympy_printer=None,
-                 signature_only=False):
-        if not sympy_printer:
-            sympy_printer = OpenClSympyPrinter()
-
-        super().__init__(sympy_printer, signature_only)
-        self._dialect = 'opencl'
-
-    def _print_Type(self, node):
-        code = super()._print_Type(node)
-        if isinstance(node, pystencils.data_types.PointerType):
-            return "__global " + code
-        else:
-            return code
-
-    def _print_ThreadBlockSynchronization(self, node):
-        raise NotImplementedError()
-
-    def _print_TextureDeclaration(self, node):
-        raise NotImplementedError()
-
-
-class OpenClSympyPrinter(CudaSympyPrinter):
-    language = "OpenCL"
-
-    DIMENSION_MAPPING = {
-        'x': '0',
-        'y': '1',
-        'z': '2'
-    }
-    INDEXING_FUNCTION_MAPPING = {
-        'blockIdx': 'get_group_id',
-        'threadIdx': 'get_local_id',
-        'blockDim': 'get_local_size',
-        'gridDim': 'get_global_size'
-    }
-
-    def __init__(self):
-        CustomSympyPrinter.__init__(self)
-        self.known_functions = OPENCL_KNOWN_FUNCTIONS
-
-    def _print_Type(self, node):
-        code = super()._print_Type(node)
-        if isinstance(node, pystencils.data_types.PointerType):
-            return "__global " + code
-        else:
-            return code
-
-    def _print_ThreadIndexingSymbol(self, node):
-        symbol_name: str = node.name
-        function_name, dimension = tuple(symbol_name.split("."))
-        dimension = self.DIMENSION_MAPPING[dimension]
-        function_name = self.INDEXING_FUNCTION_MAPPING[function_name]
-        return f"(int64_t) {function_name}({dimension})"
-
-    def _print_TextureAccess(self, node):
-        raise NotImplementedError()
-
-    # For math functions, OpenCL is more similar to the C++ printer CustomSympyPrinter
-    # since built-in math functions are generic.
-    # In CUDA, you have to differentiate between `sin` and `sinf`
-    try:
-        _print_math_func = CustomSympyPrinter._print_math_func
-    except AttributeError:
-        pass
-    _print_Pow = CustomSympyPrinter._print_Pow
-
-    def _print_Function(self, expr):
-        if isinstance(expr, fast_division):
-            return "native_divide(%s, %s)" % tuple(self._print(a) for a in expr.args)
-        elif isinstance(expr, fast_sqrt):
-            return f"native_sqrt({tuple(self._print(a) for a in expr.args)})"
-        elif isinstance(expr, fast_inv_sqrt):
-            return f"native_rsqrt({tuple(self._print(a) for a in expr.args)})"
-        return CustomSympyPrinter._print_Function(self, expr)

pystencils/cache.py

-import os
-from collections.abc import Hashable
-from functools import partial
-from itertools import chain
-
-try:
-    from functools import lru_cache as memorycache
-except ImportError:
-    from backports.functools_lru_cache import lru_cache as memorycache
-
-
-try:
-    from joblib import Memory
-    from appdirs import user_cache_dir
-    if 'PYSTENCILS_CACHE_DIR' in os.environ:
-        cache_dir = os.environ['PYSTENCILS_CACHE_DIR']
-    else:
-        cache_dir = user_cache_dir('pystencils')
-    disk_cache = Memory(cache_dir, verbose=False).cache
-    disk_cache_no_fallback = disk_cache
-except ImportError:
-    # fallback to in-memory caching if joblib is not available
-    disk_cache = memorycache(maxsize=64)
-
-    def disk_cache_no_fallback(o):
-        return o
-
-
-def _wrapper(wrapped_func, cached_func, *args, **kwargs):
-    if all(isinstance(a, Hashable) for a in chain(args, kwargs.values())):
-        return cached_func(*args, **kwargs)
-    else:
-        return wrapped_func(*args, **kwargs)
-
-
-def memorycache_if_hashable(maxsize=128, typed=False):
-
-    def wrapper(func):
-        return partial(_wrapper, func, memorycache(maxsize, typed)(func))
-
-    return wrapper
-
-# Disable memory cache:
-# disk_cache = lambda o: o
-# disk_cache_no_fallback = lambda o: o

pystencils/datahandling/pycuda.py

-try:
-    import pycuda.gpuarray as gpuarray
-except ImportError:
-    gpuarray = None
-import numpy as np
-
-import pystencils
-
-
-class PyCudaArrayHandler:
-
-    def __init__(self):
-        import pycuda.autoinit  # NOQA
-
-    def zeros(self, shape, dtype=np.float64, order='C'):
-        cpu_array = np.zeros(shape=shape, dtype=dtype, order=order)
-        return self.to_gpu(cpu_array)
-
-    def ones(self, shape, dtype=np.float64, order='C'):
-        cpu_array = np.ones(shape=shape, dtype=dtype, order=order)
-        return self.to_gpu(cpu_array)
-
-    def empty(self, shape, dtype=np.float64, layout=None):
-        if layout:
-            cpu_array = pystencils.field.create_numpy_array_with_layout(shape=shape, dtype=dtype, layout=layout)
-            return self.to_gpu(cpu_array)
-        else:
-            return gpuarray.empty(shape, dtype)
-
-    @staticmethod
-    def to_gpu(array):
-        return gpuarray.to_gpu(array)
-
-    @staticmethod
-    def upload(array, numpy_array):
-        array.set(numpy_array)
-
-    @staticmethod
-    def download(array, numpy_array):
-        array.get(numpy_array)
-
-    def randn(self, shape, dtype=np.float64):
-        cpu_array = np.random.randn(*shape).astype(dtype)
-        return self.to_gpu(cpu_array)
-
-    from_numpy = to_gpu
-
-
-class PyCudaNotAvailableHandler:
-    def __getattribute__(self, name):
-        raise NotImplementedError("Unable to initiaize PyCuda! "
-                                  "Try to run `import pycuda.autoinit` to check whether PyCuda is working correctly!")

pystencils/datahandling/pyopencl.py

-try:
-    import pyopencl.array as gpuarray
-except ImportError:
-    gpuarray = None
-
-import numpy as np
-
-import pystencils
-
-
-class PyOpenClArrayHandler:
-
-    def __init__(self, queue):
-        if not queue:
-            from pystencils.opencl.opencljit import get_global_cl_queue
-            queue = get_global_cl_queue()
-        assert queue, "OpenCL queue missing!\n" \
-            "Use `import pystencils.opencl.autoinit` if you want it to be automatically created"
-        self.queue = queue
-
-    def zeros(self, shape, dtype=np.float64, order='C'):
-        cpu_array = np.zeros(shape=shape, dtype=dtype, order=order)
-        return self.to_gpu(cpu_array)
-
-    def ones(self, shape, dtype=np.float64, order='C'):
-        cpu_array = np.ones(shape=shape, dtype=dtype, order=order)
-        return self.to_gpu(cpu_array)
-
-    def empty(self, shape, dtype=np.float64, layout=None):
-        if layout:
-            cpu_array = pystencils.field.create_numpy_array_with_layout(shape=shape, dtype=dtype, layout=layout)
-            return self.to_gpu(cpu_array)
-        else:
-            return gpuarray.empty(self.queue, shape, dtype)
-
-    def to_gpu(self, array):
-        return gpuarray.to_device(self.queue, array)
-
-    def upload(self, gpuarray, numpy_array):
-        gpuarray.set(numpy_array, self.queue)
-
-    def download(self, gpuarray, numpy_array):
-        gpuarray.get(self.queue, numpy_array)
-
-    def randn(self, shape, dtype=np.float64):
-        cpu_array = np.random.randn(*shape).astype(dtype)
-        return self.from_numpy(cpu_array)
-
-    from_numpy = to_gpu

pystencils/gpucuda/texture_utils.py

-# -*- coding: utf-8 -*-
-#
-# Copyright © 2019 Stephan Seitz <stephan.seitz@fau.de>
-#
-# Distributed under terms of the GPLv3 license.
-
-"""
-
-"""
-from typing import Union
-
-import numpy as np
-
-try:
-    import pycuda.driver as cuda
-    from pycuda import gpuarray
-    import pycuda
-except Exception:
-    pass
-
-
-def ndarray_to_tex(tex_ref,  # type: Union[cuda.TextureReference, cuda.SurfaceReference]
-                   ndarray,
-                   address_mode=None,
-                   filter_mode=None,
-                   use_normalized_coordinates=False,
-                   read_as_integer=False):
-
-    if isinstance(address_mode, str):
-        address_mode = getattr(pycuda.driver.address_mode, address_mode.upper())
-    if address_mode is None:
-        address_mode = cuda.address_mode.BORDER
-    if filter_mode is None:
-        filter_mode = cuda.filter_mode.LINEAR
-
-    if isinstance(ndarray, np.ndarray):
-        cu_array = cuda.np_to_array(ndarray, 'C')
-    elif isinstance(ndarray, gpuarray.GPUArray):
-        cu_array = cuda.gpuarray_to_array(ndarray, 'C')
-    else:
-        raise TypeError(
-            'ndarray must be numpy.ndarray or pycuda.gpuarray.GPUArray')
-
-    tex_ref.set_array(cu_array)
-
-    tex_ref.set_address_mode(0, address_mode)
-    if ndarray.ndim >= 2:
-        tex_ref.set_address_mode(1, address_mode)
-    if ndarray.ndim >= 3:
-        tex_ref.set_address_mode(2, address_mode)
-    tex_ref.set_filter_mode(filter_mode)
-
-    if not use_normalized_coordinates:
-        tex_ref.set_flags(tex_ref.get_flags() & ~cuda.TRSF_NORMALIZED_COORDINATES)
-
-    if not read_as_integer:
-        tex_ref.set_flags(tex_ref.get_flags() & ~cuda.TRSF_READ_AS_INTEGER)

pystencils/include/PyStencilsField.h

-#pragma once
-
-extern "C++" {
-#ifdef __CUDA_ARCH__
-template <typename DTYPE_T, std::size_t DIMENSION> struct PyStencilsField {
-  DTYPE_T *data;
-  DTYPE_T shape[DIMENSION];
-  DTYPE_T stride[DIMENSION];
-};
-#else
-#include <array>
-
-template <typename DTYPE_T, std::size_t DIMENSION> struct PyStencilsField {
-  DTYPE_T *data;
-  std::array<DTYPE_T, DIMENSION> shape;
-  std::array<DTYPE_T, DIMENSION> stride;
-};
-#endif
-}

pystencils/include/opencl_stdint.h

-#ifndef OPENCL_STDINT
-#define OPENCL_STDINT
-
-typedef unsigned int      uint_t;
-
-typedef signed char       int8_t;
-typedef signed short      int16_t;
-typedef signed int        int32_t;
-typedef signed long int   int64_t;
-typedef unsigned char     uint8_t;
-typedef unsigned short    uint16_t;
-typedef unsigned int      uint32_t;
-typedef unsigned long int uint64_t;
-
-#endif

pystencils/interpolation_astnodes.py

-# -*- coding: utf-8 -*-
-#
-# Copyright © 2019 Stephan Seitz <stephan.seitz@fau.de>
-#
-# Distributed under terms of the GPLv3 license.
-
-"""
-
-"""
-
-import hashlib
-import itertools
-from enum import Enum
-from typing import Set
-
-import sympy as sp
-from sympy.core.cache import cacheit
-
-import pystencils
-from pystencils.astnodes import Node
-from pystencils.data_types import TypedSymbol, cast_func, create_type
-
-try:
-    import pycuda.driver
-except Exception:
-    pass
-
-_hash = hashlib.md5
-
-
-class InterpolationMode(str, Enum):
-    NEAREST_NEIGHBOR = "nearest_neighbour"
-    NN = NEAREST_NEIGHBOR
-    LINEAR = "linear"
-    CUBIC_SPLINE = "cubic_spline"
-
-
-class _InterpolationSymbol(TypedSymbol):
-
-    def __new__(cls, name, field, interpolator):
-        obj = cls.__xnew_cached_(cls, name, field, interpolator)
-        return obj
-
-    def __new_stage2__(cls, name, field, interpolator):
-        obj = super().__xnew__(cls, name, 'dummy_symbol_carrying_field' + field.name)
-        obj.field = field
-        obj.interpolator = interpolator
-        return obj
-
-    def __getnewargs__(self):
-        return self.name, self.field, self.interpolator
-
-    def __getnewargs_ex__(self):
-        return (self.name, self.field, self.interpolator), {}
-
-    # noinspection SpellCheckingInspection
-    __xnew__ = staticmethod(__new_stage2__)
-    # noinspection SpellCheckingInspection
-    __xnew_cached_ = staticmethod(cacheit(__new_stage2__))
-
-
-class Interpolator(object):
-    """
-    Implements non-integer accesses on fields using linear interpolation.
-
-    On GPU, this interpolator can be implemented by a :class:`.TextureCachedField` for hardware acceleration.
-
-    Address modes are different boundary handlings possible choices are like for CUDA textures
-
-        **CLAMP**
-
-        The signal c[k] is continued outside k=0,...,M-1 so that c[k] = c[0] for k < 0, and c[k] = c[M-1] for k >= M.
-
-        **BORDER**
-
-        The signal c[k] is continued outside k=0,...,M-1 so that c[k] = 0 for k < 0and for k >= M.
-
-        Now, to describe the last two address modes, we are forced to consider normalized coordinates,
-        so that the 1D input signal samples are assumed to be c[k / M], with k=0,...,M-1.
-
-        **WRAP**
-
-        The signal c[k / M] is continued outside k=0,...,M-1 so that it is periodic with period equal to M.
-        In other words, c[(k + p * M) / M] = c[k / M] for any (positive, negative or vanishing) integer p.
-
-        **MIRROR**
-
-        The signal c[k / M] is continued outside k=0,...,M-1 so that it is periodic with period equal to 2 * M - 2.
-        In other words, c[l / M] = c[k / M] for any l and k such that (l + k)mod(2 * M - 2) = 0.
-
-    Explanations from https://stackoverflow.com/questions/19020963/the-different-addressing-modes-of-cuda-textures
-    """
-
-    required_global_declarations = []
-
-    def __init__(self,
-                 parent_field,
-                 interpolation_mode: InterpolationMode,
-                 address_mode='BORDER',
-                 use_normalized_coordinates=False,
-                 allow_textures=True):
-        super().__init__()
-
-        self.field = parent_field
-        self.field.field_type = pystencils.field.FieldType.CUSTOM
-        self.address_mode = address_mode
-        self.use_normalized_coordinates = use_normalized_coordinates
-        self.interpolation_mode = interpolation_mode
-        self.hash_str = hashlib.md5(
-            f'{self.field}_{address_mode}_{self.field.dtype}_{interpolation_mode}'.encode()).hexdigest()
-        self.symbol = _InterpolationSymbol(str(self), parent_field, self)
-        self.allow_textures = allow_textures
-
-    @property
-    def ndim(self):
-        return self.field.ndim
-
-    @property
-    def _hashable_contents(self):
-        return (str(self.address_mode),
-                str(type(self)),
-                self.hash_str,
-                self.use_normalized_coordinates)
-
-    def at(self, offset):
-        return InterpolatorAccess(self.symbol, *[sp.S(o) for o in offset])
-
-    def __getitem__(self, offset):
-        return InterpolatorAccess(self.symbol, *[sp.S(o) for o in offset])
-
-    def __str__(self):
-        return f'{self.field.name}_interpolator_{self.reproducible_hash}'
-
-    def __repr__(self):
-        return self.__str__()
-
-    def __hash__(self):
-        return hash(self._hashable_contents)
-
-    def __eq__(self, other):
-        return hash(self) == hash(other)
-
-    @property
-    def reproducible_hash(self):
-        return _hash(str(self._hashable_contents).encode()).hexdigest()
-
-
-class LinearInterpolator(Interpolator):
-
-    def __init__(self,
-                 parent_field: pystencils.Field,
-                 address_mode='BORDER',
-                 use_normalized_coordinates=False):
-        super().__init__(parent_field,
-                         InterpolationMode.LINEAR,
-                         address_mode,
-                         use_normalized_coordinates)
-
-
-class NearestNeightborInterpolator(Interpolator):
-
-    def __init__(self,
-                 parent_field: pystencils.Field,
-                 address_mode='BORDER',
-                 use_normalized_coordinates=False):
-        super().__init__(parent_field,
-                         InterpolationMode.NN,
-                         address_mode,
-                         use_normalized_coordinates)
-
-
-class InterpolatorAccess(TypedSymbol):
-    def __new__(cls, field, *offsets):
-        obj = InterpolatorAccess.__xnew_cached_(cls, field, *offsets)
-        return obj
-
-    def __new_stage2__(cls, symbol, *offsets):
-        assert offsets is not None
-        obj = super().__xnew__(cls, '%s_interpolator_%s' %
-                               (symbol.field.name, _hash(str(tuple(offsets)).encode()).hexdigest()),
-                               symbol.field.dtype)
-        obj.offsets = offsets
-        obj.symbol = symbol
-        obj.field = symbol.field
-        obj.interpolator = symbol.interpolator
-        return obj
-
-    def _hashable_contents(self):
-        return super()._hashable_content() + ((self.symbol, self.field, tuple(self.offsets), self.symbol.interpolator))
-
-    def __str__(self):
-        return f"{self.field.name}_interpolator({', '.join(str(o) for o in self.offsets)})"
-
-    def __repr__(self):
-        return self.__str__()
-
-    def _latex(self, printer, *_):
-        n = self.field.latex_name if self.field.latex_name else self.field.name
-        foo = ", ".join(str(printer.doprint(o)) for o in self.offsets)
-        return f'{n}_{{interpolator}}\\left({foo}\\right)'
-
-    @property
-    def ndim(self):
-        return len(self.offsets)
-
-    @property
-    def is_texture(self):
-        return isinstance(self.interpolator, TextureCachedField)
-
-    def atoms(self, *types):
-        if self.offsets:
-            offsets = set(o for o in self.offsets if isinstance(o, types))
-            if isinstance(self, *types):
-                offsets.update([self])
-            for o in self.offsets:
-                if hasattr(o, 'atoms'):
-                    offsets.update(set(o.atoms(*types)))
-            return offsets
-        else:
-            return set()
-
-    def neighbor(self, coord_id, offset):
-        offset_list = list(self.offsets)
-        offset_list[coord_id] += offset
-        return self.interpolator.at(tuple(offset_list))
-
-    @property
-    def free_symbols(self):
-        symbols = set()
-        if self.offsets is not None:
-            for o in self.offsets:
-                if hasattr(o, 'free_symbols'):
-                    symbols.update(set(o.free_symbols))
-                # if hasattr(o, 'atoms'):
-                    # symbols.update(set(o.atoms(sp.Symbol)))
-
-        return symbols
-
-    @property
-    def required_global_declarations(self):
-        required_global_declarations = self.symbol.interpolator.required_global_declarations
-        if required_global_declarations:
-            required_global_declarations[0]._symbols_defined.add(self)
-        return required_global_declarations
-
-    @property
-    def args(self):
-        return [self.symbol, *self.offsets]
-
-    @property
-    def symbols_defined(self) -> Set[sp.Symbol]:
-        return {self}
-
-    @property
-    def interpolation_mode(self):
-        return self.interpolator.interpolation_mode
-
-    @property
-    def _diff_interpolation_vec(self):
-        return sp.Matrix([DiffInterpolatorAccess(self.symbol, i, *self.offsets)
-                          for i in range(len(self.offsets))])
-
-    def diff(self, *symbols, **kwargs):
-        if symbols == (self,):
-            return 1
-        rtn = self._diff_interpolation_vec.T * sp.Matrix(self.offsets).diff(*symbols, **kwargs)
-        if rtn.shape == (1, 1):
-            rtn = rtn[0, 0]
-        return rtn
-
-    def implementation_with_stencils(self):
-        field = self.field
-
-        default_int_type = create_type('int64')
-        use_textures = isinstance(self.interpolator, TextureCachedField)
-        if use_textures:
-            def absolute_access(x, _):
-                return self.symbol.interpolator.at((o for o in x))
-        else:
-            absolute_access = field.absolute_access
-
-        sum = [0, ] * (field.shape[0] if field.index_dimensions else 1)
-
-        offsets = self.offsets
-        rounding_functions = (sp.floor, lambda x: sp.floor(x) + 1)
-
-        for channel_idx in range(field.shape[0] if field.index_dimensions else 1):
-            if self.interpolation_mode == InterpolationMode.NN:
-                if use_textures:
-                    sum[channel_idx] = self
-                else:
-                    sum[channel_idx] = absolute_access([sp.floor(i + 0.5) for i in offsets], channel_idx)
-
-            elif self.interpolation_mode == InterpolationMode.LINEAR:
-                # TODO optimization: implement via lerp: https://devblogs.nvidia.com/lerp-faster-cuda/
-                for c in itertools.product(rounding_functions, repeat=field.spatial_dimensions):
-                    weight = sp.Mul(*[1 - sp.Abs(f(offset) - offset) for (f, offset) in zip(c, offsets)])
-                    index = [f(offset) for (f, offset) in zip(c, offsets)]
-                    # Hardware boundary handling on GPU
-                    if use_textures:
-                        weight = sp.Mul(*[1 - sp.Abs(f(offset) - offset) for (f, offset) in zip(c, offsets)])
-                        sum[channel_idx] += \
-                            weight * absolute_access(index, channel_idx if field.index_dimensions else ())
-                    # else boundary handling using software
-                    elif str(self.interpolator.address_mode).lower() == 'border':
-                        is_inside_field = sp.And(
-                            *itertools.chain([i >= 0 for i in index],
-                                             [idx < field.shape[dim] for (dim, idx) in enumerate(index)]))
-                        index = [cast_func(i, default_int_type) for i in index]
-                        sum[channel_idx] += sp.Piecewise(
-                            (weight * absolute_access(index, channel_idx if field.index_dimensions else ()),
-                                is_inside_field),
-                            (sp.simplify(0), True)
-                        )
-                    elif str(self.interpolator.address_mode).lower() == 'clamp':
-                        index = [sp.Min(sp.Max(0, cast_func(i, default_int_type)), field.spatial_shape[dim] - 1)
-                                 for (dim, i) in enumerate(index)]
-                        sum[channel_idx] += weight * \
-                            absolute_access(index, channel_idx if field.index_dimensions else ())
-                    elif str(self.interpolator.address_mode).lower() == 'wrap':
-                        index = [sp.Mod(cast_func(i, default_int_type), field.shape[dim] - 1)
-                                 for (dim, i) in enumerate(index)]
-                        index = [cast_func(sp.Piecewise((i, i > 0),
-                                                        (sp.Abs(cast_func(field.shape[dim] - 1 + i, default_int_type)),
-                                                         True)), default_int_type)
-                                 for (dim, i) in enumerate(index)]
-                        sum[channel_idx] += weight * \
-                            absolute_access(index, channel_idx if field.index_dimensions else ())
-                        # sum[channel_idx] = 0
-                    elif str(self.interpolator.address_mode).lower() == 'mirror':
-                        def triangle_fun(x, half_period):
-                            saw_tooth = cast_func(sp.Abs(cast_func(x, 'int32')), 'int32') % (
-                                cast_func(2 * half_period, create_type('int32')))
-                            return sp.Piecewise((saw_tooth, saw_tooth < half_period),
-                                                (2 * half_period - 1 - saw_tooth, True))
-                        index = [cast_func(triangle_fun(i, field.shape[dim]),
-                                           default_int_type) for (dim, i) in enumerate(index)]
-                        sum[channel_idx] += weight * \
-                            absolute_access(index, channel_idx if field.index_dimensions else ())
-                    else:
-                        raise NotImplementedError()
-            elif self.interpolation_mode == InterpolationMode.CUBIC_SPLINE:
-                raise NotImplementedError("only works with HW interpolation for float32")
-
-            sum = [sp.factor(s) for s in sum]
-
-            if field.index_dimensions:
-                return sp.Matrix(sum)
-            else:
-                return sum[0]
-
-    # noinspection SpellCheckingInspection
-    __xnew__ = staticmethod(__new_stage2__)
-    # noinspection SpellCheckingInspection
-    __xnew_cached_ = staticmethod(cacheit(__new_stage2__))
-
-    def __getnewargs__(self):
-        return (self.symbol, *self.offsets)
-
-    def __getnewargs_ex__(self):
-        return (self.symbol, *self.offsets), {}
-
-
-class DiffInterpolatorAccess(InterpolatorAccess):
-    def __new__(cls, symbol, diff_coordinate_idx, *offsets):
-        if symbol.interpolator.interpolation_mode == InterpolationMode.LINEAR:
-            from pystencils.fd import Diff, Discretization2ndOrder
-            return Discretization2ndOrder(1)(Diff(symbol.interpolator.at(offsets), diff_coordinate_idx))
-        obj = DiffInterpolatorAccess.__xnew_cached_(cls, symbol, diff_coordinate_idx, *offsets)
-        return obj
-
-    def __new_stage2__(self, symbol: sp.Symbol, diff_coordinate_idx, *offsets):
-        assert offsets is not None
-        obj = super().__xnew__(self, symbol, *offsets)
-        obj.diff_coordinate_idx = diff_coordinate_idx
-        return obj
-
-    def __hash__(self):
-        return hash((self.symbol, self.field, self.diff_coordinate_idx, tuple(self.offsets), self.interpolator))
-
-    def __str__(self):
-        return '%s_diff%i_interpolator(%s)' % (self.field.name, self.diff_coordinate_idx,
-                                               ', '.join(str(o) for o in self.offsets))
-
-    def __repr__(self):
-        return str(self)
-
-    @property
-    def args(self):
-        return [self.symbol, self.diff_coordinate_idx, *self.offsets]
-
-    @property
-    def symbols_defined(self) -> Set[sp.Symbol]:
-        return {self}
-
-    @property
-    def interpolation_mode(self):
-        return self.interpolator.interpolation_mode
-
-    # noinspection SpellCheckingInspection
-    __xnew__ = staticmethod(__new_stage2__)
-    # noinspection SpellCheckingInspection
-    __xnew_cached_ = staticmethod(cacheit(__new_stage2__))
-
-    def __getnewargs__(self):
-        return (self.symbol, self.diff_coordinate_idx, *self.offsets)
-
-    def __getnewargs_ex__(self):
-        return (self.symbol, self.diff_coordinate_idx, *self.offsets), {}
-
-
-##########################################################################################
-# GPU-specific fast specializations (for precision GPUs can also use above nodes/symbols #
-##########################################################################################
-
-
-class TextureCachedField(Interpolator):
-
-    def __init__(self, parent_field,
-                 address_mode=None,
-                 filter_mode=None,
-                 interpolation_mode: InterpolationMode = InterpolationMode.LINEAR,
-                 use_normalized_coordinates=False,
-                 read_as_integer=False
-                 ):
-        super().__init__(parent_field, interpolation_mode, address_mode, use_normalized_coordinates)
-
-        if address_mode is None:
-            address_mode = 'border'
-        if filter_mode is None:
-            filter_mode = pycuda.driver.filter_mode.LINEAR
-
-        self.read_as_integer = read_as_integer
-        self.required_global_declarations = [TextureDeclaration(self)]
-
-    @property
-    def ndim(self):
-        return self.field.ndim
-
-    @classmethod
-    def from_interpolator(cls, interpolator: LinearInterpolator):
-        if (isinstance(interpolator, cls)
-                or (hasattr(interpolator, 'allow_textures') and not interpolator.allow_textures)):
-            return interpolator
-        obj = cls(interpolator.field, interpolator.address_mode, interpolation_mode=interpolator.interpolation_mode)
-        return obj
-
-    def __str__(self):
-        return f'{self.field.name}_texture_{self.reproducible_hash}'
-
-    def __repr__(self):
-        return self.__str__()
-
-    @property
-    def reproducible_hash(self):
-        return _hash(str(self._hashable_contents).encode()).hexdigest()
-
-
-class TextureDeclaration(Node):
-    """
-    A global declaration of a texture. Visible both for device and host code.
-
-    .. code:: cpp
-
-        // This Node represents the following global declaration
-        texture<float, cudaTextureType2D, cudaReadModeElementType> x_texture_5acc9fced7b0dc3e;
-
-        __device__ kernel(...) {
-            // kernel acceses x_texture_5acc9fced7b0dc3e with tex2d(...)
-        }
-
-        __host__ launch_kernel(...) {
-            // Host needs to bind the texture
-            cudaBindTexture(0, x_texture_5acc9fced7b0dc3e, buffer, N*sizeof(float));
-        }
-
-    This has been deprecated by CUDA in favor of :class:`.TextureObject`.
-    But texture objects are not yet supported by PyCUDA (https://github.com/inducer/pycuda/pull/174)
-    """
-
-    def __init__(self, parent_texture):
-        self.texture = parent_texture
-        self._symbols_defined = {self.texture.symbol}
-
-    @property
-    def symbols_defined(self) -> Set[sp.Symbol]:
-        return self._symbols_defined
-
-    @property
-    def args(self) -> Set[sp.Symbol]:
-        return set()
-
-    @property
-    def headers(self):
-        headers = ['"pycuda-helpers.hpp"']
-        if self.texture.interpolation_mode == InterpolationMode.CUBIC_SPLINE:
-            headers.append('"cubicTex%iD.cu"' % self.texture.ndim)
-        return headers
-
-    def __str__(self):
-        from pystencils.backends.cuda_backend import CudaBackend
-        return CudaBackend()(self)
-
-    def __repr__(self):
-        return str(self)
-
-
-class TextureObject(TextureDeclaration):
-    """
-    A CUDA texture object. Opposed to :class:`.TextureDeclaration` it is not declared globally but
-    used as a function argument for the kernel call.
-
-    Like :class:`.TextureDeclaration` it defines :class:`.TextureAccess` symbols.
-    Just the printing representation is a bit different.
-    """
-    pass
-
-
-def dtype_supports_textures(dtype):
-    """
-    Returns whether CUDA natively supports texture fetches with this numpy dtype.
-
-    The maximum word size for a texture fetch is four bytes.
-
-    With this trick also larger dtypes can be fetched:
-    https://github.com/inducer/pycuda/blob/master/pycuda/cuda/pycuda-helpers.hpp
-
-    """
-    if hasattr(dtype, 'numpy_dtype'):
-        dtype = dtype.numpy_dtype
-
-    if isinstance(dtype, type):
-        return dtype().itemsize <= 4
-
-    return dtype.itemsize <= 4

pystencils/kerncraft_coupling/__init__.py

-from .generate_benchmark import generate_benchmark, run_c_benchmark
-from .kerncraft_interface import KerncraftParameters, PyStencilsKerncraftKernel
-
-__all__ = ['PyStencilsKerncraftKernel', 'KerncraftParameters', 'generate_benchmark', 'run_c_benchmark']

pystencils/kerncraft_coupling/generate_benchmark.py

-import subprocess
-import warnings
-import tempfile
-from pathlib import Path
-
-from jinja2 import Environment, PackageLoader, StrictUndefined
-
-from pystencils.astnodes import PragmaBlock
-from pystencils.backends.cbackend import generate_c, get_headers
-from pystencils.cpu.cpujit import get_compiler_config, run_compile_step
-from pystencils.data_types import get_base_type
-from pystencils.include import get_pystencils_include_path
-from pystencils.integer_functions import modulo_ceil
-from pystencils.sympyextensions import prod
-
-import numpy as np
-
-
-def generate_benchmark(ast, likwid=False, openmp=False, timing=False):
-    """Return C code of a benchmark program for the given kernel.
-
-    Args:
-        ast: the pystencils AST object as returned by create_kernel
-        likwid: if True likwid markers are added to the code
-        openmp: relevant only if likwid=True, to generated correct likwid initialization code
-        timing: add timing output to the code, prints time per iteration to stdout
-
-    Returns:
-        C code as string
-    """
-    accessed_fields = {f.name: f for f in ast.fields_accessed}
-    constants = []
-    fields = []
-    call_parameters = []
-    for p in ast.get_parameters():
-        if not p.is_field_parameter:
-            constants.append((p.symbol.name, str(p.symbol.dtype)))
-            call_parameters.append(p.symbol.name)
-        else:
-            assert p.is_field_pointer, "Benchmark implemented only for kernels with fixed loop size"
-            field = accessed_fields[p.field_name]
-            dtype = str(get_base_type(p.symbol.dtype))
-            np_dtype = get_base_type(p.symbol.dtype).numpy_dtype
-            size_data_type = np_dtype.itemsize
-
-            dim0_size = field.shape[-1]
-            dim1_size = np.prod(field.shape[:-1])
-            elements = prod(field.shape)
-
-            if ast.instruction_set:
-                align = ast.instruction_set['width'] * size_data_type
-                padding_elements = modulo_ceil(dim0_size, ast.instruction_set['width']) - dim0_size
-                padding_bytes = padding_elements * size_data_type
-                ghost_layers = max(max(ast.ghost_layers))
-
-                size = dim1_size * padding_bytes + np.prod(field.shape) * size_data_type
-
-                assert align % np_dtype.itemsize == 0
-                offset = ((dim0_size + padding_elements + ghost_layers) % ast.instruction_set['width']) * size_data_type
-
-                fields.append((p.field_name, dtype, elements, size, offset, align))
-                call_parameters.append(p.field_name)
-            else:
-                size = elements * size_data_type
-                fields.append((p.field_name, dtype, elements, size, 0, 0))
-                call_parameters.append(p.field_name)
-
-    header_list = get_headers(ast)
-    includes = "\n".join(["#include %s" % (include_file,) for include_file in header_list])
-
-    # Strip "#pragma omp parallel" from within kernel, because main function takes care of that
-    # when likwid and openmp are enabled
-    if likwid and openmp:
-        if len(ast.body.args) > 0 and isinstance(ast.body.args[0], PragmaBlock):
-            ast.body.args[0].pragma_line = ''
-
-    jinja_context = {
-        'likwid': likwid,
-        'openmp': openmp,
-        'kernel_code': generate_c(ast, dialect='c'),
-        'kernelName': ast.function_name,
-        'fields': fields,
-        'constants': constants,
-        'call_argument_list': ",".join(call_parameters),
-        'includes': includes,
-        'timing': timing,
-    }
-
-    env = Environment(loader=PackageLoader('pystencils.kerncraft_coupling'), undefined=StrictUndefined)
-
-    return env.get_template('benchmark.c').render(**jinja_context)
-
-
-def run_c_benchmark(ast, inner_iterations, outer_iterations=3, path=None):
-    """Runs the given kernel with outer loop in C
-
-    Args:
-        ast: pystencils ast which is used to compile the benchmark file
-        inner_iterations: timings are recorded around this many iterations
-        outer_iterations: number of timings recorded
-        path: path where the benchmark file is stored. If None a tmp folder is created
-
-    Returns:
-        list of times per iterations for each outer iteration
-    """
-    import kerncraft
-
-    benchmark_code = generate_benchmark(ast, timing=True)
-
-    if path is None:
-        path = tempfile.mkdtemp()
-
-    if isinstance(path, str):
-        path = Path(path)
-
-    with open(path / 'bench.c', 'w') as f:
-        f.write(benchmark_code)
-
-    kerncraft_path = Path(kerncraft.__file__).parent
-
-    extra_flags = ['-I' + get_pystencils_include_path(),
-                   '-I' + str(kerncraft_path / 'headers')]
-
-    compiler_config = get_compiler_config()
-    compile_cmd = [compiler_config['command']] + compiler_config['flags'].split()
-    compile_cmd += [*extra_flags,
-                    str(kerncraft_path / 'headers' / 'timing.c'),
-                    str(kerncraft_path / 'headers' / 'dummy.c'),
-                    str(path / 'bench.c'),
-                    '-o', str(path / 'bench'),
-                    ]
-    run_compile_step(compile_cmd)
-
-    time_pre_estimation_per_iteration = float(subprocess.check_output(['./' / path / 'bench', str(10)]))
-    benchmark_time_limit = 20
-    if benchmark_time_limit / time_pre_estimation_per_iteration < inner_iterations:
-        warn = (f"A benchmark run with {inner_iterations} inner_iterations will probably take longer than "
-                f"{benchmark_time_limit} seconds for this kernel")
-        warnings.warn(warn)
-
-    results = []
-    for _ in range(outer_iterations):
-        benchmark_time = float(subprocess.check_output(['./' / path / 'bench', str(inner_iterations)]))
-        results.append(benchmark_time)
-    return results