\label{sec:hpl}

For this work, we use the freely-available reference-implementation of the High-Performance Linpack benchmark\cite{HPL}, HPL, which is used to benchmark systems for the TOP500\cite{top500} list. HPL requires MPI to be available and implements a LU decomposition, \ie a factorization of a square matrix \(A\) as the product of a lower triangular matrix \(L\) and an upper triangular matrix \(U\). HPL checks the correctness of this factorization by solving a linear system \(A\cdot{}x=b\), but only the factorization step is benchmarked. The factorization is based on a right-looking variant of the LU factorization with row partial pivoting and allows multiple look-ahead depths. The working principle of the factorization is depicted in Figure\ref{fig:hpl_overview} and consists of a series of panel factorizations followed by an update of the trailing sub-matrix. HPL uses a two-dimensional block-cyclic data distribution of \(A\) and implements several custom collective communication algorithms to efficiently overlap communication with computation. The main parameters of HPL are listed subsequently:

• \(N\) is the order of the square matrix \(A\).
• NB is the ``blocking factor'', \ie the granularity at which HPL operates when panels are distributed or worked on.
• \(P\) and \(Q\) denote the number of process rows and the number of process columns, respectively.
• RFACT determines the panel factorization algorithm. Possible values are Crout, left- or right-looking.
• SWAP specifies the swapping algorithm used while pivoting. Two algorithms are available: one based on binary exchange (along a virtual tree topology) and the other one based on a spread-and-roll (with a higher number of parallel communications). HPL also provides a panel-size threshold triggering a switch from one variant to the other.
• BCAST sets the algorithm used to broadcast the panel of columns to the other process columns. Legacy versions of the MPI standard only supported non-blocking point-to-point communications but did not support non-blocking collective communications, which is why HPL ships with in total 6 self-implemented variants to efficiently overlap the time spent waiting for an incoming panel with updates to the trailing matrix: ring, ring-modified, 2-ring, 2-ring-modified, long, and long-modified. The modified versions guarantee that the process right after the root (\ie the process that will become the root in the next iteration) receives data first and does not participate further in the broadcast. This process can thereby start working on the panel as soon as possible. The ring and 2-ring versions correspond to the name-giving two virtual topologies while the long version is a spread and roll algorithm where messages are chopped into \(Q\) pieces. This generally leads to better bandwidth exploitation. The ring and 2-ring variants rely on MPI_Iprobe, meaning they return control if no message has been fully received yet and hence facilitate partial overlapping of communication with computations. In HPL 2.2 and 2.1, this capability has been deactivated for the long and long-modified algorithms. A comment in the source code states that some machines apparently get stuck when there are too many ongoing messages.
• DEPTH controls how many iterations of the outer loop can overlap with each other.

The sequential complexity of this factorization is \(\mathrm{flop}(N) = \frac{2}{3}N^3 + 2N^2 + \O(N)\) where \(N\) is the order of the matrix to factorize. The time complexity can be approximated by \[T(N) \approx \frac{\left(\frac{2}{3}N^3 + 2N^2\right)}{P\cdot{}Q\cdot{}w} + \Theta((P+Q)\cdot{}N^2),\] where \(w\) is the flop rate of a single node and the second term corresponds to the communication overhead which is influenced by the network capacity and by the previously listed parameters (RFACT, SWAP, BCAST, DEPTH, \ldots). After each run, HPL reports the overall flop rate \(\mathrm{flop}(N)/T(N)\) (expressed in \si{\giga\flops}) for the given configuration. See Figure\ref{fig:hpl_output} for a (shortened) example output.

A large-scale execution of HPL on a real machine in order to submit to the TOP500 can therefore be quite time consuming as all the BLAS kernels, the MPI runtime, and HPL's numerous parameters need to be tuned carefully in order to reach optimal performance.

