Sparse Solver Research Articles

A sparse convex AC OPF solver and convex iteration implementation based on 3-node cycles

In this paper, we exploit sparsity technique and strengthen second-order cone programming (SOCP) relaxation of alternating current optimal power flow (AC OPF) through enforcing submatrices corresponding to maximal cliques and minimal chordless cycles in the cycle basis positive semi-definite (PSD). The proposed method adds virtual lines in minimal chordless cycles to decompose each of them into 3-node cycles. By enforcing the submatrices related to 3-node cycles PSD, the resulting convex relaxation has a tight gap. For majority of the test instances, the resulting gap is as tight as that of a semi-definite programming (SDP) relaxation, yet the computing efficiency is much higher. Further, convex iteration is implemented on the proposed solver to achieve exactness by enforcing all submatrices corresponding to lines and virtual lines rank-1. Test results from small size to large size with thousands of buses demonstrate the capability of the proposed solver and the feasibility of the convex iteration implementation.

Parallelization of the solve phase in a task-based Cholesky solver using a sequential task flow model

We describe the parallelization of the solve phase in the sparse Cholesky solver SpLLT when using a sequential task flow model. In the context of direct methods, the solution of a sparse linear system is achieved through three main phases: the analyse, the factorization and the solve phases. In the last two phases, which involve numerical computation, the factorization corresponds to the most computationally costly phase, and it is therefore crucial to parallelize this phase in order to reduce the time-to-solution on modern architectures. As a consequence, the solve phase is often not as optimized as the factorization in state-of-the-art solvers, and opportunities for parallelism are often not exploited in this phase. However, in some applications, the time spent in the solve phase is comparable to or even greater than the time for the factorization, and the user could dramatically benefit from a faster solve routine. This is the case, for example, for a conjugate gradient (CG) solver using a block Jacobi preconditioner. The diagonal blocks are factorized once only, but their factors are used to solve subsystems at each CG iteration. In this study, we design and implement a parallel version of a task-based solve routine for an OpenMP version of the SpLLT solver. We show that we can obtain good scalability on a multicore architecture enabling a dramatic reduction of the overall time-to-solution in some applications.

Elastodiffusion and cluster mobilities using kinetic Monte Carlo simulations: Fast first-passage algorithms for reversible diffusion processes

The microstructural evolution of metals and alloys is governed by the diffusion of defects over complex energy landscapes. Whenever metastability occurs in atomistic simulations, well-separated timescales emerge making it necessary to implement event-based kinetic models at larger scales. The crucial task then involves characterizing the important events contributing to mass transport. We herein describe fast first-passage algorithms based on the theory of absorbing Markov chains assuming that defects undergo reversible diffusion. We show that the absorbing transition rate matrix can be transformed into a symmetric definite-positive matrix enabling us to implement direct and iterative sparse solvers. The efficiency of the approach is demonstrated with direct computations of elastodiffusion properties around a cavity in aluminum and Monte Carlo computations of cluster diffusivity in low-alloyed manganese steels.

A parallel multithreaded sparse triangular linear system solver

We propose a parallel sparse triangular linear system solver based on the Spike algorithm. Sparse triangular systems are required to be solved in many applications. Often, they are a bottleneck due to their inherently sequential nature. Furthermore, typically many successive systems with the same coefficient matrix and with different right hand side vectors are required to be solved. The proposed solver decouples the problem at the cost of extra arithmetic operations as in the banded case. Compared to the banded case, there are extra savings due to the sparsity of the triangular coefficient matrix. We show the parallel performance of the proposed solver against the state-of-the-art parallel sparse triangular solver in Intel’s Math Kernel Library (MKL) on a multicore architecture. We also show the effect of various sparse matrix reordering schemes. Numerical results show that the proposed solver outperforms MKL’s solver in ∼80% of cases by a factor of 2.47, on average.

Efficient use of sparsity by direct solvers applied to 3D controlled-source EM problems

