\section{Introduction}

The communication performance of applications running on High-Performance Computing (HPC) systems
depends on a variety of factors like the
capability and topology of the underlying communication system, the
required communication (patterns, frequencies, volumes, and
dependencies) between processes, and the
software and algorithms used to realize the communication.  If the
communication pattern is known, and if a hardware topology description
is given, it is natural to attempt to find a good mapping of the
application processes onto the hardware processors such that pairs of
processes that frequently communicate large amounts of data
become~located~closely.

Many important scientific computing applications involve stencil
computations.  For example, stencil computations are used for climate
and ocean modeling \cite{haralick1992computer}, in computational
electromagnetic codes
\cite{taflove2005computational,DBLP:journals/jcphy/UstyugovPKN09}, for
image-processing \cite{sawdey1995design}, in Jacobi or multigrid
solvers~\cite{DBLP:conf/ipps/RenganarayananHDR07}, for earthquake
simulations \cite{DBLP:conf/sc/ChristenSC12} or in general in
simulations systems such as
OpenLB~\cite{DBLP:conf/europar/FietzKSSH12}.  In most cases, elements
of a $d$-dimensional matrix are repeatedly updated using the values of
fixed \emph{stencil} pattern of neighboring elements. When run on a
parallel computer, this yields communication patterns that are very
regular and more or less symmetric depending on the organization of
the processors. More precisely, each processing element exchanges data
repeatedly with a small set of neighboring processing elements, and
all processing element neighborhoods have the same structure. In this
situation, in each exchange step, all processes communicate with
other processes and all follow the same pattern determined by the
computational stencil and the organization of the processes.
\begin{figure}[t]
	\center
	\includegraphics[trim = 0 0 0 0, clip, height=4.5cm]{./images/motivational_example}
	\caption{Motivational example: Given a set of compute nodes
          with possibly different numbers of processes per node, and a
          computational grid, find a mapping of the processes on the nodes
          to the grid, \textrm{s.t.}\xspace the number of communication edges between
          compute nodes is minimized.}
	\label{fig:motivation}
\end{figure}

The Message Passing Interface (MPI)~\cite{MPI-3.1} supports complex
communication patterns by providing functions to specify virtual
process topologies by process neighborhoods.  Using for instance
Cartesian topologies the user can refer to processes by rank or by
coordinate vectors. Moreover, MPI supports neighborhood collective
operations such as \texttt{MPI\_\-Neighbor\_\-alltoall}\xspace which make it possible for the
MPI library to exploit (regular) communication patterns to provide
more efficient data exchange operations.  MPI also
defines functionality to reorder processes in order to optimize the
communication performance, however, at the moment most MPI libraries
do not actually perform such remapping (in the general case). \\

\noindent \textbf{Contribution.} We make the following contributions:
\begin{itemize}
\item We show that the general Cartesian mapping problem under stencil
  patterns is NP-hard.
  This result motivates
  our work on heuristic algorithms for the problem.
\item We present new algorithms for the process-mapping problem for
  stencil patterns which in contrast to previous solutions are also
  applicable to cases where 1) the number of \texttt{MPI}\xspace processes per node is
  different and 2) where the number of processes is not factorizable
  or divisible by the number of processes per node, and 3) consider
  the case of arbitrary stencil patterns, not only the nearest-neighbor stencils
  implied by the \texttt{MPI}\xspace specification.
\item We perform an extensive experimental evaluation and benchmark
  the time needed for an \texttt{MPI\_\-Neighbor\_\-alltoall}\xspace operation. The results
  show that our algorithms significantly outperform previous solutions in terms of
  communication performance as well as initialization time.
 For example, our algorithms are up to two orders of magnitude faster in running time than the (sequential) 
	high-quality general graph mapping tool Vienna Mapping (\texttt{VieM}\xspace), while obtaining similar results in
	communication performance. Moreover, our algorithms are up to three times faster than other 
	Cartesian grid mapping algorithms, and achieving significantly better mapping quality.  \\
\end{itemize} 

\noindent \textbf{Organization.} The rest of the paper is organized as follows.  We start by
introducing the process-to-node mapping problem in
Section~\ref{sec:preliminaries} and discuss related work in
Section~\ref{sec:relatedwork}.  In Section~\ref{sec:hardness}, we look
at the mapping problem for Cartesian graphs, for which we show that
mapping problem for specific graph types is NP-hard.  For that reason,
we introduce three different, efficient algorithms to solve the
process-to-node mapping problem for Cartesian graphs in
Section~\ref{sec:neighborhoodaware}.  We present the results of an
extensive experimental evaluation of our novel algorithms in
Section~\ref{sec:experiments} before we conclude in Section~\ref{sec:conclusion}.

\section{Notation and Problem Formulation}
\label{sec:preliminaries}

We are considering the traditional setup of HPC system architectures,
where several compute nodes are interconnected via a high-speed
network.  The compute nodes usually comprise multiple processor-cores,
often on two or more CPU sockets. In order to solve a computational
problem on such systems, data needs to be exchanged among the
distributed processes.  We call communication between processes
residing on different compute nodes \emph{inter-node communication}
and communication between processes residing on the same compute node
\emph{intra-node communication}.  We follow the common assumption that
intra-node communication (on a compute node) is (much) faster than
inter-node communication with higher cumulated bandwidth. We
assume homogeneous communication performance between the computation nodes,
and also within the nodes~\cite{Gropp2019,Niethammer2019}.

We denote by $N$ the \emph{number of
  (compute) nodes} allocated for the application to run.  Let $p$ be
the total \emph{number of processes} of an application and $n_i$ with
$i \in \{0, 1, \dots, N-1\}$ the \emph{number of processes per node
  $i$}, \textrm{i.e.}\xspace, $\sum_{i=0}^{N-1} n_i = p$. If all the nodes have the
same number of processes (homogeneous node sizes), we denote by $n$
the number of processes on each node, \textrm{i.e.}\xspace, $p = N n$.
     
We assume that the processes are organized in a $d$-dimensional
Cartesian grid with dimension sizes $\mathcal{D}=[d_0, \dots, d_{d-1}]$, and thus,
the \emph{size of a grid} is the number of processes it comprises,
$p=\prod_{i=0}^{d-1}d_i$.  Each process with rank~$r$,
$0 \leq r < p$, is associated with a vector
$\vec{r}=[r_0, \dots, r_{d-1}]$, where $0 \leq r_i < d_i$ for
$i \in \{0, 1, \dots, d-1\}$, uniquely determining the position of
the process in the grid. W.l.o.g., processes are assigned in row-major
order to the grid.

\paragraph*{Target Stencils}

We now define three different stencils which will be used in the
remainder of the article. A $2$-d example of the considered stencils is depicted in
Figure~\ref{fig:stencils}.
To that end, we consider a
$k$-\emph{neighborhood} of a process to be a set of communication
targets which can be described as a list of relative coordinates
$\mathcal{S}=\{\vec{R_0}, \vec{R_1}, \dots, \vec{R_{k-1}}\}$.
Every~$\vec{R_i}=[R_{i,0}, R_{i,1}, \dots, R_{i,d-1}]$ with
$0 \leq i < k$ describes the relative offset along the dimensions
to the target process.  Let $\mathds{1}_i$ be a vector with the only
non-zero component being one at index $i$.
Then, we define the following stencils:
\begin{enumerate}[label=(\alph*)]
\item \emph{nearest neighbor stencil}:
  $\mathcal{S}=\{\mathds{1}_i, -\mathds{1}_i \mid  0 \leq i <
  d\}$,
\item \emph{component stencil}: $\mathcal{S}=\{\mathds{1}_i, -\mathds{1}_i \mid  0 \leq i <
  d-1\}$, and
\item \emph{nearest neighbor with hops}: $\mathcal{S}=\{\mathds{1}_i, - \mathds{1}_i \mid 0~\leq~i~<~d\} \cup \{a\mathds{1}_0, -a\mathds{1}_0 \mid \forall a \in \{2, 3\}\}$.
\end{enumerate}
\begin{figure}[t]
	\begin{subfigure}{0.20\textwidth}
		\center
	\includegraphics[trim=0 110 250 60, clip, scale=0.55]{./images/nearestneighbor}
	\caption{nearest neighbor stencil}
	\label{fig:nearestneighbor}
	\vfill
	\center
	\includegraphics[trim=0 110 250 60, clip, scale=0.55]{./images/component}
	\caption{component stencil}
	\label{fig:component}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.30\textwidth}
		\center
                \vspace*{.5cm}
	\includegraphics[trim=0 35 250 20, clip, scale=0.55]{./images/nearestneighborhops}
	\caption{nearest neighbor with hops stencil}
	\label{fig:nearestneighborhops}
	\end{subfigure}
	\caption{Examples of the three $2$-d stencils.}
	\label{fig:stencils}
\end{figure}


\paragraph*{Optimization Problem}
By defining the $k$-neighborhood communication neighbors of each process in the
Cartesian process grid with dimension sizes $\mathcal{D}$,
we induce a Cartesian communication graph $\mathnormal{C}=(V, E)$ (Cartesian graph)
where $V$ is the vertex set representing the processes, \textrm{i.e.}\xspace, $|V|=p$
and $E$ is the set of communication edges between the processes. We
assume unit edge-weight and sparse communication, \textrm{i.e.}\xspace, the number of
communication neighbors is much smaller than the total number of
processes~($k \ll p$).

