Sparse Linear Systems Research Articles

Fast parallel-in-time quasi-boundary value methods for backward heat conduction problems

In this paper we propose two new quasi-boundary value methods for regularizing the ill-posed backward heat conduction problems. With a standard finite difference discretization in space and time, the obtained all-at-once nonsymmetric sparse linear systems have the desired block ω-circulant structure, which can be utilized to design an efficient parallel-in-time (PinT) direct solver that is built upon an explicit FFT-based diagonalization of the time discretization matrix. Convergence analysis is presented to justify the optimal choice of the regularization parameter under suitable assumptions. Numerical examples are reported to validate our analysis and illustrate the superior computational efficiency of our proposed PinT methods.

Balanced incomplete factorization preconditioner with pivoting

In this work we study pivoting strategies for the preconditioner presented in Bru (SIAM J Sci Comput 30(5):2302–2318, 2008) which computes the LU factorization of a matrix A. This preconditioner is based on the Inverse Sherman Morrison (ISM) decomposition [Preconditioning sparse nonsymmetric linear systems with the Sherman–Morrison formula. Bru (SIAM J Sci Comput 25(2):701–715, 2003), that using recursion formulas derived from the Sherman-Morrison formula, obtains the direct and inverse LU factors of a matrix. We present a modification of the ISM decomposition that allows for pivoting, and so the computation of preconditioners for any nonsingular matrix. While the ISM algorithm at a given step computes only a new pair of vectors, the new pivoting algorithm in the k-th step also modifies all the remaining vectors from k+1 to n. Thus, it can be seen as a right looking version of the ISM decomposition. The results of numerical experiments with ill-conditioned and highly indefinite matrices arising from different applications show the robustness of the new algorithm, since it is able to solve problems that are not possible to solve otherwise.

On the evaluation of general sparse hybrid linear solvers

AbstractGeneral sparse hybrid solvers are commonly used kernels for solving wide range of scientific and engineering problems. This work addresses the current problems of efficiently solving general sparse linear equations with direct/iterative hybrid solvers on many core distributed clusters. We briefly discuss the solution stages of Maphys, HIPS, and PDSLin hybrid solvers for large sparse linear systems with their major algorithmic differences. In this category of solvers, different methods with sophisticated preconditioning algorithms are suggested to solve the trade off between memory and convergence. Such solutions require a certain hierarchical level of parallelism more suitable for modern supercomputers that allow to scale for thousand numbers of processors using Schur complement framework. We study the effect of reordering and analyze the performance, scalability as well as memory for each solve phase of PDSLin, Maphys, and HIPS hybrid solvers using large set of challenging matrices arising from different actual applications and compare the results with SuperLU_DIST direct solver. We specifically focus on the level of parallel mechanisms used by the hybrid solvers and the effect on scalability. Tuning and Analysis Utilities (TAU) is employed to assess the efficient usage of heap memory profile and measuring communication volume. The tests are run on high performance large memory clusters using up to 512 processors.

Adaptive parameter based matrix splitting iteration method for the large and sparse linear systems

For the large and sparse linear systems, we utilize the efficient splittings of the system matrix and introduce an intermediate variable. The main contribution of this paper is that the adaptive parameter based matrix splitting iteration method is constructed from the view of numerical optimization and hybrid procedure to solve the derived equation instead. The novel method adopts the prediction and correction two-step iteration and the global convergence results are established under some suitable conditions. Compared with some well-known methods, numerical experiments show the efficiency and effectiveness of the new proposal with application to the three-dimensional convection-diffusion equation and the image restoration problems.

Data-Driven Fuzzy Transform

The Fuzzy transform is applied mainly to 1-D signals and 2-D data organized as a regular grid (e.g., 2-D images), thus, limiting its potential application to arbitrary data in terms of dimensionality and structure. This article defines and analyzes the properties of the data-driven F-transform, with a focus on the construction of the class of data-driven membership functions, which are multiscale, local, linearly independent, intrinsic, and robust to data discretization. Data-driven membership functions are defined by applying a 1-D filter to the Laplace–Beltrami operator, which encodes the geometric and topological properties of the input data. Then, we address the efficient computation of the data-driven F-transform through a polynomial or a rational polynomial approximation of the input filter. In this way, the computation of the data-driven F-transform is independent of the evaluation of the membership functions at any point of the input domain and reduces to the solution of a small set of sparse and symmetric linear systems. Finally, the data-driven F-transform is efficiently evaluated on large and arbitrary data, in terms of dimensionality, structure, and size.

A time-splitting local meshfree approach for time-fractional anisotropic diffusion equation: application in image denoising

Image denoising approaches based on partial differential modeling have attracted a lot of attention in image processing due to their high performance. The nonlinear anisotropic diffusion equations, specially Perona–Malik model, are powerful tools that improve the quality of the image by removing noise while preserving details and edges. In this paper, we propose a powerful and accurate local meshless algorithm to solve the time-fractional Perona–Malik model which has an adjustable fractional derivative making the control of the diffusion process more convenient than the classical one. In order to overcome the complexities of the problem, a suitable combination of the compactly supported radial basis function method and operator splitting technique is proposed to convert a complex time-fractional partial differential equation into sparse linear algebraic systems that standard solvers can solve. The numerical results of classical and fractional models are explored in different metrics to demonstrate the proposed scheme’s effectiveness. The numerical experiments confirm that the method is suitable to denoise digital images and show that the fractional derivative increases the model’s ability to remove noise in images.

Second‐order accurate hierarchical approximate factorizations for solving sparse linear systems

AbstractWe describe a second‐order accurate approach to sparsifying the off‐diagonal matrix blocks in the hierarchical approximate factorization methods for solving sparse linear systems. These methods repeatedly sparsify the fill‐in matrix blocks that arise in block Gaussian elimination, to compute an approximate factorization of the given matrix, assuming that the fill‐in blocks are low‐rank. The factorization is then used for preconditioning in a Krylov subspace method, such as the conjugate gradient method (CG), BiCGSTAB or GMRES. However, to achieve fast convergence on ill‐conditioned systems, sparsifications can only introduce a small error, in which case they may inefficiently restore sparsity, and consequently the factorization can take a lot of time to compute. In the novel approach to sparsification, the 2‐norm of the incurred error in the sparsification of a matrix block is squared compared to previous approaches, with no additional computations. We incorporate the new approach into the recent sparsified nested dissection algorithm and test it on a wide range of symmetric positive definite problems. The new approach halves the number of CG iterations needed for convergence, significantly improving the overall performance of the algorithm. Our approach can be incorporated into other solvers that exploit the low‐rank property of matrix blocks.

Cross-Interactive Residual Smoothing for Global and Block Lanczos-Type Solvers for Linear Systems with Multiple Right-Hand Sides

Global and block Krylov subspace methods are efficient iterative solvers for large sparse linear systems with multiple right-hand sides. However, global or block Lanczos-type solvers often exhibit large oscillations in the residual norms and may have a large residual gap relating to the loss of attainable accuracy of the approximations. Conventional residual smoothing schemes suppress these oscillations but cannot improve the attainable accuracy, whereas a recent residual smoothing scheme enables the improvement of the attainable accuracy for single right-hand side Lanczos-type solvers. The underlying concept of this scheme is that the primary and smoothed sequences of the approximations and residuals influence one another, thereby avoiding the severe propagation of rounding errors. In the present study, we extend this cross-interactive residual smoothing to the case of solving linear systems with multiple right-hand sides. The resulting smoothed methods can reduce the residual gap with a low additional cost compared to their original counterparts. We demonstrate the effectiveness of the proposed approach through rounding error analysis and numerical experiments.

Compressed basis GMRES on high-performance graphics processing units

Krylov methods provide a fast and highly parallel numerical tool for the iterative solution of many large-scale sparse linear systems. To a large extent, the performance of practical realizations of these methods is constrained by the communication bandwidth in current computer architectures, motivating the investigation of sophisticated techniques to avoid, reduce, and/or hide the message-passing costs (in distributed platforms) and the memory accesses (in all architectures). This article leverages Ginkgo’s memory accessor in order to integrate a communication-reduction strategy into the (Krylov) GMRES solver that decouples the storage format (i.e., the data representation in memory) of the orthogonal basis from the arithmetic precision that is employed during the operations with that basis. Given that the execution time of the GMRES solver is largely determined by the memory accesses, the cost of the datatype transforms can be mostly hidden, resulting in the acceleration of the iterative step via a decrease in the volume of bits being retrieved from memory. Together with the special properties of the orthonormal basis (whose elements are all bounded by 1), this paves the road toward the aggressive customization of the storage format, which includes some floating-point as well as fixed-point formats with mild impact on the convergence of the iterative process. We develop a high-performance implementation of the “compressed basis GMRES” solver in the Ginkgo sparse linear algebra library using a large set of test problems from the SuiteSparse Matrix Collection. We demonstrate robustness and performance advantages on a modern NVIDIA V100 graphics processing unit (GPU) of up to 50% over the standard GMRES solver that stores all data in IEEE double-precision.

Parallel implementation of the SHYFEM (System of HydrodYnamic Finite Element Modules) model

Abstract. This paper presents the message passing interface (MPI)-based parallelization of the three-dimensional hydrodynamic model SHYFEM (System of HydrodYnamic Finite Element Modules). The original sequential version of the code was parallelized in order to reduce the execution time of high-resolution configurations using state-of-the-art high-performance computing (HPC) systems. A distributed memory approach was used, based on the MPI. Optimized numerical libraries were used to partition the unstructured grid (with a focus on load balancing) and to solve the sparse linear system of equations in parallel in the case of semi-to-fully implicit time stepping. The parallel implementation of the model was validated by comparing the outputs with those obtained from the sequential version. The performance assessment demonstrates a good level of scalability with a realistic configuration used as benchmark.

Space-time generalized finite difference method for solving the thin elastic plate bending under dynamic loading

In this paper, the space-time generalized finite difference method (ST-GFDM) with supplementary nodes is applied to solve the thin elastic plate bending under dynamic loading. The proposed method treats the temporal dimension as another spatial dimension, thus the d-dimensional problem in space can be viewed as a new (d+1)-dimensional problem in the space-time domain. This method effectively reduces the complexity of numerical solution of the problem by simultaneously discretizing the temporal and spatial dimensions in the space-time domain. It avoids the temporal-difference measurement while maintaining all the advantages of the GFDM. The ST-GFDM leads a sparse linear system due to its local meshless feature, which is suitable for solving large-scale and long-time bending problems. Moreover, the supplementary nodes on the boundary are introduced to make the system well-determined. Several numerical examples are illustrated to demonstrate the accuracy and stability of the proposed method.

A GPU implementation of the PCG method for large-scale image-based finite element analysis in heterogeneous periodic media

Image-based physics simulations in heterogeneous media at a microscopic scale are a growing trend in various fields within and around Scientific Computing. We consider the scope of Numerical Homogenization, where micro-computed tomography is used to obtain digital models of physical samples, naturally leading to pixel and voxel-based solutions. In this context, the Finite Element Method (FEM) is commonly employed to solve the governing differential equations, via a system of algebraic equations. As image dimensions increase, the memory allocation due to the matrix associated with the FEM quickly becomes unfeasible, even in sparse format. Assembly-free strategies are adopted to reduce memory usage, with the caveat of increased computational cost. The Preconditioned Conjugate Gradient (PCG) method is widely employed to solve this sort of large-scale sparse linear systems , and is fitting to be adapted for assembly-free implementations. This work focuses on a massively parallel PCG solver applied to finite element analyses of heat conduction and linear elasticity on image-based models. Memory-efficiency is one of our main concerns, in an attempt to make feasible the employment of personal-use GPUs for large-scale simulations. The resulting solver is validated with an analytical benchmark, and by comparing the obtained results for a microtomographic model of a cast iron sample against experimental values found in the literature. Time and memory metrics are presented and discussed. It is shown that the developed program allows for homogenization studies of nearly 500 million degrees-of-freedom to be conducted in personal computers equipped with CUDA-enabled devices of 8 GB RAM, taking seconds or a few minutes per system solution with the PCG method. Up to 400× speed-up was observed in comparison to an analogous solver running in a 16-thread CPU. Our GPU implementation makes it possible to conduct, in a matter of minutes, homogenization studies that would take hours, or even days, in personal CPUs.

Computer big data modeling system based on finite element mathematical equation simulation

Abstract In the traditional computer analysis system, the equipment is cumbersome and takes up space, and is gradually being eliminated by the market. In view of the characteristics of small size, large throughput and high precision of computer systems in the era of big data, in order to solve the phenomenon of swell disasters induced by high-speed moving landslides, debris flows and glaciers impacting water bodies in enclosed waters, a mathematical equation based on finite element was developed. The sparse linear multi-thread parallel solving system of the group, using the sparse linear equation solver implemented on the shared memory machine provided by Intel MKL, completed the design of the high-performance finite element analysis program for surge. The theoretical solution and numerical calculation results of vertical dam failure and laboratory experimental results show that the normalized error between the theoretical value calculated by the calculation system designed in this paper and the actual value obtained from the experiment is only 5.2%~6.8% In between, it is verified that the system has rapid reliability and practical engineering significance.

AAQAL: A Machine Learning-Based Tool for Performance Optimization of Parallel SPMV Computations Using Block CSR

The sparse matrix–vector product (SpMV), considered one of the seven dwarfs (numerical methods of significance), is essential in high-performance real-world scientific and analytical applications requiring solution of large sparse linear equation systems, where SpMV is a key computing operation. As the sparsity patterns of sparse matrices are unknown before runtime, we used machine learning-based performance optimization of the SpMV kernel by exploiting the structure of the sparse matrices using the Block Compressed Sparse Row (BCSR) storage format. As the structure of sparse matrices varies across application domains, optimizing the block size is important for reducing the overall execution time. Manual allocation of block sizes is error prone and time consuming. Thus, we propose AAQAL, a data-driven, machine learning-based tool that automates the process of data distribution and selection of near-optimal block sizes based on the structure of the matrix. We trained and tested the tool using different machine learning methods—decision tree, random forest, gradient boosting, ridge regressor, and AdaBoost—and nearly 700 real-world matrices from 43 application domains, including computer vision, robotics, and computational fluid dynamics. AAQAL achieved 93.47% of the maximum attainable performance with a substantial difference compared to in practice manual or random selection of block sizes. This is the first attempt at exploiting matrix structure using BCSR, to select optimal block sizes for the SpMV computations using machine learning techniques.

Three-dimensional finite-difference finite-element frequency-domain wave simulation with multi-level optimized additive Schwarz domain-decomposition preconditioner: A tool for FWI of sparse node datasets

Efficient frequency-domain full-waveform inversion (FWI) of long-offset node data can be designed with a few discrete frequencies, which lead to modest data volumes to be managed during the inversion process. Moreover, attenuation effects can be straightforwardly implemented in the forward problem without the computational overhead. However, 3D frequency-domain seismic modeling is challenging because it requires solving a large and sparse linear indefinite system for each frequency with multiple right-hand sides (RHSs). This linear system can be solved by direct or iterative methods. The former allows efficient processing of multiple RHSs but may suffer from limited scalability for very large problems. Iterative methods equipped with a domain-decomposition preconditioner provide a suitable alternative to process large computational domains for sparse-node acquisition. We have investigated the domain-decomposition preconditioner based on the optimized restricted additive Schwarz (ORAS) method, in which a Robin or perfectly matched layer condition is implemented at the boundaries between the subdomains. The preconditioned system is solved by a Krylov subspace method, whereas a block low-rank lower-upper decomposition of the local matrices is performed at a preprocessing stage. Multiple sources are processed in groups with a pseudoblock method. The accuracy, the computational cost, and the scalability of the ORAS solver are assessed against several realistic benchmarks. In terms of discretization, we compare a compact wavelength-adaptive 27-point finite-difference stencil on a regular Cartesian grid with a P3 finite-element method on h-adaptive tetrahedral mesh. Although both schemes have comparable accuracy, the former is more computationally efficient, the latter being beneficial to comply with known boundaries such as bathymetry. The scalability of the method, the block processing of multiple RHSs, and the straightforward implementation of attenuation, which further improves the convergence of the iterative solver, make the method a versatile forward engine for large-scale 3D FWI applications from sparse node data sets.

A General Alternating-Direction Implicit Framework with Gaussian Process Regression Parameter Prediction for Large Sparse Linear Systems

This paper proposes an efficient general alternating-direction implicit (GADI) framework for solving large sparse linear systems. The convergence property of the GADI framework is discussed. Most of the existing ADI methods can be viewed as particular schemes of the developed framework. Meanwhile the GADI framework can derive new ADI methods. Moreover, as the algorithm efficiency is sensitive to the splitting parameters, we offer a data-driven approach, the Gaussian process regression (GPR) method based on the Bayesian inference, to predict the GADI framework's relatively optimal parameters. The GPR method requires a small training data set to learn the regression prediction mapping, which has sufficient accuracy and high generalization capability. It allows us to efficiently solve linear systems with a one-shot computation, and does not require any repeated computations to obtain relatively optimal splitting parameters. Finally, we use the three-dimensional convection-diffusion equation and continuous Sylvester matrix equation to examine the performance of our proposed methods. Numerical results demonstrate that the proposed framework is faster tens to thousands of times than the existing ADI methods, such as (inexact) Hermitian and skew-Hermitian splitting type methods in which the consumption of obtaining relatively optimal splitting parameters is ignored. Due to the efficiency of the developed methods, we can solve much larger linear systems which these existing ADI methods have been not reached.

Local knot method for solving inverse Cauchy problems of Helmholtz equations on complicated two‐ and three‐dimensional domains

AbstractThis article presents a local knot method (LKM) to solve inverse Cauchy problems of Helmholtz equations in arbitrary 2D and 3D domains. The Moore–Penrose pseudoinverse using the truncated singular value decomposition is employed in the local approximation of supporting domain instead of the moving least squares method. The developed approach is a semi‐analytical and local radial basis function collocation method using the non‐singular general solution as the basis function. Like the traditional boundary knot method, the LKM is simple, accurate and easy‐to‐program in solving inverse Cauchy problems associated with Helmholtz equations. Unlike the boundary knot method, the new scheme can directly reconstruct the unknowns both inside the physical domain and along its boundary by solving a sparse linear system, and can achieve a satisfactory solution. Numerical experiments, involving the complicated geometry and the high noise level, confirm the accuracy and reliability of the proposed method for solving inverse Cauchy problems of Helmholtz equations.

Comparative study of inner–outer Krylov solvers for linear systems in structured and high-order unstructured CFD problems

Advanced Krylov subspace methods are investigated for the solution of large sparse linear systems arising from stiff adjoint-based aerodynamic shape optimization problems. A special attention is paid to the flexible inner–outer GMRES strategy combined with most relevant preconditioning and deflation techniques. The choice of this specific class of Krylov solvers for challenging problems is based on its outstanding convergence properties. Typically in our implementation the efficiency of the preconditioner is enhanced with a domain decomposition method with overlapping. However, maintaining the performance of the preconditioner may be challenging since scalability and efficiency of a preconditioning technique are properties often antagonistic to each other. In this paper we demonstrate how flexible inner–outer Krylov methods are able to overcome this critical issue. A numerical study is performed considering either a Finite Volume (FV), or a high-order Discontinuous Galerkin (DG) discretization which affect the arithmetic intensity and memory-bandwidth of the algebraic operations. We consider test cases of transonic turbulent flows with RANS modeling over the two-dimensional supercritical OAT15A airfoil and the three-dimensional ONERA M6 wing. Benefits in terms of robustness and convergence compared to standard GMRES solvers are obtained. Strong scalability analysis shows satisfactory results. Based on these representative problems a discussion of the recommended numerical practices is proposed.

Strategies of Preconditioner Updates for Sequences of Linear Systems Associated with the Neutron Diffusion

The time-dependent neutron diffusion equation approximates the neutronic power evolution inside a nuclear reactor core. Applying a Galerkin finite element method for the spatial discretization of these equations leads to a stiff semi-discrete system of ordinary differential equations. For time discretization, an implicit scheme is used, which implies solving a large and sparse linear system of equations for each time step. The GMRES method is used to solve these systems because of its fast convergence when a suitable preconditioner is provided. This work explores several matrix-free strategies based on different updated preconditioners, which are constructed by low-rank updates of a given initial preconditioner. They are two tuned preconditioners based on the bad and good Broyden’s methods, initially developed for nonlinear equations and optimization problems, and spectral preconditioners. The efficiency of the resulting preconditioners under study is closely related to the selection of the subspace used to construct the update. Our numerical results show the effectiveness of these methodologies in terms of CPU time and storage for different nuclear benchmark transients, even if the initial preconditioner is not good enough.

Adaptive identification of sparse underwater acoustic channels with a mix of static and time-varying parameters

We consider identification of sparse linear systems with a mix of static and time-varying parameters. Such systems are typical in underwater acoustics (UWA), for instance, in applications requiring identification of the acoustic channel, such as UWA communications, navigation and continuous-wave sonar. The recently proposed fast local basis function (fLBF) algorithm provides high performance when identifying time-varying systems. In this paper, we further improve the performance of the fLBF algorithm by exploiting properties of the system. Specifically, we propose an adaptive time-invariance test to identify whether a particular system tap is static or time-varying and exploit this knowledge for choosing the number of basis functions. We also propose a regularization scheme that exploits the system sparsity and an adaptive technique for estimating the regularization parameter. Finally, a debiasing technique is proposed to reduce an inherent bias of fLBF estimates. The high performance of the fLBF algorithm with the proposed techniques is demonstrated in scenarios of UWA communications, using numerical and real experiments.