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 os
-import subprocess
-
-from jinja2 import Template
-
-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.sympyextensions import prod
-
-benchmark_template = Template("""
-#include "kerncraft.h"
-#include <stdlib.h>
-#include <stdint.h>
-#include <stdbool.h>
-#include <math.h>
-#include <stdio.h>
-
-{{ includes }}
-
-{%- if likwid %}
-#include <likwid.h>
-{%- endif %}
-
-#define RESTRICT __restrict__
-#define FUNC_PREFIX
-void dummy(void *);
-void timing(double* wcTime, double* cpuTime);
-extern int var_false;
-
-
-{{kernel_code}}
-
-
-int main(int argc, char **argv)
-{
-  {%- if likwid %}
-  likwid_markerInit();
-  {%- endif %}
-
-  {%- for field_name, dataType, size in fields %}
-
-  // Initialization {{field_name}}
-  double * {{field_name}} = (double *) aligned_malloc(sizeof({{dataType}}) * {{size}}, 64);
-  for (unsigned long long i = 0; i < {{size}}; ++i)
-    {{field_name}}[i] = 0.23;
-
-  if(var_false)
-    dummy({{field_name}});
-
-  {%- endfor %}
-
-
-
-  {%- for constantName, dataType in constants %}
-
-  // Constant {{constantName}}
-  {{dataType}} {{constantName}};
-  {{constantName}} = 0.23;
-  if(var_false)
-      dummy(& {{constantName}});
-
-  {%- endfor %}
-
-  {%- if likwid and openmp %}
-  #pragma omp parallel
-  {
-  likwid_markerRegisterRegion("loop");
-  #pragma omp barrier
-  {%- elif likwid %}
-  likwid_markerRegisterRegion("loop");
-  {%- endif %}
-
-  for(int warmup = 1; warmup >= 0; --warmup) {
-    int repeat = 2;
-    if(warmup == 0) {
-      repeat = atoi(argv[1]);
-      {%- if likwid %}
-      likwid_markerStartRegion("loop");
-      {%- endif %}
-    }
-    
-    {%- if timing %}
-    double wcStartTime, cpuStartTime, wcEndTime, cpuEndTime;
-    timing(&wcStartTime, &cpuStartTime);
-    {%- endif %}
-    
-    for (; repeat > 0; --repeat)
-    {
-      {{kernelName}}({{call_argument_list}});
-
-      // Dummy calls
-      {%- for field_name, dataType, size in fields %}
-      if(var_false) dummy((void*){{field_name}});
-      {%- endfor %}
-      {%- for constantName, dataType in constants %}
-      if(var_false) dummy((void*)&{{constantName}});
-      {%- endfor %}
-    }
-    {%- if timing %}
-    timing(&wcEndTime, &cpuEndTime);
-    if( warmup == 0)
-        printf("%e\\n", (wcEndTime - wcStartTime) / atoi(argv[1]) );
-    {%- endif %}
-
-  }
-
-  {%- if likwid %}
-  likwid_markerStopRegion("loop");
-  {%- if openmp %}
-  }
-  {%- endif %}
-  {%- endif %}
-
-  {%- if likwid %}
-  likwid_markerClose();
-  {%- endif %}
-}
-""")
-
-
-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))
-            fields.append((p.field_name, dtype, prod(field.shape)))
-            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 = ''
-
-    args = {
-        '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,
-    }
-    return benchmark_template.render(**args)
-
-
-def run_c_benchmark(ast, inner_iterations, outer_iterations=3):
-    """Runs the given kernel with outer loop in C
-
-    Args:
-        ast:
-        inner_iterations: timings are recorded around this many iterations
-        outer_iterations: number of timings recorded
-
-    Returns:
-        list of times per iterations for each outer iteration
-    """
-    import kerncraft
-
-    benchmark_code = generate_benchmark(ast, timing=True)
-    with open('bench.c', 'w') as f:
-        f.write(benchmark_code)
-
-    kerncraft_path = os.path.dirname(kerncraft.__file__)
-
-    extra_flags = ['-I' + get_pystencils_include_path(),
-                   '-I' + os.path.join(kerncraft_path, 'headers')]
-
-    compiler_config = get_compiler_config()
-    compile_cmd = [compiler_config['command']] + compiler_config['flags'].split()
-    compile_cmd += [*extra_flags,
-                    os.path.join(kerncraft_path, 'headers', 'timing.c'),
-                    os.path.join(kerncraft_path, 'headers', 'dummy.c'),
-                    'bench.c',
-                    '-o', 'bench',
-                    ]
-    run_compile_step(compile_cmd)
-
-    results = []
-    for _ in range(outer_iterations):
-        benchmark_time = float(subprocess.check_output(['./bench', str(inner_iterations)]))
-        results.append(benchmark_time)
-    return results

