Preconditioned Conjugate Gradient Research Articles

Acceleration of unstructured implicit low-order finite-element earthquake simulation using OpenACC on Pascal GPUs

We accelerate CPU-based unstructured implicit low-order finite-element simulations by porting to a GPU-CPU heterogeneous compute environment using OpenACC. We modified performance-sensitive parts of the code, such as sparse matrix-vector multiplication and MPI communication, so that computations would be suitable for GPUs. Other parts of the earthquake simulation code are ported by directly inserting OpenACC directives into the CPU code. This porting approach enables high performance with relatively low development costs. When comparing eight K computer nodes and eight NVIDIA Pascal P100 GPUs, we achieve 20.8 times speedup for the 3 × 3 block Jacobi preconditioned conjugate gradient finite-element solver. We show the effectiveness of the proposed method through many-case crust-deformation simulations and a large-scale computation using finite element model with 109 degrees-of-freedom on a GPU cluster.

Gravitational gradients derived from GOCE and density structures beneath the North China Craton

In this study, the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) gravitational gradients in the North China Craton (NCC) are processed. The normal gravitational gradient correction, the downward continuation from the mean orbit altitude to the geoid, and the gravitational gradient corrections for the topography, sediment, Moho undulations and the deep-source field are taken into account. After processing, the residual gravitational gradient anomalies corresponding to density heterogeneities are acquired. The preconditioned conjugate gradient method is introduced to invert the density structures using the roughness matrix to constrain the density distribution and the depth-weighting function to prevent the density distribution from migrating to the top of the model due to the decay of the gravitational gradient kernel with depth. The density structure derived from the seismic tomography is used as the initial model, and the GOCE gravitational gradient anomalies are used to invert the density structures in the NCC using the preconditioned conjugate gradient method. The results show that the density structures of the NCC lithosphere reflect obvious lateral and vertical heterogeneities. The density structure in the eastern part of the NCC is different from that in the western part of the NCC. The distribution of the density anomalies corresponds to the tectonic characteristics. The density structures in Taihang Mountain exhibit segmentation features. The thinning of the lithosphere observed in the NCC may extend into the northern part of Taihang Mountain. The density changes significantly within the area between the Ordos Block and the Qilian Block, which corresponds with the different lithospheric features in these two blocks.

An Operator Smoothing with ILU(0) for Aggregation-based Algebraic Multigrid

Due to the potential optimal convergence, algebraic multigrid is widely used in solving large sparse linear systems. The aggregation-based version is one of the widely used methods for its cheap cost and easy implementation. But its convergence is often slow compared to other versions. In this paper, a smoothing technique based on incomplete LU factorization without fill-in is presented to the operators on each level. Each operator is approximately factorized, and the derived lower and upper triangular factors are approximately inverted, which are applied to the operator from both sides to improve the diagonal dominance, and then the effectiveness of the smoothing process and the accuracy of the grid transfer operators. The numerical results show that when incorporated into the preconditioned conjugate gradient iteration, the convergence rate is greatly improved, and though the time used for setup is larger, the time used for iteration and the overall time can also be reduced for large-scale systems.

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.

Adaptive BEM with inexact PCG solver yields almost optimal computational costs

We consider the preconditioned conjugate gradient method (PCG) with optimal preconditioner in the frame of the boundary element method for elliptic first-kind integral equations. Our adaptive algorithm steers the termination of PCG as well as the local mesh-refinement. Besides convergence with optimal algebraic rates, we also prove almost optimal computational complexity. In particular, we provide an additive Schwarz preconditioner which can be computed in linear complexity and which is optimal in the sense that the condition numbers of the preconditioned systems are uniformly bounded. As model problem serves the 2D or 3D Laplace operator and the associated weakly-singular integral equation with energy space widetilde{H}^{-1/2}(Gamma ). The main results also hold for the hyper-singular integral equation with energy space H^{1/2}(Gamma ).

Learning target-aware correlation filters for visual tracking

