Sparse general matrix-matrix multiplication (SpGEMM) is a critical operation in many applications.Current multithreaded implementations are based on Gustavson's algorithm and often perform poorly on large matrices due to limited cache reuse by the accumulators.We present MAGNUS (Matrix Algebra for Gigantic NUmerical Systems), a novel algorithm to maximize data locality in SpGEMM.To generate locality, MAGNUS reorders the intermediate product into discrete cache-friendly chunks using a two-level hierarchical approach.The accumulator is applied to each chunk, where the chunk size is chosen such that the accumulator is cacheefficient.MAGNUS is input-and system-aware: based on the matrix characteristics and target system specifications, the optimal number of chunks is computed by minimizing the storage cost of the necessary data structures.MAGNUS allows for a hybrid accumulation strategy in which each chunk uses a different accumulator based on an input threshold.We consider two accumulators: an AVX-512 vectorized bitonic sorting algorithm and classical dense accumulation.An OpenMP implementation of MAGNUS is compared with several baselines, including Intel MKL, for a variety of different matrices on three Intel architectures.For matrices from the SuiteSparse collection, MAGNUS is faster than all the baselines in most cases and is often an order of magnitude faster than at least one baseline.For massive random matrices, MAGNUS scales to the largest matrix sizes, while the baselines do not.Furthermore, MAGNUS is close to the optimal bound for these matrices, regardless of the matrix size, structure, and density.

Reducing synchronization in iterative methods for solving large sparse linear systems may become one of the most important goals for such solvers on exascale computers. Research in asynchronous iterative methods has primarily considered basic iterative methods. In this paper, we examine how multigrid methods can be executed asynchronously. We present models of asynchronous additive multigrid methods, and use these models to study the convergence properties of these methods. We also introduce two parallel algorithms for implementing asynchronous additive multigrid, the global-res and local-res algorithms. These two algorithms differ in how the fine grid residual is computed, where local-res requires less computation than global-res but converges more slowly. We compare two types of asynchronous additive multigrid methods: the asynchronous fast adaptive composite grid method with smoothing (AFACx) and additive variants of the classical multiplicative method (Multadd). We implement asynchronous versions of Multadd and AFACx in OpenMP and generate the prolongation and coarse grid matrices using the BoomerAMG package. Our experimental results show that asynchronous multigrid can exhibit grid-size independent convergence and can be faster than classical multigrid in terms of solve wall-clock time. We also show that asynchronous smoothing is the best choice of smoother for our test cases, even when only one smoothing sweep is used.

Asynchronous iterative methods for solving linear systems have been gaining attention due to the high cost of synchronization points in massively parallel codes. Since future parallel computers will likely achieve exascale performance, synchronization may become the primary bottleneck. Historically, theory on asynchronous iterative methods has focused on proving that an asynchronous version of a fixed point method will converge. Additionally, some theory has shown that asynchronous methods can be faster, which has been supported by shared memory experiments. In this paper, we introduce a new way to model asynchronous Jacobi using propagation matrices, which are similar in concept to iteration matrices. With this model, we show that asynchronous Jacobi can reduce the error even if some processes are delayed for a long period of time, which could be due to a hardware malfunction or a large imbalance. We also show that asynchronous Jacobi can converge when synchronous Jacobi does not. We compare model results to shared and distributed memory implementation results, and show that in practice, asynchronous Jacobi's convergence rate improves as we increase the number of processes.

We present a new algorithm, the Distributed Southwell method, as a competitor to Block Jacobi for preconditioning and multigrid smoothing. It is based on the Southwell iterative method, which is sequential, where only the equation with the largest residual is relaxed per iteration. The Parallel Southwell method extends this idea by relaxing equation i if it has the largest residual among all the equations coupled to variable i. Since communication is required for processes to exchange residuals, this method in distributed memory can be expensive. Distributed Southwell uses a novel scheme to reduce this communication of residuals while avoiding deadlock. Using test problems from the SuiteSparse Matrix Collection, we show that Distributed Southwell requires less communication to reach the same accuracy when compared to Parallel Southwell. Additionally, we show that the convergence of Distributed Southwell does not degrade like that of Block Jacobi when the number of processes is increased.