# S-Step BiCGStab Algorithms for Geoscience Dynamic Simulations

Ani Anciaux-Sedrakian, Laura Grigori, Sophie Moufawad, Soleiman Yousef

<figure><img src="image_1.png" /></figure>

## To cite this version:

Ani Anciaux-Sedrakian, Laura Grigori, Sophie Moufawad, Soleiman Yousef. S-Step BiCGStab Algorithms for Geoscience Dynamic Simulations. Oil & Gas Science and Technology - Revue d'IFP Energies nouvelles, 2016, 71 (6), 10.2516/ogst/2016021. [hal-01436729](https://hal-01436729)

HAL Id: hal-01436729

https://hal.sorbonne-universite.fr/hal-01436729

Submitted on 16 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# S-Step BiCGStab Algorithms for Geoscience Dynamic Simulations

Ani Anciaux-Sedrakian<fnref n="1" />, Laura Grigori<fnref n="2" />,<fnref n="3" />, Sophie Moufawad<fnref n="1" />* and Soleiman Yousef<fnref n="1" />

**Abstract** — In basin and reservoir simulations, the most expensive and time consuming phase is solving systems of linear equations using Krylov subspace methods such as BiCGStab. For this reason, we explore the possibility of using communication avoiding Krylov subspace methods (s-step BiCGStab), that speedup of the convergence time on modern-day architectures, by restructuring the algorithms to reduce communication. We introduce some variants of s-step BiCGStab with better numerical stability for the targeted systems.

**Résumé** — Méthodes s-step BiCGStab appliquées en Géosciences — Dans les simulateurs d'écoulement en milieu poreux, comme les simulateurs de réservoir et de bassin, la résolution de système linéaire constitue l'étape la plus consommatrice en temps de calcul et peut même représenter jusqu'à 80 % du temps de la simulation. Ceci montre que la performance de ces simulateurs dépend fortement de l'efficacité des solveurs linéaires. En même temps, les machines parallèles modernes disposent d'un grand nombre de processeurs et d'unités de calcul massivement parallèle. Dans cet article, nous proposons de nouveaux algorithmes BiCGStab, basés sur l'algorithme à moindre communication nommé s-step, permettant d'éviter un certain nombre de communication afin d'exploiter pleinement les architectures hautement parallèles.

## INTRODUCTION

Many scientific problems require the solution of systems of linear equations of the form $Ax = b$, where the input matrix $A$ is very large and sparse. These systems arise mainly from the discretization of Partial Differential Equations (PDE), and are usually solved using Krylov subspace methods, such as Generalized Minimal RESidual (GMRES) [1], Conjugate Gradient (CG) [2] and Bi-Conjugated Gradient Stabilized (BiCGStab) [3].

In the case of basin modeling or reservoir simulations with highly heterogeneous data and complex geometries,

complex non-linear systems of PDE are solved. These PDE are discretized with a cell-centered finite volume scheme in space, leading to a non-linear system which is solved with an iterative Newton solver. At each Newton step, the system is linearized. Then, the generated large, sparse and unstructured linear system is solved using preconditioned GMRES, BiCG-Stab, CG, Orthomin or other preconditioned iterative methods. Some of the most commonly used preconditioners are ILU(k), ILUT, AMG and CPR-AMG. This resolution phase constitutes the most expensive part of the simulation. Thus we focus on linear solvers, since their efficiency is a key point for the simulator's performance.

<footnote n="1">¹ IFP Energies nouvelles, 1-4 avenue de Bois-Préau, 92852 Rueil-Malmaison Cedex - France</footnote>

<footnote n="2">² INRIA Paris, Alpines, 2 Rue Simone IFF, 75012 Paris - France</footnote>

<footnote n="3">³ UPMC - Univ Paris 6, CNRS UMR 7598, Laboratoire Jacques-Louis Lions, 4 Place Jussieu, 75005 Paris - France</footnote>

<footnote>e-mail: ani.anciaux-sedrakian@ifpen.fr - laura.grigori@inria.fr - sm101@aub.edu.lb - soleiman.yousef@ifpen.fr</footnote>

<footnote>* Corresponding author</footnote>
Furthermore, modern parallel computing resources are based on complex hardware architecture. They are composed of several multi-core processors and massively parallel processing units such as many-cores or General-Purpose GPU (GPGPU) cards. Most of the current algorithms are not able to fully exploit the highly parallel architectures. In fact, a severe degradation of performance is detected when the number of processing units is increased. This is due to the difference between the required time to perform floating point operations (flops) by processing units and the time to communicate the obtained results, where flops have become much cheaper than data communication.

Thus, recent research has focused on reformulating dense and sparse linear algebra algorithms with the aim of reducing and avoiding communication. These methods are referred to as communication avoiding methods, whereby communication refers to data movement between different processing units in parallel, and different levels of memory hierarchy. In the case of Krylov subspace methods, the introduced communication avoiding Krylov subspace methods [4-8] are based on s-step methods [9-11]. The goal is to restructure the algorithms to perform s iterations at a time by using kernels that avoid or reduce communication, such as the matrix powers kernel [12], Tall and Skinny QR (TSQR) [13], and Block Gram Schmidt (BGS) methods.

Our aim is to reduce the overall cost of the linear solver resolution phase in geoscience simulations, specifically basin and reservoir simulations, using a parallel implementation of BiCGStab that avoids communication on multi-core hardware (CA-BiCGStab) [5] and has a similar convergence behavior as the classical BiCGStab method. Communication Avoiding BiCGStab (CA-BiCGStab), which was introduced in [5], is a reformulation of BiCGStab into s-step BiCGStab that avoids communication. Thus, in this paper we study the convergence behavior of a sequential version of the unpreconditioned s-step BiCGStab [5], on matrices obtained from reservoir simulations, with different s values. The obtained results show that, for most of the tested matrices, s-step BiCGStab requires more iterations to converge than BiCGStab. Thus, we design new variants of s-step BiCGStab, that have the same convergence rate as BiCGStab for s values between 2 and 6, and reduce communication similarly to s-step BiCGStab.

In Section 1, we introduce BiCGStab, its reformulation to s-step BiCGStab [5], and we discuss the performance of s-step BiCGStab in geoscience applications, specifically reservoir simulations. Then, in Section 2, we introduce the new s-step BiCGStab variants that we call orthonormalized s-step BiCGStab, split orthonormalized s-step BiCGStab, and modified split orthonormalized s-step BiCGStab. In Section 3, we present the convergence results of the newly introduced s-step BiCGStab variants and compare them to that of s-step BiCGStab. Finally, we conclude.

# 1 FROM BICGSTAB TO S-STEP BICGSTAB

In this section we briefly introduce BiCGStab (Sect. 1.1) and s-step BiCGStab (Sect. 1.2). We show the relation between both methods and their convergence in reservoir simulations (Sect. 1.3).

## 1.1 BiCGStab

The Bi-Conjugate Gradient Stabilized method (BiCGStab), introduced by van der Vorst in 1992 [3], is an iterative Krylov subspace method that solves the general systems $Ax = b$. It is a variant of the Bi-Conjugate Gradient (BiCG) method that aims at smoothing BiCG's erratic convergence. At each iteration $m \ge 0$, $r_{m+1} = P_{m+1}(A)r_0$ is replaced by $r_{m+1} = Q_{m+1}(A)P_{m+1}(A)r_0$ where $Q_{m+1}(z) \in P_{m+1}$ and $P_{m+1}(z) \in P_{m+1}$ are polynomials of degree $m+1$. $Q_{m+1}(z)$ is chosen to be

$$Q_{m+1}(z) = \prod_{j=1}^{m+1} (1 - \omega_{j-1}z) = (1 - \omega_m z)Q_m(z)$$

where $\omega_m$ minimizes the norm of $r_{m+1}$.

BiCGStab, being a variant of BiCG, has a similar form. But the recurrence relations of $x_{m+1}$, $r_{m+1}$, $p_{m+1}$, $\alpha_m$ and $\beta_m$ are different.

$$x_{m+1} = x_m + \alpha_m p_m + \omega_m [r_m - \alpha_m A p_m] \quad (1)$$

