Symbolic Polyhedral-Based Energy Analysis for
Nested Loop Programs

TCPAs [6, 4] are parallel processor array architectures designed to execute multidimensional, loop applications represented as polyhedral recurrence equations as piecewise regular algorithms (PRAs). A TCPA works with a host, such as a CPU or FPGA. In this section, we provide an overview of the TCPA architecture. Next, we introduce the PRA representation for expressing loop nests and explain how iteration spaces are partitioned, scheduled, and mapped for spatial execution.

TCPAs are two-dimensional arrays of programmable processing elements (PEs) connected by a configurable circuit-switched neighbor-to-neighbor interconnect, as shown in Figure 1. The PE grid is bordered by four input/output (I/O) buffers, each comprising multiple dedicated address generators (AG). The architecture also includes peripheral control units, a global controller (GC), and a configuration manager (CM), which orchestrates loop-program execution on the array. The TCPA operates as a memory-coupled coprocessor alongside a host system. The loop I/O controller schedules direct memory access (DMA) to transfer data between the host’s external memory and the I/O buffers. After the TCPA is configured by the host, it executes the loop program independently, while the loop I/O controller schedules DMA transfers during execution to support dynamic refilling.

TCPAs avoid costly direct PE-to-DRAM communication by transferring data between the host and I/O buffers rather than having the PEs access them directly. During execution, all active data movement is restricted to the I/O buffers and the PE array [20]. This ensures predictable, high-bandwidth on-chip communication.

Each PE can be configured with multiple parallel functional units (FUs). Each FU has its own instruction memory, branch unit, and program counter. Input and output data streams between the I/O buffers and the PEs, while all intermediate values remain within the PE’s internal register hierarchy. Within each PE, the architecture supports multiple register classes, namely general-purpose, feedback, and input/output registers. The compiler [23] maps data to these register classes based on dependence structures and schedules. During execution, values produced and consumed by operations are stored and forwarded through these registers according to architectural constraints, and the circuit-switched interconnect supports inter-PE communication when dependencies span multiple PEs, while the exact mapping rules are defined in the binding phase and discussed later.

Almost any standard programming language supports loop nests syntactically, e.g., C, C++, Java, and Python. In this paper, we assume loop nests described using a polyhedral notation called Piecewise Linear Algorithms (PLAs) [14, 15]. Normal sequential for-loop specifications can be converted to such a polyhedral description for parallelization. A PLA describes an $n$-dimensional loop nest by an iteration space $\mathcal{I}\subseteq\mathbb{Z}^{n}$, where each element $\mathbf{i}=(i_{0},i_{1},\dots,i_{n-1})^{\mathrm{T}}\in\mathcal{I}$, called iteration vector, corresponds to one loop iteration. In the following, we assume $\mathcal{I}$ is described by a polyhedral description $\mathcal{I}=\{\,\mathbf{i}\mid\mathbf{A}\mathbf{i}\geq\mathbf{b}\,\}$, where $A\in\mathbb{Z}^{m\times n},\;b\in\mathbb{Z}^{n}$. The computations of a loop nest are described by a set of quantified statements $S=\{\ldots,S_{q}\ldots\}$, with $S_{q}$ given by

$\displaystyle S_{q}:$

$\displaystyle x_{q}[\mathbf{P}_{q}\mathbf{i}+\mathbf{f}_{q}]=F_{q}(\ldots,x_{q,r}[\mathbf{Q}_{q,r}\mathbf{i}-\mathbf{d}_{q,r}],\ldots)\;\text{if }\mathbf{i}\in\mathcal{I}_{q}.$

(1)

For all $\mathbf{i}\in\mathcal{I}\cap\mathcal{I}_{q}$ (with $\mathcal{I}_{q}$ called condition space), the variable $x_{q}$ at index (vector) $\mathbf{P}_{q}\mathbf{i}+\mathbf{f}_{q}$ is defined as value of the function $F_{q}$ evaluated on variables $x_{q,r}$ at $\mathbf{Q}_{q,r}\mathbf{i}-\mathbf{d}_{q,r}$. We finally call a PLA Piecewise Regular Algorithm (or PRA) for the special case that $\mathbf{P}_{q}$ and $\mathbf{Q}_{q,r}$ are identity matrices and $\mathbf{f_{q}}$ is zero. As a result, a statement of a PRA has the following form:

$\displaystyle S_{q}:$

$\displaystyle x_{q}[\mathbf{i}]=F_{q}(\ldots,x_{q,r}[\mathbf{i}-\mathbf{d}_{q,r}],\ldots)\;\text{if }\mathbf{i}\in\mathcal{I}_{q}.$

(2)