pystencils/kerncraft_coupling/kerncraft_interface.py

-import warnings
-from collections import defaultdict
-from tempfile import TemporaryDirectory
-from typing import Optional
-
-import kerncraft
-import sympy as sp
-from kerncraft.kerncraft import KernelCode
-from kerncraft.machinemodel import MachineModel
-
-from pystencils.astnodes import (
-    KernelFunction, LoopOverCoordinate, ResolvedFieldAccess, SympyAssignment)
-from pystencils.field import get_layout_from_strides
-from pystencils.kerncraft_coupling.generate_benchmark import generate_benchmark
-from pystencils.sympyextensions import count_operations_in_ast
-from pystencils.transformations import filtered_tree_iteration
-from pystencils.utils import DotDict
-
-
-class PyStencilsKerncraftKernel(KernelCode):
-    """
-    Implementation of kerncraft's kernel interface for pystencils CPU kernels.
-    Analyses a list of equations assuming they will be executed on a CPU
-    """
-    LIKWID_BASE = '/usr/local/likwid'
-
-    def __init__(self, ast: KernelFunction, machine: Optional[MachineModel] = None,
-                 assumed_layout='SoA', debug_print=False, filename=None):
-        """Create a kerncraft kernel using a pystencils AST
-
-        Args:
-            ast: pystencils ast
-            machine: kerncraft machine model - specify this if kernel needs to be compiled
-            assumed_layout: either 'SoA' or 'AoS' - if fields have symbolic sizes the layout of the index
-                    coordinates is not known. In this case either a structures of array (SoA) or
-                    array of structures (AoS) layout is assumed
-        """
-        kerncraft.kernel.Kernel.__init__(self, machine)
-
-        # Initialize state
-        self.asm_block = None
-        self._filename = filename
-
-        self.kernel_ast = ast
-        self.temporary_dir = TemporaryDirectory()
-        self._keep_intermediates = debug_print
-
-        # Loops
-        inner_loops = [l for l in filtered_tree_iteration(ast, LoopOverCoordinate, stop_type=SympyAssignment)
-                       if l.is_innermost_loop]
-        if len(inner_loops) == 0:
-            raise ValueError("No loop found in pystencils AST")
-        else:
-            if len(inner_loops) > 1:
-                warnings.warn("pystencils AST contains multiple inner loops. "
-                              "Only one can be analyzed - choosing first one")
-            inner_loop = inner_loops[0]
-
-        self._loop_stack = []
-        cur_node = inner_loop
-        while cur_node is not None:
-            if isinstance(cur_node, LoopOverCoordinate):
-                loop_counter_sym = cur_node.loop_counter_symbol
-                loop_info = (loop_counter_sym.name, cur_node.start, cur_node.stop, 1)
-                # If the correct step were to be provided, all access within that step length will
-                # also need to be passed to kerncraft: cur_node.step)
-                self._loop_stack.append(loop_info)
-            cur_node = cur_node.parent
-        self._loop_stack = list(reversed(self._loop_stack))
-
-        # Data sources & destinations
-        self.sources = defaultdict(list)
-        self.destinations = defaultdict(list)
-
-        def get_layout_tuple(f):
-            if f.has_fixed_shape:
-                return get_layout_from_strides(f.strides)
-            else:
-                layout_list = list(f.layout)
-                for _ in range(f.index_dimensions):
-                    layout_list.insert(0 if assumed_layout == 'SoA' else -1, max(layout_list) + 1)
-                return layout_list
-
-        reads, writes = search_resolved_field_accesses_in_ast(inner_loop)
-        for accesses, target_dict in [(reads, self.sources), (writes, self.destinations)]:
-            for fa in accesses:
-                coord = [sp.Symbol(LoopOverCoordinate.get_loop_counter_name(i), positive=True, integer=True) + off
-                         for i, off in enumerate(fa.offsets)]
-                coord += list(fa.idx_coordinate_values)
-                layout = get_layout_tuple(fa.field)
-                permuted_coord = [sp.sympify(coord[i]) for i in layout]
-                target_dict[fa.field.name].append(permuted_coord)
-
-        # Variables (arrays)
-        fields_accessed = ast.fields_accessed
-        for field in fields_accessed:
-            layout = get_layout_tuple(field)
-            permuted_shape = list(field.shape[i] for i in layout)
-            self.set_variable(field.name, str(field.dtype), tuple(permuted_shape))
-
-        # Scalars may be safely ignored
-        # for param in ast.get_parameters():
-        #     if not param.is_field_parameter:
-        #         # self.set_variable(param.symbol.name, str(param.symbol.dtype), None)
-        #         self.sources[param.symbol.name] = [None]
-
-        # data type
-        self.datatype = list(self.variables.values())[0][0]
-
-        # flops
-        operation_count = count_operations_in_ast(inner_loop)
-        self._flops = {
-            '+': operation_count['adds'],
-            '*': operation_count['muls'],
-            '/': operation_count['divs'],
-        }
-        for k in [k for k, v in self._flops.items() if v == 0]:
-            del self._flops[k]
-        self.check()
-
-        if debug_print:
-            from pprint import pprint
-            print("-----------------------------  Loop Stack --------------------------")
-            pprint(self._loop_stack)
-            print("-----------------------------  Sources -----------------------------")
-            pprint(self.sources)
-            print("-----------------------------  Destinations ------------------------")
-            pprint(self.destinations)
-            print("-----------------------------  FLOPS -------------------------------")
-            pprint(self._flops)
-
-    def as_code(self, type_='iaca', openmp=False, as_filename=False):
-        """
-        Generate and return compilable source code.
-
-        Args:
-            type_: can be iaca or likwid.
-            openmp: if true, openmp code will be generated
-            as_filename:
-        """
-        code = generate_benchmark(self.kernel_ast, likwid=type_ == 'likwid', openmp=openmp)
-        if as_filename:
-            fp, already_available = self._get_intermediate_file('kernel_{}.c'.format(type_),
-                                                                machine_and_compiler_dependent=False)
-            if not already_available:
-                fp.write(code)
-            return fp.name
-        else:
-            return code
-
-
-class KerncraftParameters(DotDict):
-    def __init__(self, **kwargs):
-        super(KerncraftParameters, self).__init__(**kwargs)
-        self['asm_block'] = 'auto'
-        self['asm_increment'] = 0
-        self['cores'] = 1
-        self['cache_predictor'] = 'SIM'
-        self['verbose'] = 0
-        self['pointer_increment'] = 'auto'
-        self['iterations'] = 10
-        self['unit'] = 'cy/CL'
-        self['ignore_warnings'] = True
-
-
-# ------------------------------------------- Helper functions ---------------------------------------------------------
-
-
-def search_resolved_field_accesses_in_ast(ast):
-    def visit(node, reads, writes):
-        if not isinstance(node, SympyAssignment):
-            for a in node.args:
-                visit(a, reads, writes)
-            return
-
-        for expr, accesses in [(node.lhs, writes), (node.rhs, reads)]:
-            accesses.update(expr.atoms(ResolvedFieldAccess))
-
-    read_accesses = set()
-    write_accesses = set()
-    visit(ast, read_accesses, write_accesses)
-    return read_accesses, write_accesses

