Preconditioned Conjugate Gradient Solver Research Articles

A matching Schur complement preconditioning technique for inverse source problems

Numerical discretization of a regularized inverse source problem leads to a non-symmetric saddle point linear system. Interestingly, the Schur complement of the non-symmetric saddle point system is Hermitian positive definite (HPD). Then, we propose a preconditioner matching the Schur complement (MSC). Theoretically, we show that the preconditioned conjugate gradient (PCG) method for a linear system with the preconditioned Schur complement as coefficient has a linear convergence rate independent of the matrix size and value of the regularization parameter involved in the inverse problem. Fast implementations are proposed for the matrix-vector multiplication of the preconditioned Schur complement so that the PCG solver requires only quasi-linear operations. To the best of our knowledge, this is the first solver with guarantee of linear convergence for the inversion of Schur complement arising from the discrete inverse problem. Combining the PCG solver for inversion of the Schur complement and the fast solvers for the forward problem in the literature, the discrete inverse problem (the saddle point system) is solved within a quasi-linear complexity. Numerical results are reported to show the performance of the proposed matching Schur complement (MSC) preconditioning technique.

A low-rank isogeometric solver based on Tucker tensors

We propose an isogeometric solver for Poisson problems that combines (i) low-rank tensor techniques to approximate the unknown solution and the system matrix, as a sum of a few terms having Kronecker product structure, (ii) a Truncated Preconditioned Conjugate Gradient solver to keep the rank of the iterates low, and (iii) a novel low-rank preconditioner, based on the Fast Diagonalization method where the eigenvector multiplication is approximated by the Fast Fourier Transform. Although the proposed strategy is written in arbitrary dimension, we focus on the three-dimensional case and adopt the Tucker format for low-rank tensor representation, which is well suited in low dimension. We show by numerical tests that this choice guarantees significant memory saving compared to the full tensor representation. We also extend and test the proposed strategy to linear elasticity problems.

A kernel-free boundary integral method for the nonlinear Poisson-Boltzmann equation

This work proposes a boundary integral formulation based Cartesian grid method for the nonlinear Poisson-Boltzmann (PB) interface problem in biophysics. The method (a) does not have the limitation associated with the standard finite difference method for complex interfaces, (b) avoids the generation of any unstructured volume or surface grids needed by the finite element method, and (c) solves the nonlinear and variable coefficient PDE in the framework of boundary integral equations. The method solves the nonlinear PB equation with the Newton iterative method. It first reformulates the linearized PB equation in each Newton iteration as a Fredholm system of two boundary integral (BI) equations of the second kind, which is well-conditioned, and then solves the BI system with a Krylov subspace method. The evaluation of boundary and volume integrals involved in the solution of the BI system is done with a kernel-free method, which does not need to know Green's function associated with the variable coefficient PB equation. The kernel-free method evaluates volume and boundary integrals by solving equivalent simple interface problems on Cartesian grids with either a fast Fourier transform based Poisson solver or a geometric multigrid preconditioned conjugate gradient solver, either of which involves computational work essentially linearly proportional to the number of nodes on the Cartesian grid. Numerical examples in both two and three space dimensions are presented to demonstrate the efficiency and accuracy of the proposed numerical method.

Convergence acceleration of preconditioned conjugate gradient solver based on error vector sampling for a sequence of linear systems

AbstractIn this article, we focus on solving a sequence of linear systems that have identical (or similar) coefficient matrices. For this type of problem, we investigate subspace correction (SC) and deflation methods, which use an auxiliary matrix (subspace) to accelerate the convergence of the iterative method. In practical simulations, these acceleration methods typically work well when the range of the auxiliary matrix contains eigenspaces corresponding to small eigenvalues of the coefficient matrix. We develop a new algebraic auxiliary matrix construction method based on error vector sampling in which eigenvectors with small eigenvalues are efficiently identified in the solution process. We use the generated auxiliary matrix for convergence acceleration in the following solution step. Numerical tests confirm that both SC and deflation methods with the auxiliary matrix can accelerate the solution process of the iterative solver. Furthermore, we examine the applicability of our technique to the estimation of the condition number of the coefficient matrix. We also present the algorithm of the preconditioned conjugate gradient method with condition number estimation.

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.