Note that in a PLA or PRA, there exists neither an explicit order of execution of iterations like in imperative loop programs, nor any order of execution of statements. Instead, each quantified statement only implies implicit schedule restrictions by so-called data dependencies. For example, in a PRA with a statement given in Eq. (2), the left-hand-side variable $x_{q}[\mathbf{i}]$ can only be evaluated after all variables $x_{q,r}[\mathbf{i}-\mathbf{d}_{q,r}]$ (arguments of $F_{q}$) on the right-hand side have been computed. We also say $x_{q}[\mathbf{i}]$ depends on $x_{q,r}[\mathbf{i}-\mathbf{d}_{q,r}]$, and call the constant vector $\mathbf{d}_{q,r}$ the dependence vector. Finally, instances of variables appearing only on the right-hand side of statements are called input variables. Similarly, instances of variables appearing only on the left-hand side of statements are called output variables. All other instances of variables are called internal variables.

Now, in order to map a polyhedral loop nest to a 1- or 2-dimensional array of processing elements, the given iteration space $\mathcal{I}$ is partitioned into as many tiles as available processing elements. For each tile, then a feasible set of iterations within a tile and between tiles needs to be found. This is explained in the following.

According to [23, 14, 13], we partition the $n$-dimensional iteration space $\mathcal{I}$ into congruent rectangular tiles $\mathcal{J}$ by using a partitioning matrix $P={\mathrm{diag}}(p_{0},\ p_{1},\ldots,p_{n-1})$ such that $\mathcal{I}\subseteq\mathcal{J}\boldsymbol{\oplus}P\mathcal{K}$. A tile $\mathcal{J}$ can be described as follows:

$$\mathcal{J}=\bigl\{\,\mathbf{j}=(j_{0},\ldots,j_{\ell},\ldots,j_{n-1})^{\mathrm{T}}\mid 0\leq j_{\ell}<p_{\ell}\,\bigr\}.$$

(3)

Similarly, the set of non-empty tile origins $\mathcal{K}$ can be described as:

$$\mathcal{K}=\bigl\{\,\mathbf{k}=(k_{0},\ldots,k_{\ell},\ldots,k_{n-1})^{\mathrm{T}}\mid 0\leq k_{\ell}<t_{\ell}\,\bigr\}.$$

(4)

The vector $T=(t_{0},\ldots,t_{\ell}\ldots,t_{n-1})$ describes the number of tiles along each dimension $\ell$. After tiling, the original dependencies in statements as shown in Eq. (2) $\mathbf{d}\in\mathbb{Z}^{n}$ are decomposed into an intra-tile dependence vector ($\mathbf{d}_{J}$) and an inter-tile dependence vector ($\mathbf{d}_{K}$) and all instances of variables are embedded into the $2n$-dimensional space with dependence vector ($\mathbf{d}_{J},\mathbf{d}_{K})^{\mathrm{T}}$. Without loss of generality, we can split the statement $S_{q}$ in Eq. (2) into two equations below [14]:

$$S_{q}:\;x_{q}[\mathbf{j},\mathbf{k}]=F_{q}(\ldots,\,x_{q,r}^{\ast}[\mathbf{j},\mathbf{k}],\,\ldots)\;\text{if }\mathbf{j}+P\mathbf{k}\in\mathcal{I}_{q}$$

(5)

	$\displaystyle S_{q}^{\ast}:\;x_{q,r}^{\ast}[\mathbf{j},\,\mathbf{k}]=$	$\displaystyle x_{q,r}[\mathbf{j}-d-P\gamma,\,\mathbf{k}+\gamma]\;\text{if }\mathbf{j}+P\mathbf{k}\in\mathcal{I}_{q}\;$
		$\displaystyle\wedge\mathbf{j}-d-P\gamma\in\mathcal{J}.$		(6)

Whereas the dependencies of the transformed statement $S_{q}$ are all zero, we need to create one statement of the form $S_{q}^{\ast}$ for each solution $\gamma$ that satisfies the set of inequalities

$$\{\gamma\in\mathbb{Z}^{n}:-e<\gamma+P^{-1}d<e\}$$

(7)

and $e=(1\;1\,\cdots 1)^{\mathrm{T}}.$

Now, in order to execute the iterations within each tile and to respect data dependencies both within a tile and across tiles, a feasible schedule must be determined. This is explained in the following.

