# Developing a Multi-GPU-Enabled Preconditioned GMRES with Inexact Triangular Solves for Block Sparse Matrices

Authors: Wenpeng Ma; Yiwen Hu; Wu Yuan; Xiazhen Liu
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Solving triangular systems is the building block for preconditioned GMRES algorithm. Inexact preconditioning becomes attractive because of the feature of high parallelism on accelerators. In this paper, we propose and implement an iterative, inexact block triangular solve on multi-GPUs based on PETSc’s framework. In addition, by developing a distributed block sparse matrix-vector multiplication procedure and investigating the optimized vector operations, we form the multi-GPU-enabled preconditioned GMRES with the block Jacobi preconditioner. In the implementation, the GPU-Direct technique is employed to avoid host-device memory copies. The preconditioning step used by PETSc’s structure and the cuSPARSE library are also investigated for performance comparisons. The experiments show that the developed GMRES with inexact preconditioning on 8 GPUs can achieve up to 4.4x speedup over the CPU-only implementation with exact preconditioning using 8 MPI processes.


---

## Body

## 1. Introduction

Solving a large sparse linear system of equations is always necessary in scientific applications. For unsymmetric matrices, the Krylov subspace-based Generalized Minimal Residual (GMRES) algorithm [1] is widely used as an effective linear solver. Since the general purpose computing on GPUs was introduced in 2007 by NVIDIA, various GPU-enabled GMRES algorithms were designed using NVIDIA’s CUDA (Compute Unified Device Architecture). Li and Saad [2] investigated preconditioned Conjugate Gradient (CG) and GMRES methods on TESLA C1060, and 4.0x speedup can be obtained for GMRES algorithm. Khodja et al. [3] implemented the GMRES algorithm on a GPU cluster by exploiting Message Passing Interface (MPI) and CUDA, and they also focused on minimizing the communication between processes using the compressed storage and hypergraph partitioning techniques. Yamazaki et al. [4, 5] developed various procedures to perform the orthogonalization steps of GMRES efficiently and introduced a serial of preconditioners within domain decomposition methods for communication avoiding GMRES (CA-GMRES) on a hybrid CPU + GPU cluster. Their CPU + GPU implementation could obtain a speedup of 7.4x and 1.7x over CA-GMRES without preconditioning and with preconditioning, respectively. Yang [6, 7] developed the preconditioned GMRES algorithm by parallelizing the ILU(0), ILUT, block ILU(k), and triangular solves on GPUs. Gao et al. [8] proposed an efficient GPU kernel on the sparse matrix-vector multiplication (SpMV) in GMRES and applied the optimized GMRES to solving the two-dimensional Maxwell's equations. He et al. [9] presented an efficient GPU implementation of the GMRES with ILU preconditioners for solving large linear dynamic systems. They obtained 3.0-12.0x speedup over the CPU implementation. Other studies on the Krylov subspace methods on GPUs can be found in [10–19].Preconditioning techniques are effective for accelerating the convergence of GMRES algorithm. ILU factorization-based preconditioning scheme [20] is one of the most widely used methods. In this method, the ILU factorizations and triangular solves need to be performed. However, they are inherently sequential procedures and have strong data dependency. Many parallel approaches have been studied, for example, the level-scheduling algorithms [2, 6, 21] and inexact methods [22–27]. The level-scheduling schemes offer limited parallelism on GPU, and the performance highly depends on the sparsity of the matrix. In contrast, inexact methods can oﬀer sufficient parallelism and become attractive on accelerators. Chow et al. [23, 24] presented an iterative, ﬁne-grained strategy for computing ILU factorizations, and their GPU implementation achieved up to 27x over the level-scheduling approach. Anzt et al. [26] used the iterative and inexact idea to solve triangular systems on GPUs, and they showed an advantage over exact triangular solves in terms of total computing time of GMRES although inexact implementation resulted in more GMRES iterations.Block sparse matrices are very common in scientiﬁc computing areas, especially in multiphysics simulations. Most of the existing parallel GMRES and preconditioning algorithms mentioned above were designed for pointwise matrices. Motivated by the block matrix applications and the highly parallelized inexact preconditioning technique, in this paper, we focus on developing a multi-GPU-enabled preconditioned GMRES with the inexact block preconditioning algorithm for block sparse matrices. Our contributions are as follows:(1)
Based on PETSc’s framework, we develop the distributed block sparse matrix-vector multiplication for GMRES on multi-GPUs using the cuSPARSE library [28] and the GPU-Direct technique. And two optimization strategies are investigated for vector operations in the GMRES method.(2)
We proposed a strategy that implements the iterative, inexact sparse triangular solves for block matrices on GPUs and integrated it into PETSc. In order to make performance comparisons, we also develop other two versions of the preconditioning technique for block sparse matrices.(3)
By employing the block Jacobi method, we investigate the preconditioned GMRES with various block preconditioning methods on a multi-GPU system. Our experiments show that the developed multi-GPU-based GMRES can achieve a speedup of up to 4.4x using 8 GPUs over the CPU-only implementations using 8 MPI processes.The remainder of this paper is structured as follows. Section2 introduces the definition of the block sparse matrix with the storage format and describes a block-based right-preconditioned GMRES algorithm. The optimization strategies, iterative inexact block sparse triangular solves, and three implementations for the preconditioning step are presented and discussed in Section 3. In Section 4, we conduct many experiments to make detailed performance comparisons between different strategies and implementations. Section 5 gives a summary of this work.


## 2. Background

### 2.1. Block Sparse Matrix

Given a sparsity patternSp, we consider a general n×n block sparse matrix with block size s as(1)An×nr,s=Ai,j,i,j=1,2,3,…,r,where r is the number of block rows and columns, Ai,j is a s×s dense block for i,j∈Sp, and Ai,j=0 for i,j∉Sp.To conduct numerical operations on block sparse matrices, an eﬀective and eﬃcient storage algorithm is required. The well-known method for storing a block sparse matrix is the Block Compressed Sparse Row (BCSR) format [28, 29]. Currently, the popular numerical libraries such as Intel MKL [30], NVIDIA’s cuSPARSE [28], and PETSc [31] support BCSR format. In this format, the block sparse matrix An×nr,s with nnzb nonzero blocks is represented by block rows using three arrays rowptr, colval, and blkval. We assume the indexing starts from 0 in C programming language:(1)
rowptr is of size r+1, and all the column indices for the ithi<r block row is recorded from rowptri to rowptri+1 (not including) in colval array.(2)
colval is of size nnzb, and stores all column indices by block rows.(3)
blkval is of size nnzb×s2, and stores all the values of nonzero blocks by block rows. The values within a s×s block are stored consecutively.However, the BCSR format cloud be speciﬁed in row-major or column-major order depending on which index is firstly ordered within a block. The following matrix:(2)14151300171157000052027200033151015360327121643161319116×63,2,rowptr: 0 2 4 7, colval: 0 1 1 2 0 1 2, blkval: {14, 17, 15, 11}, {1, 5, 3, 7}, {5, 33, 20, 15}, {27, 10, 20, 15}, {36, 4, 0, 31}, {32, 6, 7, 13}, {12, 19, 16, 11} illustrates an example of the column-major BCSR format where four values within blocks are stored by columns. In this paper, we employ column-major BCSR format in accordance with PETSc.


### 2.2. Preconditioned GMRES Algorithm

The GMRES method, proposed by Saad and Schultz [1] for solving linear systems, is an iterative process where the residual vector of a linear system is minimized over a Krylov subspace at every iteration. In this paper, we consider the computation of a block sparse linear system as(3)ABx⟶=b⟶,where B=n,r,s represents the number of rows and columns, the number of block rows and columns, and block size of the block sparse matrix, respectively.To accelerate the convergence of GMRES, we apply right-preconditioning technique to the linear system as(4)ABMB−1x⟶d=b⟶,x⟶=MB−1x⟶d,where MB is the block preconditioning matrix.Similar to the GMRES for pointwise matrices, it is straightforward to write the preconditioned GMRES algorithm in the block format. The process is shown in Algorithm1. Compared to the pointwise case, the algorithm performs block sparse vector multiplication (line 7 in Algorithm 1) and block preconditioning steps (line 6 and line 21 in Algorithm 1) instead.Algorithm 1:GMRES algorithm with right preconditioning for block sparse matrices.
Input: (1) a block sparse matrix, AB; (2) the preconditioning matrix, MB; (3) the right-hand side, b⟶; (4) the initial guess, x⟶0; (5) the relative tolerance, ϵ; (6) the maximum number of iterations, maxit; (7) the restart number, m.Output: the solution vector x⟶.(1)
continue = 1; it = 0; allocate Hm=Hp,q,0≤p≤m,0≤q<m;(2)
whilecontinuedo(3)
r⟶0=b⟶−ABx⟶0; β=r⟶02; v⟶0=r⟶0/β;(4)
if (it == 0) thennorm0=β;(5)
fori=0tom−1do(6)
z⟶=MB−1v⟶i;(7)
w⟶=ABz⟶;(8)
fork=0toido(9)
Hk,i=w⟶,v⟶k;(10)
w⟶=w⟶−Hk,iv⟶k;(11)
end(12)
normw=w⟶2; Hi+1,i=normw; t=i;(13)
updateHessenbergH,norm;(14)
ifnormw==0orlog10norm/norm0<log10ε or ++it==maxitthen(15)
continue=0; break;(16)
end(17)
v⟶i+1=w⟶/normw;(18)
end(19)
t=t+1;(20)
solvey⟶ from argminy⟶βe⟶1−Hty⟶2, where Ht=Hp,q,0≤p≤t,0≤q<t(21)
x⟶0=x⟶0+MB−1Vty⟶, where Vt=v⟶0,v⟶1,…,v⟶t−1;(22)
end(23)
x=x0When equation (4) is solved in a whole domain, the preconditioning step can be realized by conducting the incomplete LU (ILU) factorization [6, 20, 23] of AB and performing two triangular solves [20, 21, 26] as(5)AB≈LBUB=MB,(6)MB−1v⟶=z⟶⟶LBf⟶=v⟶,UBz⟶=f⟶,where LB is the lower triangular factor and UB is the upper triangular factor. However, the ILU factorizations are used as subdomain preconditioners within the domain decomposition methods such as the Additive Schwarz method (ASM) [32] which becomes the block Jacobi method without overlaps. And, the block Jacobi preconditioner is used for the GMRES solver on a distributed system in this work.


## 2.1. Block Sparse Matrix

Given a sparsity patternSp, we consider a general n×n block sparse matrix with block size s as(1)An×nr,s=Ai,j,i,j=1,2,3,…,r,where r is the number of block rows and columns, Ai,j is a s×s dense block for i,j∈Sp, and Ai,j=0 for i,j∉Sp.To conduct numerical operations on block sparse matrices, an eﬀective and eﬃcient storage algorithm is required. The well-known method for storing a block sparse matrix is the Block Compressed Sparse Row (BCSR) format [28, 29]. Currently, the popular numerical libraries such as Intel MKL [30], NVIDIA’s cuSPARSE [28], and PETSc [31] support BCSR format. In this format, the block sparse matrix An×nr,s with nnzb nonzero blocks is represented by block rows using three arrays rowptr, colval, and blkval. We assume the indexing starts from 0 in C programming language:(1)
rowptr is of size r+1, and all the column indices for the ithi<r block row is recorded from rowptri to rowptri+1 (not including) in colval array.(2)
colval is of size nnzb, and stores all column indices by block rows.(3)
blkval is of size nnzb×s2, and stores all the values of nonzero blocks by block rows. The values within a s×s block are stored consecutively.However, the BCSR format cloud be speciﬁed in row-major or column-major order depending on which index is firstly ordered within a block. The following matrix:(2)14151300171157000052027200033151015360327121643161319116×63,2,rowptr: 0 2 4 7, colval: 0 1 1 2 0 1 2, blkval: {14, 17, 15, 11}, {1, 5, 3, 7}, {5, 33, 20, 15}, {27, 10, 20, 15}, {36, 4, 0, 31}, {32, 6, 7, 13}, {12, 19, 16, 11} illustrates an example of the column-major BCSR format where four values within blocks are stored by columns. In this paper, we employ column-major BCSR format in accordance with PETSc.