Controlled-source electromagnetic (CSEM) surveying becomes a widespread method for oil and gas exploration, which requires fast and efficient software for inverting large-scale EM datasets. In this context, one often needs to solve sparse systems of linear equations with a large number of sparse right-hand sides, each corresponding to a given transmitter position. Sparse direct solvers are very attractive for these problems, especially when combined with low-rank approximations which significantly reduce the complexity and the cost of the factorization. In the case of thousands of right-hand sides, the time spent in the sparse triangular solve tends to dominate the total simulation time, and here we propose several approaches to reduce it. A significant reduction is demonstrated for marine CSEM application by utilizing the sparsity of the right-hand sides (RHS) and of the solutions that results from the geometry of the problem. Large gains are achieved by restricting computations at the forward substitution stage to exploit the fact that the RHS matrix might have empty rows (vertical sparsity) and/or empty blocks of columns within a non-empty row (horizontal sparsity). We also adapt the parallel algorithms that were designed for the factorization to solve-oriented algorithms and describe performance optimizations particularly relevant for the very large numbers of right-hand sides of the CSEM application. We show that both the operation count and the elapsed time for the solution phase can be significantly reduced. The total time of CSEM simulation can be divided by approximately a factor of 3 on all the matrices from our set (from 3 to 30 million unknowns, and from 4 to 12 thousands RHSs).

Computational kinematics of multibody systems: Two formulations for a modular approach based on natural coordinates

Multibody systems can be divided into an ordered set of kinematically determined modules, known as structural groups, in order to compute their kinematics more efficiently. In this work a procedure for the kinematic analysis of any kind of structural group is introduced, and two different methods for their solution in natural coordinates are presented: the time derivative (TD) and the third-order tensor (3OT) approaches. Moreover, the newly derived methods are compared in terms of efficiency with a global formulation, consisting in solving the kinematics of the multibody system as a whole using dense and sparse solvers. Two scalable case studies have been considered: a 2D four-bar linkage and a 3D slider-crank mechanism with an increasing number of constraint equations. The results show that the TD approach performs better in all cases with speed ups in a range of 27 to 61 times faster in 2D, and of 2.3 to 3.7 times faster in 3D with respect to the global sparse solution.

An evidence based approach to evaluating flood adaptation effectiveness including climate change considerations for coastal cities: City of Port Phillip, Victoria, Australia

AbstractCoastal cities provide a modelling challenge as surface flow is strongly affected by urban drainage networks and there is interaction between coastal and inland flooding. We present a graphics processing unit (GPU)‐based hydrodynamic model coupled to a hydraulic network that integrates adaptation analysis in the context of current and future flooding. The hydrodynamic model is based on a finite volume implementation of the shallow water equations formulated for overland flow. The hydraulic network is based on a pressure relaxation method, and uses a GPU‐based sparse matrix solver for computational speed. The integrated model is used for modelling potential combined coastal and catchment inundation and climate adaptation analysis for the City of Port Phillip, Victoria, Australia. The key outcome from the adaptation study was that resources spent towards adaptation infrastructure should be investigated in the context of sea level rise (SLR) at least for the next 50 years. The adaptation analysis identified “a tipping point” beyond a SLR of around 0.4 m (expected in the next 30 years) where conventional adaptation approaches will fail. This outcome has resulted in the City realising that significant changes in infrastructure for the region will be necessary rather than just incremental adaptation approaches to deal with future flooding.

Parallel accelerated matting method based on local learning

Learning based (LB) matting is an effective matting algorithm but its usability is greatly limited by the heavy computations. In this paper, we cast some new insights into this algorithm and provide following contributions. First, we present an explanation for LB algorithm from manifold learning perspective and unify LB into standard two-stage matting theory. Second, based on the features of two stages, we propose an acceleration scheme utilizing both CPU and GPU parallelism to speed up LB matting up to 15X. Third, we propose an image partition method which provides optimized loading balance and precision for CPU-based block-level parallelism. Finally, we analyze the performance of sparse linear solver used by general matting problems and provide default optimized choice for solver selection. The experiments on the latest parallel framework and mathematical library prove that our scheme performs well evaluated by both acceleration effects and precision.

High-Dimensional Uncertainty Quantification of Electronic and Photonic IC With Non-Gaussian Correlated Process Variations