$$r_{m+1} = (I - \omega_m A)[r_m - \alpha_m A p_m] \quad (2)$$

$$p_{m+1} = r_{m+1} + \beta_m (I - \omega_m A) p_m \quad (3)$$

where $r_0 = b - Ax_0$, and $p_0 = r_0$.

The scalars $\alpha_m$ and $\beta_m$ are defined as follows:

$$\alpha_m = \frac{\langle \tilde{r}_0, r_m \rangle}{\langle \tilde{r}_0, A p_m \rangle}, \quad \beta_m = \frac{\alpha_m \langle \tilde{r}_0, r_{m+1} \rangle}{\omega_m \langle \tilde{r}_0, r_m \rangle} \quad (4)$$

where $\tilde{r}_0$ is chosen such that $\langle \tilde{r}_0, r_0 \rangle \neq 0$ as shown in Algorithm 1. In general $\tilde{r}_0$ is set equal to $r_0$. As for $\omega_m$, it is defined by minimizing the norm of the residual $r_{m+1}$, i.e. $||(I - \omega_m A)(r_m - \alpha_m A p_m)|| = \min_{\omega \in \mathbb{R}} ||(I - \omega A)(r_m - \alpha_m A p_m)||$, where

$$\omega_m = \frac{\langle Ar_m - \alpha_m A^2 p_m, r_m - \alpha_m A p_m \rangle}{\langle Ar_m - \alpha_m A^2 p_m, Ar_m - \alpha_m A^2 p_m \rangle} \quad (5)$$

In addition, we have that for $m > 0$