Bernstein-Bézier Galerkin-Characteristics Finite Element Method for Convection-Diffusion Problems

A class of Bernstein-Bézier basis based high-order finite element methods is developed for the Galerkin-characteristics solution of convection-diffusion problems. The Galerkin-characteristics formulation is derived using a semi-Lagrangian discretization of the total derivative in the considered problems. The spatial discretization is performed using the finite element method on unstructured meshes. The Lagrangian interpretation in this approach greatly reduces the time truncation errors in the Eulerian methods. To achieve high-order accuracy in the Galerkin-characteristics solver, the semi-Lagrangian method requires high-order interpolating procedures. In the present work, this step is carried out using the Bernstein-Bézier basis functions to evaluate the solution at the departure points. Triangular Bernstein-Bézier patches are constructed in a simple and inherent manner over finite elements along the characteristics. An efficient preconditioned conjugate gradient solver is used for the linear systems of algebraic equations. Several numerical examples including advection-diffusion equations with known analytical solutions and the viscous Burgers problem are considered to illustrate the accuracy, robustness and performance of the proposed approach. The computed results support our expectations for a stable and highly accurate Bernstein-Bézier Galerkin-characteristics finite element method for convection-diffusion problems.

A semi-Lagrangian Bernstein–Bézier finite element method for two-dimensional coupled Burgers’ equations at high Reynolds numbers

This paper aims to develop a semi-Lagrangian Bernstein–Bézier high-order finite element method for solving the two-dimensional nonlinear coupled Burgers’ equations at high Reynolds numbers. The proposed method combines the semi-Lagrangian scheme for the time integration and the high-order Bernstein–Bézier functions for the space discretization in the finite element framework. Unstructured triangular Bernstein–Bézier patches are reconstructed in a simple and inherent manner over finite elements along the characteristic curves defined by the material derivative. A fourth-order Runge–Kutta scheme is used for the approximation of departure points along with a local L2-projection to compute the solution at the semi-Lagrangian stage. By using these techniques, the nonlinear problem is decoupled and two linear diffusion problems are solved separately for each velocity component. An implicit time-steeping scheme is used and a preconditioned conjugate gradient solver is used for the resulting linear systems of algebraic equations. The proposed method is investigated through several numerical examples including convergence studies. It is found that the proposed method is stable, highly accurate and efficient in solving two-dimensional coupled Burgers’ equations at high Reynolds numbers.

An efficient matrix-free preconditioned conjugate gradient based multigrid method for phase field modeling of fracture in heterogeneous materials from 3D images

In this paper, we present a new and efficient strategy to perform 3D computational fracture modeling from computed tomography (CT) images. The image-based fracture modeling creates a new way to more accurately predict fracture in realistic structures. However, it is currently complex and expensive to perform such simulations. In this work, we propose to use a matrix-free type preconditioned conjugate gradient solver based multigrid method to achieve the fastest numerical performance for phase field modeling of fracture. Several specific methods are investigated to enhance the robustness and to improve the efficiency of the proposed strategy. An automatic load control strategy to control crack propagation is investigated. High performance computing is applied using a hybrid MPI/OpenMP strategy. With all of the proposed procedures, the phase field modeling of fracture can be performed in the real 3D microstructure of heterogeneous materials using tomographic images.

Robust Discretization and Solvers for Elliptic Optimal Control Problems with Energy Regularization