Discriminative Correlation Filters (DCF) have achieved enormous popularity in the tracking community. Generally, DCF based trackers assume that the target can be well shaped by an axis-aligned bounding box. Therefore, in terms of irregularly shaped objects, the learned correlation filter is unavoidably deteriorated by the background pixels inside the bounding box. To tackle this problem, we propose Target-Aware Correlation Filters (TACF) for visual tracking. A target likelihood map is introduced to impose discriminative weight on filter values according to the probability of this location belonging to the foreground target. According to the TACF formulation, we further propose an optimization strategy based on the Preconditioned Conjugate Gradient method for efficient filter learning. With hand-crafted features (HOG), our approach achieves state-of-the-art performance (62.8% AUC) on OTB100 while running in real-time (24 fps) on a single CPU. With shallow convolutional features, our approach achieves 66.7% AUC on OTB100 and the top rank in EAO on the VOT2016 challenge.

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.

Overlapping Schwarz methods with adaptive coarse spaces for multiscale problems in 3D

We propose two variants of the overlapping additive Schwarz method for the finite element discretization of scalar elliptic problems in 3D with highly heterogeneous coefficients. The methods are efficient and simple to construct using the abstract framework of the additive Schwarz method, and an idea of adaptive coarse spaces. In one variant, the coarse space consists of finite element functions associated with the wire basket nodes and functions based on solving some generalized eigenvalue problems on the faces. In the other variant, it contains functions associated with the vertex nodes with functions based on solving some generalized eigenvalue problems on subdomain faces and subdomain edges. The functions that constitute the coarse spaces are chosen adaptively, and they correspond to the eigenvalues that are smaller than a given threshold. The convergence rate of the preconditioned conjugate gradients method in both cases is shown to be independent of the variations in the coefficients for the sufficient number of eigenfunctions in the coarse space. Numerical results are given to support the theory.

Inversion of the density structure of the lithosphere in the North China Craton from GOCE satellite gravity gradient data

Lithospheric thinning in the North China Craton (NCC), one of the oldest cratons on Earth, has been a focus of geoscientists across the globe for a long time. In this study, seismic tomographic P-wave velocity variations are converted to density perturbations as the initial constraint values, and the four components of the GOCE Level-2 gravity gradient product are subjected to topographic effect correction, relief correction on the Moho and sedimentary layer boundary, and long-wavelength correction. In addition, the anomalous gravity gradient effect resulting from the uncertainty of the Moho and sedimentary layer depth is considered. Considering the effects of temperature variations in the lithosphere, which leads to the non-uniform density distribution, the anomalous gravity gradient effect resulting from the temperature variations is corrected at the first time in this area. Inverse computation is performed based on the preconditioned conjugate gradient algorithm. Because the Lagrange empirical parameter in the algorithm exhibits uncertainty in data inversion, the regularization parameter, which is typically an empirical parameter, is replaced with the value of the inflection point of the L-curve. The results show the followings. The non-uniform density distribution in the lithosphere is related not only to the composition but also to the internal temperature variations in the lithosphere. The density of the lithosphere in the NCC is significantly non-uniform in both the horizontal and vertical directions and exhibits a notable segment-wise spatial distribution pattern. Two tectonic units, namely the Taihangshan tectonic zone and the Linfen–Weihe graben, constitute a central gravity gradient transition zone. There is a significant difference in the density distribution on the two sides of this central transition zone.

Deflated preconditioned conjugate gradient method for solving single-step BLUP models efficiently

