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:
 
 
 

pre-push

-#!/usr/bin/env bash
-
-# Checks run before every push
-# has to be copied to .git/hooks
-
-echo "Running pre-push hook"
-echo "Running flake8 check"
-flake8 --append-config=.flake8 pystencils
-
-# $? stores exit value of the last command
-if [ $? -ne 0 ]; then
- echo "flake8 failed"
- exit 1
-fi
-
-python3 setup.py quicktest
-
-# $? stores exit value of the last command
-if [ $? -ne 0 ]; then
- echo "quicktest failed"
- exit 1
-fi

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