## 2.2. Preconditioned GMRES Algorithm

The GMRES method, proposed by Saad and Schultz [1] for solving linear systems, is an iterative process where the residual vector of a linear system is minimized over a Krylov subspace at every iteration. In this paper, we consider the computation of a block sparse linear system as(3)ABx⟶=b⟶,where B=n,r,s represents the number of rows and columns, the number of block rows and columns, and block size of the block sparse matrix, respectively.To accelerate the convergence of GMRES, we apply right-preconditioning technique to the linear system as(4)ABMB−1x⟶d=b⟶,x⟶=MB−1x⟶d,where MB is the block preconditioning matrix.Similar to the GMRES for pointwise matrices, it is straightforward to write the preconditioned GMRES algorithm in the block format. The process is shown in Algorithm1. Compared to the pointwise case, the algorithm performs block sparse vector multiplication (line 7 in Algorithm 1) and block preconditioning steps (line 6 and line 21 in Algorithm 1) instead.Algorithm 1:GMRES algorithm with right preconditioning for block sparse matrices.
Input: (1) a block sparse matrix, AB; (2) the preconditioning matrix, MB; (3) the right-hand side, b⟶; (4) the initial guess, x⟶0; (5) the relative tolerance, ϵ; (6) the maximum number of iterations, maxit; (7) the restart number, m.Output: the solution vector x⟶.(1)
continue = 1; it = 0; allocate Hm=Hp,q,0≤p≤m,0≤q<m;(2)
whilecontinuedo(3)
r⟶0=b⟶−ABx⟶0; β=r⟶02; v⟶0=r⟶0/β;(4)
if (it == 0) thennorm0=β;(5)
fori=0tom−1do(6)
z⟶=MB−1v⟶i;(7)
w⟶=ABz⟶;(8)
fork=0toido(9)
Hk,i=w⟶,v⟶k;(10)
w⟶=w⟶−Hk,iv⟶k;(11)
end(12)
normw=w⟶2; Hi+1,i=normw; t=i;(13)
updateHessenbergH,norm;(14)
ifnormw==0orlog10norm/norm0<log10ε or ++it==maxitthen(15)
continue=0; break;(16)
end(17)
v⟶i+1=w⟶/normw;(18)
end(19)
t=t+1;(20)
solvey⟶ from argminy⟶βe⟶1−Hty⟶2, where Ht=Hp,q,0≤p≤t,0≤q<t(21)
x⟶0=x⟶0+MB−1Vty⟶, where Vt=v⟶0,v⟶1,…,v⟶t−1;(22)
end(23)
x=x0When equation (4) is solved in a whole domain, the preconditioning step can be realized by conducting the incomplete LU (ILU) factorization [6, 20, 23] of AB and performing two triangular solves [20, 21, 26] as(5)AB≈LBUB=MB,(6)MB−1v⟶=z⟶⟶LBf⟶=v⟶,UBz⟶=f⟶,where LB is the lower triangular factor and UB is the upper triangular factor. However, the ILU factorizations are used as subdomain preconditioners within the domain decomposition methods such as the Additive Schwarz method (ASM) [32] which becomes the block Jacobi method without overlaps. And, the block Jacobi preconditioner is used for the GMRES solver on a distributed system in this work.


## 3. Multi-GPU Implementations

### 3.1. Chart of Development

PETSc [31] offers a wide range of Application Programming Interfaces (APIs) for users to manipulate low-level data structures. Furthermore, it uses function pointers in standard C language instead of an Object-Oriented language such as C++ to achieve data encapsulation and polymorphism [33]. We take advantage of function pointers and structures in C language to develop a user friendly multi-GPU version of the preconditioned GMRES algorithm.The chart of our code development is shown in Figure1. Data structure GMRESInfo is designed to store all necessary data ﬁelds on the GPU. The memory spaces of the fields are allocated when an instance of GMRESInfo is created. Before the GMRES solver is executed, the fields in GMRESInfo are initialized on GPUs by launching the memory copy operations provided by the basic CUDA library from the host to device. The encapsulated data on the host that needs to be transferred to GPUs can be accessed using the low-level APIs in PETSc. To make full use of the existing numerical toolkits, we employ the efficient functions in cuBLAS and cuSPARSE libraries to perform the local sparse matrix and vector operations in Algorithm 1. On a multi-GPU platform, a portion of the computation inevitably requires data exchange among GPUs. This is implemented and optimized using the CUDA-aware MPI library that supports the GPU-Direct technology and transfers data directly among GPU memories through network adapters.Figure 1
The development chart of the multi-GPU version of GMRES algorithm based on PETSc.To adopt the multi-GPU version of GMRES solver in an application, a user needs to call three additional functions. Two of them creates and destroys aGMRESInfo instance, respectively, while the other (SetUserDeﬁnedGMRES) sets the function pointer of GMRES to our developed procedure.To implement Algorithm1 on a multi-GPU system, the workload for the host and device should be proper assigned. The matrix-vector operations (line 3 and line 7) and the basic vector operations such as inner product (line 9 and line 12) and the linear combinations (line 10 and line 17) are of high parallelism, so they are fully performed on the GPUs. To conduct convergence monitoring, the conditional statements (line 14–16) are inevitable. And, this part is executed on the host. Because cuSPARSE library allows users to pass a reference on the host when an API is executed on the GPU, the scalar parameter normw is declared on the host to receive the result of w⟶2 and can be directly used in line 14. The size of the Hessenberg matrix is (m + 1) × m, where m (m is set to 30 in PETSc by default) is usually much smaller than the matrix size, so the updating of Hessenberg matrix (line 13) and solving of the least square problem (line 20) are performed on the host.To distributedly store the data on a multi-GPU system, block matrices and vectors are partitioned by rows. On a GPU cluster withN GPUs, the ith GPU stores ri block rows of AB which is represented in BCSR format. In addition, there are several vectors that need to be allocated on the ith GPU. r⟶0i, z⟶i, x⟶i, and w⟶i are allocated with the size of ri×s, and Vt stores up to m vectors each of which is of size ri×s. In the preconditioning step (line 6), the preconditioning matrix on the ith GPU, denoted as MBi, is stored in separate factors, i.e., LBi and UBi, and the size of memory requirements for ILU(0) level equals to the size of ABi, where ABi is the main diagonal part of the partitioned matrix on the ith GPU.


### 3.2. Implementation and Optimization

The major computation in Algorithm1 consists of block sparse matrix-vector multiplication (BSpMV), the preconditioning step, and the vector operations. For a block sparse matrix on a distributed system, PETSc partitions the matrix by consecutive block rows and stores each partition in two separate structures. One structure is for the main diagonal of the partitioned matrix, the other is for the oﬀ diagonals. Figure 2 shows the strategy of partitioning and storing a 18s×18s (s is the block size) block sparse matrix obtained by the 5-point finite diﬀerencing scheme in PETSc. Each partition is responsible for storing 6 block rows and assigned to a MPI process. Take the second partition (dashed box in red) as an example, PETSc stores the 6s×6s square block sparse matrix in the main diagonal as the BCSR format by renumbering the block rows and columns starting from zero. The two oﬀ-diagonal parts are merged and renumbered as a local 6s×12s block sparse matrix in the BCSR format. The main advantage of this strategy is that it separates the main diagonal part (MDP, communication-free) from the oﬀ-diagonal part (ODP, communication-required), which can be used for overlapping the computation and communication. Therefore, each partition needs three steps to perform the multiplication between a block sparse matrix and a vector. The MDP firstly multiplies the local part of the vector and stores the result temporarily. All partitions then launch MPI sending and receiving calls to exchange data of the vector globally. Lastly, the ODP in each partition multiplies the received vector and accumulates the results to the former temporary vector.Figure 2
Strategy of partitioning, storing a block sparse matrix, and performing BSpMV on a multi-GPU system based on PETSc.On a multi-GPU system, we consider using a MPI process to control a GPU card. In each MPI process, the arrays expressing the MDP and ODP in the BCSR format for the partitioned matrix are copied to the corresponding GPU memories whose data pointers are encapsulated in aGMRESInfo instance. This makes the MDP and ODP capable of conducting the local BSpMV on GPUs by using cusparseDbsrmv from the cuSPARSE library [28]. However, the multiplication between the ODP and the received vector cannot be proceeded until the MPI communication has completed. In our implementation, we perform the MPI communication from GPUs to GPUs directly using the GPU-Direct technique. Note that although in the latest version of PETSc, the GPU-Direct communication is supported in the VecScatter object, the VecScatterBegin, and VecScatterEnd functions that manage the data exchange between vectors (with either ViennaCL or CUDA type in PETSc) cannot be applied directly because the GPU vectors in our implementation are designed in the GMRESInfo structure. Efforts are made to extract the low-level communication map from the VectorScatter object that is used for distributed BSpMV on CPUs, and additional kernels are developed. The BSpMV with GPU-Direct communication is summarized in Algorithm 2.Algorithm 2:BSpMV on multi-GPUs with GPU-Direct communication.
(1)
SetGPUSendbuffer(…); ⟵ from the local vector to theith GPU sending buffer.(2)
cusparseDbsrmv(…); ⟵ MDP BSpMV on theith GPU in a CUDA stream.(3)
SpMVGPUDirectComm(…); ⟵ exchange data between GPUs and wait until completed.(4)
cusparseDbsrmv(…); ⟵ ODP BSpMV on theith GPU.According to the indices in the VecScatter object that record the partial data required to be sent to remote GPUs, the GPU kernel SetGPUSendbuﬀer is developed to move the data from the local vector to the GPU buﬀer. cusparseDbsrmv is called asynchronously to perform the MDP BSpMV in a separate CUDA stream, which makes the computation overlap with the following MPI communication. Then, the CUDA-aware MPI can automatically identify all GPU buﬀers that are requested to send and receive and transfer data between GPUs directly. Once the data is received successfully in each GPU’s memory, the BSpMV can be conducted between the ODP and the received buﬀer on GPUs.The vector operations in Algorithm1 including vector inner products (lines 3 and 12), scaling (lines 3 and 17), and linear operations (line 21) are implemented using efficient subroutines from the cuBLAS [34] library. When these functions are called, scalar values are passed by reference on the host instead of the device to receive the computing results from GPUs because some of the results are used by the calculations on CPUs. For example, lines 13–16 and 20 are assigned to CPUs due to very light workload, and they need normw and H for computing and logic judgements on CPUs. Therefore, passing host reference in the call of cuBLAS functions avoids manually copying data back to the host memories.The vector operations in lines 9-10 are time consuming because they are in a nested loop. For eachi0≤i<m, they require i+1 times of MPIAllreduce communication in total because each computation of Hk,i requires one MPIAllreduce call and Hk,i has to be used to update w⟶ immediately. In fact, since v⟶k0≤k≤i are orthogonal with each other, PETSc splits the single loop into two separate ones each of which is applied to one line of lines 9-10. Therefore, the communication in obtaining all Hk,i can be reduced by launching only one call of MPIAllreduce. Furthermore, PETSc optimized the loop of computing w⟶,v⟶k0≤k≤i by using batch computing. For example, the method conducts 8 inner products from w⟶,v⟶0 to w⟶,v⟶7 by launching two CUDA kernels each of which performs 4 inner products in a batch. In each batch, w⟶ is fetched from global memory only once and reused for all four inner products. This technique reduces the visiting of global memory for w⟶ repeatedly. Another method where the inner product is transformed to the multiplication between a dense matrix (constructed by all v⟶k) and a vector w⟶ was studied in [4]. However, this method is not covered in this paper since all of our implementations and optimizations are based on PETSc’s framework.The classical inner product on GPUs [35] employed by PETSc requires the host to copy the reductions of thread blocks from the GPU and conduct a final reduction over all thread blocks on the CPU. In the case of batch computing, the memory copy from the GPU to the host and the ﬁnal reduction on the host for one batch can be overlapped with the GPU computing for another batch. Figure 3(b) illustrates the idea of asynchronization on diﬀerent batches, whereas Figure 3(a) shows the synchronized case. To implement this idea, batches are assigned to two CUDA streams alternately so that the overhead of memory copies and reductions on the CPU can be reduced.Figure 3
Asynchronous batch computing for inner products on GPUs.
(a)(b)The method of batch computing can also be applied to the updating ofw⟶ in the loop. In PETSc’s implementation, it calls cublasDaxpy one by one to update w⟶ for the CUDA vector type or updates w⟶ with a batch size of 2 for the ViennaCL [36] vector type. We develop several CUDA kernels to use a large batch size instead of calling cublasDaxpy in our implementation. The template of the kernel is shown in Algorithm 3. VecAXPYKernelj is developed to compute a linear combination of xj using coefficients aj. Compared to visiting w⟶ once using cublasDaxpy at a time, calling VecAXPYKernelj in a batch reduces the number of times of fetching w⟶ from the global memory. For a given number of loops, the kernels need to be called in a sequential order since w⟶ that is updated in one batch has to be used in the following batch.Algorithm 3:__global__ void VecAXPYKernelj (w, len, a1, ..., aj , x1, x2, ..., xj).
(1)
inttid = threadIdx.x + blockDim.x×blockIdx.x;(2)
whiletid<lendo(3)
wtid+=a1x1tid+⋯+aixjtid;(4)
tid+=gridDim.x×blockDim.x;(5)
end