Let $\sigma: V \times V \rightarrow \{0,1\}$ be a cost function that
determines whether the communication between two processes $v_i$ and
$v_j$ involves two different compute nodes, \textrm{i.e.}\xspace, inter-node
communication is required.
Let $\mathcal{N}$ be the set of compute nodes and let
$M : V \rightarrow \mathcal{N}$ be a function that maps a process
$u \in V$ to exactly one compute node $\eta \in \mathcal{N}$.
For all $u, v \in E$, let~$\sigma(u,v)=0$ if $M(u) = M(v)$ and $\sigma(u,v)=1$ otherwise.
The total cost (amount) of inter-node communication operations is defined as
$J_\text{sum} := \sum_{ (u, v) \in E} \sigma(u, v)$.
We define the bottleneck node as the node with the largest number of
outgoing communication edges, \textrm{i.e.}\xspace,
$N_b := \operatorname*{argmax}_{\eta \in \mathcal{N}}\{\sum_{ (u,v) \in E}
\sigma(u, v) \mid M(u)=\eta\}$.  Let~$\mathcal{B} := \{u \mid \forall u \in V :
M(u)=N_b\}$ be the set of vertices assigned to the bottleneck node $N_b$.
Then the cost of this bottleneck node is~$J_\text{max} :=\sum_{ u\in
\mathcal{B}} \sum_{ (u,v) \in E} \sigma(u,v)$.

Our objective is to find a mapping function $M$ of processes to nodes
that minimizes $J_\text{sum}$. We use $J_\text{max}$ to distinguish
cases with similar values of $J_\text{sum}$, especially in the
experimental evaluation. Note that the original allocation
($N \times n$) given by the scheduler needs to be respected, \textrm{i.e.}\xspace, for
each node $N_i$ it must hold that $|\{u \in V \mid M(u)=N_i\}| = n_i$.

\section{Related Work}
\label{sec:relatedwork}

There has been an immense amount of research on partitioning and process
mapping -- we refer
to~\cite{GPOverviewBook,SPPGPOverviewPaper,DBLP:reference/bdt/0003S19} for
extensive material.  The problem of process reordering for different topologies
has been an
active field of research since the beginning of \texttt{MPI}\xspace
\cite{Brandfass2013,Traeff2002,Mercier2011,Hatazaki1998,YuChungMoreira06}.
Many reordering algorithms take an arbitrary unstructured graph as
input topology, making it difficult to perform efficient, scalable
mappings on a structured grid, where communication is implicitly implied
through the grid structure. In this paper we aim to exploit both stencil and
grid structure, \textrm{i.e.}\xspace, aim for specialized algorithms.

Gropp~\cite{Gropp2019} pointed out that many \texttt{MPI}\xspace implementations have not
implemented the \texttt{MPI\_\-Cart\_\-create}\xspace reordering method. As a response, he
proposed an algorithm (\texttt{Nodecart}\xspace) for homogeneous node sizes
$n$, based on the prime factorization of $n$. \texttt{Nodecart}\xspace decomposes
the dimensions into a grid spanning the nodes and a grid describing
the layout of the processes within a node. From this decomposition,
every process can calculate its new coordinate $\vec{r}$ from which it
can obtain its new rank. As a result, Gropp was able to show
significant improvements in the time needed for a nearest neighbor
message exchange in comparison to a blocked mapping of processes to nodes.
\texttt{Nodecart}\xspace was specifically designed for the implied nearest neighbor stencil of
Cartesian communicators in \texttt{MPI}\xspace.

Niethammer and Rabenseifner \cite{Niethammer2019} used a different approach to
assign processes to nodes. They point out that the \texttt{MPI\_\-Dims\_\-create}\xspace routine only
considers the total number of processes $p$ in an application from which it
finds a grid decomposition where the dimension sizes are as close as possible.
This algorithm can lead to bad domain decompositions in terms of inter-domain
communication, if the underlying data mesh is not shaped cubically.  Thus, they
propose to solve the task of grid dimension creation and process mapping
simultaneously. This is done by finding a factorization of the number nodes $N$, with the
aim to minimize the weighted communication over the domain boundaries (the
weights represent the expansion of the application mesh and a user defined
communication cost factor). Their algorithm can be extended to handle
hierarchical systems, where the weighted inter-domain communication is
minimized at each level, but it requires symmetric hierarchies.  With this
approach, they can achieve significant performance gains for a
nearest neighbor message exchange in comparison to the blocked assignment of 
processes to nodes.

Schulz et al. \cite{Schulz2017, Schulz2017a}
developed an algorithm (VieM, Vienna Mapping) for general process mapping with the objective
of minimizing the total weighted communication.
Their algorithm takes as input an unstructured communication graph
and maps it onto a hierarchical hardware graph. This is done in a 
recursive manner with perfectly balanced graph partitioning techniques
and randomized local search for improving found solutions. Even though
the approach is costly in terms of runtime and memory,  the communication cost in comparison
to the state-of-the-art has been significantly~reduced.

\section{NP-Hardness of Cartesian Mapping Problem}
\label{sec:hardness}
In general, graph embedding problems are NP-hard, as shown by
\cite{Bokhari1981}. However, the structure of the Cartesian graph
mapping problem induced by a $k$-neighborhood pattern could make the
problem easier to solve, in terms of NP-hardness or complexity of
approximation algorithms. We propose the following formal definition of the
Cartesian partitioning problem.

\begin{definition}
	Let $\mathcal{C}=(V, E)$ be a Cartesian graph with dimension sizes $\mathcal{D}$ and 
	$k$-neighborhood $\mathcal{S}$, as defined in
	Section~\ref{sec:preliminaries} and let $\mathcal{N}$ be a set of
	partition sizes (number of cores per compute node) $\mathcal{N} = n_0,
	\dots, n_{N-1}$, \textrm{s.t.}\xspace, $\sum_{i=0}^{N-1} n_i = |V|$.  Let $M: V
	\rightarrow \mathcal{N}$ be a mapping function that assigns each vertex
	$v \in V$ of the Cartesian graph to a distinct partition~$n \in
	\mathcal{N}$.
  The \texttt{GRID-\-PARTITION}\xspace problem answers the question whether there exists
	a mapping $M$ such that $J_\text{sum} \le Q$.
  
\end{definition}

\begin{definition}
  The \texttt{3-WAY-\-PARTITION}\xspace problem consists of dividing a
  multi-set of integers into three subsets, such that the sum of each
  subset is equal.  Formally, given a multi-set~$\mathcal{I}$ of
  integers, we ask whether~$\mathcal{I}$ can be partitioned into $3$
  disjoint sub-sets $I_{1, 2, 3}$, where
  $\sum_{x \in I_{1}} x = \sum_{y \in I_2} y = \sum_{z \in I_3} z$
	and~$I_1 \cup I_2 \cup I_3 = \mathcal{I}$ holds.
\end{definition}
It is well-known that the \texttt{3-WAY-\-PARTITION}\xspace problem is
NP-complete \cite{Korf2009}.
Now, we show that the \texttt{GRID-\-PARTITION}\xspace problem is already
NP-hard for two dimensions and a simple, one-dimensional
component stencil.
To that end, we reduce \texttt{3-WAY-\-PARTITION}\xspace to
\texttt{GRID-\-PARTITION}\xspace which leads to the following theorem.

\begin{theorem}
The \texttt{GRID-\-PARTITION}\xspace problem is NP-hard, when restricted to two
dimensions and a one-dimensional component stencil
	$\mathcal{S}=\{-\mathds{1}_1, \mathds{1}_1\}$.
\end{theorem}

