Conjugate Gradient Solver Research Articles

An adapted deflated conjugate gradient solver for robust extended/generalised finite element solutions of large scale, 3D crack propagation problems

An adapted deflation preconditioner is employed to accelerate the solution of linear systems resulting from the discretisation of fracture mechanics problems with well-conditioned extended/generalised finite elements. The deflation space typically used for linear elasticity problems is enriched with additional vectors, accounting for the enrichment functions used, thus effectively removing low frequency components of the error. To further improve performance, deflation is combined, in a multiplicative way, with a block-Jacobi preconditioner, which removes high frequency components of the error as well as near-linear dependencies introduced by enrichment. The resulting scheme is tested on a series of non-planar crack propagation problems and compared to alternative linear solvers in terms of performance.

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.

Scalable and accurate multi-GPU-based image reconstruction of large-scale ptychography data

While the advances in synchrotron light sources, together with the development of focusing optics and detectors, allow nanoscale ptychographic imaging of materials and biological specimens, the corresponding experiments can yield terabyte-scale volumes of data that can impose a heavy burden on the computing platform. Although graphics processing units (GPUs) provide high performance for such large-scale ptychography datasets, a single GPU is typically insufficient for analysis and reconstruction. Several works have considered leveraging multiple GPUs to accelerate the ptychographic reconstruction. However, most of these works utilize only the Message Passing Interface to handle the communications between GPUs. This approach poses inefficiency for a hardware configuration that has multiple GPUs in a single node, especially while reconstructing a single large projection, since it provides no optimizations to handle the heterogeneous GPU interconnections containing both low-speed (e.g., PCIe) and high-speed links (e.g., NVLink). In this paper, we provide an optimized intranode multi-GPU implementation that can efficiently solve large-scale ptychographic reconstruction problems. We focus on the maximum likelihood reconstruction problem using a conjugate gradient (CG) method for the solution and propose a novel hybrid parallelization model to address the performance bottlenecks in the CG solver. Accordingly, we have developed a tool, called PtyGer (Ptychographic GPU(multiple)-based reconstruction), implementing our hybrid parallelization model design. A comprehensive evaluation verifies that PtyGer can fully preserve the original algorithm’s accuracy while achieving outstanding intranode GPU scalability.

Efficient distributed approach for density-based topology optimization using coarsening and h-refinement

This work presents an efficient parallel implementation of density-based topology optimization using Adaptive Mesh Refinement (AMR) schemes to reduce the computational burden of the bottleneck of the process, the evaluation of the objective function using Finite Element Analysis (FEA). The objective is to obtain an equivalent design to the one generated on a uniformly fine mesh using distributed memory computing but at a much cheaper computational cost. We propose using a fine mesh for the optimization and a coarse mesh for the analysis using coarsening and refinement criteria based on the thresholding of design variables. We evaluate the functional on the coarse mesh using a distributed conjugate gradient solver preconditioned by an algebraic multigrid (AMG) method showing its computational advantages in some cases by comparing with geometric multigrid (GMG) and AMG methods in two- and three-dimensional problems. We use different computational resources with small regularization distances for such comparisons. We also evaluate the performance and scalability of the proposal using a different number of computing cores and distributed computing hosts. The numerical results show a significant increment of the computing performance for the overall computing time of the proposal combining dynamic coarsening, adaptive mesh refinement, and distributed memory computing architectures.

Learning in high-dimensional feature spaces using ANOVA-based fast matrix-vector multiplication

<p style='text-indent:20px;'>Kernel matrices are crucial in many learning tasks such as support vector machines or kernel ridge regression. The kernel matrix is typically dense and large-scale. Depending on the dimension of the feature space even the computation of all of its entries in reasonable time becomes a challenging task. For such dense matrices the cost of a matrix-vector product scales quadratically with the dimensionality <inline-formula><tex-math id="M1">\begin{document}$ N $\end{document}</tex-math></inline-formula>, if no customized methods are applied. We propose the use of an ANOVA kernel, where we construct several kernels based on lower-dimensional feature spaces for which we provide fast algorithms realizing the matrix-vector products. We employ the non-equispaced fast Fourier transform (NFFT), which is of linear complexity for fixed accuracy. Based on a feature grouping approach, we then show how the fast matrix-vector products can be embedded into a learning method choosing kernel ridge regression and the conjugate gradient solver. We illustrate the performance of our approach on several data sets.</p>

Acceleration strategies for Tridimensional Coupled hydromechanical problems based on CPU and GPU programming in MATLAB.

