Randomized projection methods have been shown to be very efficient at computing low-rank approximations (LRA) of large matrices. In this work, we investigate the design and development of such methods capable of exploiting recent mixed precision accelerators like GPUs equipped with tensor core units. We combine three new ideas to exploit mixed precision arithmetic in randomized LRA. The first is to perform the matrix multiplication with mixed precision fp16/fp32 tensor cores. The second is to use CholeskyQR orthonormalization, which is much faster on GPUs, while mitigating its numerical instability by using fp64 arithmetic. The third is to use a recently proposed iterative refinement method for LRA to improve the accuracy of the LRA by calling it twice. We implement the proposed approach on various GPU architectures and analyze its performance and accuracy. We compare with a standard randomized LRA entirely in fp32 arithmetic, which achieves an average accuracy of order 10^-4 . Our results show that our approach without refinement is up to 8× faster, with an average accuracy of order 10^-2 , which may be acceptable for some applications. Otherwise, we show that using refinement significantly improves the accuracy to an average of order 10^-5 , while remaining up to 2.2× faster than the standard fp32 randomized LRA. This work illustrates the convergence of approximate computing techniques by combining low-rank approximations, randomization, mixed precision arithmetic, and GPU acceleration.

Stencil calculations and matrix-free Krylov subspace solvers represent important components of many scientific computing applications. In these solvers, stencil applications are often the dominant part of the computation; an efficient parallel implementation of the kernel is therefore crucial to reduce the time to solution. Inspired by polynomial preconditioning, we remove upper bounds on the arithmetic intensity of the Krylov subspace building block by replacing the matrix with a higher-degree matrix polynomial. Using the latest state-of-the-art stencil compiler programs with temporal blocking, reduced memory bandwidth usage and, consequently, better utilization of SIMD vectorization and thus speedup on modern hardware, we are able to obtain performance improvements for higher polynomial degrees than simpler cache-blocking approaches have yielded in the past, demonstrating the new appeal of polynomial techniques on emerging architectures. We present results in a shared-memory environment and an extension to a distributed-memory environment with local shared memory.

Gaussian elimination is a canonical linear algebra procedure for solving linear systems of equations. In the last few years, the algorithm has received a lot of attention in an attempt to improve its parallel performance. This article surveys recent developments in parallel implementations of Gaussian elimination for shared memory architecture. Five different flavors are investigated. Three of them are based on different strategies for pivoting: partial pivoting, incremental pivoting, and tournament pivoting. The fourth one replaces pivoting with the Partial Random Butterfly Transformation, and finally, an implementation without pivoting is used as a performance baseline. The technique of iterative refinement is applied to recover numerical accuracy when necessary. All parallel implementations are produced using dynamic, superscalar, runtime scheduling and tile matrix layout. Results on two multisocket multicore systems are presented. Performance and numerical accuracy is analyzed. Copyright © 2014 John Wiley & Sons, Ltd.

Graphics processing units (GPUs) brought huge performance improvements in the scientific and numerical fields. We present an efficient hybrid CPU/GPU approach that is portable, dynamically and efficiently balances the workload between the CPUs and the GPUs, and avoidsdata transfer bottlenecks that are frequently present in numerical algorithms. Our approach determines the amount of initial work to assign to the CPUs before the execution, and then dynamically balances workloads during the execution. Then, we present a theoretical model to guide the choice of the initial amount of work for the CPUs. The validation of our model allows our approach to self-adapt on any architecture using the manufacturer's characteristics of the underlying machine. We illustrate our method for the LU factorization. For this case, we show that the use of our approach combined with a communication avoiding LU algorithm is efficient. For example, our experiments on a 24 cores AMD opteron 6172 show that by adding one GPU (Tesla S2050) we accelerate LU up to 2.4× compared to the corresponding routine in MKL using 24 cores. The comparisons with MAGMA also show significant improvements.