Whereas tiling defines the mapping of iterations and their corresponding computations to PEs, we still need to find execution schedules that minimize execution time. A schedule assigns a variable $x_{q}$ defined in a statement $S_{q}$ a start time $t_{q}(J,K)$ when the operation $F_{q}$ is computed for each iteration $(J,K)$. According to [12, 23], iterations within a tile must be executed in either a sequential or a pipelined (modulo-scheduled) order with an initiation interval $\pi$ between two iterations. Different tiles are scheduled in parallel across a PE array, such an execution corresponds to a locally sequential, globally parallel (LSGP) modulo schedule. Moreover, the regular computations of loop nests can be efficiently scheduled using linear schedules in which the schedule of iterations within a tile is described by an intra-tile schedule vector $\boldsymbol{\lambda}^{J}$. Similarly, the start time of tiles is also described by a so-called inter-tile schedule vector $\boldsymbol{\lambda}^{K}$. According to [24], schedules $(\boldsymbol{\lambda}^{J},\boldsymbol{\lambda}^{K})$ that minimize the global latency $L$ of a loop nest for a given initiation interval $\pi$ can be determined efficiently as follows:

$$L\;=\;\boldsymbol{\lambda}^{J}\begin{pmatrix}p_{0}-1\\ \vdots\\ p_{n-1}-1\end{pmatrix}\;+\;\boldsymbol{\lambda}^{K}\begin{pmatrix}t_{0}-1\\ \vdots\\ t_{n-1}-1\end{pmatrix}\;+\;L_{c}.$$

(8)

The first term determines the maximal difference of start times between any two tiles. The second term determines the maximum difference of any pair of iterations within a tile. Finally, $L_{c}=\max_{1\leq q\leq|S|}(\tau_{q}+w_{q})$ describes the latency of a single iteration within a tile, where $\tau_{q}$ is the start time and $w_{q}$ the execution latency of operation $F_{q}$ in statement $S_{q}$.

Without loss of generality, we can assume that any PRA loop statement as described in the most general form in Eq. (1) can be split into one statement with only zero dependence vector right-hand side (RHS) arguments $S_{q}$ but containing computations expressed by a function ${\cal F}_{q}$ as described in Eq. (5) and statements that just relate each RHS variable with variables at displacements, i.e., statements of type $S_{q}^{*}$ in Eq. (6). This eases the description of the following energy analysis by splitting the analysis of computational energy and the analysis of memory-related energy. In the following, we denote $C$ the set of computational statements, e.g., and let $T$ denote the set of memory statements.

A processor array, as shown in Figure 2, comprises a memory system consisting of 6 different memory types $\tau\in\mathcal{T}$ as classified in the above table: on-chip I/O buffers ($\mathrm{IO}$) at the periphery, through which data transfers for fetching and storing tensor I/O data from/to a host DRAM ($\mathrm{DR}$) are performed explicitly via DMA operations. Because DRAM accesses are the most expensive, the loop mapping strategy should avoid moving any input or output variable instance of a loop nest across the chip boundary more than once. Moreover, within each PE, we distinguish four classes of registers with the following use: general-purpose registers ($\mathrm{RD}$) for handling intra-iteration dependencies, feedback registers ($\mathrm{FD}$) for temporarily storing variables with inter-iteration dependencies that are processed locally within a PE, and finally input/output registers ($\mathrm{ID}$ and $\mathrm{OD}$) used to receive data from, or transmit data to, neighboring PEs via communication ports. A breakdown of typical energies per access for different types of memory accesses ${\mathrm{T}}=\{\mathrm{RD},\mathrm{FD},\mathrm{ID},\mathrm{OD},\mathrm{IOb},\mathrm{DR}\}$ and for some basic arithmetic operations is given in Table I for a \qty45\nano technology [9].

With the above analysis of energy estimates per statement and iteration, we can finally derive estimates of the total energy $E_{\mathrm{tot}}$ for execution of a parametric loop nest by summing up the energies per statement after tiling multiplied by the volume ${S}_{q}$ (number of integer points of) the polyhedral space where statement $S_{q}$ is defined:

$$E_{\mathrm{tot}}=\sum_{S_{q}\in C}\operatorname{Vol}({S}_{q})\cdot E_{q}^{C}+\sum_{S_{q}\in M}\operatorname{Vol}({S}_{q})\cdot E_{q}^{M}.$$

(11)

In the previous section, we provided a generic energy analysis per statement of a loop nest and a formula for the total energy for executing a given loop nest by multiplying these estimates by the number of integer points in the corresponding parametric polyhedra. For computational statements $S_{q}\in C$ as described in Eq. (5), the number of iterations after tiling is obtained from the tiled iteration space.

$$\operatorname{Vol}({S}_{q})=|\big\{\mathbf{i}\mid\mathbf{i}=\mathbf{j}+P\mathbf{k}\land\mathbf{i}\in\mathcal{I}_{q}\cap\mathcal{I}\big\}|$$

(12)