Abstract We consider elliptic distributed optimal control problems with energy regularization. Here the standard L 2 {L_{2}} -norm regularization is replaced by the H - 1 {H^{-1}} -norm leading to more focused controls. In this case, the optimality system can be reduced to a single singularly perturbed diffusion-reaction equation known as differential filter in turbulence theory. We investigate the error between the finite element approximation u ϱ ⁢ h {u_{\varrho h}} to the state u and the desired state u ¯ {\overline{u}} in terms of the mesh-size h and the regularization parameter ϱ. The choice ϱ = h 2 {\varrho=h^{2}} ensures optimal convergence the rate of which only depends on the regularity of the target function u ¯ {\overline{u}} . The resulting symmetric and positive definite system of finite element equations is solved by the conjugate gradient (CG) method preconditioned by algebraic multigrid (AMG) or balancing domain decomposition by constraints (BDDC). We numerically study robustness and efficiency of the AMG preconditioner with respect to h, ϱ, and the number of subdomains (cores) p. Furthermore, we investigate the parallel performance of the BDDC preconditioned CG solver.

Topology optimization of Stokes flow with traction boundary conditions using low-order finite elements

We consider topology optimization of Stokes flow with traction boundary conditions using finite elements with low-order velocity-approximation and an element-wise constant hydrostatic pressure. The finite element formulation is stabilized using a penalty on the jump in pressure between adjacent elements. Convergence of solutions to the finite element-discretized topology optimization problem is shown, and several optimization problems are solved using a preconditioned conjugate gradient solver for the finite element matrix problem. Stable convergence to high-quality designs without an excessive number of linear solver iterations is observed, and it is seen that the finite element formulation is not particularly sensitive to the choice of the pressure jump penalty parameter, thus making it a practically useful method.

Three-dimensional magnetotelluric modeling in general anisotropic media using nodal-based unstructured finite element method

Magnetotelluric sounding is an efficient and economical method for detecting Earth's deep structure because of its large surveying depth. We use a finite element approach to solve magnetotelluric fields in three-dimensional general anisotropic media based on the vector–scalar potentials and unstructured meshes. The implementation of a three-dimensional unstructured mesh generator TetGen with high-quality local mesh refinement can easily handle complicated models with high performance. The electromagnetic field equations are expressed in terms of Coulomb-gauged vector–scalar potentials, and the nodal-based finite element method can be used to compute the electromagnetic fields with high efficiency. The symmetric successive over-relaxation preconditioned conjugate gradient solver is used to solve the resulting linear equation system. The magnetotelluric responses in one-, two-, and three-dimensional medias are computed through our proposed approach. These results are then compared with the analytical solutions, two- and three-dimensional adaptive finite element solutions. Numerical examples demonstrate that these results are in good agreement. Because an unstructured mesh is incorporated in this approach, it has great potential for simulating the magnetotelluric responses in geometrically complicated situations.

Time-Stepped Finite-Element Modeling of Three-Phase Transformer for Electromagnetic Transient Emulation on FPGA

The finite-element analysis is a powerful method to obtain detailed insight into the operation of any electromagnetic equipment. However, the required computational power to solve a finite-element modeled power equipment is so heavy that most Newton-Raphson-method-based algorithms can barely achieve real-time simulation. The low latency and hardware parallelism of the field-programmable gate array (FPGA) provides a path forward. In this paper, a parallel and deeply pipelined adaptive transmission-line modeling method with preconditioned conjugate gradient solver is designed in hardware and implemented on two Xilinx <sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">®</sup> XCVU37P FPGAs for the finite-element modeling of a three-phase transformer. The accuracy of the transformer solver under both current-excited and voltage-excited conditions of the transformer was validated against the commercial FE simulation tool.

High-fidelity nonlinear low-order unstructured implicit finite-element seismic simulation of important structures by accelerated element-by-element method

