A new parallel sparse direct solver: Presentation and numerical experiments in large‐scale structural mechanics parallel computing

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

International audience

Cooling process on a run-out table by the simulation method

SummarySequentially linear analysis (SLA), an event‐by‐event procedure for finite element (FE) simulation of quasi‐brittle materials, is based on sequentially identifying a critical integration point in the FE model, to reduce its strength and stiffness, and the corresponding critical load multiplier (λcrit), to scale the linear analysis results. In this article, two strategies are proposed to efficiently reuse previous stiffness matrix factorisations and their corresponding solutions in subsequent linear analyses, since the global system of linear equations representing the FE model changes only locally. The first is based on a direct solution method in combination with the Woodbury matrix identity, to compute the inverse of a low‐rank corrected stiffness matrix relatively cheaply. The second is a variation of the traditional incomplete LU preconditioned conjugate gradient method, wherein the preconditioner is the complete factorisation of a previous analysis step's stiffness matrix. For both the approaches, optimal points at which the factorisation is recomputed are determined such that the total analysis time is minimised. Comparison and validation against a traditional parallel direct sparse solver, with regard to a two‐dimensional (2D) and three‐dimensional (3D) benchmark study, illustrates the improved performance of the Woodbury‐based direct solver over its counterparts, especially for large 3D problems.

Accurate simulation of seismic waves is essential for achieving high‐precision full‐waveform inversion (FWI). Within the Cartesian coordinate system‐based frequency‐domain finite‐difference (FDFD) framework, we propose a one‐direction composition average‐derivative optimal method for the 3D heterogeneous isotropic elastic‐wave equation, referred to as the 45‐point scheme. The results of dispersion analysis and weighted coefficient optimization demonstrate that the 45‐point scheme achieves higher dispersion accuracy than the existing 27‐point average‐derivative scheme. More importantly, by constructing the impedance matrix along the ‘composition’ direction, the bandwidth of the sparse impedance matrix increases only slightly, with nonzero elements compactly distributed in strips. On the basis of the multifrontal massively parallel sparse direct solver (MUMPS) on a supercomputer platform, the 45‐point scheme does not significantly increase computational complexity compared to the 27‐point scheme. To further test the performance of the 45‐point scheme, we provide several numerical experiments, including simple homogeneous and complex SEG/EAGE overthrust models. In comparison with the 27‐point scheme, the 45‐point scheme yields a notable improvement in computational accuracy, particularly for large grid ratios, while imposing only a modest increase in computational cost. These findings thus strongly suggest that the 45‐point scheme holds promise as a viable option for the forward part of frequency‐domain FWI in practical high‐accuracy seismic imaging applications.

Summary form only given. Fast direct integral equations methods for analyzing time-harmonic scattering phenomena constitute an active area of research in the CEM community. These solvers present an interesting alternative to FMM-based iterative schemes for problems that are inherently ill-conditioned and/or involve many excitations. Present fast direct methods leverage the low-rank (LR) nature of blocks of the LU factors or inverse of the method of moment (MoM) matrix (if and when this property manifests itself). In essence, they leverage the limited number of degrees of freedom of the fields radiated by sources in the presence of the scatterer (or portions thereof), observed across a finite region of the scatterer (whenever the concept makes sense). LR compression strategies do not lead to low-complexity solvers for electrically large concave structures. This fact was demonstrated in [H. Guo et al., Antennas and Wireless Propagation Letters, IEEE, vol. 11, 2012], where blocks of LU factors that proved incompressible by LR methods were recompressed by MLMDA/butterfly methods. Somewhat unexpectedly, this study showed that MLMDA/butterfly scheme produced compression above and beyond that provided by LR schemes, requiring only O(N log2 N) final storage irrespective of the nature or electrical size of the scatterer. The practical applicability of the scheme however was severely hindered by the fact that the MLMDA/butterfly scheme acted as an add-on to a LR solver. Here, we report on a new direct integral equation solver by that utilizes MLMDA/butterfly schemes for compressing the LU factors of a MoM matrix. When compared to its predecessor, the new solver has three important features of note. (i) The solver entirely bypasses the LR compression step of its predecessor and operates directly on butterfly-compressed blocks. The latter is achieved using new schemes for rapidly adding and multiplying butterfly-compressed operators. (ii) The solver executes in parallel, using a hybrid OpenMP-MPI strategy to accelerate the hierarchical inverse/decomposition of the MoM matrix. Despite the inherently sequential nature of this operation, the solver exhibits excellent scaling properties up to several hundred processors. (iii) The solver applies to large-scale 3D (as opposed to 2D) analysis and was implemented to invert a combined field (as opposed to an electric field) integral operator. The solver is capable of inverting MoM matrices that discretize combined field integral equations modeling scattering from electrically large structures involving millions of unknowns on a small cluster using just a few hours of CPU time.

Sparse matrix is used in many modern industrial applications, to provide low cost and high performance solution. Unlike a number of previous works in the field of distribution networks, which analyzed the sparse matrix calculations as integral part of its calculation engine, this paper presents the considerations and analysis for the utilization of commercial of-the-shelf sparse matrix solvers (COTS SMS). Analysis is given for the parallel direct sparse solver utilized in unbalanced and unsymmetrical Distribution System Power Flow solution, which maximizes parallelization on multi-core processors. Optimal matrix population algorithm and container class from the power network object model is proposed for the fast population, search and access to the matrix elements. Optimal use cases of symmetrical and unsymmetrical matrices for the given method are analyzed. Test results for the Current Injection—Power Flow algorithm on a Three-Phase large scale system are presented which show the performance scales for symmetric and unsymmetrical matrices over a scope of system sizes and for different network topologies.

