Biconjugate Gradient Research Articles

SpMV and BiCG-Stab optimization for a class of hepta-diagonal-sparse matrices on GPU

The abundant data parallelism available in many-core GPUs has been a key interest to improve accuracy in scientific and engineering simulation. In many cases, most of the simulation time is spent in linear solver involving sparse matrix–vector multiply. In forward petroleum oil and gas reservoir simulation, the application of a stencil relationship to structured grid leads to a family of generalized hepta-diagonal solver matrices with some regularity and structural uniqueness. We present a customized storage scheme that takes advantage of generalized hepta-diagonal sparsity pattern and stencil regularity by optimizing both storage and matrix–vector computation. We also present an in-kernel optimization for implementing sparse matrix–vector multiply (SpMV) and biconjugate gradient stabilized (BiCG-Stab) solver. In-kernel is intended to avoid the multiple kernels invocation associated with the use of the numerical library operators. To keep in-kernel, a lock-free inter-block synchronization is used in which completing thread blocks are assigned some independent computations to avoid repeatedly polling the global memory. Other optimizations enable combining reductions and collective write operations to memory. The in-kernel optimization is particularly useful for the iterative structure of BiCG-Stab for preserving vector data locality and to avoid saving vector data back to memory and reloading on each kernel exit and re-entry. Evaluation uses generalized hepta-diagonal matrices that derives from a range of forward reservoir simulation’s structured grids. Results show the profitability of proposed generalized hepta-diagonal custom storage scheme over standard library storage like compressed sparse row, hybrid sparse, and diagonal formats. Using proposed optimizations, SpMV and BiCG-Stab have been noticeably accelerated compared to other implementations using multiple kernel exit–re-entry when the solver is implemented by invoking numerical library operators.

A Novel and Fast SimRank Algorithm

SimRank is a widely adopted similarity measure for objects modeled as nodes in a graph, based on the intuition that two objects are similar if they are referenced by similar objects. The recursive nature of SimRank definition makes it expensive to compute the similarity score even for a single pair of nodes. This defect limits the applications of SimRank. To speed up the computation, some existing works replace the original model with an approximate model to seek only rough solution of SimRank scores. In this work, we propose a novel solution for computing all-pair SimRank scores. In particular, we propose to convert SimRank to the problem of solving a linear system in matrix form, and further prove that the system is non-singular, diagonally dominate, and symmetric definite positive (for undirected graphs). Those features immediately lead to the adoption of Conjugate Gradient (CG) and Bi-Conjugate Gradient (BiCG) techniques for efficiently computing SimRank scores. As a result, a significant improvement on the convergence rate can be achieved; meanwhile, the sparsity of the adjacency matrix is not damaged all the time. Inspired by the existing common neighbor sharing strategy, we further reduce the computational complexity of the matrix multiplication and resolve the scalable issues. The experimental results show our proposed algorithms significantly outperform the state-of-the-art algorithms.

Fast FOCUSS method based on bi-conjugate gradient and its application to space-time clutter spectrum estimation

The focal underdetermined system solver (FOCUSS) is a powerful tool for sparse representation in complex underdetermined systems. This paper presents the fast FOCUSS method based on the bi-conjugate gradient (BICG), termed BICG-FOCUSS, to speed up the convergence rate of the original FOCUSS. BICGFOCUSS was specifically designed to reduce the computational complexity of FOCUSS by solving a complex linear equation using the BICG method according to the rank of the weight matrix in FOCUSS. Experimental results show that BICG-FOCUSS is more efficient in terms of computational time than FOCUSS without losing accuracy. Since FOCUSS is an efficient tool for estimating the space-time clutter spectrum in sparse recoverybased space-time adaptive processing (SR-STAP), we propose BICG-FOCUSS to achieve a fast estimation of the space-time clutter spectrum in mono-static array radar and in the mountaintop system. The high performance of the proposed BICG-FOCUSS in the application is demonstrated with both simulated and real data.

Through-Casing Hydraulic Fracture Evaluation by Induction Logging II: The Inversion Algorithm and Experimental Validations