We enable large-scale high-fidelity finite-element seismic response simulations of important structures, that are expected to contribute towards improvement in seismic design verification, by reducing cost of the nonlinear dynamic unstructured low-order implicit finite-element method. Most of the computational cost of this method is involved in the element-by-element (EBE) method, which is a typical example of a “low computation/(data load or store)” kernel that appears in many applications that is not straightforward to attain performance on current computer systems. Therefore, special care based on computer science is required to make use of the potential of computer architecture and achieve fast analysis. In this study, we developed a kernel algorithm and implementation suitable for the target Arm v8.2-A scalable vector extension (SVE) CPU-based supercomputer Fugaku. 5.11- and 8.69-fold speedup was attained by using the developed EBE kernel in a standard preconditioned conjugate gradient solver and a state-of-the-art SC14 massively parallel solver algorithm, respectively. Furthermore, by using the developed EBE kernel in the state-of-the-art solver, a 49 billion degrees-of-freedom high-fidelity seismic response analysis can be conducted in practical speed corresponding to 60,000 time-steps in half a day using Fugaku. The obtained insights are expected to be useful for accelerating other scientific computing methods with “low computation/(data load or store)” kernels.

A GPU-based algorithm for efficient LES of high Reynolds number flows in heterogeneous CPU/GPU supercomputers

Αn optimized MPI+OpenACC implementation model that performs efficiently in CPU/GPU systems using large-eddy simulation is presented. The code was validated for the simulation of wave boundary-layer flows against numerical and experimental data in the literature. A direct Fast-Fourier-Transform-based solver was developed for the solution of the Poisson equation for pressure taking advantage of the periodic boundary conditions. This solver was optimized for parallel execution in CPUs and outperforms by 10 times in computational time a typical iterative preconditioned conjugate gradient solver in GPUs. In terms of parallel performance, an overlapping strategy was developed to reduce the overhead of performing MPI communications using GPUs. As a result, the weak scaling of the algorithm was improved up to 30%. Finally, a large-scale simulation (Re = 2 × 105) using a grid of 4 × 108 cells was executed, and the performance of the code was analyzed. The simulation was launched using up to 512 nodes (512 GPUs + 6144 CPU-cores) on one of the current top 10 supercomputers of the world (Piz Daint). A comparison of the overall computational time showed that the GPU version was 4.2 times faster than the CPU one. The parallel efficiency of this strategy (47%) is competitive compared with the state-of-the-art CPU implementations, and it has the potential to take advantage of modern supercomputing capabilities.

An Out-of-Core Eigen-Solver with OpenMP Parallel Scheme for Large Spare Damped System

An out-of-core block Lanczos method with the OpenMP parallel scheme was developed to solve large spare damped eigenproblems. The symmetric generalized eigenproblem is first solved using the block Lanczos method with the preconditioned conjugate gradient (PCG) method, and the condensed damped eigenproblem is then solved to obtain the complex eigenvalues. Since the PCG solvers and out-of-core schemes are used, a large-scale eigenproblem can be solved using minimal computer memory. The out-of-core arrays only need to be read once in each Lanczos iteration, so the proposed method requires little extra CPU time. In addition, the second-level OpenMP parallel computation in the PCG solver is suggested to avoid using a large block size that often increases the number of iterations needed to achieve convergence.

Efficient basis updating for parametric nonlinear model order reduction

AbstractNonlinear model reduction is used to speed up the numerical solution of finite element models for vibration analysis of structures undergoing large deflections. This speed up is highly desired in design and optimization applications where parametric models are considered and the outcoming high‐dimensional differential equation must be solved multiple times. A first step is to approximate the solution vector i.e. the displacements of the nodes by a linear combination of some basis vectors. One common choice for these basis vectors is a combination of vibration modes and static modal derivatives. However, these vectors depend on parameter values of the parameterized system. This contribution shows how these basis vectors can be updated in an efficient manner. The vibration modes are updated by an inverse free preconditioned Krylov subspace method while the static derivatives are updated by a preconditioned conjugate gradient solver. A case study with a parametric beam gives a first insight into the performance of the proposed method.

Solving linear equations with messenger-field and conjugate gradient techniques: An application to CMB data analysis