pystencils/llvm/__init__.py

-from .kernelcreation import create_kernel
-from .llvmjit import make_python_function
-
-__all__ = ['create_kernel', 'make_python_function']

pystencils/llvm/control_flow.py

-import llvmlite.ir as ir
-
-
-class Loop(object):
-    def __init__(self, builder, start_val, stop_val, step_val=1, loop_name='loop', phi_name="_phi"):
-        self.builder = builder
-        self.start_val = start_val
-        self.stop_val = stop_val
-        self.step_val = step_val
-        self.loop_name = loop_name
-        self.phi_name = phi_name
-
-    def __enter__(self):
-        self.loop_end, self.after, phi = self._for_loop(self.start_val, self.stop_val, self.step_val, self.loop_name,
-                                                        self.phi_name)
-        return phi
-
-    def _for_loop(self, start_val, stop_val, step_val, loop_name, phi_name):
-        # TODO size of int??? unisgned???
-        integer = ir.IntType(64)
-
-        # Loop block
-        pre_loop_bb = self.builder.block
-        loop_bb = self.builder.append_basic_block(name='loop_' + loop_name)
-
-        self.builder.branch(loop_bb)
-
-        # Insert an explicit fall through from the current block to loop_bb
-        self.builder.position_at_start(loop_bb)
-
-        # Add phi
-        phi = self.builder.phi(integer, name=phi_name)
-        phi.add_incoming(start_val, pre_loop_bb)
-
-        loop_end_bb = self.builder.append_basic_block(name=loop_name + "_end_bb")
-        self.builder.position_at_start(loop_end_bb)
-
-        next_var = self.builder.add(phi, step_val, name=loop_name + '_next_it')
-        cond = self.builder.icmp_unsigned('<', next_var, stop_val, name=loop_name + "_cond")
-
-        after_bb = self.builder.append_basic_block(name=loop_name + "_after_bb")
-
-        self.builder.cbranch(cond, loop_bb, after_bb)
-        phi.add_incoming(next_var, loop_end_bb)
-
-        self.builder.position_at_end(loop_bb)
-
-        return loop_end_bb, after_bb, phi
-
-    def __exit__(self, exc_type, exc, exc_tb):
-        self.builder.branch(self.loop_end)
-        self.builder.position_at_end(self.after)