Evaluation of hydraulic fractures has been under intensive study since the last decade. Among published works, only a few have included casing effects. This paper focuses on the through-casing electromagnetic (EM) induction imaging method that evaluates hydraulic fractures with enhanced contrasts. An experimental system and an inverse scattering algorithm are presented. The laboratory scaled experimental system is built for the feasibility study of the EM induction imaging method for the through-casing hydraulic fracture evaluation. To develop the inverse scattering algorithm, fractures outside boreholes with metallic casing are modeled by a novel hybrid approximation method. This method combines the distorted Born approximation and the mixed order stabilized bi-conjugate gradient fast Fourier transform method to solve the forward scattering problem. The variational Born iterative method is applied to solve the nonlinear inverse problem iteratively. Experimental results show that the inverse scattering algorithm is effective for EM contrast enhanced through-casing hydraulic fracture evaluation.

Through-Casing Hydraulic Fracture Evaluation by Induction Logging I: An Efficient EM Solver for Fracture Detection

Hydraulic fracturing is an essential way to improve the production of unconventional shale oil and gas. It is important to characterize the produced fractures using either acoustic or electromagnetic (EM) methods. Conventional EM solvers in the low-frequency range face significant challenges by such multiscale problems where the fracture width is orders of magnitude smaller than its diameters. Furthermore, the cased borehole environment is extremely difficult to simulate with conventional EM solvers due to meshing difficulties and the multiscale nature of the problem. In this paper, we develop a hybrid distorted Born approximation and 3-D mixed ordered stabilized biconjugate gradient fast Fourier transform (DBA-BCGS-FFT) method to simulate the very challenging 2-D, 2-D-axisymmetric, and 3-D hydraulic fracture models under both open and cased borehole environments. Numerical examples show that this method has orders of magnitude higher efficiency than the finite element method. The capabilities of the DBA-BCGS-FFT method for the induction tool fracture mapping are demonstrated by comparing with laboratory experimental results and other reference results.

Preconditioned Multishift BiCG for $\mathcal{H}_2$-Optimal Model Reduction

We propose the use of a multishift biconjugate gradient method (BiCG) in combination with a suitable chosen polynomial preconditioning, to efficiently solve the two sets of multiple shifted linear ...

Fast Simulation of Scattering Problem for Magnetodielectric Materials With General Anisotropy in Layered Media

In this paper, the mixed-order stabilized biconjugate gradient fast Fourier transform (mixed-order BCGS-FFT) method is presented to solve the scattering problem of magnetodielectric materials with general anisotropy in layered media. While the volumetric roof-top functions are used as the testing functions for the coupled field volume integral equation and the basis functions for flux densities, the second-order curl conforming basis functions are applied to expand the vector potentials with the aim of both preserving the continuity of their tangential components and avoiding the zero terms that might otherwise be caused by the divergence operator. The layered medium Green’s function (LMGF) is efficiently evaluated through the recursive matrix method along with an interpolation technique. Several numerical experiments are presented to demonstrate the high accuracy and efficiency of the method. Different from the previously published works that aim to solve the similar problem, the new contribution of this work is to extend the mixed-order BCGS-FFT method to accommodate the layered background medium. Therefore, the 3-D FFT acceleration for integral kernels associated with the LMGF as well as the interpolation technique has been implemented and combined with the mixed-order BCGS-FFT method.

Numerical simulation of a Finite Moment Log Stable model for a European call option

Compared to the classical Black-Scholes model for pricing options, the Finite Moment Log Stable (FMLS) model can more accurately capture the dynamics of the stock prices including large movements or jumps over small time steps. In this paper, the FMLS model is written as a fractional partial differential equation and we will present a new numerical scheme for solving this model. We construct an implicit numerical scheme with second order accuracy for the FMLS and consider the stability and convergence of the scheme. In order to reduce the storage space and computational cost, we use a fast bi-conjugate gradient stabilized method (FBi-CGSTAB) to solve the discrete scheme. A numerical example is presented to show the efficiency of the numerical method and to demonstrate the order of convergence of the implicit numerical scheme. Finally, as an application, we use the above numerical technique to price a European call option. Furthermore, by comparing the FMLS model with the classical B-S model, the characteristics of the FMLS model are also analyzed.

On the solution of the Neumann Poisson problem arising from a compact differencing scheme using the full multi-grid method