### 3.3. Preconditioning

To implement equations (5) and (6) on a parallel system with distributed GPUs, we consider the block Jacobi preconditioner as(7)ABi≈LBiUBi=MBi,i=0,1,2,…,l−1,(8)MBi−1v⟶i=z⟶i⟶LBif⟶i=v⟶i,UBiz⟶i=f⟶i,i=0,1,2,…,l−1,where ABi and MBi represents the main diagonal part of the partitioned matrix and the preconditioning matrix in the ith process, respectively, z⟶i, f⟶i, and v⟶i are the partitioned vectors in the ith process, and l is the total number of processes.For pointwise matrices, the level-scheduling-based parallel methods [2, 21] for the preconditioning step depend heavily on the matrix pattern. Therefore, some cases failed to gain speedups on GPUs over the CPU implementation. To keep a flexible choice of computing the preconditioning step, we investigate three implementations of computing equation (8) for block sparse matrix based on PETSc in this section and make performance comparisons in Section 5.PETSc implements several efficient subroutines to conduct the ILU factorizations and preconditioning operations for block sparse matrices with various block sizes. In the first implementation, we consider performing (8) using PETSc’s API on CPUs and copying the preconditioned vector to the device memory at every iteration of Algorithm 1. The developed function is listed in Algorithm 4. The function accepts four input variables as ksp, vecv, vecz, and dv, where ksp is an instance of the PETSc’s KSP structure, vecv and vecz, encapsulated by the PETSc’s Vec structure in the host memory, are the left- and right-side vectors in equation (8), and dv is the device pointer to the left side vector on the GPU. The function outputs the pointer to the right-side vector. To obtain an explicit pointer to the low-level array data of the Vec object in PETSc, we employ two APIs, VecGetArray and VecRetoreArray, before and after using the object. Then, the preconditioning procedure can be divided into three steps. The data transfer is firstly launched to copy the left-side vector from the device to host using cudaMempcy from the CUDA library. This is followed by a call of PCApply in PETSc to apply the declared preconditioner to vecv. Inside PCApply, it performs two block sparse triangular solves in which the calculations are conducted in blockwise. Figure 4 shows the pattern of solving a lower block sparse triangular system with block size s in PETSc. The preconditioned vector, vecz, is then copied to the device space which dz points to.Algorithm 4:PETScPrecond(in: ksp, vecv, vecz, dv, out: dz).
(1)
VecGetArray(vecv, & hv);(2)
cudaMemcpy(hv, dv, sizeof(vecv), cudaMemcpyDeviceToHost);(3)
VecRestoreArray(vecv, & hv);(4)
PCApply(ksp ⟶ pc, vecv, vecz);(5)
VecGetArray(vecz, & hz);(6)
cudaMemcpy(dz, hz, sizeof(vecz), cudaMemcpyHostToDevice);(7)
VecRestoreArray(vecz, & hz);Figure 4
Pattern of the lower block sparse triangular solver in PETSc.In the second implementation, we utilize the efficient triangular solvers in the BCSR format provided by the cuSPARSE library. Algorithm5 lists the main steps of the function (cuSPARSEPrecond) which we develop and integrate into PETSc. The first step (lines 2–7), also called the preprocessing phase, is to estimate the memory requirement, allocate adequate spaces and extract possible parallelism for the subsequent solve phase. The preprocessing step is required to be executed only once on the local GPUs because the lower LBi and upper UBi factors remain unchanged during the iterative process of GMRES once the factors has been constructed by the ILU factorization. The second step is to perform two block sparse triangular systems by calling cusparseDbsrsv2_solve twice with the parameters of LBi and UBi.Algorithm 5:cuSPARSEPrecond(in: dv, df, firstcall, out: dz).
(1)
iffirstcallthen(2)
cusparseDbsrsv2_bufferSize(…,Ldescr, Ablkval, Arowptr, Acolval, s, …, & Lbufsz);(3)
cusparseDbsrsv2_bufferSize(…,Udescr, Ablkval, Arowptr, Acolval, s, …, & U bufsz);(4)
cudaMalloc(&pLbuf, Lbufsz);(5)
cudaMalloc(&pUbuf, Ubufsz);(6)
cusparseDbsrsv2_analysis(…,Ldescr, Ablkval, Arowptr, Acolval, s, …, pLbuf);(7)
cusparseDbsrsv2_analysis(…,Udescr, Ablkval, Arowptr, Acolval, s, …, pUbuf);(8)
end(9)
cusparseDbsrsv2_solve(…,Ldescr, Ablkval, Arowptr, Acolval, s, …, dv, df, …, pLbuf);(10)
cusparseDbsrsv2_solve(…,Udescr, Ablkval, Arowptr, Acolval, s, …, df, dz, …, pUbuf);Since the ILU factorization ofABi equation (7) is conducted outside Algorithm 1, the structure compatibility between the ILU factorization and triangular solves has to be studied. Specifically, when the block ILU factorization using PETSc on CPUs is preferred, the solving of the upper triangular system has to be transformed into(9)U˜Biz⟶i=diagUBi−1f⟶i,where U˜Bi=diagUBi−1UBi. This is because cusparseDbsrsv2_solve ignores all the values in the lower part of the input matrix according to the pointwise diagonal. And, UBiz⟶i=f⟶i is equivalent to the modified system above with all identity matrices in the diagonal blocks. The process can be illustrated by Figure 5(a).Figure 5
Two different patterns of the upper block sparse triangular solver using cusparseDbsrsv2_solve: (a)U˜Bi; (b) the upper part of (b) ABi after the cuSPARSE’s ILU factorization.
(a)(b)By contrast, the block-based ILU factorization using cuSPARSE can acceptABi in the BCSR format as the input matrix and replace ABi with the lower and upper factors. In this case, it conducts a pointwise ILU factorization by treating s2 values in each block as independent scalars, even though ABi is expressed in the BCSR format. Therefore, the replaced ABi can be directly used as the input matrix for the two calls of cusparseDbsrsv2_solve (lines 9-10 in Algorithm 5) where the pointwise lower and upper factors can be automatically extracted. This strategy is shown in Figure 5(b).The asynchronous techniques [22, 23, 25–27] where the ILU factorizations and triangular solves are calculated iteratively and inexactly have been attractive. In these methods, although a linear system may need more iterations to converge, the overall performance can be improved because the overhead of the preconditioning step at every step is reduced. In the last implementation, we investigate the inexact triangular solves for block sparse systems on a multi-GPU system.In theith process, by considering the upper block triangular system UBiz⟶i=f⟶i, we explicitly express the block vector z⟶ji in an iterative way as(10)z⟶jiit+1=UBij,j−1f⟶ji−∑k>jUBij,k−1z⟶kiit,k,j∈Sp,j=r−1,r−2,…,0,where z⟶jiit+1 represents the jth block vector at the it+1th iteration. And the lower block triangular system can be written in a similar way. By removing the data dependencies in the r constraints, we can update z⟶ji iteratively without the synchronization of data to approximate the exact solutions. Compared to the scalar version for pointwise matrices, equation (10) contains totally different computations. The scalar multiplication is changed to the small dense matrixs×s vector multiplication, and the inverse of a scalar is replaced with the inverse of a s×s matrix.It is straightforward to implement the scalar version of equation (10) on GPUs since each thread can be assigned to the calculation of a scalar solution. However, assigning the computation of z⟶ji to a single CUDA thread is not preferred on GPUs because the global memory accesses are far from coalesced when a dense block matrix multiplies a block vector [29]. Inspired by the fine-grained idea for the block sparse matrix-vector multiplication on GPUs [29], we develop a CUDA kernel shown in Algorithm 6 to update z⟶ji in an iteration.Algorithm 6:KernelBtrisvUpper(urowptr, ucolval, ublkval, df, dz, r, stride).
(1)
tid=blockDim.x×blockIdx.x+threadIdx.x;(2)
__shared__ doubleshwarps_per_blockwarp_size;(3)
__shared__ doubledgwarps_per_blocks2;(4)
irow_forward=tid/warp_size;(5)
ifirow_forward<rthen(6)
irow=r−1−irow_forward;(7)
sidx=rowptrirow; eidx=rowptrirow+1;(8)
rpos=threadIdx.x%warp_size; warp_id=threadIdx.x/warp_size;(9)
cl=rpos/s%s; rw=rpos%s;(10)
range=warp_size/s2×s2;(11)
idx=sidx+1+rpos/s2;(12)
sum=0.0; tmp=0.0; shwarp_idrpos=0.0;(13)
ifrpos<rangethen(14)
ifrpos<s2then(15)
dgwarp idrpos=blkvalsidx×s2+rpos;(16)
end(17)
whileidx<eidxdo(18)
col=col_validx;(19)
sum=sum−blkvalidx×s2+cl×s+rw×zcol×s+cl;(20)
idx=idx+warp_size/s2;(21)
end(22)
swarp_idrpos=sum;(23)
forstep≥1;step/=2do(24)
ifrpos<step×sand rpos+step×s<warp_sizethen(25)
shwarp_idrpos+=shwarp_idrpos+step×s;(26)
end(27)
end(28)
ifrpos<sthen(29)
tmp=0.0;shwarp_idrpos+=dfirow×s+rpos;(30)
fori=0;i<s;i++do(31)
tmp+=dgwarp_idrpos+i×s×shwarp_idi;(32)
end(33)
dzirow×s+rpos=tmp;(34)
end(35)
end(36)
endAs illustrated in Figure6, we use 32 threads within a warp to accumulate UBij,kz⟶kiitover all blocks at the jth row of UBi that satisfy k>j and j,k∈Sp. Specifically, each thread in the warp computes the multiplication of only one element in UBij,k and the corresponding element in the block vector z⟶kiit (lines 18–20 in Algorithm 6) and stores the scalar result in a register variable. However, 32 threads in a warp may be unable to cover all blocks simultaneously; then, each thread accumulates the results one by one via performing several cycles on the remaining blocks. To make the thread private variable visible to each other, the threads in each warp share a space which is declared as the type of shared memory (line 2 in Algorithm 6). After that, the scalar result in each thread is copied into the shared memory, and a serial of reduction operations (lines 23–27 in Algorithm 6) are launched to summarize the results for the corresponding rows of the resulting block vector of ∑k>jUBij,kz⟶kiit. This is followed by the operation of vector addition on f⟶iand sh with the first s elements. The last task of a warp is to conduct a matrix-vector multiplication using the inverse of UBij,j and the resulting vector in the shared memory from the last step. The inverse of UBij,j can be obtained from the block ILU factorization in PETSc because PETSc replaces all diagonal blocks of UBi with their inversions to prepare for the subsequent block triangular solves. We take advantage of this fact by copying the lower and upper factors in PETSc from the host to device before the kernel is executed. And, the first s2 threads in a warp load the diagonal block into the shared memory for the matrix-vector multiplication at line 31 in Algorithm 6.Figure 6
Illustration of the CUDA kernel.For the lower block triangular system, the iterative process is written as(11)f⟶jiit+1=v⟶ji−∑k<jLBij,kf⟶kiit,k,j∈Sp,j=0,1,2,…,r−1,and the kernel to iterate equation (11) for one step is developed in a similar way. Slightly different from the upper one, the kernel maps r rows of LBi from top to bottom. In addition, the matrix-vector multiplication can be discarded because the diagonal blocks of LBi are set to identity matrices in PETSc. We summarize the preconditioning step in Algorithm 7 by performing the two kernels several times.Algorithm 7:IterativePrecond(in: lrowptr, lcolval, lblkval, urowptr, ucolval, ublkval, dv, df, r, stride, itmax out: dz).
(1)
for it = 0; it<itmax; it++ do(2)
KernelBtrisvLower<<< nblocksl, nthreadsl >>> (lrowptr, lcolval, lblkval, dv, df, r, stride);(3)
end(4)
for it = 0;it<itmax; it + + do(5)
KernelBtrisvUpper<<< nblocksu, nthreadsu >>> (urowptr, ucolval, ublkval, df, dz, r, stide);(6)
end


