## 1 Introduction  
Preconditioning techniques are of great importance for solving linear systems because they always make the system easier to solve in terms of the number of iterations in a given iterative method. Incomplete LU (ILU) [1, 2] based preconditioner in which the lower and upper matrices are factorized to approximate the preconditioning matrix is known as an effective and stable preconditioner for many problems. In this method, the preconditioning step is accomplished by solving two sparse triangular systems.  
GPU technique has been very attractive for the recent decade since CUDA (Compute Unified Device Architecture) [3] was released by NVIDIA in 2007. Many linear solvers [4, 5, 6] were accelerated by the fine-grained programming model. However, the solve of sparse triangular systems [1] where traditional forward and backward substitutions are employed is regards as the performance impediment in a preconditioned linear solver on GPUs because the forward (backward) substitutions are designed from the completely sequential programming idea. As the parallel techniques are very popular in today's parallel computers consisting of CPUs and GPUs, some strategies for parallelizing the triangular solves has been proposed. The popular one is the level-scheduling algorithm [1, 2, 7] that extracts the parallelism from the factors as a group of levels and data in each level can be operated in parallel. NVIDIA's cuSPAMSE [8] implements a GPU version of this algorithm. However, the parallelism from this algorithm has a strong connection to the sparsity pattern of the original matrix, thus it fails on GPUs for acceleration sometimes. A synchronization-free algorithm for the parallelization of sparse triangular solves [9, 10] is proposed without estimating the possible levels by employing the compressed sparse column format for the storing of the given matrix.  
In recent years, the idea of asynchronous and inexact computing [11, 12, 13, 14, 15, 16, 17] for is attractive because it aims to balance exactness and performance in today's parallel architecture.  
For preconditioners, the representative works are the iterative ILU factorization [12; 13], asynchronous triangular solves [15; 16] and incomplete sparse approximate inverses (ISAI) for triangular solves [11; 17]. The advantage and common ground of these algorithms is that they eliminate the parallel bottleneck by removing the strong dependency that ILU factorizations and triangular solves hold, which gains enough parallelism on the fine-grained parallel accelerators such as GPUs. As a result, exactness becomes the main challenge in these algorithms because removing data dependency may make the ILU factorization or triangular solves inexact. However, as preconditioners, these methods are always allowed and just affect the convergence of the linear system in terms of number of iterations. Although the number of iterations increases due to the degradation of the quality for the preconditioner, the total wall clock time for the linear solver could be reduced because the preconditioning step in every iteration is greatly accelerated.  
Motivated by the GPU ISAI for scalar matrices [17], in our previous work [18], we proposed a warp-based GPU algorithm for block matrices. In this work, we investigate the algorithm on multiple GPUs. Our contributions are  
* We extend the single GPU implementation to multiple GPUs version using the Block-Jacobi preconditioner.
* A detailed multi-GPU implementation based on the open source scientific computing framework, PETSc, is introduced and investigated.
* We make detailed performance comparisons and analysis on the proposed algorithm using NVIDIA's Tesla K20 and Tesla V100 GPUs.  
For the convenience of our statement in this paper, we use many mathematical symbols and abbreviations that are shown in Table 1.