BackgroundThe single-step single nucleotide polymorphism best linear unbiased prediction (ssSNPBLUP) method, such as single-step genomic BLUP (ssGBLUP), simultaneously analyses phenotypic, pedigree, and genomic information of genotyped and non-genotyped animals. In contrast to ssGBLUP, SNP effects are fitted explicitly as random effects in the ssSNPBLUP model. Similarly, principal components associated with the genomic information can be fitted explicitly as random effects in a single-step principal component BLUP (ssPCBLUP) model to remove noise in genomic information. Single-step genomic BLUP is solved efficiently by using the preconditioned conjugate gradient (PCG) method. Unfortunately, convergence issues have been reported when solving ssSNPBLUP by using PCG. Poor convergence may be linked with poor spectral condition numbers of the preconditioned coefficient matrices of ssSNPBLUP. These condition numbers, and thus convergence, could be improved through the deflated PCG (DPCG) method, which is a two-level PCG method for ill-conditioned linear systems. Therefore, the first aim of this study was to compare the properties of the preconditioned coefficient matrices of ssGBLUP and ssSNPBLUP, and to document convergence patterns that are obtained with the PCG method. The second aim was to implement and test the efficiency of a DPCG method for solving ssSNPBLUP and ssPCBLUP.ResultsFor two dairy cattle datasets, the smallest eigenvalues obtained for ssSNPBLUP (ssPCBLUP) and ssGBLUP, both solved with the PCG method, were similar. However, the largest eigenvalues obtained for ssSNPBLUP and ssPCBLUP were larger than those for ssGBLUP, which resulted in larger condition numbers and in slow convergence for both systems solved by the PCG method. Different implementations of the DPCG method led to smaller condition numbers, and faster convergence for ssSNPBLUP and for ssPCBLUP, by deflating the largest unfavourable eigenvalues.ConclusionsPoor convergence of ssSNPBLUP and ssPCBLUP when solved by the PCG method are related to larger eigenvalues and larger condition numbers in comparison to ssGBLUP. These convergence issues were solved with a DPCG method that annihilates the effect of the largest unfavourable eigenvalues of the preconditioned coefficient matrix of ssSNPBLUP and of ssPCBLUP on the convergence of the PCG method. It resulted in a convergence pattern, at least, similar to that of ssGBLUP.

A Summary and Analysis of Low Complexity Methods for Massive MIMO Systems

The massive multiple-input multiple-output (MIMO) systems help increase the wireless channel capacity and spectral efficiency manifold, and the coverage as well as power consumption. In MIMO systems, the linear minimum mean-square error (MMSE) can achieve the near-performance. However, along with the huge number of the BS-to-user antennas at the BS station, the computational complexity becomes a serious problem due to computing matrix inversion. Many methods were proposed to improve complexity, for instance, Neumann series expansion method, Jacobi methods and conjugate gradient (CG) method, but it is still costly. And recently many other methods on account of the above methods were proposed, respectively. Such as Gauss-Seidel method, Symmetric Successive Over Relaxation (SSOR) method, further improved Jacobi method based on soft-output massive MIMO detection scheme, pre-conditioner conjugate gradient (PCG) based on SSOR, split pre-conditioned conjugate gradient (SPCG) method, etc, which all show the reduced complexity, improved convergence and BER performance.

Communication Optimization for Multiphase Flow Solver in the Library of OpenFOAM

Multiphase flow solvers are widely-used applications in OpenFOAM, whose scalability suffers from the costly communication overhead. Therefore, we establish communication-optimized multiphase flow solvers in OpenFOAM. In this paper, we first deliver a scalability bottleneck test on the typical multiphase flow case damBreak and reveal that the Message Passing Interface (MPI) communication in a Multidimensional Universal Limiter for Explicit Solution (MULES) and a Preconditioned Conjugate Gradient (PCG) algorithm is the short slab of multiphase flow solvers. Furthermore, an analysis of the communication behavior is carried out. We find that the redundant communication in MULES and the global synchronization in PCG are the performance limiting factors. Based on the analysis, we propose our communication optimization algorithm. For MULES, we remove the redundant communication and obtain optMULES. For PCG, we import several intermediate variables and rearrange PCG to reduce the global communication. We also overlap the computation of matrix-vector multiply and vector update with the non-blocking computation. The resulting algorithms are respectively referred to as OFPiPePCG and OFRePiPePCG. Extensive experiments show that our proposed method could dramatically increase the parallel scalability and solving speed of multiphase flow solvers in OpenFOAM approximately without the loss of accuracy.

An Improved Domain Decomposition Method for Drop Impact Reliability Analysis of 3-D ICs