\begin{figure}[t]
	\center
	\includegraphics[trim=0 0 0 0, clip, height=2.5cm]{./images/np_hardness_example_bw}
	\caption{Example of \texttt{3-WAY-\-PARTITION}\xspace to \texttt{GRID-\-PARTITION}\xspace transformation. Given 
	$\mathcal{I'}=\{6,3,3,2,2,2\}$, create a grid with $\mathcal{D}=[6,2]$ for the Cartesian graph
	and find a mapping \textrm{s.t.}\xspace $J_\text{sum} \leq 2|\mathcal{I'}|-6$.}
	\label{fig:hardness}
\end{figure}

\begin{proof}
	Given an arbitrary instance $\mathcal{I'}$ of \texttt{3-WAY-\-PARTITION}\xspace,
        we construct an instance of
        \texttt{GRID-\-PARTITION}\xspace with a Cartesian graph, composed of $\mathcal{S}$ and~$\mathcal{D}$,
	and $\mathcal{N}, Q$
        as follows: $\mathcal{S} = \{-\mathds{1}_1, \mathds{1}_1\}$, $\mathcal{D}=[3,
          \frac{\sum_{x \in \mathcal{I'}} x}{3}]$, $\mathcal{N} =
	\{x \mid x \in \mathcal{I'}\}$, $Q = 2 |\mathcal{I'}| - 6$.

An optimal mapping of a two-dimensional
\texttt{GRID-\-PARTITION}\xspace problem with the component stencil always
traverses the vertices in the grid along the communicating dimension
given by the stencil, assigning them to a partition until it is
full. Thus, for $|\mathcal{N}| \geq 3$, each partition has at most two
outgoing communication edges (the first vertex assigned to the
partition and the last), except partitions at the border of the grid
which have one outgoing communication edge, \textrm{i.e.}\xspace, we can say \textrm{w.l.o.g}\xspace
that $Q= 2|\mathcal{N}| - 6$.

A \emph{yes} instance of \texttt{3-WAY-\-PARTITION}\xspace now corresponds to a
\emph{yes} instance of \texttt{GRID-\-PARTITION}\xspace, since we can assign every
vertex in the first, second and third column to the partitions that correlate
to the values in of $I_1$, $I_2$, $I_3$, respectively.
\end{proof}

An example for an instance with $\mathcal{I'}=\{6,3,3,2,2,2\}$ is shown in
Figure~\ref{fig:hardness}. By encoding solutions with the first and last vertex of each
partition, one can easily show that two-dimensional \texttt{GRID-\-PARTITION}\xspace
with the component stencil is NP-complete (that is, is also in
NP). With some adjustments, one can also show that the problem is
NP-complete if we allow periodicity along the dimension of
communication. The more interesting case is for which fixed stencils
the problem remains NP-complete.

\section{Rank Reordering Algorithms for $k$-Neighborhoods}
\label{sec:neighborhoodaware}

In this section, we propose three algorithms for a $k$-neighborhood
aware process reordering for Cartesian grids. The goal is to find
reordering schemes that are a) fully distributed, that is, each process
can compute its new rank independently of the other processes based on
the input alone (grid and stencil), and b) efficient, that is, not
polynomially dependent on $p$, preferably dependent only on the size
of the (sparse, compared to $p$) stencil and the number of dimensions
(we will tolerate polylog $p$).

\subsection{Hyperplane Algorithm}
The \texttt{Hyperplane}\xspace algorithm is a variation of recursive bisection. 
The main idea consists in recursively finding a split of a suitable dimension
$d_i$ of the Cartesian grid $g$ into $d_i', d_i''$ \textrm{s.t.}\xspace $d_i= d_i' + d_i''$.
This induces two grids $g'$ and $g''$, where the $i$th dimension size of $g'$
is $d_i'$ and of $g''$ is $d_i''$. The split is chosen \textrm{s.t.}\xspace the sizes of the two
new grids is a multiple of $n$, \textrm{i.e.}\xspace, $n\mid |g'|$ and $n\mid|g''|$.
Note that if we have heterogeneous nodes, one can use the mean, minimum
or maximum of the node sizes as an input for the algorithm.
This produces $N$ grids each of
size $n$ which can be mapped to the nodes. The cuts should be chosen \textrm{s.t.}\xspace the
minimal possible amount of communication between the grids is induced, since
those correspond to
inter-node communication. For that purpose, we calculate how parallel each
vector~$\vec{R} \in \mathcal{S}$ in the stencil is to a grid dimension $j$
using the cosine.
\begin{equation}\label{angle}
	\cos(\alpha_{\vec{R}, \vec{e_j}}) = \frac{\vec{R}\vec{e_j}}{||\vec{R}||||\vec{e_j}||} \in [-1, 1].
\end{equation}
Here, $\vec{e_j}$ is the unit vector along dimension~$j$ with~$0 \leq j < d$,~$\alpha_{\vec{R_i}, \vec{e_j}}$ 
is the angle between the relative coordinate vector~$\vec{R_i}$ and dimension $j$'s unit vector $\vec{e_j}$.
In order to have a monotonic increasing function, we 
square each of the values in Equation~\eqref{angle} and sum them over all relative coordinate
vectors $\vec{R} \in \mathcal{S}$, giving us the following list
\begin{equation}\label{cosines}
	\left[\sum_{i=0}^{k-1} \cos^2(\alpha_{\vec{R_i}, \vec{e_0}}),\ 
	\dots,\
	\sum_{i=0}^{k-1} \cos^2(\alpha_{\vec{R_i}, \vec{e_{d-1}}})\right].
\end{equation}
The dimension with the minimal value in Equation~\eqref{cosines} is the most orthogonal
to all $\vec{R} \in \mathcal{S}$, thus, we try to partition the grid alongside of it. 
Ties are broken by size, \textrm{i.e.}\xspace, we want to partition along the bigger dimension.
By sorting the dimensions according to their value in Equation~\eqref{cosines},
we define a preferred dimension order along which we try~to~perform~the~cuts.

The pseudo-code can be found in Algorithm~\ref{hyperplanealgo}. The input consists
of the dimension sizes~$\mathcal{D}$, the~$k$-neighborhood $\mathcal{S}$,
the number of processes per node~$n$, the rank~$r$ of the calling process and
it outputs the new position~$\vec{r}_\text{new}$ of the calling rank on the
Cartesian grid.

In each recursive step, we check if the grid is smaller than~$2 n$. 
We do this, for it is not needed to find an explicit cut for this grid size,
rather, we can directly calculate the new coordinates with the preferred
dimension order. This avoids bad splitting along very skewed grids, \textrm{e.g.}\xspace, a
nearest neighbor stencil on a two-dimensional grid with dimensions~$[2, n]$
where~$n$ is large and odd. Instead of being forced to cut along the first
dimension of size~$2$ to obtain two~$[1, n]$ partitions, we obtain two
partitions, each with~$3$~outgoing~communication~edges.

Otherwise, the algorithm finds the best possible split of the current
grid into two new grids. For that purpose, it traverses the
current dimensions~$\mathcal{D}$ sorted in increasing order of the values
defined in Equation~\eqref{cosines} (sorting in each recursive step is necessary,
because of changing dimension
sizes) and tries to position the splitting hyperplane in the current dimension~$d_i$.  The
hyperplane is initially placed at the center~$\frac{d_i}{2}$ of the candidate
dimension~$d_i$. If the initial split is not
suitable, the position of the hyperplane is incremented/decremented, respectively, \textrm{s.t.}\xspace
the position is as close as possible to the original grid's border in
an effort to reduce~$J_\text{max}$. If it cannot find a suitable split along the 
candidate dimension, it will proceed to the next, until it finds a split.
With the following proof, we show that it is always possible to find such a split.
\begin{theorem}\label{hyperplaneproof}
	Let $C, n \in \mathbb{N}$ and $C \geq 2$, let the dimension sizes
	be given by $\mathcal{D}=[d_0, \dots, d_{d-1}]$ and $\forall d
	\in \mathcal{D} : d \in \mathbb{N}$, \textrm{s.t.}\xspace $C  n =
	\prod_{i=0}^{d-1} d_i$. Then, it is always possible to split a dimension
	$d'$ \textrm{s.t.}\xspace the two induced grids are of size which is a multiple of
	$n$.
\end{theorem}

\begin{proof}
	Let $F(x) =
	\{ f_{1}, \dots, f_{l} \} $ be the multi-set of all prime factors of $x$.
	Then,
	\begin{equation}\label{prime_factor_product}
		F(\prod_{i=0}^{d-1} d_i) = \prod_{i=0}^{d-1} F(d_i) = F(C) F(p) = F(C p).
	\end{equation}
	Since $C \geq 2$ it must hold that, $\exists d' \in \mathcal{D}$, \textrm{s.t.}\xspace
	$F(d') \cap F(C) \neq \emptyset$
\end{proof}
In Line~$5$ of Algorithm~\ref{hyperplanealgo}, the subroutine \texttt{find\_split}
returns the index $i$ of the dimension to be split, and split sizes $d'$ and $d''$.
When a suitable split is found, two new grids $g'$, $g''$ are created where the $i$th dimension
size is replaced with $d'$ and $d''$ for $g'$ and $g''$, respectively. If the calling
rank is located on the left-hand side of the split, it will call \texttt{Hyperplane}\xspace with
$g'$ ($LHS$) as input and else, it calls \texttt{Hyperplane}\xspace with $g''$ ($RHS$) as new~input.

The number
of recursions executed by the \texttt{Hyperplane}\xspace algorithm is logarithmic
in the number of compute nodes~$N$, although the two grids per
recursion step can be very imbalanced in terms of size.
\begin{theorem}
  \label{proof:hyperplane_runtime}
	Let $g$ be a grid with dimension sizes $\mathcal{D}=[d_0, \dots,\
	d_{d-1}]$ and $\forall d \in \mathcal{D} : d \in \mathbb{N}$,
	and $\prod_{i=0}^{d-1} d_i = C p$ for some $C \in \mathbb{N}$.
	Then, the \texttt{Hyperplane}\xspace algorithm will always partition $g$ into two
	grids $g'$ and $g''$ \textrm{s.t.}\xspace $\frac{1}{2} \leq \frac{|g'|}{|g''|} \leq
	1$.
\end{theorem}

\begin{proof}
	Let $F(x)$ be as defined in the proof of Theorem~\ref{hyperplaneproof}.
	Let $d'$ be the candidate dimension of the algorithm with $F(d') \cap
	F(C) \neq \emptyset$. Let $F(d') \cap F(C) = \{d'_1, \dots,
	d'_m\}$ be ordered in ascending order.
	Then, the algorithm will surely find a suitable split at $d'_1$
	with $\left\lfloor \frac{d'_{1}}{2} \right\rfloor$ and $\left\lceil \frac{d'_{1}}{2} \right\rceil$. If $d'_{1}=2$, the resulting split will yield two partitions of exactly the same size.
	
	If $d'_{1} \geq 3$, then a split yields two partitions with a bigger or
	equal difference than a split of any other prime factor $d'_k \in F(d') \cap F(C)$.
	With $\left\lfloor
	\frac{d'_{1}}
	{2} \right\rfloor = \frac{d'_{1}-1}{2}$ and $\left\lceil \frac{d'_{1}}{2} \right\rceil = \frac{d'_{1}+1}{2}$,
	\begin{equation}
		\begin{aligned}
			\frac{|g'|}{|g''|} =
			\frac{\frac{d'_{1}-1}{2} \prod_{j=2}^m d'_j \prod_{d_i \neq d'} d_i}
			{\frac{d'_{1}+1}{2} \prod_{j=2}^m d'_j \prod_{d_i \neq d'} d_i} &\leq
			\frac{\frac{d'_{k}-1}{2} \prod_{j \neq k }^m d'_j \prod_{d_i \neq d'} d_i}
			{\frac{d'_{k}+1}{2} \prod_{j\neq k}^m d'_j \prod_{d_i \neq d'} d_i} \\
			\implies 
			\frac{d'_{1}-1}{d'_{1}+1} &\leq \frac{d'_{k}-1}{d'_{k}+1}.\\
		\end{aligned}
	\end{equation}	
	This always holds, since this is a strictly monotonic increasing function that converges to $1$ for growing $d'_k$ values.
		To see this suppose $x,y \geq 0$
		\begin{equation}
			\begin{aligned}
				\frac{x-1}{x+1} &\leq \frac{y-1}{y+1}\\
				(x-1)(y+1) &\leq (y-1)(x+1)\\
				xy +x -y -1 &\leq xy + y -x -1\\
				x &\leq y.
			\end{aligned}
		\end{equation}
	Note that if the Algorithm~\ref{hyperplanealgo} first positions the hyperplane at $\left\lfloor \frac{d'}{2} \right\rfloor$ it will eventually find a
	suitable split, latest at $\left\lfloor \frac{d'_{1}}{2} \right\rfloor \prod_{j=2}^{m} d'_j$.
	\begin{equation}
		\begin{aligned}
			\left\lfloor \frac{d'_1}{2} \right\rfloor \prod_{j=1}^m d'_j &\leq \left\lfloor \frac{d'}{2} \right\rfloor\\
			\frac{d'_1 - 1}{2} \prod_{j=1}^m d'_j &\leq \frac{d' - 1}{2}\\
			d' - \prod_{j=2}^m d'_j &\leq d' -1\\
			1 &\leq  \prod_{j=2}^m d'_j
		\end{aligned}
	\end{equation}
	
	We can bound the ratio of the two grid sizes $|g'|$ and $|g''|$ from below, since $d'_1 \geq 3$
	\begin{equation}
		\begin{aligned}
			\frac{|g'|}{|g''|} =  
			\frac{\frac{d'_{1}-1}{2} \prod_{j=2}^m d'_j \prod_{d_i \neq d'} d_i}
			{\frac{d'_{1}+1}{2} \prod_{j=2}^m d'_j \prod_{d_i \neq d'} d_i} &\geq \frac{d'_1 - 1}{d'_1 +1}
			\geq \frac{3-1}{3+1} &\geq \frac{1}{2}.
		\end{aligned}
	\end{equation}
\end{proof}

It follows that the running time 
of \texttt{Hyperplane}\xspace is bounded by~$\mathcal{O}\left( \log N  \sum_{i=0}^{d-1} d_i  \right)$.

\begin{algorithm}[t!]
	\SetKw{break}{break}
	\SetKw{return}{return}
	\SetAlgoLined
	\SetKwFunction{hyperplane}{hyperplane}
	\SetKwFunction{newcoordinate}{new\_coordinate}
	\SetKwFunction{sort}{sorted}
	\SetKwFunction{adapt}{push\_to\_boundary}
	\SetKwFunction{findsplit}{find\_split}
	\SetKwFunction{inducedgrids}{induced\_grids}
	\SetKwFunction{calcremain}{n\_remain\_parts}
	\KwIn{Grid with dimension sizes~$\mathcal{D}=[d_0, \dots, d_{d-1}]$,
	set of relative vectors~$\mathcal{S}$ describing the $k$-neighborhood,
	number of processes per node~$n$ and rank $r$ of calling process.}
	\KwResult{New coordinate $\vec{r}_\text{new}$ of calling rank.}
	\Indm\renewcommand{\nl}{\let\nl\oldnl}\hyperplane{$\mathcal{D}, \mathcal{S}, n, r, \vec{r}_\text{new}$}\\
	\Indp
	\eIf{$\prod_{d \in \mathcal{D}} d \leq 2 n$}{
		$\vec{r}_\text{new} \gets \newcoordinate{$\mathcal{D}, \mathcal{S}, n, r$}$\\
		\return
	}{
		$i, d', d'' \gets \findsplit{$\mathcal{D}, \mathcal{S}, n$}$\\
		$LHS, RHS~\gets~\inducedgrids{$\mathcal{D}, i, d',d''$}$\\
		\eIf{$r \in LHS$}{
			$\hyperplane{$LHS, \mathcal{S}, n, r_\text{new}$}$\\
		}{
			$\hyperplane{$RHS, \mathcal{S}, n, r_\text{new}$}$\\
		}
	}
	\caption{\texttt{Hyperplane}\xspace Algorithm}
	\label{hyperplanealgo}
\end{algorithm}
\begin{figure}[t]
	\begin{subfigure}{0.15\textwidth}
	\includegraphics[trim = 100 70 100 45, clip, scale=0.45]{./images/cartesian_paper_5x4}
	\caption{}
	\end{subfigure}
	\begin{subfigure}{0.15\textwidth}
	\includegraphics[trim = 100 70 100 45, clip, scale=0.45]{./images/cartesian_paper_5x4_1}
	\caption{}
	\end{subfigure}
	\begin{subfigure}{0.15\textwidth}
	\includegraphics[trim = 100 70 100 45, clip, scale=0.45]{./images/cartesian_paper_5x4_2}
	\caption{}
	\end{subfigure}
	\caption{\texttt{Hyperplane}\xspace algorithm example on $5 \times 4$ grid with a 
	nearest neighbor stencil and $N=5$, $n=4$:
	(a) The first split is along the largest dimension;
	(b) Two new splits found on both grids; (c) Final split of the grid.}
	\label{hyperplanefig}
\end{figure}

\subsection{$k$-$d$ Tree Algorithm}
\label{sec:kdtree}

Similar to the \texttt{Hyperplane}\xspace algorithm, but inspired by the $k$-$d$ tree data
structure, the $k$-$d$~\texttt{tree}\xspace algorithm works also by recursively splitting
the grid along the dimensions. The main difference to the \texttt{Hyperplane}\xspace algorithm
is that the recursion continues until there is only one vertex left in the Cartesian
grid. The advantage of this approach is that the algorithm is no
longer bound to the number of processes per node $n$. In fact, it is oblivious
to it and only tries to find \emph{dense} mappings, \textrm{i.e.}\xspace, localize communicating
vertices.

Instead of partitioning the grid dimensions in a round-robin manner, like in
the $k$-$d$ tree data-structure, we choose the largest dimension, weighted by
the inverse amount of communication that is performed across it, so that we avoid
splitting a dimension along which is communicated intensively. 
To be more precise, we define the amount of communication across a dimension
$j$ to be 
$
	f_j := |\{\vec{R} \mid \forall \vec{R} \in \mathcal{S} : R_{j} \neq
0\}|,
$
where $R_{j}$ is the $j$th component of $\vec{R}$, then we want to find the
dimension index $i := \operatorname*{argmax}_{0 \leq i < |\mathcal{D}|} \frac{d_i}{f_i}$
over the current grid dimensions sizes $\mathcal{D}$ and split equally along $d_i$.
All ranks $r$ with $r \leq \frac{1}{2}\prod_{j=0}^{d-1} d_j$ are assigned to
the left-hand side, and the others to the right-hand side of the split.

The pseudo-code
is given in Algorithm~\ref{algo:kdtree}. If there is only one vertex left in
the calling grid, we enter
the base-case of the recursive function, in which the vector $\vec{r}_\text{new}$ already
holds the coordinates of the calling process in the grid. Otherwise, 
we split the best suited dimension equally and recurse accordingly. 
\begin{algorithm}[t!]
	\SetKw{break}{break}
	\SetKw{return}{return}
	\SetAlgoLined
	\SetKwFunction{findsplit}{find\_split\_index}
	\SetKwFunction{kdtreecoord}{kd\_tree}
	\SetKwFunction{newcoordinate}{new\_coordinate}
	\KwIn{Dimension sizes $\mathcal{D}$, set of relative vectors
	describing the~$k$-neighborhood $\mathcal{S}$ and 
	rank~$r$ of calling process.}
	\KwResult{New coordinate $\vec{r}_\text{new}$ of calling rank}
	\Indm\renewcommand{\nl}{\let\nl\oldnl}\kdtreecoord{$\mathcal{D}, \mathcal{S}, r, \vec{r}_\text{new}$}\\
	\Indp
	\eIf{$\prod_{d \in \mathcal{D}} d = 1$}{ \return\\ }{
		$k \gets \findsplit{$\mathcal{D}, \mathcal{S}$}$\\
		\eIf{$r \leq \frac{\prod_{d \in \mathcal{D}}d}{2}$}{
			$d_k \gets \lfloor \frac{d}{2} \rfloor$\\
			$\kdtreecoord{$\mathcal{D}, \mathcal{S}, r, \vec{r}_\text{new}$}$\\
		}{
			$d_k \gets \lceil \frac{d}{2} \rceil$\\
			$\vec{r}_\text{new}[k] \gets \vec{r}_\text{new} + d_k$\\
			$\kdtreecoord{$\mathcal{D}, \mathcal{S}, r, \vec{r}_\text{new}$}$\\
		}
	}
	\caption{$k$-$d$~\texttt{tree}\xspace Algorithm}
	\label{algo:kdtree}
\end{algorithm}

For the run-time analysis, it is not difficult to see that depth of the 
recursion tree is $\log_2 p$, where $p$ is the total number of processes
in the Cartesian grid. At each step of the recursion, we have to find the
best dimension for the split. Using a priority queue, this can be achieved
in $\mathcal{O}\left(\log d\right)$ steps. Since the remaining computations
are constant, the runtime for calculating the new reordering is
$\mathcal{O}\left(\log p \log d \right)$.

\subsection{Stencil Strips Algorithm}
\label{sec:stencilstrips}
\begin{algorithm}[b!]
	\SetKw{break}{break}
	\SetAlgoLined
	\SetKwFunction{stencilstrips}{stencil\_strips}
	\SetKwFunction{distortionfactors}{get\_distortion\_factors}
	\SetKwFunction{striplengths}{get\_strip\_lengths}
	\SetKwFunction{stripcoord}{get\_strip\_id}
	\SetKwFunction{newcoordinate}{new\_coordinate}
	\SetKwFunction{setassigndirection}{set\_assignment\_direction}
	\KwIn{Dimension sizes $\mathcal{D}$, set of relative vectors 
	$\mathcal{S}$ describing the~$k$-neighborhood,
	number of processes per compute node $n$
	and rank $r$ of calling process.}
	\KwResult{New coordinate $\vec{r}_\text{new}$ of calling rank.}
	\Indm\renewcommand{\nl}{\let\nl\oldnl}\stencilstrips{$\mathcal{D}, \mathcal{S}, n, r, \vec{r}_\text{new}$}\\
	\Indp
	$\alpha \gets \distortionfactors{$\mathcal{S}, |\mathcal{D}|, n$}$\\
	$s \gets \striplengths{$|\mathcal{D}|, n, \alpha$}$\\
	$v \gets \prod_{i \gets 0}^{|\mathcal{D}|-1} s[i]$\\ 
	\For{$i \gets 0$ \KwTo $|\mathcal{D}|-1$}{
		$v \gets \frac{v}{s[i]}$\\
		$s_\text{coord}[i] \gets \lfloor \frac{r}{v} \rfloor$\\
		$r \gets r - (s_\text{coord}[i]-1)v$\\
	}
	$\vec{r}_\text{new} \gets \newcoordinate{$\mathcal{D}, s, s_\text{coord}, r$}$\\
	\caption{\texttt{Stencil Strips}\xspace Algorithm}
	\label{stencilstripscoord}
\end{algorithm}

During early experimentation, we noticed that a consecutive assignment of processes to 
compute nodes for grids where the dimension sizes were close to the $d$th root
of $n$ (\textrm{i.e.}\xspace, optimal side length) and the nearest neighbor stencil
resulted in lower values of $J_\text{sum}$ and
$J_\text{max}$ in comparison to recursive bisection.
This observation was the
motivation for this algorithm. The main idea is to partition the grid into
strips, where strip lengths are chosen such that they are close to the scaled
length of an optimal bounding rectangle of $\mathcal{S}$ (\textrm{e.g.}\xspace, $\sqrt[\leftroot{-3}\uproot{3}d]{n}$ for
the nearest neighbor stencil). Processes in
coherent strips are~assigned~to~partitions/nodes.

For that purpose, let $R_{i}$ be the $i$th component of $\vec{R}$, then we define
$n_{i, \text{max}} := \max \{R_{i} \mid \forall \vec{R} \in \mathcal{S}\}$ and
$n_{i, \text{min}} := \min \{R_{i}\mid \forall \vec{R} \in \mathcal{S}\}$
to be the maximal and minimal value along dimension $i$ in the $k$-neighborhood
$\mathcal{S}$, respectively. The \emph{extension} $e_i$ of $\mathcal{S}$ along dimension $i$ is
$e_i := n_{i, \text{max}} - n_{i, \text{min}}$.
We use the extensions to find a bounding $n$-dimensional rectangle of
$\mathcal{S}$, with dimension sizes $\mathcal{E} = [e_0, \dots,
e_{d-1}]$.  We define the volume $V_b$ of the bounding $n$-dimensional
rectangle to be \begin{equation*} 
	V_b := \prod_{i=0}^{d-1} \epsilon_i, \quad \epsilon_i=
	\begin{cases}
		1 & \text{if } e_i = 0\\
		e_i & \text{else}.
	\end{cases}
\end{equation*}
With $d_b:= |\{e_i \mid \forall e_i \in \mathcal{E} \land e_i \neq 0\}|$ being
the number of non-zero expansions, we define the \emph{distortion factor} $\alpha_i$ to 
a~$d_b$-dimensional cube along dimension $i$ as
\begin{equation*} 
	\alpha_i := \frac{e_i}{\sqrt[\leftroot{-3}\uproot{3}d_b]{V_b}}. 
\end{equation*} 
As input for the algorithm serves the dimension sizes~$\mathcal{D}$ of the
grid, the~$k$-neighborhood $\mathcal{S}$, the number of processes per node~$n$
and the rank $r$ of the calling process. The pseudo-code is presented in
Algorithm~\ref{stencilstripscoord}. The algorithm outputs the new position
vector $\vec{r}_\text{new}$ of the calling rank. The algorithm itself works as
follows: we calculate a list of distortion factors $\alpha$ for every dimension
(Line~$1$). With the distortion factors, we can calculate a list $s$, containing
the optimal strip length~$s_i$ along dimension~$i$ (Line~$2$), with
consideration to already computed strip lengths, the distortion
factors~$\alpha_i$ and the node size~$n$:
\begin{equation*}
	s_i := \sqrt[\leftroot{-3}\uproot{3}d-i]{\frac{\alpha_i n}{\prod_{j=0}^{i-1} s_j}}.
\end{equation*}
This is done for every dimension except the largest one, for we iteratively position
the strips along the largest dimension. To be more precise, we assume that
along the largest dimension the strip length is one. In every other dimension
$d_i$, we fit~$\lfloor \frac{d_i}{s_i} \rfloor$ strips (the last strips is of
size~$s_i + \frac{d_i}{s_i}$). 
Every process can locally compute the number of small and large
strips along the dimensions, how many processes fit into each strip
and, with the calling process' rank~$r$ the coordinates $s_\text{coord}$ of
the strip in which it is contained (Line~$3$-$5$),
from which it can calculate its new position in the strip and thus, in the grid
$\vec{r}_\text{new}$~(Line~$9$). In order to obtain cohesive partitions, we flip the
strip assignment direction accordingly as depicted in
Figure~\ref{fig:stencilstrips}.

The strip widths, strip coordinates and the ranks position in the strip can
each be calculated in~$\mathcal{O}(d)$. The most expensive part of the 
algorithm is calculating the distortion factors, for which we need the maximal
and minimal expansion along the dimensions. This amounts in $\mathcal{O}(k d)$
steps. Thus, the run-time of the algorithm is bounded by $\mathcal{O}(k d)$.
\begin{figure}[t]
	\begin{subfigure}{0.2\textwidth}
		\center
	\includegraphics[trim = 0 0 0 0, clip, height=2cm]{./images/strip_assignment_correct} \\
	\caption{}
	\end{subfigure}    
	\hfill	
	\begin{subfigure}{0.2\textwidth}
		\center
	\includegraphics[trim = 0 0 0 0, clip, height=2cm]{./images/strip_assignment_incorrect} \\
	\caption{}
	\end{subfigure}
	\caption{Assignment direction of the \texttt{Stencil Strips}\xspace algorithm.
	(a): Alternating strips assignment direction accordingly guarantees 
	coherent partitions.
	(b): Imprudent assignment direction can lead to incoherent partitions.
	}
	\label{fig:stencilstrips}
\end{figure}

\section{Experimental evaluation}
\label{sec:experiments}
We have implemented the presented algorithms for Cartesian rank
reordering and the algorithm presented by Gropp~\cite{Gropp2019}, 
in accordance with the detailed pseudo-code of his paper.
We aim to show the advantage of approaches that do not rely on factorization of
the number of processes per node $n$. For that purpose, we compare the presented 
algorithms to a blocked assignment of ranks to nodes (henceforth referred as
to blocked), \texttt{Nodecart}\xspace, and \texttt{VieM}\xspace, which are described in
Section~\ref{sec:relatedwork}.  We start by describing the
systems used to run the experiments in Section~\ref{sec:systems}, implementation and
benchmarking details are found in Section~\ref{sec:methodology}. In
Section~\ref{sec:reduction}, we visualize the distribution of the reduction in
inter-node communication achieved by the different algorithms over the blocked
mapping for a wide set of instances. We proceed in
Section~\ref{sec:throughput}, where we present the throughput gained by our
algorithms for a neighbor all-to-all exchange and different message sizes on
different machines. Finally, we compare the algorithmic running time in
Section~\ref{sec:instantiation}.  
\begin{figure*}[t]
	\captionsetup[sub]{labelformat=empty, skip=0pt}
	\subcaptionbox{Nearest neighbor}[\textwidth]['u']{
		\begin{subfigure}{0.24\textwidth}
			\center
			\caption{}
			\includegraphics[trim=0 20 0 10, clip, height=3.3cm]{./images/five_point_50_48}
		\end{subfigure}
		\hfill
		\begin{subfigure}{0.24\textwidth}
			\centering
			\caption{\emph{VSC}$4$\xspace}
			\includegraphics[trim = 0 20 0 30, clip, height=3.1cm]{./images/bandwidth_VSC4/barplot_five_point_50_48}
			\end{subfigure}
		\begin{subfigure}{0.24\textwidth}
			\centering
			\caption{\emph{SuperMUC-NG}\xspace}
			\includegraphics[trim = 0 20 0 30, clip, height=3.1cm]{./images/bandwidth_superMUC/barplot_five_point_50_48}
		\end{subfigure}
		\begin{subfigure}{0.24\textwidth}
			\centering
			\caption{\emph{JUWELS}\xspace}
                         \includegraphics[trim = 0 20 0 30, clip, height=3.1cm]{./images/bandwidth_JUWELS/barplot_five_point_50_48}
		\end{subfigure}
	}
	\captionsetup[sub]{labelformat=empty, skip=0pt}
	\subcaptionbox{Nearest neighbor with hops}[\textwidth]['u']{
		\begin{subfigure}{0.24\textwidth}
			\center
			\caption{}
			\includegraphics[trim=0 20 0 10, clip, height=3.3cm]{./images/five_point_hops_first_50_48}
		\end{subfigure}
		\hfill
		\begin{subfigure}{0.24\textwidth}
			\centering
			\includegraphics[trim = 0 20 0 30, clip, height=3.1cm]{./images/bandwidth_VSC4/barplot_five_point_hops_first_50_48}
			\end{subfigure}
		\begin{subfigure}{0.24\textwidth}
			\centering
			\includegraphics[trim = 0 20 0 30, clip, height=3.1cm]{./images/bandwidth_superMUC/barplot_five_point_hops_first_50_48}
		\end{subfigure}
		\begin{subfigure}{0.24\textwidth}
			\centering
			\includegraphics[trim = 0 20 0 30, clip, height=3.1cm]{./images/bandwidth_JUWELS/barplot_five_point_hops_first_50_48}
		\end{subfigure}
	}
	\captionsetup[sub]{labelformat=empty, skip=0pt}
	\subcaptionbox{Component stencil}[\textwidth]['u']{
		\begin{subfigure}{0.24\textwidth}
			\center
			\caption{}
			\includegraphics[trim=0 20 0 10, clip, height=3.3cm]{./images/component_50_48}
			\label{fig:relnn50}
		\end{subfigure}
		\hfill
		\begin{subfigure}{0.24\textwidth}
			\centering
			\includegraphics[trim = 0 0 0 30, clip, height=3.2cm]{./images/bandwidth_VSC4/barplot_component_50_48}
			\end{subfigure}
		\begin{subfigure}{0.24\textwidth}
			\centering
			\includegraphics[trim = 0 0 0 30, clip, height=3.2cm]{./images/bandwidth_superMUC/barplot_component_50_48}
		\end{subfigure}
		\begin{subfigure}{0.24\textwidth}
			\centering
			\includegraphics[trim = 0 0 0 30, clip, height=3.2cm]{./images/bandwidth_JUWELS/barplot_component_50_48}
		\end{subfigure}
	}
	\caption{Left column: scores of the algorithms (smaller is better).
	Right three columns: Speedup over blocked mapping (higher is better) with
	$N=50$ number of nodes and $p=48$ processes per node and grid sizes
	$50\times 48$ for the nearest neighbor stencil,
	the nearest neighbor stencil with hops and the
	component stencil on the \emph{VSC}$4$\xspace, \emph{SuperMUC-NG}\xspace, and \emph{JUWELS}\xspace. Absolute times 
	are written in the corresponding bar.
	}\label{fig:bandwidth50}
\end{figure*}
\begin{figure*}[ht]
	\captionsetup[sub]{labelformat=empty, skip=0pt}
	\subcaptionbox{Nearest neighbor}[\textwidth]['l']{
		\begin{subfigure}{0.24\textwidth}
			\center
			\caption{}
			\includegraphics[trim=0 20 0 10, clip, height=3.3cm]{./images/five_point_100_48}
		\end{subfigure}
		\hfill
		\begin{subfigure}{0.24\textwidth}
			\centering
			\caption{\emph{VSC}$4$\xspace}
			\includegraphics[trim = 0 20 0 30, clip, height=3.1cm]{./images/bandwidth_VSC4/barplot_five_point_100_48}
			\end{subfigure}
		\begin{subfigure}{0.24\textwidth}
			\centering
			\caption{\emph{SuperMUC-NG}\xspace}
			\includegraphics[trim = 0 20 0 30, clip, height=3.1cm]{./images/bandwidth_superMUC/barplot_five_point_100_48}
		\end{subfigure}
		\begin{subfigure}{0.24\textwidth}
			\centering
			\caption{\emph{JUWELS}\xspace}
			\includegraphics[trim = 0 20 0 30, clip, height=3.1cm]{./images/bandwidth_JUWELS/barplot_five_point_100_48}
		\end{subfigure}
	}
	\captionsetup[sub]{labelformat=empty, skip=0pt}
	\subcaptionbox{Nearest neighbor with hops}[\textwidth]['l']{
		\begin{subfigure}{0.24\textwidth}
			\center
			\caption{}
			\includegraphics[trim=0 20 0 10, clip, height=3.3cm]{./images/five_point_hops_first_100_48}
		\end{subfigure}
		\hfill
		\begin{subfigure}{0.24\textwidth}
			\centering
			\includegraphics[trim = 0 20 0 30, clip, height=3.1cm]{./images/bandwidth_VSC4/barplot_five_point_hops_first_100_48}
			\end{subfigure}
		\begin{subfigure}{0.24\textwidth}
			\centering
			\includegraphics[trim = 0 20 0 30, clip, height=3.1cm]{./images/bandwidth_superMUC/barplot_five_point_hops_first_100_48}
		\end{subfigure}
		\begin{subfigure}{0.24\textwidth}
			\centering
			\includegraphics[trim = 0 20 0 30, clip, height=3.1cm]{./images/bandwidth_JUWELS/barplot_five_point_hops_first_100_48}
		\end{subfigure}
	}
	\captionsetup[sub]{labelformat=empty, skip=0pt}
	\subcaptionbox{Component stencil}[\textwidth]['l']{
		\begin{subfigure}{0.24\textwidth}
			\center
			\caption{}
			\includegraphics[trim=0 20 0 10, clip, height=3.3cm]{./images/component_100_48}
		\end{subfigure}
		\hfill
		\begin{subfigure}{0.24\textwidth}
			\centering
			\includegraphics[trim = 0 0 0 30, clip, height=3.2cm]{./images/bandwidth_VSC4/barplot_component_100_48}
			\end{subfigure}
		\begin{subfigure}{0.24\textwidth}
			\centering
			\includegraphics[trim = 0 0 0 30, clip, height=3.2cm]{./images/bandwidth_superMUC/barplot_component_100_48}
		\end{subfigure}
		\begin{subfigure}{0.24\textwidth}
			\centering
			\includegraphics[trim = 0 0 0 30, clip, height=3.2cm]{./images/bandwidth_JUWELS/barplot_component_100_48}
		\end{subfigure}
	}
	\caption{Left column: scores of the algorithms (smaller is better). 
	Right three columns: Speedup over blocked mapping (higher is better) with
	$N=100$ number of nodes and $p=48$ processes per node and grid sizes
	$75\times 64$ for the nearest neighbor stencil,
	the nearest neighbor stencil with hops and the
	component stencil on the \emph{VSC}$4$\xspace, \emph{SuperMUC-NG}\xspace, and \emph{JUWELS}\xspace. Absolute times 
	are written in the corresponding bar.
	}\label{fig:bandwidth100}
\end{figure*}

\subsection{Machine Description}
\label{sec:systems}
\begin{table}
	\caption{Machines used in the Experiments.}
	\label{table:machines}
	\begin{scriptsize}
		\begin{tabularx}{\the\columnwidth}{l l c r}
		\hline
		Name & Processor & \texttt{MPI}\xspace libraries & Compiler \\
		\hline
		\emph{VSC}$4$\xspace & Intel Skylake Platinum 8174 & Intel\,MPI & icc 19.0.5\\
		\emph{SuperMUC-NG}\xspace & Intel Skylake Platinum 8174 & Intel\,MPI & icc 19.0.5\\
			\emph{JUWELS}\xspace & Intel Xeon Platinum 8168 & Intel\,MPI & icc 19.0.3 \\
			\hline
	\end{tabularx}
	\end{scriptsize}
\end{table}
We perform the experiments in Section~\ref{sec:throughput} 
on the Vienna Scientific Cluster $4$ (\emph{VSC}$4$\xspace), \emph{SuperMUC-NG}\xspace and
\emph{JUWELS}\xspace. \emph{VSC}$4$\xspace is composed of $790$ nodes, where each node is equipped with two Intel Skylake Platinum $8174$ processors running at
\SI{3.1}{\GHz}, \textrm{i.e.}\xspace, each node has $48$ cores. 
The nodes are connected with a two-level fat-tree (blocking factor 2:1) via
OmniPath with a capacity of \SI[per-mode=symbol]{100}{\giga\bit\per\second}.
\emph{SuperMUC-NG}\xspace consists of $\num{6336}$ compute nodes, also equipped with two Intel Skylake Platinum $8174$ processors running at
\SI{3.1}{\GHz}.
The nodes are bundled into islands, with a pruned OmniPath connection
between the islands (pruning factor $1$:$4$) and nodes within an island are connected
via a OmniPath fat-tree.
 \emph{JUWELS}\xspace is composed of $\num{2271}$ compute
nodes, with two Intel Xeon Platinum $8168$, running at \SI{2.7}{\GHz}. The
nodes are connected with a two-level fat-tree InfiniBand network (pruning factor 2:1)
with a capacity of \SI[per-mode=symbol]{100}{\giga\bit\per\second}. A brief summary of the
systems and libraries is given in
Table~\ref{table:machines}.

\subsection{Experimental Setup}
\label{sec:methodology}
\begin{lstlisting}[frame=lines, caption={Interface used for a $k$-neighborhood aware \texttt{MPI}\xspace
Cartesian communicator}, label={stencil_cart}]
int MPIX_Cart_stencil_comm(MPI_Comm oldcomm,
    const int ndims, const int dims[],
    const int periods[], const int reorder,
    const int stencil[], const int k,
    MPI_Comm *cartcomm);
\end{lstlisting}

All algorithms were implemented in \mbox{C\texttt{++}}\xspace$11$.  The code was compiled
with the Intel compiler and full optimization flags (-O$3$) and Intel
MPI.


\begin{figure*}[ht]
	\begin{subfigure}{0.31\textwidth}
		\centering
		\includegraphics[trim = 0 15 0 25, clip, scale=0.35]{./images/NearestNeighbor}
	\caption{}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.31\textwidth}
		\includegraphics[trim = 0 15 0 25, clip, scale=0.35]{./images/NearestNeighborHops}
	\caption{}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.31\textwidth}
		\includegraphics[trim = 0 15 0 25, clip, scale=0.35]{./images/component_}
	\caption{}
	\end{subfigure}
	\caption{Reduction for
	nearest neighbor stencil (a), for the nearest neighbor stencil with
	additional two hops along the first dimension (b) and for the
	artificial component stencil (c). Lower is better. For each algorithm \texttt{sum}
	defines the total reduction, whereas \texttt{max} defines the reduction
	of the bottleneck node. \texttt{Hyperplane}\xspace depicted in blue, $k$-$d$~\texttt{tree}\xspace in
	yellow, \texttt{Stencil Strips}\xspace in green, \texttt{Nodecart}\xspace in red, and \texttt{VieM}\xspace as
	comparison in gray.} 
	\label{fig:distribution}
\end{figure*}
To remap for arbitrary stencils that cannot be expressed with the
\texttt{MPI}\xspace Cartesian interfaces, we use an interface similar to
\texttt{MPI\_\-Cart\_\-create}\xspace as shown in Listing~\ref{stencil_cart}. The array
\texttt{stencil[]} consists of a flattened list of
relative offsets along each dimension per neighbor. The number of
neighbors is~\texttt{k}, thus, the array
\texttt{stencil[]} is of length~$\texttt{k}\cdot \texttt{ndims}$.

As for the $k$-neighborhoods, we choose to test with the presented 
$2$-dimensional stencils from Section~\ref{sec:preliminaries}, \textrm{i.e.}\xspace, 
nearest neighbor stencil, the nearest neighbor stencil with hops and
the component stencil, as depicted in Figure~\ref{fig:stencils}. 

All grids were created according to the \texttt{MPI\_\-Dims\_\-create}\xspace specifications,
that is with the sizes of the dimensions being as close as possible to
each other~\cite{MPI-3.1,Traff15:dimscreate}.

For the stencil exchange, we instantiated a distributed graph communicator from
the Cartesian communicator and the $k$-neighborhood in order to call the
\texttt{MPI\_\-Neighbor\_\-alltoall}\xspace routine. Synchronization between the processes before each
collective message exchange was done with an \texttt{MPI\_\-Barrier}\xspace and we define the time
needed for the operation, as the maximal time any process spent in the
\texttt{MPI\_\-Neighbor\_\-alltoall}\xspace routine.  As we assumed unit-weighted communication edges,
every process sends and receives the same amount of data to its communication neighbors.

\subsection{Inter-Node Communication Analysis}
\label{sec:reduction}
\begin{figure}[b]
	\center
	\includegraphics[trim= 0 70 0 30, clip, scale=0.35]{./images/Instantiation_time_100_48_five_point}
	\caption{Instantiation time in \si{\milli\second} on a $100\times 48$
	nearest neighbor stencil instance, run on \emph{VSC}$4$\xspace (lower is better). Each algorithm was
	instantiated $200$ times and after outlier removal, we plot the
	mean with a $95\%$ confidence interval. \texttt{VieM}\xspace was omitted since it 
	took two orders of magnitude longer and would distort~the~scale.}
	\label{fig:instantiation} 
\end{figure}

In this section, we investigate the
reduction of~$J_\text{sum}$ and~$J_\text{max}$ of different algorithms over a blocked mapping.
Note that the evaluation of the objective function is machine independent.
We evaluate the performance of the algorithms on a varying input
number of compute nodes, processes per compute node, and number of dimensions.  To
be more precise, the number of nodes is given by the set~$\mathcal{N}=\{10, 13,
16, \dots, 33\}$ and the number of processes per node by the
set~$\mathcal{P}=\{10, 13, 16, \dots, 31\} \cup \{32\}$, while we restrict the set
of dimensions to be $\mathcal{D}=\{2, 3\}$.  The set of instances is therefore
the Cartesian product~$\mathcal{I}=\mathcal{N} \times \mathcal{P} \times
\mathcal{D}$, with $|\mathcal{I}|=144$. We define the reduction to
be~$\frac{C_X}{C_b}$, where~$C_X$ is the amount of inter-node communication
induced by algorithm~$X$ and~$C_b$ the amount of inter-node communication
induced by a blocked mapping. We compare both the
total reduction of inter-node communication~$J_\text{sum}$, as well as the
reduction over the bottleneck node~$J_\text{max}$, as defined in
Section~\ref{sec:preliminaries}.  We evaluate the presented algorithms and
compare ourselves to Gropp's algorithm (\texttt{Nodecart}\xspace) and \texttt{VieM}\xspace, both described
in Section~\ref{sec:relatedwork}.  The configuration for \texttt{VieM}\xspace was set to the
strongest setting, prioritizing quality in terms of reduction of inter-node
communication instead of speed. As for the local-search neighborhood, we
allowed swaps between any connected pair of vertices, \textrm{i.e.}\xspace, we considered the
largest search space. We set the \texttt{hierarchy\_parameter\_string}
and the \texttt{distance\_parameter\_string} to the $n$:$N$ and $0$:$1$, respectively
so that \texttt{VieM}\xspace optimizes our objective function. 

The results for the three different $k$-neighborhood communication patterns can
be seen in Figure~\ref{fig:distribution}. We plot the median improvement (orange)
with a~$95\%$ confidence interval (indicated by the notches), calculated with a
Gaussian-based asymptotic approximation. Since the confidence interval of
the medians do not overlap, we can say with statistical evidence~\cite{murphy2019probability} that the median
reduction improvement of \texttt{Hyperplane}\xspace and \texttt{Stencil Strips}\xspace is better than
\texttt{Nodecart}\xspace in all communication patterns.
Only for the nearest neighbor with hops pattern, there
seems to be no statistical reduction
between the $k$-$d$~\texttt{tree}\xspace and \texttt{Nodecart}\xspace. For the
nearest neighbor and component stencil, the reduction distribution between the
\texttt{Stencil Strips}\xspace and \texttt{VieM}\xspace looks very similar and there is no statistical
difference in the median reduction, indicating that the \texttt{Stencil Strips}\xspace
algorithm, which is not bound to factorization in any way, can greatly exploit
the structure of the grid partitioning problem.

\subsection{Throughput Analysis}
\label{sec:throughput}
We continue by examining the influence of rank reordering on the time
needed for an \texttt{MPI\_\-Neighbor\_\-alltoall}\xspace exchange and different message sizes to be
sent to each communication partner. The experiment was conducted on
\emph{VSC}$4$\xspace, \emph{SuperMUC-NG}\xspace, and \emph{JUWELS}\xspace described in Section~\ref{sec:systems}. In order to 
see how the reordering affects the communication performance of larger instances,
we perform the experiments with $50$ and $100$ nodes and
for each used the full number of processes per node $48$. Each message was sent
$200$ times to obtain large samples.
From each run, we capture the maximum time needed over all processes. Before
each exchange, we synchronize the processes with \texttt{MPI\_\-Barrier}\xspace. In
Figure~\ref{fig:bandwidth50} and Figure~\ref{fig:bandwidth100}, we plot the
mean speedup over the blocked assignments, after removing outliers beyond $1.5$
inter-quartile range from the third and first quartile, respectively.
Along with the speedup, we plot the scores of the
algorithms sorted in increasing order of $J_\text{sum}$ and $J_\text{max}$.
The message size is the number of bytes send per neighbor.
In the appendix we include tables with the absolute
running times and the communication performance of a random assignment of
processes to compute nodes (\texttt{Random}\xspace). Due to constraint space, we omitted
\texttt{Random}\xspace in the speedup plots.

For the nearest neighbor stencil and $N=50$, Figure~\ref{fig:bandwidth50},
\texttt{Hyperplane}\xspace and \texttt{Stencil Strips}\xspace are up to three and four times faster than the
blocked assignment on the \emph{VSC}$4$\xspace, up to two and almost
up to three times faster on \emph{SuperMUC-NG}\xspace and more than five times faster on
\emph{JUWELS}\xspace. All of the presented algorithms are able to outperform \texttt{Nodecart}\xspace in
terms of the mapping metric $J_\text{sum}$, $J_\text{max}$ and communication
performance on \emph{VSC}$4$\xspace and mostly on \emph{SuperMUC-NG}\xspace and \emph{JUWELS}\xspace.  As for \texttt{VieM}\xspace,
\texttt{Stencil Strips}\xspace outperforms it consistently on all machines, whereas
\texttt{Hyperplane}\xspace and $k$-$d$~\texttt{tree}\xspace attain similar or better performance.
Surprisingly, the communication performance of \texttt{VieM}\xspace on \emph{JUWELS}\xspace is worse than
predicted by $J_\text{sum}$, $J_\text{max}$.  In the case of $N=100$,
Figure~\ref{fig:bandwidth100}, \texttt{Hyperplane}\xspace and \texttt{Stencil Strips}\xspace obtain speedups
of up to three over the blocked assignment on \emph{VSC}$4$\xspace, up to
$2.5$ on \emph{SuperMUC-NG}\xspace and up to $2$ on \emph{JUWELS}\xspace. For message sizes
$\num{65536}$, $\num{262144}$ the performance of the blocked assignment,
\texttt{Hyperplane}\xspace, \texttt{Stencil Strips}\xspace, and \texttt{Nodecart}\xspace performed very badly,
resulting in speedups of more than forty, which is why we cut the axis at $4$. 

In the case of the nearest neighbor with hops stencil for $N=50$, Figure~\ref{fig:bandwidth50},
\texttt{Hyperplane}\xspace, \texttt{Stencil Strips}\xspace obtained speedups of up to four on the \emph{VSC}$4$\xspace,
up to $3.5$ on \emph{SuperMUC-NG}\xspace and up to three on \emph{JUWELS}\xspace. \texttt{VieM}\xspace outperforms the presented 
algorithms and both \emph{VSC}$4$\xspace and \emph{SuperMUC-NG}\xspace, but \texttt{Stencil Strips}\xspace obtained similar
to better performance for large message sizes on \emph{JUWELS}\xspace. In comparison to \texttt{Nodecart}\xspace,
the presented algorithms seemed to be up to two and three times faster. 
For $N=100$, Figure~\ref{fig:bandwidth100}, the presented algorithms are up to a factor
$4.5$ faster than the blocked assignment on \emph{VSC}$4$\xspace, up to $3.5$ faster 
on \emph{SuperMUC-NG}\xspace and more than a factor $2$ faster on \emph{JUWELS}\xspace. It seems that the 
performance of the blocked assignment for message size
$\num{65536}$ performed relatively bad in comparison to the other message
sizes. On the \emph{VSC}$4$\xspace, our algorithms are more than a factor two faster than
\texttt{Nodecart}\xspace, whereas on \emph{SuperMUC-NG}\xspace and \emph{JUWELS}\xspace we outperform \texttt{Nodecart}\xspace for
most message sizes by a factor of more than $1.3$.

As for the synthetic component stencil and $N=50$, only $k$-$d$~\texttt{tree}\xspace and \texttt{Stencil Strips}\xspace 
managed to find an optimal mapping, where each compute node has two outgoing communication
edges, as described in Section~\ref{sec:hardness}. This reordering results in speedups
of up to ten over the blocked assignment on \emph{VSC}$4$\xspace, up to five on \emph{SuperMUC-NG}\xspace
and up twelve on \emph{JUWELS}\xspace. Interestingly, \texttt{VieM}\xspace seems to have similar performance on 
\emph{VSC}$4$\xspace and \emph{SuperMUC-NG}\xspace, and outperforms our algorithms for small message sizes on 
\emph{JUWELS}\xspace. 
In the case of $N=100$, again only $k$-$d$~\texttt{tree}\xspace and \texttt{Stencil Strips}\xspace found optimal mappings,
leading to speedups of up to fourteen on \emph{VSC}$4$\xspace, five on \emph{SuperMUC-NG}\xspace, and sixteen on
\emph{JUWELS}\xspace. Surprisingly, \texttt{VieM}\xspace seems to match and outperform the performance
of our algorithms, even though it has a factor $1.6$ more $J_\text{sum}$
and a factor $8.5$ greater $J_\text{max}$.

\subsection{Instantiation Time}
\label{sec:instantiation}

Since the theoretical complexity of the three presented algorithms, \texttt{Nodecart}\xspace and \texttt{VieM}\xspace
all dependent on different parameters, we benchmark the algorithmic runtime
needed to calculate the new ranks only on the largest nearest neighbor stencil
instance, described in Section~\ref{sec:throughput} ($N=100$) on \emph{VSC}$4$\xspace. The algorithmic runtime
does not include the instantiation time of the reordered Cartesian
communicator, rather, it only includes the necessary operations to calculate
the new ranks. For the presented algorithms, this includes creating a sorted
communicator, with precise knowledge over which ranks are assigned to which
nodes. This custom communicator consists of a \emph{node communicator} (grouped
by processes having shared memory) and a \emph{leader communicator}, which
consists of one process per compute node. In this benchmark, the $k$-$d$~\texttt{tree}\xspace
algorithm was implemented with a linear search, for finding the dimension along
which to split, resulting in a theoretical run-time of $\mathcal{O}(d\log p)$.
Each algorithm was instantiated~$200$ times for large sample sizes, and we measured the longest time
needed over all processes. Before each instantiation, the processes were
synchronized using \texttt{MPI\_\-Barrier}\xspace. In Figure~\ref{fig:instantiation}, we plot the
mean instantiation time, after outlier removal (beyond $1.5$ inter-quartile
range from the first and third quartile) with a $95\%$ confidence interval. We
did not include \texttt{VieM}\xspace, since it took on average over \SI{7.95}{\second} 
(a factor of more than $400$) and it would distort the scale of the plot.

We can see, that for this instance the \texttt{Hyperplane}\xspace and $k$-$d$~\texttt{tree}\xspace algorithm 
are the two fastest, and statistically equivalent. \texttt{Nodecart}\xspace seems to 
need $28\%$ longer, whereas the convoluted process of calculating the
strips and strip positions results in \texttt{Stencil Strips}\xspace being the slowest 
algorithm, more than a factor $2$ slower than \texttt{Hyperplane}\xspace and $k$-$d$~\texttt{tree}\xspace.

\section{Conclusion}
\label{sec:conclusion}
We introduced three new efficient algorithms for process to compute node assignment on Cartesian grids and stencils communication patterns. 
By thoroughly exploiting the inherently present structure of the problem, we arrive at algorithms that outperform the state-of-the-art in terms of running time and communication performance.
We implemented the algorithms for \texttt{MPI\_\-Cart\_\-create}\xspace as reordering functions and performed
extensive benchmarks. 
 An intensive experimental evaluation shows that our algorithms are up to two
 orders of magnitude faster in running time than a (sequential) high-quality
 general graph mapping tool \texttt{VieM}\xspace, while obtaining similar results in
 communication performance. Furthermore, our algorithms are three times
 faster in an \texttt{MPI\_\-Neighbor\_\-alltoall}\xspace exchange than a state-of-the-art Cartesian grid mapping
 algorithm (\texttt{Nodecart}\xspace) by achieving a significantly better mapping quality.
 Considering the good results, we plan to release the algorithms and integrate
 them into publicly available \texttt{MPI}\xspace implementations.

\subsection*{Acknowledgments}
This work was partially supported by the Austrian Science Fund (FWF): project P
31763-N31.  The computational results presented have been achieved in part
using the Vienna Scientific Cluster (VSC).  We acknowledge PRACE for awarding
us access to JUWELS at GCS@FZJ, Germany and SuperMUC-NG at GCS@LRZ, Germany.

\FloatBarrier
\bibliographystyle{IEEEtran}
