Incomplete Cholesky Research Articles

On cost-efficient parallel iterative solvers for 3D frequency-domain seismic multisource viscoelastic anisotropic wave modeling

Solving large sparse linear systems in 3D frequency-domain seismic wave modeling, especially in viscoelastic anisotropic media, poses significant challenges due to the increasing number of discrete moduli and nonzero elements in the linear system matrix. The computational load surpasses that of acoustic or viscoacoustic media, making it even more challenging when dealing with multisource problems. Popular scientific tools for solving a linear system, such as multifrontal massively parallel sparse direct solver (MUMPS), structured matrix package (STRUMPACK), and portable extensible toolkit for scientific computation (PETSc), can be used, but their applicability to our specific problem has not been comprehensively evaluated. Our study aims to tackle the challenges in solving large, sparse, complex-valued symmetric linear systems with multiple right-hand side vectors for 3D frequency-domain seismic wave modeling. We have adopted the preconditioned conjugate gradient iterative algorithms as the foundation for our research, introducing two highly cost-effective parallel iterative solvers: the parallel symmetric successive overrelaxation conjugate gradient (P-SSORCG) and the parallel incomplete Cholesky conjugate gradient (P-ICCG). These novel solvers are subjected to a comprehensive comparative analysis against well-established scientific tools, such as MUMPS, STRUMPACK, and PETSc, in the context of 3D frequency-domain seismic wave modeling. We show their promising performances in a practical 3D SEG/EAGE overthrust model and demonstrate that the grouped P-SSORCG offers an efficient alternative to parallel direct solvers, particularly in situations wherein computational resources are limited.

A non-conformal domain decomposition method utilizing rotating subdomains and non-matching grids for periodic metamaterial simulation

Electromagnetic metamaterials possess characteristics such as having large-scale, periodic, multi-media, and complex geometric structures, making them highly suitable for simulation using the finite element domain decomposition method (DDM). This work proposes a non-conformal DDM for simulating finite-period metamaterials, achieving simulation with a minimal number of domains while significantly reducing computational resource consumption. First, the computational domain of the periodic metamaterial is partitioned into several non-conformal hexahedral subdomains. Subsequently, a limited number of non-repetitive subdomain models with non-matching grids are constructed, leveraging the three-dimensional rotational and translational properties of these subdomains to facilitate simulation. By eliminating the necessity to construct models for every subdomain, a substantial reduction in memory consumption is achieved. Second, an iteration method for solving the matrix equations of subdomains is enhanced by introducing a multifrontal block incomplete Cholesky decomposition preconditioner, thereby enhancing the computational efficiency of matrix equations with a large number of unknowns. Meanwhile, parallel computing techniques are employed to accelerate the proposed method. Finally, we integrate the aforementioned method into a solver and leverage it to develop an electromagnetic simulation tool tailored for metamaterials. The tool is employed to simulate metamaterial structures of varying scales, resulting in notable reductions in both memory and time consumption while maintaining accuracy comparable to commercial software.

Avoiding breakdown in incomplete factorizations in low precision arithmetic

The emergence of low precision floating-point arithmetic in computer hardware has led to a resurgence of interest in the use of mixed precision numerical linear algebra. For linear systems of equations, there has been renewed enthusiasm for mixed precision variants of iterative refinement. We consider the iterative solution of large sparse systems using incomplete factorization preconditioners. The focus is on the robust computation of such preconditioners in half precision arithmetic and employing them to solve symmetric positive definite systems to higher precision accuracy; however, the proposed ideas can be applied more generally. Even for well-conditioned problems, incomplete factorizations can break down when small entries occur on the diagonal during the factorization. When using half precision arithmetic, overflows are an additional possible source of breakdown. We examine how breakdowns can be avoided and we implement our strategies within new half precision Fortran sparse incomplete Cholesky factorization software. Results are reported for a range of problems from practical applications. These demonstrate that, even for highly ill-conditioned problems, half precision preconditioners can potentially replace double precision preconditioners, although unsurprisingly this may be at the cost of additional iterations of a Krylov solver.

Coupled Incomplete Cholesky and Jacobi Preconditioned Conjugate Gradient on the New Generation of Sunway Many-Core Architecture

Coupled Incomplete Cholesky and Jacobi Preconditioned Conjugate Gradient on the New Generation of Sunway Many-Core Architecture

Unleashing the Power of Graph Spectral Sparsification for Power Grid Analysis via Incomplete Cholesky Factorization