## 3.1. Chart of Development

PETSc [31] offers a wide range of Application Programming Interfaces (APIs) for users to manipulate low-level data structures. Furthermore, it uses function pointers in standard C language instead of an Object-Oriented language such as C++ to achieve data encapsulation and polymorphism [33]. We take advantage of function pointers and structures in C language to develop a user friendly multi-GPU version of the preconditioned GMRES algorithm.The chart of our code development is shown in Figure1. Data structure GMRESInfo is designed to store all necessary data ﬁelds on the GPU. The memory spaces of the fields are allocated when an instance of GMRESInfo is created. Before the GMRES solver is executed, the fields in GMRESInfo are initialized on GPUs by launching the memory copy operations provided by the basic CUDA library from the host to device. The encapsulated data on the host that needs to be transferred to GPUs can be accessed using the low-level APIs in PETSc. To make full use of the existing numerical toolkits, we employ the efficient functions in cuBLAS and cuSPARSE libraries to perform the local sparse matrix and vector operations in Algorithm 1. On a multi-GPU platform, a portion of the computation inevitably requires data exchange among GPUs. This is implemented and optimized using the CUDA-aware MPI library that supports the GPU-Direct technology and transfers data directly among GPU memories through network adapters.Figure 1
The development chart of the multi-GPU version of GMRES algorithm based on PETSc.To adopt the multi-GPU version of GMRES solver in an application, a user needs to call three additional functions. Two of them creates and destroys aGMRESInfo instance, respectively, while the other (SetUserDeﬁnedGMRES) sets the function pointer of GMRES to our developed procedure.To implement Algorithm1 on a multi-GPU system, the workload for the host and device should be proper assigned. The matrix-vector operations (line 3 and line 7) and the basic vector operations such as inner product (line 9 and line 12) and the linear combinations (line 10 and line 17) are of high parallelism, so they are fully performed on the GPUs. To conduct convergence monitoring, the conditional statements (line 14–16) are inevitable. And, this part is executed on the host. Because cuSPARSE library allows users to pass a reference on the host when an API is executed on the GPU, the scalar parameter normw is declared on the host to receive the result of w⟶2 and can be directly used in line 14. The size of the Hessenberg matrix is (m + 1) × m, where m (m is set to 30 in PETSc by default) is usually much smaller than the matrix size, so the updating of Hessenberg matrix (line 13) and solving of the least square problem (line 20) are performed on the host.To distributedly store the data on a multi-GPU system, block matrices and vectors are partitioned by rows. On a GPU cluster withN GPUs, the ith GPU stores ri block rows of AB which is represented in BCSR format. In addition, there are several vectors that need to be allocated on the ith GPU. r⟶0i, z⟶i, x⟶i, and w⟶i are allocated with the size of ri×s, and Vt stores up to m vectors each of which is of size ri×s. In the preconditioning step (line 6), the preconditioning matrix on the ith GPU, denoted as MBi, is stored in separate factors, i.e., LBi and UBi, and the size of memory requirements for ILU(0) level equals to the size of ABi, where ABi is the main diagonal part of the partitioned matrix on the ith GPU.


## 3.2. Implementation and Optimization

The major computation in Algorithm1 consists of block sparse matrix-vector multiplication (BSpMV), the preconditioning step, and the vector operations. For a block sparse matrix on a distributed system, PETSc partitions the matrix by consecutive block rows and stores each partition in two separate structures. One structure is for the main diagonal of the partitioned matrix, the other is for the oﬀ diagonals. Figure 2 shows the strategy of partitioning and storing a 18s×18s (s is the block size) block sparse matrix obtained by the 5-point finite diﬀerencing scheme in PETSc. Each partition is responsible for storing 6 block rows and assigned to a MPI process. Take the second partition (dashed box in red) as an example, PETSc stores the 6s×6s square block sparse matrix in the main diagonal as the BCSR format by renumbering the block rows and columns starting from zero. The two oﬀ-diagonal parts are merged and renumbered as a local 6s×12s block sparse matrix in the BCSR format. The main advantage of this strategy is that it separates the main diagonal part (MDP, communication-free) from the oﬀ-diagonal part (ODP, communication-required), which can be used for overlapping the computation and communication. Therefore, each partition needs three steps to perform the multiplication between a block sparse matrix and a vector. The MDP firstly multiplies the local part of the vector and stores the result temporarily. All partitions then launch MPI sending and receiving calls to exchange data of the vector globally. Lastly, the ODP in each partition multiplies the received vector and accumulates the results to the former temporary vector.Figure 2
Strategy of partitioning, storing a block sparse matrix, and performing BSpMV on a multi-GPU system based on PETSc.On a multi-GPU system, we consider using a MPI process to control a GPU card. In each MPI process, the arrays expressing the MDP and ODP in the BCSR format for the partitioned matrix are copied to the corresponding GPU memories whose data pointers are encapsulated in aGMRESInfo instance. This makes the MDP and ODP capable of conducting the local BSpMV on GPUs by using cusparseDbsrmv from the cuSPARSE library [28]. However, the multiplication between the ODP and the received vector cannot be proceeded until the MPI communication has completed. In our implementation, we perform the MPI communication from GPUs to GPUs directly using the GPU-Direct technique. Note that although in the latest version of PETSc, the GPU-Direct communication is supported in the VecScatter object, the VecScatterBegin, and VecScatterEnd functions that manage the data exchange between vectors (with either ViennaCL or CUDA type in PETSc) cannot be applied directly because the GPU vectors in our implementation are designed in the GMRESInfo structure. Efforts are made to extract the low-level communication map from the VectorScatter object that is used for distributed BSpMV on CPUs, and additional kernels are developed. The BSpMV with GPU-Direct communication is summarized in Algorithm 2.Algorithm 2:BSpMV on multi-GPUs with GPU-Direct communication.
(1)
SetGPUSendbuffer(…); ⟵ from the local vector to theith GPU sending buffer.(2)
cusparseDbsrmv(…); ⟵ MDP BSpMV on theith GPU in a CUDA stream.(3)
SpMVGPUDirectComm(…); ⟵ exchange data between GPUs and wait until completed.(4)
cusparseDbsrmv(…); ⟵ ODP BSpMV on theith GPU.According to the indices in the VecScatter object that record the partial data required to be sent to remote GPUs, the GPU kernel SetGPUSendbuﬀer is developed to move the data from the local vector to the GPU buﬀer. cusparseDbsrmv is called asynchronously to perform the MDP BSpMV in a separate CUDA stream, which makes the computation overlap with the following MPI communication. Then, the CUDA-aware MPI can automatically identify all GPU buﬀers that are requested to send and receive and transfer data between GPUs directly. Once the data is received successfully in each GPU’s memory, the BSpMV can be conducted between the ODP and the received buﬀer on GPUs.The vector operations in Algorithm1 including vector inner products (lines 3 and 12), scaling (lines 3 and 17), and linear operations (line 21) are implemented using efficient subroutines from the cuBLAS [34] library. When these functions are called, scalar values are passed by reference on the host instead of the device to receive the computing results from GPUs because some of the results are used by the calculations on CPUs. For example, lines 13–16 and 20 are assigned to CPUs due to very light workload, and they need normw and H for computing and logic judgements on CPUs. Therefore, passing host reference in the call of cuBLAS functions avoids manually copying data back to the host memories.The vector operations in lines 9-10 are time consuming because they are in a nested loop. For eachi0≤i<m, they require i+1 times of MPIAllreduce communication in total because each computation of Hk,i requires one MPIAllreduce call and Hk,i has to be used to update w⟶ immediately. In fact, since v⟶k0≤k≤i are orthogonal with each other, PETSc splits the single loop into two separate ones each of which is applied to one line of lines 9-10. Therefore, the communication in obtaining all Hk,i can be reduced by launching only one call of MPIAllreduce. Furthermore, PETSc optimized the loop of computing w⟶,v⟶k0≤k≤i by using batch computing. For example, the method conducts 8 inner products from w⟶,v⟶0 to w⟶,v⟶7 by launching two CUDA kernels each of which performs 4 inner products in a batch. In each batch, w⟶ is fetched from global memory only once and reused for all four inner products. This technique reduces the visiting of global memory for w⟶ repeatedly. Another method where the inner product is transformed to the multiplication between a dense matrix (constructed by all v⟶k) and a vector w⟶ was studied in [4]. However, this method is not covered in this paper since all of our implementations and optimizations are based on PETSc’s framework.The classical inner product on GPUs [35] employed by PETSc requires the host to copy the reductions of thread blocks from the GPU and conduct a final reduction over all thread blocks on the CPU. In the case of batch computing, the memory copy from the GPU to the host and the ﬁnal reduction on the host for one batch can be overlapped with the GPU computing for another batch. Figure 3(b) illustrates the idea of asynchronization on diﬀerent batches, whereas Figure 3(a) shows the synchronized case. To implement this idea, batches are assigned to two CUDA streams alternately so that the overhead of memory copies and reductions on the CPU can be reduced.Figure 3
Asynchronous batch computing for inner products on GPUs.
(a)(b)The method of batch computing can also be applied to the updating ofw⟶ in the loop. In PETSc’s implementation, it calls cublasDaxpy one by one to update w⟶ for the CUDA vector type or updates w⟶ with a batch size of 2 for the ViennaCL [36] vector type. We develop several CUDA kernels to use a large batch size instead of calling cublasDaxpy in our implementation. The template of the kernel is shown in Algorithm 3. VecAXPYKernelj is developed to compute a linear combination of xj using coefficients aj. Compared to visiting w⟶ once using cublasDaxpy at a time, calling VecAXPYKernelj in a batch reduces the number of times of fetching w⟶ from the global memory. For a given number of loops, the kernels need to be called in a sequential order since w⟶ that is updated in one batch has to be used in the following batch.Algorithm 3:__global__ void VecAXPYKernelj (w, len, a1, ..., aj , x1, x2, ..., xj).
(1)
inttid = threadIdx.x + blockDim.x×blockIdx.x;(2)
whiletid<lendo(3)
wtid+=a1x1tid+⋯+aixjtid;(4)
tid+=gridDim.x×blockDim.x;(5)
end