Uncertainty quantification based on generalized polynomial chaos has been used in many applications. It has also achieved great success in variation-aware design automation. However, almost all existing techniques assume that the parameters are mutually independent or Gaussian correlated, which is rarely true in real applications. For instance, in chip manufacturing, many process variations are actually correlated. Recently, some techniques have been developed to handle non-Gaussian correlated random parameters, but they are time-consuming for high-dimensional problems. We present a new framework to solve uncertainty quantification problems with many non-Gaussian correlated uncertainties. First, we propose a set of smooth basis functions to well capture the impact of non-Gaussian correlated process variations. We develop a tensor approach to compute these basis functions in a high-dimension setting. Second, we investigate the theoretical aspect and practical implementation of a sparse solver to compute the coefficients of all basis functions. We provide some theoretical analysis for the exact recovery condition and error bound of this sparse solver in the context of uncertainty quantification. We present three adaptive sampling approaches to improve the performance of the sparse solver. Finally, we validate our methods by synthetic and practical electronic/photonic ICs with 19 to 57 non-Gaussian correlated variation parameters. Our approach outperforms Monte Carlo by thousands of times in terms of efficiency. It can also accurately predict the output density functions with multiple peaks caused by non-Gaussian correlations, which are hard to capture by existing methods.

Auto-Gauging of Vector Potential by Parallel Sparse Direct Solvers—Numerical Observations

When using magnetic vector potential (MVP)-based formulations for magnetostatic or eddy-current problems, either gauge conditions specifying the divergence of the MVP or tree gauging by eliminating redundant degrees of freedom of the MVP is usually imposed to ensure uniqueness of solutions. Explicit gauging of the MVP is not always necessary since classical iterative solvers can automatically and implicitly fix the gauge as long as the right-hand side vectors are consistent. Besides iterative solvers, implicit gauging is also observed when using state-of-the-art parallel sparse direct solvers (PSDSs), thanks to the built-in functions of handling null-spaces of either real symmetric positive semi-definite matrix systems or those complex symmetric systems from eddy-current problems. Both static and eddy-current examples are solved by PSDS to demonstrate results of local physical quantities or global quantities such as magnetic energy or joule losses. High-order edge/nodal elements are also considered in our numerical examples and it is observed that PSDS can also easily and correctly handle the delicate discrete null spaces.

Improved Equilibrated Error Estimates for Open Boundary Magnetostatic Problems Based on Dual A and H Formulations

Calculating the bounds of global energy is an important issue in computational electromagnetism, which can provide guaranteed results when extracting inductance parameters. In this paper, an improved equilibrated type a posteriori error estimate for open boundary magnetostatic problems is proposed. We derive our error estimator based on vector dual formulations, which can be efficiently solved using parallel sparse direct solvers. The new estimator can provide a sharp and guaranteed estimate of the finite-element spatial discretization error. Moreover, the computational cost is cheaper than using existing equilibrated error estimators. Numerical experiments are carried out to showcase the performance of our error estimator, including the modified TEAM workshop problem 13 and the benchmark IEEJ problem.

An order N numerical method to efficiently calculate the transport properties of large systems: An algorithm optimized for sparse linear solvers

We present a self-contained description of the wave-function matching (WFM) method to calculate electronic quantum transport properties of nanostructures using the Landauer-Büttiker approach. The method is based on a partition of the system between a central region (“conductor”) containing NS sites and an asymptotic region (“leads”) characterized by NP open channels. The two subsystems are linearly coupled and solved simultaneously using an efficient sparse linear solver. Invoking the sparsity of the Hamiltonian matrix representation of the central region, we show that the number of operations required by the WFM method in conductance calculations scales linearly with the number of sites, more precisely with ∼NS×NP for large NS, faster than previously claimed.

Analysis of quasi-Monte Carlo methods for elliptic eigenvalue problems with stochastic coefficients

We consider the forward problem of uncertainty quantification for the generalised Dirichlet eigenvalue problem for a coercive second order partial differential operator with random coefficients, motivated by problems in structural mechanics, photonic crystals and neutron diffusion. The PDE coefficients are assumed to be uniformly bounded random fields, represented as infinite series parametrised by uniformly distributed i.i.d. random variables. The expectation of the fundamental eigenvalue of this problem is computed by (a) truncating the infinite series which define the coefficients; (b) approximating the resulting truncated problem using lowest order conforming finite elements and a sparse matrix eigenvalue solver; and (c) approximating the resulting finite (but high dimensional) integral by a randomly shifted quasi-Monte Carlo lattice rule, with specially chosen generating vector. We prove error estimates for the combined error, which depend on the truncation dimension s, the finite element mesh diameter h, and the number of quasi-Monte Carlo samples N. Under suitable regularity assumptions, our bounds are of the particular form $${\mathcal {O}}(h^2 + N^{-1 + \delta })$$ , where $$\delta > 0$$ is arbitrary and the hidden constant is independent of the truncation dimension, which needs to grow as $$h\rightarrow 0$$ and $$N \rightarrow \infty $$ . As for the analogous PDE source problem, the conditions under which our error bounds hold depend on a parameter $$p \in (0, 1)$$ representing the summability of the terms in the series expansions of the coefficients. Although the eigenvalue problem is nonlinear, which means it is generally considered harder than the source problem, in almost all cases ( $$p \ne 1$$ ) we obtain error bounds that converge at the same rate as the corresponding rate for the source problem. The proof involves a detailed study of the regularity of the fundamental eigenvalue as a function of the random parameters. As a key intermediate result in the analysis, we prove that the spectral gap (between the fundamental and the second eigenvalues) is uniformly positive over all realisations of the random problem.

