Sparse Direct Solver Research Articles

Optimal numerical solvers for transient simulations of ice flow using the Ice Sheet System Model (ISSM versions 4.2.5 and 4.11)

Abstract. Identifying fast and robust numerical solvers is a critical issue that needs to be addressed in order to improve projections of polar ice sheets evolving in a changing climate. This work evaluates the impact of using advanced numerical solvers for transient ice-flow simulations conducted with the JPL–UCI Ice Sheet System Model (ISSM). We identify optimal numerical solvers by testing a broad suite of readily available solvers, ranging from direct sparse solvers to preconditioned iterative methods, on the commonly used Ice Sheet Model Intercomparison Project for Higher-Order ice sheet Models benchmark tests. Three types of analyses are considered: mass transport, horizontal stress balance, and incompressibility. The results of the fastest solvers for each analysis type are ranked based on their scalability across mesh size and basal boundary conditions. We find that the fastest iterative solvers are ∼ 1.5–100 times faster than the default direct solver used in ISSM, with speed-ups improving rapidly with increased mesh resolution. We provide a set of recommendations for users in search of efficient solvers to use for transient ice-flow simulations, enabling higher-resolution meshes and faster turnaround time. The end result will be improved transient simulations for short-term, highly resolved forward projections (10–100 year time scale) and also improved long-term paleo-reconstructions using higher-order representations of stresses in the ice. This analysis will also enable a new generation of comprehensive uncertainty quantification assessments of forward sea-level rise projections, which rely heavily on ensemble or sampling approaches that are inherently expensive.

Reordering Strategy for Blocking Optimization in Sparse Linear Solvers

Solving sparse linear systems is a problem that arises in many scientific applications, and sparse direct solvers are a time-consuming and key kernel for those applications and for more advanced solvers such as hybrid direct-iterative solvers. For this reason, optimizing their performance on modern architectures is critical. The preprocessing steps of sparse direct solvers---ordering and block-symbolic factorization---are two major steps that lead to a reduced amount of computation and memory and to a better task granularity to reach a good level of performance when using BLAS kernels. With the advent of GPUs, the granularity of the block computation has become more important than ever. In this paper, we present a reordering strategy that increases this block granularity. This strategy relies on block-symbolic factorization to refine the ordering produced by tools such as Metis or Scotch, but it does not impact the number of operations required to solve the problem. We integrate this algorithm in the PaSt...

Solving Mixed Sparse-Dense Linear Least-Squares Problems by Preconditioned Iterative Methods

The efficient solution of large linear least-squares problems in which the system matrix $A$ contains rows with very different densities is challenging. Previous work has focused on direct methods for problems in which $A$ has a few relatively dense rows. These rows are initially ignored, a factorization of the sparse part is computed using a sparse direct solver, and then the solution is updated to take account of the omitted dense rows. In some practical applications the number of dense rows can be significant, and for very large problems, using a direct solver may not be feasible. We propose processing rows that are identified as dense separately within a conjugate gradient method using an incomplete factorization preconditioner combined with the factorization of a dense matrix of size equal to the number of dense rows. Numerical experiments on large-scale problems from real applications are used to illustrate the effectiveness of our approach. The results demonstrate that we can efficiently solve prob...

Analysis of a Splitting Approach for the Parallel Solution of Linear Systems on GPU Cards

We discuss an approach for solving sparse or dense banded linear systems $Ax = b$ on a graphics processing unit (GPU) card. The matrix ${A} \in {\mathbb{R}}^{N \times N}$ is possibly nonsymmetric and moderately large, i.e., ${10,000}{} \leq N \leq {500,000}{}$. The split and parallelize (SaP) approach seeks to partition the matrix $A$ into diagonal subblocks $A_i$, $i=1,\ldots,P$, which are independently factored in parallel. The solution may choose to consider or to ignore the matrices that couple the diagonal subblocks $A_i$. This approach, along with the Krylov-subspace-based iterative method that it preconditions, are implemented in a solver called SaP::GPU, which is compared in terms of efficiency with three commonly used sparse direct solvers: PARDISO, SuperLU, and MUMPS. SaP::GPU, which runs entirely on the GPU except for several stages involved in preliminary row and column permutations, is robust and compares well in terms of efficiency with the aforementioned direct solvers. In a comparison against Intel's MKL, SaP::GPU also fared well when used to solve dense banded systems that are close to being diagonally dominant. SaP::GPU is publicly available and distributed as open source under a permissive BSD-3 license.