pystencils/llvm/kernelcreation.py

-from pystencils.llvm.llvmjit import make_python_function
-from pystencils.transformations import insert_casts
-
-
-def create_kernel(assignments, function_name="kernel", type_info=None, split_groups=(),
-                  iteration_slice=None, ghost_layers=None):
-    """
-    Creates an abstract syntax tree for a kernel function, by taking a list of update rules.
-
-    Loops are created according to the field accesses in the equations.
-
-    Args:
-    assignments: list of sympy equations, containing accesses to :class:`pystencils.field.Field`.
-                 Defining the update rules of the kernel
-    function_name: name of the generated function - only important if generated code is written out
-                   type_info: a map from symbol name to a C type specifier. If not specified all symbols are assumed to
-                   be of type 'double' except symbols which occur on the left hand side of equations where the
-                   right hand side is a sympy Boolean which are assumed to be 'bool' .
-    split_groups: Specification on how to split up inner loop into multiple loops. For details see
-                  transformation :func:`pystencils.transformation.split_inner_loop`
-    iteration_slice: if not None, iteration is done only over this slice of the field
-    ghost_layers: a sequence of pairs for each coordinate with lower and upper nr of ghost layers
-                 if None, the number of ghost layers is determined automatically and assumed to be equal for a
-                 all dimensions
-
-    :return: :class:`pystencils.ast.KernelFunction` node
-    """
-    from pystencils.cpu import create_kernel
-    code = create_kernel(assignments, function_name, type_info, split_groups, iteration_slice, ghost_layers)
-    code.body = insert_casts(code.body)
-    code._compile_function = make_python_function
-    code._backend = 'llvm'
-    return code

