## 4 Experiments
### Results
We observe from Figure 5 that solving SBTS achieved over 9\(\times\) speedup on Tesla V100 GPU compared to Tesla K20 GPU, and it shows nearly linear scalability up to four GPUs. The scalability is expected, because each warp independently solved a small block triangular system without any data dependency on other warps, which made maximum use of the fine-grained feature of GPUs. On K20 GPUs, because the L1 cache was not used for global memory loads but only for register spills to local memory [28], we preferred a configuration of 48 KB shared memory and 16 KB L1 cache to guarantee that there were enough shared memory spaces for storing temporary variables whose sizes varied with block sizes. According to NVIDIA's profiler, the achieved_occupancy for warps was about 90% on Tesla V100 GPUs, whereas it was 73% on Tesla K20 GPUs. This is because warps had overlap with loading data from the global memory, and using L1 cache on Tesla V100 with low-latency was effective for latency hiding. \(T_{ESP}\) and \(T_{cs2 sc}\) implemented using cuSPARSE library were much reduced due to the performance improvement in both GPU hardware and the software library. Figure 5: atmosmodd case: overhead breakdown of gaining \(L_{i}^{inv}\) and \(U_{i}^{inv}\) using Block–ISAI with \(S(|X_{i}|)\), \(S(|X_{i}|^{2})\) and \(S(|X_{i}|^{3})\). (**a**) \(T_{ESP}\): time used for estimating the sparsity pattern. (**b**) \(T_{SBTS}\): time used for solving small block triangular systems. (**c**) \(T_{cs2 sc}\): time used for format transformations. Figure 3b shows the BjPB-ISAI implementation had an advantage over the cuSPARSE implementation on both K20 and V100 GPUs in terms of \(T_{gmres}\). By counting \(T_{BISAI}\) for BjPB-ISAI algorithm, we still observed GMRES(30) with BjPB-ISAI(\(|X_{i}|\)) could achieve 1.6x over that with cuSPARSE on Tesla K20 platform and GMRES(30) with BjPB-ISAI(\(|X_{i}|^{2}\)) using four V100 GPUs could achieve 1.87x over that with cuSPARSE on Tesla V100 platform. In the second case, we solved a non-linear system from the 2D nonlinear driven cavity problem expressed as  
\[-\triangle u-\frac{\partial\omega}{\partial y} = 0,\] \[-\triangle v-\frac{\partial\omega}{\partial x} = 0,\] \[-\triangle\omega+u\frac{\partial\omega}{\partial x}+v\frac{ \partial\omega}{\partial y} = 0.\]  
For the velocity variables \((u,v)\), the no-slip, rigid-wall Dirichlet conditions were applied. For vorticity \(\omega\), the boundary condition is given by \(\omega=-\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\). The Reynolds number was set to 1 in this case. The problem was solved on a 600 \(\times\) 600 mesh in a unit square domain by using a five-point stencil finite difference method. Since each mesh cell was associated with three variables, the resulting linear system at each Newton step was formed by the Jacobian matrix with a block size of three. The relative tolerance for each linear system was set to \(10^{-5}\). The convergence criteria was set to \(10^{-8}\) for the non-linear process. In each subdomain, the block ILU(0) factorization was employed. We note that the linear system could not converge without using the ILU preconditioner. The non-linear system required three Newton steps to converge to the given tolerance. At each Newton step, a linear system consisting of 1,080,000 block rows (columns) and 1,797,600 non-zero blocks needed to be solved. The sparsity pattern for the linear system stayed unchanged at different Newton steps, therefore the data dependency and computational workload for the preconditioning step was also unchanged even though the data in the linear system changed at every Newton step. Figure 6 shows a comparison of \(T_{p}\) between Tesla K20 and Tesla V100 GPUs. Figures 7 and 8 show the comparisons of \(T_{gmres}\) and \(I_{gmres}\) in each Newton step, respectively.