\label{sec:stampede} In June 2013, the Stampede supercomputer at TACC was ranked 6th in the TOP500 by achieving \SI{5168.1}{\tera\flops} and was still ranked 20th in June 2017. In 2017, this machine got upgraded and renamed Stampede2. The Stampede platform consisted of 6400 Sandy Bridge nodes, each with two 8-core Xeon E5-2680 and one Intel Xeon Phi KNC MIC coprocessor. The nodes were interconnected through a \SI{56}{\giga\bit\per\second} FDR InfiniBand 2-level Clos fat-tree topology built on Mellanox switches. As can be seen in Figure\ref{fig:fat_tree_topology}, the 6400 nodes are divided into groups of 20, with each group being connected to one of the 320 36-port switches (\SI{4}{\tera\bit\per\second} capacity), which are themselves connected to 8 648-port ``core switches'' (each with a capacity of \SI{73}{\tera\bit\per\second}). The peak performance of the 2 Xeon CPUs per node was approximately \SI{346}{\giga\flops}, while the peak performance of the KNC co-processor was about \SI{1}{\tera\flops}. The theoretical peak performance of the platform was therefore \SI{8614}{\tera\flops}. However, in the TOP500, Stampede was ranked with \SI{5168}{\tera\flops}. According to the log submitted to the TOP500 (see Figure\ref{fig:hpl_output}) that was provided to us, this execution took roughly two hours and used \(77\times78 = 6,006\) processes. The matrix of order \(N = 3,875,000\) occupied approximately \SI{120}{\tera\byte} of memory, \ie \SI{20}{\giga\byte} per node. One MPI process per node was used and each node's computational resources (the 16 CPU-cores and the Xeon Phi) must have been controlled by OpenMP and/or Intel's MKL.

The performance achieved by Stampede, \SI{5168}{\tera\flops}, needs to be compared to the peak performance of the 6,006 nodes, \ie \SI{8084}{\tera\flops}. This difference may be attributed to the node usage (\eg the MKL), to the MPI library, to the network topology that may be unable to deal with the very intensive communication workload, to load imbalance among nodes because some node happens to be slower for some reason (defect, system noise, \ldots), to the algorithmic structure of HPL, etc. All these factors make it difficult to know precisely what performance to expect without running the application at scale.