pystencils/rng.py

-import numpy as np
-import sympy as sp
-
-from pystencils import TypedSymbol
-from pystencils.astnodes import LoopOverCoordinate
-from pystencils.backends.cbackend import CustomCodeNode
-
-philox_two_doubles_call = """
-{result_symbols[0].dtype} {result_symbols[0].name};
-{result_symbols[1].dtype} {result_symbols[1].name};
-philox_double2({parameters}, {result_symbols[0].name}, {result_symbols[1].name});
-"""
-
-philox_four_floats_call = """
-{result_symbols[0].dtype} {result_symbols[0].name};
-{result_symbols[1].dtype} {result_symbols[1].name};
-{result_symbols[2].dtype} {result_symbols[2].name};
-{result_symbols[3].dtype} {result_symbols[3].name};
-philox_float4({parameters},
-              {result_symbols[0].name}, {result_symbols[1].name}, {result_symbols[2].name}, {result_symbols[3].name});
-
-"""
-
-
-def _get_philox_code(template, dialect, vector_instruction_set, time_step, offsets, keys, dim, result_symbols):
-    parameters = [time_step] + [LoopOverCoordinate.get_loop_counter_symbol(i) + offsets[i]
-                                for i in range(dim)] + list(keys)
-
-    while len(parameters) < 6:
-        parameters.append(0)
-    parameters = parameters[:6]
-
-    assert len(parameters) == 6
-
-    if dialect == 'cuda' or (dialect == 'c' and vector_instruction_set is None):
-        return template.format(parameters=', '.join(str(p) for p in parameters),
-                               result_symbols=result_symbols)
-    else:
-        raise NotImplementedError("Not yet implemented for this backend")
-
-
-class PhiloxTwoDoubles(CustomCodeNode):
-
-    def __init__(self, dim, time_step=TypedSymbol("time_step", np.uint32), offsets=(0, 0, 0), keys=(0, 0)):
-        self.result_symbols = tuple(TypedSymbol(sp.Dummy().name, np.float64) for _ in range(2))
-        symbols_read = [s for s in keys if isinstance(s, sp.Symbol)]
-        super().__init__("", symbols_read=symbols_read, symbols_defined=self.result_symbols)
-        self._time_step = time_step
-        self._offsets = offsets
-        self.headers = ['"philox_rand.h"']
-        self.keys = tuple(keys)
-        self._args = sp.sympify((dim, time_step, keys))
-        self._dim = dim
-
-    @property
-    def args(self):
-        return self._args
-
-    @property
-    def undefined_symbols(self):
-        result = {a for a in (self._time_step, *self._offsets, *self.keys) if isinstance(a, sp.Symbol)}
-        loop_counters = [LoopOverCoordinate.get_loop_counter_symbol(i)
-                         for i in range(self._dim)]
-        result.update(loop_counters)
-        return result
-
-    def fast_subs(self, _):
-        return self  # nothing to replace inside this node - would destroy intermediate "dummy" by re-creating them
-
-    def get_code(self, dialect, vector_instruction_set):
-        return _get_philox_code(philox_two_doubles_call, dialect, vector_instruction_set,
-                                self._time_step, self._offsets, self.keys, self._dim, self.result_symbols)
-
-    def __repr__(self):
-        return "{}, {} <- PhiloxRNG".format(*self.result_symbols)
-
-
-class PhiloxFourFloats(CustomCodeNode):
-
-    def __init__(self, dim, time_step=TypedSymbol("time_step", np.uint32), offsets=(0, 0, 0), keys=(0, 0)):
-        self.result_symbols = tuple(TypedSymbol(sp.Dummy().name, np.float32) for _ in range(4))
-        symbols_read = [s for s in keys if isinstance(s, sp.Symbol)]
-        super().__init__("", symbols_read=symbols_read, symbols_defined=self.result_symbols)
-        self._time_step = time_step
-        self._offsets = offsets
-        self.headers = ['"philox_rand.h"']
-        self.keys = tuple(keys)
-        self._args = sp.sympify((dim, time_step, offsets, keys))
-        self._dim = dim
-
-    @property
-    def args(self):
-        return self._args
-
-    @property
-    def undefined_symbols(self):
-        result = {a for a in (self._time_step, *self._offsets, *self.keys) if isinstance(a, sp.Symbol)}
-        loop_counters = [LoopOverCoordinate.get_loop_counter_symbol(i)
-                         for i in range(self._dim)]
-        result.update(loop_counters)
-        return result
-
-    def fast_subs(self, _):
-        return self  # nothing to replace inside this node - would destroy intermediate "dummy" by re-creating them
-
-    def get_code(self, dialect, vector_instruction_set):
-        return _get_philox_code(philox_four_floats_call, dialect, vector_instruction_set,
-                                self._time_step, self._offsets, self.keys, self._dim, self.result_symbols)
-
-    def __repr__(self):
-        return "{}, {}, {}, {} <- PhiloxRNG".format(*self.result_symbols)
-
-
-def random_symbol(assignment_list, seed=TypedSymbol("seed", np.uint32), rng_node=PhiloxTwoDoubles, *args, **kwargs):
-    counter = 0
-    while True:
-        node = rng_node(*args, keys=(counter, seed), **kwargs)
-        inserted = False
-        for symbol in node.result_symbols:
-            if not inserted:
-                assignment_list.insert(0, node)
-                inserted = True
-            yield symbol

