Research on the conjugate gradient algorithm with a modified incomplete Cholesky preconditioner on GPU

Highlights

• We present a parallel method of the forward/backward substitutions on the GPU.

• A new kernel of SpMV with high optimization on the GPU is proposed.

• We suggest a new kernel of inner product on the GPU.

• An efficient MIC preconditioned conjugate gradient algorithm on the GPU is presented.

Abstract

In this study, we discover the parallelism of the forward/backward substitutions (FBS) for two cases and thus propose an efficient preconditioned conjugate gradient algorithm with the modified incomplete Cholesky preconditioner on the GPU (GPUMICPCGA). For our proposed GPUMICPCGA, the following are distinct characteristics: (1) the vector operations are optimized by grouping several vector operations into single kernels, (2) a new kernel of inner product and a new kernel of the sparse matrix–vector multiplication with high optimization are presented, and (3) an efficient parallel implementation of FBS on the GPU (GPUFBS) for two cases are suggested. Numerical results show that our proposed kernels outperform the corresponding ones presented in CUBLAS or CUSPARSE, and GPUFBS is almost 3 times faster than the implementation of FBS using the CUSPARSE library. Furthermore, GPUMICPCGA has better behavior than its counterpart implemented by the CUBLAS and CUSPARSE libraries.

Introduction

Graphic Processing Units (GPUs), originally developed for computer games, now provide computational power for scientific applications and, with their many core structures, provide a serious opportunity for scientific computations. Processing large data sets in parallel on GPUs has drawn much attentions in recent years. Following the introduction of CUDA (Compute Unified Device Architecture) by NVIDIA in 2007  [16], a programming model designed to support joint CPU/GPU execution of applications, GPUs have become increasingly serious competitors among the general-purpose parallel programming systems.

The conjugate gradient (CG) algorithm is one of the best known iterative methods. With a suitable preconditioner, the performance of the CG algorithm can be dramatically improved. The preconditioned conjugate gradient (PCG) algorithm has proven its efficiency and robustness in a wide range of applications. Recently, many researchers have attempted to develop the suitable and flexible PCG algorithm for the GPU architecture. In  [3], an efficient implementation of a Jacobi preconditioned conjugate gradient algorithm on GPUs is developed. However, the Jacobi preconditioner has a limited impact on the efficiency and robustness of the PCG algorithm on the whole. In  [10], the authors give a parallel implementation of the SSOR (Symmetric Successive Over-Relaxation) PCG algorithm by utilizing a first order approximation. Numerical results show as compared to the CPU implementation of the CG algorithm, their GPU PCG implementation is up to 10 times faster. A new parallel computational technique is proposed for the parallelization of the explicit inverse preconditioned biconjugate conjugate gradient type method on GPUs in  [9]. Other related work can be found in  [7], [13], [14].

From the studies mentioned above, it can be found that for each PCG algorithm, one of the key factors is how to parallel solve the equation $Mz=r$ ($M$ is the preconditioned matrix) on GPUs when shifting it to the GPU platform. For the modified incomplete Cholesky factorization (MIC) preconditioning, it has proven its effectiveness in improving the performance of the CG algorithm, and $M=LD{L}^T$, where $L$ is a lower triangular matrix. However, for solving the equation $LD{L}^{T}z=r$, the following two steps are required. First, $LDv=r$ is solved for $v$ by the forward substitution, and then $D{L}^{T}z=Dv$ is solved for $z$ by the backward substitution. It is obviously observed that the forward/backward substitutions are not easy to implement on GPUs, which to a great extent prevents us from shifting the PCG algorithm with the MIC preconditioner to the GPU platform.

Therefore, in this study, we will discuss how to implement the MIC PCG algorithm on the GPU and emphasize whether the forward/backward substitutions can be parallelized on the GPU. Fortunately, in our recent studies, we observe that because $L$ may have the same non-zero entries as the coefficient matrix of the linear system, the forward/backward substitutions can be better parallelized on the GPU if the coefficient matrix of the linear system satisfies some conditions.

Next, we take the following partial differential equation (PDE) for example to represent our recent findings and thus propose an efficient parallel implementation of the MIC PCG algorithm on the GPU:

$$u_t=\Delta u\text{,}\text{in}\Omega \times [0,T]\text{,}$$

where $\Omega$ is a regular domain in ${\mathbb{R}}^d$, $d=2,3$, and $\Delta$ is the Laplace operator, and $u$ is the unknown solution. A solution $u$ of the PDE (1) must satisfy the boundary conditions (Dirichlet, Neumann or Cauchy) and the initial conditions.

For $d=2$, we utilize the discrete variational derivative method (DVDM) to discretize the PDE (1). For readers who are interested in the DVDM, please refer to the literature  [6]. Assuming that a 3×3 two-dimensional grid, shown in Fig. 1, is defined, the scheme obtained above can be written as the following linear system for a certain time $t\in [0,T]$:

$$A_{\left(2\right)}{u}_{\left(2\right)}=b_{\left(2\right)}\text{,}$$

where

