In this tutorial we will show how to access and navigate the Iteration/Expression Tree (IET) rooted in an Operator.

Part I - Top Down

Let’s start with a fairly trivial example. First of all, we disable all performance-related optimizations, to maximize the simplicity of the created IET as well as the readability of the generated code.

from devito import configuration
configuration['opt'] = 'noop'
configuration['language'] = 'C'

Then, we create a TimeFunction with 3 points in each of the space Dimensions x and y.

from devito import Grid, TimeFunction

grid = Grid(shape=(3, 3))
u = TimeFunction(name='u', grid=grid)

We now create an Operator that increments by 1 all points in the computational domain.

from devito import Eq, Operator

eq = Eq(u.forward, u+1)
op = Operator(eq)

An Operator is an IET node that can generate, JIT-compile, and run low-level code (e.g., C). Just like all other types of IET nodes, it’s got a number of metadata attached. For example, we can query an Operator to retrieve the input/output Functions.

op.input
(u(t, x, y),)
op.writes
(u(t, x, y),)

If we print op, we can see how the generated code looks like.

print(op)
#define _POSIX_C_SOURCE 200809L
#define START(S) struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);
#define STOP(S,T) gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;

#include "stdlib.h"
#include "math.h"
#include "sys/time.h"

struct dataobj
{
  void *restrict data;
  unsigned long * size;
  unsigned long * npsize;
  unsigned long * dsize;
  int * hsize;
  int * hofs;
  int * oofs;
  void * dmap;
} ;

struct profiler
{
  double section0;
} ;


int Kernel(struct dataobj *restrict u_vec, const int time_M, const int time_m, const int x_M, const int x_m, const int y_M, const int y_m, struct profiler * timers)
{
  float (*restrict u)[u_vec->size[1]][u_vec->size[2]] __attribute__ ((aligned (64))) = (float (*)[u_vec->size[1]][u_vec->size[2]]) u_vec->data;

  for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))
  {
    START(section0)
    for (int x = x_m; x <= x_M; x += 1)
    {
      for (int y = y_m; y <= y_M; y += 1)
      {
        u[t1][x + 1][y + 1] = u[t0][x + 1][y + 1] + 1;
      }
    }
    STOP(section0,timers)
  }

  return 0;
}

An Operator is the root of an IET that typically consists of several nested Iterations and Expressions – two other fundamental IET node types. The user-provided SymPy equations are wrapped within Expressions. Loop nest embedding such expressions are constructed by suitably nesting Iterations.

The Devito compiler constructs the IET from a collection of Clusters, which represent a higher-level intermediate representation (not covered in this tutorial).

The Devito compiler also attaches to the IET key computational properties, such as sequential, parallel, and affine, which are derived through data dependence analysis.

We can print the IET structure of an Operator, as well as the attached computational properties, using the utility function pprint.

from devito.tools import pprint
pprint(op)
<Callable Kernel>
  <CallableBody <unpacks=0, allocs=0, casts=1, maps=0, objs=0> <unmaps=0, frees=0>>
    <List (0, 1, 0)>

      <[affine,sequential] Iteration time::time::(time_m, time_M, 1)>
        <Section (section0)>

          <TimedList (1, 1, 1)>
            <[affine,parallel,parallel=] Iteration x::x::(x_m, x_M, 1)>
              <[affine,parallel,parallel=] Iteration y::y::(y_m, y_M, 1)>
                <ExpressionBundle (1)>

                  <Expression u[t1, x + 1, y + 1] = u[t0, x + 1, y + 1] + 1>

In this example, op is represented as a <Callable Kernel>. Attached to it are metadata, such as _headers and _includes, as well as the body, which includes the children IET nodes. Here, the body is the concatenation of an PointerCast and a List object.

op._headers
OrderedSet([('_POSIX_C_SOURCE', '200809L'), ('START(S)', 'struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);'), ('STOP(S,T)', 'gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;')])
op._includes
OrderedSet(['stdlib.h', 'math.h', 'sys/time.h'])
op.body
<CallableBody <unpacks=0, allocs=0, casts=1, maps=0, objs=0> <unmaps=0, frees=0>>

We can explicitly traverse the body until we locate the user-provided SymPy equations.

print(op.body.casts[0])  # Printing the PointerCast
float (*restrict u)[u_vec->size[1]][u_vec->size[2]] __attribute__ ((aligned (64))) = (float (*)[u_vec->size[1]][u_vec->size[2]]) u_vec->data;
print(op.body.body[0])  # Printing the actual body
for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))
{
  START(section0)
  for (int x = x_m; x <= x_M; x += 1)
  {
    for (int y = y_m; y <= y_M; y += 1)
    {
      u[t1][x + 1][y + 1] = u[t0][x + 1][y + 1] + 1;
    }
  }
  STOP(section0,timers)
}

Below we access the Iteration representing the time loop.

t_iter = op.body.body[0].body[0]
t_iter
<WithProperties[affine,sequential]::Iteration time[t0,t1]; (time_m, time_M, 1)>

We can for example inspect the Iteration to discover what its iteration bounds are.

t_iter.limits
(time_m, time_M, 1)

And as we keep going down through the IET, we can eventually reach the Expression wrapping the user-provided SymPy equation.

expr = t_iter.nodes[0].body[0].body[0].nodes[0].nodes[0].body[0]
expr.view
'<Expression u[t1, x + 1, y + 1] = u[t0, x + 1, y + 1] + 1>'

Of course, there are mechanisms in place to, for example, find all Expressions in a given IET. The Devito compiler has a number of IET visitors, among which FindNodes, usable to retrieve all nodes of a particular type. So we easily can get all Expressions within op as follows

from devito.ir.iet import Expression, FindNodes
exprs = FindNodes(Expression).visit(op)
exprs[0].view
'<Expression u[t1, x + 1, y + 1] = u[t0, x + 1, y + 1] + 1>'
Back to top