## 3.3. Preconditioning

To implement equations (5) and (6) on a parallel system with distributed GPUs, we consider the block Jacobi preconditioner as(7)ABi≈LBiUBi=MBi,i=0,1,2,…,l−1,(8)MBi−1v⟶i=z⟶i⟶LBif⟶i=v⟶i,UBiz⟶i=f⟶i,i=0,1,2,…,l−1,where ABi and MBi represents the main diagonal part of the partitioned matrix and the preconditioning matrix in the ith process, respectively, z⟶i, f⟶i, and v⟶i are the partitioned vectors in the ith process, and l is the total number of processes.For pointwise matrices, the level-scheduling-based parallel methods [2, 21] for the preconditioning step depend heavily on the matrix pattern. Therefore, some cases failed to gain speedups on GPUs over the CPU implementation. To keep a flexible choice of computing the preconditioning step, we investigate three implementations of computing equation (8) for block sparse matrix based on PETSc in this section and make performance comparisons in Section 5.PETSc implements several efficient subroutines to conduct the ILU factorizations and preconditioning operations for block sparse matrices with various block sizes. In the first implementation, we consider performing (8) using PETSc’s API on CPUs and copying the preconditioned vector to the device memory at every iteration of Algorithm 1. The developed function is listed in Algorithm 4. The function accepts four input variables as ksp, vecv, vecz, and dv, where ksp is an instance of the PETSc’s KSP structure, vecv and vecz, encapsulated by the PETSc’s Vec structure in the host memory, are the left- and right-side vectors in equation (8), and dv is the device pointer to the left side vector on the GPU. The function outputs the pointer to the right-side vector. To obtain an explicit pointer to the low-level array data of the Vec object in PETSc, we employ two APIs, VecGetArray and VecRetoreArray, before and after using the object. Then, the preconditioning procedure can be divided into three steps. The data transfer is firstly launched to copy the left-side vector from the device to host using cudaMempcy from the CUDA library. This is followed by a call of PCApply in PETSc to apply the declared preconditioner to vecv. Inside PCApply, it performs two block sparse triangular solves in which the calculations are conducted in blockwise. Figure 4 shows the pattern of solving a lower block sparse triangular system with block size s in PETSc. The preconditioned vector, vecz, is then copied to the device space which dz points to.Algorithm 4:PETScPrecond(in: ksp, vecv, vecz, dv, out: dz).
(1)
VecGetArray(vecv, & hv);(2)
cudaMemcpy(hv, dv, sizeof(vecv), cudaMemcpyDeviceToHost);(3)
VecRestoreArray(vecv, & hv);(4)
PCApply(ksp ⟶ pc, vecv, vecz);(5)
VecGetArray(vecz, & hz);(6)
cudaMemcpy(dz, hz, sizeof(vecz), cudaMemcpyHostToDevice);(7)
VecRestoreArray(vecz, & hz);Figure 4
Pattern of the lower block sparse triangular solver in PETSc.In the second implementation, we utilize the efficient triangular solvers in the BCSR format provided by the cuSPARSE library. Algorithm5 lists the main steps of the function (cuSPARSEPrecond) which we develop and integrate into PETSc. The first step (lines 2–7), also called the preprocessing phase, is to estimate the memory requirement, allocate adequate spaces and extract possible parallelism for the subsequent solve phase. The preprocessing step is required to be executed only once on the local GPUs because the lower LBi and upper UBi factors remain unchanged during the iterative process of GMRES once the factors has been constructed by the ILU factorization. The second step is to perform two block sparse triangular systems by calling cusparseDbsrsv2_solve twice with the parameters of LBi and UBi.Algorithm 5:cuSPARSEPrecond(in: dv, df, firstcall, out: dz).
(1)
iffirstcallthen(2)
cusparseDbsrsv2_bufferSize(…,Ldescr, Ablkval, Arowptr, Acolval, s, …, & Lbufsz);(3)
cusparseDbsrsv2_bufferSize(…,Udescr, Ablkval, Arowptr, Acolval, s, …, & U bufsz);(4)
cudaMalloc(&pLbuf, Lbufsz);(5)
cudaMalloc(&pUbuf, Ubufsz);(6)
cusparseDbsrsv2_analysis(…,Ldescr, Ablkval, Arowptr, Acolval, s, …, pLbuf);(7)
cusparseDbsrsv2_analysis(…,Udescr, Ablkval, Arowptr, Acolval, s, …, pUbuf);(8)
end(9)
cusparseDbsrsv2_solve(…,Ldescr, Ablkval, Arowptr, Acolval, s, …, dv, df, …, pLbuf);(10)
cusparseDbsrsv2_solve(…,Udescr, Ablkval, Arowptr, Acolval, s, …, df, dz, …, pUbuf);Since the ILU factorization ofABi equation (7) is conducted outside Algorithm 1, the structure compatibility between the ILU factorization and triangular solves has to be studied. Specifically, when the block ILU factorization using PETSc on CPUs is preferred, the solving of the upper triangular system has to be transformed into(9)U˜Biz⟶i=diagUBi−1f⟶i,where U˜Bi=diagUBi−1UBi. This is because cusparseDbsrsv2_solve ignores all the values in the lower part of the input matrix according to the pointwise diagonal. And, UBiz⟶i=f⟶i is equivalent to the modified system above with all identity matrices in the diagonal blocks. The process can be illustrated by Figure 5(a).Figure 5
Two different patterns of the upper block sparse triangular solver using cusparseDbsrsv2_solve: (a)U˜Bi; (b) the upper part of (b) ABi after the cuSPARSE’s ILU factorization.
(a)(b)By contrast, the block-based ILU factorization using cuSPARSE can acceptABi in the BCSR format as the input matrix and replace ABi with the lower and upper factors. In this case, it conducts a pointwise ILU factorization by treating s2 values in each block as independent scalars, even though ABi is expressed in the BCSR format. Therefore, the replaced ABi can be directly used as the input matrix for the two calls of cusparseDbsrsv2_solve (lines 9-10 in Algorithm 5) where the pointwise lower and upper factors can be automatically extracted. This strategy is shown in Figure 5(b).The asynchronous techniques [22, 23, 25–27] where the ILU factorizations and triangular solves are calculated iteratively and inexactly have been attractive. In these methods, although a linear system may need more iterations to converge, the overall performance can be improved because the overhead of the preconditioning step at every step is reduced. In the last implementation, we investigate the inexact triangular solves for block sparse systems on a multi-GPU system.In theith process, by considering the upper block triangular system UBiz⟶i=f⟶i, we explicitly express the block vector z⟶ji in an iterative way as(10)z⟶jiit+1=UBij,j−1f⟶ji−∑k>jUBij,k−1z⟶kiit,k,j∈Sp,j=r−1,r−2,…,0,where z⟶jiit+1 represents the jth block vector at the it+1th iteration. And the lower block triangular system can be written in a similar way. By removing the data dependencies in the r constraints, we can update z⟶ji iteratively without the synchronization of data to approximate the exact solutions. Compared to the scalar version for pointwise matrices, equation (10) contains totally different computations. The scalar multiplication is changed to the small dense matrixs×s vector multiplication, and the inverse of a scalar is replaced with the inverse of a s×s matrix.It is straightforward to implement the scalar version of equation (10) on GPUs since each thread can be assigned to the calculation of a scalar solution. However, assigning the computation of z⟶ji to a single CUDA thread is not preferred on GPUs because the global memory accesses are far from coalesced when a dense block matrix multiplies a block vector [29]. Inspired by the fine-grained idea for the block sparse matrix-vector multiplication on GPUs [29], we develop a CUDA kernel shown in Algorithm 6 to update z⟶ji in an iteration.Algorithm 6:KernelBtrisvUpper(urowptr, ucolval, ublkval, df, dz, r, stride).
(1)
tid=blockDim.x×blockIdx.x+threadIdx.x;(2)
__shared__ doubleshwarps_per_blockwarp_size;(3)
__shared__ doubledgwarps_per_blocks2;(4)
irow_forward=tid/warp_size;(5)
ifirow_forward<rthen(6)
irow=r−1−irow_forward;(7)
sidx=rowptrirow; eidx=rowptrirow+1;(8)
rpos=threadIdx.x%warp_size; warp_id=threadIdx.x/warp_size;(9)
cl=rpos/s%s; rw=rpos%s;(10)
range=warp_size/s2×s2;(11)
idx=sidx+1+rpos/s2;(12)
sum=0.0; tmp=0.0; shwarp_idrpos=0.0;(13)
ifrpos<rangethen(14)
ifrpos<s2then(15)
dgwarp idrpos=blkvalsidx×s2+rpos;(16)
end(17)
whileidx<eidxdo(18)
col=col_validx;(19)
sum=sum−blkvalidx×s2+cl×s+rw×zcol×s+cl;(20)
idx=idx+warp_size/s2;(21)
end(22)
swarp_idrpos=sum;(23)
forstep≥1;step/=2do(24)
ifrpos<step×sand rpos+step×s<warp_sizethen(25)
shwarp_idrpos+=shwarp_idrpos+step×s;(26)
end(27)
end(28)
ifrpos<sthen(29)
tmp=0.0;shwarp_idrpos+=dfirow×s+rpos;(30)
fori=0;i<s;i++do(31)
tmp+=dgwarp_idrpos+i×s×shwarp_idi;(32)
end(33)
dzirow×s+rpos=tmp;(34)
end(35)
end(36)
endAs illustrated in Figure6, we use 32 threads within a warp to accumulate UBij,kz⟶kiitover all blocks at the jth row of UBi that satisfy k>j and j,k∈Sp. Specifically, each thread in the warp computes the multiplication of only one element in UBij,k and the corresponding element in the block vector z⟶kiit (lines 18–20 in Algorithm 6) and stores the scalar result in a register variable. However, 32 threads in a warp may be unable to cover all blocks simultaneously; then, each thread accumulates the results one by one via performing several cycles on the remaining blocks. To make the thread private variable visible to each other, the threads in each warp share a space which is declared as the type of shared memory (line 2 in Algorithm 6). After that, the scalar result in each thread is copied into the shared memory, and a serial of reduction operations (lines 23–27 in Algorithm 6) are launched to summarize the results for the corresponding rows of the resulting block vector of ∑k>jUBij,kz⟶kiit. This is followed by the operation of vector addition on f⟶iand sh with the first s elements. The last task of a warp is to conduct a matrix-vector multiplication using the inverse of UBij,j and the resulting vector in the shared memory from the last step. The inverse of UBij,j can be obtained from the block ILU factorization in PETSc because PETSc replaces all diagonal blocks of UBi with their inversions to prepare for the subsequent block triangular solves. We take advantage of this fact by copying the lower and upper factors in PETSc from the host to device before the kernel is executed. And, the first s2 threads in a warp load the diagonal block into the shared memory for the matrix-vector multiplication at line 31 in Algorithm 6.Figure 6
Illustration of the CUDA kernel.For the lower block triangular system, the iterative process is written as(11)f⟶jiit+1=v⟶ji−∑k<jLBij,kf⟶kiit,k,j∈Sp,j=0,1,2,…,r−1,and the kernel to iterate equation (11) for one step is developed in a similar way. Slightly different from the upper one, the kernel maps r rows of LBi from top to bottom. In addition, the matrix-vector multiplication can be discarded because the diagonal blocks of LBi are set to identity matrices in PETSc. We summarize the preconditioning step in Algorithm 7 by performing the two kernels several times.Algorithm 7:IterativePrecond(in: lrowptr, lcolval, lblkval, urowptr, ucolval, ublkval, dv, df, r, stride, itmax out: dz).
(1)
for it = 0; it<itmax; it++ do(2)
KernelBtrisvLower<<< nblocksl, nthreadsl >>> (lrowptr, lcolval, lblkval, dv, df, r, stride);(3)
end(4)
for it = 0;it<itmax; it + + do(5)
KernelBtrisvUpper<<< nblocksu, nthreadsu >>> (urowptr, ucolval, ublkval, df, dz, r, stide);(6)
end