Graph spectral sparsification based preconditioning technique has shown promising results for power grid analysis. However, the conventional methods converge slowly for high accuracy requirement. In this work, we propose an efficient approach to address this issue. Instead of using Cholesky factorization, we employ incomplete Cholesky factorization to factorize the spectral sparsifier. We also propose a concept of graph spectral pattern, which can further reduce the preconditioned conjugate gradient (PCG) iterations using less number of nonzeros. Experiments show that under 10-6 relative tolerance, our proposed preconditioning technique achieves 1.17× speedup compared to AMGPCG in average; compared to the conventional spectral sparsification based preconditioning techniques, our proposed approach achieves up to 8.53× speedup of the factorization, 8.74× speedup of the PCG iteration, and 5.6× speedup of the total time. Moreover, the speedup of the total time continues to enlarge for higher accuracy requirement, e.g., 10-12. Last, but not the least, our method is compatible with existing graph spectral sparsification algorithms for power grid analysis.

GPU implementation of an incomplete Cholesky conjugate gradient solver for a FEM-generated system using full kernel consolidation

GPU implementation of an incomplete Cholesky conjugate gradient solver for a FEM-generated system using full kernel consolidation

An Incomplete Cholesky Preconditioner Based on Orthogonal Approximations

An Incomplete Cholesky Preconditioner Based on Orthogonal Approximations

Fast iterative solver for 3D finite-element modeling of controlled-source electromagnetic responses with vector and scalar potentials

The controlled-source electromagnetic (CSEM) forward modeling with vector and scalar potentials could provide insights into the inductive and galvanic effects. We have developed a finite-element (FE) forward modeling algorithm for the CSEM with vector and scalar potentials. The first-order vector and nodal shape functions are used for the vector and scalar potentials, respectively. To mitigate the null space and thus the nonuniqueness of the potentials caused by the curl operator, we introduce a new preconditioner constructed from an incomplete Cholesky decomposition with zero fill-ins of the Laplacian approximation to the resulted matrix after the FE discretization. We implemented a preconditioned quasiminimal residual method to iteratively solve the resulting linear system with the new preconditioner. We first verified the accuracy of the algorithm through comparison against the analytic solutions obtained for a three-layered earth model. The new iterative modeling algorithm converges much faster compared with the modeling algorithm with the preconditioner constructed from the original stiffness matrix. The efficiency and accuracy are further illustrated by comparison with the modeling scheme based on coupled potentials with explicit enforcement of the Coulomb gauge condition for the vector potential by modeling three 3D models. In addition, the accuracy of our algorithm is demonstrated via a comparison of our numerical solution for a 3D model with the numerical solution obtained from the integral equation method. We also present the inductive and galvanic components of the horizontal secondary electric field for numerical solutions obtained with direct and iterative solvers. The numerical test demonstrated that, aside from the efficiency improved with the new iterative solver, the unstable and nonuniqueness problem for the potentials is eliminated by the new preconditioner with the Laplacian approximation involved.

Fast finite element electrostatic analysis with domain decomposition method

Accurate simulation of electrostatic field in electronic devices is very important to improve their performance. When the traditional finite element method is used to analyze the electrostatic field of large and complex structures, the calculation speed is slow and the computing resources are greatly consumed. Therefore, in this paper, we investigate electrostatic analysis with a non-overlapping domain decomposition method based on the scalar finite element method under a novel transmission condition. In order to improve the efficiency of matrix solving and reduce the memory consumption, the block Jacobi preconditioner is introduced. In addition, the multifrontal block incomplete Cholesky preconditioner is utilized to further improve the computational performance based on the block Jacobi preconditioner. Several numerical examples are simulated and compared with the analytical solutions and the commercial software Maxwell 3D. The results demonstrate the accuracy and efficiency of the proposed domain decomposition method.

Neural-network preconditioners for solving the Dirac equation in lattice gauge theory

This work develops neural-network--based preconditioners to accelerate solution of the Wilson-Dirac normal equation in lattice quantum field theories. The approach is implemented for the two-flavor lattice Schwinger model near the critical point. In this system, neural-network preconditioners are found to accelerate the convergence of the conjugate gradient solver compared with the solution of unpreconditioned systems or those preconditioned with conventional approaches based on even-odd or incomplete Cholesky decompositions, as measured by reductions in the number of iterations and/or complex operations required for convergence. It is also shown that a preconditioner trained on ensembles with small lattice volumes can be used to construct preconditioners for ensembles with many times larger lattice volumes, with minimal degradation of performance. This volume-transferring technique amortizes the training cost and presents a pathway towards scaling such preconditioners to lattice field theory calculations with larger lattice volumes and in four dimensions.