We discuss linear system solvers invoking a messenger-field and compare them with (preconditioned) conjugate gradient approaches. We show that the messenger-field techniques correspond to fixed point iterations of an appropriately preconditioned initial system of linear equations. We then argue that a conjugate gradient solver applied to the same preconditioned system, or equivalently a preconditioned conjugate gradient solver using the same preconditioner and applied to the original system, will in general ensure at least a comparable and typically better performance in terms of the number of iterations to convergence and time-to-solution. We illustrate our conclusions with two common examples drawn from the cosmic microwave background (CMB) data analysis: Wiener filtering and map-making. In addition, and contrary to the standard lore in the CMB field, we show that the performance of the preconditioned conjugate gradient solver can depend significantly on the starting vector. This observation seems of particular importance in the cases of map-making of high signal-to-noise ratio sky maps and therefore should be of relevance for the next generation of CMB experiments.

Parallel computing of three-dimensional discontinuous deformation analysis based on OpenMP

The computing efficiency of three-dimensional discontinuous deformation analysis (3D-DDA) needs to be improved for large-scale simulations. Among all the subroutines of 3D-DDA, the equation solver is very time-consuming. To accelerate the equation-solving process, this paper proposes implementing the parallel block Jacobi (BJ) and preconditioned conjugate gradient (PCG) iterative solvers into the original 3D-DDA based on OpenMP. The calculation accuracy and computational efficiency are studied by several numerical examples, demonstrating that the modified 3D-DDA with parallel BJ or PCG solver exhibits much higher execution efficiency with satisfactory correctness. The maximum speedup ratio is up to 5.1 for the cases studied.

Weighted total variation using split Bregman fast quantitative susceptibility mapping reconstruction method**Project supported by the National Natural Science Foundation of China (Grant Nos. 11474236, 81671674, and 11775184) and the Science and Technology Project of Fujian Province, China (Grant No. 2016Y0078).

An ill-posed inverse problem in quantitative susceptibility mapping (QSM) is usually solved using a regularization and optimization solver, which is time consuming considering the three-dimensional volume data. However, in clinical diagnosis, it is necessary to reconstruct a susceptibility map efficiently with an appropriate method. Here, a modified QSM reconstruction method called weighted total variation using split Bregman (WTVSB) is proposed. It reconstructs the susceptibility map with fast computational speed and effective artifact suppression by incorporating noise-suppressed data weighting with split Bregman iteration. The noise-suppressed data weighting is determined using the Laplacian of the calculated local field, which can prevent the noise and errors in field maps from spreading into the susceptibility inversion. The split Bregman iteration accelerates the solution of the L1-regularized reconstruction model by utilizing a preconditioned conjugate gradient solver. In an experiment, the proposed reconstruction method is compared with truncated k-space division (TKD), morphology enabled dipole inversion (MEDI), total variation using the split Bregman (TVSB) method for numerical simulation, phantom and in vivo human brain data evaluated by root mean square error and mean structure similarity. Experimental results demonstrate that our proposed method can achieve better balance between accuracy and efficiency of QSM reconstruction than conventional methods, and thus facilitating clinical applications of QSM.

Adaptive precision in block‐Jacobi preconditioning for iterative sparse linear system solvers

SummaryWe propose an adaptive scheme to reduce communication overhead caused by data movement by selectively storing the diagonal blocks of a block‐Jacobi preconditioner in different precision formats (half, single, or double). This specialized preconditioner can then be combined with any Krylov subspace method for the solution of sparse linear systems to perform all arithmetic in double precision. We assess the effects of the adaptive precision preconditioner on the iteration count and data transfer cost of a preconditioned conjugate gradient solver. A preconditioned conjugate gradient method is, in general, a memory bandwidth‐bound algorithm, and therefore its execution time and energy consumption are largely dominated by the costs of accessing the problem's data in memory. Given this observation, we propose a model that quantifies the time and energy savings of our approach based on the assumption that these two costs depend linearly on the bit length of a floating point number. Furthermore, we use a number of test problems from the SuiteSparse matrix collection to estimate the potential benefits of the adaptive block‐Jacobi preconditioning scheme.