## 4. Experiments

### 4.1. Test Platform

In our experiments, the tests were conducted on a heterogeneous cluster at Computer Network Information Center, Chinese Academy of Sciences. The cluster consists of 270 blades nodes, 30 GPGPU nodes, and 40 MIC nodes. Each GPGPU node is configured with 2 Intel E5-2680 V2 (Ivy Bridge | 10C | 2.8 GHz) CPUs, 128 GB DDR3 ECC 1866 MHz memory, and two NVIDIA Tesla K20 GPGPU cards. Each card is configured by 320 bit GDDR5 5 GB memory with a clock rate of 2.6 GHz and a bandwidth of 208 GB/sec. In each card, there are 14 stream multiprocessors (SMs) containing 2496 processors in total. The peak single and double precision floating point performance is 3.52 and 1.17 TFLOPS, respectively.All algorithms introduced in the present paper were developed and tested based on PETSc-3.10.2 that was compiled with “- -with-debugging = 0” option by GCC-4.4.7 and CUDA-aware OpenMPI-3.1.4 [37]. As the optimized PETSc is generally 2 times faster than the debug version of PETSc, we use the optimized PETSc on CPUs as reference when we analyze speedups between CPUs and GPUs. The developed CUDA kernels are compiled by nvcc from CUDA 6.5.14. And, all the floating point operations were performed in double precision.


### 4.2. Test Matrices

The test matrices were selected from the SuiteSparse Matrix Collection [38], which are shown in Figure 7. Each pointwise sparse matrix was transformed into the block sparse matrix with the BCSR format using cusparseDcsr2bsr routine from the NVIDIA's cuSPARSE library [28] with a given block size. Table 1 shows the descriptions before and after the transformation. The strategy we use to partition a matrix is(12)nrowsi=rl+1modr,l≠0,0≤i≤modr,l,rl,modr,l≤i≤l,where r is the total number of block rows and columns of the matrix, l is the number of MPI processes (GPUs), and rowsi is the number of block rows of the partitioned matrix assigned to the ith process (GPU). All processes call MatCreateMPIAIJWithArrays simultaneously to construct a parallel object of the PETSc matrix using their local parts.Figure 7
The patterns of the selected matrices.Table 1
Descriptions of the test matrices.
MatrixnnnzrnnzbsG3_circuit1585478766082679273952300492thermal21228045858031324560937915455atmosmodd1270432881488031760821908444atmosmodl14897521031976049658434298883venkat25624241717763156061073624


### 4.3. Results

We use 4 matrices fromG3_circuit to atmosmodl in Table 1 to investigate the performance of GPU-Direct communication in BSpMV and batch computing on vector operations in GMRES restarted with 30.Figure8 shows a comparison between the synchronized and asynchronized inner products with two batch sizes. All cases are conducted using 8 MPI processes and GPUs. Each MPI process is associated with a CPU core and controls a GPU card. Each process accumulates the wall clock time of computing the local inner products (line 9 in Algorithm 1) on the GPU until the case has converged. And, the wall clock times are averaged over all processes and shown in the y-axis of Figure 8. In both the synchronized and asynchronized cases, results using a batch size of 4 outperform that using a batch size of 2. This is to be expected because the larger batch size can reduce the global memory accesses in the loop. We notice that the asynchronized inner products with batch size 2 perform better than the synchronized ones using batch size 4. This demonstrates that overlapping the frequent memory copies and the computing on the host is effective for a smaller batch size. And, around 2.0x speedup can be obtained using asynchronized inner products with batch size 4 compared to using the synchronized method with batch size 2.Figure 8
Comparison between the synchronized and asynchronized inner products with different batch sizes.For updatingw⟶ at line 10 of Algorithm 1, Figure 9 shows the wall clock comparison between the implementation using cublasDaxpy from the cuSPARSE library and that using VecAXPYKerneljj≤8. Since the operations are memory bound, reducing the number of global access to w⟶ improves the performance. We see from Figure 9 that over 2.5x speedup can be obtained for all four cases.Figure 9
Comparison between cublasDaxpy andVecAXPYKerneljj≤8.Figure10 shows the performance for the BSpMV on the multi-GPUs. For each matrix, two strategies are studied and compared. One is using the BSpMV on the MDP to overlap with the GPU-Direct communication that is prepared for the BSpMV on the ODP. The other is to perform the BSpMV on the MDP, the communication, and the BSpMV on the ODP in a sequential order. The wall clock time for each configuration is estimated as the maximum value over all GPUs within 30 GMRES iterations (line 7 during the nested loop (lines 5–18) with m = 30) because the performance is determined by the slowest one. We observe that the overlapping strategy has shown a significant improvement in wall clock times for G3_circuit and thermal2. This is because these two matrices have much larger bandwidth that leads to more overhead in MPI communication than the other two matrices. Especially, when the number of GPUs increases to 8, the MPI overhead accounts for the major proportion of the overall time for these two cases and the decreasing computational workload is unable to overlap the overhead effectively. This is why we observe that 8 GPUs costs more time than 6 GPUs in the two cases. When less than 6 GPUs are used for atmosmodd and atmosmodl, the communication accounts for around 20% of the total GPU time due to the small bandwidth; therefore, we see relatively small gaps between the overlap bars and the nonoverlap bars. However, the proportions of the communication increase to nearly 30% and 40% when 6 GPUs and 8 GPUs are used, which makes the overlapping strategy has an obvious advantage over the nonoverlapping method.Figure 10
Performance of BSpMV with and without overlapping strategies.By employing the optimization techniques discussed above, we report the strong scalability for all the four cases in Figure11. The wall clock times are estimated from lines 3 to 21 (without the preconditioning step) in Algorithm 1 and averaged over three runs. For G3_circuit and thermal2 cases, we observe that only 33.3% and 35.8% parallel efficiencies can be obtained on 16 GPUs compared to 2 GPUs, respectively. This can be explained by the fact that the two matrices have large bandwidths and irregular sparsity patterns which results in unbalanced workload of matrix-vector multiplication and unbalanced MPI communication. Although the overlapping strategy is employed, the decreasing computational workload fails to hide the increasing communication overhead. The parallel efficiency for both atmosmodd and atmosmodl cases is over 50% on 16GPUs, which benefits from the smaller bandwidth and structured sparsity pattern.Figure 11
Strong scaling for all the four cases.Then, we investigate the performance of preconditioned GMRES with different preconditioning algorithms. Since the existing works of preconditioned GMRES focused either on scalar matrices or on single-GPU implementation, we compare our algorithm to the implementations based on PETSc and cuSPARSE libraries where block matrices are efficiently structured.Thermal2, mosmodd, atmosmodl, and venkat25 are selected for the experiments, and the results for the first three matrices are shown from Tables 2–4 and the result for venkat25 is shown in Figure 12. In each table, the first two row sections show the performance results for Algorithms 4 and 5. And, the remaining sections are the results for Algorithm 7 with different choices of the maximum number of iterations (itmax) used in Algorithm 7. In each section, the first three rows list the minimum, maximum, and average wall clock times spent on performing the preconditioning step in a restarted loop (lines 5–18 in Algorithm 1). The fourth row shows the total wall clock times spent on executing GMRES until the given relative tolerance is satisfied. And, the fifth row shows the total number of iterations required for the convergence tolerance.Table 2
thermal2: comparisons among three implementations on the preconditioning step of GMRES(30) with the relative tolerance, ϵ=10−3.
Algorithms2GPUs4GPUs6GPUs8GPUsPETScPrecondTprecond_minms1085523327196Tprecond_maxms1125569408295Tprecon_avgms1105550358260Tgmresms2760203115241168GMRES iterations58838686cuSPARSEPrecondTprecond_minms28391081689490Tprecond_maxms3429196614311254Tprecon_avgms313415691053789Tgmresms7378603645644012GMRES iterations58838686IterativePrecond (itmax=3)Tprecond_minms35617211072Tprecond_maxms36518512893Tprecon_avgms36117711887Tgmresms13851080800644GMRES iterations669699100IterativePrecond (itmax=4)Tprecond_minms47422914696Tprecond_maxms485246170124Tprecon_avgms480236157116Tgmresms15051192883705GMRES iterations60889294IterativePrecond (itmax=5)Tprecond_minms572286188120Tprecond_maxms586308212155Tprecon_avgms579294196145Tgmresms17271329996760GMRES iterations59868989Table 3
atmosmodd: comparisons among three implementations on the preconditioning step of GMRES(30) with the relative tolerance, ϵ=10−5.
Algorithms2GPUs4GPUs6GPUs8GPUsPETScPrecondTprecond_minms526262180127Tprecond_maxms580293203142Tprecon_avgms553282191135Tgmresms331114531085765GMRES iterations136116124120cuSPARSEPrecondTprecond_minms285.3173148.3135.7Tprecond_maxms285.5173.4149137.3Tprecon_avgms285.4173.2148.5136.7Tgmresms1886969854741GMRES iterations136116124120IterativePrecond (itmax=3)Tprecond_minms215.2108.372.455Tprecond_maxms216.7108.87355.4Tprecon_avgms216108.672.755.2Tgmresms1814824661469GMRES iterations160137157143IterativePrecond (itmax  = 4)Tprecond_minms286143.914696Tprecond_maxms286.4144.6170124Tprecon_avgms286.2144.2157116Tgmresms1900934.3883705GMRES iterations1381289294IterativePrecond (itmax = 5)Tprecond_minms356.9179.5120.591Tprecond_maxms357.4180.2121.291.4Tprecon_avgms357.2179.9120.891.2Tgmresms20431056897646GMRES iterations125124154148Table 4
atmosmodl: comparisons among three implementations on the preconditioning step of GMRES (30) with the relative tolerance, ϵ=10−5.
Algorithms2GPUs4GPUs6GPUs8GPUsPETScPrecondTprecond_minms539276183132Tprecond_maxms583293199146Tprecon_avgms561283193141Tgmresms1628944690559GMRES iterations63707479cuSPARSEPrecondTprecond_minms408.8230.8191.2175.0Tprecond_maxms409231.4192.5176.2Tprecon_avgms408.9231.1191.8175.7Tgmresms1236768656621GMRES iterations63707479IterativePrecond (itmax=2)Tprecond_minms182.292.162.147.1Tprecond_maxms182.992.662.547.6Tprecon_avgms182.592.362.247.4Tgmresms887497382327GMRES iterations77828794IterativePrecond (itmax=3)Tprecond_minms272.1137.492.470.8Tprecond_maxms272.5138.392.870.4Tprecon_avgms272.3137.792.770.3Tgmresms1001580422355GMRES iterations68757882IterativePrecond (itmax=4)Tprecond_minms361.5182.8123.193.3Tprecond_maxms362.6183.4123.493.6Tprecon_avgms362.1183.1123.393.4Tgmresms1172667494402GMRES iterations65727680Figure 12
venkat25: comparison ofTgmres with increasing itmax (relative tolerance: 10−5).In Table2, the thermal2 case shows that the preconditioning step using cuSPARSEPrecond is much slower than PETScPrecond. This demonstrates that the matrix pattern of thermal2 is not suitable for level-scheduling-based algorithms. To estimate the degree of load balance for the preconditioning step, we introduce the load balance factor as(13)LBF=Tprecond_maxTprecond_avg.We observe that the LBFs for cuSPARSEPrecond increases from 1.09 to 1.58 as the number of GPUs increases, which is larger than that for PETScPrecond. This can be explained by the fact that the level-scheduling algorithm extracts much different levels of parallelism from the partitioned matricesABi with different matrix patterns. The LBFs for IterativePrecond do not exceed 1.1, which is better than that for both PETscPrecond and cuSPARSEPrecond. This is due to the fully parallelized feature of IterativePrecond. In all three implementations, the total number of GMRES iterations used to converge to the given tolerance increases with increasing the number of GPUs. This is to be expected because the number of iterations usually increases when the block Jacobi method is applied to increasing partitions. Compared to the exact triangular solver using PETScPrecond and cuSPARSEPrecond, the iterative, inexact solver using IterativePrecond makes GMRES(30) require more iterations to converge to the same tolerance. And, the number of iterations for GMRES(30) can be reduced by increasing the number of times the iterative process is executed. Despite more GMRES (30) iterations are required, IterativePrecond shows an advantage over the other two methods in terms of wall clock times. When IterativePrecond with itmax=3 is used for comparisons, we obtain 1.81x and 6.23x on 8 GPUs over the implementations of GMRES using PETScPrecond and cuSPARSEPrecond, respectively.We also compare our CPU + GPU-based GMRES (30) with IterativePrecond to the corresponding CPU only implementation in PETSc, which is shown in Figure13. For a number of N in the x-axis of Figure 13, the baseline bar in blue is obtained by running the PETSc’s CPU-only (pure MPI) implementation on N CPU cores each of which is associated with a MPI process, while the other bar is obtained by running the CPU + GPU implementation on N GPUs each of which is controlled by a MPI process. By using 8 GPUs, we obtain over 4.4x speedups compared to 8 CPU cores, which greatly reduces the computing time of GMRES(30) algorithm for the problem with a fixed size.Figure 13
thermal2: comparisons of the wall clock times of GMRES(30) between the CPU-only implementation using PETSc and out CPU + GPU implementation.For theatmosmodd case, the right-hand vector comes from the accompanying data of the case. In Table 3, the atmosmodd case shows that both cuSPARSEPrecond and IterativePrecond are faster than PETScPrecond. The pointwise matrix pattern of atmosmodd has been reported to be suitable for level-scheduling algorithms [2], and we notice in this experiment that the benefit of the pattern is still kept when atmosmodd is transformed into the BCSR format. We observe that Tgmres using cuSPARSEPrecond is about 1.76x faster than that using PETScPrecond when 2 GPUs are used. However, the speedup reduces as the number of GPUs increases, and Tgmres with cuSPARSEPrecond does not show persuasive advantage over Tgmres with PETScPrecond when 8 GPUs are used. By contrast, although IterativePrecond results in more GMRES iterations, Tgmres using IterativePrecond with itmax=3 maintains more than 1.6x speedups compared to Tgmres using PETScPrecond when the number of GPUs increases from 2 to 8. In terms of LBFs for the preconditioning step, both cuSPARSEPrecond and IterativePrecond performs better than PETScPrecond because of the parallelism and the relatively balanced partitions. As we expected, the number of GMRES iterations is gradually reduced with the increasing value of itmax, but leading to more overhead in every 30 preconditioning steps, and thus makes the total wall clock time of GMRES increase. Due to the stochastic feature of updating solutions in IterativePrecond, the GMRES iterations could be less than that using the exact preconditioners. In this case, for example, when IterativePrecond with itmax=5 is employed, only 125 iterations are required to make GMRES(30) converge, which is 11 less than that using PETScPrecond. Figure 14 shows the speedups obtained by the CPU + GPU implementation over the PETSc implementation on CPUs. And, above 3.0 x speedup can be obtained for all configurations.Figure 14
atmosmodd: comparisons of the wall clock times of GMRES (30) between the CPU-only implementation using PETSc and out CPU + GPU implementation.The right-hand vector for theatmosmodl case also comes with the matrix. The results are shown in Table 4. We find that IterativePrecond using itmax=4 is sufficient for GMRES (30) to converge with nearly the same number of iterations as that obtained by the exact preconditioners. The scalability of IterativePrecond is much better than cuSPARSEPrecond because there is no data dependencies when IterativePrecond is applied, and cuSPARSEPrecond is slower than PETScPrecond when 8 GPUs are used. Tgmres using IterativePrecond with itmax=2 shows a speedup of 1.71x and 1.90x over that using PETScPrecond and cuSPARSEPrecond, respectively. Compared to the CPU-only implementation using 8 MPI processes, 2.97x speedup can be obtained using 8 GPUs, which is shown in Figure 15.Figure 15
atmosmodl: comparisons of the wall clock times of GMRES(30) between the CPU-only implementation using PETSc and out CPU + GPU implementation.Since the scale ofvenkat25 is not large enough to gain speedups on GPUs over the pure MPI implementation, we use this case only to demonstrate the advantage of the inexact triangular solve over the cuSPARSE-based triangular solve. Figure 12 shows a comparison of Tgmres with increasing itmax. We observe a optimal value of 6 for itmax because Tgmres attains the minimum value at itmax=6. Although IterativePrecond results in 69 more GMRES iterations to converge, it still outperforms cuSPARSEPrecond. When 4 GPUs are used, GMRES with IterativePrecond performs 4.37x faster than that with cuSPARSEPrecond.