Drop test is usually adopted in the integrated circuit (IC) testing to estimate the shock resistance capability of IC packaging. Generally, it is very time-consuming for the drop test simulation, and therefore the fast numerical approach is needed to reduce the computational cost. In this paper, we propose an efficient drop test simulator for through-silicon-via (TSV)-based 3-D IC to simulate its mechanical behaviors under drop impact. The proposed simulator is based on the idea of domain decomposition (DD) to improve the condition number of the coefficient matrix of the solver. We further develop two efficient numerical techniques in the simulator to improve its efficiency. First, a second-order formulation is proposed for the initial solution selection in preconditioned conjugate gradient (PCG) solver, which can efficiently reduce the number of PCG iterations at each time point during simulation. Second, an equivalent structure-embedded model is proposed and applied in the first Schwarz iteration to efficiently reduce the number of Schwarz iterations in DD method. Numerical experiments show that the proposed drop test simulator can achieve $7.03\times $ speedup in comparison with the conventional finite-element method-based solver. It is demonstrated in this paper through several examples with multichip layers, of which each chip layer consists of an $8\times 8$ TSV array, that the proposed simulator can be widely applied to reliability analysis of 3-D ICs under drop impact.

Iterative map-making with two-level preconditioning for polarized cosmic microwave background data sets

Context. An estimation of the sky signal from streams of time ordered data (TOD) acquired by the cosmic microwave background (CMB) experiments is one of the most important steps in the context of CMB data analysis referred to as the map-making problem. The continuously growing CMB data sets render the CMB map-making problem progressively more challenging in terms of computational cost and memory in particular in the context of ground-based experiments with their operational limitations as well as the presence of contaminants. Aims. We study a recently proposed, novel class of the Preconditioned Conjugate Gradient (PCG) solvers which invoke two-level preconditioners in the context of the ground-based CMB experiments. We compare them against the PCG solvers commonly used in the map-making context considering their precision and time-to-solution. Methods. We compare these new methods on realistic, simulated data sets reflecting the characteristics of current and forthcoming CMB ground-based experiments. We develop a divide-and-conquer implementation of the approach where each processor performs a sequential map-making for a subset of the TOD. Results. We find that considering the map level residuals, the new class of solvers permits us to achieve a tolerance that is better than the standard approach by up to three orders of magnitude, where the residual level often saturates before convergence is reached. This often corresponds to an important improvement in the precision of the recovered power spectra in particular on the largest angular scales. The new method also typically requires fewer iterations to reach a required precision and therefore shorter run times are required for a single map-making solution. However, the construction of an appropriate two-level preconditioner can be as costly as a single standard map-making run. Nevertheless, if the same problem needs to be solved multiple times, for example, as in Monte Carlo simulations, this cost is incurred only once, and the method should be competitive, not only as far as its precision is concerned but also its performance.

Simulation of Fluid-Thermal Field in Oil-Immersed Transformer Winding Based on Dimensionless Least-Squares and Upwind Finite Element Method

In order to study the coupling fluid and thermal problems of the local winding in oil-immersed power transformers, the least-squares finite element method (LSFEM) and upwind finite element method (UFEM) are adopted, respectively, to calculate the fluid and thermal field in the oil duct. When solving the coupling problem by sequential iterations, the effect of temperature on the material property and the loss density of the windings should be taken into account. In order to improve the computation efficiency for the coupling fields, an algorithm, which adopts two techniques, the dimensionless LSFEM and the combination of Jacobi preconditioned conjugate gradient method (JPCGM) and the two-side equilibration method (TSEM), is proposed in this paper. To validate the efficiency of the proposed algorithm, a local winding model of a transformer is built and the fluid field is computed by the conventional LSFEM, dimensionless LSFEM, and the Fluent software. While the fluid and thermal computation results of the local winding model of a transformer obtained by the two LSFEMs are basically consistent with those of the Fluent software, the stiffness matrix, which is formed by the dimensionless scheme of LSFEM and preconditioned by the JPCGM and TSEM, has a smaller condition number and a faster convergence rate of the equations. Thus, it demonstrates a broader applicability.

BestSF

The Sparse Matrix-Vector Multiplication (SpMV) kernel dominates the computing cost in numerous scientific applications. Many implementations based on different sparse formats were proposed to improve this kernel on the recent GPU architectures. However, it has been widely observed that there is no “best-for-all” sparse format for the SpMV kernel on GPU. Indeed, serious performance degradation of an order of magnitude can be observed without a careful selection of the sparse format to use. To address this problem, we propose in this article BestSF (Best Sparse Format), a new learning-based sparse meta-format that automatically selects the most appropriate sparse format for a given input matrix. To do so, BestSF relies on a cost-sensitive classification system trained using Weighted Support Vector Machines (WSVMs) to predict the best sparse format for each input sparse matrix. Our experimental results on two different NVIDIA GPU architectures using a large number of real-world sparse matrices show that BestSF achieved a noticeable overall performance improvement over using a single sparse format. While BestSF is trained to select the best sparse format in terms of performance (GFLOPS), our further experimental investigations revealed that using BestSF also led, in most of the test cases, to the best energy efficiency (MFLOPS/W). To prove its practical effectiveness, we also evaluate the performance and energy efficiency improvement achieved when using BestSF as a building block in a GPU-based Preconditioned Conjugate Gradient (PCG) iterative solver.