A methodology for the numerical solution of the Neumann–Poisson problem for pressure that arises during the simulation of the incompressible Navier–Stokes equations on non-staggered grids is presented in this study. A sixth order compact differencing scheme has been used for discretizing the governing equation. A general procedure for implementing the discretized form of the integral constraint is proposed. Furthermore, different strategies for handling the corners in a rectangular domain are proposed and evaluated. Solutions for a model problem with a known analytical solution have been obtained on different grids using the full multi-grid technique with Bi-Conjugate Gradient Stabilized method as the smoother. Systematic order studies have been carried out. These bring out the fact that the overall order of the numerical solution is determined by the order of the discretization used for the boundary condition in the case of the Neumann–Poisson problem.

Fast Frequency-Domain Forward and Inverse Methods for Acoustic Scattering from Inhomogeneous Objects in Layered Media

The fast scattering and inverse scattering algorithms for acoustic wave propagation and scattering in a layered medium with buried objects are an important research topic, especially for large-scale geophysical applications and for target detection. There have been increasing efforts in the development of practical, accurate, and efficient means of imaging subsurface target anomalies. In this work, the acoustic scattering problem in layered media is formulated as a volume integral equation and is solved by the stabilized bi-conjugate gradient fast Fourier transform (BCGS-FFT) method. By splitting the layered medium Green’s function interacting with the induced source into a convolution and a correlation, the acoustic fields can be calculated efficiently by the FFT algorithm. This allows both the forward solution and inverse solution to be computed with only [Formula: see text] computation time per iteration, where [Formula: see text] is the number of degrees of freedom. The inverse scattering is solved using a simultaneous multiple frequency contrast source inversion (CSI). The stable convergence of this inversion process makes the multiple frequency simultaneous CSI reconstruction practical for large acoustic problems. Some representative examples are shown to demonstrate the effectiveness of the forward and inverse solvers for acoustic applications.

Variants of the groupwise update strategy for short-recurrence Krylov subspace methods

Krylov subspace methods often use short-recurrences for updating approximations and the corresponding residuals. In the bi-conjugate gradient (Bi-CG) type methods, rounding errors arising from the matrix---vector multiplications used in the recursion formulas influence the convergence speed and the maximum attainable accuracy of the approximate solutions. The strategy of a groupwise update has been proposed for improving the convergence of the Bi-CG type methods in finite-precision arithmetic. In the present paper, we analyze the influence of rounding errors on the convergence properties when using alternative recursion formulas, such as those used in the bi-conjugate residual (Bi-CR) method, which are different from those used in the Bi-CG type methods. We also propose variants of a groupwise update strategy for improving the convergence speed and the accuracy of the approximate solutions. Numerical experiments demonstrate the effectiveness of the proposed method.

Two dimension magnetotelluric modeling using finite element method, incomplete lu preconditioner and biconjugate gradient stabilized technique

Magnetotelluric (MT) method is a passive geophysical exploration technique utilizing natural electromagnetic source to obtain variation of the electric field and magnetic field on the surface of the earth. The frequency range used in this modeling is 10-4 Hz to 102 Hz. The two-dimensional (2D) magnetotelluric modeling is aimed to determine the value of electromagnetic field in the earth, the apparent resistivity, and the impedance phase. The relation between the geometrical and physical parameters used are governed by the Maxwell's equations. These equations are used in the case of Transverse Electric polarization (TE) and Transverse Magnetic polarization (TM). To calculate the solutions of electric and magnetic fields in the entire domain, the modeling domain is discretized into smaller elements using the finite element method, whereas the assembled matrix of equation system is solved using the Biconjugate Gradient Stabilized (BiCGStab) technique combined with the Incomplete Lower - Upper (ILU) preconditioner. This scheme can minimize the iteration process (computational cost) and is more effective than the Biconjugate Gradient (BiCG) technique with LU preconditions and Conjugate Gradient Square (CGS).

GPU-accelerated iterative solution of complex-entry systems issued from 3D edge-FEA of electromagnetics in the frequency domain

We present a performance analysis of a parallel implementation for both preconditioned conjugate gradient and preconditioned bi-conjugate gradient solvers running on graphic processing units (GPUs) with CUDA programming model. The solvers were mainly optimized for the solution of sparse systems of algebraic equations at complex entries, arising from the three-dimensional edge-finite element analysis of the electromagnetic phenomena involved in the open-bound earth diffusion of currents under time-harmonic excitation. We used a shifted incomplete Cholesky (IC) factorization as preconditioner. Results show a significant speedup by using either a single-GPU or a multi-GPU device, compared to a serial central processing unit (CPU) implementation, thereby allowing the simulations of large-scale problems in low-cost personal computers. Additional experiments of the optimized solvers show that its use can be extended successfully to other complex systems of equations arising in electrical engineering, such as those obtained in power–system analysis.