Benefits from using mixed precision computations in the ELPA-AEO and ESSEX-II eigensolver projects

We first briefly report on the status and recent achievements of the ELPA-AEO (Eigen value Solvers for Petaflop Applications—Algorithmic Extensions and Optimizations) and ESSEX II (Equipping Sparse Solvers for Exascale) projects. In both collaboratory efforts, scientists from the application areas, mathematicians, and computer scientists work together to develop and make available efficient highly parallel methods for the solution of eigenvalue problems. Then we focus on a topic addressed in both projects, the use of mixed precision computations to enhance efficiency. We give a more detailed description of our approaches for benefiting from either lower or higher precision in three selected contexts and of the results thus obtained.

Single stage DOA-frequency representation of the array data with source reconstruction capability

In this paper, a new signal processing framework is proposed, in which the array time samples are represented in DOA-frequency domain through a single stage problem. It is shown that concatenated array data is well represented in a G dictionary atoms space, where columns of G correspond to pixels in the DOA-frequency image. We present two approaches for the G formation and compare the benefits and disadvantages of them. A mutual coherence guaranteed G construction technique is also proposed. Furthermore, unlike most of the existing methods, the proposed problem is reversible into the time domain, therefore, source recovery from the resulted DOA-frequency image is possible. The proposed representation in DOA-frequency domain can be simply transformed into a group sparse problem, in the case of non-multitone sources in a given bandwidth. Therefore, it can also be utilized as an effective wideband DOA estimator. In the simulation part, two scenarios of multitone sources with unknown frequency and DOA locations and non-multitone wideband sources with assumed frequency region are examined. In multitone scenario, sparse solvers yield more accurate DOA-frequency representation compared to some noncoherent approaches. For non-multitone wideband source scenario, the proposed method with group sparse solver outperforms some existing wideband DOA estimators in low SNR regime. In addition, simultaneous sources’ recovery and DOA estimation shows significant improvement compared to the conventional delay and sum beamformer and without prerequisites required in sophisticated wideband beamformers.

Computational Efficiency on Frequency Domain Analysis of Large-scale Finite Element Model by Combination of Iterative and Direct Sparse Solver

대규모 유한요소 모델을 빠르게 해석하기는 위해서 병렬 희소 솔버를 필수적으로 적용해야 한다. 이 논문에서는 미세하게 변화하는 시스템 행렬을 대상으로 연속적으로 해를 구해야 하는 문제에서 효율적으로 적용가능한 반복-직접 희소 솔버 조 합 기법을 소개한다. 반복-직접 희소 솔버 조합 기법은 병렬 희소 솔버 패키지인 PARDISO에 제안 및 구현된 기법으로 새 롭게 행렬값이 갱신된 선형 시스템의 해를 구할 때 이전 선형 시스템에 적용된 직접 희소 솔버의 행렬 분해(factorization) 결과를 Krylov 반복 희소 솔버의 preconditioner로 활용하는 방법을 의미한다. PARDISO에서는 미리 설정된 반복 회수까지 해가 수렴하지 않으면 직접 희소 솔버로 해를 구하며, 이후 이어지는 갱신된 선형 시스템의 해를 구할 때는 최종적으로 사용 된 직법 희소 솔버의 행렬 분해 결과를 preconditioner로 사용한다. 이 연구에서는 첫 번째 Krylov 반복 단계에서 소요되는 시간을 동적으로 계산하여 최대 반복 회수를 설정하는 기법을 제안하였으며, 주파수 영역 해석에 적용하여 그 효과를 검증 하였다.

