## 3 Block-Jacobi with Block-ISAI on Multi-GPUs
### Multi-GPU Implementation
```
1:if\(myid==root\)then
2:\(\text{cusparseXcsr2 bsrNnz(handle,dir,\)\(cstrinfo,\)\(\ldots,\)\(bs,\)\(brinfo,\)\(\ldots,\)\(out:\)\(nnzb)\);
3:allocate memory for\(browptr,\)\(bcolval,\)\(bblkval\) on the device memory
4:\(\text{cusparseDcsr2 bsr(handle,dir,\)\(cstrinfo,\)\(\ldots,\)\(bs,\)\(bsrinfo,\)\(\ldots)\);
5:\(\text{cudaMemcpy}\)\(browptr,\)\(bcolval,\)\(bblkval\) back into the host memory
6:partition the block matrices by block rows (Figure 1a)
7:for\(iproc=1;iproc<p;iproc++\)do
8:// send the local parts\(browptr_{iproc},\)\(bcolval_{iproc},\) and \(bblkval_{iproc}\) to\(iproc\)
9:MPI_Send(...,to_iproc,\(\ldots\));
10:endfor
11:else
12:// receive from \(root\) as \(browptr_{iproc},\)\(bcolval_{iproc},\) and \(bblkval_{iproc}\)
13:MPI_Recv(...,\(from\_root,\)\(\ldots\));
14:endif
15:MatCreateMPIBAIJWithArrays(...,\(bs,\)\(...,\)\(browptr_{iproc},\)\(bcolval_{iproc},\)\(bblkval_{iproc},\)\(\&A\));
16:VecCreate(\(comm\_world,\&x\));
17:VecCreate(\(comm\_world,\&b\));
```  
**List 1** Procedure for transforming user provided matrices to PETSc's BCSR matrices. ```
1:if\(myid==root\)then
2:\(\text{cusparseXcsr2 bsrNnz(handle,dir,\)\(cstrinfo,\)\(\ldots,\)\(bs,\)\(brinfo,\)\(\ldots,\)\(out:\)\(nnzb)\);
3:allocate memory for\(browptr,\)\(bcolval,\)\(bblkval\) on the device memory
4:\(\text{cusparseDcsr2 bsr(handle,dir,\)\(cstrinfo,\)\(\ldots,\)\(bs,\)\(bsrinfo,\)\(\ldots)\);
5:\(\text{cudaMemcpy}\)\(browptr,\)\(bcolval,\)\(bblkval\) back into the host memory
6:partition the block matrices by block rows (Figure 1a)
7:for\(iproc=1;iproc<p;iproc++\)do
8:// send the local parts\(browptr_{iproc},\)\(bcolval_{iproc},\) and \(bblkval_{iproc}\) to\(iproc\)
9:MPI_Send(...,to_iproc,\(\ldots\));
10:endfor
11:else
12:// receive from \(root\) as \(browptr_{iproc},\)\(bcolval_{iproc},\) and \(bblkval_{iproc}\)
13:MPI_Recv(...,\(from\_root,\)\(\ldots\));
14:endif
15:MatCreateMPIBAIJWithArrays(...,\(bs,\)\(...,\)\(browptr_{iproc},\)\(bcolval_{iproc},\)\(bblkval_{iproc},\)\(\&A\));
16:VecCreate(\(comm\_world,\&x\));
17:VecCreate(\(comm\_world,\&b\));
```  
**List 2** Locally accessing block lower and upper factors. The key part of computation on GPU is the numerous small block triangular systems (SBTS) (line 9 and line 14) with a degree of parallelism of \(m\), where \(m\) is the number of block rows (columns) for \(L_{i}\) (\(U_{i}\)). Different from the large triangular systems in (8), one of the SBTS is viewed as a sub-task for a thread or several threads of a GPU rather than the whole task for a single GPU. How to map the sub-task to threads depends on the algorithmic technique. The SBTS are independent with each other that provides the possibility of exploiting the fine-grained parallelism on the GPU. The proposed warp-based algorithm on a single GPU for solving the SBTS [18] is applied to a sub-domain in the multi-GPU case. Figure 2 shows an example of solving the first four block columns of \(L_{i}^{inv}\) using a thread block containing 128 threads. From the hardware point of view, the basic unit of execution and dispatching on a streaming multiprocessor (SM) of the GPU are warps and a warp groups 32 consecutive threads that are executed in Single Instruction Multiple Thread (SIMT) model. Therefore, the threads in the thread block are divided into 4 warps each of which is implicitly issued by the warp scheduler on the GPU. A total of 32 threads in a warp solve one of the four SBTS with executing exactly the same instruction but accessing different data. The right part of Figure 2 illustrates the process of solving a \(5\times 5\) lower triangular system with block size 3 in \(warp_{3}\). To explain how the \(5\times 5\) lower triangular system is solved by \(warp_{3}\), we simply denote the system as  
\[WSol=rhs \tag{9}\]  
where \(W=L_{i}(R_{L}^{3},R_{L}^{3})\), \(Sol=L_{i}^{inv}(R_{L},3)\), and \(rhs=I_{3}^{R_{L}^{3}}\). The solution (\(Sol\), marked in black triangles in blue squares) is solved by looping over five block columns. In the process of handing the first block column for example, three steps are required.