Integer Set Library (ISL) - A Primer

Polyhedral optimization is a tool used in compilers for optimizing loop nests. While the major compilers that use this implement polyhedral optimizations from scratch,¹ there is a generally-applicable open source C library called the Integer Set Library (ISL) that implements the core algorithms used in polyhedral optimization.

This article gives an overview of a subset of ISL, mainly focusing on the representation of sets and relations and basic manipulations on them. Because my personal use of ISL is related to MLIR, I will also include some notes about the MLIR port of ISL (called the Fast Presburger Library, or FPL), as well as an interoperability layer I wrote between ISL and FPL. This article does not cover the details of the algorithms used in ISL.

All the code written for this article is available on GitHub.

Integer sets and quasi-affine formulas

The core data object in ISL is a set of integer points, and the core of polyhedral optimization is to analyze and manipulate such sets. Because the integer sets of interest can be large or infinite, the points are represented indirectly as the set of integer points satisfying a system of equalities and inequalities with a restricted form called quasi-affine formulas (defined formally after the following motivation).

One guiding principle in polyhedral optimization is to construct an integer set corresponding to the “iteration space” of a loop nest. That is, there is one point in the integer set for each execution of a statement within a loop. I’ll work with a simplified loop nest that is a “perfect nest” that contains a single inner-most statement abstracted as S.

for (int i = 0; i < 100; i++) {
    for (int j = 0; j < 200; j++) {
        S(i, j);
    }
}

We want to construct an integer set that contains all pairs (i, j) for which S(i, j) is executed. More generally, a loop nest can have different offsets, bounds, and step sizes, as well as an if statement guarding the execution of S.