On Using Cholesky-Based Factorizations and Regularization for Solving Rank-Deficient Sparse Linear Least-Squares Problems

By examining the performance of modern parallel sparse direct solvers and exploiting our knowledge of the algorithms behind them, we perform numerical experiments to study how they can be used to efficiently solve rank-deficient sparse linear least-squares problems arising from practical applications. The Cholesky factorization of the normal equations breaks down when the least-squares problem is rank-deficient, while applying a symmetric indefinite solver to the augmented system can give an unacceptable level of fill in the factors. To try to resolve these difficulties, we consider a regularization procedure that modifies the diagonal of the unregularized matrix. This leads to matrices that are easier to factorize. We consider both the regularized normal equations and the regularized augmented system. We employ the computed factors of the regularized systems as preconditioners with an iterative solver to obtain the solution of the original (unregularized) problem. Furthermore, we look at using limited-me...

Three-Dimensional Modeling of the Casing Effect in Onshore Controlled-Source Electromagnetic Surveys

The presence of steel-cased wells and other infrastructure causes a significant change in the electromagnetic fields that has to be taken into consideration in modeling and interpretation of field data. A realistic and accurate simulation requires the borehole casing to be incorporated into the modeling scheme, which is numerically challenging. Due to the huge conductivity contrast between the casing and surrounding media, a spatial discretization that provides accurate results at different spatial scales ranging from millimeters to hundreds of meters is required. In this paper, we present a full 3D frequency-domain electromagnetic modeling based on a parallel finite-difference algorithm considering the casing effect and investigate its applicability on the borehole-to-surface configuration of the Hontomin CO2 storage site. To guarantee a robust solution of linear systems with highly ill-conditioned matrices caused by huge conductivity contrasts and multiple spatial scales in the model, we employ direct sparse solvers. Different scenarios are simulated in order to study the influence of the source position, conductivity model, and the effect of the steel casing on the measured data. Several approximations of the real hollow casing that allow for a large reduction in the number of elements in the resulting meshes are studied. A good agreement between the modeled responses and the real field data demonstrates the feasibility of simulating casing effects in complex geological areas. The steel casing of the well greatly increases the amplitude of the surface electromagnetic fields and thus improves the signal-to-noise ratio and the sensitivity to deep targets.

On the improvement of a scalable sparse direct solver for unsymmetrical linear equations