pystencils/sympy_gmpy_bug_workaround.py

-# Disable gmpy backend until this bug is resolved if joblib serialize
-# See https://github.com/sympy/sympy/pull/13530
-import os
-import warnings
-
-os.environ['MPMATH_NOGMPY'] = '1'
-try:
-    import mpmath.libmp
-    # In case the user has imported sympy first, then pystencils
-    if mpmath.libmp.BACKEND == 'gmpy':
-        warnings.warn("You are using the gmpy backend. You might encounter an error 'argument is not an mpz sympy'. "
-                      "This is due to a known bug in sympy/gmpy library. "
-                      "To prevent this, import pystencils first then sympy or set the environment variable "
-                      "MPMATH_NOGMPY=1")
-except ImportError:
-    pass
-
-__all__ = []

pystencils_tests/benchmark/main.c

-#include "iacaMarks.h"
-
-int main(int argc, char * argv[]){
-	int a = 0;
-	for(int i = 0; i < argc+100000; i++){
-		IACA_START
-		a += i;
-	}
-	IACA_END
-	return a;
-}

pystencils_tests/kerncraft_inputs/2d-5pt.c

-double a[30][50][3];
-double b[30][50][3];
-double s;
-
-for(int j=1; j<30-1; ++j)
-    for(int i=1; i<50-1; ++i)
-        b[j][i] = ( a[j][i-1] + a[j][i+1]
-                  + a[j-1][i] + a[j+1][i]) * s;

pystencils_tests/kerncraft_inputs/3d-7pt.c

-double a[M][N][N];
-double b[M][N][N];
-double s;
-
-for(int k=1; k<M-1; ++k)
-    for(int j=1; j<N-1; ++j)
-        for(int i=1; i<N-1; ++i)
-            b[k][j][i] = ( a[k][j][i-1] + a[k][j][i+1]
-                         + a[k][j-1][i] + a[k][j+1][i]
-                         + a[k-1][j][i] + a[k+1][j][i]) * s;