# Notes on stencils

Numba provides the @stencil decorator to represent stencil computations. This document explains how this feature is implemented in the several different modes available in Numba. Currently, calls to the stencil from non-jitted code is supported as well as calls from jitted code, either with or without the parallel=True option.

## The stencil decorator

The stencil decorator itself just returns a `StencilFunc` object. This object encapsulates the original stencil kernel function as specified in the program and the options passed to the stencil decorator. Also of note is that after the first compilation of the stencil, the computed neighborhood of the stencil is stored in the `StencilFunc` object in the `neighborhood` attribute.

## Handling the three modes

As mentioned above, Numba supports the calling of stencils from inside or outside a `@jit` compiled function, with or without the parallel=True option.

### Outside jit context

`StencilFunc` overrides the `__call__` method so that calls to `StencilFunc` objects execute the stencil:

    def __call__(self, *args, **kwargs):
        result = kwargs.get('out')

        new_stencil_func = self._stencil_wrapper(result, None, *args)

        if result is None:
            return new_stencil_func.entry_point(*args)
        else:
            return new_stencil_func.entry_point(*args, result)

First, the presence of the optional out parameter is checked. If it is present then the output array is stored in `result`. Then, the call to `_stencil_wrapper` generates the stencil function given the result and argument types and finally the generated stencil function is executed and its result returned.

### Jit without `parallel=True`

When constructed, a `StencilFunc` inserts itself into the typing context’s set of user functions and provides the `_type_me` callback. In this way, the standard Numba compiler is able to determine the output type and signature of a `StencilFunc`. Each `StencilFunc` maintains a cache of previously seen combinations of input argument types and keyword types. If previously seen, the `StencilFunc` returns the computed signature. If not previously computed, the `StencilFunc` computes the return type of the stencil by running the Numba compiler frontend on the stencil kernel and then performing type inference on the Numba IR (IR) to get the scalar return type of the kernel. From that, a Numpy array type is constructed whose element type matches that scalar return type.

After computing the signature of the stencil for a previously unseen combination of input and keyword types, the `StencilFunc` then creates the stencil function itself. `StencilFunc` then installs the new stencil function’s definition in the target context so that jitted code is able to call it.

Thus, in this mode, the generated stencil function is a stand-alone function called like a normal function from within jitted code.

### Jit with `parallel=True`

When calling a `StencilFunc` from a jitted context with `parallel=True`, a separate stencil function as generated by Creating the stencil function is not used. Instead, parfors (Stage 5b: Perform Automatic Parallelization) are created within the current function that implements the stencil. This code again starts with the stencil kernel and does a similar kernel size computation but then rather than standard Python looping syntax, corresponding parfors are created so that the execution of the stencil will take place in parallel.

The stencil to parfor translations can also be selectively disabled by setting `parallel={'stencil': False}`, among other sub-options described in Stage 5b: Perform Automatic Parallelization.

## Creating the stencil function

Conceptually, a stencil function is created from the user-specified stencil kernel by adding looping code around the kernel, transforming the relative kernel indices into absolute array indices based on the loop indices, and replacing the kernel’s `return` statement with a statement to assign the computed value into the output array.

To accomplish this transformation, first, a copy of the stencil kernel IR is created so that subsequent modifications of the IR for different stencil signatures will not effect each other.

Then, an approach similar to how GUFunc’s are created for parfors is employed. In a text buffer, a Python function is created with a unique name. The input array parameter is added to the function definition and if the `out` argument type is present then an `out` parameter is added to the stencil function definition. If the `out` argument is not present then first an output array is created with `numpy.zeros` having the same shape as the input array.

The kernel is then analyzed to compute the stencil size and the shape of the boundary (or the `neighborhood` stencil decorator argument is used for this purpose if present). Then, one `for` loop for each dimension of the input array is added to the stencil function definition. The range of each loop is controlled by the stencil kernel size previously computed so that the boundary of the output image is not modified but instead left as is. The body of the innermost `for` loop is a single `sentinel` statement that is easily recognized in the IR. A call to `exec` with the text buffer is used to force the stencil function into existence and an `eval` is used to get access to the corresponding function on which `run_frontend` is used to get the stencil function IR.

Various renaming and relabeling is performed on the stencil function IR and the kernel IR so that the two can be combined without conflict. The relative indices in the kernel IR (i.e., `getitem` calls) are replaced with expressions where the corresponding loop index variables are added to the relative indices. The `return` statement in the kernel IR is replaced with a `setitem` for the corresponding element in the output array. The stencil function IR is then scanned for the sentinel and the sentinel replaced with the modified kernel IR.

Next, `compile_ir` is used to compile the combined stencil function IR. The resulting compile result is cached in the `StencilFunc` so that other calls to the same stencil do not need to undertake this process again.

## Exceptions raised

Various checks are performed during stencil compilation to make sure that user-specified options do not conflict with each other or with other runtime parameters. For example, if the user has manually specified a `neighborhood` to the stencil decorator, the length of that neighborhood must match the dimensionality of the input array. If this is not the case, a `ValueError` is raised.

If the neighborhood has not been specified then it must be inferred and a requirement to infer the kernel is that all indices are constant integers. If they are not, a `ValueError` is raised indicating that kernel indices may not be non-constant.

Finally, the stencil implementation detects the output array type by running Numba type inference on the stencil kernel. If the return type of this kernel does not match the type of the value passed to the `cval` stencil decorator option then a `ValueError` is raised.