Large-scale fluid flow in porous media demands intense computations and occurs in the most diverse applications, including groundwater flow and oil recovery. This article presents novel computational strategies applied to reservoir geomechanics. Advances are proposed for the efficient assembly of finite element matrices and the solution of linear systems using highly vectorized code in MATLAB. In the CPU version, element matrix assembly is performed using conventional vectorization procedures, based on two strategies: the explicit matrices, and the multidimensional products. Further assembly of the global sparse matrix is achieved using the native sparse function. For the GPU version, computation of the complete set of element matrices is performed with the same strategies as the CPU approach, using gpuArray structures and the native CUDA support provided by MATLAB Parallel Computing Toolbox. Solution of the resulting linear system in CPU and GPU versions is obtained with two strategies using a one-way approach: the native conjugate gradient solver (pcg), and the one provided by the Eigen library. A broad discussion is presented in a dedicated benchmark, where the different strategies using CPU and GPU are compared in processing time and memory requirements. These analyses present significant speedups over serial codes.

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.

An Inexact Optimal Hybrid Conjugate Gradient Method for Solving Symmetric Nonlinear Equations

This article presents an inexact optimal hybrid conjugate gradient (CG) method for solving symmetric nonlinear systems. The method is a convex combination of the optimal Dai–Liao (DL) and the extended three-term Polak–Ribiére–Polyak (PRP) CG methods. However, two different formulas for selecting the convex parameter are derived by using the conjugacy condition and also by combining the proposed direction with the default Newton direction. The proposed method is again derivative-free, therefore the Jacobian information is not required throughout the iteration process. Furthermore, the global convergence of the proposed method is shown using some appropriate assumptions. Finally, the numerical performance of the method is demonstrated by solving some examples of symmetric nonlinear problems and comparing them with some existing symmetric nonlinear equations CG solvers.

Multiparameter optimization of nonuniform passive diffusion properties for creating coarse-grained equivalent models of cardiac propagation

The arrhythmogenic role of discrete cardiac propagation may be assessed by comparing discrete (fine-grained) and equivalent continuous (coarse-grained) models. We aim to develop an optimization algorithm for estimating the smooth conductivity field that best reproduces the diffusion properties of a given discrete model. Our algorithm iteratively adjusts local conductivity of the coarse-grained continuous model by simulating passive diffusion from white noise initial conditions during 3–10 ms and computing the root mean square error with respect to the discrete model. The coarse-grained conductivity field was interpolated from up to 300 evenly spaced control points. We derived an approximate formula for the gradient of the cost function that required (in two dimensions) only two additional simulations per iteration regardless of the number of estimated parameters. Conjugate gradient solver facilitated simultaneous optimization of multiple conductivity parameters. The method was tested in rectangular anisotropic tissues with uniform and nonuniform conductivity (slow regions with sinusoidal profile) and random diffuse fibrosis, as well as in a monolayer interconnected cable model of the left atrium with spatially-varying fibrosis density. Comparison of activation maps served as validation. The results showed that after convergence the errors in activation time were < 1 ms for rectangular geometries and 1–3 ms in the atrial model. Our approach based on the comparison of passive properties (<10 ms simulation) avoids performing active propagation simulations (>100 ms) at each iteration while reproducing activation maps, with possible applications to investigating the impact of microstructure on arrhythmias.

Variable-Metric Localization of Occupied and Virtual Orbitals.

The key idea of the variable-metric approach to orbital localization is to allow nonorthogonality between orbitals while, at the same time, preventing them from becoming linearly dependent. The variable-metric localization has been shown to improve the locality of occupied nonorthogonal orbitals relative to their orthogonal counterparts. In this work, numerous localization algorithms are designed and tested to exploit the conceptual simplicity of the variable-metric approach with the goal of creating a straightforward and reliable localization procedure for virtual orbitals. The implemented algorithms include the steepest descent, conjugate gradient (CG), limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS), and hybrid procedures as well as trust-region (TR) methods based on the CG and Cauchy-point subproblem solvers. Comparative analysis shows that the CG-based TR algorithm is the best overall method to obtain nonorthogonal localized molecular orbitals (NLMOs), occupied or virtual. The L-BFGS and CG algorithms can also be used to obtain NLMOs reliably but often at higher computational cost. Extensive tests demonstrate that the implemented methods allow us to obtain well-localized Boys-Foster (i.e., maximally localized Wannier functions) and Pipek-Mezey, orthogonal and nonorthogonal, and occupied and virtual orbitals for a variety of gas-phase molecules and periodic materials. The tests also show that virtual NLMOs, which have not been described before, are, on average, 13% (Boys-Foster) and 18% (Pipek-Mezey) more localized than their orthogonal counterparts.

Multiscale Cholesky preconditioning for ill-conditioned problems