Performance irregularities on massively parallel processors lead to load imbalances and a significant loss of performance. Multi-core nodes suggest a promising way to redistribute work within a node, thus mitigating performance irregularities. However, there exists a non-trivial cost to redistributing work, and associated data, across cores. We investigate how work can be equitably distributed across cores without significantly disturbing data locality, and without incurring significant scheduling overhead. Towards this end, we design a series of scheduling strategies and tuning mechanisms; our foundational technique is an intelligent blending of static and dynamic scheduling. We also implement a basic runtime system and library to minimize programmer effort in applying these strategies. Our techniques provide 28.16% performance gains over static scheduling and 17.13% gains over guided scheduling for a widely used regular mesh benchmark, and 44.45% gains over static scheduling and 13.06% gains over guided scheduling for an n-body simulation, both on 1024 nodes. I. MAINTAINING LOAD BALANCE AND LOCALITY Load imbalances come from the application (typically coarse-grained, persistent) or the architecture (typically finegrained, transient). Figure 1 shows the timeline of a single execution of an already highly optimized dense matrix factorization code, which we refer to as CALU, on a multicore node. In the timeline, load imbalances across cores occur due to the white areas of the timeline, which indicate idle times on a core. Work could be re-distributed across cores during execution time to reduce these idle times. The purpose of Figure 1 is to show that load balancing is an important and non-trivial problem for HPC applications, not just due to increasing complexity of applications, but also due to increasing complexity of architectures; load balancing is important even for computations that have conventionally been considered to not need load balancing, with this dense Communication-avoiding LU factorization being one example. A basic load balancing technique, e.g., dynamic scheduling [1], can incur a cost of data movement to re-distribute work across cores, and this cost can significantly degrade performance of an application [2]. We also note that the significant costs of cache misses only increase for a multi-core node with a larger number of cores [3]). To minimize this cost, localityaware scheduling [1], [2], [4] can be used. Figure 2a shows the application execution timeline for CALU when localityaware scheduling (called “scheduling with locality tags”) is used, and Figure 2b shows the application execution timeline when CALU uses a different type of locality-aware scheduling (called “scheduling with locality tags and stealing”). The basic locality-aware scheduling strategy and the optimized locality-aware scheduling strategy is explained in [2], and its application to the CALU code is discussed in [4]. The green rectangular block is the execution timeline of CALU. The cache misses during execution of CALU code (with the locality-aware scheduling technique applied) are in the red line graph below this execution timeline. In this case, there are several cache misses, indicating loss of locality and causing large performance degradations due to the significant costs of cache misses. The purpose of these graphs is to explain the problem that the strategies developed for this particular example may not be useful in another application or architecture, where the balance needed between locality and load balance may be different. With this, to balance the tradeoff between load balance and locality, our approach is to use an intelligent blend of static and dynamic scheduling, with the proportion of static and dynamic scheduling used being carefully tuned. The two diagrams in the bottom of this section, which each contain three rectangular solids, show how our strategies are applied to an MPI+OpenMP 3D stencil code. The diagram on the left shows a single timestep of a 3D stencil code written in MPI+OpenMP, where work within each MPI process is statically scheduled across cores. The diagram on the right shows the application of our lightweight scheduling technique to this 3D stencil code, where the basic static/dynamic scheduling strategy is used. The strategy is defined to the left of these two diagrams. We note that the static fraction is the scheduling parameter of mixed static/dynamic scheduling. We call this strategy mixed static/dynamic scheduling, or Lightweight Scheduling. II. IMPLEMENTATION AND SCHEDULER OPTIMIZATION Optimizations over the mixed static/dynamic scheduling strategy can help to further balance the tradeoff between load balance and locality. Examples, as shown in the bottom-left quadrant, are constrained staggered static/dynamic scheduling and mixed static/dynamic scheduling with variable-sized tasks. Each strategy allows for new opportunities to help find the best balance between load balance and locality. Additionally, we must define the methodology for tuning the balance between load balance and locality. We can do this through experimental tuning, where we run the scheduler with different scheduler parameter values, and use the best-performing parameter values during application execution. We can also adjust the scheduler parameters at runtime, based on information from previous application timesteps. Finally, we can use modelguided optimization, developed through performance modeling and theoretical analysis, during application execution to prune the search space for efficiently finding the best-performing scheduler parameters. An illustration of the slack-conscious scheduling strategy, along with its corresponding formula used in its mechanism of model-guided optimization to tune the scheduler parameters, is shown in the right column. The formula shown is used to tune the static fraction per MPI process, based on MPI slack. We note that we assume one MPI process is assigned to a node of a cluster, as is done in [2], [3]. The slack is shown in orange in the figure, and is denoted by S in the formula shown. Figure 5 shows the code for a conjugate gradient (CG) computation transformed to use our technique. The static fraction used for the threaded computation region shown, enclosed by “FORALL BEGIN(...)” and “FORALL END(...)”, is calculated through the use of the function call predict_static_fraction(). III. WITHIN-NODE RESULTS AND ANALYSIS The purpose of Figure 3 is to show the competitive performance of the mixed static/dynamic scheduled version of CALU for widely used implementations of LU factorization such as MKL and PLASMA. Our mixed static/dynamic approach provides 30.34% gains over the PLASMA library and 34.46% gains over the MKL library, and is further discussed in [3]. The timeline shown in Figure 4a shows the CALU computation as it is applied to static/dynamic scheduling. This timeline should be compared to the timeline in Figure 1. This shows how mixed static/dynamic scheduling improves both load imbalance. The histograms in Figure 4b show the distribution of execution times for 5000 independent executions of the CALU code on a single node of an Intel Westmere cluster. The distribution of execution times show impact to scalability, as discussed in [3]; having both absolute performance and performance variations are necessary for obtaining performance gains at large-scale. The performance of the executions in fully static scheduling is multi-modal, showing the impact of different load imbalances. The load imbalances come from both the application and the architecture. The performance variations are small for dynamic scheduling, but because of the overheads that dynamic scheduling incurs, performance degradation increases. Both the performance variations and absolute performance are least when mixed static/dynamic scheduling is used. Table 1 and Table 2 show the performance improvement of our scheduling technique over OpenMP static scheduling on the NAS benchmarks for an Intel Westmere 16-core node and IBM BG/Q 16-core node, respectively. As can be seen in Table 1, the gains for the CG benchmark are large on the Intel Westmere machine due to the scheduler’s handling of application load imbalance of NAS CG along with load imbalance due to performance irregularities sourcing from OS noise. As seen in Table 2, while the BG/Q machine has low noise [3], the scheduler still acheives significant gains for CG due to its ability to handle the application load imbalance of CG. IV. ACROSS-NODE RESULTS AND SCALABILITY Figure 6a shows the performance for different strategies as applied to the MPI+OpenMP regular mesh code SNAP [5]. Figure 6b shows performance for different scheduling strategies applied to an MPI+OpenMP n-body particle galaxy simulation code from the Rebound application [6], as we increase the number of nodes used in a cluster of SMPs. The speedups over static scheduling are shown. We focus on the results for figure 6b. The OpenMP dynamic scheduling strategy does only slightly better at small node counts, but helps at larger node counts where load imbalance across cores is more. The OpenMP guided scheduling strategy [7] does better than the dynamic scheduling strategy. The uSched strategy, which is the basic mixed static/dynamic scheduling strategy illustrated in the end of the top-right quadrant of the poster, greatly reduces the scheduling cost. This strategy helps both at small and large node counts. We note that uSched performs better than guided scheduling, due to guided scheduling not being able to maintain locality across timesteps, and due to its significant dequeue overheads for the large tasks. Additional results are in [3]. The uSched strategy provides 25.52% performance gains over static scheduling at large node counts. The callsite strategy is the optimization to the uSched strategy that takes into account MPI slack [8] to reduce scheduling overheads, and is illustr