This paper focuses on the application level improvements in a sparse direct solver specifically used for large-scale unsymmetrical linear equations resulting from unstructured mesh discretization of coupled elliptic/hyperbolic PDEs. Existing sparse direct solvers are designed for distributed server systems taking advantage of both distributed memory and processing units. We conducted extensive numerical experiments with three state-of-the-art direct linear solvers that can work on distributed-memory parallel architectures; namely, MUMPS (MUMPS solver website, http://graal.ens-lyon.fr/MUMPS), WSMP (Technical Report TR RC-21886, IBM, Watson Research Center, Yorktown Heights, 2000), and SUPERLU_DIST (ACM Trans Math Softw 29(2):110---140, 2003). The performance of these solvers was analyzed in detail, using advanced analysis tools such as Tuning and Analysis Utilities (TAU) and Performance Application Programming Interface (PAPI). The performance is evaluated with respect to robustness, speed, scalability, and efficiency in CPU and memory usage. We have determined application level issues that we believe they can improve the performance of a distributed-shared memory hybrid variant of this solver, which is proposed as an alternative solver [SuperLU_MCDT (Many-Core Distributed)] in this paper. The new solver utilizing the MPI/OpenMP hybrid programming is specifically tuned to handle large unsymmetrical systems arising in reservoir simulations so that higher performance and better scalability can be achieved for a large distributed computing system with many nodes of multicore processors. Two main tasks are accomplished during this study: (i) comparisons of public domain solver algorithms; existing state-of-the-art direct sparse linear system solvers are investigated and their performance and weaknesses based on test cases are analyzed, (ii) improvement of direct sparse solver algorithm (SuperLU_MCDT) for many-core distributed systems is achieved. We provided results of numerical tests that were run on up to 16,384 cores, and used many sets of test matrices for reservoir simulations with unstructured meshes. The numerical results showed that SuperLU_MCDT can outperform SuperLU_DIST 3.3 in terms of both speed and robustness.

On some practical issues concerning the implementation of Cahn–Hilliard–Navier–Stokes type models

Suspensions, solutions and colloids are physical systems that can be generally denoted as “multicomponent systems” or “mixtures”. Physical systems of this type are frequently met in many industrial applications which rises the need for numerical simulations of the behaviour of such complex systems. A particular practical example of such system is glass/tin/nitrogen system that must be studied in modelling of the float glass process (Pilkington process) that partly motivated this research. The aim of this paper is to discuss the numerical challenges along with some practical issues concerning the implementation of a chosen Cahn–Hilliard–Navier–Stokes type model that showed up to be particularly suitable for the description of three-component systems. The discretization of the system of governing partial differential equations is based on the finite element method using the FEniCS Project. Numerical experiments carried out in two and three spatial dimensions verify our straightforward implementation that uses parallel direct sparse solvers to resolve the intermediate linear systems of algebraic equations. Possible improvements of the current implementation are briefly outlined within the paper.

PSPIKE+: A family of parallel hybrid sparse linear system solvers

We present PSPIKE+, as a family of parallel hybrid linear system solvers that are more robust than other available preconditioned iterative methods, and more scalable than parallel sparse direct solvers. A critical step in using PSPIKE+ is the extraction of a preconditioner that is in the form of a generalized central band which consists of a set of overlapping diagonal blocks that encapsulate as many of the heaviest elements of the coefficient matrix. To address this issue, we propose the G-PAVER reordering scheme that produces an effective preconditioner in the form of overlapping diagonal blocks while simultaneously accommodating load balancing on distributed-memory architectures and controlling overlap sizes. In each outer Krylov subspace iteration for solving the linear system, we need to solve linear systems involving the preconditioner. This is accomplished using a parallel tearing-based method, which may be regarded as an algebraic domain decomposition scheme. The tearing-based method requires the solution of a much smaller linear system (we call the balance system) whose size is equal to the sum of the overlaps. In this paper, we form the balance system explicitly and solve it directly. Our parallel implementation of PSPIKE+ based on G-PAVER, PARDISO 5.0.0, and solving the explicit balance system results in impressive robustness and good parallel scalability on computing platforms with a limited number of multi-core nodes.

Fast 3D frequency-domain full-waveform inversion with a parallel block low-rank multifrontal direct solver: Application to OBC data from the North Sea

Wide-azimuth long-offset ocean bottom cable (OBC)/ocean bottom node surveys provide a suitable framework to perform computationally efficient frequency-domain full-waveform inversion (FWI) with a few discrete frequencies. Frequency-domain seismic modeling is performed efficiently with moderate computational resources for a large number of sources with a sparse multifrontal direct solver (Gauss-elimination techniques for sparse matrices). Approximate solutions of the time-harmonic wave equation are computed using a block low-rank (BLR) approximation, leading to a significant reduction in the operation count and in the volume of communication during the lower upper (LU) factorization as well as offering great potential for reduction in the memory demand. Moreover, the sparsity of the seismic source vectors is exploited to speed up the forward elimination step during the computation of the monochromatic wavefields. The relevance and the computational efficiency of the frequency-domain FWI performed in the viscoacoustic vertical transverse isotropic (VTI) approximation was tested with a real 3D OBC case study from the North Sea. The FWI subsurface models indicate a dramatic resolution improvement relative to the initial model built by reflection traveltime tomography. The amplitude errors introduced in the modeled wavefields by the BLR approximation for different low-rank thresholds have a negligible footprint in the FWI results. With respect to a standard multifrontal sparse direct factorization, and without compromise of the accuracy of the imaging, the BLR approximation can bring a reduction of the LU factor size by a factor of up to three. This reduction is not yet exploited to reduce the effective memory usage (ongoing work). The flop reduction can be larger than a factor of 10 and can bring a factor of time reduction of around three. Moreover, this reduction factor tends to increase with frequency, namely with the matrix size. Frequency-domain viscoacoustic VTI FWI can be viewed as an efficient tool to build an initial model for elastic FWI of 4C OBC data.

Accelerating sparse Cholesky factorization on GPUs

Sparse direct factorization is a critical component of scientific computingGPU acceleration of it is difficult due to algorithmic irregularity and PCIe perf.We present an algorithm for GPU acceleration which resolves these issues.The algorithm shows 2x speedup vs. best alternative CPU-only perf.Aspects of this Cholesky case should be useful in other algorithms (LDL, LU, ... ). Sparse factorization is a fundamental tool in scientific computing. As the major component of a sparse direct solver, it represents the dominant computational cost for many analyses. For factorizations which involve sufficient dense math, the substantial computational capability provided by GPUs (Graphics Processing Units) can help alleviate this cost. However, for many other cases, the prevalence of small/irregular dense math and the relatively slow communication between the host and device over the PCIe bus, make it challenging to significantly accelerate sparse factorization using the GPU.In this paper we describe a left-looking supernodal Cholesky factorization algorithm which permits improved utilization of the GPU when factoring sparse matrices. The central idea is to stream subtrees of the elimination tree through the GPU and perform the factorization of each subtree entirely on the GPU. This avoids the majority of the PCIe communication without the need for a complex task scheduler. Importantly, within these subtrees, many independent, small, dense operations are batched to minimize kernel launch overhead and many of these batched kernels are executed concurrently to maximize device utilization.Performance results for commonly studied matrices are presented along with suggested actions for further optimization.

Laplace-domain waveform modeling and inversion for the 3D acoustic–elastic coupled media

Laplace-domain waveform inversion reconstructs long-wavelength subsurface models by using the zero-frequency component of damped seismic signals. Despite the computational advantages of Laplace-domain waveform inversion over conventional frequency-domain waveform inversion, an acoustic assumption and an iterative matrix solver have been used to invert 3D marine datasets to mitigate the intensive computing cost. In this study, we develop a Laplace-domain waveform modeling and inversion algorithm for 3D acoustic–elastic coupled media by using a parallel sparse direct solver library (MUltifrontal Massively Parallel Solver, MUMPS). We precisely simulate a real marine environment by coupling the 3D acoustic and elastic wave equations with the proper boundary condition at the fluid-solid interface. In addition, we can extract the elastic properties of the Earth below the sea bottom from the recorded acoustic pressure datasets. As a matrix solver, the parallel sparse direct solver is used to factorize the non-symmetric impedance matrix in a distributed memory architecture and rapidly solve the wave field for a number of shots by using the lower and upper matrix factors. Using both synthetic datasets and real datasets obtained by a 3D wide azimuth survey, the long-wavelength component of the P-wave and S-wave velocity models is reconstructed and the proposed modeling and inversion algorithm are verified. A cluster of 80CPU cores is used for this study.

3D CSEM data inversion using Newton and Halley class methods

For the first time in 3D controlled source electromagnetic data inversion, we explore the use of the Newton and the Halley optimization methods, which may show their potential when the cost function has a complex topology. The inversion is formulated as a constrained nonlinear least-squares problem which is solved by iterative optimization. These methods require the derivatives up to second order of the residuals with respect to model parameters. We show how Green's functions determine the high-order derivatives, and develop a diagrammatical representation of the residual derivatives. The Green's functions are efficiently calculated on-the-fly, making use of a finite-difference frequency-domain forward modelling code based on a multi-frontal sparse direct solver. This allow us to build the second-order derivatives of the residuals keeping the memory cost in the same order as in a Gauss–Newton (GN) scheme. Model updates are computed with a trust-region based conjugate-gradient solver which does not require the computation of a stabilizer. We present inversion results for a synthetic survey and compare the GN, Newton, and super-Halley optimization schemes, and consider two different approaches to set the initial trust-region radius. Our analysis shows that the Newton and super-Halley schemes, using the same regularization configuration, add significant information to the inversion so that the convergence is reached by different paths. In our simple resistivity model examples, the convergence speed of the Newton and the super-Halley schemes are either similar or slightly superior with respect to the convergence speed of the GN scheme, close to the minimum of the cost function. Due to the current noise levels and other measurement inaccuracies in geophysical investigations, this advantageous behaviour is at present of low consequence, but may, with the further improvement of geophysical data acquisition, be an argument for more accurate higher-order methods like those applied in this paper.

2.5D forward modeling and inversion of frequency-domain airborne electromagnetic data

Frequency-domain airborne electromagnetics is a proven geophysical exploration method. Presently, the interpretation is mainly based on resistivity—depth imaging and one-dimensional layered inversion; nevertheless, it is difficult to obtain satisfactory results for two- or three-dimensional complex earth structures using 1D methods. 3D forward modeling and inversion can be used but are hampered by computational limitations because of the large number of data. Thus, we developed a 2.5D frequency-domain airborne electromagnetic forward modeling and inversion algorithm. To eliminate the source singularities in the numerical simulations, we split the fields into primary and secondary fields. The primary fields are calculated using homogeneous or layered models with analytical solutions, and the secondary (scattered) fields are solved by the finite-element method. The linear system of equations is solved by using the large-scale sparse matrix parallel direct solver, which greatly improves the computational efficiency. The inversion algorithm was based on damping least-squares and singular value decomposition and combined the pseudo forward modeling and reciprocity principle to compute the Jacobian matrix. Synthetic and field data were used to test the effectiveness of the proposed method.

A Sparse Symmetric Indefinite Direct Solver for GPU Architectures

In recent years, there has been considerable interest in the potential for graphics processing units (GPUs) to speed up the performance of sparse direct linear solvers. Efforts have focused on symmetric positive-definite systems for which no pivoting is required, while little progress has been reported for the much harder indefinite case. We address this challenge by designing and developing a sparse symmetric indefinite solver SSIDS. This new library-quality LDL T factorization is designed for use on GPU architectures and incorporates threshold partial pivoting within a multifrontal approach. Both the factorize and the solve phases are performed using the GPU. Another important feature is that the solver produces bit-compatible results. Numerical results for indefinite problems arising from a range of practical applications demonstrate that, for large problems, SSIDS achieves performance improvements of up to a factor of 4.6 × compared with a state-of-the-art multifrontal solver on a multicore CPU.

Evaluation of parallel direct sparse linear solvers in electromagnetic geophysical problems

High performance computing is absolutely necessary for large-scale geophysical simulations. In order to obtain a realistic image of a geologically complex area, industrial surveys collect vast amounts of data making the computational cost extremely high for the subsequent simulations. A major computational bottleneck of modeling and inversion algorithms is solving the large sparse systems of linear ill-conditioned equations in complex domains with multiple right hand sides. Recently, parallel direct solvers have been successfully applied to multi-source seismic and electromagnetic problems. These methods are robust and exhibit good performance, but often require large amounts of memory and have limited scalability. In this paper, we evaluate modern direct solvers on large-scale modeling examples that previously were considered unachievable with these methods. Performance and scalability tests utilizing up to 65,536 cores on the Blue Waters supercomputer clearly illustrate the robustness, efficiency and competitiveness of direct solvers compared to iterative techniques. Wide use of direct methods utilizing modern parallel architectures will allow modeling tools to accurately support multi-source surveys and 3D data acquisition geometries, thus promoting a more efficient use of the electromagnetic methods in geophysics.

An Efficient Multicore Implementation of a Novel HSS-Structured Multifrontal Solver Using Randomized Sampling

We present a sparse linear system solver that is based on a multifrontal variant of Gaussian elimination and exploits low-rank approximation of the resulting dense frontal matrices. We use hierarchically semiseparable (HSS) matrices, which have low-rank off-diagonal blocks, to approximate the frontal matrices. For HSS matrix construction, a randomized sampling algorithm is used together with interpolative decompositions. The combination of the randomized compression with a fast ULV HSS factorization leads to a solver with lower computational complexity than the standard multifrontal method for many applications, resulting in speedups up to sevenfold for problems in our test suite. The implementation targets many-core systems by using task parallelism with dynamic runtime scheduling. Numerical experiments show performance improvements over state-of-the-art sparse direct solvers. The implementation achieves high performance and good scalability on a range of modern shared memory parallel systems, including ...

Frequency-domain elastic full-waveform multiscale inversion method based on dual-level parallelism

The complexity of an elastic wavefield increases the nonlinearity of inversion. To some extent, multiscale inversion decreases the nonlinearity of inversion and prevents it from falling into local extremes. A multiscale strategy based on the simultaneous use of frequency groups and layer stripping method based on damped wave field improves the stability of inversion. A dual-level parallel algorithm is then used to decrease the computational cost and improve practicability. The seismic wave modeling of a single frequency and inversion in a frequency group are computed in parallel by multiple nodes based on multifrontal massively parallel sparse direct solver and MPI. Numerical tests using an overthrust model show that the proposed inversion algorithm can effectively improve the stability and accuracy of inversion by selecting the appropriate inversion frequency and damping factor in lowfrequency seismic data.

Parallel direct solver for solving systems of linear equations resulting from finite element method on multi-core desktops and workstations

This paper presents the PARFES—parallel sparse direct solver for solving systems of linear equations resulting from the finite element method applied to problems of structural and solid mechanics. PARFES is intended for use on shared-memory multi-core desktops, laptops and workstations. The solver’s high performance relies on the use of a specific algorithm for the procedure of block matrix multiplication using AVX instructions for the vectorization of computations, and the parallelizing approach that ensures a balanced processor load. The block looking-left factorization algorithm allows for the use of disk memory if the dimensionality of problem exceeds the capacity of RAM.

Direct Finite-Element Solver of Linear Complexity for Large-Scale 3-D Electromagnetic Analysis and Circuit Extraction

In this paper, we develop a linear-complexity direct finite-element solver for the electromagnetic analysis of general 3-D problems containing arbitrarily shaped lossy or lossless conductors in inhomogeneous materials. Both theoretical analysis and numerical experiments have demonstrated the solver’s linear complexity in CPU time and memory consumption with prescribed accuracy satisfied. The proposed direct solver has successfully analyzed an industry product-level full package involving over 22.8488 million unknowns in approximately 16 h on a single core running at 3 GHz. It has also rapidly solved large-scale antenna arrays of over 73 wavelengths with 3600 antenna elements whose number of unknowns is over 10 million. The proposed direct solver has been compared with the finite-element methods that utilize the most advanced direct sparse solvers and a widely used commercial iterative finite-element solver. Clear advantages of the proposed solver in time and memory complexity, as well as computational efficiency, have been demonstrated.