## 3 Block-Jacobi with Block-ISAI on Multi-GPUs
### Preconditioning on Multi-GPUs
The first part provides an approximation to the inverses of \(L_{i}\) and \(U_{i}\) requires to be preformed only once during the iterative solver because we consider fixed preconditioners in this work and the ILU factors do not change once the factorization is completed. Specifically, the sparsity pattern estimation for the inverses of \(L_{i}\) and \(U_{i}\), denoted as \(L_{i}^{inv}\) and \(U_{i}^{inv}\), is firstly performed using the symbolic matrix-matrix multiplication on \(S(|L_{i}|)\) and \(S(|U_{i}|)\) for \(k\) times, where \(|L_{i}|\) and \(|U_{i}|\) represent the absolute values of \(L_{i}\) and \(U_{i}\). The sparsity patterns become thicker and with increasing \(k\), but resulting in more memory and computation requirement. The parallel idea of ISAI [17; 18] is to solve \(L_{i}^{inv}\) by columns, so we transform \(L_{i}^{inv}\) and \(U_{i}^{inv}\) into BCSC formats for the convenience of accessing row indices in each column of \(L_{i}^{inv}\) or \(U_{i}^{inv}\). More importantly, BCSC formats make the access of \(L_{i}^{inv}\) or \(U_{i}^{inv}\) in the global memory in a coalesced manner. Then two GPU kernels that solve \(L_{i}^{inv}\) and \(U_{i}^{inv}\) in block columns [18] in parallel are launched. The major feature of the GPU algorithm is making full use of the fine-grained programming model of GPU and solves each block column by a warp consisting of 32 consecutive threads (CUDA threads are issued by warps implicitly [28]). The implementation issues are discussed in the next section. Once \(L_{i}^{inv}\) and \(U_{i}^{inv}\) are ready, the second part of Algorithm 1 performs block-based sparse matrix-vector multiplication instead of exactly solving block triangular systems. ```
1:(1) MPI process identifier \(i\);
2:(2) the lower and upper block triangular factors: \(L_{i}\) and \(U_{i}\);
3:(3) the partitioned vector \(v_{i}\);
4:the preconditioned vector \(z_{i}\);
5:map the \(i\)th process to a GPU with identifier \(\mathit{idve}\) on a node;
6:cudaSetDevice(\(\mathit{idve}\));
7:ifirstcallthen
8: use \(S(|L_{i}|^{k})\) and \(S(|U_{i}|^{k})\) to give an estimation of sparsity pattern of the inverses of \(L_{i}\) and \(U_{i}\): \(S(L_{i}^{inv})=S(|L_{i}|^{k})\) and \(S(U_{i}^{inv})=S(|U_{i}|^{k})\);
9: transform \(L_{i}^{inv}\) and \(U_{i}^{inv}\) to BCSC format;
10: // unroll the \(\mathit{for}\) loop on GPU to form a warp-based GPU kernel [18]
11:for\(\mathit{icol}=0\) to \(m-1\)do
12:\(R_{L}=\) row indices for \(L_{i}^{inv}(:,\mathit{icol})\);
13: locally and temporarily form and solve: \(L_{i}(R_{L},R_{L})L_{i}^{inv}(R_{L},\mathit{icol})=I_{\mathit{icol}}^{R_{L}}\);
14:endfor
15: // unroll the \(\mathit{for}\) loop on GPU to form a warp-based GPU kernel [18]
16:for\(\mathit{icol}=0\) to \(m-1\)do
17:\(R_{U}=\) row indices for \(U_{i}^{inv}(:,\mathit{icol})\);
18: locally and temporarily form and solve: \(U_{i}(R_{U},R_{U})U_{i}^{inv}(R_{U},\mathit{icol})=I_{\mathit{icol}}^{R_{U}}\);
19:endfor
20: // transform \(L_{i}^{inv}\) and \(U_{i}^{inv}\) back into BCSR format;
21:endif
22:\(y_{i}=L_{i}^{inv}v_{i}\);
23:\(z_{i}=U_{i}^{inv}y_{i}\);
```  
**Algorithm 1** Block-Jacobi preconditioning with Block-ISAI on multi-GPUs: \(M^{-1}v=z\).