Two state‐of‐the‐art parallel software packages for the direct solution of sparse linear systems based on LU‐decomposition, MUMPS and SuperLU_DIST have been tested as black‐box solvers on problems derived from finite difference discretizations of the Helmholtz equation. The target architecture has been Linux clusters, for which no consistent set of tests of the algorithms implemented in these packages has been published. The investigation consists of series of memory and time scalability checks and has focused on examining the applicability of the algorithms when processing very large sparse matrices on Linux cluster platforms. Special emphasis has been put on monitoring the behaviour of the packages when the equation systems need to be solved for multiple right‐hand sides, which is the case, for instance, when modelling a seismic survey. The outcome of the tests points at poor efficiency of the tested algorithms during application of the LU‐factors in the solution phase on this type of architecture, where the communication acts as an impasse. Copyright © 2005 John Wiley & Sons, Ltd.

Research is on-going that examines parallel direct block-diagonal-bordered sparse linear solvers for irregular sparse matrix problems derived from electrical power system applications. Parallel block-diagonal-bordered sparse linear solvers exhibit distinct advantages when compared to current general parallel direct sparse matrix solvers. Our research shows that actual power system matrices can be readily ordered into block-diagonal-bordered form, although load imbalance becomes excessive beyond 16 processors, limiting scalability for a single parallel linear solver within an application. Nevertheless, other dimensions exist in electrical power system applications that can be exploited to efficiently make use of large-scale multi-processors. >

The effectiveness of sparse matrix techniques for directly solving large-scale linear least-squares problems is severely limited if the system matrix A has one or more nearly dense rows. In this paper, we partition the rows of A into sparse rows and dense rows (As and Ad) and apply the Schur complement approach. A potential difficulty is that the reduced normal matrix AsTAs is often rank-deficient, even if A is of full rank. To overcome this, we propose explicitly removing null columns of As and then employing a regularization parameter and using the resulting Cholesky factors as a preconditioner for an iterative solver applied to the symmetric indefinite reduced augmented system. We consider complete factorizations as well as incomplete Cholesky factorizations of the shifted reduced normal matrix. Numerical experiments are performed on a range of large least-squares problems arising from practical applications. These demonstrate the effectiveness of the proposed approach when combined with either a sparse parallel direct solver or a robust incomplete Cholesky factorization algorithm.

The goal of this paper is to develop a parallel algorithm for the direct solution of large sparse linear systems and integrate it into domain decomposition methods. The computational effort for these linear systems, often encountered in numerical simulation of structural mechanics problems by finite element codes, is very significant in terms of run-time and memory requirements.In this paper, a two-level parallelism is exploited. The exploitation of the lower level of parallelism is based on the development of a parallel direct solver with a nested dissection algorithm and to introduce it into the FETI methods. This direct solver has the advantage of handling zero-energy modes in floating structures automatically and properly. The upper level of parallelism is a coarse-grain parallelism between substructures of FETI. Some numerical tests are carried out to evaluate the performance of the direct solver.

This paper presents a parallelized version of the edge-based finite element method with a novel post-processing approach for numerical modeling of an electromagnetic field in complex media. The method uses an unstructured tetrahedral mesh which can reduce the number of degrees of freedom significantly. The linear system of finite element equations is solved using parallel direct solvers which are robust for ill-conditioned systems and efficient for multiple source electromagnetic (EM) modeling. We also introduce a novel approach to compute the scalar components of the electric field from the tangential components along each edge based on field redatuming. The method can produce a more accurate result as compared to conventional approach. We have applied the developed algorithm to compute the EM response for a typical 3D anisotropic geoelectrical model of the off-shore HC reservoir with complex seafloor bathymetry. The numerical study demonstrates that the modeling algorithm is capable of simulating the complex topography and bathymetry that is commonly encountered in controlled source electromagnetic problems.

This report concerns the generation of support vector machine classifiers for solving the pattern recognition problem in machine learning. Several methods are proposed based on interior point methods for convex quadratic programming. Software implementations are developed by adapting the object-oriented packaging OOQP to the problem structure and by using the software package PETSc to perform time-intensive computations in a distributed setting. Linear systems arising from classification problems with moderately large numbers of features are solved by using two techniques--one a parallel direct solver, the other a Krylov-subspace method incorporating novel preconditioning strategies. Numerical results are provided, and computational experience is discussed.

Robust and scalable adaptive BDDC preconditioners for virtual element discretizations of elliptic partial differential equations in mixed form

Simulators of semiconductor devices have to solve systems of equations generated by the discretization of partial differential equations, which are the most time-consuming part of the simulation process. Therefore, the use of an effective method to solve these linear systems is essential. In this work we have evaluated the efficiency of different parallel direct and iterative solvers used for the solution of the drift–diffusion equations in semiconductor device simulation. Several preconditioning techniques have been applied in order to minimize the execution times. We have found that FGMRES and BCGSTAB solvers preconditioned with Additive Schwarz are the most suitable for these types of problems. The results were obtained in an HP Superdome cluster with 128 Itanium2 1.5 GHz. Copyright © 2005 John Wiley & Sons, Ltd.