We introduce a data distribution scheme for H-matrices and a distributedmemory algorithm for H-matrix-vector multiplication. Our data distribution scheme avoids an expensive Ω(P2) scheduling procedure used in previous work, where P is the number of processes, while data balancing is well-preserved. Based on the data distribution, our distributed-memory algorithm evenly distributes all computations among P processes and adopts a novel tree-communication algorithm to reduce the latency cost. The overall complexity of our algorithm is O ( N logN P +αlogP+βlog 2 P ) for H-matrices under weak admissibility condition, where N is the matrix size, α denotes the latency, and β denotes the inverse bandwidth. Numerically, our algorithm is applied to address both twoand three-dimensional problems of various sizes among various numbers of processes. On thousands of processes, good parallel efficiency is still observed. AMS subject classifications: 65F99, 65Y05

Determinantal point processes (DPPs) were introduced by Macchi (Macchi 1975 Adv. Appl. Probab. 7, 83-122) as a model for repulsive (fermionic) particle distributions. But their recent popularization is largely due to their usefulness for encouraging diversity in the final stage of a recommender system (Kulesza & Taskar 2012 Found. Trends Mach. Learn. 5, 123-286). The standard sampling scheme for finite DPPs is a spectral decomposition followed by an equivalent of a randomly diagonally pivoted Cholesky factorization of an orthogonal projection, which is only applicable to Hermitian kernels and has an expensive set-up cost. Researchers Launay et al. 2018 (http://arxiv.org/abs/1802.08429); Chen & Zhang 2018 NeurIPS (https://papers.nips.cc/paper/7805-fast-greedy-map-inference-for-determinantal-point-process-to-improve-recommendation-diversity.pdf) have begun to connect DPP sampling to LDLH factorizations as a means of avoiding the initial spectral decomposition, but existing approaches have only outperformed the spectral decomposition approach in special circumstances, where the number of kept modes is a small percentage of the ground set size. This article proves that trivial modifications of LU and LDLH factorizations yield efficient direct sampling schemes for non-Hermitian and Hermitian DPP kernels, respectively. Furthermore, it is experimentally shown that even dynamically scheduled, shared-memory parallelizations of high-performance dense and sparse-direct factorizations can be trivially modified to yield DPP sampling schemes with essentially identical performance. The software developed as part of this research, Catamari (hodgestar.com/catamari) is released under the Mozilla Public License v.2.0. It contains header-only, C++14 plus OpenMP 4.0 implementations of dense and sparse-direct, Hermitian and non-Hermitian DPP samplers. This article is part of a discussion meeting issue 'Numerical algorithms for high-performance computational science'.

We expose a systematic approach for developing distributed-memory parallel matrix-matrix multiplication algorithms. The journey starts with a description of how matrices are distributed to meshes of nodes (e.g., MPI processes), relates these distributions to scalable parallel implementation of matrix-vector multiplication and rank-1 update, continues on to reveal a family of matrix-matrix multiplication algorithms that view the nodes as a two-dimensional (2D) mesh, and finishes with extending these 2D algorithms to so-called three-dimensional (3D) algorithms that view the nodes as a 3D mesh. A cost analysis shows that the 3D algorithms can attain the (order of magnitude) lower bound for the cost of communication. The paper introduces a taxonomy for the resulting family of algorithms and explains how all algorithms have merit depending on parameters such as the sizes of the matrices and architecture parameters. The techniques described in this paper are at the heart of the Elemental distributed-memory linear algebra library. Performance results from implementation within and with this library are given on a representative distributed-memory architecture, the IBM Blue Gene/P supercomputer.

Multi-shift triangular solves are basic linear algebra calculations with applications in eigenvector and pseudospectra computation. We propose blocked algorithms that efficiently exploit Level 3 BLAS to perform multi-shift triangular solves and safe multi-shift triangular solves. Numerical experiments indicate that computing triangular eigenvectors with a safe multi-shift triangular solve achieves speedups by a factor of 60 relative to LAPACK. This algorithm accelerates the calculation of general eigenvectors threefold. When using multi-shift triangular solves to compute pseudospectra, we report ninefold speedups relative to EigTool.

The butterfly algorithm is a fast algorithm which approximatelyevaluates a discrete analogue of the integral transform$\int_{\mathbb{R}^d} K(x,y) g(y) dy$ at large numbers of target pointswhen the kernel, $K(x,y)$, is approximately low-rank when restricted tosubdomains satisfying a certain simple geometric condition.In $d$ dimensions with $O(N^d)$ quasi-uniformly distributed source and target points,when each appropriate submatrix of $K$ is approximately rank-$r$, the running time ofthe algorithm is at most $O(r^2 N^d \log N)$.A parallelization of the butterfly algorithm is introduced which, assuming amessage latency of $\alpha$ and per-process inverse bandwidth of $\beta$,executes in at most$O(r^2 \frac{N^d}{p} \log N + (\beta r\frac{N^d}{p}+\alpha)\log p)$time using $p$ processes.This parallel algorithm was then instantiated in the form of the open-source\textttDistButterfly library for the special case where$K(x,y)=\exp(i \Phi(x,y))$, where $\Phi(x,y)$ is a black-box, sufficientlysmooth, real-valued phase function .Experiments on Blue Gene/Q demonstrate impressive strong-scaling results forimportant classes of phase functions. Using quasi-uniform sources,hyperbolic Radon transforms, and an analogue of a three-dimensional generalized Radon transformwere, respectively, observed to strong-scale from 1-node/16-cores up to1024-nodes/16,384-cores with greater than 90% and82% efficiency, respectively.

We present a parallel preconditioning method for the iterative solution of the time-harmonic elastic wave equation which makes use of higher-order spectral elements to reduce pollution error. In particular, the method leverages perfectly matched layer boundary conditions to efficiently approximate the Schur complement matrices of a block LDLT factorization. Roth sequential and parallel versions of the algorithm are discussed and results for large-scale problems from exploration geophysics are presented.

This paper introduces a parallel directional fast multipole method (FMM) for solving N-body problems with highly oscillatory kernels, with a focus on the Helmholtz kernel in three dimensions. This class of oscillatory kernels requires a more restrictive low-rank criterion than that of the low-frequency regime, and thus effective parallelizations must adapt to the modified data dependencies. We propose a simple partition at a fixed level of the octree and show that, if the partitions are properly balanced between p processes, the overall runtime is essentially O(N log N/p+ p). By the structure of the low-rank criterion, we are able to avoid communication at the top of the octree. We demonstrate the effectiveness of our parallelization on several challenging models.

Parallelizing dense matrix computations to distributed memory architectures is a well-studied subject and generally considered to be among the best understood domains of parallel computing. Two packages, developed in the mid 1990s, still enjoy regular use: ScaLAPACK and PLAPACK. With the advent of many-core architectures, which may very well take the shape of distributed memory architectures within a single processor, these packages must be revisited since the traditional MPI-based approaches will likely need to be extended. Thus, this is a good time to review lessons learned since the introduction of these two packages and to propose a simple yet effective alternative. Preliminary performance results show the new solution achieves competitive, if not superior, performance on large clusters.

Parallelizing dense matrix computations to distributed memory architectures is a well-studied subject and generally considered to be among the best understood domains of parallel computing. Two packages, developed in the mid 1990s, still enjoy regular use: ScaLAPACK and PLAPACK. With the advent of many-core architectures, which may very well take the shape of distributed memory architectures within a single processor, these packages must be revisited since the traditional MPI-based approaches will likely need to be extended. Thus, this is a good time to review lessons learned since the introduction of these two packages and to propose a simple yet effective alternative. Preliminary performance results show the new solution achieves competitive, if not superior, performance on large clusters.

A parallelization of a sweeping preconditioner for three-dimensional Helmholtz equations without large cavities is introduced and benchmarked for several challenging velocity models. The setup and application costs of the sequential preconditioner are shown to be $O(\gamma^2 N^{4/3})$ and $O(\gamma N \log N)$, where $\gamma(\omega)$ denotes the modestly frequency-dependent number of grid points per perfectly matched layer. Several computational and memory improvements are introduced relative to using black-box sparse-direct solvers for the auxiliary problems, and competitive runtimes and iteration counts are reported for high-frequency problems distributed over thousands of cores. Two open-source packages are released along with this paper: Parallel Sweeping Preconditioner (PSP) and the underlying distributed multifrontal solver, Clique.

To implement dense linear algebra algorithms for distributed-memory computers, an expert applies knowledge of the domain, the target architecture, and how to parallelize common operations. This is often a rote process that becomes tedious for a large collection of algorithms. We have developed a way to encode this expert knowledge such that it can be applied by a system to generate mechanically the same (and sometimes better) highly-optimized code that an expert creates by hand. This paper illustrates how we have encoded a subset of this knowledge and how our system applies it and searches a space of generated implementations automatically.

The efforts of an expert to parallelize and optimize a dense linear algebra algorithm for distributed-memory targets are largely mechanical and repetitive. We demonstrate that these efforts can be encoded and automatically applied to obviate the manual implementation of many algorithms in high-performance code.

PreviousNext No AccessSEG Technical Program Expanded Abstracts 2012A parallel sweeping preconditioner for frequency-domain seismic wave propagationAuthors: Jack PoulsonBjörn EngquistSiwei LiLexing YingJack PoulsonThe University of Texas at AustinSearch for more papers by this author, Björn EngquistThe University of Texas at AustinSearch for more papers by this author, Siwei LiThe University of Texas at AustinSearch for more papers by this author, and Lexing YingThe University of Texas at AustinSearch for more papers by this authorhttps://doi.org/10.1190/segam2012-0765.1 SectionsSupplemental MaterialAboutPDF/ePub ToolsAdd to favoritesDownload CitationsTrack CitationsPermissions ShareFacebookTwitterLinked InRedditEmail Abstract We present a parallel implementation of Engquist and Ying's sweeping preconditioner, which exploits radiation boundary conditions in order to form an approximate block LDLT factorization of the Helmholtz operator with only O(N4/3) work and an application (and memory) cost of only O(N logN). The approximate factorization is then used as a preconditioner for GMRES, and we show that essentially O(1) iterations are required for convergence, even for the full SEG/EAGE over-thrust model at 30 Hz. In particular, we demonstrate the solution of said problem in a mere 15 minutes on 8192 cores of TACC's Lonestar, which may be the largest-scale 3D heterogeneous Helmholtz calculation to date. Generalizations of our parallel strategy are also briefly discussed for time-harmonic linear elasticity and Maxwell's equations. Permalink: https://doi.org/10.1190/segam2012-0765.1FiguresReferencesRelatedDetails SEG Technical Program Expanded Abstracts 2012ISSN (print):1052-3812 ISSN (online):1949-4645Copyright: 2012 Pages: 4609 Publisher:Society of Exploration Geophysicists HistoryPublished: 25 Oct 2012 CITATION INFORMATION Jack Poulson, Björn Engquist, Siwei Li, and Lexing Ying, (2012), "A parallel sweeping preconditioner for frequency-domain seismic wave propagation," SEG Technical Program Expanded Abstracts : 1-6. https://doi.org/10.1190/segam2012-0765.1 Plain-Language Summary PDF DownloadLoading ...

In spite of an extensive literature on fast algorithms for synthetic aperture radar (SAR) imaging, it is not currently known if it is possible to accurately form an image from $N$ data points in provable near-linear time complexity. This paper seeks to close this gap by proposing an algorithm which runs in complexity $O(N \log N \log(1/\epsilon))$ without making the far-field approximation or imposing the beam pattern approximation required by time-domain backprojection, with $\epsilon$ the desired pixelwise accuracy. It is based on the butterfly scheme, which unlike the FFT works for vastly more general oscillatory integrals than the discrete Fourier transform. A complete error analysis is provided: the rigorous complexity bound has additional powers of $\log N$ and $\log(1/\epsilon)$ that are not observed in practice.

SUMMARYA message passing, distributed‐memory parallel computer on a chip is one possible design for future, many‐core architectures. We discuss initial experiences with the Intel Single‐chip Cloud Computer research processor, which is a prototype architecture that incorporates 48 cores on a single die that can communicate via a small, shared, on‐die buffer. The experiment is to port a state‐of‐the‐art, distributed‐memory, dense matrix library, Elemental, to this architecture and gain insight from the experience. We show that programmability addressed by this library, especially the proper abstraction for collective communication, greatly aids the porting effort. This enables us to support a wide range of functionality with limited changes to the library code. Copyright © 2011 John Wiley & Sons, Ltd.

A message passing, distributed-memory parallel computer on a chip is one possible design for future, many-core architectures. We discuss initial experiences with the Intel Single-chip Cloud Computer research processor, which is a prototype architecture that incorporates 48 cores on a single die that can communicate via a small, shared, ondie buffer. The experiment is to port a state-of-the-art, distributed-memory, dense matrix library, Elemental, to this architecture and gain insight from the experience. We show that programmability addressed by this library, especially the proper abstraction for collective communication, greatly aids the porting effort. This enables us to support a wide range of functionality with limited changes to the library code.

PreviousNext No AccessSEG Technical Program Expanded Abstracts 2011Sweeping preconditioner for the 3D Helmholtz equationAuthors: Björn EngquistJack PoulsonLexing YingBjörn EngquistThe University of Texas at AustinSearch for more papers by this author, Jack PoulsonThe University of Texas at AustinSearch for more papers by this author, and Lexing YingThe University of Texas at AustinSearch for more papers by this authorhttps://doi.org/10.1190/1.3627851 SectionsSupplemental MaterialAboutPDF/ePub ToolsAdd to favoritesDownload CitationsTrack CitationsPermissions ShareFacebookTwitterLinked InRedditEmail Abstract The Helmholtz equation describes wave propagation in the frequency domain and, as such, can be used for seismic imaging and full waveform inversion. We present two novel preconditioners for the efficient solution of the Helmholtz equation in three dimensions. Both methods follow the general structure of constructing an approximate LDLt factorization by eliminating the unknowns layer by layer starting from an absorbing layer or boundary condition. In the first approach, we represent the Schur complement matrices of the factorization in the hierarchical matrix framework. In the second approach, applying each Schur complement matrix is equivalent to solving a quasi‐2D problem with the multifrontal method. These preconditioners have linear application cost, and the preconditioned iterative solvers converge in a number of iterations that is essentially independent of the number of unknowns or the frequency. Numerical results on realistic 3D seismic models to confirm the efficiency of these methods.Permalink: https://doi.org/10.1190/1.3627851FiguresReferencesRelatedDetailsCited ByA Unified 2D/3D Large-Scale Software Environment for Nonlinear Inverse ProblemsACM Transactions on Mathematical Software, Vol. 45, No. 1Restrictions on the Use of Sweeping Type Preconditioners for Helmholtz Problems5 January 2019Sweeping preconditioners for elastic wave propagation with spectral element methods20 February 2014 | ESAIM: Mathematical Modelling and Numerical Analysis, Vol. 48, No. 2A parallel sweeping preconditioner for frequency-domain seismic wave propagationJack Poulson, Björn Engquist, Siwei Li, and Lexing Ying25 October 2012Efficient 3D frequency response modeling with spectral accuracy by the rapid expansion methodChunlei Chu and Paul L. Stoffa19 June 2012 | GEOPHYSICS, Vol. 77, No. 4Why it is Difficult to Solve Helmholtz Problems with Classical Iterative Methods24 August 2011Fast Algorithms for High Frequency Wave Propagation24 August 2011Analysis of a fast method for solving the high frequency Helmholtz equation in one dimension1 February 2011 | BIT Numerical Mathematics, Vol. 51, No. 3The least squares AMG Solver for the one-dimensional Helmholtz operator26 August 2011 | Computing and Visualization in Science, Vol. 14, No. 1 SEG Technical Program Expanded Abstracts 2011ISSN (print):1052-3812 ISSN (online):1949-4645Copyright: 2011 Pages: 4424 Publisher:Society of Exploration Geophysicists HistoryPublished: 25 May 2012 CITATION INFORMATION Björn Engquist, Jack Poulson, and Lexing Ying, (2011), "Sweeping preconditioner for the 3D Helmholtz equation," SEG Technical Program Expanded Abstracts : 3158-3163. https://doi.org/10.1190/1.3627851 Plain-Language Summary PDF DownloadLoading ...