## 3 Block-Jacobi with Block-ISAI on Multi-GPUs
### Preconditioning on Multi-GPUs  
On a distributed system with GPU accelerators, the matrix is always partitioned into sub-matrices each of which is assigned to an accelerator that a process maps to, which is shown in Figure 1a. The popular partitioning strategy is to divide a square matrix in successive rows. A straightforward way to construct \(M\) in a distributed system is to conduct incomplete LU factorization (ILU) on the diagonal part of the partitioned matrix residing in each process, which can be expressed as  
\[D(A_{i})\approx L_{i}U_{i}=M_{i} \tag{4}\]  
where \(A_{i}\) is the partitioned matrix in \(i\)th process, \(D(A_{i})\) is the main diagonal part of \(A_{i}\), \(L_{i}\) and \(U_{i}\) is the block lower and upper triangular factors obtained by the ILU factorization, respectively. The process is illustrated in Figure 1b. \[D(A_{i})\approx L_{i}U_{i}=M_{i} \tag{5}\]  
where \(A_{i}\) is the partitioned matrix in \(i\)th process, \(D(A_{i})\) is the main diagonal part of \(A_{i}\), \(L_{i}\) and \(U_{i}\) is the block lower and upper triangular factors obtained by the ILU factorization, respectively. The process is illustrated in Figure 1b. Figure 1: Matrix partition and the process of constructing preconditioners. We apply the Block-Jacobi preconditioning method [1], which also can be regarded as the Additive Schwarz Method (ASM) with zero-overlap [25], to perform the preconditioning step as  
\[M^{-1}v=\left[ \begin{array}{cccc}(M_{0})^{-1}&&\\ &(M_{1})^{-1}&&\\ &&\ldots&\\ &&&(M_{p-1})^{-1}\end{array} \right]\left[ \begin{array}{cccc}v_{0}\\ v_{1}\\ \ldots\\ v_{p-1}\end{array} \right]=\left[ \begin{array}{cccc}z_{0}\\ z_{1}\\ \ldots\\ z_{p-1}\end{array} \right] \tag{6}\]  
where \(M_{i}\) is defined as (5), \(v_{i}\) the partitioned vector of a global vector \(v\), \(z_{i}\) is the partitioned vector of the preconditioned vector \(z\), and \(p\) is the number of processes or accelerators. The task of the \(i\)th \((0\leq i<p)\) process is to compute \(z_{i}\) by solving  
\[M_{i}^{-1}v_{i}=z_{i} \tag{7}\]  
independently without any communication. The aim of performing ILU factorization in (5) is to form (7) two triangular systems  
\[L_{i}y_{i}=v_{i},U_{i}z_{i}=y_{i} \tag{8}\]  
which can be solved by forward and backward substitutions. Equation (8) is hard to parallelize in the ILU preconditioned krylov solvers [24] on parallel architectures because forward and backward substitutions are of inherently sequential feature. Applying ISAI idea [11; 17] to (8), we can transform the forward (backward) substitutions into highly parallelized operations consisting of sparse matrix-vector multiplications by computing the inverse of \(L_{i}\) and \(U_{i}\) approximately. In this way, the quality of the preconditioning may degrade and then results in more iterations to converge, but the overall time for the solver could be reduced due to the performance improvement at each preconditioning step. In this work, to seek a balance between performance and exactness, we combine Block-Jacobi preconditioner with Block-ISAI (BJPB-ISAI) in a multi-GPU scenario. We list BJPB-ISAI method in Algorithm 1. As the GPU workload needs to be launched from hosts, the algorithm runs by a hybrid of Message Passing Interface (MPI) programming model [26] and Compute Unified Development Architecture (CUDA) [27]. MPI is used to launch \(p\) host processes each of which is corresponding to a GPU card. Based on the partition strategy illustrated in Figure 1a, the local block triangular factors \(L_{i}\) and \(U_{i}\) can be obtained from the block ILU factorizations on \(D(A_{i})\), and they together with a partitioned vector \(v_{i}\) are used as the input parameters. The algorithm outputs the partitioned, preconditioned vector \(z_{i}\) on each GPU. We divide the content of the algorithm into two parts.