## 4.1. Test Platform

In our experiments, the tests were conducted on a heterogeneous cluster at Computer Network Information Center, Chinese Academy of Sciences. The cluster consists of 270 blades nodes, 30 GPGPU nodes, and 40 MIC nodes. Each GPGPU node is configured with 2 Intel E5-2680 V2 (Ivy Bridge | 10C | 2.8 GHz) CPUs, 128 GB DDR3 ECC 1866 MHz memory, and two NVIDIA Tesla K20 GPGPU cards. Each card is configured by 320 bit GDDR5 5 GB memory with a clock rate of 2.6 GHz and a bandwidth of 208 GB/sec. In each card, there are 14 stream multiprocessors (SMs) containing 2496 processors in total. The peak single and double precision floating point performance is 3.52 and 1.17 TFLOPS, respectively.All algorithms introduced in the present paper were developed and tested based on PETSc-3.10.2 that was compiled with “- -with-debugging = 0” option by GCC-4.4.7 and CUDA-aware OpenMPI-3.1.4 [37]. As the optimized PETSc is generally 2 times faster than the debug version of PETSc, we use the optimized PETSc on CPUs as reference when we analyze speedups between CPUs and GPUs. The developed CUDA kernels are compiled by nvcc from CUDA 6.5.14. And, all the floating point operations were performed in double precision.


## 4.2. Test Matrices

The test matrices were selected from the SuiteSparse Matrix Collection [38], which are shown in Figure 7. Each pointwise sparse matrix was transformed into the block sparse matrix with the BCSR format using cusparseDcsr2bsr routine from the NVIDIA's cuSPARSE library [28] with a given block size. Table 1 shows the descriptions before and after the transformation. The strategy we use to partition a matrix is(12)nrowsi=rl+1modr,l≠0,0≤i≤modr,l,rl,modr,l≤i≤l,where r is the total number of block rows and columns of the matrix, l is the number of MPI processes (GPUs), and rowsi is the number of block rows of the partitioned matrix assigned to the ith process (GPU). All processes call MatCreateMPIAIJWithArrays simultaneously to construct a parallel object of the PETSc matrix using their local parts.Figure 7
The patterns of the selected matrices.Table 1
Descriptions of the test matrices.
MatrixnnnzrnnzbsG3_circuit1585478766082679273952300492thermal21228045858031324560937915455atmosmodd1270432881488031760821908444atmosmodl14897521031976049658434298883venkat25624241717763156061073624


## 4.3. Results

