Stateful Dataflow Multigraphs: A Data-Centric Model for Performance Portability on Heterogeneous Architectures
Tal Ben-Nun Johannes de Fine Licht Alexandros Nikolaos Ziogas Timo Schneider Torsten Hoefler










 
	language=Python,
	frame=single, 
	framerule=0pt,
	basicstyle=, 
	backgroundcolor=,
	keywordstyle=,
	commentstyle=,
	rulecolor=,
	stringstyle=,
	lineskip=0pt,
	keywords=def, @dace, program, map, lambda, dace, float32, with, as, True, False, return, for, in, class, @staticmethod,
	numbers=left,
	numbersep=5pt,
	showstringspaces=false

















Stateful Dataflow Multigraphs]Stateful Dataflow Multigraphs: A Data-Centric Model for Performance Portability on Heterogeneous Architectures 


Ben-Nun et al.]Tal Ben-Nun, Johannes de Fine Licht, Alexandros N. Ziogas, Timo Schneider, Torsten Hoefler


  Department of Computer Science, ETH Zurich, Switzerland    

talbn,definelicht,alziogas,timos,htor@inf.ethz.ch    



The ubiquity of accelerators in high-performance computing has driven programming complexity beyond the skill-set of the average domain scientist. To maintain performance portability in the future, it is imperative to decouple architecture-specific programming paradigms from the underlying scientific computations.
We present the Stateful DataFlow multiGraph (SDFG), a data-centric intermediate representation that enables separating program definition from its optimization. By combining fine-grained data dependencies with high-level control-flow, SDFGs are both expressive and amenable to program transformations, such as tiling and double-buffering. These transformations are applied to the SDFG in an interactive process, using extensible pattern matching, graph rewriting, and a graphical user interface. We demonstrate SDFGs on CPUs, GPUs, and FPGAs over various motifs — from fundamental computational kernels to graph analytics. We show that SDFGs deliver competitive performance, allowing domain scientists to develop applications naturally and port them to approach peak hardware performance without modifying the original scientific code.






	<ccs2012>
	<concept>
	<concept_id>10011007.10011006.10011008.10011009.10010175</concept_id>
	<concept_desc>Software and its engineering Parallel programming languages</concept_desc>
	<concept_significance>500</concept_significance>
	</concept>
	<concept>
	<concept_id>10011007.10011006.10011008.10011009.10011016</concept_id>
	<concept_desc>Software and its engineering Data flow languages</concept_desc>
	<concept_significance>500</concept_significance>
	</concept>
	<concept>
	<concept_id>10003120.10003121.10003129</concept_id>
	<concept_desc>Human-centered computing Interactive systems and tools</concept_desc>
	<concept_significance>500</concept_significance>
	</concept>	
	<concept>
	<concept_id>10011007.10011006.10011041.10011044</concept_id>
	<concept_desc>Software and its engineering Just-in-time compilers</concept_desc>
	<concept_significance>300</concept_significance>
	</concept>
	</ccs2012>