Many computer graphics applications boil down to solving sparse systems of linear equations. While the current arsenal of numerical solvers available in various specialized libraries and for different computer architectures often allow efficient and scalable solutions to image processing, modeling and simulation applications, an increasing number of graphics problems face large-scale and ill-conditioned sparse linear systems --- a numerical challenge which typically chokes both direct factorizations (due to high memory requirements) and iterative solvers (because of slow convergence). We propose a novel approach to the efficient preconditioning of such problems which often emerge from the discretization over unstructured meshes of partial differential equations with heterogeneous and anisotropic coefficients. Our numerical approach consists in simply performing a fine-to-coarse ordering and a multiscale sparsity pattern of the degrees of freedom, using which we apply an incomplete Cholesky factorization. By further leveraging supernodes for cache coherence, graph coloring to improve parallelism and partial diagonal shifting to remedy negative pivots, we obtain a preconditioner which, combined with a conjugate gradient solver, far exceeds the performance of existing carefully-engineered libraries for graphics problems involving bad mesh elements and/or high contrast of coefficients. We also back the core concepts behind our simple solver with theoretical foundations linking the recent method of operator-adapted wavelets used in numerical homogenization to the traditional Cholesky factorization of a matrix, providing us with a clear bridge between incomplete Cholesky factorization and multiscale analysis that we leverage numerically.

Efficient second-order ADI difference schemes for three-dimensional Riesz space-fractional diffusion equations

In this paper, a three-dimensional time-dependent Riesz space-fractional diffusion equation is considered, and an alternating direction implicit (ADI) difference scheme is proposed, in which a weighted and shifted Grünwald difference scheme (see Tian et al. (2015) [33]) is utilized for the discretizations of space-fractional derivatives, and a fractional-Douglas-Gunn type ADI method is utilized for the discretization of time derivative. The method is proved to be unconditionally stable and convergent with second-order accuracy both in time and space with respect to a weighted discrete energy norm. Efficient implementation of the method is carefully discussed, and then based on fast matrix-vector multiplications, a fast conjugate gradient (FCG) solver for the resulting symmetric positive definite linear algebraic system is developed. Numerical experiments support the theoretical analysis and show strong effectiveness and efficiency of the method for large-scale modeling and simulations. Finally, a linearized ADI scheme based on second-order extrapolation method is developed and tested for the nonlinear Riesz space-fractional diffusion equation.

Solution of large-scale 3D controlled-source electromagnetic modeling problem using efficient iterative solvers

With geophysical surveys evolving from traditional 2D to 3D models, the large volume of data adds challenges to inversion, especially when aiming to resolve complex 3D structures. An iterative forward solver for the controlled-source electromagnetic (CSEM) method requires less memory than that for a direct solver; however, it is not easy to iteratively solve an ill-conditioned linear system of equations arising from finite-element discretization of Maxwell’s equations. To solve this problem, we have developed efficient and robust iterative solvers for frequency- and time-domain CSEM modeling problems. For the frequency-domain problem, we first transform the linear system into its equivalent real-number format, and then introduce an optimal block-diagonal preconditioner. Because the condition number of the preconditioned linear equation system has an upper bound of [Formula: see text], we can achieve fast solution convergence when applying a flexible generalized minimum residual solver. Applying the block preconditioner further results in solving two smaller linear systems with the same coefficient matrix. For the time-domain problem, we first discretize the partial differential equation for the electric field in time using an unconditionally stable backward Euler scheme. We then solve the resulting linear equation system iteratively at each time step. After the spatial discretization in the frequency domain, or space-time discretization in the time domain, we exploit the conjugate-gradient solver with auxiliary-space preconditioners derived from the Hiptmair-Xu decomposition to solve these real linear systems. Finally, we check the efficiency and effectiveness of our iterative methods by simulating complex CSEM models. The most significant advantage of our approach is that the iterative solvers we adopt have almost the same accuracy and robustness as direct solvers but require much less memory, rendering them more suitable for large-scale 3D CSEM forward modeling and inversion.

Multi-GPU acceleration of large-scale density-based topology optimization

This work presents a parallel implementation of density-based topology optimization using distributed GPU computing systems. The use of multiple GPU devices allows us accelerating the computing process and increasing the device memory available for GPU computing. This increment of device memory enables us to address large models that commonly do not fit into one GPU device. The most modern scientific computers incorporate these devices to design energy-efficient, low-cost, and high-computing power systems. However, we should adopt the proper techniques to take advantage of the computational resources of such high-performance many-core computing systems. It is well-known that the bottleneck of density-based topology optimization is the solving of the linear elasticity problem using Finite Element Analysis (FEA) during the topology optimization iterations. We solve the linear system of equations obtained from FEA using a distributed conjugate gradient solver preconditioned by a smooth aggregation-based algebraic multigrid (AMG) using GPU computing with multiple devices. The use of aggregation-based AMG reduces memory requirements and improves the efficiency of the interpolation operation. This fact is rewarding for GPU computing. We evaluate the performance and scalability of the distributed GPU system using structured and unstructured meshes. We also test the performance using different 3D finite elements and relaxing operators. Besides, we evaluate the use of numerical approaches to increase the topology optimization performance. Finally, we present a comparison between the many-core computing instance and one efficient multi-core implementation to highlight the advantages of using GPU computing in large-scale density-based topology optimization problems.