High performance computing in modelling and simulation

With the arrival of high performance computing (HPC) technologies, one of the most well-known scientific fields that succeeded in taking advantage of the new skyrocketing HPC opportunities was the modelling and simulation of computationally intensive scientific applications. It is now clear that the development of such models through which computers can extensively simulate the evolution of artificial and natural systems is fundamental for the advancement of science. As a result, different HPC approaches have allowed researchers nowadays to considerably extend the application of computing methodologies in research and industry, but also extend it to the quantitative study of complex phenomena. Moreover, this has permitted a broader and more effective application of numerical methods for differential equation systems (e.g. finite element method (FEM), finite difference method (FDM), etc.) on the one hand, and, most importantly, the application of alternative computational paradigms, such as cellular automata (CA), genetic algorithms, neural networks, swarm intelligence, and so on, on the other. The latter have demonstrated their effectiveness for modelling purposes when traditional simulation methodologies have proven to be impracticable. In conjunction with high parallel computer architectures, like the ones found in graphics processing units (GPUs), with thousands of computational threads on each of them, they are applied as fine candidates to overcome the computational burden imposed by vastly different scientific applications in diverse fields, such as engineering, physics, chemistry, biology, geology, medicine, ecology, sociology, traffic control, economy, and so on. Nevertheless, the resulting new parallelization paradigms and novel parallel computational approaches for all the aforementioned applications are considered one of the hot topics of today’s research. Having all this in mind, this special issue aims to provide a platform for a multidisciplinary community composed of scholars, researchers, developers, educators, practitioners and experts from world-leading universities, institutions, agencies and companies in computational science, and thus in the high performance computing for modelling and simulation field. As a sequel to the special session on high performance computing in modelling and simulation (HPCMS) within the 22nd Euromicro international conference on parallel, distributed and network-based processing (PDP), held in Turin, Italy on 12–14 February 2014, our goal in organizing this special issue of the ‘International Journal of High Performance Computing Applications’ was to gather the most suitable publications advancing HPCMS from both a theoretical and an application point of view and which were presented in an earlier reduced version in PDP 2014. We are pleased to say that this goal was achieved by the four excellent contributions presented in this issue, briefly summarized in the following. More specifically, this special issue addresses inherent parallel models and complex numerical models and their implementations by using GPUs as a fundamental computational primitive in various application fields like electromagnetics, environmental modelling and aircraft evacuation. Moreover, the evaluation of a one-instruction-per-cycle (one-IPC) high-level simulator of the microthreaded many-core systems is provided in detail while integrating the architecture model with the application model by abstracting the details of instruction execution. The first paper, by Camargos et al., provides a performance analysis of a parallel implementation for both preconditioned conjugate gradient and preconditioned bi-conjugate gradient solvers running on GPUs with a compute unified device architecture (CUDA) programming model. As a result, they focused on a GPU-accelerated iterative solution of complex-entry systems issued from 3D edge-Finite Element Analysis (FEA) of electromagnetics in the frequency domain. Their results

Fast permutation preconditioning for fractional diffusion equations.

In this paper, an implicit finite difference scheme with the shifted Grünwald formula, which is unconditionally stable, is used to discretize the fractional diffusion equations with constant diffusion coefficients. The coefficient matrix possesses the Toeplitz structure and the fast Toeplitz matrix-vector product can be utilized to reduce the computational complexity from {mathcal {O}}{(N^{3})} to {mathcal {O}}{(N log N)}, where N is the number of grid points. Two preconditioned iterative methods, named bi-conjugate gradient method for Toeplitz matrix and bi-conjugate residual method for Toeplitz matrix, are proposed to solve the relevant discretized systems. Finally, numerical experiments are reported to show the effectiveness of our preconditioners.

2D Zero-Inertia Model for Solution of Overland Flow Problems in Flexible Meshes