[500]Software and its engineering Parallel programming languages
[500]Software and its engineering Data flow languages
[500]Human-centered computing Interactive systems and tools
[300]Software and its engineering Just-in-time compilers


 
2019 
2019 
[SC '19]The International Conference for High Performance Computing, Networking, Storage, and AnalysisNovember 17–22, 2019Denver, CO, USA
The International Conference for High Performance Computing, Networking, Storage, and Analysis (SC '19), November 17–22, 2019, Denver, CO, USA
15.00
10.1145/3295500.3356173
978-1-4503-6229-0/19/11


 
[
    [
    December 30, 2023
=====================





§ MOTIVATION




HPC programmers have long sacrificed ease of programming and portability
for achieving better performance. This mindset was established at a time
when computer nodes had a single processor/core and were programmed with
C/Fortran and MPI. 

The last decade, witnessing the end of Dennard scaling and Moore's law,
brought a flurry of new technologies into the compute nodes. Those range
from simple multi-core and manycore CPUs to heterogeneous GPUs and
specialized FPGAs. 

To support those architectures, the complexity of OpenMP's specification grew by more than
an order of magnitude from 63 pages in OpenMP 1.0 to 666 pages in OpenMP
5.0.

This one example illustrates how (performance) programming complexity
shifted from network scalability to node utilization. 

Programmers would now not only worry about communication (fortunately,
the MPI specification grew by less than 4x from MPI-1.0 to 3.1) but also
about the much more complex on-node heterogeneous programming. 

The sheer number of new approaches, such as OpenACC, OpenCL, or CUDA demonstrate
the difficult situation in on-node programming. 

This increasing complexity makes it nearly impossible for domain
scientists to write portable and performant code today.

The growing complexity in performance programming led to a
specialization of roles into domain scientists and performance
engineers.

Performance engineers typically optimize codes by moving functionality
to performance libraries such as BLAS or LAPACK.

If this is insufficient, they translate 
the user-code to optimized versions, often in different languages such as
assembly code, CUDA, or tuned OpenCL.

Both libraries and manual tuning reduce code maintainability, because the optimized
versions are not only hard to understand for the original author (the
domain scientist) but also cannot be changed without major effort.

Code annotations as used by OpenMP or OpenACC do not change the original
code that then remains understandable to the domain programmer.

However, the annotations must re-state (or modify) some of the semantics
of the annotated code (e.g., data placement or reduction operators).

This means that a (domain scientist) programmer who modifies the code,
must modify some annotations or she may introduce hard-to-find
bugs.

With heterogeneous target devices, it now becomes common that the
complexity of annotations is higher than the code they
describe <cit.>.

Thus, scientific programmers can barely manage the complexity of the
code targeted at heterogeneous devices.

The main focus of the community thus moved from scalability to
performance portability as a major research
target <cit.>. 

We call a code-base performance-portable if the domain
scientist's view (“what is computed”) does not change while the code
is optimized to different target architectures, achieving consistently high performance.
The execution should be approximately as performant (e.g., attaining similar ratio of peak performance) as the best-known implementation or theoretical best performance on the target architecture <cit.>.

As discussed before, hardly any existing programming model that supports
portability to different accelerators satisfies this definition.

Our Data-centric Parallel Programming (DAPP) concept addresses
performance portability. It uses a data-centric viewpoint of an
application to separate the roles of domain scientist and
performance programmer, as shown in Fig. <ref>.

DAPP relies on Stateful DataFlow multiGraphs (SDFGs) to represent code
semantics and transformations, and supports modifying
them to tune for particular target architectures.

It bases on the observation that data-movement dominates time and
energy in today's computing systems <cit.> and pioneers the
necessary fundamental change of view in parallel programming.

As such, it builds on ideas of data-centric mappers and schedule annotations 
such as Legion <cit.> and Halide <cit.>
and extends them with a multi-level visualization of data movement, code transformation and compilation for heterogeneous targets,
and strict separation of concerns for programming roles.

The domain programmer thus works in a convenient and well-known language
such as (restricted) Python or MATLAB. 

The compiler transforms the code into an SDFG, on which the performance engineer
solely works on, specifying transformations that match
certain data-flow structures on all levels (from registers to inter-node communication) and modify them.

Our transformation language can implement arbitrary changes to the SDFG
and supports creating libraries of transformations to optimize
workflows.

Thus, SDFGs separate the concerns of the domain scientist and the
performance engineers through a clearly defined interface, enabling
highest productivity of both roles.

We provide a full implementation of this concept in our Data-Centric 
(DaCe) programming environment, which supports (limited) Python, MATLAB, and
TensorFlow as frontends, as well as support for selected DSLs.

DaCe is easily extensible to other frontends through an SDFG builder
interface.

Performance engineers develop potentially domain-specific transformation
libraries (e.g., for stencil-patterns) and can tune them through DaCe's
Interactive Optimization Environment (DIODE).

The current implementation focuses on on-node parallelism as the most
challenging problem in scientific computing today. However, it is
conceivable that the principles can be extended beyond node-boundaries
to support large-scale parallelism using MPI as a backend.



The key contributions of our work are as follows:

  
  * We introduce the principle of Data-centric Parallel Programming,
    in which we use Stateful Dataflow Multigraphs, a data-centric
    Intermediate Representation that enables separating code definition
    from its optimization.
  
  * We provide an open-source implementation[<https://www.github.com/spcl/dace>] of the data-centric environment and its 
    performance-optimization IDE.
  
  * We demonstrate performance portability on fundamental kernels, graph algorithms, and a real-world quantum transport simulator — results are
    competitive with and faster than expert-tuned libraries from Intel and NVIDIA,
    approaching peak hardware performance, and up to five orders of
    magnitude faster than naïve FPGA code written with High-Level
    Synthesis, all from the same program source code.








§ DATA-CENTRIC PROGRAMMING



Current approaches in high-performance computing optimizations <cit.>
revolve around improving data locality. Regardless of the underlying
architecture, the objective is to keep information as close as possible
to the processing elements and promote memory reuse. Even a simple
application, such as matrix multiplication, requires multiple stages of
transformations, including data layout modifications (packing) and
register-aware caching <cit.>. 

Because optimizations do not modify computations and differ for each
architecture, maintaining performance portability of scientific
applications requires separating computational semantics from data
movement.

SDFGs enable separating application development into two stages, as
shown in Fig. <ref>. The problem is formulated as a
high-level program (Fig. <ref>), and is then transformed into a human-readable
SDFG as an Intermediate Representation (IR, Fig. <ref>). The SDFG can then
be modified without changing the original code, and as long as the dataflow aspects do not change, the original code can be updated while keeping SDFG transformations intact. What differentiates the SDFG from other IRs is the ability to hierarchically and parametrically view data movement, where scopes in the graph contain overall data requirements. This enables reusing transformations (e.g., tiling) at different levels of the memory hierarchy, as well as performing cross-level optimizations.

The modifications to the SDFG are not completely
automatic. Rather, they are made by the performance engineer as a
result of informed decisions based on the program structure, hardware
information, and intermediate performance results. 

To support this, a transformation interface and 
common optimization libraries
should be at the performance engineer's disposal,
enabling modification of the IR in a verifiable manner (i.e., without
breaking semantics), either programmatically or interactively.  

The domain scientist, in turn, writes an entire
application once for all architectures, and can freely update the
underlying calculations without undoing optimizations on the SDFG.

Conceptually, we perform the separation of computation from data movement logic by viewing programs as data flowing between operations, much like Dataflow and Flow-Based Programming <cit.>.

One key difference between dataflow and data-centric parallel programming,
however, is that in a pure dataflow model execution is stateless, which
means that constructs such as loops have to be unrolled.
At the other extreme, traditional, control-centric programs revolve
around statements that are executed in order. 

Data-centric parallel programming promotes the use of stateful
	dataflow, in which execution order depends first on data dependencies,
but also on a global execution state. The former fosters the expression
of concurrency, whereas the latter increases expressiveness and
compactness by enabling concepts such as loops and data-dependent
execution.

The resulting concurrency works in several granularities, from
utilizing processing elements on the same chip, to ensuring overlapped
copy and execution of programs on accelerators in clusters.

A data-centric model combines the following concepts:

	
  * Separating Containers from Computation: Data-holding constructs with volatile or non-volatile information are defined as separate entities from computations, which consist of stateless functional units that perform arithmetic or logical operations in any granularity. 
	
  * Dataflow: The concept of information moving from one container or computation to another. This may be translated to copying, communication, or other forms of movement.
	
  * States: Constructs that provide a mechanism to introduce execution order independent of data movement.
	
  * Coarsening: The ability to view parallel patterns in a hierarchical manner, e.g., by grouping repeating computations.


The resulting programming interface should thus enable these concepts without drastically modifying development process, both in terms of languages and integration with existing codebases.



 §.§ Domain Scientist Interface





  
Languages
Scientific applications typically employ different programming models
and Domain-Specific Languages (DSLs) to solve problems. To cater to the
versatile needs of the domain scientists, SDFGs should be easily
generated from various languages. We thus implement SDFG frontends in
high-level languages (Python, MATLAB, TensorFlow), and provide a
low-level (builder) API to easily map other DSLs to SDFGs. In the rest of this section, we focus on the Python <cit.> interface, which is the most extensible.




  
Interface
The Python interface creates SDFGs from restricted Python code, supporting numpy operators and functions, as well as the option to explicitly specify dataflow. In Fig. <ref>, we demonstrate the data-centric interface on a one-dimensional Laplace operator.

DaCe programs exist as decorated, strongly-typed functions in the application ecosystem, so that they can interact with existing codes using array-based interfaces (bottom of figure).

The Python interface contains primitives such as  and  (which translate directly into SDFG components), allows programmers to use multi-dimensional arrays, and implements an extensible subset of operators from numpy <cit.> on such arrays to ease the use of linear algebra operators. For instance, the code  generates the dataflow of a matrix multiplication.



  
Extensibility
For operators and functions that are not implemented, a user can easily provide dataflow implementations using decorated functions () that describe the SDFG. Otherwise, unimplemented functions fall-back into Python, casting the array pointers (which may be defined internally in the DaCe program) into numpy arrays and emitting a “potential slowdown” warning. If the syntax is unsupported (e.g., dynamic dictionaries), an error is raised.



  
Explicit Dataflow
If the programmer does not use predefined operators (e.g., for custom element-wise computation), dataflow “intrinsics” can be explicitly defined separately from code, in constructs which we call Tasklets. Specifically, tasklet functions cannot access data unless it was explicitly moved in or out using pre-declared operators (, ) on arrays, as shown in the code.



Data movement operations (memlets) can be versatile, and the Python syntax of explicit memlets is defined using the syntax shown in Fig. <ref>. First, a local variable (i.e., that can be used in computation) is defined, whether it is an input or an output. After the direction of the movement, the data container is specified, along with an optional range (or index). In some applications (e.g., with indirect or data-dependent access), it is a common occurrence that the subset of the accessed data is known, but not exact indices; specifying memory access constraints both enables this behavior and facilitates access tracking for decomposition, e.g., which data to send to an accelerator. Finally, the two optional values in parentheses govern the nature of the access — the number of data elements moved, used for performance modeling, and a lambda function that is called when write-conflicts may occur. 




Using explicit dataflow is beneficial when defining nontrivial data accesses. Fig. <ref> depicts a full implementation of Sparse Matrix-Vector multiplication (SpMV). In the implementation, the access  is translated into an indirect access subgraph (see Appendix F) that can be identified and used in transformations.






  
External Code
Supporting scientific code, in terms of performance and productivity, requires the ability to call previously-defined functions or invoke custom code (e.g., intrinsics or assembly). In addition to falling back to Python, the frontend enables defining tasklet code in the generated code language directly. In Fig. <ref> we see a DaCe program that calls a BLAS function directly. The semantics of such tasklets require that memlets are defined separately (for correctness); the code can in turn interact with the memory directly (memlets that are larger than one element are pointers). With this feature, users can use existing codes and benefit from concurrent scheduling that the SDFG provides.




  
Parametric Dimensions
To support parametric sizes (e.g., of arrays and maps) in DaCe, we utilize symbolic math evaluation. In particular, we extend the SymPy <cit.> library to support our expressions and strong typing. The code can thus define symbolic sizes and use complex memlet subset expressions, which will be analyzed during SDFG compilation.

The separation of access and computation, flexible interface, and symbolic sizes are the core enablers of data-centric parallel programming, helping domain scientists create programs that are amenable to efficient hardware mapping.









§ STATEFUL DATAFLOW MULTIGRAPHS




We define an SDFG as a directed graph of directed acyclic
multigraphs, whose components are summarized in Table <ref>.
Briefly, the SDFG is composed of acyclic dataflow multigraphs,
in which nodes represent containers or computation, and edges (memlets) represent data movement.
To support cyclic data dependencies and control-flow, these multigraphs reside in State nodes at the top-level graph. Following complete execution of the dataflow in a state, state transition edges on the top-level graph specify conditions and assignments, forming a state machine. For complete operational semantics of SDFGs, we refer to Appendix A.






 §.§ Containers

As a data-centric model, SDFGs offer two forms of data containers: Data and Stream nodes. 

Data nodes represent a location in memory that is mapped to a multi-dimensional array, whereas Stream nodes are defined as multi-dimensional arrays of concurrent queues, which can be accessed using push/pop semantics.

Containers are tied to a specific storage location (as a node property), which may be on a GPU or even a file. In the generated code, memlets between containers either generate appropriate memory copy operations or fail with illegal accesses (for instance, when trying to access paged CPU memory within a GPU kernel).

In FPGAs, Stream nodes instantiate FIFO interfaces that can be used to connect hardware modules.

Another property of containers is whether they are transient, i.e., only allocated for the duration of SDFG execution. This allows transformations and performance engineers to distinguish between buffers that interact with external systems, and ones that can be manipulated (e.g., data layout) or eliminated entirely, even across devices. This feature is advantageous, as standard compilers cannot make this distinction, especially in the presence of accelerators.






 §.§ Computation

Tasklet nodes contain stateless, arbitrary computational functions of any granularity. 

The SDFG is designed, however, for fine-grained tasklets, so as to enable performance engineers to analyze and optimize the most out of the code, leaving computational semantics intact.

Throughout the process of data-centric transformations and compilation, the tasklet code remains immutable.

This code, provided that it cannot access external memory without memlets, can be written in any source language that can compile to the target platform, and is implemented in Python by default.


In order to support Python as a high-level language for tasklets, we implement a Python-to-C++ converter. The converter traverses the Python Abstract Syntax Tree (AST), performs type and shape inference, tracks local variables for definitions, and uses features from C++14 (such as lambda expressions and s) to create the corresponding code.

Features that are not supported include dictionaries, dynamically-sized lists, exceptions, and other Python high-level constructs.

Given that tasklets are meant to be fine-grained, and that our DaCe interface is strongly typed ( <ref>), this feature-set is sufficient for HPC kernels and real-world applications. 






 §.§ Concurrency



Expressing parallelism is inherent in SDFGs by design, supported by the Map and Consume scopes.

Extending the traditional task-based model, SDFGs expose concurrency by grouping parallel subgraphs (computations, local data, movement) into one symbolic instance, enclosed within two “scope” nodes.

Formally, we define an enclosed subgraph as nodes dominated by a scope entry node and post-dominated by an exit node.

The subgraphs are thus connected to external data only through scope nodes, which enables analysis of their overall data requirements (useful, e.g., for automatically transforming a map to GPU code).



Map scopes represent parallel computation on all levels, and can be nested hierarchically. This feature consolidates many parallel programming concepts, including multi-threading, GPU kernels, multi-GPU synchronization, and multiple processing elements on FPGAs.

The semantics of a Map are illustrated in Fig. <ref> — a symbolic integer set attribute of the scope entry/exit nodes called range (Fig. <ref>) defines how the subgraph should be expanded on evaluation (Fig. <ref>).

Like containers, Maps are tied to schedules that determine how they translate to code. When mapped to multi-core CPUs, Map scopes generate OpenMP  loops; for GPUs, device schedules generate CUDA kernels (with the map range as thread-block indices), whereas thread-block schedules determine the dimensions of blocks, emitting synchronization calls () as necessary; for FPGAs, Maps synthesize different hardware modules as processing elements.

Streams can also be used in conjunction with Maps to compactly represent systolic arrays, constructs commonly used in circuit design to represent efficient pipelines, as can be seen in Fig. <ref>. Note that no data is flowing in or out of the Map scope (using empty memlets for the enclosed subgraph) — this would replicate the scope's contents as separate connected components.







Consume scopes enable producer/consumer relationships via dynamic processing of streams. Consume nodes are defined by the number of processing elements, an input stream to consume from, and a quiescence condition that, when evaluated to true, stops processing. An example is shown in Fig. <ref>, which computes the Fibonacci recurrence relation of an input N without memoization. In the SDFG, the value is first pushed into the stream S and asynchronously processed by P workers, with the memlet annotated as  for dynamic number of accesses. The tasklet adds the result to  and pushes two more values to S for processing. The consume scope then operates until the number of elements in the stream is zero, which terminates the program.


Consume scopes are implemented using batch stream dequeue and atomic operations to asynchronously pop and process elements. The potential to encompass complex parallel patterns like work stealing schedulers using high-performance implementations of this node dramatically reduces code complexity.





In order to handle concurrent memory writes from scopes, we define Write-Conflict Resolution memlets. As shown in Fig. <ref>, such memlets are visually highlighted for the performance engineer using dashed lines. Implementation-wise, such memlets can be implemented as atomic operations, critical sections, or accumulator modules, depending on the target architecture and the function.

Reduce nodes complement conflict resolution by implementing target-optimized reduction procedures on data nodes. An example can be seen with a map-reduce implementation of matrix multiplication (Fig. <ref>), where a tensor with multiplied pairs of the input matrices is reduced to the resulting matrix. As we shall show in the next section, this inefficient representation can be easily optimized using data-centric transformations.


Different connected components within an SDFG multigraph also run concurrently (by definition).

Thus, they are mapped to  in OpenMP, different CUDA streams on GPUs, or different command queues on FPGAs. 

These concepts are notoriously cumbersome to program manually for all platforms, where synchronization mistakes, order of library calls, or less-known features (e.g., , non-blocking CUDA streams) may drastically impact performance or produce wrong results. 

Therefore, the SDFG's automatic management of concurrency, and configurable fine-tuning of synchronization aspects by the performance engineer (e.g., maximum number of concurrent streams, nested parallelism) make the IR attractive for HPC programming on all platforms.










 §.§ State


Sequential operation in SDFGs either implicitly occurs following data dependencies, or explicitly specified using multiple states. 

State transition edges define a condition, which can depend on data in containers, and a list of assignments to inter-state symbols (e.g., loop iteration). 

The concept of a state machine enables both complex control flow patterns, such as imperfectly nested loops, and data-dependent execution, as shown in Fig. <ref>.


To enable control-flow within data-flow (e.g., a loop in a map), or a parametric number of state machines, SDFGs can be nested via the Invoke node. 

The semantics of Invoke are equivalent to a tasklet, thereby disallowing access to external memory without memlets.

The Mandelbrot example (Fig. <ref>) demonstrates nested SDFGs. In the program, each pixel requires a different number of iterations to converge. In this case, an invoke node calls another SDFG within the map to manage the convergence loop. 

Recursive calls to the same SDFG are disallowed, as the nested state machine may be inlined or transformed by the performance engineer. 







§ PERFORMANCE ENGINEER WORKFLOW


The Stateful Dataflow Multigraph is designed to expose application data movement motifs, regardless of the underlying computations.
As such, the optimization process of an SDFG consists of finding and leveraging such motifs, in order to mutate program dataflow.
Below, we describe the two principal tools we provide the performance engineer to guide optimization, followed by the process of compiling an SDFG into an executable library.



 §.§ Graph Transformations



Informally, we define a graph transformation on an SDFG as a “find and
replace” operation, either within one state or between several, which
can be performed if all of the conditions match. For general
optimization strategies (e.g., tiling), we provide a standard library of
such transformations, which is meant to be used as a baseline for
performance engineers. Transformations can be designed with symbolic expressions, or specialized for given sizes in order to fine-tune applications. A list of 16 transformations implemented in DaCe can be found in Appendix B. 

Transformations consist of a pattern subgraph and a replacement subgraph. 
A transformation also includes a matching function,
used to identify instances of the pattern in an SDFG, and programmatically verify that requirements are met. To find matching subgraphs in SDFGs, we use the VF2 algorithm <cit.> to find isomorphic subgraphs.

Transformations can be applied interactively, or using a Python API for matching and applying transformations. An example of a full source code of a transformation is found in Appendix D.

Using the transformation infrastructure, we enable seamless knowledge transfer of optimizations across applications.



Two examples of SDFG transformations can be found in Figures <ref> and <ref>. In Fig. <ref>, the transformation detects a pattern L where Reduce is invoked immediately following a Map, reusing the computed values. In this case, the transient array  can be removed (if unused later) and computations can be fused with a conflict resolution, resulting in the replacement R. 
In the second example (Fig. <ref>), a local array, which can be assigned to scratch-pad memory, is added between two map nodes. As a result, the relative indices are changed in all subsequent memlets. 




 §.§ DIODE


SDFGs are intended to be malleable and interactive, which we realize with an Integrated Development Environment (IDE). The Data-centric Integrated Optimization Development Environment, or DIODE (Fig. <ref>), is a specialized IDE for performance engineers that enables editing SDFGs and applying transformations in a guided manner, in addition to the programmatic interface. In particular, performance engineers can:





	
  * interactively modify attributes of SDFG nodes and memlets;
	
  * transform and tune transformation parameters;
	
  * inspect an integrated program view that maps between lines in the domain code, SDFG components, and generated code;
	
  * run and compare historical performance of transformations;
	
  * save transformation chains to files, a form of “optimization version control” that can be used when tuning to different architectures (diverging from a mid-point in the chain); and
	
  * hierarchically interact with large-scale programs, collapsing irrelevant parts and  focusing on bottleneck subgraphs.




We demonstrate the process of interactively optimizing the matrix multiplication SDFG (Fig. <ref>) in a video[<https://www.vimeo.com/301317247>]. Using the IDE and the SDFG representation yields nearly the same performance as Intel's optimized Math Kernel Library (MKL) <cit.> in minutes (<ref>).





 §.§ Compilation Pipeline


Compiling an SDFG consists of three steps:   data dependency inference,   code generation, and   compiler invocation.

In step , data dependencies on all levels are resolved. First, a validation pass is run on the graph to ensure that scopes are correctly structured, memlets are connected properly, and map schedules and data storage locations are feasible (failing when, e.g., FPGA code is specified in a GPU map). 

Then, memlet ranges are propagated from tasklets and containers outwards (through scopes) to obtain the overall data dependencies of each scope, using the image of the scope function (e.g., Map range) on the union of the internal memlet subsets.
This information is later used to generate exact memory copy calls to/from accelerators.

The code generation process of an SDFG (step ) is hierarchical, starting from top-level states and traversing into scopes.
It begins by emitting external interface code and the top-level state machine. Within each state, nodes are traversed in topological order, and a platform-specific dispatcher is assigned to generate the respective code, depending on the node's storage or schedule type. The process continues recursively via map/consume scopes and nested SDFGs, generating heterogeneous codes using several dispatchers. Between states, transitions are generated by emitting for-loops and branches when detected, or using conditional  statements as a fallback. 

In step , we invoke compiler(s) for the generated code according to the used dispatchers. 
The compiled library can then be used directly by way of inclusion in existing applications, or through Python/DaCe.






§ ASSESSING PERFORMANCE WITHOUT TRANSFORMATIONS



To understand how the inherently-concurrent representation of the SDFG creates reasonably performing naïve code, 
we reimplement and run the Polybench <cit.> benchmark over SDFGs, without any optimizing transformations, using the experimental setup of Section <ref>. We show that the representation itself exposes enough parallelism to compete with state-of-the-art polyhedral compilers, outperform them on GPUs, and provide the first complete set of placed-and-routed Polybench kernels on an FPGA.

To demonstrate the wide coverage provided by SDFGs, we
apply the  automatic transformation to offload each Polybench application to the FPGA during runtime, use our simulation flow to verify correctness, and finally
perform the full placement and routing process. The same applies for . We execute all kernels on the accelerators, including potentially unprofitable ones (e.g., including tasklets without parallelism).

The results are shown in Fig. <ref>, comparing SDFGs both with general-purpose compilers (green bars in the figure), and with pattern-matching and polyhedral compilers (blue bars). We use the default Polybench flags, the  dataset size, and select the best performance of competing compilers from the flags specified in the Appendix C. On the CPU, we see that for most kernels, the performance of unoptimized SDFGs is closer to that of the polyhedral compilers than to the general-purpose compilers. The cases where SDFGs are on the order of standard compilers are solvers (e.g., , ) and simple linear algebra (e.g., ). In both cases, data-centric transformations are necessary to optimize the computations, which we exclude from this test in favor of demonstrating SDFG expressibility.

On the GPU, in most cases SDFGs generate code that outperforms PPCG, a tool specifically designed to transform polyhedral applications to GPUs. In particular, the  GPU SDFG is 11.8× faster than PPCG. We attribute these speedups to avoiding unnecessary array copies due to explicit data dependencies, as well as the inherent parallel construction of the data-centric representation.





§ PERFORMANCE EVALUATION




We evaluate the performance of SDFGs on a set of fundamental kernels, followed by three case studies: analysis of matrix multiplication, a graph algorithm, and a real-world physics application.



  
Experimental Setup
We run all of our experiments on a server that contains an Intel 12-core Xeon E5-2650 v4 CPU (clocked at 2.20GHz, HyperThreading disabled, 64 GB DDR4 RAM) and a Tesla P100 GPU (16 GB HBM2 RAM) connected over PCI Express. For FPGA results, we use a Xilinx VCU1525 board, hosting an XCVU9P FPGA and four DDR4 banks at 2400 MT/s. We run the experiments 30 times and report the median running time (including memory copy), where the error-bars indicate 95% confidence interval of all runs, and points are outliers (Tukey fences, 1.5 IQR). For Polybench running times, we use the provided measurement tool, which reports the average time of five runs. All reported results were executed in hardware, including the FPGA.



  
Compilation
Generated code from SDFGs is compiled using  8 for CPU results, CUDA 9.2 for GPU, and Xilinx SDAccel 2018.2 for FPGAs. Compilation flags:  for CPU,  for GPU, and  for FPGA. Fundamental kernels use single-precision floating point types, Polybench uses the default experiment data-types (mostly double-precision), and graphs use integers.




 §.§ Fundamental Computational Kernels



We begin by evaluating 5 core scientific computing kernels, implemented over SDFGs:

	
  * Matrix Multiplication (MM) of two 2,048×2,048 matrices.
	
  * Jacobi Stencil: A 5-point stencil repeatedly computed on a 2,048 square 2D domain for T=1,024 iterations, with constant (zero) boundary conditions.
	
  * Histogram: Outputs the number of times each value (evenly binned) occurs in a 8,192 square 2D image.
	
  * Query: Filters a column of 67,108,864 elements according to a given condition (filters roughly 50%).
	
  * Sparse Matrix-Vector Multiplication (SpMV) of a CSR matrix (8,192 square; 33,554,432 non-zero values).





Our results, shown in Fig. <ref>, are compared with naïve implementations of the code, compiled with GCC, Clang, NVCC, and ICC; Intel MKL <cit.> and HPX <cit.> corresponding library calls for CPU; NVIDIA CUBLAS <cit.>, CUSPARSE <cit.>, and CUB <cit.>, as well as Hybrid Hexagonal Tiling over PPCG <cit.> for GPU; Halide <cit.> (auto- and manually-tuned) for CPU and GPU; and Xilinx Vivado HLS/SDAccel <cit.> and Spatial <cit.> for FPGA.

On all applications, our SDFG results only employ data-centric transformations, keeping tasklets intact (<ref>). We highlight key results for all platforms below.





In MM, a highly tuned kernel, SDFGs achieve ∼98.6% of the performance of
MKL, ∼70% of CUBLAS, and ∼90% of CUTLASS, which is the upper bound
of CUDA-based implementations of MM. On FPGA, SDFGs yield a result
4,992× faster than naïve HLS over SDAccel. We also run the
FPGA kernel for 1024×1024 matrices and compare to the runtime of
878 ms reported for Spatial <cit.> on the same chip. We
measure 45 ms, yielding a speedup of 19.5×.




We observe similar results in SpMV, which is more complicated to optimize due to its irregular memory access characteristics. SDFGs are on par with MKL (99.9% performance) on CPU, and are successfully vectorized on GPUs.





For Histogram, SDFGs enable vectorizing the program, achieving 8× the performance of gcc, where other compilers fail due to the kernel's data-dependent accesses. 
We implement a two-stage kernel for the FPGA, where the first stage reads
16 element vectors and scatters them to 16 parallel processing elements
generated with map unrolling (similar to Fig. <ref>), each
computing a separate histogram. In the second stage, the histograms are merged
on the FPGA before copying back the final result. This yields a 10×
speedup in hardware. 




In Query, SDFGs are able to use streaming data access to parallelize the
operation automatically, achieving significantly better results than HPX
and STL. On FPGA we read wide vectors, then use a deep pipeline
to pack the sparse filtered vectors into dense vectors. This scheme yields a 10× speedup, similar to Histogram. 




For Jacobi on CPU, we use a domain-specific transformation (). We see that it outperforms standard implementations
by up to two orders of magnitude, performing 90× faster than Polly and marginally outperforming
Pluto, which uses a similar tiling technique. 
In Halide, when all timesteps are hard-coded in the pipeline (which takes ∼68 minutes to compile due to the stateless dataflow representation), its auto-tuner yields the best result, which is 20% faster than SDFGs. For FPGAs, Jacobi is
mapped to a systolic array of processing elements, allowing it to scale up with
FPGA logic to 139 GOp/s.
Overall, the results indicate that data-centric transformations can yield competitive performance across both architectures and memory access patterns.




 §.§ Case Study I: Matrix Multiplication


The transformation chain and performance results from Fig. <ref> to the MM CPU SDFG in the previous section are shown in Fig. <ref>.

After fusing the map and reduction into a write-conflict resolution memlet ( transformation), we largely follow the approach of Goto and van de Geijn <cit.>. Manually reordering the map dimensions in the SDFG reduces data movement and yields a marginal improvement. We then tile for the L3 cache () and tile again for mapping to vector registers. To use caches efficiently, we pack the tiles of B and store tiles of C using the  transformation (twice), followed by  to generate code that uses vector extensions. Lastly, we apply a transformation to convert the internal write-conflict resolution memlet into a state machine, peeling the loop ().

The figure shows that not all transformations yield immediate speedups, yet they are necessary to expose the next steps. Moreover, after only 7 steps the performance increases by ∼536x (75% of MKL), and further increases to 98.6% of MKL after tuning transformation parameters for a specific size (Fig. <ref>).








 §.§ Case Study II: Graph Algorithms



We implement an irregular computation problem on multi-core CPUs: Breadth-First Search (BFS). 
We use the data-driven push-based algorithm <cit.>, and test five graphs with different characteristics as shown in Appendix E.

Due to the irregular nature of the algorithm, BFS is not a trivial problem to optimize. However, SDFGs inherently support constructing the algorithm using streams and data-dependent map ranges. The primary state of the optimized SDFG is shown in Fig. <ref>, which contains only 14 nodes (excluding input/output data) and is the result of three transformations from the base Python program. In particular, the initial generated OpenMP code does not map well to cores due to the dynamic scheduling imposed by the frontier. We mitigate this effect by applying   with T tiles. Since the accesses to the next frontier are performed concurrently through , major overhead was incurred due to atomic operations. Using  , the results are accumulated to  and bulk updates are performed to the global frontier. Similarly, we use   for frontier size accumulation.






We compare our results with two state-of-the-art CPU graph processing frameworks: Galois <cit.> and Gluon <cit.>. We use the default settings ( for Gluon,  for Galois) and use 12 threads (1 thread per core).

In Fig. <ref>, we see that SDFGs perform on-par with the frameworks on all graphs, where Galois is marginally faster on social networks  (∼1.53× on ) and the Kronecker graph. However, in road maps (, ) SDFGs are up to 2× faster than Galois. This result could stem from our fine-grained scheduling imposed by the map scopes.










 §.§ Case Study III: Quantum Transport


Quantum Transport (QT) Simulation is used for estimating heat dissipation in nanoscale transistors. OMEN <cit.> is an efficient (two-time Gordon Bell finalist) QT simulator based on a nonlinear solver, written in C++ and CUDA using MKL, CUBLAS, and CUSPARSE.

We use SDFGs and transformations to optimize the computation of OMEN, the full details of which are described by Ziogas et al. <cit.>. A significant portion of the OMEN runtime is spent in computing Scattering Self-Energies (SSE), which is given by the formula in Fig. <ref> (top-left).

Here we focus on the computation of Σ^≷. Upon generating the SDFG from the formula, we see that there are multiple parallel dimensions that compute small matrix multiplications (denoted as ) and Hadamard products (), reducing the result with a summation. 

In step , we split the parallel computation into several maps, creating transient 5D and 6D arrays. Steps  and  transform the data layout by reorganizing the map dimensions and transient arrays such that the multiplications can be implemented with one “batched-strided” optimized CUBLAS GEMM operation. Lastly, step  replaces the CUBLAS call with a specialized nested SDFG that performs small-scale batched-strided matrix multiplication (SBSMM), transformed w.r.t. matrix dimensions and cache sizes in hardware.





Table <ref> lists the overall SSE runtime simulating a 4,864 atom nanostructure over OMEN, numpy (using MKL for dense/sparse LA), and DaCe. Using data-centric transformations, the under-utilization resulting from the multitude of small matrix multiplications is mitigated, translating into a 32.26× speedup for SDFGs over manually-tuned implementations, and 1,021× over Python. 

Breaking down the speedup, the specialized SDFG implementation of SBSMM (Table <ref>) outperforms CUBLAS by up to 4.76×. 
This demonstrates that performance engineers can use the data-centric view to easily tune dataflow for specific problems that are not considered by existing library function implementations.




§ RELATED WORK


Multiple previous works locally address issues posed in this paper. We discuss those papers below, and summarize them in Fig. <ref>.




  
Performance Portability
Performance-portable programming models consist of high-performance libraries and programming languages. Kokkos <cit.> and RAJA <cit.> are task-based HPC libraries that provide parallel primitives (e.g., forall, reduce) and offer execution policies/spaces for CPUs and GPUs. The Legion <cit.> programming model defines hierarchical parallelism and asynchronous tasking by inferring data dependencies from access sets, called logical regions, which are composed of index spaces and accessed fields. Logical regions are similar to subsets in memlets, but differ by implicitly representing data movement. Legate <cit.> is an implementation of the numpy interface that uses Legion to target distributed GPU systems from Python. Language and directive-based standards such as OpenCL <cit.>, OpenACC, and OpenMP <cit.> also provide portability, where some introduce FPGA support through extensions <cit.>. SYCL <cit.> is an embedded DSL standard extending OpenCL to enable single-source (C++) heterogeneous computing, basing task scheduling on data dependencies. However, true performance portability cannot be achieved with these standards, as optimized code/directives vastly differ on each platform (especially in the case of FPGAs). Other frameworks mentioned below <cit.> also support imperative and massively parallel architectures (CPUs, GPUs), where Halide and Tiramisu have been extended <cit.> to target FPGA kernels. As opposed to SDFGs, none of the above models were designed to natively support both load/store architectures and reconfigurable hardware.



  
Separation of Concerns
Multiple frameworks explicitly separate the computational algorithm from subsequent optimization schemes. In CHiLL <cit.>, a user may write high-level transformation scripts for existing code, describing sequences of loop transformations. These scripts are then combined with C/C++, FORTRAN or CUDA <cit.> programs to produce 
optimized code using the polyhedral model. Image processing pipelines written in the Halide <cit.> embedded DSL are defined as operations, whose schedule is separately generated in the code by invoking commands such as , , and . Tiramisu <cit.> operates in a similar manner, enabling loop and data manipulation. In SPIRAL <cit.>, high-level specifications of computational kernels are written in a DSL, followed by using breakdown and rewrite rules to lower them to optimized algorithmic implementations.
SDFGs, along with DIODE and the data-centric transformation workflow, streamlines such approaches and promotes a solution that enables knowledge transfer of optimization techniques across applications.





  
Dataflow Representations
Several IRs combine dataflow with control flow in graph form. The LLVM IR <cit.> is a control-flow graph composed of basic blocks of statements. Each block is given in Single Static Assignment (SSA) form and can be transformed into a dataflow DAG. Program Dependence Graphs (PDGs) <cit.> represent statements and predicate expressions with nodes, and the edges represent both data dependencies and control flow conditions for execution.

The PDG is a control-centric representation in which statements are assumed to execute dynamically, a model that fits instruction-fetching architectures well. SDFGs, on the other hand, define explicit state machines of dataflow execution, which translate natively to reconfigurable hardware.
Additionally, in PDGs data-dependent access (e.g., resulting from an indirection) creates the same edge as an access of the full array. This inhibits certain transformations that rely on such accesses, and does not enable inferring the total data movement volume, as opposed to the memlet definition.

SDFGs can be trivially converted to the SSA and PDG representations. In the other direction, however, the parametric concurrency that the SDFG scopes offer would be lost (barring specific application classes, e.g., polyhedral).
As the representations operate on a low level, with statements as the unit element, they do not encapsulate the multi-level view as in SDFGs.
HPVM <cit.> extends the LLVM IR by introducing hierarchical dataflow graphs for mapping to accelerators, yet still lacks a high-level view and explicit state machines that SDFGs offer.

Other representations include Bamboo <cit.>, an object-oriented dataflow model that tracks state locally through data structure mutation over the course of the program; Dryad <cit.> and Naiad <cit.>, parametric graphs intended for coarse-grained distributed data-parallel applications, where Naiad extends Dryad with definition of loops in a nested context; simplified data dependency graphs for optimization of GPU applications <cit.>; deterministic producer/consumer graphs <cit.>; and other
combinations of task DAGs with data movement <cit.>. As the SDFG provides general-purpose state machines with dataflow, all the above models can be fully represented within it, where SDFGs have the added benefit of encapsulating fine-grained data dependencies.



  
Data-Centric Transformations
Several representations <cit.> provide a fixed set of high-level program transformations, similar to those presented on SDFGs. In particular, Halide's schedules are by definition data-centric, and the same applies to polyhedral loop transformations in CHiLL. HPVM also applies a number of optimization passes on a higher level than LLVM, such as tiling, node fusion, and mapping of data to GPU constant memory. Lift <cit.> programs are written in a high-level functional language with predefined primitives (e.g., map, reduce, split), while a set of rewrite rules is used to optimize the expressions and map them to OpenCL constructs.
Loop transformations and the other aforementioned coarse-grained optimizations are all contained within our class of data-centric graph-based transformations, which can express arbitrary data movement patterns.



  
Decoupling Data Access and Computation
The Chapel <cit.> language supports controlling and reasoning about locality by defining object locations and custom iteration spaces. Charm++ <cit.> is a parallel programming framework, organized around message passing between collections of objects that perform tasks in shared- or distributed-memory environments.
In the Sequoia <cit.> programming model, all communication is made explicit by having tasks exchange data through argument passing and calling subtasks.  MAPS <cit.> separates data accesses from computation by coupling data with their memory access patterns. This category also includes all frameworks that implement the polyhedral model, including CHiLL, PENCIL <cit.>, Pluto <cit.>, Polly <cit.> and Tiramisu. Furthermore, the concept of decoupling data access can be found in models for graph analytics <cit.>, stream processing <cit.>, machine learning <cit.>, as well as other high-throughput <cit.> and high-performance <cit.> libraries. Such models enable automatic optimization by reasoning about accessed regions. However, it is assumed that the middleware will carry most of the burden of optimization, and thus frameworks are tuned for existing memory hierarchies and architectures. This does not suffice for fine-tuning kernels to approach peak performance, nor is it portable to new architectures. In these cases, a performance engineer typically resorts to a full re-implementation of the algorithm, as opposed to the workflow proposed here, where SDFG transformations can be customized or extended.




  
Summary In essence, SDFGs provide the expressiveness of a general-purpose programming language, while enabling performance portability without interfering with the original code.
Differing from previous models, the SDFG is not limited to specific application classes or hardware, and the extensible data-centric transformations generalize existing code optimization approaches.





§ CONCLUSION

In this paper, we present a novel data-centric model for producing high-performance computing applications from scientific code. Leveraging dataflow tracking and graph rewriting, we enable the role of a performance engineer — a developer who is well-versed in program optimization, but does not need to comprehend the underlying domain-specific mathematics. We show that by performing transformations on an SDFG alone, i.e., without modifying the input code, it is possible to achieve performance comparable to the state-of-the-art on three fundamentally different platforms.

The IR proposed in this paper can be extended in several directions. Given a collection of transformations, research may be conducted into their systematic application, enabling automatic optimization with reduced human intervention. More information about tasklets, such as arithmetic intensity, can be recovered and added to such automated cost models to augment dataflow captured by memlets. Another direction is the application of SDFGs to distributed systems, where data movement minimization is akin to communication-avoiding formulation.


We thank Hussein Harake, Colin McMurtrie, and the whole CSCS team granting
access to the Greina and Daint machines, and for their excellent technical
support.

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 programme (grant agreement DAPP, No. 678880).

T.B.N. is supported by the ETH Zurich Postdoctoral Fellowship and Marie Curie Actions for People COFUND program. 


	
















§ SDFG OPERATIONAL SEMANTICS



We begin by defining the different elements of our IR, which is a graph
of graphs.  We follow by defining how a function expressed in our IR must be
called and then give semantic rules of how an SDFG is
evaluated. We precede each operational semantic rule by a less formal
description of the intent of the rule.



 §.§ Elements of an SDFG


An SDFG is a directed multigraph defined by the tuple (S, T, s_0), whose vertex set S
represent states, its edges T represent transitions between
them, and has one start-state s_0. It is a multigraph since there can be multiple transitions between a pair of states. 

An SDFG state s∈ S is a named, directed acyclic multigraph defined by the tuple (V, E, name).
Each node v_i ∈ V is of one of the
following types, as shown in Table 1 of the paper: data, tasklet, mapentry, mapexit, stream,
reduce,  consume-entry,
consume-exit, invoke. Each node type defines connectors — attachment points for edges defining the nature of the connection, and the edges E indicate dependencies which constrain the execution order.
Each edge carries a memlet, an annotation that specifies dataflow, as well as the connectors on its source and destination nodes. Below, we describe the node types in SDFG states and the exact definition of memlets.

As opposed to inherently sequential representations (cf. <cit.>), in SDFGs the execution order
is mainly constrained by explicit dependencies. 

To parameterize some aspects of the graph, there exists a global namespace for
symbols. 
Symbols can be used in all symbolic expressions mentioned below, and at
runtime they evaluate to scalar values. 
( <ref>).

A data node represents a typed location in memory. The memory itself
can be allocated outside of the SDFG and passed as a pointer upon invocation,
or allocated at runtime if transient. A data node is a tuple (name, basetype,
dimensions, transient).  The name is an identifier, the basetype
a basic type (e.g., , , ), and dimensions a list of symbolic
integer expressions.  Memory pointed to by differently named data nodes must
not alias. A data node has an implicit connector for each adjacent edge.

A tasklet node represents computation. It is defined by the tuple
(inputs, outputs, code), where inputs and outputs are sets of elements of the form
(name, basetype, dimensions) that define connectors with the same name, representing the only external memory the code can read or write. Code is a string, which can be written in Python or languages supported by the target architecture, such as C++. 

A stream node represents an array of typed locations in memory, accessed
using queue semantics, i.e., using  and  operations. A stream has the input connectors data and push, as well as the output
connectors data and pop. The data connectors allow
initialization/copying of the data stored in the queues from/to an array,
push and pop enqueue and dequeue elements, respectively.

The memlet annotation represents dataflow,
defined by the tuple (src_conn, dst_conn, subset, reindex, accesses, wcr). 
The subset function selects which
subset of elements visible at the source connector will flow to the destination
connector. The reindex function specifies at which indices the data will
be visible at the destination node. We express subset and reindex as functions on integer sets, and implement them as lists of exclusive ranges, where each range refers to one data dimension and defined by start:end:stride:tilesize. 
tilesize is used to propagate multiple elements at a time, e.g., vector types. Fig. <ref> extends Fig. 5b from the paper by showing the same memlets in subset/reindex notation rather than array index labels.




The wcr optional attribute (write-conflict resolution) is a function S
× S → S for a data type S.
The function receives the old value currently present at the destination and the new
value, and produces the value that will be written. 
This attribute is normally defined when data flows concurrently from multiple sources, e.g., in a parallel dot product.

Map entry/exit nodes are defined by a tuple (range, inputs, outputs), creating a scope.
The subgraph within the scope of a map is defined by the nodes that are dominated by the map-entry and post-dominated by the map-exit. This subgraph is expanded at runtime and executed concurrently according to the symbolic range of the mapentry, which takes the form identifier=begin:end:step. Input and output connectors of map nodes are either defined as / ( is any string) to map memlets outside the scope to memlets within it; or other identifiers, which can be used to create data-dependent map ranges.
The range identifier becomes a scope-wide symbol, and can used to define memlet attributes (subset, reindex).  


The invoke node allows to call an SDFG within an SDFG state. Its
semantics are equivalent to a tasklet that has input connectors for each
data-node and undefined symbol used in the invoked SDFG, and an output node for
each data node. 

The reduce alias node is used to implement target-optimized reduction procedures. It is defined by axes and a wcr function, consisting of an input and an output connector (of the same type).
It is semantically equivalent to a map whose range spans the incoming and outgoing memlet data, containing an identity tasklet () and an output memlet annotated by the given wcr for the reduced axes. 


Consume entry/exit nodes complement maps, enabling producer/consumer relationships via dynamic parallel processing of streams. The two nodes form a scope, similarly to a map, but also accept a stream_in connector (connected to a stream), and a quiescence condition after which processing stops. Thus, a consume entry node is a tuple (range, cond, inputs, outputs). Semantically, a consume scope is equivalent to a map using the same range, with an internal invoked SDFG containing the scope. This state connects back to itself without any assignments and the condition cond.

A state transition is defined as a tuple (source, destination,
	condition, assignments). Condition can depend on 
symbols and data from data nodes in the source state, whereas assignments take the form symbol = expression. Once a state finishes execution, all outgoing state transitions of
that state are evaluated in an arbitrary order, and the destination of the first
transition whose condition is true is the next state which will be executed. 
If no transition evaluates to true, the program terminates.
Before starting the execution of the next state, all assignments are performed,
and the left-hand side of assignments become symbols.




 §.§ Operational Semantics





  §.§.§ Initialization




  
Notation
We denote collections (sets/lists) as capital
letters and their members with the corresponding lowercase letter and a
subscript, i.e., in an SDFG G=(S, T, s_0) the set of states S as s_i, with 0 ≤ i < |S|.
Without loss of generality we assume s_0 to be the start state.  We denote
the value stored at memory location a as M[a], and assume all basic types
are size-one elements to simplify address calculations.

The state of execution is denoted by ρ. Within the state we carry several
sets: loc, which maps names of data nodes and transients to memory addresses;
sym, which maps symbol names (identifiers) to their current value;
and vis, which maps connectors to the data visible at that connector in the current state of execution.

We define a helper function size(), which returns the product of all
dimensions of the data node or element given as argument (using ρ to resolve symbolic values).
Furthermore, id() returns the name property of a data or transient node, and offs() the offset of a data element relative to the start of the memory region it is stored in.
The function copy() creates a copy of the object given as argument, i.e., when we modify the copy, the original object remains the same.



  
Invocation
When an SDFG G is called with the data arguments A ≡ [a_i = p_i] (a_i is an identifier, p_i is
an address/pointer) and symbol arguments Z ≡ [z_i = v_i] (z_i is an
identifier, v_i an integer) we initialize the configuration ρ:

	
  * For all symbols z_i in Z: sym[z_i] ← v_i.
	
  * For all data and stream nodes d_i ∈ G without incoming edges s.t. id(d_i) = a_i: loc(d_i) ← p_i, vis(d_i.data) ← M[p_i, ..., p_i + size(d_i)].
	
  * Set current to a copy of the start state  of G, s_0.
	
  * Set state to id(s_0).
	
  * Set qsize[f_i] to zero for all stream nodes f_i ∈ G.
	

This can be expressed as the following rule:




    [                G=(S,T),;   𝑠𝑡𝑎𝑟𝑡_𝑠𝑡𝑎𝑡𝑒(G) → s_0,;  D: data nodes in  s_0,; F: stream nodes in  s_0;                         ]/[                                (𝑐𝑎𝑙𝑙(G, A, Z), ρ) →ρ[                                        𝑠𝑡𝑎𝑡𝑒↦𝑖𝑑(s_0),;                                                                                          𝑐𝑢𝑟𝑟𝑒𝑛𝑡↦𝑐𝑜𝑝𝑦(s_0),;                                                                              ∀ a_i=p_i ∈ A: 𝑙𝑜𝑐[a_i ↦ p_i],;                                                                              ∀ z_i=v_i ∈ Z: 𝑠𝑦𝑚[z_i] ↦ v_i,;                                                       ∀ d_i ∈ D: 𝑣𝑖𝑠[d_i.𝑑𝑎𝑡𝑎] ↦ M[p_i, .., p_i+size(d_i)],;                                                                            ∀ f_i ∈ G.S: 𝑞𝑠𝑖𝑧𝑒[𝑖𝑑(f_i)] ↦ 0] ]





  §.§.§ Propagating Data in a State



Execution of a state entails propagating data along edges, governed by the rules defined below.

Data is propagated within a state in an arbitrary order, constrained only by
dependencies expressed by edges. We begin by considering all nodes in the current state, which we add to current. We gradually remove nodes as they are processed, until current is empty and we can proceed to the next state or termination.

To keep the description concise, we omit details of address translation and indexing. We also use the subset and reindex properties as functions that resolve symbols and return data elements directly.




  
Element Processing
In each step, we take one element q (either a memlet or a node) of 𝑐𝑢𝑟𝑟𝑒𝑛𝑡,
for which all input connectors have visible data, then:

If q is a memlet (src, dst, subset, reindex, wcr), update 𝑣𝑖𝑠[𝑑𝑠𝑡] to 𝑤𝑐𝑟(𝑟𝑒𝑖𝑛𝑑𝑒𝑥(𝑠𝑢𝑏𝑠𝑒𝑡(𝑣𝑖𝑠[𝑠𝑟𝑐]))):


    [          q = 𝑚𝑒𝑚𝑙𝑒𝑡(𝑠𝑟𝑐, 𝑑𝑠𝑡, 𝑠𝑢𝑏𝑠𝑒𝑡, 𝑟𝑒𝑖𝑛𝑑𝑒𝑥, 𝑤𝑐𝑟),;                                    (𝑣𝑖𝑠[𝑠𝑟𝑐], ρ) ≠∅,; (𝑤𝑐𝑟(𝑟𝑒𝑖𝑛𝑑𝑒𝑥(𝑠𝑢𝑏𝑠𝑒𝑡(𝑣𝑖𝑠[𝑠𝑟𝑐]))), ρ) → [d_0, .., d_n] ]/
    	(q, ρ) →ρ[ 𝑣𝑖𝑠[𝑑𝑠𝑡] ↦ [d_0, .., d_n] ]



If q is a data node, update its referenced memory for an input connector c_i,
 M[𝑙𝑜𝑐(𝑖𝑑(q)), .., 𝑙𝑜𝑐(𝑖𝑑(q))+𝑠𝑖𝑧𝑒(𝑣𝑖𝑠[𝑞.𝑑𝑎𝑡𝑎])] = 𝑣𝑖𝑠[𝑞.𝑑𝑎𝑡𝑎]:


    [                                 q = 𝑑𝑎𝑡𝑎(𝑖𝑑, 𝑑𝑖𝑚𝑠, 𝑏𝑡, 𝑡𝑟𝑎𝑛𝑠𝑖𝑒𝑛𝑡),;                            (∀ x,q.c_i ∈𝑐𝑢𝑟𝑟𝑒𝑛𝑡: 𝑣𝑖𝑠[q.c_i], ρ) ≠∅,; (∀ x,q.c_i ∈𝑐𝑢𝑟𝑟𝑒𝑛𝑡: 𝑣𝑖𝑠[q.c_i], ρ) → [d^i_0, d^i_1, .., d^i_k_i],;                         ∀ d^i_j : loc(id(q)) + 𝑜𝑓𝑓𝑠(d^i_j) = l^i_j ]/
    	(q, ρ) →ρ [∀ i,k_i: M[l^i_j, .., l^i_j+𝑠𝑖𝑧𝑒(d^i_j)] = d^i_j ]



If q is a tasklet, generate a prologue that allocates local variables for
all input connectors c_i of q, initialized to vis[c_i] (P_1), as well
as output connectors (P_2).  Generate an epilogue Ep which updates
ρ[vis[c_i] ↦ v_i ] for each output connector c_i of q with the
contents of the appropriate variable (declared in P_2). Execute the
concatenation of (P_1; P_2; code; Ep;).  


    [                              q = 𝑡𝑎𝑠𝑘𝑙𝑒𝑡(𝐶𝑖𝑛, 𝐶𝑜𝑢𝑡, 𝑐𝑜𝑑𝑒),;                    (∀ x,q.c_i ∈𝑐𝑢𝑟𝑟𝑒𝑛𝑡: vis[q.c_i], ρ) ≠∅,; (
    	P_1 = [∀ c_i ∈ Cin: type(c_i) id(c_i) = vis[c_i]];, ρ),;                  P_2 = [∀ c_i ∈ Cout: type(c_i) id(c_i);],;            (Ep = [∀ c_i ∈ Cout: &vis(c_i) = id(c_i);], ρ),;                                                            ]/[ (q, ρ) →ρ [ exec(P_1; P_2; code; Ep;) ] ]



If q is a mapentry node with range y = R (y is the identifier) and scope o⊂ V: Remove o
from current. Remove q and the corresponding map exit node from o. For
each element in r_i ∈ R, replicate o, resolve any occurrence of y to
r_i, connect incoming connectors of q and p in state.


    [                 q = 𝑚𝑎𝑝𝑒𝑛𝑡𝑟𝑦(𝐶𝑖𝑛, 𝐶𝑜𝑢𝑡, 𝑅),;                    ∀ cin_i : vis[cin_i] ≠∅,; o = 𝑠𝑐𝑜𝑝𝑒(q), o' = 𝑠𝑐𝑜𝑝𝑒(q) ∖{q, 𝑚𝑒𝑥𝑖𝑡(q)},;            𝑁𝑒𝑤𝑆𝑦𝑚𝑠 = [cin_i : ∄ (cin_i, y)] ]/[                                     (q, ρ) →ρ[ 𝑐𝑢𝑟𝑟𝑒𝑛𝑡↦ o' ∪ [∀ r_i : 𝑟𝑒𝑠𝑠𝑦𝑚(𝑐𝑜𝑝𝑦(o'), r_i)],;                                                       ∀ n_i ∈𝑁𝑒𝑤𝑆𝑦𝑚: sym[n_i] ↦ vis[cin_i]  ] ]



If q is a consume-entry node, defined by (range, cond, cin, cout), replace q with a mapentry 
and do the same
for the corresponding consume exit node. Then we create a new SDFG new, which contains the contents of the original consume scope scope(q). new consists of one state s_0, and a single state transition to the same state with the condition cond, defined by (s_0, s_0, cond, []). Finally, we replace scope(q) in current with an invoke node for
new and reconnect the appropriate edges between the entry and exit nodes. 


    [                    q = 𝑐𝑜𝑛𝑠𝑢𝑚𝑒-𝑒𝑛𝑡𝑟𝑦(𝑟𝑎𝑛𝑔𝑒, 𝑐𝑜𝑛𝑑, 𝑐𝑖𝑛, 𝑐𝑜𝑢𝑡); 𝑛𝑒𝑤𝑠𝑑𝑓𝑔 = 𝑆𝐷𝐹𝐺(𝑠𝑐𝑜𝑝𝑒(𝑞) ∖{𝑞, 𝑐𝑒𝑥𝑖𝑡(𝑞)}, (𝑠_0,𝑠_0, 𝑐𝑜𝑛𝑑, []));                                         iv = 𝑖𝑛𝑣𝑜𝑘𝑒(𝑛𝑒𝑤𝑠𝑑𝑓𝑔);                             men = 𝑚𝑎𝑝𝑒𝑛𝑡𝑟𝑦(𝑟𝑎𝑛𝑔𝑒, 𝑐𝑖𝑛, 𝑐𝑜𝑢𝑡);                   mex = 𝑚𝑎𝑝𝑒𝑥𝑖𝑡(𝑐𝑒𝑥𝑖𝑡(𝑞).𝑐𝑖𝑛, 𝑐𝑒𝑥𝑖𝑡(𝑞).𝑐𝑜𝑢𝑡) ]/
    	(q, ρ) →ρ[ 𝑐𝑢𝑟𝑟𝑒𝑛𝑡↦ (𝑐𝑢𝑟𝑟𝑒𝑛𝑡∖q, cexit(q)) ∪{ iv, 𝑚𝑒𝑛, 𝑚𝑒𝑥} ]



If q is a reduce node defined by the tuple (cin, cout, range), we create a mapentry
node men with the same range, a mapexit node mex, and a tasklet . We add these nodes to the node set of current, nd(current). We
connect them by adding edges to the edge set of current.


    [      q = 𝑟𝑒𝑑𝑢𝑐𝑒(𝑐𝑖𝑛, 𝑐𝑜𝑢𝑡, 𝑟𝑎𝑛𝑔𝑒, 𝑤𝑐𝑟),;                            vis[cin] ≠∅,; men = mapentry({IN_1}, {OUT_1}, range),;         mex = mapexit({IN_1}, {OUT_1}),;            t = tasklet({i},{o}, “o=i;") ]/=1.4pt
    	[                                                      (q, ρ) →ρ [                         nd(current) ↦ nd(current) ∪{men, mex, t};                                                                  ed(current) ↦ ed(current) ∪{ (men.OUT_1, t.i), (t.o, mex.IN_1) };                                                                                                       vis[men.IN_1] ↦ vis[q.cin],;                                                                                                                    vis[q.IN] ↦∅ ] ]



If q is a stream, we add the data visible at the push connector to the
appropriate memory location (indicated by qsize[id(q)]), and increment it.
A stream node with an edge connected to its pop connector can only be evaluated if qsize[id(q)] ≥ 1. When propagating through streams, there are three cases for one step:

  Data is both pushed into and popped from the stream:


    [ q = 𝑠𝑡𝑟𝑒𝑎𝑚(𝑝𝑢𝑠ℎ, 𝑝𝑜𝑝, 𝑠𝑖𝑧𝑒),;        vis[push] ≠ emptyset,;    ∃ (q.pop, dst) ∈ current,;          ρ[qsize(id(q)) ≥ 1] ]/[                                     (q, ρ) →ρ [   vis[q.pop] ↦ M[𝑜𝑓𝑓𝑠(q), .., 𝑜𝑓𝑓𝑠(q)+size(q)],;                                                               M[𝑜𝑓𝑓𝑠(q)+size(q) × qsize[q], ..,;                                                 𝑜𝑓𝑓𝑠(q) + size(q) × (qsize[q]+1)] ↦ vis[push],] ]



  Data is only pushed to the stream node:


    [ q = 𝑠𝑡𝑟𝑒𝑎𝑚(𝑝𝑢𝑠ℎ, 𝑝𝑜𝑝, 𝑠𝑖𝑧𝑒),;        vis[push] ≠ emptyset,;     ∄ (q.pop, dst) ∈ current ]/[                                  (q, ρ) →ρ [              M[𝑜𝑓𝑓𝑠(q)+size(q) qsize[q], ..,;                                              𝑜𝑓𝑓𝑠(q) + size(q) (qsize[q]+1)] ↦ vis[push],;                                                                   qsize[q] ↦ qsize[q]+1 ] ]



  Data is popped from the stream but not pushed into it:


    [ q = 𝑠𝑡𝑟𝑒𝑎𝑚(𝑝𝑢𝑠ℎ, 𝑝𝑜𝑝, 𝑠𝑖𝑧𝑒),;   ∄ (src, q.push) ∈ S_state,;     ∃ (q.pop, dst) ∈ current;                              ]/[                                   (q, ρ) →ρ [ vis[q.pop] ↦ M[𝑜𝑓𝑓𝑠(q), .., 𝑜𝑓𝑓𝑠(q)+size(q)],;                                                                     qsize[q] ↦ qsize[q]-1 ] ]



Following the element processing step, we remove q from current, repeating the above step until current is empty. 



  §.§.§ Evaluating State Transitions


Once current is empty, we evaluate all outgoing state
transitions of the current state:
(state, next, cond, assigs) ∈ T. For
each transition, we resolve all symbols in cond and the right-hand sides of assigs using ρ, then evaluate arithmetic and boolean expressions using standard semantic rules, which we omit here.  If no condition evaluates to true, signal the completion of G to the caller and stop the evaluation of G:


    [                                           𝑐𝑢𝑟𝑟𝑒𝑛𝑡 = ∅,; ∄ (ρ[state], next, cond, assigs) ∈ T: (cond, ρ) → True;                                                        ]/ρ[ state ↦∅ ]



Otherwise, we choose an arbitrary transition for which cond → true and update ρ:
Set state to next, set curr to a copy of S[next]. For each left-hand side of an assignment
z_i, update sym[z_i] with the value of the
corresponding right-hand side v_i.  Data propagation then follows Section <ref>:


    [                                            current = ∅,; ∃ (ρ[state], next, cond, assigs) ∈ T: (cond, ρ) → True,;                               (assigs, ρ) → [z_i = v_i] ]/[                    ρ[         state ↦ next,;                       current ↦ copy(next),;                            sym[z_i] ↦ v_i ] ]





§ DATA-CENTRIC GRAPH TRANSFORMATIONS







Table <ref> lists the transformations used in the Performance Evaluation section of the paper.



The implementation of graph transformations is based on
algebraic graph rewriting <cit.>. Each transformation is
implemented by a rule p: L  R, which consists of the
two sub-graphs L and R, while r is a relation that maps the vertices and
edges of L to elements of the same kind in R. Moreover, a specific
matching of L in the SDFG G is represented by a relation
m: L → G. Applying the optimization on G produces the transformed
graph H, which can be constructed as part of the pushout
<H, m^', r^'> of p, G and m. A visualization of the above
method, also known as the single-pushout approach
<cit.>, is shown in Fig. <ref>.






§ POLYBENCH FLAGS



We compile Polybench code for all compilers using the following base flags: . Each compiler was also tested with different variants of flags, where we report the best-performing result in the paper. The flags are (variants separated by semicolons): 

	
  * gcc 8.2.0: ; Base flags
	
  * clang 6.0: Base flags
	
  * icc 18.0.3: Base flags; 
	
  * Polly (over clang 6.0): ;
 
	
  * Pluto 0.11.4:  and
 among best of ; ;
 ;
 
	
  * PPCG 0.8: Base flags


All variants were also tested with compile-time size specialization.




§ TRANSFORMATION EXAMPLE: REDUNDANT ARRAY REMOVAL


A simple transformation implemented in DaCe is , whose code can be found below.

The transformation removes a transient array that is used directly before another array, creating a copy, and not used anywhere else (making the copy redundant). This situation often happens after transformations and due to the strict nature of some language frontends (e.g., TensorFlow). The subgraph expression to match (lines 6–14) is a path graph of size two, connecting the two access nodes. The nodes have no restriction on their content (see lines 6–7). The function  is called on a matching subgraph for further programmatic checks (lines 16–58), and the function  (lines 60–78) applies the transformation using the SDFG builder API. The checks ensure that the array is indeed transient and not used in other instances of data access nodes. To avoid recomputing subsets (which may not be feasible to compute symbolically), if the transformation operates in strict mode, it only matches two arrays of the same shape (lines 51–56). The transformation then operates in a straightforward manner, renaming the memlets to point to the second (not removed) array (lines 66–70) and redirecting dataflow edges to that data access node (lines 73–74). Lastly, the  node is removed from the SDFG state (line 78).

In line 81, the transformation is registered with a global registry. This process enables users to import custom transformations and integrate them into DaCe easily.
However, because  is a strict transformation (i.e., can only improve performance), it is also hardcoded in a set of such transformations within the SDFG class. Strict transformations are applied automatically after processing a DaCe program, and other such transformations include  and .


class RedundantArray(pm.Transformation):
    """ Implements the redundant array removal transformation, 
        applied when a transient array is copied to and from (to 
        another array), but never used anywhere else. """

    _in_array = nodes.AccessNode('_')
    _out_array = nodes.AccessNode('_')

    @staticmethod
    def expressions():
        return [
            nxutil.node_path_graph(RedundantArray._in_array,
                                   RedundantArray._out_array),
        ]

    @staticmethod
    def can_be_applied(graph, candidate, expr_index, sdfg, 
                       strict=False):
        in_array =   
          graph.nodes()[candidate[RedundantArray._in_array]]
        out_array =           graph.nodes()[candidate[RedundantArray._out_array]]

        # Ensure out degree is one (only one target, out_array)
        if graph.out_degree(in_array) != 1:
            return False

        # Make sure that the candidate is a transient variable
        if not in_array.desc(sdfg).transient:
            return False

        # Make sure both arrays are using the same storage location
        if (in_array.desc(sdfg).storage !=  
              out_array.desc(sdfg).storage):
            return False

        # Find occurrences in this and other states
        occurrences = []
        for state in sdfg.nodes():
            occurrences.extend([
                n for n in state.nodes() 
                if isinstance(n, nodes.AccessNode)
                and n.desc(sdfg) == in_array.desc(sdfg)
            ])

        if len(occurrences) > 1:
            return False

        # If strict, only apply if arrays are of same shape 
        # (no need to modify memlet subset)
        if (strict and (len(in_array.desc(sdfg).shape) != 
                len(out_array.desc(sdfg).shape))
                or any(i != o for i, o in zip(
                    in_array.desc(sdfg).shape,
                    out_array.desc(sdfg).shape))):
            return False

        return True

    def apply(self, sdfg):
        graph = sdfg.nodes()[self.state_id]
        in_array = gnode(RedundantArray._in_array)
        out_array = gnode(RedundantArray._out_array)

        for e in graph.in_edges(in_array):
            # Modify all incoming edges to point to out_array
            path = graph.memlet_path(e)
            for pe in path:
                if pe.data.data == in_array.data:
                    pe.data.data = out_array.data

            # Redirect edge to out_array
            graph.remove_edge(e)
            graph.add_edge(e.src, e.src_conn, out_array, 
                           e.dst_conn, e.data)

        # Finally, remove in_array node
        graph.remove_node(in_array)

# Register transformation in global registry
pm.Transformation.register_pattern(RedundantArray)




§ GRAPH DATASETS


Table <ref> lists the graph datasets used in the BFS case study of the paper.





§ INDIRECT MEMORY ACCESS


Indirect memory access, i.e., , is an important characteristic
of sparse data structure and pointer jumping algorithms.
Indirection cannot be directly represented by a single memlet. As shown in Fig. <ref> (an excerpt of Sparse Matrix-Vector Multiplication), such statements are converted to a subgraph,
where the internal access is given exactly (
in the figure), and the indirect memory is given by a memlet
with one access (). The memory is then accessed dynamically
and copied using a tasklet.