Gaussian elimination is a canonical linear algebra procedure for solving linear systems of equations. In the last few years, the algorithm received a lot of attention in an attempt to improve its parallel performance. This article surveys recent developments in parallel implementations of the Gaussian elimination. Five different flavors are investigated. Three of them are based on different strategies for pivoting: partial pivoting, incremental pivoting, and tournament pivoting. The fourth one replaces pivoting with the Random Butterfly Transformation, and finally, an implementation without pivoting is used as a performance baseline. The technique of iterative refinement is applied to recover numerical accuracy when necessary. All parallel implementations are produced using dynamic, superscalar, runtime scheduling and tile matrix layout. Results on two multi-socket multicore systems are presented. Performance and numerical accuracy is analyzed.

We study several solvers for the solution of general linear systems where the main objective is to reduce the communication overhead due to pivoting. We first describe two existing algorithms for the LU factorization on hybrid CPU/GPU architectures. The first one is based on partial pivoting and the second uses a random preconditioning of the original matrix to avoid pivoting. Then we introduce a solver where the panel factorization is performed using a communication-avoiding pivoting heuristic while the update of the trailing submatrix is performed by the GPU. We provide performance comparisons and tests on accuracy for these solvers on current hybrid multicore-GPU parallel machines.

Due to the evolution of massively parallel computers towards deeper levels of parallelism and memory hierarchy, and due to the exponentially increasing ratio of the time required to transfer data, either through the memory hierarchy or between different compute units, to the time required to compute floating point operations, the algorithms are confronted with two challenges. They need not only to be able to exploit multiple levels of parallelism, but also to reduce the communication between the compute units at each level of the hierarchy of parallelism and between the different levels of the memory hierarchy. In this paper we present an algorithm for performing the LU factorization of dense matrices that is suitable for computer systems with two levels of parallelism. This algorithm is able to minimize both the volume of communication and the number of messages transferred at every level of the two-level hierarchy of parallelism. We present its implementation for a cluster of multicore processors based on MPI and Pthreads. We show that this implementation leads to a better performance than routines implementing the LU factorization in well-known numerical libraries. For matrices that are tall and skinny, that is they have many more rows than columns, our algorithm outperforms the corresponding algorithm from ScaLAPACK by a factor of 4.5 on a cluster of 32 nodes, each node having two quad-core Intel Xeon EMT64 processors.