Параллельная реализация метода сопряженных градиентов с предобусловливателем IC1 на основе использования переупорядочения узлов сетки

The paper proposes a method for using MPI+OpenMP technology for constructing and inverting a first-order incomplete triangular Cholesky expansion preconditioner IC1(τ) for solving system of linear algebraic equations with an arbitrary symmetric positive definite matrix on a small number of processors. Methods of using MPI and MPI+OpenMP technologies are based on the use of grid node orderings consistent with the division of the calculation area. When constructing the preconditioning matrix IC1(τ) using MPI, cutting is performed at positions in some of its rows. The use of OpenMP technology in the construction and inversion of the preconditioner is carried out for most rows of the matrix. The time taken to solve problems using the conjugate gradient method with the IC1(τ) preconditioner using MPI and hybrid MPI+OpenMP technology is compared using the example of a model problem and a number of problems from the SuiteSparse collection of sparse matrices.

State Estimation for Closed-Loop LPV System Identification via Kernel Methods⋆

This paper proposes a state-estimation method for closed-loop identification of linear-parameter-varying (LPV) systems by extending an approach for linear-time-invariant (LTI) systems and applying kernel canonical correlation analysis (KCCA). The proposed method estimates the state term from the one-step-ahead prediction via a kernel approach. The incomplete Cholesky decomposition (ICD) is introduced to reduce the complexity of the KCCA. A simple numerical simulation shows the effectiveness of the proposed approach.

Linear programming with a feasible direction interior point technique for smooth optimization

We propose an adaptation of the Feasible Direction Interior Points Algorithm (FDIPA) of J. Herskovits, for solving large-scale linear programs. At each step, the solution of two linear systems with the same coefficient matrix is determined. This step involves a significant computational effort. Reducing the solution time of linear systems is, therefore, a way to improve the performance of the method. The linear systems to be solved are associated with definite positive symmetric matrices. Therefore, we use Split Preconditioned Conjugate Gradient (SPCG) method to solve them, together with an Incomplete Cholesky preconditioner using Matlab’s ICHOL function. We also propose to use the first iteration of the conjugate gradient, and to presolve before applying the algorithm, in order to reduce the computational cost. Following, we then provide mathematica proof that show that the iterations approach Karush–Kuhn–Tucker points of the problem under reasonable assumptions. Finally, numerical evidence show that the method not only works in theory but is also competitive with more advanced methods.

Adaptive rigidification of elastic solids

We present a method for reducing the computational cost of elastic solid simulation by treating connected sets of non-deforming elements as rigid bodies. Non-deforming elements are identified as those where the strain rate squared Frobenius norm falls below a threshold for several frames. Rigidification uses a breadth first search to identify connected components while avoiding connections that would form hinges between rigid components. Rigid elements become elastic again when their approximate strain velocity rises above a threshold, which is fast to compute using a single iteration of conjugate gradient with a fixed Laplacian-based incomplete Cholesky preconditioner. With rigidification, the system size to solve at each time step can be greatly reduced, and if all elastic element become rigid, it reduces to solving the rigid body system. We demonstrate our results on a variety of 2D and 3D examples, and show that our method is likewise especially beneficial in contact rich examples.

A New Approach to the Splitting Factor Preconditioner Applied to Linear Programming Problems

In this paper, we present the results of a new approach to the splitting factor preconditioner, which is a preconditioner based on the Incomplete Cholesky factorization and the splitting preconditioner. In previous work for small linear programming problems, the preconditioner was applied in all iterations of the interior point method and compared with the splitting preconditioner also applied in all iterations. In this paper, we will do a hybrid approach, in which in the first iterations the preconditioner is the Incomplete Cholesky Factorization, and in the last iterations, the preconditioner used is the splitting factor preconditioner or the splitting preconditioner. The results obtained show that even in the hybrid approach, the splitting factor preconditioner achieved better performance.

Using incomplete Cholesky factorization to increase the time step in molecular dynamics simulations