AbstractA study of the efficiency of a zero-inertia model (ZI) for two-dimensional (2D) overland flow simulation is presented in this work. An upwind numerical scheme is used for the spatial discretization in the frame of finite-volume methods and an implicit formulation is chosen to avoid numerical instability. The scheme is applied in both structured and unstructured meshes, focusing in the latter ones due to their good adaptability. The ZI equation has a nonlinear character; hence, a linearization is required in the implicit procedure. This is carried out by means of Picard iterations method as a previous step to the system matrix resolution, characteristic of implicit techniques. The BiConjugate Gradient Stabilized (BiCGStab) method combined with sparse storage strategies is selected for the system resolution. A dual-threshold incomplete lower upper factorization (ILUT) is chosen as matrix preconditioner. Computational efficiency of the implicit temporal discretization for ZI model is explored under b...

Bi-conjugate gradient based computation of weight vector in space-time adaptive processing

在空时自适应信号处理中, 权重向量是通过杂波协方差矩阵的逆和一个导向向量的乘积得到的。为了减少计算权重向量的复杂度, 很多研究都是为了降低计算杂波协方差矩阵逆的计算量。本文将计算该权重向量的计算归结为求解一个线性方程组解, 然后用双共轭梯度法求解该权重向量。分析得到该方法的计算复杂度比直接求逆的方法降低了一个级别。

A time-integration method for stable simulation of extremely deformable hyperelastic objects

This paper presents a time integration method for realtime simulation of extremely deformable objects subject to geometrically nonlinear hyperelasticity. In the presented method, the equation of motion of the system is discretized by the backward Euler method, and linearly approximated through the first-order Taylor expansion. The approximate linear equation is solved with the quasi-minimal residual method (QMR), which is an iterative linear equation solver for non-symmetric or indefinite matrices. The solution is then corrected considering the nonlinear term that is omitted at the Taylor expansion. The method does not demand the constitutive law to guarantee the positive definiteness of the stiffness matrix. Experimental results show that the presented method realizes stable behavior of the simulated model under such deformation that the tetrahedral elements are almost flattened. It is also shown that QMR outperforms the biconjugate gradient stabilized method (BiCGStab) in this application.

BiCGCR2: A new extension of conjugate residual method for solving non-Hermitian linear systems

In the present paper, we introduce a new extension of the conjugate residual (CR) method for solving non-Hermitian linear systems with the aim of developing an alternative basic solver to the established biconjugate gradient (BiCG), biconjugate residual (BiCR) and biconjugate A-orthogonal residual (BiCOR) methods. The proposed Krylov subspace method, referred to as the BiCGCR2 method, is based on short-term vector recurrences and is mathematically equivalent to both BiCR and BiCOR. We demonstrate by extensive numerical experiments that the proposed iterative solver has often better convergence performance than BiCG, BiCR and BiCOR. Hence, it may be exploited for the development of new variants of non-optimal Krylov subspace methods.

Parallel FETI‐DP algorithm for efficient simulation of large‐scale EM problems

SummaryAn efficient parallelization of the dual‐primal finite‐element tearing and interconnecting (FETI‐DP) algorithm is presented for large‐scale electromagnetic simulations. As a nonoverlapping domain decomposition method, the FETI‐DP algorithm formulates a global interface problem, whose iterative solution is accelerated with a solution of a global corner problem. To achieve a good load balance for parallel computation, the original computational domain is decomposed into subdomains with similar sizes and shapes. The subdomains are then distributed to processors based on their close proximity to minimize inter‐processor communication. The parallel generalized minimal residual method, enhanced with the iterative classical Gram‐Schmidt orthogonalization scheme to reduce global communication, is adopted to solve the global interface problem with a fast convergence rate. The global corner‐related coarse problem is solved iteratively with a parallel communication‐avoiding biconjugate gradient stabilized method to minimize global communication, and its convergence is accelerated by a diagonal preconditioner constructed from the coarse system matrix. To alleviate neighboring communication overhead, the non‐blocking communication approach is employed in both generalized minimal residual and communication‐avoiding biconjugate gradient stabilized iterative solutions. Three numerical examples are presented to demonstrate the accuracy, scalability, and capability of the proposed parallel FETI‐DP algorithm for electromagnetic modeling of general objects and antenna arrays. Copyright © 2016 John Wiley & Sons, Ltd.