Accelerating the task/data-parallel version of ILUPACK’s BiCG in multi-CPU/GPU configurations

ILUPACK is a valuable tool for the solution of sparse linear systems via iterative Krylov subspace-based methods. Its relevance for the solution of real problems has motivated several efforts to enhance its performance on parallel machines. In this work we focus on exploiting the task-level parallelism derived from the structure of the BiCG method, in addition to the data-level parallelism of the internal matrix computations, with the goal of boosting the performance of a GPU (graphics processing unit) implementation of this solver. First, we revisit the use of dual-GPU systems to execute independent stages of the BiCG concurrently on both accelerators, while leveraging the extra memory space to improve the data access patterns. In addition, we extend our ideas to compute the BiCG method efficiently in multicore platforms with a single GPU. In this line, we study the possibilities offered by hybrid CPU-GPU computations, as well as a novel synchronization-free sparse triangular linear solver. The experimental results with the new solvers show important acceleration factors with respect to the previous data-parallel CPU and GPU versions.

Improving resilience of scientific software through a domain-specific approach

In this paper we present research on improving the resilience of the execution of scientific software, an increasingly important concern in High Performance Computing (HPC). We build on an existing high-level abstraction framework, the Oxford Parallel library for Structured meshes (OPS), developed for the solution of multi-block structured mesh-based applications, and implement an algorithm in the library to carry out checkpointing automatically, without the intervention of the user. The target applications are a hydrodynamics benchmark application from the Mantevo Suite, CloverLeaf 3D, the sparse linear solver proxy application TeaLeaf, and the OpenSBLI compressible Navier–Stokes direct numerical simulation (DNS) solver. We present (1) the basic algorithm that OPS relies on to determine the optimal checkpoint in terms of size and location, (2) improvements that supply additional information to improve the decision, (3) techniques that reduce the cost of writing the checkpoints to non-volatile storage, (4) a performance analysis of the developed techniques on a single workstation and on several supercomputers, including ORNL’s Titan. Our results demonstrate the utility of the high-level abstractions approach in automating the checkpointing process and show that performance is comparable to, or better than the reference in all cases.

Moving mesh finite difference solution of non-equilibrium radiation diffusion equations

A moving mesh finite difference method based on the moving mesh partial differential equation is proposed for the numerical solution of the 2T model for multi-material, non-equilibrium radiation diffusion equations. The model involves nonlinear diffusion coefficients and its solutions stay positive for all time when they are positive initially. Nonlinear diffusion and preservation of solution positivity pose challenges in the numerical solution of the model. A coefficient-freezing predictor-corrector method is used for nonlinear diffusion while a cutoff strategy with a positive threshold is used to keep the solutions positive. Furthermore, a two-level moving mesh strategy and a sparse matrix solver are used to improve the efficiency of the computation. Numerical results for a selection of examples of multi-material non-equilibrium radiation diffusion show that the method is capable of capturing the profiles and local structures of Marshak waves with adequate mesh concentration. The obtained numerical solutions are in good agreement with those in the existing literature. Comparison studies are also made between uniform and adaptive moving meshes and between one-level and two-level moving meshes.

Frame Theory and Fractional Programming for Sparse Recovery-Based mmWave Channel Estimation

We consider the estimation of millimeter wave (mmWave) channel gains following the now well-established approach of sparse signal processing. Within that context, we offer two major contributions. The first is a complete frame-theoretical treatment of the sensing matrix design (beam management) problem, compactly described by a pair of Lemmas that together with efficient low-coherence frame construction algorithms offer a general solution for the optimal transmitter (Tx) and receiver (Rx) beamforming components of the sparse mmWave channel estimation problem. The second contribution is a pair of novel sparse recovery algorithms, which unlike the majority of sparse solvers found in the literature, is not based on the relaxation of l 0 -norm into an l 1 -norm, but rather on a smooth approximation of the l 0 -norm handled further via fractional programming (FPG). The two algorithms differ from each other in that in the first (slightly more accurate), the resulting convex problem is solved via interior point methods, while the second (stand-alone) makes use of the alternating direction method of multipliers (ADMM). As a bonus, an original ADMM variation of the well-known basis pursuit denoising (BPDN)-l 1 -reweighted sparse solver is also given. Simulation results confirm the channel estimation accuracy improvements obtained by both contributions.