Migration deconvolution using domain decomposition

Least-squares migration (LSM) is an effective technique for mitigating blurring effects and migration artifacts generated by limited data frequency bandwidth, incomplete coverage of geometry, source signature, and unbalanced amplitudes caused by complex wavefield propagation in the subsurface. Migration deconvolution (MD) is an image-domain approach for LSM that approximates the Hessian operator using a set of precomputed point spread functions. We have developed a new workflow by integrating the MD and domain decomposition (DD) methods. DD techniques aim to solve large and complex linear systems by splitting problems into smaller parts, facilitating parallel computing, and providing a higher convergence in iterative algorithms. We suggest that instead of solving the problem in a unique domain, as conventionally performed, splitting the problem into subdomains that overlap and solve each of them independently. We accelerate the convergence rate of the conjugate-gradient solver by applying the DD methods to retrieve better reflectivity, which is mainly visible in regions with low amplitudes. Moreover, using the pseudo-Hessian operator, the convergence of the algorithm is accelerated, suggesting that the inverse problem becomes better conditioned. Experiments using the synthetic Pluto model demonstrate that our algorithm dramatically reduces the required number of iterations while providing a considerable enhancement in image resolution and better continuity of poorly illuminated events.

Iterative Approximation of Preconditioning Matrices Through Krylov-Type Solver Iterations

Linear equation solvers require considerable computational time in many computer simulation methods, such as the stress analysis using the finite element method (FEM), or the fluid dynamics using an implicit time integration scheme. Because of their favorable nature to modern supercomputers, Krylov-type linear equation solvers have become dominant. Krylov-type solvers are usually used with preconditioners, which accelerate the convergence of Krylov solvers, and therefore developing robust preconditioners draws attention from researchers. The basic idea of preconditioning is to transform the original system to a system which can be solved more easily by using a matrix which resembles the original system but the associated linear system can be solved without any difficulties. In this study, we propose a new methodology of constructing such a matrix by updating the matrix using the information obtained through the Krylov solver computation. More precisely, we employ the limited memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) Hessian matrix update scheme. The proposed scheme is combined with the orthogonal minimum (ORTHOMIN) residual method, which accepts a variable preconditioner. As a result of performance tests, in the best case scenario, (1) with the Matrix Market matrices, our new method combined with Algebraic Multigrid (AMG) shows the successful convergence in 12 cases out of 13 problems, whilst the conventional AMG converged only in four cases and (2) with FEM problems, we obtained [Formula: see text] acceleration in the rate of converge and [Formula: see text] in computational time over the AMG preconditioned conjugate gradient (CG) solver. With those results, our newly developed algorithm provides robustness to the conventional preconditioning, whilst the computational speed is superior to conventional ones.

Replicated Computational Results (RCR) Report for “Adaptive Precision Block-Jacobi for High Performance Preconditioning in the Ginkgo Linear Algebra Software”

The article by Flegar et al. titled “Adaptive Precision Block-Jacobi for High Performance Preconditioning in the Ginkgo Linear Algebra Software” presents a novel, practical implementation of an adaptive precision block-Jacobi preconditioner. Performance results using state-of-the-art GPU architectures for the block-Jacobi preconditioner generation and application demonstrate the practical usability of the method, compared to a traditional full-precision block-Jacobi preconditioner. A production-ready implementation is provided in the Ginkgo numerical linear algebra library. In this report, the Ginkgo library is reinstalled and performance results are generated to perform a comparison to the original results when using Ginkgo’s Conjugate Gradient solver with either the full or the adaptive precision block-Jacobi preconditioner for a suite of test problems on an NVIDIA GPU accelerator. After completing this process, the published results are deemed reproducible.

Low occupancy high performance elemental products in assembly free FEM on GPU

Assembly free FEM bypasses the assembly step and solves the system of linear equations at the element level using Conjugate Gradient (CG) type iterative solver. The smaller dense Matrix-vector Products (MvPs) are encapsulated within the CG solver and are computed either at element level or degree of freedom (DoF) level. Both these strategies exploit the computing power of GPU effectively, but the performance is lagging due to the uncoalesced global memory access on GPU. This paper proposes an improved MvP strategy in assembly free FEM, which improves the performance by coalesced global memory access using on-chip faster shared memory and using the texture cache memory on GPU. Since GPU has limited shared memory (in few KBs), the proposed technique suffers from a problem known as low occupancy. Despite the low occupancy issue, the proposed strategy outperforms both element based and DoF based MvP strategies on GPU. Numerical experiments compared with element level and DoF level strategies on GPU and found that, GPU instance of proposed MvP outperforms both strategies approximately by factor of 7 and 1.5 respectively.