A Class of Symmetric Factored Approximate Inverses and Hybrid Two-Level Solver

A new class of symmetric factored approximate inverses is proposed and used in conjunction with the Preconditioned Conjugate Gradient method for solving sparse symmetric linear systems. Additionally, a new hybrid two-level solver is proposed utilizing a block independent set reordering, in order to create the two level hierarchy. The Schur complement is formed explicitly by inverting the blocks created by reordering. Then, the preconditioned conjugate gradient method is used in conjunction with the symmetric factored approximate inverse to solve the reduced order linear system. Furthermore, numerical results on the performance and convergence behavior for solving various model problems are presented.

PIBLUP: High-Performance Software for Large-Scale Genetic Evaluation of Animals and Plants.

Today, the rapid increase in phenotypic and genotypic information is leading to larger mixed model equations (MMEs) and rendering genetic evaluation more time-consuming. It has been demonstrated that a preconditioned conjugate gradient (PCG) algorithm via an iteration on data (IOD) technique is the most efficient method of solving MME at a low computing cost. Commonly used software applications implementing PCG by IOD merely employ functions from the Intel Math Kernel Library (MKL) to accelerate numerical computations and have not taken full advantage of the multicores or multiprocessors of computer systems to reduce the execution time. Making the most of multicore/multiprocessor systems, we propose PIBLUP, a parallel, shared memory implementation of PCG by IOD to minimize the execution time of genetic evaluation. In addition to functions in MKL, PIBLUP uses Message Passing Interface (MPI) shared memory programming to parallelize code in the entire workflow where possible. Results from the analysis of the two datasets show that the execution time was reduced by more than 80% when solving MME using PIBLUP with 16 processes in parallel, compared to a serial program using a single process. PIBLUP is a high-performance tool for users to efficiently perform genetic evaluation. PIBLUP with its user manual is available at https://github.com/huiminkang/PIBLUP.

A Folding Calculation Method Based on the Preconditioned Conjugate Gradient Inversion Algorithm of Gravity Gradient Tensor

To solve the non-uniqueness problem of gravity gradient inversion, we proposed a folding calculation method based on preconditioned conjugate gradient inversion. Compared with the original algorithm in which the entire study area is taken as the research subject and all grids are used simultaneously in the inversion, the proposed folding method divides the research area into several sub-areas. A prism unit from any of the four corner grids is selected for the first iteration, whose density anomaly result is taken as the initial density anomaly for the next iteration of the same sub-area. The folding in the left–right and up–down directions takes turns during the calculation until the inversion calculation has covered the entire research area. This folding algorithm demonstrates strong regularity. The inversion results of multiple synthetic models show that the folding calculation method performs multiple parameter corrections in the initial model. Meanwhile, on the basis that the relative fitting of standard errors of the observed values satisfies the convergence condition, the model errors are constrained, and thus, the model errors of the inversion results are consequently reduced. The 3-D density anomaly pattern over the Craton area in North China was obtained via a joint inversion of the four components ( $$\varvec{T}_{\text{xx}} ,$$ $$\varvec{T}_{\text{xz}} ,$$ $$\varvec{T}_{\text{yy}} ,$$ and $$\varvec{T}_{\text{zz}}$$ ) of the Gravity field and steady-state Ocean Circulation Explorer L2 gravity gradient after preprocessing. We compared the inversion results from the folding calculation method and the original method, and performed detailed analysis and discussion on the inversion results with existing geological and geophysical data. Our analysis shows that the improved calculation method is effectively applicable to the inversion of measured gravity gradient data, and the inversion results provide more detailed and reliable pattern information for the density anomaly.