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.



In [1]:

    
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.



In [2]:

    
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.



In [3]:

    
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.



In [4]:

    
op.input









    Out[4]:





(u(t, x, y),)



In [5]:

    
op.output









    Out[5]:





(u(t, x, y),)

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



In [6]:

    
print(op)









    



#define _POSIX_C_SOURCE 200809L
#include "stdlib.h"
#include "math.h"
#include "sys/time.h"

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

struct profiler
{
  double section0;
} ;


int Kernel(struct dataobj *restrict u_vec, const int time_M, const int time_m, struct profiler * timers, const int x_M, const int x_m, const int y_M, const int y_m)
{
  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))
  {
    struct timeval start_section0, end_section0;
    gettimeofday(&start_section0, NULL);
    /* Begin 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;
      }
    }
    /* End section0 */
    gettimeofday(&end_section0, NULL);
    timers->section0 += (double)(end_section0.tv_sec-start_section0.tv_sec)+(double)(end_section0.tv_usec-start_section0.tv_usec)/1000000;
  }
  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.



In [7]:

    
from devito.tools import pprint
pprint(op)









    



<Callable Kernel>
  <ArrayCast>
  <List (0, 1, 0)>

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

          <[affine,parallel,parallel=,tilable] Iteration x::x::(x_m, x_M, 1)>
            <[affine,parallel,parallel=,tilable] 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 ArrayCast and a List object.



In [8]:

    
op._headers









    Out[8]:





['#define _POSIX_C_SOURCE 200809L']



In [9]:

    
op._includes









    Out[9]:





['stdlib.h', 'math.h', 'sys/time.h']



In [10]:

    
op.body









    Out[10]:





(<ArrayCast(u(t, x, y))>, <List (0, 1, 0)>)

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



In [11]:

    
print(op.body[0])  # Printing the ArrayCast









    



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;



In [12]:

    
print(op.body[1])  # Printing the List









    



for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))
{
  struct timeval start_section0, end_section0;
  gettimeofday(&start_section0, NULL);
  /* Begin 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;
    }
  }
  /* End section0 */
  gettimeofday(&end_section0, NULL);
  timers->section0 += (double)(end_section0.tv_sec-start_section0.tv_sec)+(double)(end_section0.tv_usec-start_section0.tv_usec)/1000000;
}

Below we access the Iteration representing the time loop.



In [13]:

    
t_iter = op.body[1].body[0]
t_iter









    Out[13]:





<WithProperties[affine,sequential,wrappable]::Iteration time[t0,t1]; (time_m, time_M, 1)>

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



In [14]:

    
t_iter.limits









    Out[14]:





(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.



In [15]:

    
expr = t_iter.nodes[0].body[0].body[0].nodes[0].nodes[0].body[0]
expr.view









    Out[15]:





'<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



In [16]:

    
from devito.ir.iet import Expression, FindNodes
exprs = FindNodes(Expression).visit(op)
exprs[0].view









    Out[16]:





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