diff --git a/example/06_array_data_dependency_detection.ipynb b/example/06_array_data_dependency_detection.ipynb
new file mode 100644
index 000000000..445f9e254
--- /dev/null
+++ b/example/06_array_data_dependency_detection.ipynb
@@ -0,0 +1,625 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# How to use the data dependency checker for static array access\n",
+ "\n",
+ "This is a small introduction notebook on the use of dependency detection inside of loop enviroments for static [(1)](#fn1) array acesses. The intention is to show the required steps to setup the dependency problem, how to solve it and interpret the result.\n",
+ "\n",
+ "Lets start by parsing the file `src/phys_mod.F90` from the `example` directory and let us consider the simple routine `phys_kernel_LITE_LOOP`.\n",
+ "\n",
+ "[(1)](#fn1-back) \"static access is an array reference at a particular location in a programm\" - Compilers: Principles, Techniques, and Tools "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "[Loki::Sourcefile] Constructed from src/phys_mod.F90 in 1.40s\n"
+ ]
+ }
+ ],
+ "source": [
+ "from loki import Sourcefile\n",
+ "\n",
+ "source = Sourcefile.from_file(\"src/phys_mod.F90\", preprocess=True)\n",
+ "routine = source[\"phys_kernel_LITE_LOOP\"]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Let us examine that subroutine, it consists of a nested loop of depth two with two assignments in the inner most loop body. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ " \n",
+ " \n",
+ " \n",
+ " "
+ ]
+ }
+ ],
+ "source": [
+ "routine.body.view()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Before we consider the actual iteration space (i.e. space spanned by the loop indicies) we define some small utility functions. One simple interpreting anything as a column vector and one pretty printing a system of equations for better viewing experience."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from numpy import array, isscalar\n",
+ "\n",
+ "\n",
+ "def as_column(element):\n",
+ " return array(element).reshape(-1, 1)\n",
+ "\n",
+ "\n",
+ "def pprint_equation(*args):\n",
+ " def make_strings_equal_length(strings):\n",
+ " max_length = max(len(s) for s in strings)\n",
+ " padded_strings = [s.ljust(max_length) for s in strings]\n",
+ " return padded_strings\n",
+ "\n",
+ " collection = [\n",
+ " str(element).split(\"\\n\") if not isscalar(element) else str(element)\n",
+ " for element in args\n",
+ " ]\n",
+ " count_rows_per_element = [\n",
+ " len(element) if isinstance(element, list) else 1 for element in collection\n",
+ " ]\n",
+ "\n",
+ " max_number_lines = max(count_rows_per_element)\n",
+ " collection = [\n",
+ " [\" \"] * ((max_number_lines - len(element) + 1) // 2)\n",
+ " + element\n",
+ " + [\" \"] * ((max_number_lines - len(element)) // 2)\n",
+ " if isinstance(element, list)\n",
+ " else element\n",
+ " for index, element in enumerate(collection)\n",
+ " ]\n",
+ "\n",
+ " collection = [\n",
+ " make_strings_equal_length(element) if isinstance(element, list) else element\n",
+ " for element in collection\n",
+ " ]\n",
+ "\n",
+ " for i in range(max_number_lines):\n",
+ " for element in collection:\n",
+ " if isinstance(element, list):\n",
+ " print(element[i], end=\" \")\n",
+ " else:\n",
+ " if i == max_number_lines // 2:\n",
+ " print(element, end=\" \")\n",
+ " else:\n",
+ " print(\" \" * len(element), end=\" \")\n",
+ " print(\"\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Additionaly we create a small loop extractor, selecting the loops based on their nested structure."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from loki import FindNodes, Loop\n",
+ "\n",
+ "\n",
+ "def simple_loop_extractor(start_node):\n",
+ " \"\"\"Find all loops in the AST and structure them depending on their nesting level\"\"\"\n",
+ " start_loops = FindNodes(Loop, greedy=True).visit(start_node)\n",
+ " return [FindNodes(Loop).visit(node) for node in start_loops]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Extracting the iteration space is an easy process. We know that all our assignments which are worth comparing are inside the deepest loop body therfore simply collecting all nested loops is enough and transforming them via the polyhedron utility into the required inequality relation. For a more complex situation view further down \"Collecting access in different loop bodies\""
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from loki.analyse.util_polyhedron import Polyhedron\n",
+ "\n",
+ "nested_loops = list(simple_loop_extractor(routine.body)[0])\n",
+ "poly = Polyhedron.from_nested_loops(nested_loops)\n",
+ "B, b = poly.A, as_column(poly.b)\n",
+ "iteration_space_required_variables = list(str(v) for v in poly.variables)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ " [['k'] \n",
+ "[[-1 0 0 0 0] ['i'] [[-1] \n",
+ " [ 1 0 -1 0 0] * ['dim2'] ≤ [ 0] \n",
+ " [ 0 -1 0 1 0] ['i1'] [ 0] \n",
+ " [ 0 1 0 0 -1]] ['i2']] [ 0]] \n"
+ ]
+ }
+ ],
+ "source": [
+ "pprint_equation(B, \"*\", as_column(iteration_space_required_variables), \"≤\", b)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Next is the collection of the actual array access that we would like to compare. With the knowledge of the subroutine we deal with only two assignments, by traversing the IR and collecting all assignments."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Write into out1(i,k): \n",
+ "\tAssignment:: out1(i, k) = (in1(i, k) + in2(i, k) + in3(i, k) + in4(i, k) + in5(i, k) + in6(i, k) + in7(i, k) + in8(i, k) + in9(i, k) + in10(i, k))*0.1\n",
+ "Read of out1(i,k): \n",
+ "\tAssignment:: in1(i, k) = out1(i, k)\n"
+ ]
+ }
+ ],
+ "source": [
+ "from loki import Assignment\n",
+ "\n",
+ "assignments = FindNodes(Assignment).visit(routine.body)\n",
+ "\n",
+ "first_assignement, second_assignment = (*assignments,)\n",
+ "\n",
+ "print(\n",
+ " f\"Write into out1(i,k): \\n\\t{first_assignement}\",\n",
+ " f\"Read of out1(i,k): \\n\\t{second_assignment}\",\n",
+ " sep=\"\\n\",\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The task of how to select which array access to compare with which is currently not automatized, therefore the user is required to produce a reasonable. Here we compare the write into `out1` of the first assignment with the read of it in the second assignment, which represents a **True Dependency**, i.e. a write which is folled by a read of the same location. Even though the overall relationship is more complex, since the read of `out1` is stored in `in1` which then may be used in the some future iteration.\n",
+ "\n",
+ "Now we need to describe this access in their affine array access function representation, i.e. as a matrix vector product with an additional added vector. This is done by using the `construct_affine_array_access_function_representation` function which takes a tuple describing the dimensions access of an array and a list of to be considered variables."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from loki.analyse.analyse_array_data_dependency_detection import (\n",
+ " construct_affine_array_access_function_representation,\n",
+ ")\n",
+ "\n",
+ "(\n",
+ " first_access_function,\n",
+ " access_variables,\n",
+ ") = construct_affine_array_access_function_representation(\n",
+ " first_assignement.lhs.dimensions, iteration_space_required_variables\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Resulting in this representation:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ " [['k'] \n",
+ " ['i'] \n",
+ "[[0 1 0 0 0] * ['dim2'] + [[0] = 0 \n",
+ " [1 0 0 0 0]] ['i1'] [0]] \n",
+ " ['i2']] \n"
+ ]
+ }
+ ],
+ "source": [
+ "pprint_equation(first_access_function.matrix, \"*\", as_column(access_variables), \"+\", first_access_function.vector, \"=\", 0)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The same is done for the second access."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "(\n",
+ " second_access_function,\n",
+ " access_variables_dash,\n",
+ ") = construct_affine_array_access_function_representation(\n",
+ " second_assignment.rhs.dimensions, iteration_space_required_variables\n",
+ ")\n",
+ "\n",
+ "assert access_variables == access_variables_dash"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "With the representation:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ " [['k'] \n",
+ " ['i'] \n",
+ "[[0 1 0 0 0] * ['dim2'] + [[0] = 0 \n",
+ " [1 0 0 0 0]] ['i1'] [0]] \n",
+ " ['i2']] \n"
+ ]
+ }
+ ],
+ "source": [
+ "pprint_equation(second_access_function.matrix, \"*\", as_column(access_variables_dash), \"+\", second_access_function.vector, \"=\", 0)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now we have a iteration space per function access, here it is the same iteration space since all function accesses happen in the inner most loop. Passing this combination to the `has_data_dependency` function we now emply various techniques (gcd test, independet variable test, integer linear programming) to confirm or deny a data dependency."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "There is a data dependency between the two accesses.\n"
+ ]
+ }
+ ],
+ "source": [
+ "from loki.analyse.analyse_array_data_dependency_detection import has_data_dependency\n",
+ "\n",
+ "first_access_represenetation = ((B, b), first_access_function)\n",
+ "second_access_represenetation = ((B, b), second_access_function)\n",
+ "\n",
+ "if has_data_dependency(first_access_represenetation, second_access_represenetation):\n",
+ " print(\"There is a data dependency between the two accesses.\")\n",
+ "else:\n",
+ " print(\"There is no data dependency between the two accesses.\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "# Collecting access in different loop bodies"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "For this we consider a different subroutine, with two independent loops one iterating over all even integers and one iterating over all odd, clearly showing no data dependency."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "[Loki::Sourcefile] Constructed from ../tests/sources/data_dependency_detection/loop_carried_dependencies.f90 in 0.19s\n"
+ ]
+ }
+ ],
+ "source": [
+ "source = Sourcefile.from_file(\n",
+ " \"../tests/sources/data_dependency_detection/loop_carried_dependencies.f90\",\n",
+ " preprocess=True,\n",
+ ")\n",
+ "routine = source[\"NoDependency\"]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " \n",
+ " "
+ ]
+ }
+ ],
+ "source": [
+ "routine.body.view()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Here the two different loops are collected greedily and then **two** different iteration spaces are constructed in the same way as above."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "outer_loops = FindNodes(Loop, greedy=True).visit(routine.body)\n",
+ "\n",
+ "first_iteration_space = Polyhedron.from_nested_loops([outer_loops[0]])\n",
+ "second_iteration_space = Polyhedron.from_nested_loops([outer_loops[1]])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[-1] [[-1] \n",
+ " [ 1]] * [['i']] ≤ [10]] \n"
+ ]
+ }
+ ],
+ "source": [
+ "B, b = first_iteration_space.A, as_column(first_iteration_space.b)\n",
+ "iteration_space_required_variables = [str(v) for v in first_iteration_space.variables]\n",
+ "\n",
+ "pprint_equation(B, \"*\", as_column(iteration_space_required_variables), \"≤\", b)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[-1] [[-1] \n",
+ " [ 1]] * [['i']] ≤ [ 5]] \n"
+ ]
+ }
+ ],
+ "source": [
+ "B_dash, b_dash = second_iteration_space.A, as_column(second_iteration_space.b)\n",
+ "iteration_space_required_variables_dash = [\n",
+ " str(v) for v in second_iteration_space.variables\n",
+ "]\n",
+ "\n",
+ "assert iteration_space_required_variables == iteration_space_required_variables_dash\n",
+ "\n",
+ "pprint_equation(\n",
+ " B_dash, \"*\", as_column(iteration_space_required_variables_dash), \"≤\", b_dash\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Again the assignments are collected, then the access on the `data` array is represented for the two different assignments."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Write into data(even values): \n",
+ "\tAssignment:: data(2*i) = 10\n",
+ "Write into data(odd values): \n",
+ "\tAssignment:: data(2*i + 1) = 20\n"
+ ]
+ }
+ ],
+ "source": [
+ "assignments = FindNodes(Assignment).visit(routine.body)\n",
+ "\n",
+ "first_assignement, second_assignment = (*assignments,)\n",
+ "\n",
+ "print(\n",
+ " f\"Write into data(even values): \\n\\t{first_assignement}\",\n",
+ " f\"Write into data(odd values): \\n\\t{second_assignment}\",\n",
+ " sep=\"\\n\",\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[2]] * [['i']] + [[0]] = 0 \n"
+ ]
+ }
+ ],
+ "source": [
+ "first_access, access_variables = construct_affine_array_access_function_representation(\n",
+ " first_assignement.lhs.dimensions, iteration_space_required_variables\n",
+ ")\n",
+ "\n",
+ "pprint_equation(first_access.matrix, \"*\", as_column(access_variables), \"+\", first_access.vector, \"=\", 0)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[2]] * [['i']] + [[1]] = 0 \n"
+ ]
+ }
+ ],
+ "source": [
+ "(\n",
+ " second_access,\n",
+ " access_variables_dash,\n",
+ ") = construct_affine_array_access_function_representation(\n",
+ " second_assignment.lhs.dimensions, iteration_space_required_variables\n",
+ ")\n",
+ "\n",
+ "assert access_variables == access_variables_dash\n",
+ "\n",
+ "pprint_equation(second_access.matrix, \"*\", as_column(access_variables_dash), \"+\", second_access.vector, \"=\", 0)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "And last but not least the question of a data dependency is answered."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "There is no data dependency between the two accesses.\n"
+ ]
+ }
+ ],
+ "source": [
+ "first_access_represenetation = ((B, b), first_access)\n",
+ "second_access_represenetation = ((B_dash, b_dash), second_access)\n",
+ "\n",
+ "if has_data_dependency(first_access_represenetation, second_access_represenetation):\n",
+ " print(\"There is a data dependency between the two accesses.\")\n",
+ "else:\n",
+ " print(\"There is no data dependency between the two accesses.\")"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "loki_env",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.12"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/loki/analyse/analyse_array_data_dependency_detection.py b/loki/analyse/analyse_array_data_dependency_detection.py
new file mode 100644
index 000000000..7dc5c74aa
--- /dev/null
+++ b/loki/analyse/analyse_array_data_dependency_detection.py
@@ -0,0 +1,480 @@
+# (C) Copyright 2018- ECMWF.
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+
+from math import gcd as math_gcd
+from warnings import warn
+from dataclasses import dataclass
+from typing import Any, List
+import numpy as np
+import numpy.typing as npt
+from loki.expression import (
+ FindVariables,
+ accumulate_polynomial_terms,
+ simplify,
+ is_constant,
+ symbols as sym,
+)
+from loki.analyse.util_linear_algebra import (
+ generate_row_echelon_form,
+ back_substitution,
+ is_independent_system,
+ yield_one_d_systems,
+)
+
+try:
+ _ = math_gcd(4, 3, 2)
+ gcd = math_gcd
+except TypeError: # Python 3.8 can only handle two arguments
+ from functools import reduce
+
+ def gcd(*args):
+ return reduce(math_gcd, args)
+
+
+try:
+ from ortools.linear_solver import pywraplp
+
+ HAVE_ORTOOLS = True
+except ImportError:
+ HAVE_ORTOOLS = False
+
+__all__ = [
+ "has_data_dependency",
+ "HAVE_ORTOOLS",
+ "construct_affine_array_access_function_representation",
+]
+
+NDArrayInt = npt.NDArray[np.int_]
+
+
+@dataclass
+class IterationSpaceRepresentation:
+ """
+ Represents an iteration space as an inequality system. The matrix A, vector b and unkowns x
+ represent the inequality system Ax <= b, where the inequality is to be interpreted element wise.
+ """
+
+ matrix: NDArrayInt
+ vector: NDArrayInt
+
+
+# abbrivation for readability
+IterSpacRepr = IterationSpaceRepresentation
+
+
+@dataclass
+class AffineArrayAccessFunctionRepresentation:
+ """
+ Represents an affine array access function as a linear system. The matrix F and vector f
+ represent the linear system Fx + f, where x is the iteration space vector of unkowns.
+ """
+
+ matrix: NDArrayInt
+ vector: NDArrayInt
+
+
+# abbrivation for readability
+AffiAcceRepr = AffineArrayAccessFunctionRepresentation
+
+
+@dataclass
+class StaticArrayAccessInstanceRepresentation:
+ """
+ Represents a static array access instance using an iteration space representation and
+ an affine array access function representation.
+ """
+
+ iteration_space: IterSpacRepr
+ access_function: AffiAcceRepr
+
+
+# abbrivation for readability
+StatInstRepr = StaticArrayAccessInstanceRepresentation
+
+
+def _assert_correct_dimensions(
+ first_array_access: StatInstRepr,
+ second_array_access: StatInstRepr,
+) -> None:
+ def assert_compatibility_of_iteration_space(
+ first: IterSpacRepr, second: IterSpacRepr
+ ):
+ (number_inequalties_1, number_variables_1) = first.matrix.shape
+ (rhs_number_inequalities_1, column_count) = first.vector.shape
+ assert column_count == 1, "Right hand side is a single column vector!"
+ assert (
+ number_inequalties_1 == rhs_number_inequalities_1
+ ), "System matrix and right hand side vector require same amount of rows!"
+
+ (number_inequalties_2, number_variables_2) = second.matrix.shape
+ (rhs_number_inequalities_2, column_count) = first.vector.shape
+ assert column_count == 1, "Right hand side is a single column vector!"
+ assert (
+ number_inequalties_2 == rhs_number_inequalities_2
+ ), "System matrix and right hand side vector require same amount of rows!"
+
+ assert number_variables_1 == number_variables_2, (
+ "Same number of variables per iteration space is assumed with "
+ "each variable from the first iteration space having its equivalent "
+ "at the same place in the second iteration space"
+ )
+
+ def assert_compatibility_of_access_function(
+ first: AffiAcceRepr,
+ second: AffiAcceRepr,
+ ):
+ (number_rows_1, number_columns_1) = first.matrix.shape
+ (rhs_number_rows_1, column_count) = first.vector.shape
+ assert column_count == 1, "Right hand side is a single column vector!"
+ assert (
+ number_rows_1 == rhs_number_rows_1
+ ), "System matrix and right hand side vector require same amount of rows!"
+
+ (number_rows_2, number_columns_2) = second.matrix.shape
+ (rhs_number_rows_2, column_count) = first.vector.shape
+ assert column_count == 1, "Right hand side is a single column vector!"
+ assert (
+ number_rows_2 == rhs_number_rows_2
+ ), "System matrix and right hand side vector require same amount of rows!"
+
+ assert number_columns_1 == number_columns_2, (
+ "Same number of variables per access function is assumed with "
+ "each variable from the first access function having its equivalent "
+ "at the same place in the second access function"
+ )
+
+ def assert_compatibility_of_access_and_iteration_unkowns(
+ function: AffiAcceRepr,
+ space: IterSpacRepr,
+ ):
+ (_, number_variables_space) = space.matrix.shape
+ (_, number_variables_function) = function.matrix.shape
+
+ assert (
+ number_variables_space == number_variables_function
+ ), "Same number of variables per access function and iteration space is assumed"
+
+ assert_compatibility_of_iteration_space(
+ first_array_access.iteration_space, second_array_access.iteration_space
+ )
+ assert_compatibility_of_access_function(
+ first_array_access.access_function, second_array_access.access_function
+ )
+ assert_compatibility_of_access_and_iteration_unkowns(
+ first_array_access.access_function, second_array_access.iteration_space
+ )
+
+
+def _safe_integer_division(x, y):
+ result = x // y
+ if (x % y != 0).all():
+ raise ValueError("Division does not result in an integer.")
+ return result
+
+
+def _gaussian_eliminiation_for_diophantine_equations(
+ augmented_matrix: NDArrayInt,
+) -> (bool, NDArrayInt):
+ """
+ Compute the Row Echelon Form (REF) of an augmented matrix while ensuring integer solutions
+ to linear Diophantine equations.
+
+ Parameters
+ ----------
+ augmented_matrix : numpy.ndarray
+ The input augmented system matrix.
+
+ Returns
+ -------
+ tuple
+ HasIntegerSolution : bool
+ A boolean indicating if a solution with integer coefficients exists.
+ REF_matrix : numpy.ndarray
+ The Row Echelon Form (REF) of the input matrix.
+
+ Notes
+ -----
+ This function calculates the Row Echelon Form (REF) of a given matrix of integers while enforcing that
+ each linear Diophantine equation in the system has integer solutions. It follows Theorem 11.32
+ in "Compilers: Principles, Techniques, and Tools" by Aho, Lam, Sethi, and Ullman (2015).
+ """
+
+ def gcd_condition(A):
+ """Check that gcd condition of linear Diophantine equation is satisfied"""
+ if A[0, -1] % gcd(*A[0, :-1]) != 0:
+ raise ValueError
+
+ status = True
+ solution = None
+ try:
+ solution = generate_row_echelon_form(
+ augmented_matrix,
+ conditional_check=gcd_condition,
+ division_operator=_safe_integer_division,
+ )
+ except ValueError:
+ status = False
+
+ return (status, solution)
+
+
+def _does_independent_system_violate_bounds(
+ matrix: NDArrayInt, vector: NDArrayInt
+) -> bool:
+ """
+ Check if a system of inequalities, which can be separated into N one-dimensional problems
+ represented by a matrix and a vector, violates its bounds. The system is defined as A x <= b,
+ where A is the matrix, b is the vector, and x is the vector of unknowns.
+
+ Parameters
+ ----------
+ matrix : numpy.ndarray
+ A 2D NumPy array representing the matrix of the independent system.
+
+ vector : numpy.ndarray
+ A 1D NumPy array representing the vector of the independent system.
+
+ Returns
+ -------
+ bool
+ True if any of the 1D systems violate their bounds, False otherwise.
+ """
+
+ def get_bounds(matrix, vector):
+ """Gets x <= c and x >= d conditions"""
+ larger_zero = matrix > 0
+ upper_bounds = vector[larger_zero] / matrix[larger_zero]
+
+ smaller_zero = matrix < 0
+ lower_bounds = vector[smaller_zero] / matrix[smaller_zero]
+
+ return np.unique(lower_bounds), np.unique(upper_bounds)
+
+ for A, b in yield_one_d_systems(matrix, vector):
+ if np.all(0 == A):
+ if not np.all(0 <= b[A == 0]):
+ return True
+
+ lower_bounds, upper_bounds = get_bounds(A.astype(float), b.astype(float))
+ if np.amax(lower_bounds) > np.amin(upper_bounds):
+ return True
+ return False
+
+
+def _solve_inequalties_with_ortools(
+ constraint_coeffs: NDArrayInt, right_hand_side: NDArrayInt
+):
+ """
+ Solve a system of inequalities using ortools. With A the constraint coefficient matrix, b the right
+ hand side and x the unkowns, the system is of the form A x <= b.
+ """
+ number_constraints, number_variables = constraint_coeffs.shape
+ assert (
+ number_constraints,
+ 1,
+ ) == right_hand_side.shape, (
+ "number of number_constraints, "
+ "i.e. the number of rows in matrix must match number of columns in right hand side vector"
+ )
+
+ solver = pywraplp.Solver.CreateSolver("SCIP")
+ assert solver, "Solver could not be created."
+
+ inf = solver.infinity()
+ unkowns = [solver.IntVar(-inf, inf, "x_" + str(i)) for i in range(number_variables)]
+
+ for i in range(number_constraints):
+ constraint_expr = [
+ constraint_coeffs[i][j] * unkowns[j] for j in range(number_variables)
+ ]
+ solver.Add(sum(constraint_expr) <= right_hand_side[i, 0])
+
+ solver.Maximize(0) # arbitrary simple objective function
+
+ return solver.Solve()
+
+
+def has_data_dependency_impl(
+ first_array_access: StatInstRepr,
+ second_array_access: StatInstRepr,
+) -> bool:
+ # fmt: off
+ (has_solution, solved_system) = _gaussian_eliminiation_for_diophantine_equations(
+ np.hstack(
+ (
+ first_array_access.access_function.matrix,
+ -second_array_access.access_function.matrix,
+ first_array_access.access_function.vector - second_array_access.access_function.vector,
+ )
+ )
+ )
+ # fmt: on
+ if not has_solution:
+ return False
+
+ n, m = solved_system.shape
+
+ assert n < m, "The augmented system is assumed to be underdetermined."
+
+ R, S, g_hat = solved_system[:, :n], solved_system[:, n:-1], solved_system[:, -1:]
+
+ # fmt: off
+ K = np.block(
+ [
+ [
+ first_array_access.iteration_space.matrix, np.zeros_like(first_array_access.iteration_space.matrix),
+ ],
+ [
+ np.zeros_like(second_array_access.iteration_space.matrix), second_array_access.iteration_space.matrix,
+ ],
+ ]
+ )
+ # fmt: on
+ k = np.vstack(
+ (
+ first_array_access.iteration_space.vector,
+ second_array_access.iteration_space.vector,
+ )
+ )
+
+ n, m = R.shape
+ K_R = K[0:n, 0:m]
+
+ U = K[n:, 0:m]
+ V = K[0:n, m:]
+ Z = K[n:, m:]
+
+ tmp1 = back_substitution(R, S, _safe_integer_division)
+ tmp2 = back_substitution(R, g_hat, _safe_integer_division)
+
+ matrix_final_polygon = np.vstack((-K_R.dot(tmp1) + V, -U.dot(tmp1) + Z))
+ vector_final_polygon = k + np.vstack((K_R.dot(tmp2), U.dot(tmp2)))
+
+ if is_independent_system(
+ matrix_final_polygon
+ ) and _does_independent_system_violate_bounds(
+ matrix_final_polygon, vector_final_polygon
+ ):
+ return False
+
+ if not HAVE_ORTOOLS:
+ warn(
+ "ortools is not installed --> assuming data dependency --> will lead to alot of false positives"
+ )
+ return True
+
+ status = _solve_inequalties_with_ortools(matrix_final_polygon, vector_final_polygon)
+
+ known_status = {
+ pywraplp.Solver.OPTIMAL: True,
+ pywraplp.Solver.FEASIBLE: True,
+ pywraplp.Solver.INFEASIBLE: False,
+ }
+
+ if status in known_status:
+ return known_status[status]
+ warn("ortools produced a not considered status --> assuming data dependency.")
+
+ return True # assume data dependency if all tests are inconclusive
+
+
+def _ensure_correct_StatInstRepr(array_access: Any) -> StatInstRepr:
+ if isinstance(array_access, StatInstRepr):
+ return array_access
+
+ if len(array_access) == 2:
+ iteration_space = array_access[0]
+ access_function = array_access[1]
+
+ if isinstance(iteration_space, IterSpacRepr) and isinstance(
+ access_function, AffiAcceRepr
+ ):
+ return StatInstRepr(iteration_space, access_function)
+
+ if isinstance(iteration_space, IterSpacRepr) and len(access_function) == 2:
+ return StatInstRepr(iteration_space, AffiAcceRepr(*access_function))
+
+ if len(iteration_space) == 2 and isinstance(access_function, AffiAcceRepr):
+ return StatInstRepr(IterSpacRepr(*iteration_space), access_function)
+
+ if len(iteration_space) == 2 and len(access_function) == 2:
+ return StatInstRepr(
+ IterSpacRepr(*iteration_space), AffiAcceRepr(*access_function)
+ )
+
+ raise ValueError(f"Cannot convert {str(array_access)} to StatInstRepr")
+
+
+def has_data_dependency(
+ first_array_access: StatInstRepr, second_array_access: StatInstRepr
+):
+ first_array_access = _ensure_correct_StatInstRepr(first_array_access)
+ second_array_access = _ensure_correct_StatInstRepr(second_array_access)
+
+ _assert_correct_dimensions(first_array_access, second_array_access)
+
+ return has_data_dependency_impl(first_array_access, second_array_access)
+
+
+def construct_affine_array_access_function_representation(
+ array_dimensions_expr: tuple(), additional_variables: List[str] = None
+):
+ """
+ Create a matrix and vector representation of the array access function.
+
+ This function is used to represent array access expressions like z[i], where you should pass ("i", )
+ or y[1+3, 4-j], where you should pass ("1+3", "4-j"). If 'var' is an Array object, then 'var.dimensions'
+ should be passed to this function.
+
+ Returns:
+ AffiAcceRepr(matrix, vector): A mapping that relates a vector 'i' within the bounds Bi+b>=0 to the
+ array location Fi+f.
+ """
+
+ def generate_row(expr, variables):
+ supported_types = (sym.TypedSymbol, sym.MetaSymbol, sym.Sum, sym.Product)
+ if not (is_constant(expr) or isinstance(expr, supported_types)):
+ raise ValueError(f"Cannot derive inequality from expr {str(expr)}")
+ simplified_expr = simplify(expr)
+ terms = accumulate_polynomial_terms(simplified_expr)
+ const_term = terms.pop(1, 0) # Constant term or 0
+ row = np.zeros(len(variables), dtype=np.dtype(int))
+
+ for base, coef in terms.items():
+ if not len(base) == 1:
+ raise ValueError(f"Non-affine bound {str(simplified_expr)}")
+ row[variables.index(base[0].name.lower())] = coef
+
+ return row, const_term
+
+ def unique_order_preserving(sequence):
+ seen = set()
+ return [x for x in sequence if not (x in seen or seen.add(x))]
+
+ if additional_variables is None:
+ additional_variables = []
+
+ for variable in additional_variables:
+ assert variable.lower() == variable
+
+ variables = additional_variables.copy()
+ variables += list(
+ {v.name.lower() for v in FindVariables().visit(array_dimensions_expr)}
+ )
+ variables = unique_order_preserving(variables)
+
+ n = len(array_dimensions_expr)
+ d = len(variables)
+
+ F = np.zeros([n, d], dtype=np.dtype(int))
+ f = np.zeros([n, 1], dtype=np.dtype(int))
+
+ for dimension_index, sub_expr in enumerate(array_dimensions_expr):
+ row, constant = generate_row(sub_expr, variables)
+ F[dimension_index] = row
+ f[dimension_index, 0] = constant
+ return AffiAcceRepr(F, f), variables
diff --git a/loki/analyse/util_linear_algebra.py b/loki/analyse/util_linear_algebra.py
index 2d1485afa..2c81f13a2 100644
--- a/loki/analyse/util_linear_algebra.py
+++ b/loki/analyse/util_linear_algebra.py
@@ -5,9 +5,7 @@
# granted to it by virtue of its status as an intergovernmental organisation
# nor does it submit to any jurisdiction.
-from numpy import zeros_like, dot
-from numpy import vstack, hstack
-from numpy import all as np_all, sum as np_sum, isin as np_isin
+import numpy as np
__all__ = [
"back_substitution",
@@ -21,13 +19,19 @@ def is_independent_system(matrix):
"""
Check if a linear system of equations can be split into independent one-dimensional problems.
- Args:
- matrix (numpy.ndarray): A rectangular matrix representing coefficients.
+ Parameters
+ ----------
+ matrix : numpy.ndarray
+ A rectangular matrix representing coefficients.
- Returns:
- bool: True if the system can be split into independent one-dimensional problems, False otherwise.
+ Returns
+ -------
+ bool
+ True if the system can be split into independent one-dimensional problems, False otherwise.
- This function checks whether a linear system of equations in the form of matrix x [operator] right_hand_side
+ Notes
+ -----
+ This function checks whether a linear system of equations in the form of `matrix [operator] right_hand_side`
can be split into independent one-dimensional problems. The number of problems is determined by the
number of variables (the row number of the matrix).
@@ -35,61 +39,58 @@ def is_independent_system(matrix):
if each row of the matrix has exactly one non-zero coefficient or no non-zero coefficients.
"""
- return np_all(np_isin(np_sum(matrix != 0, axis=1), [0, 1]))
+ return np.all(np.isin(np.sum(matrix != 0, axis=1), [0, 1]))
def yield_one_d_systems(matrix, right_hand_side):
"""
- Split a linear system of equations (<=, >= or ==) into independent one-dimensional problems.
-
- Args:
- matrix (numpy.ndarray): A rectangular matrix representing coefficients.
- right_hand_side (numpy.ndarray): The right-hand side vector.
-
- Yields:
- tuple[single_dimensional_array, single_dimensional_array]:
- A tuple containing a coefficient vector and the corresponding right-hand side.
-
- This function takes a linear system of equations in the form of matrix x [operator] right_hand_side,
- where "matrix" is a rectangular matrix, x is a vector of variables, and "right_hand_side" is
+ Split a linear system of equations (<=, >=, or ==) into independent one-dimensional problems.
+
+ Parameters
+ ----------
+ matrix : numpy.ndarray
+ A rectangular matrix representing coefficients.
+ right_hand_side : numpy.ndarray
+ The right-hand side vector.
+
+ Yields
+ ------
+ tuple[numpy.ndarray, numpy.ndarray]
+ A tuple containing a coefficient vector and the corresponding right-hand side.
+
+ Notes
+ -----
+ The independence of the problems is NOT explicitly checked; call `is_independent_system` before using this
+ function if unsure.
+
+ This function takes a linear system of equations in the form of `matrix [operator] right_hand_side`,
+ where "matrix" is a rectangular matrix, "x" is a vector of variables, and "right_hand_side" is
the right-hand side vector. It splits the system into assumed independent one-dimensional problems.
Each problem consists of a coefficient vector and a right-hand side. The number of problems is equal to the
number of variables (the row number of the matrix).
- Note:
- - The independence of the problems is not explicitly checked, call is_independent_system before using this
- function if unsure.
-
- Example:
- ```
+ Example
+ -------
+ ```python
for A, b in yield_one_d_systems(matrix, right_hand_side):
# Solve the one-dimensional problem A * x = b
solution = solve_one_d_system(A, b)
```
"""
- # drop completly empty rows
- mask = ~np_all(hstack((matrix, right_hand_side)) == 0, axis=1)
- matrix = matrix[mask]
- right_hand_side = right_hand_side[mask]
-
# yield systems with empty left hand side (A) and non empty right hand side
- mask = np_all(matrix == 0, axis=1)
- if right_hand_side[mask].size == 0:
- return
-
- for A, b in zip(matrix[mask].T, right_hand_side[mask].T):
- yield A, b
+ mask = np.all(matrix == 0, axis=1)
+ if right_hand_side[mask].size != 0:
+ for A in matrix[mask].T:
+ yield A.reshape((-1,1)), right_hand_side[mask]
matrix = matrix[~mask]
right_hand_side = right_hand_side[~mask]
- if right_hand_side.size == 0:
- return
-
- for A, b in zip(matrix.T, right_hand_side.T):
- mask = A != 0
- yield A[mask], b[mask]
+ if right_hand_side.size != 0:
+ for A in matrix.T:
+ mask = A != 0
+ yield A[mask].reshape((-1,1)), right_hand_side[mask]
def back_substitution(
@@ -100,33 +101,42 @@ def back_substitution(
"""
Solve a linear system of equations using back substitution for an upper triangular square matrix.
- Args:
- upper_triangular_square_matrix (numpy.ndarray): An upper triangular square matrix (R).
- right_hand_side (numpy.ndarray): A vector (y) on the right-hand side of the equation Rx = y.
- division_operation (function, optional): A custom division operation function. Default is standard division (/).
+ Parameters
+ ----------
+ upper_triangular_square_matrix : numpy.ndarray
+ An upper triangular square matrix (R).
+
+ right_hand_side : numpy.ndarray
+ A vector (y) on the right-hand side of the equation Rx = y.
- Returns:
- numpy.ndarray: The solution vector (x) to the system of equations Rx = y.
+ division_operation : function, optional
+ A custom division operation function. Default is standard division (/).
+ Returns
+ -------
+ numpy.ndarray
+ The solution vector (x) to the system of equations Rx = y.
+
+ Notes
+ -----
The function performs back substitution to find the solution vector x for the equation Rx = y,
where R is an upper triangular square matrix and y is a vector. The division_operation
function is used for division (e.g., for custom division operations).
- Note:
- - The function assumes that the upper right element of the upper_triangular_square_matrix (R)
- is nonzero for proper back substitution.
+ The function assumes that the upper right element of the upper_triangular_square_matrix (R)
+ is nonzero for proper back substitution.
"""
R = upper_triangular_square_matrix
y = right_hand_side
- x = zeros_like(y)
+ x = np.zeros_like(y)
assert R[-1, -1] != 0
x[-1] = divison_operation(y[-1], R[-1, -1])
for i in range(len(y) - 2, -1, -1):
- x[i] = divison_operation((y[i] - dot(R[i, i + 1 :], x[i + 1 :])), R[i, i])
+ x[i] = divison_operation((y[i] - np.dot(R[i, i + 1 :], x[i + 1 :])), R[i, i])
return x
@@ -135,25 +145,32 @@ def generate_row_echelon_form(
A, conditional_check=lambda A: None, division_operator=lambda x, y: x / y
):
"""
- Calculate the Row Echelon Form (REF) of a matrix A using Gaussian elimination.
-
- Args:
- A (numpy.ndarray): The input matrix for which RREF is to be calculated.
- conditional_check (function, optional): A custom function to check conditions during the computation.
- division_operation (function, optional): A custom division operation function. Default is standard division (/).
-
- Returns:
- numpy.ndarray: The REF of the input matrix A.
-
- This function computes the Row Echelon Form (RREF) of a given matrix A using Gaussian elimination.
- You can provide a custom conditional_check function to add checks or operations during the computation.
- You can also specify a custom division_operation function for division (e.g., for custom division operations).
-
- Note:
+ Calculate the Row Echelon Form (REF) of a matrix.
+
+ Parameters
+ ----------
+ A : numpy.ndarray
+ The input matrix for which the REF is to be calculated.
+ conditional_check : function, optional
+ A custom function to check conditions during the computation.
+ division_operation : function, optional
+ A custom division operation function. Default is standard division (/).
+
+ Returns
+ -------
+ numpy.ndarray
+ The REF of the input matrix A.
+
+ Notes
+ -----
- If the input matrix has no rows or columns, it is already in REF, and the function returns itself.
- The function utilizes the specified division operation (default is standard division) for division.
- Reference: https://math.stackexchange.com/questions/3073083/how-to-reduce-matrix-into-row-echelon-form-in-numpy
+ Reference
+ ---------
+ https://math.stackexchange.com/a/3073117
+ for question:
+ https://math.stackexchange.com/questions/3073083/how-to-reduce-matrix-into-row-echelon-form-in-numpy
"""
# if matrix A has no columns or rows,
# it is already in REF, so we return itself
@@ -170,7 +187,7 @@ def generate_row_echelon_form(
# we perform REF on matrix from second column
B = generate_row_echelon_form(A[:, 1:], conditional_check, division_operator)
# and then add the first zero-column back
- return hstack([A[:, :1], B])
+ return np.hstack([A[:, :1], B])
# if non-zero element happens not in the first row,
# we switch rows
@@ -190,4 +207,4 @@ def generate_row_echelon_form(
B = generate_row_echelon_form(A[1:, 1:], conditional_check, division_operator)
# we add first row and first (zero) column, and return
- return vstack([A[:1], hstack([A[1:, :1], B])])
+ return np.vstack([A[:1], np.hstack([A[1:, :1], B])])
diff --git a/loki/analyse/util_polyhedron.py b/loki/analyse/util_polyhedron.py
new file mode 100644
index 000000000..a001622b7
--- /dev/null
+++ b/loki/analyse/util_polyhedron.py
@@ -0,0 +1,354 @@
+# (C) Copyright 2018- ECMWF.
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+
+from typing import List
+import numpy as np
+from loki.ir import Loop
+from loki.expression import (
+ symbols as sym,
+ FindVariables,
+ simplify,
+ is_constant,
+ accumulate_polynomial_terms,
+)
+from loki.tools import as_tuple
+
+__all__ = ["Polyhedron"]
+
+
+class Polyhedron:
+ """
+ Halfspace representation of a (convex) polyhedron.
+
+ A polyhedron `P c R^d` is described by a set of inequalities, in matrix form
+ ```
+ P = { x=[x1,...,xd]^T c R^d | Ax <= b }
+ ```
+ with n-by-d matrix `A` and d-dimensional right hand side `b`.
+
+ In loop transformations, polyhedrons are used to represent iteration spaces of
+ d-deep loop nests.
+
+
+ Parameters
+ ----------
+ A : numpy.array
+ The representation matrix A.
+ b : numpy.array
+ The right hand-side vector b.
+ variables : list, optional
+ List of variables representing the dimensions in the polyhedron.
+
+ Attributes
+ ----------
+ A : numpy.array
+ The representation matrix A.
+ b : numpy.array
+ The right hand-side vector b.
+ variables : list, optional, default = None
+ """
+
+ def __init__(self, A, b, variables=None):
+ A = np.array(A, dtype=np.dtype(int))
+ b = np.array(b, dtype=np.dtype(int))
+ assert A.ndim == 2 and b.ndim == 1
+ assert A.shape[0] == b.shape[0]
+ self.A = A
+ self.b = b
+
+ self.variables = None
+ self.variable_names = None
+ if variables is not None:
+ assert len(variables) == A.shape[1]
+ self.variables = variables
+ self.variable_names = [v.name.lower() for v in self.variables]
+
+ def __str__(self):
+ str_A = "[" + ", ".join([str(row) for row in self.A]) + "]"
+ str_b = f"[{', '.join(map(str, self.b))}]"
+ str_variable_names = (
+ f"[{', '.join(map(str, self.variable_names))}]"
+ if self.variable_names
+ else "[]"
+ )
+
+ return f"Polyhedron(\n\tA={str_A}, \n\tb={str_b}, \n\tvariables={str_variable_names}\n)"
+
+ def __repr__(self):
+ return str(self)
+
+ def _has_satisfiable_constant_restrictions(self):
+ """
+ Check whether the constant restrictions of the polyhedron are satisfiable.
+
+ This method checks if 0 <= b, assuming that A x = 0.
+
+ Returns:
+ bool: True if all constant restrictions are satisfiable, False otherwise.
+
+ """
+
+ return (0 <= self.b).all()
+
+ def is_empty(self):
+ """
+ Determine whether a polyhedron is empty.
+
+ A polyhedron is considered empty under the following conditions:
+ 1. It contains no inequalities.
+ 2. It spans no space, which is a nontrivial problem. The simplest case is when it has an empty
+ matrix A and does not satisfy the constant restrictions 0 <= b.
+
+ Notes
+ -----
+ An empty polyhedron implies that it has no valid solutions or feasible points within its boundaries.
+ This function is expected to be called only for polyhedrons with an empty matrix.
+
+ Returns
+ -------
+ bool
+ True if the polyhedron is empty; False if it is not.
+ """
+ if self.A.size == 0:
+ return self.b.size == 0 or not self._has_satisfiable_constant_restrictions()
+
+ raise RuntimeError(
+ """
+ Checking if a polyhedron with a non-empty matrix spans no space is a nontrivial problem.
+ This function is expected to be only called upon polyhedrons with an empty matrix!
+ """
+ )
+
+ def variable_to_index(self, variable):
+ if self.variable_names is None:
+ raise RuntimeError("No variables list associated with polyhedron.")
+ if isinstance(variable, sym.TypedSymbol):
+ variable = variable.name.lower()
+ assert isinstance(variable, str)
+ return self.variable_names.index(variable)
+
+ @staticmethod
+ def _to_literal(value):
+ if value < 0:
+ return sym.Product((-1, sym.IntLiteral(abs(value))))
+ return sym.IntLiteral(value)
+
+ def lower_bounds(self, index_or_variable, ignore_variables=None):
+ """
+ Return all lower bounds imposed on a variable.
+
+ The lower bounds for the variable `j` are given by the index set:
+ ```
+ L = {i | A_ij < 0, i in {0, ..., d-1}}
+ ```
+ Parameters
+ ----------
+ index_or_variable : int or str or sym.Array or sym.Scalar
+ The index, name, or expression symbol for which the lower bounds are produced.
+ ignore_variables : list or None, optional
+ A list of variable names, indices, or symbols for which constraints should be ignored
+ if they depend on one of them.
+
+ Returns
+ -------
+ list
+ The bounds for the specified variable.
+ """
+ if isinstance(index_or_variable, int):
+ j = index_or_variable
+ else:
+ j = self.variable_to_index(index_or_variable)
+
+ if ignore_variables:
+ ignore_variables = [
+ i if isinstance(i, int) else self.variable_to_index(i)
+ for i in ignore_variables
+ ]
+
+ bounds = []
+ for i in range(self.A.shape[0]):
+ if self.A[i, j] < 0:
+ if ignore_variables and any(
+ self.A[i, k] != 0 for k in ignore_variables
+ ):
+ # Skip constraint that depends on any of the ignored variables
+ continue
+
+ components = [
+ self._to_literal(self.A[i, k]) * self.variables[k]
+ for k in range(self.A.shape[1])
+ if k != j and self.A[i, k] != 0
+ ]
+ if not components:
+ lhs = sym.IntLiteral(0)
+ elif len(components) == 1:
+ lhs = components[0]
+ else:
+ lhs = sym.Sum(as_tuple(components))
+ bounds += [
+ simplify(
+ sym.Quotient(
+ self._to_literal(self.b[i]) - lhs,
+ self._to_literal(self.A[i, j]),
+ )
+ )
+ ]
+ return bounds
+
+ def upper_bounds(self, index_or_variable, ignore_variables=None):
+ """
+ Return all upper bounds imposed on a variable.
+
+ The upper bounds for the variable `j` are given by the index set:
+ ```
+ U = {i | A_ij > 0, i in {0, ..., d-1}}
+ ```
+ Parameters
+ ----------
+ index_or_variable : int or str or sym.Array or sym.Scalar
+ The index, name, or expression symbol for which the upper bounds are produced.
+ ignore_variables : list or None, optional
+ A list of variable names, indices, or symbols for which constraints should be ignored
+ if they depend on one of them.
+
+ Returns
+ -------
+ list
+ The bounds for the specified variable.
+ """
+ if isinstance(index_or_variable, int):
+ j = index_or_variable
+ else:
+ j = self.variable_to_index(index_or_variable)
+
+ if ignore_variables:
+ ignore_variables = [
+ i if isinstance(i, int) else self.variable_to_index(i)
+ for i in ignore_variables
+ ]
+
+ bounds = []
+ for i in range(self.A.shape[0]):
+ if self.A[i, j] > 0:
+ if ignore_variables and any(
+ self.A[i, k] != 0 for k in ignore_variables
+ ):
+ # Skip constraint that depends on any of the ignored variables
+ continue
+
+ components = [
+ self._to_literal(self.A[i, k]) * self.variables[k]
+ for k in range(self.A.shape[1])
+ if k != j and self.A[i, k] != 0
+ ]
+ if not components:
+ lhs = sym.IntLiteral(0)
+ elif len(components) == 1:
+ lhs = components[0]
+ else:
+ lhs = sym.Sum(as_tuple(components))
+ bounds += [
+ simplify(
+ sym.Quotient(
+ self._to_literal(self.b[i]) - lhs,
+ self._to_literal(self.A[i, j]),
+ )
+ )
+ ]
+ return bounds
+
+ @staticmethod
+ def generate_entries_for_lower_bound(bound, variables, index):
+ """
+ Helper routine to generate matrix and right-hand side entries for a given lower bound.
+
+ Note that this routine can only handle affine bounds, which means expressions that are
+ constant or can be reduced to a linear polynomial.
+
+ Upper bounds can be derived from this by multiplying the left-hand side and right-hand side
+ with -1.
+
+ Parameters
+ ----------
+ bound : int or str or sym.Array or sym.Scalar
+ The expression representing the lower bound.
+ variables : list of str
+ The list of variable names.
+ index : int
+ The index of the variable constrained by this bound.
+
+ Returns
+ -------
+ tuple(np.array, np.array)
+ The pair ``(lhs, rhs)`` of the row in the matrix inequality, where `lhs` is the left-hand side
+ and `rhs` is the right-hand side.
+ """
+ supported_types = (sym.TypedSymbol, sym.MetaSymbol, sym.Sum, sym.Product)
+ if not (is_constant(bound) or isinstance(bound, supported_types)):
+ raise ValueError(f"Cannot derive inequality from bound {str(bound)}")
+ summands = accumulate_polynomial_terms(bound)
+ b = -summands.pop(1, 0) # Constant term or 0
+ A = np.zeros([1, len(variables)], dtype=np.dtype(int))
+ A[0, index] = -1
+ for base, coef in summands.items():
+ if not len(base) == 1:
+ raise ValueError(f"Non-affine bound {str(bound)}")
+ A[0, variables.index(base[0].name.lower())] = coef
+ return A, b
+
+ @classmethod
+ def from_loop_ranges(cls, loop_variables, loop_ranges):
+ """
+ Create polyhedron from a list of loop ranges and associated variables.
+ """
+ assert len(loop_ranges) == len(loop_variables)
+
+ # Add any variables that are not loop variables to the vector of variables
+ variables = list(loop_variables)
+ variable_names = [v.name.lower() for v in variables]
+ for v in sorted(
+ FindVariables().visit(loop_ranges), key=lambda v: v.name.lower()
+ ):
+ if v.name.lower() not in variable_names:
+ variables += [v]
+ variable_names += [v.name.lower()]
+
+ n = 2 * len(loop_ranges)
+ d = len(variables)
+ A = np.zeros([n, d], dtype=np.dtype(int))
+ b = np.zeros([n], dtype=np.dtype(int))
+
+ for i, (loop_variable, loop_range) in enumerate(
+ zip(loop_variables, loop_ranges)
+ ):
+ assert loop_range.step is None or loop_range.step == "1"
+ j = variables.index(loop_variable.name.lower())
+
+ # Create inequality from lower bound
+ lhs, rhs = cls.generate_entries_for_lower_bound(
+ loop_range.start, variable_names, j
+ )
+ A[2 * i, :] = lhs
+ b[2 * i] = rhs
+
+ # Create inequality from upper bound
+ lhs, rhs = cls.generate_entries_for_lower_bound(
+ loop_range.stop, variable_names, j
+ )
+ A[2 * i + 1, :] = -lhs
+ b[2 * i + 1] = -rhs
+
+ return cls(A, b, variables)
+
+ @classmethod
+ def from_nested_loops(cls, nested_loops: List[Loop]):
+ """
+ Helper function, for creating a polyhedron from a list of loops.
+ """
+ return cls.from_loop_ranges(
+ [l.variable for l in nested_loops], [l.bounds for l in nested_loops]
+ )
diff --git a/loki/transform/transform_loop.py b/loki/transform/transform_loop.py
index f11b2e667..117c7405e 100644
--- a/loki/transform/transform_loop.py
+++ b/loki/transform/transform_loop.py
@@ -15,8 +15,7 @@
from loki.expression import (
symbols as sym, SubstituteExpressions, FindVariables,
- accumulate_polynomial_terms, simplify, is_constant, symbolic_op,
- TypedSymbol
+ simplify, is_constant, symbolic_op,
)
from loki.frontend.fparser import parse_fparser_expression
from loki.ir import Loop, Conditional, Comment, Pragma
@@ -33,216 +32,10 @@
dataflow_analysis_attached, read_after_write_vars, loop_carried_dependencies
)
-__all__ = ['loop_interchange', 'loop_fusion', 'loop_fission', 'Polyhedron']
+__all__ = ['loop_interchange', 'loop_fusion', 'loop_fission']
-class Polyhedron:
- """
- Halfspace representation of a (convex) polyhedron.
-
- A polyhedron `P c R^d` is described by a set of inequalities, in matrix form
- ```
- P = { x=[x1,...,xd]^T c R^d | Ax <= b }
- ```
- with n-by-d matrix `A` and d-dimensional right hand side `b`.
-
- In loop transformations, polyhedrons are used to represent iteration spaces of
- d-deep loop nests.
-
- :param np.array A: the representation matrix A.
- :param np.array b: the right hand-side vector b.
- :param list variables: list of variables representing the dimensions in the polyhedron.
- """
-
- def __init__(self, A, b, variables=None):
- A = np.array(A, dtype=np.dtype(int))
- b = np.array(b, dtype=np.dtype(int))
- assert A.ndim == 2 and b.ndim == 1
- assert A.shape[0] == b.shape[0]
- self.A = A
- self.b = b
-
- self.variables = None
- self.variable_names = None
- if variables is not None:
- assert len(variables) == A.shape[1]
- self.variables = variables
- self.variable_names = [v.name.lower() for v in self.variables]
-
- def variable_to_index(self, variable):
- if self.variable_names is None:
- raise RuntimeError('No variables list associated with polyhedron.')
- if isinstance(variable, TypedSymbol):
- variable = variable.name.lower()
- assert isinstance(variable, str)
- return self.variable_names.index(variable)
-
- @staticmethod
- def _to_literal(value):
- if value < 0:
- return sym.Product((-1, sym.IntLiteral(abs(value))))
- return sym.IntLiteral(value)
-
- def lower_bounds(self, index_or_variable, ignore_variables=None):
- """
- Return all lower bounds imposed on a variable.
-
- Lower bounds for variable `j` are given by the index set
- ```
- L = {i in {0,...,d-1} | A_ij < 0}
- ```
-
- :param index_or_variable: the index, name, or expression symbol for which the
- lower bounds are produced.
- :type index_or_variable: int or str or sym.Array or sym.Scalar
- :param ignore_variables: optional list of variable names, indices or symbols
- for which constraints should be ignored if they depend on one of them.
- :type ignore_variables: list or None
-
- :returns list: the bounds for that variable.
- """
- if isinstance(index_or_variable, int):
- j = index_or_variable
- else:
- j = self.variable_to_index(index_or_variable)
-
- if ignore_variables:
- ignore_variables = [i if isinstance(i, int) else self.variable_to_index(i)
- for i in ignore_variables]
-
- bounds = []
- for i in range(self.A.shape[0]):
- if self.A[i,j] < 0:
- if ignore_variables and any(self.A[i, k] != 0 for k in ignore_variables):
- # Skip constraint that depends on any of the ignored variables
- continue
-
- components = [self._to_literal(self.A[i,k]) * self.variables[k]
- for k in range(self.A.shape[1]) if k != j and self.A[i,k] != 0]
- if not components:
- lhs = sym.IntLiteral(0)
- elif len(components) == 1:
- lhs = components[0]
- else:
- lhs = sym.Sum(as_tuple(components))
- bounds += [simplify(sym.Quotient(self._to_literal(self.b[i]) - lhs,
- self._to_literal(self.A[i,j])))]
- return bounds
-
- def upper_bounds(self, index_or_variable, ignore_variables=None):
- """
- Return all upper bounds imposed on a variable.
-
- Upper bounds for variable `j` are given by the index set
- ```
- U = {i in {0,...,d-1} | A_ij > 0}
- ```
-
- :param index_or_variable: the index, name, or expression symbol for which the
- upper bounds are produced.
- :type index_or_variable: int or str or sym.Array or sym.Scalar
- :param ignore_variables: optional list of variable names, indices or symbols
- for which constraints should be ignored if they depend on one of them.
- :type ignore_variables: list or None
-
- :returns list: the bounds for that variable.
- """
- if isinstance(index_or_variable, int):
- j = index_or_variable
- else:
- j = self.variable_to_index(index_or_variable)
-
- if ignore_variables:
- ignore_variables = [i if isinstance(i, int) else self.variable_to_index(i)
- for i in ignore_variables]
-
- bounds = []
- for i in range(self.A.shape[0]):
- if self.A[i,j] > 0:
- if ignore_variables and any(self.A[i, k] != 0 for k in ignore_variables):
- # Skip constraint that depends on any of the ignored variables
- continue
-
- components = [self._to_literal(self.A[i,k]) * self.variables[k]
- for k in range(self.A.shape[1]) if k != j and self.A[i,k] != 0]
- if not components:
- lhs = sym.IntLiteral(0)
- elif len(components) == 1:
- lhs = components[0]
- else:
- lhs = sym.Sum(as_tuple(components))
- bounds += [simplify(sym.Quotient(self._to_literal(self.b[i]) - lhs,
- self._to_literal(self.A[i,j])))]
- return bounds
-
- @staticmethod
- def generate_entries_for_lower_bound(bound, variables, index):
- """
- Helper routine to generate matrix and right-hand side entries for a
- given lower bound.
-
- NB: This can only deal with affine bounds, i.e. expressions that are
- constant or can be reduced to a linear polynomial.
-
- Upper bounds can be derived from this by multiplying left-hand side and
- right-hand side with -1.
-
- :param bound: the expression representing the lower bound.
- :param list variables: the list of variable names.
- :param int index: the index of the variable constrained by this bound.
-
- :return: the pair ``(lhs, rhs)`` of the row in the matrix inequality.
- :rtype: tuple(np.array, np.array)
- """
- supported_types = (sym.TypedSymbol, sym.MetaSymbol, sym.Sum, sym.Product)
- if not (is_constant(bound) or isinstance(bound, supported_types)):
- raise ValueError(f'Cannot derive inequality from bound {str(bound)}')
- summands = accumulate_polynomial_terms(bound)
- b = -summands.pop(1, 0) # Constant term or 0
- A = np.zeros([1, len(variables)], dtype=np.dtype(int))
- A[0, index] = -1
- for base, coef in summands.items():
- if not len(base) == 1:
- raise ValueError(f'Non-affine bound {str(bound)}')
- A[0, variables.index(base[0].name.lower())] = coef
- return A, b
-
- @classmethod
- def from_loop_ranges(cls, loop_variables, loop_ranges):
- """
- Create polyhedron from a list of loop ranges and associated variables.
- """
- assert len(loop_ranges) == len(loop_variables)
-
- # Add any variables that are not loop variables to the vector of variables
- variables = list(loop_variables)
- variable_names = [v.name.lower() for v in variables]
- for v in sorted(FindVariables().visit(loop_ranges), key=lambda v: v.name.lower()):
- if v.name.lower() not in variable_names:
- variables += [v]
- variable_names += [v.name.lower()]
-
- n = 2 * len(loop_ranges)
- d = len(variables)
- A = np.zeros([n, d], dtype=np.dtype(int))
- b = np.zeros([n], dtype=np.dtype(int))
-
- for i, (loop_variable, loop_range) in enumerate(zip(loop_variables, loop_ranges)):
- assert loop_range.step is None or loop_range.step == '1'
- j = variables.index(loop_variable.name.lower())
-
- # Create inequality from lower bound
- lhs, rhs = cls.generate_entries_for_lower_bound(loop_range.start, variable_names, j)
- A[2*i,:] = lhs
- b[2*i] = rhs
-
- # Create inequality from upper bound
- lhs, rhs = cls.generate_entries_for_lower_bound(loop_range.stop, variable_names, j)
- A[2*i+1,:] = -lhs
- b[2*i+1] = -rhs
-
- return cls(A, b, variables)
-
+from loki.analyse.util_polyhedron import Polyhedron
def eliminate_variable(polyhedron, index_or_variable):
"""
diff --git a/pyproject.toml b/pyproject.toml
index f451c5d0b..559a3fb93 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -64,6 +64,10 @@ examples = [
"ipyparams",
]
+dataflow_analysis = [
+ "ortools"
+]
+
[project.scripts]
"loki-transform.py" = "scripts.loki_transform:cli"
"loki-lint.py" = "scripts.loki_lint:cli"
diff --git a/tests/sources/data_dependency_detection/loop_carried_dependencies.f90 b/tests/sources/data_dependency_detection/loop_carried_dependencies.f90
new file mode 100644
index 000000000..b208a6157
--- /dev/null
+++ b/tests/sources/data_dependency_detection/loop_carried_dependencies.f90
@@ -0,0 +1,58 @@
+subroutine SimpleDependency(data, n)
+ implicit none
+ integer, intent(in) :: n
+ real(8), dimension(n) :: data
+ integer :: i
+
+ ! Loop with a simple loop-carried dependency
+ do i = 1, n
+ data(i) = data(i) + data(i-1)
+ end do
+
+end subroutine SimpleDependency
+
+
+subroutine NestedDependency(data, n)
+ implicit none
+ integer, intent(in) :: n
+ real(8), dimension(n) :: data
+ integer :: i, j
+
+ ! Nested loop with a loop-carried dependency
+ do i = 2, n
+ do j = 1, i-1
+ data(i) = data(i) + data(j)
+ end do
+ end do
+
+end subroutine NestedDependency
+
+
+subroutine ConditionalDependency(data, n)
+ implicit none
+ integer, intent(in) :: n
+ real(8), dimension(n) :: data
+ integer :: i
+
+ ! Loop with a conditional loop-carried dependency
+ do i = 2, n
+ if (data(i-1) > 0.0) then
+ data(i) = data(i) + data(i-1)
+ endif
+ end do
+
+end subroutine ConditionalDependency
+
+subroutine NoDependency(data)
+ implicit none
+ real(8), dimension(20) :: data
+ integer :: i
+
+ do i = 1, 10, 1
+ data(2*i) = 10;
+ end do
+
+ do i = 1, 5, 1
+ data(2*i + 1) = 20;
+ end do
+end subroutine NoDependency
\ No newline at end of file
diff --git a/tests/sources/data_dependency_detection/various_loops.f90 b/tests/sources/data_dependency_detection/various_loops.f90
new file mode 100644
index 000000000..514d7e567
--- /dev/null
+++ b/tests/sources/data_dependency_detection/various_loops.f90
@@ -0,0 +1,79 @@
+subroutine single_loop(arr, n)
+ implicit none
+ integer, intent(inout) :: arr(:)
+ integer, intent(in) :: n
+ integer :: i
+
+ do i = 1, n
+ arr(i) = arr(i) * 2
+ end do
+end subroutine single_loop
+
+subroutine single_loop_split_access(arr, n)
+ implicit none
+ integer, intent(inout) :: arr(:)
+ integer, intent(in) :: n
+ integer :: i, nhalf
+ nhalf = n/2
+ do i = 1, nhalf
+ arr(2*i) = arr(2*i) * 2
+ arr(2*i + 1) = arr(2*i + 1) * 2
+ end do
+end subroutine single_loop_split_access
+
+subroutine single_loop_arithmetic_operations_for_access(arr, n)
+ implicit none
+ integer, intent(inout) :: arr(:)
+ integer, intent(in) :: n
+ integer :: i
+
+ do i = 1, n
+ arr(i*i) = arr(i + i) * 2
+ end do
+end subroutine single_loop_arithmetic_operations_for_access
+
+subroutine nested_loop_single_dimensions_access(arr, n)
+ implicit none
+ integer, intent(inout) :: arr(:)
+ integer, intent(in) :: n
+ integer :: i, j, nhalf
+
+ nhalf = n/2
+ do i = 1, nhalf
+ do j = 1, nhalf
+ arr(i + j) = arr(i + j) * 2
+ end do
+ end do
+end subroutine nested_loop_single_dimensions_access
+
+
+
+subroutine nested_loop_partially_used(arr, n)
+ implicit none
+ integer, intent(inout) :: arr(:)
+ integer, intent(in) :: n
+ integer :: i, j, nfourth
+
+ nfourth = n / 4
+ do i = 1, nfourth
+ do j = 1, nfourth
+ arr(i + j) = arr(i + j) * 2
+ end do
+ end do
+end subroutine nested_loop_partially_used
+
+
+subroutine partially_used_array(arr, n)
+ implicit none
+ integer, intent(out) :: arr(:)
+ integer, intent(in) :: n
+ integer :: i = 1 , j = 3, nhalf
+
+ nhalf = n/2
+ arr(1) = 1
+ do i = 2, nhalf
+ arr(i) = arr(i - 1)
+ end do
+
+ j = arr(j)
+end subroutine partially_used_array
\ No newline at end of file
diff --git a/tests/test_analyse_array_data_dependency_detection.py b/tests/test_analyse_array_data_dependency_detection.py
new file mode 100644
index 000000000..d90bbe06e
--- /dev/null
+++ b/tests/test_analyse_array_data_dependency_detection.py
@@ -0,0 +1,403 @@
+# (C) Copyright 2018- ECMWF.
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+
+import pytest
+import numpy as np
+from loki import Scope, parse_fparser_expression
+from loki.analyse.analyse_array_data_dependency_detection import (
+ has_data_dependency,
+ HAVE_ORTOOLS,
+ construct_affine_array_access_function_representation,
+)
+
+array = np.array
+
+@pytest.mark.parametrize(
+ "first_access_represenetation, second_access_represenetation, error_type",
+ [
+ # totaly empty tuple
+ (tuple(), tuple(), ValueError),
+ # tuple with two empty tuple --> correct form
+ ((tuple(), tuple()), (tuple(), tuple()), ValueError),
+ #
+ ],
+)
+def test_has_data_dependency_false_inputs(
+ first_access_represenetation, second_access_represenetation, error_type
+):
+ with pytest.raises(error_type):
+ _ = has_data_dependency(
+ first_access_represenetation, second_access_represenetation
+ )
+
+@pytest.mark.parametrize(
+ "array_dimensions_expr, expected",
+ [
+ ("i-1", ([[1, 0]], [[-1]])),
+ ("i,j", ([[1, 0], [0, 1]], [[0], [0]])),
+ ("j,j+1", ([[0, 1], [0, 1]], [[0], [1]])),
+ ("1,2", ([[0, 0], [0, 0]], [[1], [2]])),
+ ("1,i,2*i+j", ([[0, 0], [1, 0], [2, 1]], [[1], [0], [0]])),
+ ],
+)
+def test_access_function_creation(array_dimensions_expr, expected):
+ scope = Scope()
+ first = parse_fparser_expression(f"z({array_dimensions_expr})", scope)
+
+ use_these_variables = ["i", "j"]
+
+ access, variables = construct_affine_array_access_function_representation(
+ first.dimensions, use_these_variables
+ )
+
+ assert np.array_equal(access.matrix, np.array(expected[0], dtype=np.dtype(int)))
+ assert np.array_equal(access.vector, np.array(expected[1], dtype=np.dtype(int)))
+ assert np.array_equal(variables, np.array(["i", "j"], dtype=np.dtype(object)))
+
+expected_result = {
+ "Example 11.29.1. from Compilers: Principles, Techniques, and Tools": True,
+ "Example 11.30 from Compilers: Principles, Techniques, and Tools": False,
+ "Example 11.35 from Compilers: Principles, Techniques, and Tools": False,
+ "Example Anti Dependency": True,
+ "Exercise 11.6.5 a) from Compilers: Principles, Techniques, and Tools": not HAVE_ORTOOLS,
+ "Exercise 11.6.5 b) from Compilers: Principles, Techniques, and Tools": not HAVE_ORTOOLS,
+ "Exercise 11.6.5 c) from Compilers: Principles, Techniques, and Tools": True,
+ "Exercise 11.6.5 d) from Compilers: Principles, Techniques, and Tools": not HAVE_ORTOOLS,
+}
+
+
+examples = [
+ (
+ "Example 11.29.1. from Compilers: Principles, Techniques, and Tools",
+ # Comparing Z[i] and Z[i-1]
+ # for (i = 1; i <= 10; i++) {
+ # Z[i] = Z[i-1];
+ # }
+ #
+ # first access
+ (
+ # iteration space representation
+ (array([[1], [-1]]), array([[10], [-1]])),
+ # function access representation
+ (array([[1]]), array([[-1]])),
+ ),
+ # second access
+ (
+ # iteration space representation
+ (array([[1], [-1]]), array([[10], [-1]])),
+ # function access representation
+ (array([[1]]), array([[0]])),
+ ),
+ ),
+ # Example 11.29.2. from Compilers: Principles, Techniques, and Tools
+ # is ignored, since it shows a comparison between an access with itself which does not make any sense
+ #
+ (
+ "Example 11.30 from Compilers: Principles, Techniques, and Tools",
+ # Comparing Z[2*i] and Z[2*i+1]
+ # for (i = 1; i < 10; i++) {
+ # Z[2*i] = 10;
+ # }
+ # for (j = 1; j < 10; j++) {
+ # Z[2*j+1] = 20;
+ # }
+ #
+ # first access
+ (
+ # iteration space representation
+ (array([[1, 0], [-1, 0]]), array([[10], [-1]])),
+ # function access representation
+ (array([[2, 0]]), array([[0]])),
+ ),
+ # second access
+ (
+ # iteration space representation
+ (array([[0, 1], [0, -1]]), array([[10], [-1]])),
+ # function access representation
+ (array([[0, 2]]), array([[1]])),
+ ),
+ ),
+ (
+ "Example 11.35 from Compilers: Principles, Techniques, and Tools",
+ # Comparing Z[i,j] and Z[j+10,i+11]
+ # for (i = 0; i <= 10; i++)
+ # for (j = 0; j <= 10; j++)
+ # Z[i,j] = Z[j+10,i+11];
+ #
+ # first access
+ (
+ # iteration space representation
+ (
+ array([[-1, 0], [0, -1], [1, 0], [0, 1]]),
+ array([[0], [0], [10], [10]]),
+ ),
+ # function access representation
+ (array([[1, 0], [0, 1]]), array([[0], [0]])),
+ ),
+ # second access
+ (
+ # iteration space representation
+ (
+ array([[-1, 0], [0, -1], [1, 0], [0, 1]]),
+ array([[0], [0], [10], [10]]),
+ ),
+ # function access representation
+ (array([[0, 1], [1, 0]]), array([[10], [11]])),
+ ),
+ ),
+ (
+ "Example Anti Dependency",
+ # Comparing Z[i,j] and Z[j+10,i+11]
+ # for (i = 2; i <= N; i++)
+ # Z[i] = Z[i-1] + 1;
+ #
+ # first access
+ (
+ # iteration space representation
+ (
+ array([[1, -1], [-1, 0]]),
+ array([[0], [-2]]),
+ ),
+ # function access representation
+ (array([[1, 0]]), array([[0]])),
+ ),
+ # second access
+ (
+ # iteration space representation
+ (
+ array([[1, -1], [-1, 0]]),
+ array([[0], [-2]]),
+ ),
+ # function access representation
+ (array([[1, 0]]), array([[-1]])),
+ ),
+ ),
+ (
+ "Exercise 11.6.5 a) from Compilers: Principles, Techniques, and Tools",
+ # Comparing Z[i,j,k] = Z[i+100,j+100,k+100];
+ # for (i=0; i<100; i++)
+ # for (j=0; j<100; j++)
+ # for (k=0; k<100; k++)
+ # Z[i,j,k] = Z[i+100,j+100,k+100];
+ #
+ # first access
+ (
+ # iteration space representation
+ (
+ array(
+ [
+ [1, 0, 0],
+ [0, 1, 0],
+ [0, 0, 1],
+ [-1, 0, 0],
+ [0, -1, 0],
+ [0, 0, -1],
+ ]
+ ),
+ array([[99], [99], [99], [0], [0], [0]]),
+ ),
+ # function access representation
+ (
+ array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]),
+ array([[0], [0], [0]]),
+ ),
+ ),
+ # second access
+ (
+ # iteration space representation
+ (
+ array(
+ [
+ [1, 0, 0],
+ [0, 1, 0],
+ [0, 0, 1],
+ [-1, 0, 0],
+ [0, -1, 0],
+ [0, 0, -1],
+ ]
+ ),
+ array([[99], [99], [99], [0], [0], [0]]),
+ ),
+ # function access representation
+ (
+ array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]),
+ array([[100], [100], [100]]),
+ ),
+ ),
+ ),
+ (
+ "Exercise 11.6.5 b) from Compilers: Principles, Techniques, and Tools",
+ # Comparing Z[i,j,k] = Z[j+100,k+100,i+100];
+ # for (i=0; i<100; i++)
+ # for (j=0; j<100; j++)
+ # for (k=0; k<100; k++)
+ # Z[i,j,k] = Z[j+100,k+100,i+100];
+ #
+ # first access
+ (
+ # iteration space representation
+ (
+ array(
+ [
+ [1, 0, 0],
+ [0, 1, 0],
+ [0, 0, 1],
+ [-1, 0, 0],
+ [0, -1, 0],
+ [0, 0, -1],
+ ]
+ ),
+ array([[99], [99], [99], [0], [0], [0]]),
+ ),
+ # function access representation
+ (
+ array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]),
+ array([[0], [0], [0]]),
+ ),
+ ),
+ # second access
+ (
+ # iteration space representation
+ (
+ array(
+ [
+ [1, 0, 0],
+ [0, 1, 0],
+ [0, 0, 1],
+ [-1, 0, 0],
+ [0, -1, 0],
+ [0, 0, -1],
+ ]
+ ),
+ array([[99], [99], [99], [0], [0], [0]]),
+ ),
+ # function access representation
+ (
+ array([[0, 1, 0], [0, 0, 1], [1, 0, 0]]),
+ array([[100], [100], [100]]),
+ ),
+ ),
+ ),
+ (
+ "Exercise 11.6.5 c) from Compilers: Principles, Techniques, and Tools",
+ # Comparing Z[i,j,k] = Z[j-50,k-50,i-50];
+ # for (i=0; i<100; i++)
+ # for (j=0; j<100; j++)
+ # for (k=0; k<100; k++)
+ # Z[i,j,k] = Z[j-50,k-50,i-50];
+ #
+ # first access
+ (
+ # iteration space representation
+ (
+ array(
+ [
+ [1, 0, 0],
+ [0, 1, 0],
+ [0, 0, 1],
+ [-1, 0, 0],
+ [0, -1, 0],
+ [0, 0, -1],
+ ]
+ ),
+ array([[99], [99], [99], [0], [0], [0]]),
+ ),
+ # function access representation
+ (
+ array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]),
+ array([[0], [0], [0]]),
+ ),
+ ),
+ # second access
+ (
+ # iteration space representation
+ (
+ array(
+ [
+ [1, 0, 0],
+ [0, 1, 0],
+ [0, 0, 1],
+ [-1, 0, 0],
+ [0, -1, 0],
+ [0, 0, -1],
+ ]
+ ),
+ array([[99], [99], [99], [0], [0], [0]]),
+ ),
+ # function access representation
+ (
+ array([[0, 1, 0], [0, 0, 1], [1, 0, 0]]),
+ array([[-50], [-50], [-50]]),
+ ),
+ ),
+ ),
+ (
+ "Exercise 11.6.5 d) from Compilers: Principles, Techniques, and Tools",
+ # Comparing Z[i,j,k] = Z[i+99,k+100,j];
+ # for (i=0; i<100; i++)
+ # for (j=0; j<100; j++)
+ # for (k=0; k<100; k++)
+ # Z[i,j,k] = Z[i+99,k+100,j];
+ #
+ # first access
+ (
+ # iteration space representation
+ (
+ array(
+ [
+ [1, 0, 0],
+ [0, 1, 0],
+ [0, 0, 1],
+ [-1, 0, 0],
+ [0, -1, 0],
+ [0, 0, -1],
+ ]
+ ),
+ array([[99], [99], [99], [0], [0], [0]]),
+ ),
+ # function access representation
+ (
+ array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]),
+ array([[0], [0], [0]]),
+ ),
+ ),
+ # second access
+ (
+ # iteration space representation
+ (
+ array(
+ [
+ [1, 0, 0],
+ [0, 1, 0],
+ [0, 0, 1],
+ [-1, 0, 0],
+ [0, -1, 0],
+ [0, 0, -1],
+ ]
+ ),
+ array([[99], [99], [99], [0], [0], [0]]),
+ ),
+ # function access representation
+ (
+ array([[1, 0, 0], [0, 0, 1], [0, 1, 0]]),
+ array([[99], [100], [0]]),
+ ),
+ ),
+ ),
+]
+
+
+@pytest.mark.parametrize(
+ "name, first_access_represenetation, second_access_represenetation",
+ examples,
+)
+def test_has_data_dependency(
+ name, first_access_represenetation, second_access_represenetation
+):
+ assert expected_result[name] == has_data_dependency(
+ first_access_represenetation, second_access_represenetation
+ ), f"Test {name} failed!"
diff --git a/tests/test_transform_loop.py b/tests/test_transform_loop.py
index b498558d7..77afc6f38 100644
--- a/tests/test_transform_loop.py
+++ b/tests/test_transform_loop.py
@@ -11,9 +11,9 @@
import numpy as np
from conftest import jit_compile, clean_test, available_frontends
-from loki import Subroutine, OMNI, FindNodes, Loop, Conditional, Scope, Assignment
-from loki.frontend.fparser import parse_fparser_expression, HAVE_FP
-from loki.transform import loop_interchange, loop_fusion, loop_fission, Polyhedron, normalize_range_indexing
+from loki import Subroutine, OMNI, FindNodes, Loop, Conditional, Assignment
+from loki.frontend.fparser import HAVE_FP
+from loki.transform import loop_interchange, loop_fusion, loop_fission, normalize_range_indexing
from loki.expression import symbols as sym
from loki.pragma_utils import is_loki_pragma, pragmas_attached
@@ -26,90 +26,6 @@
def fixture_here():
return Path(__file__).parent
-
-@pytest.mark.parametrize('variables, lbounds, ubounds, A, b, variable_names', [
- # do i=0,5: do j=i,7: ...
- (['i', 'j'], ['0', 'i'], ['5', '7'],
- [[-1, 0], [1, 0], [1, -1], [0, 1]], [0, 5, 0, 7], ['i', 'j']),
- # do i=1,n: do j=0,2*i+1: do k=a,b: ...
- (['i', 'j', 'k'], ['1', '0', 'a'], ['n', '2*i+1', 'b'],
- [[-1, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, -1], [0, -1, 0, 0, 0, 0], [-2, 1, 0, 0, 0, 0],
- [0, 0, -1, 1, 0, 0], [0, 0, 1, 0, -1, 0]], [-1, 0, 0, 1, 0, 0], ['i', 'j', 'k', 'a', 'b', 'n']),
- # do jk=1,klev: ...
- (['jk'], ['1'], ['klev'], [[-1, 0], [1, -1]], [-1, 0], ['jk', 'klev']),
- # do JK=1,klev-1: ...
- (['JK'], ['1'], ['klev - 1'], [[-1, 0], [1, -1]], [-1, -1], ['jk', 'klev']),
- # do jk=ncldtop,klev: ...
- (['jk'], ['ncldtop'], ['klev'], [[-1, 0, 1], [1, -1, 0]], [0, 0], ['jk', 'klev', 'ncldtop']),
- # do jk=1,KLEV+1: ...
- (['jk'], ['1'], ['KLEV+1'], [[-1, 0], [1, -1]], [-1, 1], ['jk', 'klev']),
-])
-def test_polyhedron_from_loop_ranges(variables, lbounds, ubounds, A, b, variable_names):
- """
- Test converting loop ranges to polyedron representation of iteration space.
- """
- scope = Scope()
- loop_variables = [parse_fparser_expression(expr, scope) for expr in variables]
- loop_lbounds = [parse_fparser_expression(expr, scope) for expr in lbounds]
- loop_ubounds = [parse_fparser_expression(expr, scope) for expr in ubounds]
- loop_ranges = [sym.LoopRange((l, u)) for l, u in zip(loop_lbounds, loop_ubounds)]
- p = Polyhedron.from_loop_ranges(loop_variables, loop_ranges)
- assert np.all(p.A == np.array(A, dtype=np.dtype(int)))
- assert np.all(p.b == np.array(b, dtype=np.dtype(int)))
- assert p.variables == variable_names
-
-
-def test_polyhedron_from_loop_ranges_failures():
- """
- Test known limitation of the conversion from loop ranges to polyhedron.
- """
- # m*n is non-affine and thus can't be represented
- scope = Scope()
- loop_variable = parse_fparser_expression('i', scope)
- lower_bound = parse_fparser_expression('1', scope)
- upper_bound = parse_fparser_expression('m * n', scope)
- loop_range = sym.LoopRange((lower_bound, upper_bound))
- with pytest.raises(ValueError):
- _ = Polyhedron.from_loop_ranges([loop_variable], [loop_range])
-
- # no functionality to flatten exponentials, yet
- upper_bound = parse_fparser_expression('5**2', scope)
- loop_range = sym.LoopRange((lower_bound, upper_bound))
- with pytest.raises(ValueError):
- _ = Polyhedron.from_loop_ranges([loop_variable], [loop_range])
-
-
-@pytest.mark.parametrize('A, b, variable_names, lower_bounds, upper_bounds', [
- # do i=1,n: ...
- ([[-1, 0], [1, -1]], [-1, 0], ['i', 'n'], [['1'], ['i']], [['n'], []]),
- # do i=1,10: ...
- ([[-1], [1]], [-1, 10], ['i'], [['1']], [['10']]),
- # do i=0,5: do j=i,7: ...
- ([[-1, 0], [1, 0], [1, -1], [0, 1]], [0, 5, 0, 7], ['i', 'j'], [['0'], ['i']], [['5', 'j'], ['7']]),
- # do i=1,n: do j=0,2*i+1: do k=a,b: ...
- ([[-1, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, -1], [0, -1, 0, 0, 0, 0], [-2, 1, 0, 0, 0, 0],
- [0, 0, -1, 1, 0, 0], [0, 0, 1, 0, -1, 0]], [-1, 0, 0, 1, 0, 0],
- ['i', 'j', 'k', 'a', 'b', 'n'], # variable names
- [['1', '-1 / 2 + j / 2'], ['0'], ['a'], [], ['k'], ['i']], # lower bounds
- [['n'], ['1 + 2*i'], ['b'], ['k'], [], []]), # upper bounds
-])
-def test_polyhedron_bounds(A, b, variable_names, lower_bounds, upper_bounds):
- """
- Test the production of lower and upper bounds.
- """
- scope = Scope()
- variables = [parse_fparser_expression(v, scope) for v in variable_names]
- p = Polyhedron(A, b, variables)
- for var, ref_bounds in zip(variables, lower_bounds):
- lbounds = p.lower_bounds(var)
- assert len(lbounds) == len(ref_bounds)
- assert all(str(b1) == b2 for b1, b2 in zip(lbounds, ref_bounds))
- for var, ref_bounds in zip(variables, upper_bounds):
- ubounds = p.upper_bounds(var)
- assert len(ubounds) == len(ref_bounds)
- assert all(str(b1) == b2 for b1, b2 in zip(ubounds, ref_bounds))
-
-
@pytest.mark.parametrize('frontend', available_frontends())
def test_transform_loop_interchange_plain(here, frontend):
"""
diff --git a/tests/test_util_linear_algebra.py b/tests/test_util_linear_algebra.py
index fc93d4ee5..caf45dda4 100644
--- a/tests/test_util_linear_algebra.py
+++ b/tests/test_util_linear_algebra.py
@@ -4,10 +4,18 @@
# In applying this licence, ECMWF does not waive the privileges and immunities
# granted to it by virtue of its status as an intergovernmental organisation
# nor does it submit to any jurisdiction.
-from math import gcd
+from math import gcd as math_gcd
import pytest
import numpy as np
+try:
+ _ = math_gcd(4,3,2)
+ gcd = math_gcd
+except TypeError: #Python 3.8 can only handle two arguments
+ from functools import reduce
+ def gcd(*args):
+ return reduce(math_gcd, args)
+
from loki.analyse.util_linear_algebra import (
back_substitution,
generate_row_echelon_form,
@@ -149,30 +157,32 @@ def test_is_independent_system(matrix, expected_result):
(
np.array([[1, 0], [0, 1], [0, 0]]),
np.array([[1], [2], [0]]),
- [[1], [1]],
- [[1], [2]],
+ [np.array([[0]]), np.array([[0]]), np.array([[1]]), np.array([[1]])],
+ [np.array([[0]]), np.array([[0]]), np.array([[1]]), np.array([[2]])],
),
(
np.array([[1, 0], [0, 1], [0, 0]]),
np.array([[1], [2], [1]]),
- [[0], [1], [1]],
- [[1], [1], [2]],
+ [np.array([[0]]), np.array([[0]]), np.array([[1]]), np.array([[1]])],
+ [np.array([[1]]), np.array([[1]]), np.array([[1]]), np.array([[2]])],
),
(
np.array([[0, 0, 0], [0, 0, 0], [0, 0, 0]]),
np.array([[1], [2], [3]]),
- [[0]],
- [[1, 2, 3]],
+ [np.array([[0], [0], [0]])] * 3,
+ [np.array([[1], [2], [3]])] * 3,
),
( # will even split non independent systems, call is_independent_system before
np.array([[2, 1], [1, 3]]),
np.array([[3], [4]]),
- [[2], [1]],
- [[3], [4]],
+ [np.array([[2], [1]]), np.array([[1], [3]])],
+ [np.array([[3], [4]]), np.array([[3], [4]])],
),
],
)
def test_yield_one_d_systems(matrix, rhs, list_of_lhs_column, list_of_rhs_column):
- for index, (A, b) in enumerate(yield_one_d_systems(matrix, rhs)):
- assert np.allclose(A, list_of_lhs_column[index])
- assert np.allclose(b, list_of_rhs_column[index])
+ results = list(yield_one_d_systems(matrix, rhs))
+ assert len(results) == len(list_of_lhs_column) == len(list_of_rhs_column)
+ for index, (A, b) in enumerate(results):
+ assert np.array_equal(A, list_of_lhs_column[index])
+ assert np.array_equal(b, list_of_rhs_column[index])
diff --git a/tests/test_util_polyhedron.py b/tests/test_util_polyhedron.py
new file mode 100644
index 000000000..db4d716f0
--- /dev/null
+++ b/tests/test_util_polyhedron.py
@@ -0,0 +1,294 @@
+# (C) Copyright 2018- ECMWF.
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+
+from pathlib import Path
+import pytest
+import numpy as np
+
+from loki.ir import Scope, Loop
+from loki.sourcefile import Sourcefile
+from loki.frontend.fparser import parse_fparser_expression, HAVE_FP
+from loki.analyse.util_polyhedron import Polyhedron
+from loki.expression import symbols as sym
+from loki.visitors import FindNodes
+
+
+@pytest.fixture(scope="module", name="here")
+def fixture_here():
+ return Path(__file__).parent
+
+
+# Polyhedron functionality relies on FParser's expression parsing
+pytestmark = pytest.mark.skipif(not HAVE_FP, reason="Fparser not available")
+
+
+@pytest.mark.parametrize(
+ "variables, lbounds, ubounds, A, b, variable_names",
+ [
+ # do i=0,5: do j=i,7: ...
+ (
+ ["i", "j"],
+ ["0", "i"],
+ ["5", "7"],
+ [[-1, 0], [1, 0], [1, -1], [0, 1]],
+ [0, 5, 0, 7],
+ ["i", "j"],
+ ),
+ # do i=1,n: do j=0,2*i+1: do k=a,b: ...
+ (
+ ["i", "j", "k"],
+ ["1", "0", "a"],
+ ["n", "2*i+1", "b"],
+ [
+ [-1, 0, 0, 0, 0, 0],
+ [1, 0, 0, 0, 0, -1],
+ [0, -1, 0, 0, 0, 0],
+ [-2, 1, 0, 0, 0, 0],
+ [0, 0, -1, 1, 0, 0],
+ [0, 0, 1, 0, -1, 0],
+ ],
+ [-1, 0, 0, 1, 0, 0],
+ ["i", "j", "k", "a", "b", "n"],
+ ),
+ # do jk=1,klev: ...
+ (["jk"], ["1"], ["klev"], [[-1, 0], [1, -1]], [-1, 0], ["jk", "klev"]),
+ # do JK=1,klev-1: ...
+ (["JK"], ["1"], ["klev - 1"], [[-1, 0], [1, -1]], [-1, -1], ["jk", "klev"]),
+ # do jk=ncldtop,klev: ...
+ (
+ ["jk"],
+ ["ncldtop"],
+ ["klev"],
+ [[-1, 0, 1], [1, -1, 0]],
+ [0, 0],
+ ["jk", "klev", "ncldtop"],
+ ),
+ # do jk=1,KLEV+1: ...
+ (["jk"], ["1"], ["KLEV+1"], [[-1, 0], [1, -1]], [-1, 1], ["jk", "klev"]),
+ ],
+)
+def test_polyhedron_from_loop_ranges(variables, lbounds, ubounds, A, b, variable_names):
+ """
+ Test converting loop ranges to polyedron representation of iteration space.
+ """
+ scope = Scope()
+ loop_variables = [parse_fparser_expression(expr, scope) for expr in variables]
+ loop_lbounds = [parse_fparser_expression(expr, scope) for expr in lbounds]
+ loop_ubounds = [parse_fparser_expression(expr, scope) for expr in ubounds]
+ loop_ranges = [sym.LoopRange((l, u)) for l, u in zip(loop_lbounds, loop_ubounds)]
+ p = Polyhedron.from_loop_ranges(loop_variables, loop_ranges)
+ assert np.all(p.A == np.array(A, dtype=np.dtype(int)))
+ assert np.all(p.b == np.array(b, dtype=np.dtype(int)))
+ assert p.variables == variable_names
+
+
+def test_polyhedron_from_loop_ranges_failures():
+ """
+ Test known limitation of the conversion from loop ranges to polyhedron.
+ """
+ # m*n is non-affine and thus can't be represented
+ scope = Scope()
+ loop_variable = parse_fparser_expression("i", scope)
+ lower_bound = parse_fparser_expression("1", scope)
+ upper_bound = parse_fparser_expression("m * n", scope)
+ loop_range = sym.LoopRange((lower_bound, upper_bound))
+ with pytest.raises(ValueError):
+ _ = Polyhedron.from_loop_ranges([loop_variable], [loop_range])
+
+ # no functionality to flatten exponentials, yet
+ upper_bound = parse_fparser_expression("5**2", scope)
+ loop_range = sym.LoopRange((lower_bound, upper_bound))
+ with pytest.raises(ValueError):
+ _ = Polyhedron.from_loop_ranges([loop_variable], [loop_range])
+
+
+@pytest.mark.parametrize(
+ "A, b, variable_names, lower_bounds, upper_bounds",
+ [
+ # do i=1,n: ...
+ ([[-1, 0], [1, -1]], [-1, 0], ["i", "n"], [["1"], ["i"]], [["n"], []]),
+ # do i=1,10: ...
+ ([[-1], [1]], [-1, 10], ["i"], [["1"]], [["10"]]),
+ # do i=0,5: do j=i,7: ...
+ (
+ [[-1, 0], [1, 0], [1, -1], [0, 1]],
+ [0, 5, 0, 7],
+ ["i", "j"],
+ [["0"], ["i"]],
+ [["5", "j"], ["7"]],
+ ),
+ # do i=1,n: do j=0,2*i+1: do k=a,b: ...
+ (
+ [
+ [-1, 0, 0, 0, 0, 0],
+ [1, 0, 0, 0, 0, -1],
+ [0, -1, 0, 0, 0, 0],
+ [-2, 1, 0, 0, 0, 0],
+ [0, 0, -1, 1, 0, 0],
+ [0, 0, 1, 0, -1, 0],
+ ],
+ [-1, 0, 0, 1, 0, 0],
+ ["i", "j", "k", "a", "b", "n"], # variable names
+ [["1", "-1 / 2 + j / 2"], ["0"], ["a"], [], ["k"], ["i"]], # lower bounds
+ [["n"], ["1 + 2*i"], ["b"], ["k"], [], []],
+ ), # upper bounds
+ ],
+)
+def test_polyhedron_bounds(A, b, variable_names, lower_bounds, upper_bounds):
+ """
+ Test the production of lower and upper bounds.
+ """
+ scope = Scope()
+ variables = [parse_fparser_expression(v, scope) for v in variable_names]
+ p = Polyhedron(A, b, variables)
+ for var, ref_bounds in zip(variables, lower_bounds):
+ lbounds = p.lower_bounds(var)
+ assert len(lbounds) == len(ref_bounds)
+ assert all(str(b1) == b2 for b1, b2 in zip(lbounds, ref_bounds))
+ for var, ref_bounds in zip(variables, upper_bounds):
+ ubounds = p.upper_bounds(var)
+ assert len(ubounds) == len(ref_bounds)
+ assert all(str(b1) == b2 for b1, b2 in zip(ubounds, ref_bounds))
+
+
+@pytest.mark.parametrize(
+ "polyhedron,is_empty,will_fail",
+ [
+ # totaly empty polyhedron
+ (Polyhedron.from_nested_loops([]), True, False),
+ # full matrix --> non trivial problem
+ (Polyhedron([[1]], [1]), None, True),
+ # empty matrix, full and fullfiled b --> non empty polyhedron
+ (Polyhedron([[]], [1]), False, False),
+ # empty matrix, full b but not fullfiled b --> empty polyhedron
+ (Polyhedron([[]], [-1]), True, False),
+ ],
+)
+def test_check_empty_polyhedron(polyhedron, is_empty, will_fail):
+ if will_fail:
+ with pytest.raises(RuntimeError):
+ _ = polyhedron.is_empty()
+ else:
+ assert polyhedron.is_empty() == is_empty
+
+
+def simple_loop_extractor(start_node):
+ """Find all loops in the AST and structure them depending on their nesting level"""
+ start_loops = FindNodes(Loop, greedy=True).visit(start_node)
+ return [FindNodes(Loop).visit(node) for node in start_loops]
+
+
+def assert_equal_polyhedron(poly_A, poly_B):
+ assert poly_A.variables == poly_B.variables
+ assert (poly_A.A == poly_B.A).all()
+ assert (poly_A.b == poly_B.b).all()
+
+
+@pytest.mark.parametrize(
+ "filename, loop_extractor, polyhedrons_per_subroutine",
+ [
+ (
+ "sources/data_dependency_detection/loop_carried_dependencies.f90",
+ simple_loop_extractor,
+ {
+ "SimpleDependency": [
+ Polyhedron(
+ [[-1, 0], [1, -1]], [-1, 0], [sym.Scalar("i"), sym.Scalar("n")]
+ ),
+ ],
+ "NestedDependency": [
+ Polyhedron(
+ [[-1, 0, 0], [1, 0, -1], [0, -1, 0], [-1, 1, 0]],
+ [-2, 0, -1, -1],
+ [sym.Scalar("i"), sym.Scalar("j"), sym.Scalar("n")],
+ ),
+ ],
+ "ConditionalDependency": [
+ Polyhedron(
+ [[-1, 0], [1, -1]],
+ [-2, 0],
+ [sym.Scalar("i"), sym.Scalar("n")],
+ ),
+ ],
+ "NoDependency": [
+ Polyhedron(
+ [[-1], [1]],
+ [-1, 10],
+ [sym.Scalar("i")],
+ ),
+ Polyhedron(
+ [[-1], [1]],
+ [-1, 5],
+ [sym.Scalar("i")],
+ ),
+ ],
+ },
+ ),
+ (
+ "sources/data_dependency_detection/various_loops.f90",
+ simple_loop_extractor,
+ {
+ "single_loop": [
+ Polyhedron(
+ [[-1, 0], [1, -1]],
+ [-1, 0],
+ [sym.Scalar("i"), sym.Scalar("n")],
+ ),
+ ],
+ "single_loop_split_access": [
+ Polyhedron(
+ [[-1, 0], [1, -1]],
+ [-1, 0],
+ [sym.Scalar("i"), sym.Scalar("nhalf")],
+ ),
+ ],
+ "single_loop_arithmetic_operations_for_access": [
+ Polyhedron(
+ [[-1, 0], [1, -1]],
+ [-1, 0],
+ [sym.Scalar("i"), sym.Scalar("n")],
+ ),
+ ],
+ "nested_loop_single_dimensions_access": [
+ Polyhedron(
+ [[-1, 0, 0], [1, 0, -1], [0, -1, 0], [0, 1, -1]],
+ [-1, 0, -1, 0],
+ [sym.Scalar("i"), sym.Scalar("j"), sym.Scalar("nhalf")],
+ ),
+ ],
+ "nested_loop_partially_used": [
+ Polyhedron(
+ [[-1, 0, 0], [1, 0, -1], [0, -1, 0], [0, 1, -1]],
+ [-1, 0, -1, 0],
+ [sym.Scalar("i"), sym.Scalar("j"), sym.Scalar("nfourth")],
+ ),
+ ],
+ "partially_used_array": [
+ Polyhedron(
+ [[-1, 0], [1, -1]],
+ [-2, 0],
+ [sym.Scalar("i"), sym.Scalar("nhalf")],
+ ),
+ ],
+ },
+ ),
+ ],
+)
+def test_polyhedron_construction_from_nested_loops(
+ here, filename, loop_extractor, polyhedrons_per_subroutine
+):
+ source = Sourcefile.from_file(here / filename)
+
+ for subroutine in source.all_subroutines:
+ expected_polyhedrons = polyhedrons_per_subroutine[subroutine.name]
+
+ list_of_loops = loop_extractor(subroutine.body)
+
+ polyhedrons = [Polyhedron.from_nested_loops(loops) for loops in list_of_loops]
+
+ for polyhedron, expected_polyhedron in zip(polyhedrons, expected_polyhedrons):
+ assert_equal_polyhedron(polyhedron, expected_polyhedron)