for (int i = 3; i < 100; i += 5) {
    for (int j = 100; j > 0; j -= 2) {
        if ((i + j) % 2 == 0 {
            S(i, j);
        }
    }
}

The observation of polyhedral optimization is that most for loops and if statements that occur in practical code have a limited structure. The loops have constant bounds and static, constant step sizes (e.g., we don’t typically increment j as j *= 2 or j += i*i) and the conditions in if statements typically involve just a little bit more than linear arithmetic. Notably, code often uses floor/ceiling division and modulo operations to express things like operating over a tiled memory layout.

This was the inspiration for the restriction of an integer set to be defined by a system of quasi-affine formulas.

Definition: A quasi-affine formula is a multivariate formula built from the following operations:

• Integer literals
• Integer-valued variables
• addition and subtraction
• multiplication by an integer constant
• floor- and ceiling-rounded division by a nonzero integer constant
• modulus by a nonzero integer constant

Borrowing a BNF grammar from the MLIR website, we can also define it as

affine-expr ::= `(` affine-expr `)`
              | affine-expr `+` affine-expr
              | affine-expr `-` affine-expr
              | `-`? integer-literal `*` affine-expr
              | affine-expr `ceildiv` integer-literal
              | affine-expr `floordiv` integer-literal
              | affine-expr `mod` integer-literal
              | `-`affine-expr
              | bare-id
              | `-`? integer-literal

A Presburger formula is a first-order logical formula involving quasi-affine expressions and the standard binary comparison operators (=, !=, <, <=, >, >=). In particular, Presburger formulas may contain existentially and universally quantified variables.

A Presburger set (or just integer set, scoped to this article) is a subset of points in $\mathbb{Z}^n$ satisfying a Presburger formula.

Examples using isl_cat

Our first foray into using ISL will be to show some examples of Presburger sets.

The code repository for this article uses bazel with a pre-configured ISL dependency, so if you have bazel (or better bazelisk) installed, you can clone the repo and build ISL with

git clone git@github.com:j2kun/isl-primer.git
cd isl-primer
bazel build @isl//:isl

The file examples.isl contains the example Presburger sets we will show, and the Polly Labs Playground gives a UI that can be used to visualize the sets.² Otherwise, you can run them through bazel build @isl//:isl_cat to have ISL parse them, validate/simplify them, and print them back out.

bazel run @isl//:isl_cat < examples.isl

# isl_cat reads from stdin, so run without args and paste a formula to the terminal
bazel run @isl//:isl_cat
{ [i, j, a, b] : i = a + 1 and (b + j) mod 2 = 0 }

ISL splits the kinds of sets you can define into a few categories. The simplest is a basic set, which is defined as a Presburger set that happens to be a single convex integer polyhedron. As an example, this defines a triangle:

{ [i,j] : 0 < i < 10 and 0 < j < 10 and i < j }

The syntax here is close to a set-builder notation from math. The curly braces {} denote a set, the square brackets [ ] around i, j define the two variables that name the two dimensions of the underlying integer lattice $\mathbb{Z}^2$, and the colon : separates the variable definition from the formulas.

The “underlying integer lattice” is called a space in ISL terminology. A space is simply defined by a number of dimensions, with names for each dimension, and possibly a name for the overall space itself. To name the space above, we put a symbol before the opening square bracket.

{ A[i,j] : 0 < i < 10 and 0 < j < 10 and i < j }

The names of the spaces can be used to identify statements in a loop nest.

Next, ISL defines a set (as opposed to a basic set) as a Presburger set that happens to have multiple disjoint convex components. These have to be associated with the same underlying space. The Presburger set below defines a union of a triangle and a square.

{ [i, j] : 0 <= i <= 10 and 0 <= j <= 10 and j > i or (6 <= i <= 9 and 0 < j <= 4) }

Next, a union set is a finite union of Presburger sets, which may use different underlying spaces with different dimensions. A semicolon separates the space definitions.

{
    A[i, j] : 0 <= i <= 10 and 0 <= j <= 10;
    B[i, j, k] : 6 <= i <= 9 and 0 < j <= 4 and 0 <= k < 2
}

As far as I understand, the purpose of a union set is to represent two code statements (A and B) within the same structure, as they may be statements that share a common parent loop and should be optimized jointly. However, I have not worked with this in detail myself.

Finally there are map analogues of the above set types. A map type is just another Presburger set, but for which some of the variables are distinguished as domain variables, while others are codomain variables (in ISL they are called range variables). An -> separates the domain from the codomain which are otherwise defined as spaces as above.

For example, a basic map splits the variables of a basic set, e.g.,

{
    [i, j] -> [a, b] : i = a + 1 and (b + j) mod 2 = 0
    and 0 <= i < 10 and 0 <= j < 10
    and 0 <= a < 10 and 0 <= b < 10
}

A map is like a basic map but may have multiple disjoint components, and a union map is a finite union of maps.

Mathematically speaking, maps in polyhedral analysis are not really “maps,” but rather binary relations. Nothing here stops two domain points from “mapping” to different range points, or vice versa. Using relations is a strength of the formalism, because it adds more generality. But because, I suppose, someone thought typing “relation” instead of “map” was too laborious, we now have to live with this extra little cognitive overhead.

Finally, polyhedral analysis (and ISL) also supports symbolic constants by prefixing a Presburger set with a list of symbols in square brackets and an ->.

For example, the Presburger set below uses a symbolic constant N to define the upper bound of the variables i and j.

[N] -> { [i, j] : 0 < i < N and j > i and 0 < j < N }

Another supported syntax is to represent an existential local variable. Below is an example of a presburger set where a local variable ensures the domain and range form sequential pairs here the first entry in the pair is even.

{
    [i] -> [j] :
    exists q :
    0 <= i < N
    and 0 <= j < N
    and i = 2*q
    and j = 2*q + 1
}

Generating C code

The first nontrivial thing we can do is to generate C code that iterates over the points in a Presburger set and executes a statement for each point in the set.

To do this we need to specify what’s called an iteration schedule. Iteration schedules can be quite complex, and to my knowledge there is no systematic way to determine what is the “right” iteration schedule. It’s application dependent. So I will only scratch the surface here and focus on the simplest case where you specify what subset of variables are iterated over, and in what order.

This specification is defined by converting a presburger set to a map, which maps all the space variables to the subset of variables you want to iterate over.

For the example of the N-parameterized triangle above,

[N] -> { [i, j] : 0 < i < N and j > i and 0 < j < N }

we can convert the [i,j] space to A[i, j] -> [i, j] to tell it to iterate over both i and j, in that order.

[N] -> { A[i, j] -> [i, j] : 0 < i < N and j > i and 0 < j < N }

And then run isl_codegen on the result. (Ingore the second two lines of the input below, they are for advanced control of the generated code)

bazel run @isl//:isl_codegen << EOF
[N] -> { A[i, j] -> [i, j] : 0 < i < N and j > i and 0 < j < N }
{ : }
{}
EOF

for (int c0 = 1; c0 < N - 1; c0 += 1)
  for (int c1 = c0 + 1; c1 < N; c1 += 1)
    A(c0, c1);

Or reversing the order of iteration by mapping to [j, i] instead of [i, j]:

bazel run @isl//:isl_codegen << EOF
[N] -> { A[i, j] -> [j, i] : 0 < i < N and j > i and 0 < j < N }
{ : }
{}
EOF

for (int c0 = 2; c0 < N; c0 += 1)
  for (int c1 = 1; c1 < c0; c1 += 1)
    A(c1, c0);

Note how, because the constraints include $0 < i < j$, the generated code knows to start j and 2 and i at 1, as those are the smallest values satisfying the constraints. Anote note how it swapped the order the induction variables are passed to the statement A when the iteration order was switched.

One more example: the produced loop nest can include if statements, or arithmetic in the final parameters passed to the statement. To show this I use an integer relation that describes a diagonal packing from an earlier article on this blog about diagonal data packing.

bazel run @isl//:isl_codegen << EOF
{ [i0, i1] -> [ct, slot] : (i0 - i1 + ct) mod 16 = 0 and (-i0 + slot) mod 16 = 0 and 0 <= i0 <= 9 and 0 <= i1 <= 15 and 0 <= ct <= 15 and 0 <= slot <= 1023 }
{ : }
{}
EOF

for (int c0 = 0; c0 <= 15; c0 += 1)
  for (int c1 = 0; c1 <= 1023; c1 += 1)
    if ((c1 + 6) % 16 >= 6)
      (c1 % 16, (c0 + c1) % 16);

Later in this article I’ll show code that can be used to generate this using the ISL C API.

Basic manipulation with the ISL C API

Hello, ISL: extracting the domain of a map

The ISL C API is quite a beast. By my count there are over 20 thousand different function declarations, and because it’s a pure-C library, they look like this

__isl_give isl_set *isl_map_domain(__isl_take isl_map *bmap);

There are types for the isl_basic_set, isl_set, isl_union_set, isl_basic_map, isl_map, and isl_union_map.³ And the functions are named so that they start with isl_, then the type of the main argument being operated on (e.g., map), then a description of the operation.

So the isl_map_domain function takes as input an isl_map, which is frees before returning a pointer to an isl_set representing the subset of integer points that are in the domain of the map.

ISL internally does some reference counting, each object is associated with an isl_ctx context object that handles this. The macros __isl_give, __isl_take, and __isl_keep (which are no-ops as far as I can tell), are used to indicate to the caller whether the function takes ownership of the argument pointer (__isl_take), whether the caller keeps ownership (__isl_keep), and whether the caller is responsible for freeing the returned pointer (__isl_give). Even with this helpful convention, I still find myself leaking memory when using ISL, so running bazel with --config=asan is a must.

So as a hello world example, here is a C++ program that creates a context, parses a basic map from a string, and then extracts the domain and prints it.

#include "include/isl/ctx.h"
#include "include/isl/map.h"
#include "include/isl/map_type.h"
#include "include/isl/set.h"

std::string parse_map_and_extract_domain_as_string(
    isl_ctx *ctx, std::string islStr) {
  isl_set *domainMap = isl_map_domain(
        isl_map_read_from_str(ctx, islStr.c_str()));
  char *result_str = isl_set_to_str(domainMap);
  std::string result(result_str);

  free(result_str);
  isl_set_free(domainMap);
  return result;
}

And a simple binary to run it:

#include "isl_api_examples.h"
#include <iostream>
#include <string>
#include "include/isl/ctx.h"

int main() {
  std::string islMap;
  std::cout << "Enter an ISL map: ";
  std::getline(std::cin, islMap);

  isl_ctx *ctx = isl_ctx_alloc();
  std::string result = parse_map_and_extract_domain_as_string(ctx, islMap);
  std::cout << "Result: " << result << std::endl;
  isl_ctx_free(ctx);
  return 0;
}

The build file

load("@rules_cc//cc:cc_binary.bzl", "cc_binary")
load("@rules_cc//cc:cc_library.bzl", "cc_library")
load("@rules_cc//cc:cc_test.bzl", "cc_test")

cc_library(
    name = "isl_api_examples",
    srcs = ["isl_api_examples.cpp"],
    hdrs = ["isl_api_examples.h"],
    deps = ["@isl"],
)

cc_binary(
    name = "get_domain",
    srcs = ["get_domain.cpp"],
    deps = [":isl_api_examples"],
)

and running it:

$ bazel run :get_domain
Enter an ISL map: { [i] -> [j] : (-i + j) mod 7 = 0 and 0 <= i <= 20 and 0 <= j <= 20 }
Result: { [i] : 0 <= i <= 20 }

# In this next example, I replace the upper bound on j with 3, to demonstrate
# that the extracted domain really only includes the subset of i values for which
# there is a corresponding j value.
$ bazel run :get_domain
Enter an ISL map: { [i] -> [j] : (-i + j) mod 7 = 0 and 0 <= i <= 20 and 0 <= j <= 3 }
Result: { [i] : 0 <= i <= 20 and 7*floor((-1 - i)/7) <= -4 - i }

Composing maps

While the ISL API has a number of basic operations like taking unions and intersections of sets and maps, in my opinion things get interesting when you start composing maps. Mainly this is because, a general composition of two relations requires introducing existentially quantified variables, and ISL has powerful algorithms for simplifying and eliminating these to arrive at a nice representation of the composed map.

Because most of my present uses for polyhedral optimization involve analyzing memory layouts, this example will show how, if a tensor is laid out in memory according to a given integer relation, then a modification to the tensor in the program can be composed with the original layout to get a descriptor of the new layout.

For example, suppose we have a 3-dimensional tensor, and a Presburger relation that describes a column-major layout in vector registers of size $16$.

{
    S[i, j, k] -> [reg, lane] :
        0 <= i < 32
    and 0 <= j < 64
    and 0 <= k < 1024
    and 0 <= lane < 16
    and i + j * 32 + k * 32 * 64 = 16 * reg + lane
}

Now suppose the program includes an operation that implements a transposition by swapping axes i and j. Subsequent analysis after this step would need to incorporate this operation into a new layout. This can be done by precomposing the memory layout relation with the following relation:

{ S[i, j, k] -> [j, i, k] }

Which is shorthand for

{ S[i0, j0, k0] -> [i1, j1, k1] : i0 = j1 and i1 = j0 and k0 = k1 }

To do this in ISL, we will peek behind the curtain a bit to see how the underlying constraints are organized as a matrix of coefficients of equalities and inequalities.

The declaration: parse a layout from a string, then precompose it with a relation that swaps two specified domain dimensions.

std::string precompose_transposition(isl_ctx* ctx, std::string starting_layout,
                                     int dim1, int dim2);

The implementation as a whole, then line by line.

std::string precompose_transposition(isl_ctx* ctx, std::string starting_layout,
                                     int dim1, int dim2) {
  isl_map* layout = isl_map_read_from_str(ctx, starting_layout.c_str());
  isl_set* domain = isl_map_domain(isl_map_copy(layout));

  isl_map* domain_map =
      isl_map_from_domain_and_range(isl_set_copy(domain), domain);
  isl_space* transpose_space = isl_map_get_space(domain_map);
  isl_map_free(domain_map);
  unsigned num_domain_vars = isl_space_dim(transpose_space, isl_dim_in);

  isl_mat* eq_mat =
      create_empty_constraint_matrix(transpose_space, num_domain_vars);
  isl_mat* ineq_mat = create_empty_constraint_matrix(transpose_space, 0);

  // Column order is: [domain_vars, range_vars, div_vars, symbol_vars, constant]
  // so the offset between a domain variable and its corresponding range
  // variable is num_domain_vars
  for (int domain_var = 0; domain_var < num_domain_vars; ++domain_var) {
    // First constraint: domain_var dim1 - range_var dim2 = 0
    if (domain_var == dim1) {
      isl_mat_set_element_si(eq_mat, /*row=*/dim1, /*col=*/dim1, 1);
      isl_mat_set_element_si(eq_mat, /*row=*/dim1,
                             /*col=*/num_domain_vars + dim2, -1);
      continue;
    }

    if (domain_var == dim2) {
      // Second constraint: domain_var dim2 - range_var dim1 = 0
      isl_mat_set_element_si(eq_mat, /*row=*/dim2, /*col=*/dim2, 1);
      isl_mat_set_element_si(eq_mat, /*row=*/dim2,
                             /*col=*/num_domain_vars + dim1, -1);
      continue;
    }

    // Otherwise, domain_var d - range_var d = 0
    isl_mat_set_element_si(eq_mat, /*row=*/domain_var, /*col=*/domain_var, 1);
    isl_mat_set_element_si(eq_mat, /*row=*/domain_var,
                           /*col=*/num_domain_vars + domain_var, -1);
  }

  isl_map* transpose =
      isl_map_from_basic_map(isl_basic_map_from_constraint_matrices(
          transpose_space, eq_mat, ineq_mat, isl_dim_in, isl_dim_out,
          isl_dim_div, isl_dim_param, isl_dim_cst));

  isl_map* composed = isl_map_apply_range(transpose, layout);

  char* result_str = isl_map_to_str(composed);
  std::string result(result_str);

  free(result_str);
  isl_map_free(composed);
  return result;
}

After parsing the layout, the first step is to create an underlying space for the new Presburger map which uses the exact same domain variables as the input layout. Most of this step is jumping through ISL API hoops and managing memory.

isl_set* domain = isl_map_domain(isl_map_copy(layout));
isl_map* domain_map =
  isl_map_from_domain_and_range(isl_set_copy(domain), domain);
isl_space* transpose_space = isl_map_get_space(domain_map);
isl_map_free(domain_map);

Next we create new empty constraint matrices to store the constraints that will represent the transposition.

isl_mat* eq_mat = create_empty_constraint_matrix(transpose_space, num_domain_vars);
isl_mat* ineq_mat = create_empty_constraint_matrix(transpose_space, 0);

This uses a helper function that allocates a matrix and zero-initializes all entries.

__isl_give isl_mat* create_empty_constraint_matrix(__isl_keep isl_space* space,
                                                   unsigned num_constraints) {
  unsigned num_in = isl_space_dim(space, isl_dim_in);
  unsigned num_out = isl_space_dim(space, isl_dim_out);
  unsigned num_div = isl_space_dim(space, isl_dim_div);
  unsigned num_param = isl_space_dim(space, isl_dim_param);

  // The layout of the columns is: [domain_vars, range_vars, div_vars,
  // symbol_vars, constant]
  unsigned num_cols = num_in + num_out + num_div + num_param + 1;

  isl_mat* mat =
      isl_mat_alloc(isl_space_get_ctx(space), num_constraints, num_cols);

  for (int i = 0; i < isl_mat_rows(mat); ++i)
    for (int j = 0; j < isl_mat_cols(mat); ++j)
      isl_mat_set_element_si(mat, /*row=*/i, /*col=*/j, 0);
  return mat;
}

Next we define the actual constraints by setting individual coefficients. Recall, our desired relation is:

{ S[i0, j0, k0] -> [i1, j1, k1] : i0 = j1 and i1 = j0 and k0 = k1 }

And so we want one constraint per (domain, range) variable pair. Iterating over all the domain variables, we have two constraints special for the transposed dimensions:

// First constraint: domain_var dim1 - range_var dim2 = 0
isl_mat_set_element_si(eq_mat, /*row=*/0, /*col=*/dim1, 1);
isl_mat_set_element_si(eq_mat, /*row=*/0, /*col=*/num_domain_vars + dim2, -1);

// Second constraint: domain_var dim2 - range_var dim1 = 0
isl_mat_set_element_si(eq_mat, /*row=*/1, /*col=*/dim2, 1);
isl_mat_set_element_si(eq_mat, /*row=*/1, /*col=*/num_domain_vars + dim1, -1);

And one that equates the domain and range variables by default.

// Otherwise, domain_var d - range_var d = 0
isl_mat_set_element_si(eq_mat, /*row=*/domain_var, /*col=*/domain_var, 1);
isl_mat_set_element_si(eq_mat, /*row=*/domain_var,
                       /*col=*/num_domain_vars + domain_var, -1);

The constraint matrix specifies one constraint per row, one variable per column, and each row is a linear combination of the variables and coefficients, with the last column representing the constant term.

Finally, we construct the relation from the constraint matrices and call apply_range on the result. The construction allows us to specify the intended order of the columns in relation to the variables (e.g., the columns have domain variables first for isl_dim_in, then range variables (isl_dim_out), then other variables that are not relevant here.

isl_map* transpose =
  isl_map_from_basic_map(isl_basic_map_from_constraint_matrices(
      transpose_space, eq_mat, ineq_mat, isl_dim_in, isl_dim_out,
      isl_dim_div, isl_dim_param, isl_dim_cst));

isl_map* composed = isl_map_apply_range(transpose, layout);

The apply_range call takes the first argument, and composes its range with the domain of the second argument to construct a composite relation. The rest of the code converts the result to a string and frees memory as before.

While the code is mucky, it exposes one critical question to understand how ISL works internally: how does this representation of constraints as linear combinations of variables permit us to represent modular arithmetic, and floor and ceiling division? The answer is in the div dimension. In ISL, a div dimension is a special variable, always existentially quantified, that is used to represent the intermediate computation of a division.

For example, if you print the constraint matrix for a simple set like { [i] : 5 <= i } (I included a binary and API example to dump the raw constraint matrix for a basic_set), you get

bazel run dump_constraints << EOF
{ [i] : i = 5 }
EOF

Equality Matrix:

[[1,-5]]

The sole equality here is signifying $1 \cdot i - 5 = 0$, as the first column is the coefficient of i and the last column is the constant term.

However, if you dump a more complicated system that has modular arithmetic, like { [i] : (i % 2) = 1 }, you will get the following equality constraint:

bazel run dump_constraints << EOF
{ [i] : (i % 2) = 1 }
EOF

Equality Matrix:

[[-1,2,-1]]

The last column is still the constant term, and the second-to-last column is a new “div” variable $d$ representing the result of an integer division. $d$ is existentially quantified, so the formula is really

Rearranging and taking both sides of the equation mod 2 recovers the original formula. ISL does a similar trick for floor and ceiling division; all quasi-affine expressions can reduce such operations to a system of purely linear equalities and inequalities, where the intermediate divisors are represented as existentially quantified variables.

Manually adding your own div/mod variables to the constraint system is a bit of a pain. ISL has a helper API for constructing quasi-affine expressions at include/isl/aff.h, but for brevity I will defer the exploration of that API to a future article. ⁴

Enumerating maps

Next I’ll show another debugging tool: enumerating and dumping all the points in a basic map or set.

The key API is isl_set_foreach_point, which takes a callback function pointer and calls it for each point in the set. In my example, the callback inserts points into a struct that is aware of the prior domain/range breakdown and is able to split a point into domain and range points.

struct PointPairCollector {
  using Point = std::vector<int64_t>;
  std::vector<std::pair<Point, Point>> points;
  int domain_dims;
  int range_dims;

  PointPairCollector(int domain_dims, int range_dims)
      : domain_dims(domain_dims), range_dims(range_dims) {}
};

Then the function pointer to populate it:

isl_stat pointPairCallback(__isl_take isl_point* pnt, void* user) {
  PointPairCollector* collector = static_cast<PointPairCollector*>(user);

  isl_space* space = isl_point_get_space(pnt);
  isl_space_free(space);

  std::vector<int64_t> domainPoint(collector->domain_dims);
  std::vector<int64_t> rangePoint(collector->range_dims);

  // Extract domain coordinates
  for (int i = 0; i < collector->domain_dims; i++) {
    isl_val* coord = isl_point_get_coordinate_val(pnt, isl_dim_set, i);
    if (isl_val_is_int(coord)) {
      domainPoint[i] = isl_val_get_num_si(coord);
    }
    isl_val_free(coord);
  }

  // Extract range coordinates
  for (int i = 0; i < collector->range_dims; i++) {
    isl_val* coord = isl_point_get_coordinate_val(pnt, isl_dim_set,
                                                  collector->domain_dims + i);
    if (isl_val_is_int(coord)) {
      rangePoint[i] = isl_val_get_num_si(coord);
    }
    isl_val_free(coord);
  }

  collector->points.emplace_back(domainPoint, rangePoint);
  isl_point_free(pnt);
  return isl_stat_ok;
}

And the routine to combine them:

void enumeratePoints(isl_basic_map *bmap,
                     PointPairCollector& collector) {
  isl_set* set = isl_set_from_basic_set(isl_basic_map_wrap(bmap));
  isl_set_foreach_point(set, &pointPairCallback, &collector);
  isl_set_free(set);
}

AST generation

Next I’ll show how to use the C API to generate the loop nest as a C string.⁵

std::string generate_loop_nest_as_c_str(isl_basic_map *bmap) {
  auto *ctx = isl_basic_map_get_ctx(bmap);
  isl_union_map *schedule =
      isl_union_map_from_basic_map(isl_basic_map_copy(bmap));

  // Context and options are intentionally empty.
  isl_set *context = isl_set_universe(isl_space_params_alloc(ctx, 0));
  isl_union_map *options = isl_union_map_empty(isl_space_params_alloc(ctx, 0));

  // Build the AST
  isl_ast_build *build = isl_ast_build_from_context(context);
  build = isl_ast_build_set_options(build, options);
  isl_ast_node *tree = isl_ast_build_node_from_schedule_map(build, schedule);
  isl_ast_build_free(build);

  char *cStr = isl_ast_node_to_C_str(tree);
  std::string actual = std::string(cStr);
  free(cStr);
  isl_ast_node_free(tree);

  // Add a leading newline for ease of comparison with multiline strings.
  return actual.insert(0, "\n");
}

The key API is isl_ast_build_node_from_schedule_map, which takes a map whose [domain] -> [range] part specifies the iteration order as mentioned above. Then isl_ast_build_node_from_schedule_map cosntructs an isl_ast_node, which is a classical AST structure with nodes for loops, if statements, etc. The isl_ast_node_to_C_str is a helper that converts the AST to a C string, but you can walk the AST yourself if you want to, say, generate an MLIR loop nest from the AST. You can see how we do this in the HEIR project here.

Interoperating with MLIR’s FPL

The MLIR compiler project has its own internal implementation of a subset of the core algorithms of ISL, called the Fast Polyhedral Library (FPL). However, because it’s incomplete, there are some features of ISL, including code generation and advanced simplification algorithms, that our group at the HEIR project have found critical.

So until FPL is a more complete port of ISL, we added an interoperability layer between MLIR and ISL.⁶ It can convert between MLIR’s IntegerRelation and ISL’s basic_map, and it does so simply by copying the underlying constraint matrices between the two formats. This is part of why I avoided the higher-level ISL APIs in my examples above, so that this conversion layer would be easier to understand.

Adding MLIR as a dependency to the code repo for this project would incur massive build times for little benefit, so instead I’ll point to the code and show the core snippets here.

Because both libraries use the same representation internally, a set of constraint matrices for equalities and inequalities, the conversion requires simply copying the matrices back and forth.

Converting a basic_map to an IntegerRelation:

presburger::IntegerRelation convertBasicMapToRelation(
    isl_basic_map* bmap) {
  // Variables in an IntegerRelation are stored in the order
  //
  //   Domain, Range, Symbols, Locals, Constant
  //
  // In ISL this corresponds to
  //   In,     Out,   Param,   Div,    Cst
  isl_mat* eqMat = isl_basic_map_equalities_matrix(
      bmap, isl_dim_in, isl_dim_out, isl_dim_param, isl_dim_div, isl_dim_cst);
  isl_mat* ineqMat = isl_basic_map_inequalities_matrix(
      bmap, isl_dim_in, isl_dim_out, isl_dim_param, isl_dim_div, isl_dim_cst);

  PresburgerSpace fplSpace = PresburgerSpace::getRelationSpace(
      /*numDomain=*/isl_basic_map_dim(bmap, isl_dim_in),
      /*numRange=*/isl_basic_map_dim(bmap, isl_dim_out),
      /*numSymbols=*/isl_basic_map_dim(bmap, isl_dim_param),
      /*numLocals=*/isl_basic_map_dim(bmap, isl_dim_div));
  IntegerRelation result(
      /*numReservedInequalities=*/isl_mat_rows(ineqMat),
      /*numReservedEqualities=*/isl_mat_rows(eqMat),
      /*numReservedCols=*/isl_mat_cols(eqMat), fplSpace);

  populateConstraints(result, eqMat, /*eq=*/true);
  populateConstraints(result, ineqMat, /*eq=*/false);

  isl_mat_free(eqMat);
  isl_mat_free(ineqMat);

  return result;
}

Then populateConstraints iterates and copies the coefficients one by one.

void populateConstraints(IntegerRelation& rel, __isl_keep isl_mat* mat,
                         bool eq) {
  unsigned numRows = isl_mat_rows(mat);
  unsigned numCols = isl_mat_cols(mat);

  for (unsigned i = 0; i < numRows; i++) {
    SmallVector<int64_t, 8> row;
    for (unsigned j = 0; j < numCols; j++) {
      isl_val* val = isl_mat_get_element_val(mat, i, j);
      row.push_back(isl_val_get_num_si(val));
      isl_val_free(val);
    }

    if (eq) {
      rel.addEquality(row);
    } else {
      rel.addInequality(row);
    }
  }
}

Similary, the converse conversion

__isl_give isl_mat* createConstraintRows(__isl_keep isl_ctx* ctx,
                                         const IntegerRelation& rel,
                                         bool isEq) {
  unsigned numRows = isEq ? rel.getNumEqualities() : rel.getNumInequalities();
  unsigned numDimVars = rel.getNumDimVars();
  unsigned numLocalVars = rel.getNumLocalVars();
  unsigned numSymbolVars = rel.getNumSymbolVars();
  unsigned numCols = rel.getNumCols();
  isl_mat* mat = isl_mat_alloc(ctx, numRows, numCols);

  for (unsigned i = 0; i < numRows; i++) {
    auto row = isEq ? rel.getEquality(i) : rel.getInequality(i);

    // Dims stay at the same positions.
    for (unsigned j = 0; j < numDimVars; j++)
      mat = isl_mat_set_element_si(mat, i, j, (int64_t)row[j]);
    // Output locals before symbols.
    for (unsigned j = 0; j < numLocalVars; j++)
      mat = isl_mat_set_element_si(
          mat, i, j + numDimVars, (int64_t)row[j + numDimVars + numSymbolVars]);
    // Output symbols in the end.
    for (unsigned j = 0; j < numSymbolVars; j++)
      mat = isl_mat_set_element_si(mat, i, j + numDimVars + numLocalVars,
                                   (int64_t)row[j + numDimVars]);
    // Finally output the constant.
    mat = isl_mat_set_element_si(
         mat, i, numCols - 1, (int64_t)row[numCols - 1]);
  }
  return mat;
}

__isl_give isl_basic_map* convertRelationToBasicMap(const IntegerRelation& rel,
                                                    __isl_keep isl_ctx* ctx) {
  isl_mat* eqMat = createConstraintRows(ctx, rel, /*isEq=*/true);
  isl_mat* ineqMat = createConstraintRows(ctx, rel, /*isEq=*/false);
  isl_space* space =
      isl_space_alloc(ctx, rel.getNumSymbolVars(), rel.getNumDomainVars(),
                      rel.getNumRangeVars());
  return isl_basic_map_from_constraint_matrices(
      space, eqMat, ineqMat, isl_dim_in, isl_dim_out, isl_dim_div,
      isl_dim_param, isl_dim_cst);
}

Next steps

In future articles I’d like to take a closer look at some of the core algorithms underlying polyhedral optimization, and ISL in particular. Many of these algorithms utilize some deep work in convex optimization, including simplex-like methods, Gomory cuts, Fourier-Motzkin elimination, and basis reduction. Chapter 2 of the ISL reference documentation gives a taste of what is involved.

Thanks to Asra Ali for her feedback on a draft of this article.