Incomplete Cholesky (IC) factorization is widely used as a preconditioner for accelerating the convergence of the conjugate gradient iterative method. The IC factorization has been used in continuum mechanics and other applications that require solutions of elliptic partial differential equations. In this study, we propose an efficient use of the IC factorization to increase the time step and accelerate molecular dynamics simulations. Previously, we proposed the semi-implicit Hessian correction (SimHec) scheme (Washio et al., 2021) for overdamped Langevin dynamics of polymer simulations. SimHec constructs an approximation of the Hessian matrix by superposing the corrected elemental Hessian matrices associated with the interactions within a limited bandwidth along the polymer chain. The resulting narrow-bandwidth system is efficiently solved by an overlapped skyline solver in parallel. In this study, we integrate the non-local interactions in the polymer chain into the corrected Hessian matrix, and approximate the components of the Hessian matrix outside the narrow bandwidth using the IC factorization. This strategy allows us to use time steps that are 150–1000 times larger than those used in the explicit scheme without a severe increase in the computational load.

MPI+OpenMP реализация метода сопряженных градиентов с предобусловливателем IC(0) на основе использования переупорядочения узлов сетки

The paper proposes methods for using MPI and MPI+OpenMP technologies to construct and invert the preconditioner of an incomplete triangular Cholesky decomposition without filling IC(0) for solving systems of linear algebraic equations with an arbitrary symmetric positive definite matrix. Methods of using MPI and MPI+OpenMP technologies are based on the use of grid node orderings consistent with the division of the calculation area. The use of OpenMP technology in the construction and inversion of the preconditioner is carried out for most rows of the matrix. Comparative timing results for the MPI+OpenMP and MPI implementations of the proposed preconditioning used with the conjugate gradient method for a model problems and the sparse matrix collections SuiteSparse are presented.

Analysis of SAR reduction to human head with plasma photonic crystals shield using ICCG-SFDTD method

A novel and effective biological electromagnetic method using plasma photonic crystals shield to protect the human head from electromagnetic radiation of an ultra-high-field magnetic resonance imaging (MRI) frequency is explored. An improved incomplete Cholesky conjugate gradient (ICCG)-based symplectic finite-difference time-domain (ICCG-SFDTD) method, which maintains the property of the Hamilton system and ensures energy conservation numerically with high precision, is applied to a bioelectromagnetic simulation of human head models without shield, with plasma shield and with plasma photonic crystals shield. Results show that plasma photonic crystals shield is superior to plasma shield with the same thickness in reducing the specific absorption rate (SAR) and electric field levels of the human head. In addition, the relationship of SAR and electric field in the human head models with the parameters of plasma photonic crystals, such as the plasma frequency, plasma collision frequency and dielectric constant of the medium is studied in order to provide a guidance to design plasma photonic crystals protective layer for biological tissue.

Three-dimensional controlled-source electromagnetic forward modeling by edge-based finite element with a divergence correction

We have developed an edge-based finite-element (FE) modeling algorithm with a divergence correction for calculating the controlled-source electromagnetic (CSEM) responses of a 3D conductivity earth model. We solve a curl-curl equation to directly calculate the secondary electric field to eliminate the source singularity. The choice of the edge-based FE method enables us to properly handle the discontinuity of the normal component of electric fields across conductivity boundaries. Although we can solve the resulting complex-symmetric linear system of equations efficiently by a quasiminimal residual (QMR) method preconditioned with an incomplete Cholesky decomposition for the high-frequency band, the iterative solution process encounters a common problem in the field formulation and does not converge within a practically feasible number of iterations for low frequencies. To overcome this difficulty and to accelerate the iterative solution process in general, we combine a divergence correction technique with the secondary field solution using the QMR solver. We have found that applying the divergence correction intermittently during the iterative solution process ensures the calculation of sufficiently accurate electric and magnetic fields and can significantly speed up the solution process by more than an order of magnitude. We have tested the efficiency and accuracy of our algorithm with 1D and 3D models, and we have found that the divergence correction technique is able to guide the electric field to satisfy the boundary conditions across conductivity interfaces. Although there is a computational overhead required for applying the divergence correction, that cost is significantly offset by the substantial gains in the solution accuracy and speed-up. The work makes the field-based curl-curl formulation using edge elements an efficient and practical method for CSEM simulations.

Acceleration and Parallelization of a Linear Equation Solver for Crack Growth Simulation Based on the Phase Field Model

We aim to accelerate the linear equation solver for crack growth simulation based on the phase field model. As a first step, we analyze the properties of the coefficient matrices and prove that they are symmetric positive definite. This justifies the use of the conjugate gradient method with the efficient incomplete Cholesky preconditioner. We then parallelize this preconditioner using so-called block multi-color ordering and evaluate its performance on multicore processors. The experimental results show that our solver scales well and achieves an acceleration of several times over the original solver based on the diagonally scaled CG method.