It is clear that due to the level of complexity of both HPL and the underlying hardware, simple performance models (analytic expressions based on \(N, P, Q\) and estimations of platform characteristics as presented in Section\ref{sec:hpl}) may be able to provide trends but can by no means predict the performance for each configuration (\ie consider the exact effect of HPL's 6 different broadcast algorithms on network contention). Additionally, these expressions do not allow engineers to improve the performance through actively identifying performance bottlenecks. For complex optimizations such as partially non-blocking collective communication algorithms intertwined with computations, very faithful modeling of both the application and the platform is required. Given the scale of this scenario (3,785 steps on 6,006 nodes in two hours), detailed simulations quickly become intractable without significant effort.

The complex network optimizations done in real MPI implementations need to be considered when predicting performance of MPI applications. For instance, message size not only influences the network's latency and bandwidth factors but also the protocol used, such as ``eager'' or ``rendez-vous'', as they are selected based on the message size, with each protocol having its own synchronization semantics. To deal with this, SMPI relies on a generalization of the LogGPS model\cite{smpi} and supports specifying synchronization and performance modes. This model needs to be instantiated once per platform through a carefully controlled series of messages (MPI_Send and MPI_Recv) between two nodes and through a set of piece-wise linear regressions.

Modeling network topologies and contention is also difficult. SMPI relies on SimGrid's communication models where each ongoing communication is represented as a whole (as opposed to single packets) by a flow. Assuming steady-state, contention between active communications can be modeled as a bandwidth sharing problem that accounts for non-trivial phenomena (\eg RTT-unfairness of TCP, cross-traffic interference or network heterogeneity\cite{Velho_TOMACS13}). Communications that start or end trigger re-computation of the bandwidth sharing if needed. In this model, the time to simulate a message passing through the network is independent of its size, which is advantageous for large-scale applications frequently sending large messages. SimGrid does not model transient phenomena incurred by the network protocol but accounts for network topology and heterogeneity.

Finally, collective operations are also challenging, particularly since these operations often play a key factor to an application's performance. Consequently, performance optimization of these operations has been studied intensively. As a result, MPI implementations now commonly have several alternatives for each collective operation and select one at runtime, depending on message size and communicator geometry. SMPI implements collective communication algorithms and the selection logic from several MPI implementations (\eg Open MPI, MPICH), which helps to ensure that simulations are as close as possible to real executions. Although SMPI supports these facilities, they are not required in the case of HPL as it ships with its own implementation of collective operations.

In Section\ref{sec:relwork} we explained that SMPI relies on the online simulation approach. Since SimGrid is a sequential simulator, SMPI maps every MPI process of the application onto a lightweight simulation thread. These threads are then run one at a time, \ie in mutual exclusion. Every time a thread enters an MPI call, SMPI takes control and the time that was spent computing (isolated from the other threads) since the previous MPI call can be injected into the simulator as a virtual delay.

Mapping MPI processes to threads of a single process effectively folds them into the same address space. Consequently, global variables in the MPI application are shared between threads unless these variables are privatized and the simulated MPI ranks thus isolated from each other. Several technical solutions are possible to handle this issue\cite{smpi}. The default strategy in SMPI consists of making a copy of the data segment (containing all global variables) per MPI rank at startup and, when context switching to another rank, to remap the data segment via mmap to the private copy of that rank. SMPI also implements another mechanism relying on the dlopen function that saves calls to mmap when context switching.

This causes online simulation to be expensive in terms of both simulation time and memory since the whole parallel application is executed on a single node. To deal with this, SMPI provides two simple annotation mechanisms:

• Kernel sampling: Control flow is in many cases independent of the computation results. This allows computation-intensive kernels (\eg BLAS kernels for HPL) to be skipped during the simulation. For this purpose, SMPI supports annotation of regular kernels through several macros such as SMPI_SAMPLE_LOCAL and SMPI_SAMPLE_GLOBAL. The regularity allows SMPI to execute these kernels a few times, estimate their cost and skip the kernel in the future by deriving its cost from these samples, hence cutting simulation time significantly. Skipping kernels renders the content of some variables invalid but in simulation, only the behavior of the application and not the correctness of computation results are of concern.
• Memory folding: SMPI provides the SMPI_SHARED_MALLOC (SMPI_SHARED_FREE) macro to replace calls to malloc (free). They indicate that some data structures can safely be shared between processes and that the data they contain is not critical for the execution (\eg an input matrix) and that it may even be overwritten. SMPI_SHARED_MALLOC works as follows (see Figure\ref{fig:global_shared_malloc}) : a single block of physical memory (of default size \SI{1}{\mega\byte}) for the whole execution is allocated and shared by all MPI processes. A range of virtual addresses corresponding to a specified size is reserved and cyclically mapped onto the previously obtained physical address. This mechanism allows applications to obtain a nearly constant memory footprint, regardless of the size of the actual allocations.

As explained in Section\ref{sec:con:diff}, faithful prediction of HPL necessitates emulation, \ie to execute the code. HPL relies heavily on BLAS kernels such as dgemm (for matrix-matrix multiplication) or dtrsm (for solving an equation of the form \(Ax=b\)). An analysis of an HPL simulation with \(64\) processes and a very small matrix of order \(30,000\) showed that roughly \SI{96}{\percent} of the time is spent in these two very regular kernels. For larger matrices, these kernels will consume an even bigger percentage of the computation time. Since these kernels do not influence the control flow, simulation time can be reduced by substituting dgemm and dtrsm function calls with a performance model for the respective kernel. Figure\ref{fig:macro_simple} shows an example of this macro-based mechanism that allows us to keep HPL code modifications to an absolute minimum. The (1.029e-11) value represents the inverse of the flop rate for this computation kernel and was obtained through calibration. The estimated time for the real kernel is calculated based on the parameters and eventually passed on to smpi_execute_benched that advances the clock of the executing rank by this estimate by entering a sleep state. The effect on simulation time for a small scenario is depicted in Figure\ref{fig:kernel_sampling}. On the one hand, this modification speeds up the simulation by orders of magnitude, especially when the matrix order grows. On the other hand, this kernel model leads to an optimistic estimation of the floprate. This may be caused by inaccuracies in our model as well as by the fact that the initial emulation is generally more sensitive to pre-emptions, \eg by the operating system, and therefore more likely to be pessimistic compared to a real execution.

HPL uses pseudo-randomly generated matrices that need to be setup every time HPL is executed. The time spent on this just as the validation of the computed result is not considered in the reported \si{\giga\flops} performance. We skip all the computations since we replaced them by a kernel model and therefore, result validation is meaningless. Since both phases do not have an impact on the reported performance, we can safely skip them.

In addition to the main computation kernels dgemm and dtrsm, we identified seven other BLAS functions through profiling as computationally expensive enough to justify a specific handling: dgemv, dswap, daxpy, dscal, dtrsv, dger and idamax. Similarly, a significant amount of time was spent in fifteen functions implemented in HPL: HPL_dlaswp*N, HPL_dlaswp*T, HPL_dlacpy and HPL_dlatcpy.

All of these functions are called during the LU factorization and hence impact the performance measured by HPL; however, because of the removal of the dgemm and dtrsm computations, they all operate on bogus data and hence also produce bogus data. We also determined through experiments that their impact on the performance prediction is minimal and hence modeled them for the sake of simplicity as being instantaneous.

Note that HPL implements an LU factorization with partial pivoting and a special treatment of the idamax function that returns the index of the first element equaling the maximum absolute value. Although we ignored the cost of this function as well, we set its return value to an arbitrary value to make the simulation fully deterministic. We confirmed that this modification is harmless in terms of performance prediction while it speeds up the simulation by an additional factor of \(\approx3\) to \(4\) on small (\(N=30,000\)) and even more on large scenarios.

As explained in Section\ref{sec:smpi}, when emulating an application with SMPI, all MPI processes are run within the same simulation process on a single node. The memory consumption of the simulation can therefore quickly reach several \si{\tera\byte} of RAM.

Yet, as we no longer operate on real data, storing the whole input matrix \(A\) is needless. However, since only a minimal portion of the code was modified, some functions may still read or write some parts of the matrix. It is thus not possible to simply remove the memory allocations of large data structures altogether. Instead, SMPI's SHARED_MALLOC mechanism can be used to share unimportant data structures between all ranks, minimizing the memory footprint.

The largest two allocated data structures in HPL are the input matrix A (with a size of typically several \si{\giga\byte} per process) and the panel which contains information about the sub-matrix currently being factorized. This sub-matrix typically occupies a few hundred \si{\mega\byte} per process. Although using the default SHARED_MALLOC mechanism works flawlessly for A, a more careful strategy needs to be used for the panel. Indeed, the panel is an intricate data structure with both \texttt{int}s (accounting for matrix indices, error codes, MPI tags, and pivoting information) and \texttt{double}s (corresponding to a copy of a sub-matrix of A). To optimize data transfers, HPL flattens this structure into a single allocation of \texttt{double}s (see Figure\ref{fig:panel_structure}). Using a fully shared memory allocation for the panel therefore leads to index corruption that results in classic invalid memory accesses as well as communication deadlocks, as processes may not send to or receive from the correct process. Since \texttt{int}s and \texttt{double}s are stored in non-contiguous parts of this flat allocation, it is therefore essential to have a mechanism that preserves the process-specific content. We have thus introduced the macro SMPI_PARTIAL_SHARED_MALLOC that works as follows: mem = SMPI_PARTIAL_SHARED_MALLOC(500, {27,42 , 100,200}, 2). In this example, 500 bytes are allocated in mem with the elements mem[27], …, mem[41] and mem[100], …, mem[199] being shared between processes (they are therefore generally completely corrupted) while all other elements remain private. To apply this to HPL's panel data­structure and partially share it between processes, we only had to modify a few lines.

Designating memory explicitly as private, shared or partially shared helps with both memory management and overall performance. As SMPI is internally aware of the memory's visibility, it can avoid calling memcopy when large messages containing shared segments are sent from one MPI rank to another. For fully private or partially shared segments, SMPI identifies and copies only those parts that are process-dependent (private) into the corresponding buffers on the receiver side.

HPL simulation times were considerably improved in our experiments because the panel as the most frequently transferred datastructure is partially shared with only a small part being private. The additional error introduced by this technique was negligible (below \SI{1}{\percent}) while the memory consumption was lowered significantly: for a matrix of order \(40,000\) and \(64\) MPI processes, the memory consumption decreased from about \SI{13.5}{\giga\byte} to less than \SI{40}{\mega\byte}.

HPL \texttt{malloc}s/\texttt{free}s panels in each iteration, with the size of the panel strictly decreasing from iteration to iteration. As we explained above, the partial sharing of panels requires many calls to mmap and introduces an overhead that makes these repeated allocations / frees become a bottleneck. Since the very first allocation can fit all subsequent panels, we modified HPL to allocate only the first panel and reuse it for subsequent iterations (see Figure\ref{fig:panel_reuse}).

We consider this optimization harmless with respect to simulation accuracy as the maximum additional error that we observed was always less than \SI{1}{\percent}. Simulation time is reduced significantly, albeit the reached speed-up is less impressive than for previous optimizations: For a very small matrix of order \(40,000\) and \(64\) MPI processes, the simulation time decreases by four seconds, from \SI{20.5}{\sec} to \SI{16.5}{\sec}. Responsible for this is a reduction of system time, namely from \SI{5.9}{\sec} to \SI{1.7}{\sec}. The number of page faults decreased from \(2\) million to \(0.2\) million, confirming the devastating effect these allocations/deallocations would have at scale.

We already explained in Section\ref{sec:appmodeling} that SMPI supports two mechanisms to keep local static and global variables private to each rank, even though they run in the same process. In this section, we discuss the impact of the choice.

• mmap When mmap is used, SMPI copies the data segment on startup for each rank into the heap. When control is transferred from one rank to another, the data segment is mmap'ed to the location of the other rank's copy on the heap. All ranks have hence the same addresses in the virtual address space at their disposition although mmap ensures they point to different physical addresses. This also means inevitably that caches must be flushed to ensure that no data of one rank leaks into the other rank, making mmap a rather expensive operation.

• dlopen With dlopen, copies of the global variables are still made but they are stored inside the data segment as opposed to the heap. When switching from one rank to another, the starting virtual address for the storage is readjusted rather than the target of the addresses. This means that each rank has distinct addresses for global variables. The main advantage of this approach is that caches do not need to be flushed as is the case for the mmap approach, because data consistency can always be guaranteed.

\noindent Impact of choice of mmap/dlopen The choice of mmap or dlopen influences the simulation time indirectly through its impact on system/user time and page faults, \eg for a matrix of order \(80,000\) and \(32\) MPI processes, the number of minor page faults drops from \num{4412047} (with mmap) to \num{6880} (with dlopen). This results in a reduction of system time from \SI{10.64}{\sec} (out of \SI{51.47}{\sec} in total) to \SI{2.12}{\sec}. Obviously, the larger the matrix and the number of processes, the larger the number of context switch during the simulation, and thus the higher the gain.

For larger matrix orders (\ie \(N\) larger than a few hundred thousand), the performance of the simulation quickly deteriorates as the memory consumption rises rapidly.

We explained already how we fold the memory in order to reduce the physical memory usage. The virtual memory, on the other hand, is still allocated for every process since the allocation calls are still executed. Without a reduction of allocated virtual addresses, the page table rapidly becomes too large to fit in a single node. More precisely, the size of the page table containing pages of size \SI{4}{\kibi\byte} can be computed as:

This means that the addresses in the page table for a matrix of order \(N=4,000,000\) consume \(PT_{size}(4,000,000) = \num{2.5e11}\) bytes, \ie \SI{250}{\giga\byte} on a system where double-precision floating-point numbers and addresses take 8 bytes. Thankfully, the x86-64 architecture supports several page sizes, known as ``huge pages'' in Linux. Typically, these pages are around \SI{2}{\mebi\byte} (instead of \SI{4}{\kibi\byte}), although other sizes (\SIrange{2}{256}{\mebi\byte}) are possible as well. Changing the page size requires administrator (root) privileges as the Linux kernel support for hugepages needs to be activated and a hugetlbfs file system must be mounted. After at least one huge page has been allocated, the path of the allocated file system can then be passed on to SimGrid. Setting the page size to \SI{2}{\mebi\byte} reduces drastically the page table size. For example, for a matrix of order \(N=4,000,000\), it shrinks from \SI{250}{\giga\byte} to \SI{0.488}{\giga\byte}.