The growing demand for scalable and effective scientific and numerical libraries on multicore architectures forces hardware manufacturers to design solutions that improve both the processor speed and transfer rates between their memory hierarchies. Several studies show that these improvement factors are disproportionate and may vary widely from one architecture to another and then have a strong impact on the tuning and the performance prediction of numerical libraries. In this paper, we analyze the communication and performance of some routines in well known libraries on different architectures and we establish a relation model between hardware parameters and performance. We focus on the LU factorization, which is one the most popular algorithms in the scientific field, therefore also used as a benchmark, e.g., the HPL benchmark to rank the TOP500 supercomputers. Our experiments in terms of hardware counter measurements allow us to predict the performance behavior of numerical algorithms (LU in particular) on different architectures.

We present the use of a hybrid static/dynamic scheduling strategy of the task dependency graph for direct methods used in dense numerical linear algebra. This strategy provides a balance of data locality, load balance, and low dequeue overhead. We show that the usage of this scheduling in communication avoiding dense factorization leads to significant performance gains. On a 48 core AMD Opteron NUMA machine, our experiments show that we can achieve up to 64% improvement over a version of CALU that uses fully dynamic scheduling, and up to 30% improvement over the version of CALU that uses fully static scheduling. On a 16-core Intel Xeon machine, our hybrid static/dynamic scheduling approach is up to 8% faster than the version of CALU that uses a fully static scheduling or fully dynamic scheduling. Our algorithm leads to speedups over the corresponding routines for computing LU factorization in well known libraries. On the 48 core AMD NUMA machine, our best implementation is up to 110% faster than MKL, while on the 16 core Intel Xeon machine, it is up to 82% faster than MKL. Our approach also shows significant speedups compared with PLASMA on both of these systems.

In this paper we discuss a hypergraph-based unsymmetric nested dissection (HUND) ordering for reducing the fill-in incurred during Gaussian elimination. It has several important properties. It takes a global perspective of the entire matrix, as opposed to local heuristics. It takes into account the asymmetry of the input matrix by using a hypergraph to represent its structure. It is suitable for performing Gaussian elimination in parallel, with partial pivoting. This is possible because the row permutations performed due to partial pivoting do not destroy the column separators identified by the nested dissection approach. The hypergraph nested dissection approach is essentially equivalent to graph nested dissection on the matrix $A^{T}A$, but we need only the original matrix $A$ and never form the usually denser matrix $A^{T}A$. The usage of hypergraphs in our approach is fairly standard, and HUND can be implemented by calling an existing hypergraph partitioner that uses recursive bisection. Our implementation uses local reordering constrained column approximate minimum degree (CCOLAMD) to further improve the ordering. We also explain how weighted matching (HSL routine MC64) can be used in this context. Experimental results on 27 medium and large size matrices with highly unsymmetric structures compare our approach to four other well-known reordering algorithms. The results show that it provides a robust reordering algorithm, in the sense that it is the best or close to the best (often within 10%) of all the other methods, in particular on matrices with highly unsymmetric structures.

In this paper we study algorithms for performing the LU and QR factorizations of dense matrices. Recently, two communication optimal algorithms have been introduced for distributed memory architectures, referred to as communication avoiding CALU and CAQR. In this paper we discuss two algorithms based on CAQR and CALU that are adapted to multicore architectures. They combine ideas to reduce communication from communication avoiding algorithms with asynchronism and dynamic task scheduling. For matrices that are tall and skinny, that is, they have many more rows than columns, the two algorithms outperform the corresponding algorithms from Intel MKL vendor library on a dual-socket, quad-core machine based on Intel Xeon EMT64 processor and on a four-socket, quad-core machine based on AMD Opteron processor. For these matrices, multithreaded CALU outperforms the corresponding routine dgetrf from Intel MKL library up to a factor of 2.3 and the corresponding routine dgetrf from ACML library up to a factor of 5, while multithreaded CAQR outperforms by a factor of 5.3 the corresponding dgeqrf routine from MKL library.