We use 4 matrices fromG3_circuit to atmosmodl in Table 1 to investigate the performance of GPU-Direct communication in BSpMV and batch computing on vector operations in GMRES restarted with 30.Figure8 shows a comparison between the synchronized and asynchronized inner products with two batch sizes. All cases are conducted using 8 MPI processes and GPUs. Each MPI process is associated with a CPU core and controls a GPU card. Each process accumulates the wall clock time of computing the local inner products (line 9 in Algorithm 1) on the GPU until the case has converged. And, the wall clock times are averaged over all processes and shown in the y-axis of Figure 8. In both the synchronized and asynchronized cases, results using a batch size of 4 outperform that using a batch size of 2. This is to be expected because the larger batch size can reduce the global memory accesses in the loop. We notice that the asynchronized inner products with batch size 2 perform better than the synchronized ones using batch size 4. This demonstrates that overlapping the frequent memory copies and the computing on the host is effective for a smaller batch size. And, around 2.0x speedup can be obtained using asynchronized inner products with batch size 4 compared to using the synchronized method with batch size 2.Figure 8
Comparison between the synchronized and asynchronized inner products with different batch sizes.For updatingw⟶ at line 10 of Algorithm 1, Figure 9 shows the wall clock comparison between the implementation using cublasDaxpy from the cuSPARSE library and that using VecAXPYKerneljj≤8. Since the operations are memory bound, reducing the number of global access to w⟶ improves the performance. We see from Figure 9 that over 2.5x speedup can be obtained for all four cases.Figure 9
Comparison between cublasDaxpy andVecAXPYKerneljj≤8.Figure10 shows the performance for the BSpMV on the multi-GPUs. For each matrix, two strategies are studied and compared. One is using the BSpMV on the MDP to overlap with the GPU-Direct communication that is prepared for the BSpMV on the ODP. The other is to perform the BSpMV on the MDP, the communication, and the BSpMV on the ODP in a sequential order. The wall clock time for each configuration is estimated as the maximum value over all GPUs within 30 GMRES iterations (line 7 during the nested loop (lines 5–18) with m = 30) because the performance is determined by the slowest one. We observe that the overlapping strategy has shown a significant improvement in wall clock times for G3_circuit and thermal2. This is because these two matrices have much larger bandwidth that leads to more overhead in MPI communication than the other two matrices. Especially, when the number of GPUs increases to 8, the MPI overhead accounts for the major proportion of the overall time for these two cases and the decreasing computational workload is unable to overlap the overhead effectively. This is why we observe that 8 GPUs costs more time than 6 GPUs in the two cases. When less than 6 GPUs are used for atmosmodd and atmosmodl, the communication accounts for around 20% of the total GPU time due to the small bandwidth; therefore, we see relatively small gaps between the overlap bars and the nonoverlap bars. However, the proportions of the communication increase to nearly 30% and 40% when 6 GPUs and 8 GPUs are used, which makes the overlapping strategy has an obvious advantage over the nonoverlapping method.Figure 10
Performance of BSpMV with and without overlapping strategies.By employing the optimization techniques discussed above, we report the strong scalability for all the four cases in Figure11. The wall clock times are estimated from lines 3 to 21 (without the preconditioning step) in Algorithm 1 and averaged over three runs. For G3_circuit and thermal2 cases, we observe that only 33.3% and 35.8% parallel efficiencies can be obtained on 16 GPUs compared to 2 GPUs, respectively. This can be explained by the fact that the two matrices have large bandwidths and irregular sparsity patterns which results in unbalanced workload of matrix-vector multiplication and unbalanced MPI communication. Although the overlapping strategy is employed, the decreasing computational workload fails to hide the increasing communication overhead. The parallel efficiency for both atmosmodd and atmosmodl cases is over 50% on 16GPUs, which benefits from the smaller bandwidth and structured sparsity pattern.Figure 11
Strong scaling for all the four cases.Then, we investigate the performance of preconditioned GMRES with different preconditioning algorithms. Since the existing works of preconditioned GMRES focused either on scalar matrices or on single-GPU implementation, we compare our algorithm to the implementations based on PETSc and cuSPARSE libraries where block matrices are efficiently structured.Thermal2, mosmodd, atmosmodl, and venkat25 are selected for the experiments, and the results for the first three matrices are shown from Tables 2–4 and the result for venkat25 is shown in Figure 12. In each table, the first two row sections show the performance results for Algorithms 4 and 5. And, the remaining sections are the results for Algorithm 7 with different choices of the maximum number of iterations (itmax) used in Algorithm 7. In each section, the first three rows list the minimum, maximum, and average wall clock times spent on performing the preconditioning step in a restarted loop (lines 5–18 in Algorithm 1). The fourth row shows the total wall clock times spent on executing GMRES until the given relative tolerance is satisfied. And, the fifth row shows the total number of iterations required for the convergence tolerance.Table 2
thermal2: comparisons among three implementations on the preconditioning step of GMRES(30) with the relative tolerance, ϵ=10−3.
Algorithms2GPUs4GPUs6GPUs8GPUsPETScPrecondTprecond_minms1085523327196Tprecond_maxms1125569408295Tprecon_avgms1105550358260Tgmresms2760203115241168GMRES iterations58838686cuSPARSEPrecondTprecond_minms28391081689490Tprecond_maxms3429196614311254Tprecon_avgms313415691053789Tgmresms7378603645644012GMRES iterations58838686IterativePrecond (itmax=3)Tprecond_minms35617211072Tprecond_maxms36518512893Tprecon_avgms36117711887Tgmresms13851080800644GMRES iterations669699100IterativePrecond (itmax=4)Tprecond_minms47422914696Tprecond_maxms485246170124Tprecon_avgms480236157116Tgmresms15051192883705GMRES iterations60889294IterativePrecond (itmax=5)Tprecond_minms572286188120Tprecond_maxms586308212155Tprecon_avgms579294196145Tgmresms17271329996760GMRES iterations59868989Table 3
atmosmodd: comparisons among three implementations on the preconditioning step of GMRES(30) with the relative tolerance, ϵ=10−5.
Algorithms2GPUs4GPUs6GPUs8GPUsPETScPrecondTprecond_minms526262180127Tprecond_maxms580293203142Tprecon_avgms553282191135Tgmresms331114531085765GMRES iterations136116124120cuSPARSEPrecondTprecond_minms285.3173148.3135.7Tprecond_maxms285.5173.4149137.3Tprecon_avgms285.4173.2148.5136.7Tgmresms1886969854741GMRES iterations136116124120IterativePrecond (itmax=3)Tprecond_minms215.2108.372.455Tprecond_maxms216.7108.87355.4Tprecon_avgms216108.672.755.2Tgmresms1814824661469GMRES iterations160137157143IterativePrecond (itmax  = 4)Tprecond_minms286143.914696Tprecond_maxms286.4144.6170124Tprecon_avgms286.2144.2157116Tgmresms1900934.3883705GMRES iterations1381289294IterativePrecond (itmax = 5)Tprecond_minms356.9179.5120.591Tprecond_maxms357.4180.2121.291.4Tprecon_avgms357.2179.9120.891.2Tgmresms20431056897646GMRES iterations125124154148Table 4
atmosmodl: comparisons among three implementations on the preconditioning step of GMRES (30) with the relative tolerance, ϵ=10−5.
Algorithms2GPUs4GPUs6GPUs8GPUsPETScPrecondTprecond_minms539276183132Tprecond_maxms583293199146Tprecon_avgms561283193141Tgmresms1628944690559GMRES iterations63707479cuSPARSEPrecondTprecond_minms408.8230.8191.2175.0Tprecond_maxms409231.4192.5176.2Tprecon_avgms408.9231.1191.8175.7Tgmresms1236768656621GMRES iterations63707479IterativePrecond (itmax=2)Tprecond_minms182.292.162.147.1Tprecond_maxms182.992.662.547.6Tprecon_avgms182.592.362.247.4Tgmresms887497382327GMRES iterations77828794IterativePrecond (itmax=3)Tprecond_minms272.1137.492.470.8Tprecond_maxms272.5138.392.870.4Tprecon_avgms272.3137.792.770.3Tgmresms1001580422355GMRES iterations68757882IterativePrecond (itmax=4)Tprecond_minms361.5182.8123.193.3Tprecond_maxms362.6183.4123.493.6Tprecon_avgms362.1183.1123.393.4Tgmresms1172667494402GMRES iterations65727680Figure 12
venkat25: comparison ofTgmres with increasing itmax (relative tolerance: 10−5).In Table2, the thermal2 case shows that the preconditioning step using cuSPARSEPrecond is much slower than PETScPrecond. This demonstrates that the matrix pattern of thermal2 is not suitable for level-scheduling-based algorithms. To estimate the degree of load balance for the preconditioning step, we introduce the load balance factor as(13)LBF=Tprecond_maxTprecond_avg.We observe that the LBFs for cuSPARSEPrecond increases from 1.09 to 1.58 as the number of GPUs increases, which is larger than that for PETScPrecond. This can be explained by the fact that the level-scheduling algorithm extracts much different levels of parallelism from the partitioned matricesABi with different matrix patterns. The LBFs for IterativePrecond do not exceed 1.1, which is better than that for both PETscPrecond and cuSPARSEPrecond. This is due to the fully parallelized feature of IterativePrecond. In all three implementations, the total number of GMRES iterations used to converge to the given tolerance increases with increasing the number of GPUs. This is to be expected because the number of iterations usually increases when the block Jacobi method is applied to increasing partitions. Compared to the exact triangular solver using PETScPrecond and cuSPARSEPrecond, the iterative, inexact solver using IterativePrecond makes GMRES(30) require more iterations to converge to the same tolerance. And, the number of iterations for GMRES(30) can be reduced by increasing the number of times the iterative process is executed. Despite more GMRES (30) iterations are required, IterativePrecond shows an advantage over the other two methods in terms of wall clock times. When IterativePrecond with itmax=3 is used for comparisons, we obtain 1.81x and 6.23x on 8 GPUs over the implementations of GMRES using PETScPrecond and cuSPARSEPrecond, respectively.We also compare our CPU + GPU-based GMRES (30) with IterativePrecond to the corresponding CPU only implementation in PETSc, which is shown in Figure13. For a number of N in the x-axis of Figure 13, the baseline bar in blue is obtained by running the PETSc’s CPU-only (pure MPI) implementation on N CPU cores each of which is associated with a MPI process, while the other bar is obtained by running the CPU + GPU implementation on N GPUs each of which is controlled by a MPI process. By using 8 GPUs, we obtain over 4.4x speedups compared to 8 CPU cores, which greatly reduces the computing time of GMRES(30) algorithm for the problem with a fixed size.Figure 13
thermal2: comparisons of the wall clock times of GMRES(30) between the CPU-only implementation using PETSc and out CPU + GPU implementation.For theatmosmodd case, the right-hand vector comes from the accompanying data of the case. In Table 3, the atmosmodd case shows that both cuSPARSEPrecond and IterativePrecond are faster than PETScPrecond. The pointwise matrix pattern of atmosmodd has been reported to be suitable for level-scheduling algorithms [2], and we notice in this experiment that the benefit of the pattern is still kept when atmosmodd is transformed into the BCSR format. We observe that Tgmres using cuSPARSEPrecond is about 1.76x faster than that using PETScPrecond when 2 GPUs are used. However, the speedup reduces as the number of GPUs increases, and Tgmres with cuSPARSEPrecond does not show persuasive advantage over Tgmres with PETScPrecond when 8 GPUs are used. By contrast, although IterativePrecond results in more GMRES iterations, Tgmres using IterativePrecond with itmax=3 maintains more than 1.6x speedups compared to Tgmres using PETScPrecond when the number of GPUs increases from 2 to 8. In terms of LBFs for the preconditioning step, both cuSPARSEPrecond and IterativePrecond performs better than PETScPrecond because of the parallelism and the relatively balanced partitions. As we expected, the number of GMRES iterations is gradually reduced with the increasing value of itmax, but leading to more overhead in every 30 preconditioning steps, and thus makes the total wall clock time of GMRES increase. Due to the stochastic feature of updating solutions in IterativePrecond, the GMRES iterations could be less than that using the exact preconditioners. In this case, for example, when IterativePrecond with itmax=5 is employed, only 125 iterations are required to make GMRES(30) converge, which is 11 less than that using PETScPrecond. Figure 14 shows the speedups obtained by the CPU + GPU implementation over the PETSc implementation on CPUs. And, above 3.0 x speedup can be obtained for all configurations.Figure 14
atmosmodd: comparisons of the wall clock times of GMRES (30) between the CPU-only implementation using PETSc and out CPU + GPU implementation.The right-hand vector for theatmosmodl case also comes with the matrix. The results are shown in Table 4. We find that IterativePrecond using itmax=4 is sufficient for GMRES (30) to converge with nearly the same number of iterations as that obtained by the exact preconditioners. The scalability of IterativePrecond is much better than cuSPARSEPrecond because there is no data dependencies when IterativePrecond is applied, and cuSPARSEPrecond is slower than PETScPrecond when 8 GPUs are used. Tgmres using IterativePrecond with itmax=2 shows a speedup of 1.71x and 1.90x over that using PETScPrecond and cuSPARSEPrecond, respectively. Compared to the CPU-only implementation using 8 MPI processes, 2.97x speedup can be obtained using 8 GPUs, which is shown in Figure 15.Figure 15
atmosmodl: comparisons of the wall clock times of GMRES(30) between the CPU-only implementation using PETSc and out CPU + GPU implementation.Since the scale ofvenkat25 is not large enough to gain speedups on GPUs over the pure MPI implementation, we use this case only to demonstrate the advantage of the inexact triangular solve over the cuSPARSE-based triangular solve. Figure 12 shows a comparison of Tgmres with increasing itmax. We observe a optimal value of 6 for itmax because Tgmres attains the minimum value at itmax=6. Although IterativePrecond results in 69 more GMRES iterations to converge, it still outperforms cuSPARSEPrecond. When 4 GPUs are used, GMRES with IterativePrecond performs 4.37x faster than that with cuSPARSEPrecond.


## 5. Conclusion and Future Work

In this work, focusing on the block sparse matrices, we present the development of a PETSc-enabled GMRES algorithm with various optimizations and preconditioning implementations on a GPU cluster. Motivated by the fined-grained feature of the inexact preconditioners in recent years, we propose a GPU algorithm to conduct a blockwise triangular solve inexactly and iteratively. On the block matrices selected from the SuiteSparse matrix collection, our experiments have shown that the preconditioned GMRES with inexact triangular solves outperform that with exact ones implemented by the cuSPARSE library and achieve up to 4.4x speedup using 8 GPUs over the CPU-only implementation using 8 MPI processes. Efforts are undergoing to study the preconditioned GMRES algorithm using the Additive Schwarz Method with the proposed inexact triangular solves on multi-GPUs.



---

*Source: 6804723-2021-02-28.xml*