$$A_{\left(2\right)}=\left[\begin{array}{ccccccccc}\hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \ast \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \ast \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill \end{array}\right]\text{,}$$

the coefficient matrix $A_{\left(2\right)}$ is a sparse, symmetric, positive definite, five-diagonal, diagonally dominant 9×9 matrix of regular structure, $b_{\left(2\right)}$ is a vector of size 9, and $u_{\left(2\right)}$ is an unknown solution vector of size 9.

According to the definition of the MIC preconditioning, the preconditioned matrix $M_{\left(2\right)}$ of $A_{\left(2\right)}=L_{\left(2\right)}{D}_{\left(2\right)}{L}_{\left(2\right)}^T$, where

$$L_{\left(2\right)}=\left[\begin{array}{ccccccccc}\hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill \end{array}\right]\text{,}$$

a lower triangular matrix with non-zero values along the main diagonal and at off-diagonal locations only where $A_{\left(2\right)}$ has non-zero entries.

Generally, the forward/backward substitutions both need $n$ ($n$ is the element number of the solution vector) steps to finish calculations of elements. However, it can be observed that the forward substitution for this case both can be finished in 5 steps instead of 9 steps according to the progression shown in Fig. 2.

The computational grid in Fig. 2 is the same as in Fig. 1, and arrows represent dependencies between elements of the grid. From Fig. 2, we can see that the calculation of an element in each step is dependent to its two neighbors (adjacent row and column) and these neighbors must be calculated earlier than it. Also note that not all elements in Fig. 2 have these neighbor elements. For example, the element 7 of the forward substitution is only dependent to its row neighbor element 4. Moreover, the elements of each step are independent of each other and can be concurrently calculated. These observations are also suitable for the backward substitution. All these things are helpful for shifting the forward/backward substitutions to the GPU platform.

In the above case, we investigate a class of five-diagonal matrices like $A_{\left(2\right)}$ and introduce how to decrease the number of computational steps and the parallelism embraced in them for the forward/backward substitutions. On the other hand, let us consider the three-dimensional PDE ($d=3$) in order to investigate a class of seven-diagonal matrices.

Assuming that a 4×2×2 three-dimensional grid, shown in Fig. 3, is defined, the PDE (1) ($d=3$) is discretized by the DVDM and then is written as the following linear system for a certain time $t\in [0,T]$:

$$A_{\left(3\right)}{u}_{\left(3\right)}=b_{\left(3\right)}\text{,}$$

where

$$A_{\left(3\right)}=\left[\begin{array}{cccccccccccccccc}\hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \ast \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill \\ \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill & \hfill \ast \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \hfill & \hfill \hfill & \hfill \ast \hfill & \hfill \ast \hfill \end{array}\right]\text{,}$$

the coefficient matrix $A_{\left(3\right)}$ is a sparse, symmetric, positive definite, seven-diagonal, diagonally dominant 16×16 matrix of regular structure, $b_{\left(3\right)}$ is a vector of size 16, and $u_{\left(3\right)}$ is an unknown solution vector of size 16. Similarly, the MIC preconditioned matrix $M_{\left(3\right)}=L_{\left(3\right)}{D}_{\left(3\right)}{L}_{\left(3\right)}^T$ can be obtained, and

$$L_{\left(3\right)}=\left[\begin{array}{cccccccccccccccc}\hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill \hfill & \hfill \hfill & \hfill ♢\hfill & \hfill 1\hfill \end{array}\right]\text{.}$$

By analyzing the lower triangular matrix $L_{\left(3\right)}$, we conclude that the computational steps of the forward/backward substitutions both can decrease from 16 to 6 for this case. In Fig. 4, we take the forward substitution for example to show its computational progression.

As compared to the case in Fig. 2, here the calculation of an element in each step does not only depend on its adjacent row and column neighbors, but also depends on its adjacent layer neighbors. For example, the element 14 is dependent on its adjacent row and column elements 10 and 13 and its adjacent layer element 6. However, not all elements have these neighbor elements and the number of neighbors varies with the location of the element. Moreover, for any element, its dependent neighbors must be calculated earlier than it in order to start its calculation. Like the case in Fig. 2, the elements of each step are independent of each other and can be concurrently calculated.

In the following sections, we will furthermore discuss the implementation of the MIC PCG algorithm on the GPU. Compared to other related work, our researches have the following main contributions.

$•$  For the MIC preconditioning, the parallelism of the forward/backward substitutions for two cases is found, and thus an efficient parallel method on the GPU is presented.

$•$  A new method of inner product and a new method of the sparse matrix–vector multiplication (SpMV) on the GPU are proposed. For our proposed SpMV on the GPU, it is optimized by shared memory and fully coalesced global memory.

$•$  Based on the above methods, an optimized MIC PCG algorithm with double precision on the GPU is suggested.

The reminder of this paper is organized as follows. In the second section, the MIC PCG algorithm and the CUDA architecture are described. In the third section, some GPU kernels for the MIC PCG algorithm are presented. Numerical results are presented in the fourth section. The fifth section contains our conclusions and points to our future research direction.