For memory-related statements of type $S_{q}^{*}\in M$ as described in Eq. (6), the execution count of the associated memory accesses additionally depends on the dependencies $\mathbf{d}_{j}$ and $\mathbf{d}_{k}$. Therefore, for $S_{q}\in T$, the volume is given by

$$\operatorname{Vol}(S_{q}^{\ast})=|\big\{\mathbf{i}\mid\mathbf{i}=\mathbf{j}+P\mathbf{k}\land\mathbf{j}-\mathbf{d}_{j}\in\mathcal{J}\land\mathbf{i}\in\mathcal{I}_{q}\cap\mathcal{I}\big\}|$$

(13)

The computation of such volumes corresponds to counting the number of integer points in parametric polyhedra. This can be performed symbolically using Barvinok’s algorithm [2] as implemented in the Integer Set Library (ISL) [17]. The resulting volumes are returned as piecewise quasi-polynomial functions of the loop bounds $N_{i}$ and tile size parameters $p_{i}$, which can then be inserted into Eq. (11) to determine the total energy consumption $E_{\mathrm{tot}}$ of all statements as executed in a given loop nest.

We validate the accuracy of the symbolic volume computation results by comparing analytically derived integer-point counts for both data transfer and computation against reference counts obtained from a cycle-accurate simulator. The simulator operates using an XML-based architectural description that captures the entire TCPA architecture. Based on this description, the TCPA compiler [23] automatically tiles the input loop programs, maps the tiles onto the processor array, and schedules their execution. During simulation, all data transfers and computations are tracked, yielding exact reference counts. Each selected PolyBench kernel [10] is then simulated for a specified architectural configuration to obtain ground-truth memory access and computation counts.

Note that our analytical method can analyze the energy of a loop kernel symbolically without requiring concrete loop bounds. The resulting volumes are expressed as closed-form quasi-polynomials. In the following, we evaluated our approach across eight different benchmarks from [10] for varying problem sizes and architectural parameters. The analytically derived access counts and obtained total energy values match the simulation results exactly, confirming the accuracy of the symbolic formulation.

Figure 4 compares the analysis time of simulation-based counting and the proposed symbolic method for the GESUMMV benchmark mapped onto an $8\times 8$ PE array. The simulation-based approach exhibits a rapid increase in analysis time with increasing matrix size, as it explicitly executes all loop iterations, instruction activities, and memory accesses within the cycle-accurate model. As the iteration space grows quadratically with the matrix dimension, the simulation cost scales accordingly.

In contrast, the symbolic approach evaluates one-time calculated closed-form expressions for computation and memory access statements, independent of the number of loop iterations. As a result, the analysis time remains almost constant at less than $0.5\ \unit{s}$ across all evaluated problem sizes. This decoupling of analysis cost from the dynamic execution makes the symbolic method scalable and suitable for exploring large loop bounds that are impractical for simulation-based analysis.

Since the proposed analysis derives expressions, it enables a fast evaluation of bounds and exploration of architectural configurations. In the following, we study the total energy $E_{\mathrm{tot}}$ and latency $L$ of a given loop nest with increasing loop bounds. Thereby, we provide a fine-grained breakdown of energy contributions across different access locations $L(x)\in\mathrm{T}$, revealing how memory accesses and computations contribute individually. As will be shown, computation energy remains relatively small compared to data transfer energy across all configurations. The total energy and latency are computed using Eq. (11) and (8), respectively. For this analysis, we consider the GEMM kernel [10], where the iteration space grows cubically with problem size (i.e., $O(N^{3})$). As loop bounds increase, both computation and data movement grow rapidly. Figure 5 illustrates how total energy and latency increase with matrix size for GEMM, analytically evaluated for an $8\times 8$ PE grid. As expected, both metrics grow rapidly with increasing loop bounds due to the cubic growth of the iteration space.

The energy breakdown further reveals how different components scale. For smaller problem sizes, DRAM accesses dominate the total energy consumption. However, as the loop bounds increase, the relative contribution of DRAM energy decreases, while the energy associated with on-chip storage and communication—such as FD and RD registers—as well as computation, increases. This shift is primarily due to larger tile sizes, which increase intra-tile data reuse and consequently amplify activity within local storage locations. This growth is not strictly proportional to the loop bounds, as data accessed from different locations scales differently depending on tile sizes and data dependencies. In many existing energy estimation approaches, activity counts are obtained by simulating a fixed workload instance for smaller loop bounds and extrapolating to larger ones, since simulation is not efficient at handling larger bounds. Our proposed symbolic analysis does not suffer from scalability.

Since not only the energy estimates are parametric but also the schedules, performance metrics such as latency, throughput, and thus also energy efficiency can be computed analytically. This paves the way for a rapid comparison of architectural configurations and supports DSE to identify suitable accelerator architectures for a huge number of loop applications.