$$\left\{ \begin{array}{l} p_m, r_m \in K_{2m+1}(A, p_0) + K_{2m}(A, r_0) \\ x_m - x_0 \in K_{2m}(A, p_0) + K_{2m-1}(A, r_0) \end{array} \right. \quad (6)$$

and more generally for $m \ge 0$ and $j > 0$

$$\left\{ \begin{array}{l} p_{m+j}, r_{m+j} \in K_{2j+1}(A, p_m) + K_{2j}(A, r_m) \\ x_{m+j} - x_m \in K_{2j}(A, p_m) + K_{2j-1}(A, r_m) \end{array} \right. \quad (7)$$
**Algorithm 1: BiCGStab**

**Input:** A, b, x₀, mmax: the maximum allowed iterations

**Output:** x<sub>m</sub>: the mth approximate solution satisfying the stopping criteria

1: Let r₀ = b - Ax₀, p₀ = r₀, ρ₀ = ⟨r₀, r₀⟩, and m = 0

2: Choose r̃₀ such that δ₀ = ⟨r̃₀, r₀⟩ ≠ 0.

3: **While** (√ρ_m > ϵ||b||_2 and m < m_max) **Do**

4: α_m = δ_m / ⟨r̃₀, Ap_m⟩

5: s = r_m - α_m Ap_m

6: t = As

7: ω_m = ⟨t, s⟩ / ⟨t, t⟩

8: x<sub>m+1</sub> = x<sub>m</sub> + α<sub>m</sub>p<sub>m</sub> + ω<sub>m</sub>s

9: r<sub>m+1</sub> = (I - ω<sub>m</sub>A)s

10: δ<sub>m+1</sub> = ⟨r̃₀, r<sub>m+1</sub>⟩

11: β_m = (δ<sub>m+1</sub>/δ_m)(α_m/ω_m)

12: p<sub>m+1</sub> = r<sub>m+1</sub> + β_m(I - ω_m A)p_m

13: ρ<sub>m+1</sub> = ⟨r<sub>m+1</sub>, r<sub>m+1</sub>⟩, m = m + 1

14: end for

At each iteration of Algorithm 1, two sparse matrix-vector multiplications, six saxpy's, and five dot products are computed. Given that each processor has the scalars and its corresponding part of the vectors, then the saxpy's can be parallelized without communication. However, this is not the case for the sparse matrix-vector multiplications and the dot products, which require communication to obtain the desired results. Such operations cause a severe performance degradation, especially when using modern computing resources.

## 1.2 S-Step BiCGStab

To reduce the communication in parallel and sequential implementations of BiCGStab, Carson et al. [5] introduced the s-step version of BiCGStab. The reformulation is based on the computation of s BiCGStab iterations at once, and on the fact that for m ≥ 0 and 1 ≤ j ≤ s

$$ \begin{cases} p_{m+j}, r_{m+j} \in \mathcal{K}_{2s+1}(A, p_m) + \mathcal{K}_{2s}(A, r_m) \\ x_{m+j} - x_m \in \mathcal{K}_{2s}(A, p_m) + \mathcal{K}_{2s-1}(A, r_m) \end{cases} \quad (8) $$

since $\mathcal{K}_{2j+1}(A, z) \subseteq \mathcal{K}_{2s+1}(A, z)$ for any $z \neq 0$.

The goal is to perform more flops per communication, by computing 2s matrix-vector products at the beginning of each iteration of the s-step BiCGStab. This would reduce the communication cost, specifically the number of messages, by O(s) times in parallel [5]. However, this is not possible using the same formulation as BiCGStab. Therefore, at the beginning of each s-step iteration, one computes $P_{2s+1}$ and $R_{2s}$, the Krylov matrices corresponding to the $\mathcal{K}_{2s+1}(A, p_m)$ and $\mathcal{K}_{2s}(A, r_m)$ bases respectively, where m = 0, s, 2s, 3s, ... Then, by Equation (8), $p_{m+j}, r_{m+j}$,

and $x_{m+j} - x_m$ can be defined as the product of the basis vectors and a span vector, for j = 0, ..., s,

$$ p_{m+j} = [P_{2s+1}, R_{2s}]a_j \quad (9) $$

$$ r_{m+j} = [P_{2s+1}, R_{2s}]c_j \quad (10) $$

$$ x_{m+j} = x_m + [P_{2s+1}, R_{2s}]e_j \quad (11) $$

where $[P_{2s+1}, R_{2s}]$ is an $n \times (4s+1)$ matrix containing the basis vectors of $\mathcal{K}_{2s+1}(A, p_m)$ and $\mathcal{K}_{2s}(A, r_m)$, and $a_j, c_j$, and $e_j$ are span vectors of size $4s+1$. Note that $e_0 = 0$. As for $a_0$ and $c_0$, their definition depends on the type of computed basis.

One can compute a basis defined by a recurrence relation with three or less terms, such as a monomial, scaled monomial, Newton, or Chebyshev basis. Then, we have that

$$ AP_{2s} = P_{2s+1} T_{2s+1} \quad (12) $$

$$ AR_{2s-1} = R_{2s} T_{2s} \quad (13) $$

$$ A[P_{2s}, 0, R_{2s-1}, 0] = [P_{2s+1}, R_{2s}]T' \quad (14) $$

where $T_{2s+1}$ and $T_{2s}$ are change of basis matrices of size $(2s+1) \times (2s)$ and $(2s) \times (2s-1)$ respectively, and

$$ T' = \begin{bmatrix} [T_{2s+1} & 0] \\ [0 & T_{2s}] \end{bmatrix} $$

is a $(4s+1) \times (4s+1)$ matrix.

The definition of $T_{2s+1}$, $T_{2s}$ and eventually $T'$ depends on the chosen type of basis

$$ [P_{2s+1}, R_{2s}] = [\bar{p}_m, \bar{p}_{m+1}, \dots, \bar{p}_{m+2s}, \bar{r}_m, \bar{r}_{m+1}, \dots, \bar{r}_{m+2s-1}] $$

For example, in the monomial basis case, where $\bar{p}_m = p_m$, $\bar{r}_m = r_m$, $\bar{p}_{m+i} = A\bar{p}_{m+i-1}$, and $\bar{r}_{m+i} = A\bar{r}_{m+i-1}$ for $i > 0$, the matrices $T_{2s+1}$ and $T_{2s}$ are all zeros except the lower diagonal which is ones, i.e. $T_{2s+1}(i+1, i) = 1$ for $i = 1, \dots, 2s$. In the case of the scaled monomial basis, where $\bar{p}_m = \frac{p_m}{\|p_m\|}$ and $\bar{p}_{m+i} = \frac{A\bar{p}_{m+i-1}}{\|A\bar{p}_{m+i-1}\|}$, the matrices are defined as $T_{2s+1}(i+1, i) = \|A\bar{p}_{m+i-1}\|$ for $i = 1, \dots, 2s$, and $T_{2s}(i+1, i) = \|A\bar{r}_{m+i-1}\|$ for $i = 1, \dots, 2s-1$, and zero elsewhere.

The reformulation of BiCGStab into s-step BiCGStab starts by replacing the definitions (9)-(11) in Equations (1)-(3), and taking into consideration that for $j = 0$ to $s-1$.

$$ Ap_{m+j} = A[P_{2s+1}, R_{2s}]a_j = A[P_{2s}, 0, R_{2s-1}, 0]a_j \quad (15) $$

$$ = [P_{2s+1}, R_{2s}]T'a_j \quad (16) $$

$$ Ar_{m+j} = [P_{2s+1}, R_{2s}]T'c_j \quad (17) $$
$$A^2 p_{m+j} = [P_{2s+1}, R_{2s}](T')^2 a_j \quad (18)$$

Then we get the following,

$$e_{j+1} = e_j + \alpha_{m+j}a_j + \omega_{m+j}c_j - \omega_{m+j}\alpha_{m+j}T'a_j \quad (19)$$

$$a_{j+1} = c_{j+1} + \beta_{m+j}a_j - \beta_{m+j}\omega_{m+j}T'a_j \quad (20)$$

$$c_{j+1} = c_j - \alpha_{m+j}T'a_j - \omega_{m+j}T'(c_j - \alpha_{m+j}T'a_j) \quad (21)$$

with $e_0 = 0$, $a_0 = [1, 0_{4s}]^T$, and $c_0 = [0_{2s+1}, 1, 0_{2s-1}]$ in the case of monomial and Newton basis. Then, the span vectors, $a_j$, $c_j$, and $e_j$ are updated for $j = 1, \dots, s$, rather than $p_{m+j}$, $r_{m+j}$, and $x_{m+j}$ which are of size $n$ $4s+1$.

As for $\alpha_{m+j}$, $\delta_{m+j+1}$, and $\omega_{m+j}$, it is sufficient to replace $r_{m+j}$ and $p_{m+j}$ by definitions (10) and (9) to obtain

$$\begin{align*}
\delta_{m+j+1} &= \langle \tilde{r}_0, r_{m+j+1} \rangle = \langle g, c_{j+1} \rangle \\
\alpha_{m+j} &= \delta_{m+j} / \langle \tilde{r}_0, A p_{m+j} \rangle = \delta_{m+j} / \langle g, T' a_{j+1} \rangle \\
\omega_{m+j} &= \frac{\langle T' c_j - \alpha_{m+j} (T')^2 a_j, G c_j - \alpha_{m+j} G T' a_j \rangle}{\langle T' c_j - \alpha_{m+j} (T')^2 a_j, G T' c_j - \alpha_{m+j} G (T')^2 a_j \rangle}
\end{align*}$$

where $G = [P_{2s+1}, R_{2s}]^T[P_{2s+1}, R_{2s}]$ is a Gram-like matrix, $g = [P_{2s+1}, R_{2s}]^T \tilde{r}_0$. Then, $\beta_{m+j} = \frac{\delta_{m+j+1}}{\delta_{m+j}} = \frac{\alpha_{m+j}}{\omega_{m+j}}$.

$$20: \quad a_{j+1} = c_{j+1} + \beta_{m+j}a_j - \beta_{m+j}\omega_{m+j}t_a$$

21: **end for**

$$22: \quad p_{m+s} = [P_{2s+1}, R_{2s}]a_s, r_{m+s} = [P_{2s+1}, R_{2s}]z_s$$

$$23: \quad x_{m+s} = x_m + [P_{2s+1}, R_{2s}]e_s$$

$$24: \quad \rho_{m+s} = \langle r_{m+s}, r_{m+s} \rangle, k = k+1$$

25: **end While**

Algorithm 2 is a reformulation of BiCGStab that reduces communication where the matrix-vector multiplications are grouped at the beginning of the outer iteration, and the Gram-like matrix $G$ is computed once per outer iteration. Then, in the inner iterations, the vector operations of size $n$ are replaced by vector operations of size $4s+1$, where $4s+1 \ll n$. However, this reformulation alone is not sufficient to reduce communication.

For example, in the sequential case the basis computation should be done using the matrix powers kernel [12], where $\bar{p}_{m+j+1}$ is computed by parts that fit into cache memory, for $j = 0, \dots, 2s-1$. This reduces communication in the memory hierarchy of the processor and increases cache hits. In the parallel case, each processor fetches, at the beginning, the needed data from neighboring processors to compute its assigned part of the $2s$ vectors $\bar{p}_{m+j+1}$ without any communication, for $j = 0, \dots, 2s-1$. Similarly, we can compute the $2s-1$ vectors $\bar{r}_{m+j+1}$ for $j = 0, \dots, 2s-2$. Note that it is possible to compute the two bases simultaneously using a block version of the matrix powers kernel that computes a block of vectors without communication.

### Algorithm 2: s-step BiCGStab

**Input:** A, b, $x_0$, $m_{max}$, s, **Type of Basis**

**Output:** $x_m$: the mth approximate solution satisfying the stopping criteria

1: Let $r_0 = b - Ax_0$, $p_0 = r_0$, $\rho_0 = \langle r_0, r_0 \rangle$, and $k = 0$

2: Choose $\tilde{r}_0$ such that $\delta_0 = \langle \tilde{r}_0, r_0 \rangle \neq 0$.

3: **While** ($\sqrt{\rho_k} > \epsilon||b||_2$ and $k < \lfloor \frac{m_{max}}{s} \rfloor$) **Do**

4: Compute $P_{2s+1}$ and $R_{2s}$ depending on Type of

5: Basis, and output the diagonals of $T'$

6: $G = [P_{2s+1}, R_{2s}]^T[P_{2s+1}, R_{2s}]$ and $g = [P_{2s+1}, R_{2s}]^T \tilde{r}_0$

7: Initialize $a_0$, $e_0$, $c_0$ and set $m = k * s$

8: **for** ($j = 0$ to $s-1$) **Do**

9: $t_a = T'a_j$

10: $\alpha_{m+j} = \delta_{m+j} / \langle g, t_a \rangle$

11: $d = c_j - \alpha_{m+j}t_a$

12: $t_d = T'd$

13: $g_d = Gd$

14: $g_t = Gt_d$

15: $\omega_{m+j} = \langle t_d, g_d \rangle / \langle t_d, g_t \rangle$

16: $e_{j+1} = e_j + \alpha_{m+j}a_j + \omega_{m+j}d$

17: $c_{j+1} = c_j - \omega_{m+j}t_d - \alpha_{m+j}t_a$

18: $\delta_{m+j+1} = \langle g, c_{j+1} \rangle$

19: $\beta_{m+j} = (\delta_{m+j+1}/\delta_{m+j})(\alpha_{m+j}/\omega_{m+j})$

## 1.3 S-Step BiCGStab for Geoscience Applications

In geoscience applications, specifically in reservoir simulations, at each time step a new linear system of the form $Ax = b$ has to be solved. The difficulty and the ill-conditioning of the systems may vary throughout the simulation. However, in most cases, an iterative method and a preconditioner are chosen at the beginning of the simulation and are used for solving all the obtained linear systems. Since the obtained systems are not symmetric, Krylov subspace methods such as BiCGStab are used. Our aim is to implement a numerically stable version of s-step BiCGStab that has a similar convergence rate as BiCGStab for the reservoir simulations systems. The stability of s-step BiCGStab is related to the chosen s value and to the type of the basis.

Thus, we study the convergence of the s-step BiCGStab (Algorithm 2) method for different s values, using the monomial and Newton basis [7, 14, 15] and compare it to BiCG-Stab's convergence. We do not consider the scaled versions of the monomial and Newton basis, since this requires the computation of the norm of each vector at a time, which annihilates the possibility of avoiding communication in the matrix powers kernel. The test matrices, described in Section 3.1, are obtained from different reservoir simulations.
A sample of the obtained results is described in Section 3.2 (Tab. 3). For the well-conditioned matrices, the s-step BiCGStab with the monomial basis converges in fewer s-step iterations as s increases from 2 to 6. But for the ill-conditioned matrices, the convergence of the s-step BiCGStab with monomial basis is chaotic with respect to s. Note that, we focus on the consistency of the convergence behavior as s increases for the following reasons. First, as s increases, the communication cost per s-step BiCGStab iteration decreases, since more flops are performed per communication. Moreover, if the convergence of s-step BiCGStab improves as s increases, then the overall communication cost is decreased which should lead to a speedup in parallel implementations. Second, although we test the convergence of s-step BiCGStab on matrices obtained from reservoir simulations, our goal is the speedup obtained in the full simulation. As mentioned earlier, the obtained linear systems could be well-conditioned or ill-conditioned. However, at the beginning of the simulation a single s value is chosen without having any information on the systems to be solved. Thus, we would like to implement an s-step BiCGStab version for which the number of iterations needed for convergence decreases as s increases to some upper limit.

An alternative to the monomial basis is the Newton basis, which is known to be more stable [7] in the case of GMRES. However, in the case of s-step BiCGStab, computing a Newton basis is expensive. First, the 2s largest eigenvalues have to be computed once per matrix using some library such as ARPACK [16]. In general, the computed eigenvalues could be almost equal, which does not improve the stability of the basis. Thus, the eigenvalues are reordered in the Leja ordering as discussed in [7]. For the well-conditioned matrices obtained in geoscience simulations, using the Newton basis did not improve the convergence of s-step BiCGStab. Moreover, for the ill-conditioned matrices, finding the desired number of eigenvalues is time consuming. In addition, in geoscience simulations, we seek a relatively “cheap” and stable basis, since several linear systems are solved during the simulation (at least one per time step). For all these reasons, we will use the monomial basis and improve its numerical stability, as discussed in the next section.

## 2 ORTHONORMALIZED S-STEP BICGSTAB

The s-step BiCGStab method with the monomial basis, has an irregular convergence with respect to the s values, and converges slower than BiCGStab for some of the tested systems. This irregular and slow convergence might be due to the fact that the estimated residual used for stopping criterion is not equal to the exact residual [4]. However, in our case, the slow convergence is caused by the basis vectors which become numerically linearly dependent. One way to

improve the stability of the basis is by orthonormalizing it. We propose a new variant of s-step BiCGStab that orthonormalizes the basis vectors. This new version is referred to as orthonormalized s-step BiCGStab (Algorithm 3).

There are several ways of constructing the basis and orthonormalizing it. However, we derive the algorithm irrespective of the method used. We replace the 4s + 1 basis vectors [P<sub>2s+1</sub>, R<sub>2s</sub>] by an orthonormal basis Q<sub>4s+1</sub>. Then, the vectors p<sub>m+j</sub>, r<sub>m+j</sub>, and x<sub>m+j</sub> can be defined as follows,

$$p_{m+j} = Q_{4s+1}a_j \quad (22)$$

$$r_{m+j} = Q_{4s+1}c_j \quad (23)$$

$$x_{m+j} - x_m = Q_{4s+1}e_j \quad (24)$$

The n × (4s + 1) orthonormal matrix Q<sub>4s+1</sub> should satisfy AQ<sub>4s+1</sub>v = Q<sub>4s+1</sub>H<sub>4s+1</sub>v, where Q<sub>4s+1</sub>v ∈ K<sub>2s</sub>(A, p<sub>m</sub>) + K<sub>2s-1</sub>(A, r<sub>m</sub>) and H<sub>4s+1</sub> is a (4s + 1) × (4s + 1) upper Hessenberg matrix. Then, for j = 0 to s − 1 we get

$$A p_{m+j} = Q_{4s+1} H_{4s+1} a_j \quad (25)$$

$$A r_{m+j} = Q_{4s+1} H_{4s+1} c_j \quad (26)$$

$$A^2 p_{m+j} = Q_{4s+1} (H_{4s+1})^2 a_j \quad (27)$$

By replacing the definitions (22)-(27) in Equations (1)-(3), we get that

$$ \begin{aligned} e_{j+1} &= e_j + \alpha_{m+j}a_j + \omega_{m+j}c_j - \omega_{m+j}\alpha_{m+j}H_{4s+1}a_j \\ a_{j+1} &= c_{j+1} + \beta_{m+j}a_j - \beta_{m+j}\omega_{m+j}H_{4s+1}a_j \\ c_{j+1} &= c_j - \alpha_{m+j}H_{4s+1}a_j - \omega_{m+j}H_{4s+1}(c_j - \alpha_{m+j}H_{4s+1}a_j) \end{aligned} $$

with e<sub>0</sub> = 0. As for a<sub>0</sub> and c<sub>0</sub>, their definitions depend on the orthonormalization technique used. We will discuss this in Section 2.1.

As for α<sub>m+j</sub> and δ<sub>m+j+1</sub>, it is sufficient to replace r<sub>m+j</sub> and p<sub>m+j</sub> by definitions (22) and (23) to obtain

$$ \begin{aligned} \delta_{m+j+1} &= \langle \tilde{r}_0, r_{m+j+1} \rangle = \langle g, c_{j+1} \rangle \\ \alpha_{m+j} &= \frac{\delta_{m+j}}{\langle \tilde{r}_0, A p_{m+j} \rangle} = \frac{\delta_{m+j}}{\langle g, H_{4s+1} a_{j+1} \rangle} \end{aligned} $$

where g = Q<sub>4s+1</sub><sup>T</sup> r̃<sub>0</sub>. Similarly for ω<sub>m+j</sub>, we get

$$ \omega_{m+j} = \frac{\langle c_j - \alpha_{m+j}H_{4s+1}a_j, H_{4s+1}c_j - \alpha_{m+j}(H_{4s+1})^2 a_j \rangle}{\langle H_{4s+1}c_j - \alpha_{m+j}(H_{4s+1})^2 a_j, H_{4s+1}c_j - \alpha_{m+j}(H_{4s+1})^2 a_j \rangle} \quad (28) $$

$$ \text{since } Q_{4s+1}^T Q_{4s+1} = I. \text{ Finally, } \beta_{m+j} = \frac{\delta_{m+j+1}}{\delta_{m+j}} \frac{\alpha_{m+j}}{\omega_{m+j}}. $$
Algorithm 3 describes the orthonormalized s-step BiCG-Stab method except for the orthonormal basis construction phase, which we discuss in Section 2.1.

**Algorithm 3:** Orthonormalized s-step BiCGStab

**Input:** A, b, x₀, mmax, s, Type of Basis

**Output:** x<sub>m</sub>: the mth approximate solution satisfying the stopping criteria

1: Let $r_0 = b - Ax_0$, $p_0 = r_0$, $\rho_0 = \langle r_0, r_0 \rangle$, and $k = 0$

2: Choose $\tilde{r}_0$ such that $\delta_0 = \langle \tilde{r}_0, r_0 \rangle \neq 0$.

3: **While** ($\sqrt{\rho_k} > \epsilon \|b\|_2$ and $k < \lfloor \frac{m_{\max}}{s} \rfloor$) **Do**

4: Compute the orthonormal basis $Q_{4s+1}$ and

5: the upper Hessenberg matrix $H_{4s+1}$

6: Compute $g = Q_{4s+1}^T \tilde{r}$

7: Initialize $a_0, e_0, c_0$ and set $m = k * s$

8: **for** ($j = 0$ to $s - 1$) **Do**

9: $h_a = H_{4s+1}a_j$

10: $\alpha_{m+j} = \delta_{m+j} / \langle g, h_a \rangle$

11: $d = c_j - \alpha_{m+j}h_a$

12: $h_d = H_{4s+1}d$

13: $\omega_{m+j} = \langle d, h_d \rangle / \langle h_d, h_d \rangle$

14: $e_{j+1} = e_j + \alpha_{m+j}a_j + \omega_{m+j}d$

15: $c_{j+1} = c_j - \omega_{m+j}h_d - \alpha_{m+j}h_a$

16: $\delta_{m+j+1} = \langle g, c_{j+1} \rangle$

17: $\beta_{m+j} = (\delta_{m+j+1}/\delta_{m+j})(\alpha_{m+j}/\omega_{m+j})$

18: $a_{j+1} = c_{j+1} + \beta_{m+j}a_j - \beta_{m+j}\omega_{m+j}h_a$

19: **end for**

20: $p_{m+s} = Q_{4s+1}a_s, r_{m+s} = Q_{4s+1}z_s$

$x_{m+s} = x_m + Q_{4s+1}e_s$

21: $\rho_{m+s} = \langle r_{m+s}, r_{m+s} \rangle$, $k = k + 1$

22: **end While**

## 2.1 Construction of the Orthonormal Basis

The simplest parallelizable way to compute the orthonormal basis $Q_{4s+1}$ is to compute first $[P_{2s+1}, R_{2s}]$ using the matrix powers kernel. Then, orthonormalize it using a QR algorithm, such as the Tall and Skinny QR (TSQR) algorithm [13] that requires log ($p$) messages in parallel, where $p$ is the number of processors. In this case,

$$a_0 = U_{4s+1} \times [1, 0_{4s}]^T$$

and

$$c_0 = U_{4s+1} \times [0_{2s+1}, 1, 0_{2s-1}]^T$$

where $U_{4s+1}$ is the $(4s + 1) \times (4s + 1)$ upper triangular matrix obtained from the QR factorization of $[P_{2s+1}, R_{2s}]$. In addition,

$$A[P_{2s}, 0, R_{2s-1}, 0] = [P_{2s+1}, R_{2s}]T'$$

where $T'$ is the change of basis matrix. By replacing $[P_{2s+1}, R_{2s}]$ by $Q_{4s+1}U_{4s+1}$, obtained from the QR factorization, and by assuming that $U_{4s+1}$ is invertible, we get:

$$\begin{aligned} A Q_{4s+1} \nu &= A[P_{2s+1}, R_{2s}] U_{4s+1}^{-1} \nu \\ &= A[P_{2s}, 0, R_{2s-1}, 0] U_{4s+1}^{-1} \nu \\ &= [P_{2s+1}, R_{2s}] T' U_{4s+1}^{-1} \nu \\ &= Q_{4s+1} U_{4s+1} T' U_{4s+1}^{-1} \nu \\ &= Q_{4s+1} H_{4s+1} \nu \end{aligned}$$

$Q_{4s+1}$ is an $n \times (4s + 1)$ orthonormal matrix, $H_{4s+1} = U_{4s+1} T' U_{4s+1}^{-1}$ is a $(4s + 1) \times (4s + 1)$ upper Hessenberg matrix, and

$$Q_{4s+1} \nu \in K_{2s}(A, p_m) + K_{2s-1}(A, r_m)$$

Note that the matrix $H_{4s+1}$ is never constructed and multiplying $H_{4s+1}$ by a vector is equivalent to solving an upper triangular system and multiplying $T'$ and $U_{4s+1}$ by a vector.

However, there are two issues to take into consideration in the construction of the orthonormal basis. First, at iteration $k \ge 0$, we compute two bases $P_{2s+1}$ and $R_{2s}$, for two different subspaces $K_{2s+1}(A, p_k)$ and $K_{2s}(A, r_k)$. There is no guarantee that the two bases are linearly independent with respect to each others. In other words, the $4s + 1$ vectors obtained are not necessarily linearly independent. Moreover, the upper triangular matrix obtained from the orthonormalization of a linearly dependent set of vectors, is not invertible. Second, at iteration $k = 0$, we compute two bases of $K_{2s+1}(A, r_0)$ and $K_{2s}(A, r_0)$ subspaces, since $p_0$ is initialized to $r_0$, as in the BiCGStab and s-step BiCGStab algorithms.

A solution to the first problem, is to perform a split orthonormalization, where $P_{2s+1} = Q_{2s+1}U_{2s+1}$ and $R_{2s} = Q_{2s}U_{2s}$ are orthonormalized separately. Then, $Q_{4s+1} = [Q_{2s+1}, Q_{2s}]$ and $U_{4s+1} = \begin{pmatrix} U_{2s+1} & 0 \\ 0 & U_{2s} \end{pmatrix}$ still satisfy the relation

$$A Q_{4s+1} \nu = Q_{4s+1} H_{4s+1} \nu$$

where $H_{4s+1} = U_{4s+1} T' U_{4s+1}^{-1}$. But $Q_{4s+1}$ is not orthonormal, only $Q_{2s+1}$ and $Q_{2s}$ are orthonormal. Note that in the derivation of the orthonormalized s-step BiCGStab in Section 2, we only need that $Q_{4s+1}$ is orthonormal for the definition of $\omega_{m+j}$. Moreover, $\omega_{m+j}$ is obtained by minimizing the L2 norm of $r_{m+j+1}$. If instead we minimize the B norm of $r_{m+j+1}$, where B is an $n \times n$ matrix that satisfies $Q_{4s+1}^T B Q_{4s+1} = I_{4s+1}$, then $\omega_{m+j}$ would be defined as in Equation (28). We call this version split orthonormalized s-step BiCGStab (Algorithm 3).
<table><caption>**Algorithm 4: Split Orthonormalized s-step BiCGStab**</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>Input:</td><td>A, b, x<sub>0</sub>, m<sub>max</sub>, s, Type of Basis</td></tr><tr><td>Output:</td><td>x<sub>m</sub>: the mth approximate solution satisfying the stopping criteria</td></tr><tr><td>1:</td><td>Let r<sub>0</sub> = b - Ax<sub>0</sub>, p<sub>0</sub> = r<sub>0</sub>, ρ<sub>0</sub> = &lt;r<sub>0</sub>, r<sub>0</sub>&gt;, and k = 0</td></tr><tr><td>2:</td><td>Choose r̃<sub>0</sub> such that δ<sub>0</sub> = &lt;r̃<sub>0</sub>, r<sub>0</sub>&gt; ≠ 0.</td></tr><tr><td>3:</td><td>While (√ρ<sub>k</sub> &gt; ε||b||<sub>2</sub> and k &lt; ⌊(m<sub>max</sub>/s)⌋) Do</td></tr><tr><td>4:</td><td>Compute P<sub>2s+1</sub> and R<sub>2s</sub> depending on Type of Basis, and output the diagonals of T'</td></tr><tr><td>5:</td><td>Perform the QR factorization of P<sub>2s+1</sub> = Q<sub>2s+1</sub>U<sub>2s+1</sub></td></tr><tr><td>6:</td><td>and R<sub>2s</sub> = Q<sub>2s</sub>U<sub>2s</sub>.</td></tr><tr><td>7:</td><td>Let Q<sub>4s+1</sub> = [Q<sub>2s+1</sub>, Q<sub>2s</sub>], U<sub>4s+1</sub> = ( R<sub>2s+1</sub> 0 ),</td></tr><tr><td>8:</td><td>and H<sub>4s+1</sub> = U<sub>4s+1</sub>T'U<sub>4s+1</sub><sup>-1</sup></td></tr><tr><td>9:</td><td>Compute g = Q<sub>4s+1</sub>r̃</td></tr><tr><td>10:</td><td>Initialize a<sub>0</sub>, e<sub>0</sub>, c<sub>0</sub> and set m = k * s</td></tr><tr><td>11:</td><td>for (j = 0 to s - 1) Do</td></tr><tr><td>12:</td><td>h<sub>a</sub> = H<sub>4s+1</sub>a<sub>j</sub></td></tr><tr><td>13:</td><td>α<sub>m+j</sub> = δ<sub>m+j</sub> / &lt; g, h<sub>a</sub> &gt;</td></tr><tr><td>14:</td><td>d = c<sub>j</sub> - α<sub>m+j</sub>h<sub>a</sub></td></tr><tr><td>15:</td><td>h<sub>d</sub> = H<sub>4s+1</sub>d</td></tr><tr><td>16:</td><td>ω<sub>m+j</sub> = (d,h<sub>d</sub>) / (h<sub>d</sub>,h<sub>d</sub>)</td></tr><tr><td>17:</td><td>e<sub>j+1</sub> = e<sub>j</sub> + α<sub>m+j</sub>a<sub>j</sub> + ω<sub>m+j</sub>d</td></tr><tr><td>18:</td><td>c<sub>j+1</sub> = c<sub>j</sub> - ω<sub>m+j</sub>h<sub>d</sub> - α<sub>m+j</sub>h<sub>a</sub></td></tr><tr><td>19:</td><td>δ<sub>m+j+1</sub> = &lt; g, c<sub>j+1</sub> &gt;</td></tr><tr><td>20:</td><td>β<sub>m+j</sub> = (δ<sub>m+j+1</sub>/δ<sub>m+j</sub>)(α<sub>m+j</sub>/ω<sub>m+j</sub>)</td></tr><tr><td>21:</td><td>a<sub>j+1</sub> = c<sub>j+1</sub> + β<sub>m+j</sub>a<sub>j</sub> - β<sub>m+j</sub>ω<sub>m+j</sub>h<sub>a</sub></td></tr><tr><td>22:</td><td>end for</td></tr><tr><td>23:</td><td>p<sub>m+s</sub> = Q<sub>4s+1</sub>a<sub>s</sub>, r<sub>m+s</sub> = Q<sub>4s+1</sub>z<sub>s</sub>,</td></tr><tr><td>24:</td><td>x<sub>m+s</sub> = x<sub>m</sub> + Q<sub>4s+1</sub>e<sub>s</sub></td></tr><tr><td>25:</td><td>end While</td></tr></tbody></table>

For the second problem, starting with a $p_0 \neq r_0$, might improve the convergence of all the previously discussed s-step BiCGStab versions. One might pick a random $p_0$. However, its effect on the convergence of the method is unknown. Thus, to be consistent with the previously introduced BiCGStab versions, we choose to perform one iteration of BiCGStab before constructing the first $4s + 1$ basis vectors. The advantage is that all the information obtained in the BiCGStab iteration are used afterwards. Algorithm 5 could be considered as a “preprocessing” step, for all the s-step algorithms, where it is performed before the while loop and replacing the first two lines in Algorithm 2, Algorithm 3, or Algorithm 4. We refer to these versions as modified s-step BiCGStab, modified orthonormalized s-step BiCGStab, and modified split orthonormalized s-step BiCGStab.

<table><caption>**Algorithm 5: Choosing $p_0$ not equal to $r_0$**</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>1:</td><td>Let r<sub>0</sub> = b - Ax<sub>0</sub>, p<sub>0</sub> = r<sub>0</sub>, ρ<sub>0</sub> = &lt;r<sub>0</sub>, r<sub>0</sub>&gt;, and k = 0</td></tr><tr><td>2:</td><td>Choose r̃<sub>0</sub> such that δ<sub>0</sub> = &lt;r̃<sub>0</sub>, r<sub>0</sub>&gt; ≠ 0.</td></tr><tr><td>3:</td><td>α<sub>0</sub> = δ<sub>0</sub> / &lt; r̃<sub>0</sub>Ap<sub>0</sub> &gt;: s = r<sub>0</sub> - α<sub>0</sub>Ap<sub>0</sub>, t = As</td></tr></tbody></table>

<table><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>5:</td><td>ω<sub>0</sub> = &lt; t, s &gt; / &lt; t, t &gt;</td></tr><tr><td>6:</td><td>x<sub>0</sub> = x<sub>0</sub> + α<sub>0</sub>p<sub>0</sub> + ω<sub>0</sub>s</td></tr><tr><td>7:</td><td>r<sub>0</sub> = (I - ω<sub>0</sub>A)s</td></tr><tr><td>8:</td><td>β<sub>0</sub> = α<sub>0</sub>/(ω<sub>0</sub>δ<sub>0</sub>)</td></tr><tr><td>9:</td><td>δ<sub>0</sub> = &lt;r̃<sub>0</sub>, r<sub>0</sub>&gt;, β<sub>0</sub> = β<sub>0</sub> * δ<sub>0</sub></td></tr><tr><td>10:</td><td>p<sub>0</sub> = r<sub>0</sub> + β<sub>0</sub>(I - ω<sub>0</sub>A)p<sub>0</sub></td></tr><tr><td>11:</td><td>ρ<sub>0</sub> = &lt; r<sub>0</sub>, r<sub>0</sub> &gt;</td></tr></tbody></table>

## 3. RESULTS AND EXPECTED PERFORMANCE

In this Section, we show the convergence behavior of the newly introduced split orthonormalized s-step BiCGStab and modified split orthonormalized s-step BiCGStab, on the matrices defined in Section 3.1, and compare them to that of BiCGStab, s-step BiCGStab, and modified s-step BiCG-Stab in Section 3.2. Then we discuss the split orthonormalized s-step BiCGStab's computation and communication cost in Section 3.3 and compare it to that of s-step BiCGStab and BiCGStab.

### 3.1 Test Matrices

The study cases presented in this paper are obtained from different representative models for reservoir simulations at different time steps. Table 1 illustrates the characteristics of the models from which the square test matrices, GCS2K, CantaF3, SPE10 [17], and HIS, are generated.

Consequently, the obtained matrices have different profiles and varying degrees of difficulty. Table 2 describes the test matrices.

### 3.2 Convergence Results

We have implemented the s-step BiCGStab, and the split orthonormalized s-step BiCGStab (Algorithm 4) whereby the monomial bases $P_{2s+1}$ and $R_{2s}$ are first built and then orthonormalized separately using MKL's (Math Kernel Library) QR factorization. We have also implemented the corresponding modified versions, whereby one iteration of BiCGStab is performed before building the first $4s+1$ bases vectors. These algorithms are developed in MCG Solver [18, 19], a C++ based software package developed by IFPEN to provide linear solver algorithms for its industrial simulators in reservoir simulation, basin simulation or in engine combustion simulation.

Table 3 shows the convergence results for the s-step BiCG-Stab versions on the matrices introduced previously with tolerance $tol = 10^{-8}$ except for the SPE10 matrix ($tol = 10^{-4}$). The number of iterations needed for the convergence of s-step BiCGStab versions, referred to as s-Step Iterations (SI), is shown. For comparison reasons, we also show the Total number of Iterations (TI) of the s-step
<table><caption>TABLE 1
The test matrices.</caption><thead><tr><th rowspan="2">Matrix</th><th colspan="4">Block CSR format</th></tr><tr><th>Size</th><th>Block size</th><th>Nonzero blocks</th><th>Remark</th></tr></thead><tbody><tr><td>GCS2K</td><td>370 982</td><td>3</td><td>2 928 696</td><td>Well-conditioned</td></tr><tr><td>CantaF3</td><td>8016</td><td>3</td><td>65 202</td><td>Ill-conditioned</td></tr><tr><td>SPE10</td><td>1 094 421</td><td>2</td><td>6 421 171</td><td>Ill-conditioned</td></tr><tr><td>HIS</td><td>34 761</td><td>3</td><td>189 661</td><td>Ill-conditioned</td></tr></tbody></table>

BiCGStab versions, which is equal to the s-step iterations times s. Note that in exact arithmetics, each s-step BiCGStab iteration is equivalent to s iterations of the classical BiCGStab method. However, in terms of computation and communication cost, they are not equivalent, as discussed in the next section.

For the well-conditioned matrices such as GCS2K, s-step BiCGStab with monomial basis converges faster as s increases from 2 to 6, and the corresponding total iterations are in the same range as that of BiCGStab. However, this is not the case for the other ill-conditioned matrices. The s-step BiCGStab’s convergence is chaotic with respect to s and for some s values it requires more total iterations to converge than BiCGStab. As mentioned earlier, using a more stable basis such as the Newton basis could improve the convergence of the s-step BiCGStab for the ill-conditioned matrices. However, for the ill-conditioned matrices ARPACK was not able to find the requested number of eigenvalues in 3000 iterations, which is time consuming. Thus, we orthonormalize the bases for better numerical stability.

In the case of well-conditioned matrices such GCS2K, the computed monomial basis is already numerically stable. Thus, it is expected that the convergence will not improve much by orthonormalizing the basis or by starting with a $p_0 \neq r_0$. This is clear in Table 3, where the convergence of split orthonormalized s-step BiCGStab, modified s-step BiCGStab, and modified split orthonormalized s-step BiCGStab is in the same range as that of s-step BiCGStab.

On the other hand, the convergence behavior of the s-step BiCGStab methods for the ill-conditioned matrices varies. For SPE10, s-step BiCGStab converges in fewer s-step iterations, as s increases from 2 to 6. Then, orthonormalizing the basis and/or starting with $p_0 \neq r_0$ improves the numerical stability of the basis, leading to a better convergence. Note that orthonormalizing the basis (split ortho s-step BiCGStab) has a larger effect on convergence than starting with a $p_0 \neq r_0$ (modified s-step BiCGStab).

For the HIS matrix, s-step BiCGStab’s convergence fluctuates as s increases. Whereas, the other s-step BiCGStab

<table><caption>TABLE 2
The reservoir models.</caption><thead><tr><th>Matrix</th><th>Reservoir model</th><th>PVT-type</th><th>Scheme</th><th>Time step</th></tr></thead><tbody><tr><td>GCS2K</td><td>Fractured reservoir (dual-medium)</td><td>Bi component</td><td>Fully implicit</td><td>1</td></tr><tr><td>CantaF3</td><td>Fractured reservoir (dual-medium)</td><td>Multi component</td><td>IMPEX</td><td>11</td></tr><tr><td>SPE10</td><td>Classical 2 phase flow model</td><td>Black oil</td><td>Fully implicit</td><td>81</td></tr><tr><td>HIS</td><td>Classical 3 phase flow model</td><td>Black oil</td><td>Fully implicit</td><td>23</td></tr></tbody></table>

methods have a strictly decreasing convergence with respect to s. Starting with $p_0 \neq r_0$ (modified s-step BiCGStab) improved the convergence of s-step BiCGStab (except for $s = 3$). Moreover, the convergence results of modified s-step BiCGStab and modified split orthonormalized s-step BiCGStab are very similar for $s \ge 4$.

CantaF3 is a special case where the modified split orthonormalized s-step BiCGStab method was the only method that had a stable convergence as s increased from 2 to 5. All the other s-step BiCGStab methods had a chaotic convergence with respect to s.

Note that in some cases the total iterations of the s-step BiCGStab variants is more than that of BiCGStab. However, the communication cost of each s-step iteration is much less than that of s iterations of BiCGStab. Thus the s-step variants should still converge faster, in terms of runtime, in parallel implementations.

In general, we can say that for well-conditioned matrices, the s-step BiCGStab methods have a similar rate of convergence. However, for the ill-conditioned matrices, orthonormalizing the bases separately and starting with a $p_0 \neq r_0$ have positive effects on the stability of the basis which speeds up and stabilizes the convergence with respect to s.

As mentioned earlier, several linear system with different degrees of difficulty are solved throughout a basin or reservoir simulation using the same method (BiCGStab). Our goal is to speedup the convergence of BiCGStab by replacing it with an s-step version that reduces communication. Based on the presented convergence results, the modified split orthonormalized s-step BiCGStab method seems to be the most stable version with respect to s for both, well-conditioned and ill-conditioned matrices. The reason we focus on the convergence behavior as s increases, rather
<caption>TABLE 3</caption>

The convergence of the s-step BiCGStab, split orthonormalized s-step BiCGStab and their modified versions with monomial basis, for different s values, compared to BiCGStab. The s-Step Iterations (SI) and the Total Iterations (TI) are shown. The number of total iterations is equal to the number of s-step iterations times s.

<table><thead><tr><th rowspan="2">Matrix</th><th rowspan="2">BiCGStab</th><th rowspan="2">s</th><th colspan="2">s-Step</th><th colspan="2">Split ortho s-step</th><th colspan="2">Modified s-step</th><th colspan="2">Modified split ortho s-step</th></tr><tr><th>SI</th><th>TI</th><th>SI</th><th>TI</th><th>SI</th><th>TI</th><th>SI</th><th>TI</th></tr></thead><tbody><tr><td rowspan="5">GCS2K</td><td rowspan="5">193</td><td>2</td><td>99</td><td>198</td><td>91</td><td>182</td><td>89</td><td>178</td><td>90</td><td>180</td></tr><tr><td>3</td><td>61</td><td>183</td><td>62</td><td>186</td><td>59</td><td>177</td><td>61</td><td>183</td></tr><tr><td>4</td><td>45</td><td>180</td><td>47</td><td>188</td><td>46</td><td>184</td><td>50</td><td>200</td></tr><tr><td>5</td><td>41</td><td>205</td><td>41</td><td>205</td><td>38</td><td>190</td><td>38</td><td>190</td></tr><tr><td>6</td><td>34</td><td>204</td><td>33</td><td>198</td><td>36</td><td>216</td><td>35</td><td>210</td></tr><tr><td rowspan="5">CantaF3</td><td rowspan="5">4472</td><td>2</td><td>2813</td><td>5626</td><td>1266</td><td>2532</td><td>1980</td><td>3960</td><td>1566</td><td>3132</td></tr><tr><td>3</td><td>1065</td><td>3195</td><td>1240</td><td>3720</td><td>1924</td><td>5772</td><td>1013</td><td>3039</td></tr><tr><td>4</td><td>5510</td><td>22 040</td><td>801</td><td>3204</td><td>4824</td><td>19 296</td><td>947</td><td>3788</td></tr><tr><td>5</td><td>1228</td><td>6140</td><td>3612</td><td>18 060</td><td>2658</td><td>13 290</td><td>774</td><td>3870</td></tr><tr><td>6</td><td>662</td><td>3972</td><td>1079</td><td>6474</td><td>623</td><td>3738</td><td>1072</td><td>6432</td></tr><tr><td rowspan="5">SPE10</td><td rowspan="5">4744</td><td>2</td><td>2438</td><td>4876</td><td>1682</td><td>3364</td><td>2144</td><td>4288</td><td>1829</td><td>3658</td></tr><tr><td>3</td><td>1756</td><td>5268</td><td>1345</td><td>4035</td><td>1712</td><td>5136</td><td>1200</td><td>3600</td></tr><tr><td>4</td><td>1215</td><td>4860</td><td>752</td><td>3008</td><td>1137</td><td>4548</td><td>862</td><td>3448</td></tr><tr><td>5</td><td>1032</td><td>5160</td><td>752</td><td>3760</td><td>1021</td><td>5105</td><td>567</td><td>2835</td></tr><tr><td>6</td><td>758</td><td>4548</td><td>506</td><td>3036</td><td>856</td><td>5136</td><td>548</td><td>3288</td></tr><tr><td rowspan="5">HIS</td><td rowspan="5">4721</td><td>2</td><td>1867</td><td>3734</td><td>1778</td><td>3556</td><td>1721</td><td>3442</td><td>2571</td><td>5142</td></tr><tr><td>3</td><td>1249</td><td>3747</td><td>1288</td><td>3864</td><td>1568</td><td>4704</td><td>1158</td><td>3474</td></tr><tr><td>4</td><td>1631</td><td>6524</td><td>1213</td><td>4852</td><td>1107</td><td>4428</td><td>1100</td><td>4400</td></tr><tr><td>5</td><td>1013</td><td>5065</td><td>1007</td><td>5035</td><td>808</td><td>4040</td><td>869</td><td>4345</td></tr><tr><td>6</td><td>699</td><td>4194</td><td>624</td><td>3744</td><td>685</td><td>4110</td><td>676</td><td>4056</td></tr></tbody></table>
<caption>TABLE 4</caption>

The number of flops performed in one iteration of the sequential (modified) split orthonormalized s-step BiCGStab and s-step BiCGStab with monomial basis, and the corresponding flops for computing s iterations of BiCGStab.

<table><thead><tr><th></th><th>BiCGStab</th><th>s-Step BiCGStab</th><th>Split ortho s-step BiCGStab</th></tr></thead><tbody><tr><th>Flops</th><td>4 nnzs + 20 ns</td><td>(8 s − 2) nnz + 32 ns<sup>2</sup> + 44 sn +<br>11 n + 64 s<sup>3</sup> + 96 s<sup>2</sup> + 16 s − 2</td><td>(8 s − 2) nnz + 16 ns<sup>2</sup> + 32 sn + 11 n +<br><sup>32</sup>⁄<sub>3</sub> s<sup>3</sup> + <sup>112</sup>⁄<sub>3</sub> s<sup>2</sup> + <sup>308</sup>⁄<sub>3</sub> s − <sup>64</sup>⁄<sub>3</sub></td></tr></tbody></table>

than the convergence behavior for a given s value, is that the s value is fixed throughout the simulation. And the convergence effect of the chosen s value on the different linear systems is not known beforehand. Thus, we seek a robust method that will converge faster as s increases to some upper limit (5 or 6).

## 3.3 Computation and Communication Cost

In Table 4, the number of flops performed in one s-step iteration of orthonormalized s-step BiCGStab and s-step BiCGStab is presented, along with the number of flops performed in s iterations of BiCGStab.

In the (modified) split orthonormalized s-step BiCGStab there is a need for performing two QR factorizations of the matrices $P_{2s+1}$ and $R_{2s}$, however Gram-like matrix G is not computed. Thus, the computed flops in the (modified) split orthonormalized s-step BiCGStab is slightly less than the computed flops in the s-step BiCGStab as shown in Table 4. As discussed in [5], the only communication that occurs in the parallel implementation of s-step BiCGStab is in the construction of the basis and the matrix G. Similarly for the parallel implementation of the (modified) split orthonormalized s-step BiCGStab, only the construction of the basis and its orthonormalization using TSQR require communication. Note that both TSQR and the computation of G require $log(p)$ messages and sending $O((4s + 1)^2 log(p))$ words where p is the number of processors. Thus in terms of computed flops, sent messages and sent words, the (modified) split orthonormalized s-step BiCGStab and the s-step BiCGStab are equivalent. Note that the only difference between the modified split orthonormalized s-step BiCGStab and the split orthonormalized s-step BiCGStab is that Algorithm 5 is called once before the while loop. Hence, we may assume that the communication and computation cost of both methods is bounded by the same value.

On the other hand, the performed flops in one iteration of s-step BiCGStab and orthonormalized s-step BiCGStab is at least twice the flops performed in s iterations of BiCGStab. However, the number of sent messages in the s-step versions is reduced by a factor of $O(s)$, at the expense of increasing the number of sent words. Therefore, it is expected to obtain

speedup in the parallel implementations of the introduced s-step BiCGStab variants.

## CONCLUSION

In this paper, we have introduced the split orthonormalized s-step BiCGStab and the modified split orthonormalized s-step BiCGStab, variants of s-step BiCGStab where the basis vectors are orthonormalized. In addition, in the modified split orthonormalized s-step BiCGStab, we perform one iteration of BiCGStab to define a $p_0$ not equal to $r_0$.

We have studied the convergence behavior of the introduced methods with monomial basis, and compared it to that of the s-step BiCGStab for the matrices obtained from reservoir simulations. For almost all the tested matrices, the modified split orthonormalized s-step BiCGStab with monomial basis, converged faster than the s-step BiCGStab for $s = 2, \dots, 6$. Moreover, for ill-conditioned matrices, the modified split orthonormalized s-step BiCGStab has a similar convergence behavior as the BiCGStab method for $s = 2, \dots, 6$, unlike the s-step BiCGStab.

All the s-step BiCGStab versions send $O(s)$ times less messages than s iterations of BiCGStab. Moreover, the computation cost of the introduced variants is slightly less than that of the s-step BiCGStab. Hence, it is expected that the introduced s-step BiCGStab methods, specifically modified split orthonormalized s-step BiCGStab, will perform well in parallel on multi-core architectures.

As a future work, we would like to implement the split orthonormalized s-step BiCGStab and the modified split orthonormalized s-step BiCGStab, in parallel and compare its runtime to that of the parallel BiCGStab and s-step BiCGStab.

## REFERENCES

1 Saad Y., Schultz M.H. (1986) Gmres: a generalized minimal residual algorithm for solving nonsymmetric linear systems, *SIAM J. Sci. Stat. Comput.* **7**, 3, 856-869.

2 Hestenes M.R., Stiefel E. (1952) Methods of conjugate gradients for solving linear systems, *J. Res. Natl. Bur. Stand.* **49**, 409-436.
3 van der Vorst H.A. (1992) Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems, *SIAM J. Sci. Statist. Comput.* **13**, 2, 631-644.

4 Carson E., Knight N., Demmel J. (2011) Avoiding communication in two-sided Krylov subspace methods. *Technical Report UCB/EECS-2011-93*, EECS Department, University of California, Berkeley, August.

5 Carson E., Knight N., Demmel J. (2013) Avoiding communication in nonsymmetric Lanczos-based Krylov subspace methods, *SIAM J. Sci. Comput.* **35**, 5, S42-S61.

6 Grigori L, Moufawad S. (2013) Communication avoiding ILU0 preconditioner, *Technical Report*, ALPINES - INRIA Paris-Rocquencourt, March.

7 Hoemmen M. (2010) Communication-avoiding Krylov subspace methods, *PhD Thesis*, EECS Department, University of California, Berkeley.

8 Mohiyuddin M., Hoemmen M., Demmel J., Yelick K. (2009) Minimizing communication in sparse matrix solvers, *In Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis, SC'09*, New York, NY, USA, ACM, pp.1-12.

9 Chronopoulos A.T., Gear W. (1989) s-Step iterative methods for symmetric linear systems, *J. Comput. Appl. Math.* **25**, 2, 153-168.

10 Erhel J. (1995) A parallel GMRES version for general sparse matrices, *Electron. Trans. Numer. Anal.* **3**, 160-176.

11 Walker H.F. (1988) Implementation of the GMRES method using householder transformations, *SIAM J. Sci. Statist. Comput.* **9**, 1, 152-163.

12 Demmel J., Hoemmen M., Mohiyuddin M., Yelick K. (2008) Avoiding communication in sparse matrix computations,

In *Parallel and Distributed Processing, 2008. IPDPS 2008. IEEE International Symposium*, 14–18 April 2008, Held in Hyatt Regency Hotel in Miami, Florida, USA, pp. 1-12.

13 Demmel J., Grigori L., Hoemmen M., Langou J. (2012) Communication-avoiding parallel and sequential QR factorizations, *SIAM J. Sci. Comput.* **34**, 206-239.

14 Bai Z., Hu D., Reichel L. (1994) A Newton basis GMRES implementation, *IMA J. Numer. Anal.* **14**, 563-581.

15 Reichel L. (1990) Newton interpolation at Leja points, *BIT Numerical Mathematics* **30**, 332-346.

16 Lehoucq R., Sorensen D., Yang C. (1998) *ARPACK Users' Guide*, Society for Industrial and Applied Mathematics, Philadelphia, PA.

17 The 10th SPE Comparative Solution Project (2000) Retrieved from http://www.spe.org/web/csp/datasets/set02.htm.

18 Anciaux-Sedrakian A., Eaton J., Gratien J., Guignon T., Havé P., Preux C., Ricois O. (2015) Will GPGPUs be finally a credible solution for industrial reservoir simulators, *SPE Reservoir Simulation Symposium*, 23-25 February, Houston, Texas, USA, SPE-173223-MS. DOI: 10.2118/173223-MS.

19 Anciaux-Sedrakian A., Gottschling P., Gratien J., Guignon T. (2014) Survey on efficient linear solvers for porous media flow models on recent hardware architectures, *Oil Gas Sci. Technol.* - Rev. IFP **69**, 4, 753-766.

Manuscript submitted in December 2015

Manuscript accepted in October 2016

Published online in December 2016