# A note on the partial convergence management for the solution of symmetric linear systems with multiple right-hand sides

Luc Giraud, Yanfei Xiang

► To cite this version:

Luc Giraud, Yanfei Xiang. A note on the partial convergence management for the solution of symmetric linear systems with multiple right-hand sides. RR-9574, Inria Centre at the University of Bordeaux. 2025. hal-04922247

HAL Id: hal-04922247

https://inria.hal.science/hal-04922247v1

Submitted on 31 Jan 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Inria

A note on the partial convergence
management for the solution of
symmetric linear systems with
multiple right-hand sides

Luc Giraud, Yanfei Xiang

RESEARCH
REPORT

Nº 9574

January 2025

Project-Team Concace
<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# A note on the partial convergence management for the solution of symmetric linear systems with multiple right-hand sides

Luc Giraud*, Yanfei Xiang*

Project-Team Concace

Research Report n° 9574 — January 2025 — 23 pages

**Abstract:** We consider the solution of large sparse symmetric linear systems with multiple right-hand sides available simultaneously. Based on the partial convergence detection and management, described in IB-BGMRES [Linear Algebra Appl., 419 (2006), pp. 265-285] and the breakdown-free idea discussed in [BIT Numer. Math., 57 (2017), pp. 379-403], the block conjugate residual and block conjugate gradient methods with partial convergence management are proposed. It enable to select the directions to use for extending the search space from one iteration to the next by choosing the directions that contribute the most to the residual norms. We illustrate the numerical and computational benefits of these two novel block conjugate direction variants on a set of simple academic examples enabling reproducible experiments.

**Key-words:** BGC, BCR, partial convergence, numerical robustness

RESEARCH CENTRE
BORDEAUX - SUD-OUEST

200 avenue de la Vieille Tour
33405 Talence Cedex

<footnote>* Inria, Inria centre at the University of Bordeaux</footnote>
Une note sur la gestion de la convergence partielle dans les solveurs linéaires pour matrices symétriques avec second-membres multiples

**Résumé :** Nous considérons la solution de grands systèmes linéaires symétriques clairsemés avec plusieurs côtés droits disponibles simultanément. Sur la base de la détection et de la gestion de la convergence partielle, décrites dans IB-BGMRES [Linear Algebra Appl., 419 (2006), pp. 265-285] et l'idée sans rupture discutée dans [BIT Numer. Math., 57 (2017), pp. 379-403], les méthodes de résidu conjugué par bloc et de gradient conjugué par bloc avec gestion de la convergence partielle sont proposées. Elle permet de sélectionner les directions à utiliser pour étendre l'espace de recherche d'une itération à l'autre en choisissant les directions qui contribuent le plus aux normes résiduelles. Nous illustrons les avantages numériques et informatiques de ces deux nouvelles variantes de directions du gradient conjugué par blocs sur un ensemble d'exemples académiques simples permettant des expériences reproductibles.

**Mots-clés :** BGC, BCR, convergence partielle, robustesse numérique
# Contents

<table><tr><td>1</td><td>Introduction</td><td>4</td></tr><tr><td>2</td><td>Block conjugate residual and conjugate gradient variants</td><td>5</td></tr><tr><td>2.1</td><td>Conjugate residual and conjugate gradient for single right-hand side</td><td>5</td></tr><tr><td>2.2</td><td>Block extensions and their possible breakdowns and alternative remedies</td><td>6</td></tr><tr><td>2.3</td><td>Partial convergence management</td><td>8</td></tr><tr><td>3</td><td>Numerical experiments</td><td>10</td></tr><tr><td>3.1</td><td>Partial convergence with full rank and rank deficient set of right-hand sides</td><td>11</td></tr><tr><td>3.2</td><td>Influence of the value of the convergence threshold</td><td>14</td></tr><tr><td>3.3</td><td>Influence of the number of right-hand sides</td><td>15</td></tr><tr><td>3.4</td><td>Experiments with individual convergence threshold</td><td>15</td></tr><tr><td>3.5</td><td>Experiments with symmetric matrices</td><td>20</td></tr><tr><td>4</td><td>Concluding remarks</td><td>22</td></tr></table>

Distributed under a Creative Commons Attribution 4.0 International License
# 1 Introduction

For the solution of symmetric linear systems with multiple right-hand sides available simultaneously, we consider devising new block conjugate direction (residual/gradient) methods combine with a numerical component for managing the partial convergence of the multiple right-hand sides, originally referred to as the *inexact breakdown* mechanism in [14]. Partial convergence arises when one or several right-hand sides converge to the prescribe accuracy earlier than the others, or when a linear combination of some right-hand sides occurs. In such cases, the search block size can be reduced by retaining only a selected subset of directions typically used to extend the space. This adjustment ensures the robustness of the block conjugate direction methods while improving its computational efficiency.

In a matrix form the symmetric linear system with multiple right-hand sides writes

$$AX = B, \quad (1)$$

where $A \in \mathbb{R}^{n \times n}$ is a large sparse symmetric matrix, and $B = [b^{(1)}, b^{(2)}, \dots, b^{(p)}] \in \mathbb{R}^{n \times p}$ ($p > 1$) are the multiple right-hand sides. When $A$ is a symmetric positive definite matrix, the conjugate gradient method (CG) [7] is the standard choice. Alternatively, when $A$ is symmetric and non-singular but not necessarily positive definite, the conjugate residual method (CR) [11] becomes applicable. Given the multiple right-hand sides shown in Equation (1) are available simultaneously, only the block conjugate gradient method (BCG) [12], a block version of CG, and corresponding straightforward block extension of CR (refers as BCR) are considered in this work. Under the context of block methods, the search space is defined as the sum of the Krylov subspaces associated with each individual right-hand side, so that at each iteration the search space is enlarged by $p$ additional directions. Given the search space contains the individual Krylov subspaces, the convergence should be at least as fast as solving independently each right-hand side. Besides the block algorithm implementation enables the use of efficient BLAS-3 like computational kernels, the time to the solution is expected to be reduced. Unfortunately, these potential advantages come at the price of novel numerical difficulties. Such as the well-known *breakdown* [8, 12] issue caused by rank deficiency, which leads the columns of the block vectors become linear dependence. To be specific, when such breakdown occurs, at least one of the parameter matrices involved in the BCG or BCR is unavailable, leading these algorithms terminate early without finding a satisfactory approximation. Notably, Ji and Li developed a *breakdown-free* BCG method [8] with targeted formulations of the parameter matrices composed by an orthogonal search basis to avoid such breakdown issue. Another inherent challenge in block implementations is the *partial convergence* issue. Partial convergence arises when one or several right-hand sides converge to the prescribe accuracy earlier than the others, or when a linear combination of some right-hand sides occurs. Under the framework of the block version of the generalized minimum residual (GMRES) norm method [16] (refers as BGMRES), Robbé and Sadkane proposed the IB-BGMRES algorithm [14] with *inexact breakdown* (IB) mechanism for addressing the partial convergence issue. This IB mechanism ensures the robustness of IB-BGMRES and reduces the overall computational cost by focusing on the directions that contribute most to convergence. Additionally, it allows for the reintroduction of the abandoned information in next iteration when necessary, ensuring efficiency without sacrificing accuracy. Refer to IB-BGMRES-DR [1] and IB-BGCRO-DR [5] for some further applications of such IB mechanism to other block minimum residual norm methods with long-term recurrences. Although the numerical properties and advantages of the IB mechanism have been demonstrated in [1, 5, 14], to the best of our knowledge, related studies addressing partial convergence in the context of BCR and BCG methods with short-term recurrences have not yet been documented in the literature.

Based on the aforementioned discussions, this work introduces enhanced block conjugate direction methods (i.e., BCR and BCG) integrated with numerical mechanisms to address both the breakdown [8] and inexact breakdown [14] issues appeared when solving symmetric linear systems with multiple right-hand sides. As a result, two novel block conjugate direction methods with partial convergence detecting are proposed: IB-BCR and IB-BCG. These methods incorporate parameter matrices formulated using the breakdown-free strategy outlined in [8] while applying the inexact breakdown mechanism [14] to the residuals. This combination allows for adaptive adjustment of the search block size by retaining only a carefully selected subset of search directions, typically those contribute most to reduce the residual norm. This design enhances the robustness of the block conjugate direction methods while significantly improving their computational efficiency. Note that in the context of block methods with long-term recurrence, as discussed by Langou in [9], simply removing the converged parts may result in a loss of critical information, potentially slowing down convergence. To address this, the abandoned directions can be reintroduced in subsequent
iterations, as demonstrated in IB-BGMRES [14]. For block conjugate direction methods with short-term
recurrence, however, the directions abandoned in a given iteration cannot be reintroduced in later iterations
since there is no way to keep them without destroying the short-term recurrence structure.

The remainder of this manuscript is organized as follows. Section 2, the main focus of this work, details the developments aimed at devising novel block conjugate direction methods: IB-BCR and IB-BCG variants, equipped with breakdown handling [8, 12] and partial convergence detecting mechanism [14]. Section 3 demonstrates various numerical behaviors and features of the proposed block conjugate direction methods, with concluding remarks presented in Section 4.

The notations used in this work are illustrated as follows. The notation $\mathbb{R}$ refers to the real number field. The vectors are denoted by lowercase letters; matrices with multiple columns are described by uppercase letters. The symbol $|| \cdot ||$ denotes the Euclidean norm defaultly for both vectors and matrices. The superscript $^T$ denotes the transpose. For convenience of the algorithm illustration and presentation, some MATLAB notation and function are used. Without special note, a subscript $j$ for a vector (in the single right-hand case) or a matrix (in the block case) is used to indicate that the vector or matrix is obtained at iteration $j$. The span{$c_0, c_1, \dots, c_j$} refers to a space spanned by the vectors $c_0, c_1, \dots, c_j$. A matrix $C \in \mathbb{R}^{m \times l}$ consisting of $m$ rows and $l$ columns sometimes is denoted as $C_{m \times l}$ explicitly. The identity and null matrices of dimension $k$ are denoted respectively by $I_k$ and $0_k$ or just $I$ and $0$ when the dimension is evident from the context.

## 2 Block conjugate residual and conjugate gradient variants

In this section, we first retrospect the relationship and key properties of the conjugate gradient (CG) [7] and conjugate residual (CR) [11] algorithms for a single right-hand side in Section 2.1. Then, in Section 2.2, we extend these methods to their block versions and briefly review the breakdown [8, 12] issue arising from rank deficiency in the BCG method, which also manifests in the BCR context. In Section 2.3, we present an alternative approach using the inexact breakdown (IB) [14] mechanism to address partial convergence or linear combination issues. This section explores how the IB mechanism can be integrated into the BCR and BCG algorithms, enhancing their robustness and computational efficiency. This leads to two novel block conjugate direction methods, IB-BCR and IB-BCG, equipped with breakdown handling and partial convergence detection. We also introduce a simpler remedy to handle partial or individual convergence.

### 2.1 Conjugate residual and conjugate gradient for single right-hand side

The CG method [7] is a well-established approach for solving linear systems with symmetric positive definite coefficient matrices. When the coefficient matrix is symmetric but not necessarily positive definite, the CR algorithm [11] serves as a viable alternative. Notably, the CR algorithm can be viewed as a generalized form of the CG method, as elaborated by Hestenes in his work [6].

For a single right-hand side solution, the CR method is based on the Lanczos process to generate the orthonogonal basis for the Krylov subspace $K_n(A, r_0) = \text{span}\{r_0, Ar_0, \dots, A^{n-1}r_0\}$, in which $r_0$ denotes the residual vector associated with the initial guess. In the CR method, the residual vectors $r_j$ are A-conjugate, meaning they satisfy: $r_i^T A r_j = 0$ for $i \neq j$, while the vectors $Ap_j$ (where $p_j$ are the search directions) are orthogonal to each other [15, Section 6.8]: $(Ap_i)^T A p_j = 0$ for $i \neq j$. This property is in contrast to the CG method, where the search directions $p_j$ are A-conjugate, and the residuals $r_j$ are orthogonal to each other. Moreover when addressing symmetric systems, the CR method is mathematically equivalent to the minimal residual (MINRES) method [15] and the generalized minimum residual (GMRES) norm method [16], as all these methods minimize the 2-norm of the residual over the same subspaces. While the CG method minimizes the A-norm of the forward error. We refer the reader to [15, Algorithm 6.20] for the implementation of CR and to [17, Algorithm 3.1], [18, Algorithm 2] for its preconditioned variant depicted in Algorithm 1. We notice that in the preconditioned case the preconditioned residuals are A-conjugate and the $Ap_j$'s are M-orthogonal (refer to [10, Algorithm 3] and [13, Section 6] for more discussions) assuming that the preconditioner is symmetric positive definite as for MINRES.
**Algorithm 1 Preconditioned conjugate residual method for Ax = b**

**Require:** $A \in \mathbb{R}^{n \times n}$, the left-hand side of the linear systems, and a preconditioner $M \in \mathbb{R}^{n \times n}$, an approximation of the inverse of $A$

**Require:** $b \in \mathbb{R}^n$, the right-hand side, and $x_0 \in \mathbb{R}^n$, the initial guess

**Require:** $m$ maximum number of the iteration step

1: Compute $r_0 = b - Ax_0$, $z_0 = Mr_0$, $p_0 = z_0$

2: **for** $j = 0, 1, 2, \dots, m$ **do**

3: $\alpha_j = z_j^T A z_j / (M A p_j)^T A p_j$

4: $x_{j+1} = x_j + \alpha_j p_j$

5: $r_{j+1} = r_j - \alpha_j (A p_j)$

6: $z_{j+1} = M r_{j+1}$

7: $\beta_j = z_{j+1}^T A z_{j+1} / z_j^T A z_j$

8: $p_{j+1} = z_{j+1} + \beta_j p_j$

9: **end for**

**Ensure:** Return $x_{j+1}$, the computed solution

## 2.2 Block extensions and their possible breakdowns and alternative remedies

When addressing linear systems (1) with multiple right-hand sides, where $p > 1$, the block extensions of the conjugate direction methods, as discussed in Section 2.1, are generally considered. Because block methods define the search space as the sum of the Krylov subspaces associated with each individual right-hand side, so that at each iteration the search space is enlarged by $p$ additional directions rather than a single direction. Starting from a block initial guess $X_0 = [x_0^{(1)}, x_0^{(2)}, \dots, x_0^{(p)}] \in \mathbb{R}^{n \times p}$ and the associated block initial residual $R_0 = B - AX_0$, the block conjugate gradient (BCG) method was proposed by O'Leary in [12] to solve such systems with a symmetric positive definite coefficient matrix and multiple right-hand sides provided simultaneously. Its preconditioned version is outlined in Algorithm 2.

**Algorithm 2 Block preconditioned conjugate gradient method [12] for $AX = B$**

**Require:** $A \in \mathbb{R}^{n \times n}$ the left-hand side of the linear systems and a preconditioner $M \in \mathbb{R}^{n \times n}$ be an approximation of the inverse of $A$

**Require:** $B \in \mathbb{R}^{n \times p}$ the of right-hand-sides and $X_0 \in \mathbb{R}^{n \times p}$ the block initial guess

**Require:** $m$ maximum number of the block iteration step

1: Compute $R_0 = B - AX_0$, $Z_0 = MR_0$, $P_0 = Z_0$

2: **for** $j = 0, 1, 2, \dots, m$ **do**

3: $\alpha_j = (P_j^T A P_j)^{-1} (Z_j^T R_j)$

4: $X_{j+1} = X_j + P_j \alpha_j$

5: $R_{j+1} = R_j - A P_j \alpha_j$

6: $Z_{j+1} = M R_{j+1}$

7: $\beta_j = (Z_j^T R_j)^{-1} (Z_{j+1}^T R_{j+1})$

8: $P_{j+1} = Z_{j+1} + P_j \beta_j$

9: **end for**

**Ensure:** Return $X_{j+1}$ computed solution

The block algorithm implementation enables the use of efficient BLAS-3 like computational kernels, thus the time to the solution is expected to be reduced. Unfortunately, this potential advantage comes at the price of novel numerical difficulties. As mentioned in [12], in the implementation of BCG, the block parameter matrices $\alpha_j$ and $\beta_j$ involved in the $j$th iteration of BCG are respectively formulated as

$$ \alpha_j = (P_j^T A P_j)^{-1} (R_j^T Z_j) \in \mathbb{R}^{p \times p}, $$

and

$$ \beta_j = (R_j^T Z_j)^{-1} (R_{j+1}^T Z_{j+1}) \in \mathbb{R}^{p \times p}, \text{ with } Z_{j+1} = M R_{j+1}. $$

Exploiting the key orthogonal properties of the CR method [11] and the original ideas of O'Leary for
the preconditioned BCG method [12], the core loop of preconditioned BCR algorithm essentially reads

$$
\begin{align}
X_{j+1} &= X_j + P_j \alpha_j, \nonumber \\
R_{j+1} &= R_j - Q_j \alpha_j, \nonumber \\
Z_{j+1} &= M R_{j+1}, \nonumber \\
P_{j+1} &= Z_{j+1} + P_j \beta_j, \tag{2}
\end{align}
$$

$$
Q_{j+1} = AP_{j+1}, \tag{3}
$$

where $\alpha_j$ and $\beta_j$ are the parameter matrices to be determined such that $Z_{j+1}^T A Z_j = 0$ and $Q_{j+1}^T M Q_j = 0$. Specifically, these parameter matrices are defined by

$$
\alpha_j = (Q_j^T M Q_j)^{-1} (Q_j^T Z_j) \in \mathbb{R}^{p \times p},
$$

and

$$
\beta_j = -(Q_j^T M Q_j)^{-1} ((MQ_j)^T A Z_{j+1}) \in \mathbb{R}^{p \times p},
$$

in which some matrix-multiplications related to a $p \times p$ non-singular matrix may be involved during the practical implementation to improve the stability [12]. Based on the orthogonality properties $Z_{j+1}^T Q_j = 0$ and $Q_{j+1}^T M Q_j = 0$, as well as the core loop of the preconditioned BCR algorithm, we can derive an additional orthogonality property $Z_{j+1}^T A Z_j = 0$ as follow:

$$
\begin{align*}
Z_{j+1}^T A Z_j &= Z_{j+1}^T (Q_j - Q_{j-1} \beta_{j-1}) = -Z_{j+1}^T Q_{j-1} \beta_{j-1} = -R_{j+1}^T M Q_{j-1} \beta_{j-1} \\
&= (Q_j \alpha_j - R_j)^T M Q_{j-1} \beta_{j-1} = \alpha_j^T Q_j^T M Q_{j-1} \beta_{j-1} - Z_j^T Q_{j-1} \beta_{j-1} = 0.
\end{align*}
$$

This leads to the fact that BCR generates the jth approximate solution $X_j$ such that the Euclidean norm of the corresponding block residual $R_j$ is minimized over the increasing subspaces $X_0 + K_j(MA, MR_0)$ with $K_j(MA, MR_0) = \oplus_{\ell=1}^p K_j(MA, Mr_0^{(\ell)})$, i.e.,

$$
\| R_j \| = \min_{X_j \in X_0 + K_j(MA, MR_0)} \| B - A X_j \|,
$$

which is the same as the case for the block minimum residual method (Block MINRES) [12, Section 3].

In BCG, if the matrices $(P_j^T A P_j)$ or $(R_j^T Z_j)$ are singular this leads to the so-called *breakdown* issue [12] because the algorithm terminates early without finding a satisfactory approximate solution. Such a situation occurs, when some rank deficiency appears simultaneously in $R_j$ and $P_j$ as shown in [8, Proposition 2.1]. The loss of rank in $R_j$ indicates that either some solutions have converged or some linear combinations of the solutions have converged in the search space at the $j^{th}$ iteration. The same situation might occur in BCR with the parameter matrix $(Q_j^T M Q_j)$ that can become singular if some rank deficiency appear simultaneously in $Q_j$ and $R_j$ as shown in Proposition 2.1.

**Proposition 2.1.** For the preconditioned BCR algorithm, suppose A and M are non-singular symmetric matrices, $R_j$ be an $n \times p$ residual matrix with rank $p_j$ ($p_j \le p$) at the j-th iteration, then, we have

$$
\mathrm{rank}(P_j) = \mathrm{rank}(Q_j) = \mathrm{rank}(Z_j) = \mathrm{rank}(R_j) = p_j, \quad (4)
$$

where $\mathrm{rank}(\cdot)$ denotes the rank of a matrix.

*Proof.* From the core loop of the preconditioned BCR algorithm, we have $Z_j = MR_j$ and $Q_j = AP_j$. Because A and M are nonsingular, we have $\mathrm{rank}(Z_j) = \mathrm{rank}(R_j) = p_j$ and $\mathrm{rank}(P_j) = \mathrm{rank}(Q_j)$.

In the sequel we will use twice the similar arguments that is, if $C \in \mathbb{R}^{n \times n}$ is a symmetric non-singular matrix and $E \in \mathbb{R}^{n \times p}$ then

$$
\mathrm{rank}(E^T C E) = \mathrm{rank}(E). \tag{5}
$$

Because C is a symmetric matrix, it is diagonalizable in an orthonormal basis, that is $C = WDW^T$ where the columns of the unitary matrix W are the unitary eigenvectors and D the diagonal of the eigenvalues. Consequently $E^T C E = E^T W D W^T E = \tilde{E}^T D \tilde{E}$ with $\tilde{E} = W^T E$ and $\mathrm{rank}(\tilde{E}) = \mathrm{rank}(E)$. If we note $\tilde{E} = U \Sigma V^T$ the SVD decomposition of $\tilde{E}$ with U and V unitary matrices we have

$$
\mathrm{rank}(E^T C E) = \mathrm{rank}(V \Sigma^T U^T D U \Sigma V^T) = \mathrm{rank}(\Sigma^T D \Sigma) = \mathrm{rank}(\Sigma) = \mathrm{rank}(\tilde{E}) = \mathrm{rank}(E).
$$
From Equation (2)-(3), the block vector $Q_j$ is given by

$$Q_j = AZ_j + Q_{j-1}\beta_{j-1}. \qquad (6)$$

Left multiplying (6) by $Q_j^T M$ on both sides, we get

$$Q_j^T M Q_j = Q_j^T M A Z_j + Q_j^T M Q_{j-1} \beta_{j-1} = Q_j^T M A Z_j \qquad (7)$$

by definition of the parameter matrices $Q_j^T M Q_{j-1} = 0$. Because $M$ is symmetric, then using Equation (7) and (5) as well as some basic properties of matrix rank, we have

$$\mathrm{rank}(Q_j) = \mathrm{rank}(Q_j^T M Q_j) = \mathrm{rank}(Q_j^T M A Z_j) \le \mathrm{rank}(Z_j) = \mathrm{rank}(R_j). \qquad (8)$$

On the other hand, left multiplying (6) by $Z_j^T$ on both sides, we get

$$Z_j^T Q_j = Z_j^T A Z_j + Z_j^T Q_{j-1} \beta_{j-1} = Z_j^T A Z_j \qquad (9)$$

by the orthogonal property $Z_j^T Q_{j-1} = 0$.

Then, apply Equation (9) and (5) as well as some basic properties of matrix rank, we have

$$\mathrm{rank}(Z_j) = \mathrm{rank}(Z_j^T A Z_j) = \mathrm{rank}(Z_j^T Q_j) \le \mathrm{rank}(Q_j) = \mathrm{rank}(P_j). \qquad (10)$$

Based on relation (8) and (10), we eventually have $\mathrm{rank}(P_j) = \mathrm{rank}(Q_j) = \mathrm{rank}(Z_j) = \mathrm{rank}(R_j)$. $\square$

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

To overcome the singularity, Ji and Li proposed the breakdown-free BCG method [8], where the block search directions $P_j$ are replaced by $\tilde{P}_j$, an orthonormal basis of the space spanned by $P_j$, which can be computed by considering the left singular vectors of $P_j$ associated with non-zero singular values. Although, this alternative makes sense in exact arithmetic, it can be somehow relaxed when a prescribed accuracy is required. In that situation, the partial convergence detection introduced in [14] can be adapted to the two block algorithms to monitor the size block according to a prescribed threshold for the individual convergence criterion based on backward error.

## 2.3 Partial convergence management

In this section, we investigate the application of the partial convergence management, referred to as *Inexact Breakdown* (IB) mechanism [14], to effectively tackle the challenges of variable convergence speed for the different right)hand sides. This IB mechanism [14] was initially devised to address partial convergence in the block minimum residual norm methods with long-term recurrences, and has been further extended in [1, 5] to the block versions of GMRES and the generalized conjugate residual with inner orthogonalization (GCRO) [19] algorithms. Since BCR is also a type of minimum residual norm method, the IB mechanism can be adapted to enhance its robustness and efficiency in addressing these numerical challenges. However, two main differences exist. First, from a numerical perspective, since BCR involves short-term recurrences, the directions abandoned at a given iteration cannot be reintroduced later, as is possible in methods with long-term recurrences [1, 5, 14] to avoid the loss of critical information [9]. This limitation may make its application to the BCR method be less efficient compared to methods discussed in [1, 5, 14], as there is no way to retain these directions without disrupting the short-term recurrence structure of BCR. Second, from a computational point of view, the algorithm does not compute a QR-factorization of the residual block, contrarily to block GMRES and GCRO where it is a byproduct of the least squares problem solution (refer to [5, Section 4] for more details). Consequently, the partial convergence policies will require the reduced singular value decomposition (SVD) calculation of a tall and skinny matrix, namely of the residual block $R_j$ scaled by a certain diagonal matrix $D_\epsilon$ that depends on the selected stopping criterion and convergence threshold $\epsilon$, which could be described as the following form:

$$R_j D_\epsilon = \mathbb{U}_{1,L} \Sigma_1 \mathbb{V}_{1,R}^T + \mathbb{U}_{2,L} \Sigma_2 \mathbb{V}_{2,R}^T \text{ with } \sigma_{\mathrm{min}}(\Sigma_1) \ge \tau > \sigma_{\mathrm{max}}(\Sigma_2),$$

where $D_\epsilon = \epsilon^{-1} \mathrm{diag}(\|b^{(1)}\|^{-1}, \dots, \|b^{(p)}\|^{-1}) \in \mathbb{R}^{p \times p}$, with $\tau = 1$ the prescribed IB-threshold. The vectors $(\mathbb{U}_{1,L}, \mathbb{U}_{2,L})$ and $(\mathbb{V}_{1,R}, \mathbb{V}_{2,R})$ are the left and right singular vectors of $R_j D_\epsilon$, respectively. Once a partial convergence is detected, all the calculations but the update of the solution and residual blocks are
performed on blocks of lower column dimension, which lowers the number of matrix-vector products and preconditioning applications. Let us denote

$$[U_j^L, W_j] = \text{SpaceExpansion}(R_j, \epsilon), \quad (11)$$

where $U_j^L \in \mathbb{R}^{n \times p_j}$ are the $p_j$ left singular vectors computed by the selected partial convergence detection mechanism and $W_j = (U_j^L)^T R_j \in \mathbb{R}^{p_j \times p}$ the components of the residuals in the space spanned by $U_j^L$, that are the new directions to be added to the search space. The algorithm is depicted in Algorithm 3, where a bar notation is used to indicate that the dimensions of these matrices may reduce in case of partial convergence detection.

**Algorithm 3** Block preconditioned Conjugate Residual method with partial convergence detection (or Inexact-Breakdown) mechanism — IB-BCR

**Require:** $A \in \mathbb{R}^{n \times n}$ the left-hand side of the linear systems and a preconditioner $M \in \mathbb{R}^{n \times n}$ be an approximation of the inverse of A

**Require:** $B \in \mathbb{R}^{n \times p}$ the block of right-hand sides and $X_0 \in \mathbb{R}^{n \times p}$ the block initial guess

**Require:** $m$ maximum number of the block iteration step and the maximum number of matrix-vector products is set to be $\max M vps \in \mathbb{N}^+$

**Require:** $\epsilon > 0$ a threshold for the selected backward error used in stopping criterion

1: $U_0^L, W_0 = \text{SpaceExpansion}(B - AX_0, \epsilon)$ and $\bar{Z}_0 = MU_0^L$

2: $\bar{P}_0 = \bar{Z}_0$

3: **for** $j = 0, 1, 2, \dots, m$ **do**

4: $\bar{Q}_j = A\bar{P}_j$

5: $\bar{\alpha}_j = (\bar{Q}_j^T M \bar{Q}_j)^{-1}(\bar{Q}_j^T \bar{Z}_j)$

6: $X_{j+1} = X_j + \bar{P}_j \bar{\alpha}_j W_j$

7: $R_{j+1} = R_j - \bar{Q}_j \bar{\alpha}_j W_j$

8: **if** the stopping criterion related to $\epsilon$ or $\max M vps$ is met **then**

9: **return** $X_{j+1}$

10: **else**

11: $[U_{j+1}^L, W_{j+1}] = \text{SpaceExpansion}(R_{j+1}, \epsilon)$

12: $\bar{Z}_{j+1} = MU_{j+1}^L$

13: $\bar{\beta}_j = -(\bar{Q}_j^T M \bar{Q}_j)^{-1}((M \bar{Q}_j)^T A \bar{Z}_{j+1})$

14: $\bar{P}_{j+1} = \bar{Z}_{j+1} + \bar{P}_j \bar{\beta}_j$

15: **end if**

16: **end for**

**Ensure:** Return $X_{j+1}$

We notice that Algorithm 3 can be accommodated to implement a crude search space expansion that simply discards the columns of the block that correspond to individual solutions that have converged. In that case, $U_j^L$ consists of the columns of the residual that have not yet converged, and $W_j \in \mathbb{R}^{p \times p}$ is a diagonal matrix with entry equal to 0 when the corresponding right-hand side has converged. This somewhat naive alternative algorithm will be referred to as IC-BCR, for Individually Converged BCR.

Although the CG method [6] and its block version BCG [12] minimize the A-norm of the forward error, their most commonly used stopping criterion still relies on a backward error that is mostly based on the residual norm. This observation motivates us to adapt BCG with the previous developed partial convergence mechanism in the context of BCR. By exploiting the algorithmic resemblance between BCG and BCR, the development of corresponding IB-BCG counterpart is fairly straightforward and depicted in Algorithm 4. We conclude that the IB-BCG algorithm with short-term recurrences might lack efficiency and robustness for two main reasons. First, similar to the BCR case, it may be less efficient compared to algorithms with long-term recurrences, such as those discussed in [1, 5, 14], since the abandoned directions cannot be reintroduced in later iterations without disrupting the short-term recurrence structure. Second, it may lack robustness because BCG minimizes the A-norm of the errors, and therefore controlling the search space expansion through heuristics based on residuals does not fully align with this optimal numerical feature. Finally, the corresponding IC (Individually Converged) variant can also be considered; of does not
suffer from the two previously possible flaws of IB-BCG as it mostly consists in stopping iterating on the
iterates that have converged.

**Algorithm 4** Block preconditioned Conjugate Gradient method with partial convergence (or Inexact-Breakdown detection) mechanism — IB-BCG

**Require:** $A \in \mathbb{R}^{n \times n}$ the left-hand side of the linear systems and a preconditioner $M \in \mathbb{R}^{n \times n}$ that is an approximation of the inverse of $A$

**Require:** $B \in \mathbb{R}^{n \times p}$ the block of right-hand-sides and $X_0 \in \mathbb{R}^{n \times p}$ the block initial guess

**Require:** $m$ maximum number of the block iteration step and the maximum number of matrix-vector products ($\#mvps$) is set to be $\max Mvps \in \mathbb{N}^+$

**Require:** $\varepsilon > 0$ a threshold for the selected backward error used in stopping criterion

1: $U_0^L, W_0 = \text{SpaceExpansion}(B - AX_0, \varepsilon)$ and $\bar{Z}_0 = MU_0^L$

2: $\bar{P}_0 = (\bar{Z}_0)$

3: **for** $j = 0, 1, 2, \dots, m$ **do**

4: &nbsp;&nbsp;&nbsp;&nbsp; $\bar{Q}_j = A\bar{P}_j$

5: &nbsp;&nbsp;&nbsp;&nbsp; $\bar{\alpha}_j = (\bar{P}_j^T\bar{Q}_j)^{-1}(\bar{P}_j^T\bar{R}_j)$

6: &nbsp;&nbsp;&nbsp;&nbsp; $X_{j+1} = X_j + \bar{P}_j\bar{\alpha}_jW_j$

7: &nbsp;&nbsp;&nbsp;&nbsp; $R_{j+1} = R_j - \bar{Q}_j\bar{\alpha}_jW_j$

8: &nbsp;&nbsp;&nbsp;&nbsp; **if** the stopping criterion related to $\varepsilon$ or $\max Mvps$ is met **then**

9: &nbsp;&nbsp;&nbsp;&nbsp; **return** $X_{j+1}$

10: **else**

11: &nbsp;&nbsp;&nbsp;&nbsp; $[U_{j+1}^L, W_{j+1}] = \text{SpaceExpansion}(R_{j+1}, \varepsilon)$

12: &nbsp;&nbsp;&nbsp;&nbsp; $\bar{Z}_{j+1} = MU_{j+1}^L$

13: &nbsp;&nbsp;&nbsp;&nbsp; $\bar{\beta}_j = -(\bar{P}_j^T\bar{Q}_j)^{-1}(\bar{Q}_j^T\bar{Z}_{j+1})$

14: &nbsp;&nbsp;&nbsp;&nbsp; $\bar{P}_{j+1} = (\bar{Z}_{j+1} + \bar{P}_j\bar{\beta}_j)$

15: **end if**

16: **end for**

**Ensure:** Return $X_{j+1}$ for approximation of the linear systems

We illustrate in the next sections about the complementary numerical performances of the proposed IB-BCR and IB-BCG variants with breakdown-free [8] and partial convergence detecting mechanism [14], as well as their naive IC-variants for individual convergence detection.

3 Numerical experiments

Numerical experiments are carried out on a set of symmetric positive definite (SPD) and symmetric but not positive definite matrices from the University of Florida Sparse Matrix Collection [4]. The main features of these SPD and symmetric matrices are respectively described in Table 1 and Table 2. In order to illustrate the numerical benefits of the proposed IB-BCR (shown in Algorithm 3) and IB-BCG (shown in Algorithm 4) methods, the performance of IB-BCR and IB-BCG are respectively evaluated in comparison with that of the classical BCR and classical BCG and their IC variants in terms of the consumed number of matrix-vector products (referred to as #mvps) and the consumed block iteration steps (referred to as #iter). Given the error related with true solution of linear system is generally unavailable in practice, we compare the convergence histories of the related methods (i.e., BCR and BCG as well as their IB and IC variants) in terms of the smallest and largest backward errors related to Euclidean norm of the residuals at each #mvps within the following sections, in which the stopping criterion is that the *p* individual normwise backward errors satisfy $\eta_{b(i)}(x_j^{(i)}) = \frac{\|b^{(i)} - Ax_j^{(i)}\|}{\|b^{(i)}\|} < \varepsilon$ ($i = 1, \dots, p$) with respect to the approximate solution $x_j^{(i)}$, or #mvps exceeds the allowed maximum number of matrix-vector products (referred to as $\max Mvps$).

In the default setting of the experiments, the block initial guess is set to be $0 \in \mathbb{R}^{n \times p}$, where *p* represents the number of the right-hand sides. The multiple right-hand sides $B = \text{randn}(n, p) = [b^{(1)}, b^{(2)}, \dots, b^{(p)}] \in \mathbb{R}^{n \times p}$ are composed of *p* linearly independent vectors containing pseudo-random values drawn from the standard normal distribution (using the same seed when comparing these block methods). The search space expansion policy used in conjunction with the partial con-
vergence detecting is based on the backward error $\eta_b$. Without special notes, the $maxMvps$ is set to be $5000 \times p$ for each solver run, the convergence threshold is $\varepsilon = 10^{-8}$. For all the experiments involving SPD matrices, we consider the preconditioned BCR and BCG variants, where an incomplete Cholesky factorization is employed by default as the preconditioner. For the symmetric but not positive definite ones, no preconditioner is considered.

The experiments have been carried out in personal Linux (double precision floating point arithmetic) system by MATLAB (R2019a) with hardware setting as PC-Intel (R) Core (TM) i7-8665U CPU @ 1.90 GHz, 8 GB RAM. In order to evaluate the robustness and efficiency of the newly proposed IB-BCR and IB-BCG variants, we first investigate in Section 3.1 their numerical behavior when the set of right-hand sides is not full rank. Next we investigate their behavior when the convergence threshold $\varepsilon$ varies in Section 3.2 and when the number $p$ of right-hand sides varies in Section 3.3. In Section 3.4, we consider examples where all the solutions do not need to be computed using the same convergence threshold. Finally, we report on symmetric but not positive definite examples in Section 3.5.

<table><caption>Table 1: Main characteristics of the symmetric positive definite matrices</caption><thead><tr><th>Name</th><th>n</th><th>Nonzero</th><th>Origin</th><th>Cond. number</th></tr></thead><tbody><tr><td>apache1</td><td>80,800</td><td>542,184</td><td>Stru. Prob.</td><td></td></tr><tr><td>bcsstk15</td><td>3,948</td><td>117,816</td><td>Stru. Prob.</td><td>6.53e+09</td></tr><tr><td>bcsstk16</td><td>4,884</td><td>290,378</td><td>Stru. Prob.</td><td>4.94e+09</td></tr><tr><td>bcsstk17</td><td>10,974</td><td>428,650</td><td>Stru. Prob.</td><td>1.29e+10</td></tr><tr><td>bcsstk18</td><td>11,948</td><td>149,090</td><td>Stru. Prob.</td><td>3.45e+11</td></tr><tr><td>bundle1</td><td>10,581</td><td>770,811</td><td>Comp. Grap./Vis. Prob.</td><td>1.00e+03</td></tr><tr><td>cbuckle</td><td>13,681</td><td>676,515</td><td>Stru. Prob.</td><td>3.29e+07</td></tr><tr><td>crankseg_1</td><td>52,804</td><td>10,614,210</td><td>Stru. Prob.</td><td></td></tr><tr><td>crankseg_2</td><td>63,838</td><td>14,148,858</td><td>Stru. Prob.</td><td></td></tr><tr><td>gridgena</td><td>48,962</td><td>512,084</td><td>Opti. Prob.</td><td></td></tr><tr><td>gyro</td><td>17,361</td><td>1,021,159</td><td>Model Redu. Prob.</td><td>1.09e+09</td></tr><tr><td>Kuu</td><td>7,102</td><td>340,200</td><td>Stru. Prob.</td><td>1.57e+04</td></tr><tr><td>s1rmq4m1</td><td>5,489</td><td>262,411</td><td>Stru. Prob.</td><td>1.81e+06</td></tr><tr><td>s1rmt3m1</td><td>5,489</td><td>217,651</td><td>Stru. Prob.</td><td>2.54e+06</td></tr><tr><td>s2rmq4m1</td><td>5,489</td><td>263,351</td><td>Stru. Prob.</td><td>1.77e+08</td></tr><tr><td>s2rmt3m1</td><td>5,489</td><td>217,681</td><td>Stru. Prob.</td><td>2.49e+08</td></tr><tr><td>s3rmq4m1</td><td>5,489</td><td>262,943</td><td>Stru. Prob.</td><td>1.76e+10</td></tr><tr><td>s3rmt3m1</td><td>5,489</td><td>217,669</td><td>Stru. Prob.</td><td>2.48e+10</td></tr></tbody></table>

<table><caption>Table 2: Main characteristics of the symmetric but not positive definite matrices</caption><thead><tr><th>Name</th><th>n</th><th>Nonzero</th><th>Origin</th><th>Cond. number</th></tr></thead><tbody><tr><td>benzene</td><td>8,219</td><td>242,669</td><td>T/QC Prob.</td><td>1.45e+03</td></tr><tr><td>rail_5177</td><td>5,177</td><td>35,185</td><td>Model Redu. Prob.</td><td>5.33e+04</td></tr><tr><td>saylr4</td><td>3,564</td><td>22,316</td><td>CFD Prob.</td><td>6.86e+06</td></tr></tbody></table>

## 3.1 Partial convergence with full rank and rank deficient set of right-hand sides

In order to illustrate the benefit of using partial convergence detection, we consider the case where the right-hand sides are not full rank. For those experiments we select the matrix Kuu with full rank right-hand sides $B = \text{randn}(n, p)$ or rank deficient one $B = [B_{\text{pre}}, B_{\text{pre} \text{randn}(p/2, p/2)}]$ with $B_{\text{pre}} = \text{randn}(n, p/2)$ and $p = 20$.
The convergence histories of BCR variants for full rank right-hand sides are displayed in the left graph of Figure 1. Several observations can be made. The convergence histories of these three solvers overlap as long as no partial convergence is detected. After this first partial convergence, the convergence rate of IB-BCR becomes faster (in terms of #mvps) than that of BCR and IC-BCR. And the largest and smallest backward errors of IB-BCR smoothly and simultaneously decreases to the target threshold. The more naive, but cost-free, variant IC-BCR exhibits a plateau for the largest backward error when the first right-hand side has converged, then followed by a super fast convergence. Finally, the breakdown-free BCR variant illustrates the drawback of a block solver without any partial or individual convergence detection, that is, many directions are introduced in the search space that enables some backward errors to go below the convergence threshold without special attention to the right-hand sides that converge the slowest. Results of the block size $p_j$ along the iterations are displayed in the right graph of Figure 1, from which we noticed that the block size remains constant for BCR, while it reduces progressively for IB-BCR and more abruptly for IC-BCR as most of the right-hand sides converge at the same iterations for that example. Corresponding convergence histories of BCG variants are shown in Figure 2, similar comments can be made.

<figure><img src="image_5.png" /><figcaption>Figure 1: BCR variants for full rank right-hand sides with $p = 20$. Left: convergence histories of the largest/smallest backward errors $\eta_{b(i)}$ as a function of the number of matrix-vector products $\#mvps$. Right: block size $p_j$ along the iterations.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 2: BCG variants for full rank right-hand sides with $p = 20$. Left: convergence histories of the largest/smallest backward errors $\eta_{b(i)}$ as a function of the number of matrix-vector products $\#mvps$. Right: block size $p_j$ along the iterations.</figcaption></figure>

For the BCR variants, we illustrate in Figure 3 the robustness introduced by the partial convergence detection mechanism in a fake and somehow extreme case where the rank of the $p$ right-hand sides is $p/2$. As it can be seen in the right graph of Figure 3, the rank deficiency is immediately detected by the IB variant that reduces to $p/2$ block size at the very first iteration. Although no real breakdown is encountered by the other two variants, their convergence is very slow due to fact that the block size remains equal to $p$ despite the rank deficiency (except for the very last iterations in the IC variant), which also
reveals some lack of robustness. The displayed backward errors (left plot of Figure 3) are computed using the norm of the iterative residual, the true residual being computed only when the iterative one meets the convergence criterion. Corresponding numerical performance in terms of the number of matrix-vector products (#mvps) and block iterations (#iter) of these three block variants are summarized in Table 3. Corresponding numerical performance of BCG variants are summarized in Table 4 and Figure 4.

<figure><img src="image_7.png" /><figcaption>Figure 3: Same case as Figure 1 but the right-hand sides are linearly dependent as $B = [B_{pre}, B_{pre} \text{ randn}(p/2, p/2)]$ with $B_{pre} = \text{randn}(n, p/2)$.</figcaption></figure>

<table><caption>Table 3: Numerical results of BCR variants in terms of #mvps and #iter for matrices Kuu with full rank and rank deficient set of right-hand sides for Section 3.1.</caption><thead><tr><th>Columns in the right-hand sides B</th><th>Method</th><th>#mvps</th><th>#iter</th></tr></thead><tbody><tr><td rowspan="3">linearly independent (Results for Figure 1)</td><td>BCR</td><td>3200</td><td>160</td></tr><tr><td>IB-BCR</td><td><b>2944</b></td><td>166</td></tr><tr><td>IC-BCR</td><td>3121</td><td>160</td></tr><tr><td rowspan="3">linearly dependent (Results for Figure 3)</td><td>BCR</td><td>4540</td><td>227</td></tr><tr><td>IB-BCR</td><td><b>2121</b></td><td>231</td></tr><tr><td>IC-BCR</td><td>4442</td><td>227</td></tr></tbody></table>

<table><caption>Table 4: Numerical results of BCG variants in terms of both #mvps and #iter for matrices Kuu with full rank and rank deficient set of right-hand sides for Section 3.1.</caption><thead><tr><th>Columns in the right-hand sides B</th><th>Method</th><th>#mvps</th><th>#iter</th></tr></thead><tbody><tr><td rowspan="3">linearly independent (Results for Figure 2)</td><td>BCG</td><td>3240</td><td>162</td></tr><tr><td>IB-BCG</td><td><b>2987</b></td><td>166</td></tr><tr><td>IC-BCG</td><td>3162</td><td>162</td></tr><tr><td rowspan="3">linearly dependent (Results for Figure 4)</td><td>BCG</td><td>2879</td><td><b>231</b></td></tr><tr><td>IB-BCG</td><td><b>2152</b></td><td>232</td></tr><tr><td>IC-BCG</td><td>2977</td><td>235</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>Figure 4: Same case as Figure 2 but the right-hand sides are linearly dependent as $B = [B_{pre}, B_{pre} \text{ randn}(p/2, p/2)]$ with $B_{pre} = \text{randn}(n, p/2)$.</figcaption></figure>

In order to be more exhaustive, we respectively report in Table 5 and Table 6 for the numerical performances of the BCR and BCG variants in terms of #mvps and #iter for all matrices listed in Table 1. Due to potential gaps between the true residual and the iterative residual, some right-hand sides may fail to converge with respect to the backward error computed using the true residual norm, even when the stopping criterion is satisfied based on the iterative residual norm. Thus we introduce the notation "*" to indicate cases where convergence is achieved according to the iterative residual norm but not with respect to the true residual norm.

<table><caption>Table 5: Numerical results of BCR variants in terms of #mvps and #iter for all matrices listed in Table 1 with right-hand sides $B = \text{randn}(n, p)$ with $p = 20$ and $\varepsilon = 10^{-8}$ for Section 3.1.</caption><thead><tr><th rowspan="2">Matrix</th><th colspan="2">#mvps</th><th colspan="2">#iter</th></tr><tr><th>BCR / IB-BCR / IC-BCR</th><th>BCR / IB-BCR / IC-BCR</th><th>BCR / IB-BCR / IC-BCR</th><th>BCR / IB-BCR / IC-BCR</th></tr></thead><tbody><tr><td>apache1</td><td>15160 / <b>13938</b> / 14758</td><td></td><td><b>758</b> / 870 / 772</td><td></td></tr><tr><td>bcssstk15</td><td>2700 / <b>2611</b> / 2668</td><td></td><td>135 / 155 / 135</td><td></td></tr><tr><td>bundle1</td><td>700 / <b>669</b> / 685</td><td></td><td>35 / 37 / 35</td><td></td></tr><tr><td>cbuckle</td><td>14460 / <b>13999</b> / 14248</td><td></td><td><b>723</b> / 749 / 727</td><td></td></tr><tr><td>crankseg_1</td><td>6200 / <b>5535</b> / 5949</td><td></td><td><b>310</b> / 333 / 314</td><td></td></tr><tr><td>crankseg_2</td><td>6980 / <b>6480</b> / 6821</td><td></td><td><b>349</b> / 370 / 351</td><td></td></tr><tr><td>gridgena</td><td>9140 / <b>8798</b> / 9067</td><td></td><td>457 / 463 / 457</td><td></td></tr><tr><td>s1rmq4m1</td><td>2760 / <b>2445</b> / 2689</td><td></td><td><b>138</b> / 157 / 139</td><td></td></tr><tr><td>s1rmt3m1</td><td>2800 / <b>2512</b> / 2739</td><td></td><td><b>140</b> / 157 / 141</td><td></td></tr><tr><td>s1rmq4m1</td><td>2760 / <b>2445</b> / 2689</td><td></td><td><b>138</b> / 157 / 139</td><td></td></tr><tr><td>s1rmt3m1</td><td>2800 / <b>2512</b> / 2739</td><td></td><td><b>140</b> / 157 / 141</td><td></td></tr></tbody></table>

## 3.2 Influence of the value of the convergence threshold

In this section, we investigate how the value of convergence threshold affects the performance and robustness of the proposed IB-BCR and IB-BCG variants with partial convergence detecting. The convergence thresholds for $\eta_b$ are set to be $\varepsilon = 10^{-1}, 10^{-2}, 10^{-5}, 10^{-12}$ and the number of right-hand sides is set to be $p = 20$. With these numerical setting, for illustration purpose, the convergence histories for the solution involving the matrix Kuu are depicted in Figure 5. The general trends are very similar for the various convergence thresholds. We observe that, for the IB variant, all the right-hand sides converge simultaneously to
<table><caption>Table 6: Numerical results of BCG variants in terms of #mvps and #iter for all matrices listed in Table 1 with right-hand sides $B = randn(n, p)$ with $p = 20$ and $\epsilon = 10^{-8}$ for Section 3.1.</caption><thead><tr><th rowspan="2">Matrix</th><th colspan="2">#mvps</th><th colspan="2">#iter</th></tr><tr><th>BCG / IB-BCG / IC-BCG</th><th>BCG / IB-BCG / IC-BCG</th><th>BCG / IB-BCG / IC-BCG</th><th>BCG / IB-BCG / IC-BCG</th></tr></thead><tbody><tr><td>apache1</td><td>15420 / <b>13969</b> / 15017</td><td></td><td><b>771</b> / 837 / 791</td><td></td></tr><tr><td>bcsstk15</td><td>2740 / <b>2614</b> / 2693</td><td></td><td>137 / 153 / 137</td><td></td></tr><tr><td>bcsstk16</td><td>1080 / <b>993</b> / 1047</td><td></td><td>54 / 55 / 54</td><td></td></tr><tr><td>bcsstk17</td><td>13540 / <b>12085</b> / 13284</td><td></td><td><b>677</b> / 740 / 683</td><td></td></tr><tr><td>bcsstk18</td><td>6540 / <b>6303</b> / 6481</td><td></td><td><b>327</b> / 343 / 329</td><td></td></tr><tr><td>bundle1</td><td>720 / <b>682</b> / 691</td><td></td><td>36 / 39 / 36</td><td></td></tr><tr><td>cbuckle</td><td>15080 / <b>14681</b> / 14920</td><td></td><td><b>754</b> / 785 / 758</td><td></td></tr><tr><td>crankseg_1</td><td>6380 / <b>5714</b> / 6117</td><td></td><td><b>319</b> / 337 / 324</td><td></td></tr><tr><td>crankseg_2</td><td>7200 / <b>6696</b> / 7055</td><td></td><td><b>360</b> / 377 / 363</td><td></td></tr><tr><td>gridgena</td><td>9180 / <b>8846</b> / 9115</td><td></td><td>459 / 462 / 459</td><td></td></tr><tr><td>s1rmq4m1</td><td>2800 / <b>2462</b> / 2720</td><td></td><td>140 / 153 / 140</td><td></td></tr><tr><td>s1rmt3m1</td><td>2840 / <b>2522</b> / 2764</td><td></td><td>142 / 150 / 142</td><td></td></tr><tr><td>s2rmq4m1</td><td>4480 / <b>3964</b> / 4399</td><td></td><td><b>224</b> / 240 / 225</td><td></td></tr><tr><td>s2rmt3m1</td><td>5760 / <b>5190</b> / 5671</td><td></td><td><b>288</b> / 318 / 289</td><td></td></tr></tbody></table>

the required accuracy, unlike the IC variant. Although the IB variant is more effective in terms of reducing the number of #mvps, both lead to the same solution quality for all the right-hand sides, while BCR and BCG compute solutions with much smaller backward error than required.

A more exhaustive set of numerical results are reported in Table 7-8 in terms of #mvps and #iter for other examples. This illustrates that the short-term recurrence induces some residual gaps that cause trouble for stringent convergence thresholds like $10^{-12}$. The general trend is that the IB and IC variants minimize #mvps while the classical BCR and BCG ones minimize #iter.

### 3.3 Influence of the number of right-hand sides

In this section, we illustrate how the number of the right-hand sides interplays with the performance of the BCR and BCG variants; we vary $p = 5, 10, 30, 40$. The numerical experiments are displayed in Table 9-10. No significant impact on the ranking of the variants can be observed. When solving for a large number of right-hand sides, it can be seen that it is preferable (when it is affordable from a memory view point) to solve all at once rather than dividing them in chunks of smaller number to be solved in sequence. For instance solving for $p = 40$ right-hand sides with the matrix apache1 by IB-BCR does not require 4 times more #mvps but 2.8 times compared with solving a sequence of 4 block systems with $p = 10$. Corresponding convergence histories are respectively described in the left and right graph of Figure 6, from which it is easy to notice that the larger $p$, the clearer the gap between the smallest and largest backward errors.

### 3.4 Experiments with individual convergence threshold

In this section, we illustrate that the partial convergence mechanism can also be adapted to cope with different individual convergence thresholds refer to as “variable accuracy” in the VA variant. Within the BCR or BCG context, it implies to change the *SpaceExpansion* function described in Equation (11) (Step 11 of Algorithm 3 or Step 11 of Algorithm 4). We illustrate this feature in Figure 7 where we consider the solution for matrix crankseg_1 with $p = 20$ right-hand sides and set the convergence threshold to $\epsilon = 10^{-4}$
<table><caption>Table 7: Numerical results of BCR variants in terms of #mvps and #iter for parts of matrices listed in Table 1 with the right-hand sides $B = \text{randn}(n, p)$, $p = 20$, $\text{maxMvps} = 5000 \times p$, and various convergence thresholds $\epsilon = 10^{-1}, 10^{-2}, 10^{-5}, 10^{-12}$ for Section 3.2.</caption><thead><tr><th rowspan="2">Matrix</th><th rowspan="2">ε</th><th colspan="2">10<sup>-1</sup></th><th colspan="2">10<sup>-2</sup></th><th colspan="2">10<sup>-5</sup></th><th colspan="2">10<sup>-12</sup></th></tr><tr><th>Method</th><th>#mvps</th><th>#iter</th><th>#mvps</th><th>#iter</th><th>#mvps</th><th>#iter</th><th>#mvps</th><th>#iter</th></tr></thead><tbody><tr><td rowspan="3">apache1</td><td>BCR</td><td>740</td><td>37</td><td>7740</td><td><strong>387</strong></td><td>12080</td><td><strong>604</strong></td><td>18960*</td><td><strong>948</strong></td></tr><tr><td>IB-BCR</td><td><strong>725</strong></td><td>38</td><td><strong>5619</strong></td><td>887</td><td><strong>10544</strong></td><td>740</td><td><strong>17630</strong>*</td><td>1021</td></tr><tr><td>IC-BCR</td><td>726</td><td>37</td><td>6082</td><td>670</td><td>11530</td><td>624</td><td>18575*</td><td>959</td></tr><tr><td rowspan="3">bcsstk18</td><td>BCR</td><td>3160</td><td>158</td><td>4500</td><td><strong>225</strong></td><td>5780</td><td><strong>289</strong></td><td>7140*</td><td>357</td></tr><tr><td>IB-BCR</td><td><strong>2025</strong></td><td>212</td><td><strong>3986</strong></td><td>317</td><td><strong>5562</strong></td><td>343</td><td><strong>6903</strong>*</td><td>367</td></tr><tr><td>IC-BCR</td><td>2501</td><td><strong>156</strong></td><td>4183</td><td>235</td><td>5672</td><td>290</td><td>7078*</td><td>357</td></tr><tr><td rowspan="3">cbuckle</td><td>BCR</td><td>3580</td><td><strong>179</strong></td><td>6960</td><td><strong>348</strong></td><td>11640</td><td><strong>582</strong></td><td>17760*</td><td><strong>888</strong></td></tr><tr><td>IB-BCR</td><td><strong>3167</strong></td><td>205</td><td><strong>6262</strong></td><td>384</td><td><strong>11036</strong></td><td>606</td><td><strong>17231</strong>*</td><td>905</td></tr><tr><td>IC-BCR</td><td>3362</td><td>180</td><td>6655</td><td>349</td><td>11429</td><td>585</td><td>17608*</td><td>890</td></tr><tr><td rowspan="3">crankseg_2</td><td>BCR</td><td>2180</td><td><strong>109</strong></td><td>2880</td><td><strong>144</strong></td><td>4960</td><td><strong>248</strong></td><td>9480*</td><td><strong>474</strong></td></tr><tr><td>IB-BCR</td><td><strong>1063</strong></td><td>520</td><td><strong>2197</strong></td><td>204</td><td><strong>4468</strong></td><td>270</td><td><strong>9004</strong>*</td><td>496</td></tr><tr><td>IC-BCR</td><td>1143</td><td>131</td><td>2642</td><td>148</td><td>4816</td><td>251</td><td>9344*</td><td>476</td></tr><tr><td rowspan="3">gridgena</td><td>BCR</td><td>220</td><td>11</td><td>5880</td><td><strong>294</strong></td><td>8500</td><td>425</td><td>9940*</td><td>497</td></tr><tr><td>IB-BCR</td><td><strong>193</strong></td><td>11</td><td>5821</td><td>366</td><td>8846</td><td>803</td><td><strong>9595</strong>*</td><td>502</td></tr><tr><td>IC-BCR</td><td>195</td><td>11</td><td><strong>5582</strong></td><td>303</td><td><strong>8411</strong></td><td>425</td><td>9849*</td><td>497</td></tr><tr><td rowspan="3">Kuu</td><td>BCR</td><td>300</td><td>15</td><td>1640</td><td><strong>82</strong></td><td>2560</td><td>128</td><td>3840</td><td>192</td></tr><tr><td>IB-BCR</td><td>279</td><td>17</td><td><strong>1225</strong></td><td>107</td><td><strong>2280</strong></td><td>136</td><td><strong>3631</strong></td><td>197</td></tr><tr><td>IC-BCR</td><td><strong>278</strong></td><td>15</td><td>1425</td><td>85</td><td>2467</td><td>128</td><td>3785</td><td>192</td></tr></tbody></table>

<table><caption>Table 8: Numerical results of BCG variants in terms of both #mvps and #iter for a selection of the of testing matrices listed in Table 1 with the right-hand sides $B = \text{randn}(n, p)$, $p = 20$ and $\text{maxMvps} = 5000 \times p$ and various convergence thresholds $\epsilon = 10^{-1}, 10^{-2}, 10^{-5}, 10^{-12}$ for Section 3.2.</caption><thead><tr><th rowspan="2">Matrix</th><th rowspan="2">ε</th><th colspan="2">10<sup>-1</sup></th><th colspan="2">10<sup>-2</sup></th><th colspan="2">10<sup>-5</sup></th><th colspan="2">10<sup>-12</sup></th></tr><tr><th>Method</th><th>#mvps</th><th>#iter</th><th>#mvps</th><th>#iter</th><th>#mvps</th><th>#iter</th><th>#mvps</th><th>#iter</th></tr></thead><tbody><tr><td rowspan="3">apache1</td><td>BCG</td><td>7000</td><td>350</td><td>8520</td><td><strong>426</strong></td><td>12340</td><td><strong>617</strong></td><td>19040*</td><td><strong>952</strong></td></tr><tr><td>IB-BCG</td><td>3072</td><td>685</td><td><strong>6670</strong></td><td>739</td><td><strong>10752</strong></td><td>718</td><td><strong>17997</strong>*</td><td>1036</td></tr><tr><td>IC-BCG</td><td><strong>1742</strong></td><td><strong>123</strong></td><td>7727</td><td>454</td><td>11788</td><td>632</td><td>18706*</td><td>965</td></tr><tr><td rowspan="3">bcsstk18</td><td>BCG</td><td>4160</td><td><strong>208</strong></td><td>4860</td><td><strong>243</strong></td><td>5940</td><td><strong>297</strong></td><td>7660*</td><td>383</td></tr><tr><td>IB-BCG</td><td><strong>3720</strong></td><td>304</td><td><strong>4396</strong></td><td>297</td><td><strong>5728</strong></td><td>329</td><td><strong>7029</strong>*</td><td><strong>372</strong></td></tr><tr><td>IC-BCG</td><td>3962</td><td>219</td><td>4728</td><td>248</td><td>5850</td><td>299</td><td>7307*</td><td>378</td></tr><tr><td rowspan="3">cbuckle</td><td>BCG</td><td>7220</td><td><strong>361</strong></td><td>8860</td><td><strong>443</strong></td><td>12580</td><td><strong>629</strong></td><td>18380*</td><td><strong>919</strong></td></tr><tr><td>IB-BCG</td><td><strong>6541</strong></td><td>372</td><td><strong>8253</strong></td><td>456</td><td><strong>12024</strong></td><td>646</td><td><strong>17857</strong>*</td><td>928</td></tr><tr><td>IC-BCG</td><td>6994</td><td>363</td><td>8683</td><td>448</td><td>12325</td><td>634</td><td>18190*</td><td>921</td></tr><tr><td rowspan="3">crankseg_2</td><td>BCG</td><td>2040</td><td><strong>102</strong></td><td>3040</td><td><strong>152</strong></td><td>5220</td><td><strong>261</strong></td><td>9700*</td><td><strong>485</strong></td></tr><tr><td>IB-BCG</td><td><strong>1257</strong></td><td>133</td><td><strong>2503</strong></td><td>180</td><td><strong>4707</strong></td><td>281</td><td><strong>9225</strong>*</td><td>499</td></tr><tr><td>IC-BCG</td><td>1659</td><td>110</td><td>2830</td><td>156</td><td>5057</td><td>265</td><td>9567*</td><td>487</td></tr><tr><td rowspan="3">gridgena</td><td>BCG</td><td>6620</td><td>331</td><td>7740</td><td>387</td><td>8540</td><td>427</td><td>10100*</td><td><strong>505</strong></td></tr><tr><td>IB-BCG</td><td><strong>5112</strong></td><td>354</td><td>7799</td><td>643</td><td><strong>8182</strong></td><td>433</td><td><strong>9690</strong>*</td><td>525</td></tr><tr><td>IC-BCG</td><td>5761</td><td><strong>316</strong></td><td><strong>7535</strong></td><td>387</td><td>8469</td><td>427</td><td>9964*</td><td>506</td></tr><tr><td rowspan="3">Kuu</td><td>BCG</td><td>1360</td><td><strong>68</strong></td><td>1780</td><td><strong>89</strong></td><td>2620</td><td><strong>131</strong></td><td>3860</td><td><strong>193</strong></td></tr><tr><td>IB-BCG</td><td><strong>846</strong></td><td>84</td><td><strong>1414</strong></td><td>102</td><td><strong>2345</strong></td><td>138</td><td><strong>3662</strong></td><td>197</td></tr><tr><td>IC-BCG</td><td>1087</td><td>72</td><td>1639</td><td>90</td><td>2535</td><td>132</td><td>3810</td><td>194</td></tr></tbody></table>
<table><caption>Table 9: Numerical results of BCR variants with different number of right-hand sides ($p = 5, 10, 30, 40$) in terms of #mvps and #iter for parts of matrices listed in Table 1 with $B = \text{randn}(n, p)$, $\text{maxMvps} = 5000 \times p$ and $\varepsilon = 10^{-8}$ for Section 3.3.</caption><thead><tr><th rowspan="2">Matrix</th><th rowspan="2">Method</th><th colspan="2">5</th><th colspan="2">10</th><th colspan="2">30</th><th colspan="2">40</th></tr><tr><th>#mvps</th><th>#iter</th><th>#mvps</th><th>#iter</th><th>#mvps</th><th>#iter</th><th>#mvps</th><th>#iter</th></tr></thead><tbody><tr><td rowspan="3">apache1</td><td>BCR</td><td>7590</td><td><b>1518</b></td><td>10800</td><td><b>1080</b></td><td>19080</td><td><b>636</b></td><td>22200</td><td><b>555</b></td></tr><tr><td>IB-BCR</td><td><b>7082</b></td><td>1641</td><td><b>9920</b></td><td>1191</td><td><b>17022</b></td><td>752</td><td><b>19784</b></td><td>717</td></tr><tr><td>IC-BCR</td><td>7532</td><td>1520</td><td>10538</td><td>1091</td><td>18435</td><td>653</td><td>21416</td><td>569</td></tr><tr><td rowspan="3">bcsstk18</td><td>BCR</td><td>4200*</td><td><b>840</b></td><td>5370*</td><td><b>537</b></td><td>7140*</td><td><b>238</b></td><td>7760*</td><td><b>194</b></td></tr><tr><td>IB-BCR</td><td><b>3963*</b></td><td>871</td><td><b>5213*</b></td><td>569</td><td><b>6813*</b></td><td>251</td><td><b>7344*</b></td><td>206</td></tr><tr><td>IC-BCR</td><td>4070*</td><td>846</td><td>5297*</td><td>540</td><td>7053*</td><td>240</td><td>7647*</td><td>195</td></tr><tr><td rowspan="3">cbuckle</td><td>BCR</td><td>6055</td><td><b>1211</b></td><td>10010</td><td><b>1001</b></td><td>17640</td><td><b>588</b></td><td>20080</td><td><b>502</b></td></tr><tr><td>IB-BCR</td><td><b>5892</b></td><td>1236</td><td><b>9573</b></td><td>1019</td><td><b>16924</b></td><td>618</td><td><b>19255</b></td><td>531</td></tr><tr><td>IC-BCR</td><td>6000</td><td>1212</td><td>9808</td><td>1003</td><td>17491</td><td>590</td><td>19857</td><td>506</td></tr><tr><td rowspan="3">crankseg_2</td><td>BCR</td><td>2765</td><td><b>553</b></td><td>4160</td><td><b>416</b></td><td>9390</td><td><b>313</b></td><td>11760</td><td><b>294</b></td></tr><tr><td>IB-BCR</td><td><b>2581</b></td><td>579</td><td><b>3961</b></td><td>431</td><td><b>8737</b></td><td>337</td><td><b>10856</b></td><td>323</td></tr><tr><td>IC-BCR</td><td>2674</td><td>559</td><td>4060</td><td>417</td><td>9195</td><td>316</td><td>11538</td><td>297</td></tr><tr><td rowspan="3">gridgena</td><td>BCR</td><td>7130</td><td><b>1426</b></td><td>8150</td><td>815</td><td>10080</td><td>336</td><td>10880</td><td>272</td></tr><tr><td>IB-BCR</td><td>7213</td><td>1548</td><td><b>8023</b></td><td>879</td><td><b>9525</b></td><td>342</td><td><b>10203</b></td><td>279</td></tr><tr><td>IC-BCR</td><td><b>7104</b></td><td>1427</td><td>8118</td><td>815</td><td>9967</td><td>336</td><td>10752</td><td>272</td></tr><tr><td rowspan="3">Kuu</td><td>BCR</td><td>1600</td><td>320</td><td>2260</td><td>226</td><td>3900</td><td><b>130</b></td><td>4440</td><td>111</td></tr><tr><td>IB-BCR</td><td><b>1476</b></td><td>329</td><td><b>2097</b></td><td>235</td><td><b>3572</b></td><td>136</td><td><b>4075</b></td><td>118</td></tr><tr><td>IC-BCR</td><td>1582</td><td>320</td><td>2224</td><td>226</td><td>3822</td><td>131</td><td>4323</td><td>111</td></tr><tr><td rowspan="3">s2rmt3m1</td><td>BCR</td><td>3910*</td><td><b>782</b></td><td>4690*</td><td><b>469</b></td><td>6240*</td><td><b>208</b></td><td>6800*</td><td><b>170</b></td></tr><tr><td>IB-BCR</td><td><b>3785*</b></td><td>868</td><td><b>4435*</b></td><td>554</td><td><b>5620*</b></td><td>248</td><td><b>5987*</b></td><td>209</td></tr><tr><td>IC-BCR</td><td>3863*</td><td>784</td><td>4634*</td><td>471</td><td>6137*</td><td>209</td><td>6688*</td><td>172</td></tr><tr><td rowspan="3">ted_B_unscaled</td><td>BCR</td><td>135*</td><td>27</td><td>270*</td><td>27</td><td>870*</td><td>29</td><td>1120*</td><td><b>28</b></td></tr><tr><td>IB-BCR</td><td>135*</td><td>27</td><td>270*</td><td>27</td><td><b>851*</b></td><td>29</td><td>1144*</td><td>29</td></tr><tr><td>IC-BCR</td><td>135*</td><td>28</td><td>270*</td><td>28</td><td>840*</td><td>29</td><td>1120*</td><td>29</td></tr></tbody></table>
<table><caption>Table 10: Numerical results of BCG variants with different number of right-hand sides, $p = 5, 10, 30, 40$, in terms of #mvps and #iter for parts of matrices listed in Table 1 with $B = \text{randn}(n, p)$, $\text{maxMvps} = 5000 \times p$ and $\epsilon = 10^{-8}$ for Section 3.3.</caption><thead><tr><th rowspan="2">Matrix</th><th rowspan="2">Method</th><th colspan="2">5</th><th colspan="2">10</th><th colspan="2">30</th><th colspan="2">40</th></tr><tr><th>#mvps</th><th>#iter</th><th>#mvps</th><th>#iter</th><th>#mvps</th><th>#iter</th><th>#mvps</th><th>#iter</th></tr></thead><tbody><tr><td rowspan="3">apache1</td><td>BCG</td><td>7730</td><td><strong>1546</strong></td><td>10900</td><td><strong>1090</strong></td><td>19440</td><td><strong>648</strong></td><td>22320</td><td><strong>558</strong></td></tr><tr><td>IB-BCG</td><td><strong>7175</strong></td><td>1639</td><td><strong>10030</strong></td><td>1168</td><td><strong>17277</strong></td><td>756</td><td><strong>20001</strong></td><td>666</td></tr><tr><td>IC-BCG</td><td>7667</td><td>1547</td><td>10667</td><td>1103</td><td>18694</td><td>667</td><td>21695</td><td>577</td></tr><tr><td rowspan="3">bcsstk18</td><td>BCG</td><td>4395</td><td><strong>879</strong></td><td>5510</td><td><strong>551</strong></td><td>7260</td><td><strong>242</strong></td><td>7880</td><td><strong>197</strong></td></tr><tr><td>IB-BCG</td><td><strong>4130</strong></td><td>895</td><td><strong>5345</strong></td><td>582</td><td><strong>6959</strong></td><td>254</td><td><strong>7473</strong></td><td>209</td></tr><tr><td>IC-BCG</td><td>4260</td><td>883</td><td>5449</td><td>553</td><td>7178</td><td>243</td><td>7773</td><td>199</td></tr><tr><td rowspan="3">cbuckle</td><td>BCG</td><td>6525</td><td>1305</td><td>10570</td><td><strong>1057</strong></td><td>18450</td><td><strong>615</strong></td><td>20920</td><td><strong>523</strong></td></tr><tr><td>IB-BCG</td><td><strong>6376</strong></td><td>1312</td><td><strong>10209</strong></td><td>1073</td><td><strong>17755</strong></td><td>633</td><td><strong>20130</strong></td><td>547</td></tr><tr><td>IC-BCG</td><td>6469</td><td><strong>1299</strong></td><td>10423</td><td>1060</td><td>18305</td><td>617</td><td>20683</td><td>525</td></tr><tr><td rowspan="3">crankseg_2</td><td>BCG</td><td>2845</td><td><strong>569</strong></td><td>4290</td><td><strong>429</strong></td><td>9780</td><td><strong>326</strong></td><td>12160</td><td><strong>304</strong></td></tr><tr><td>IB-BCG</td><td><strong>2642</strong></td><td>586</td><td><strong>4088</strong></td><td>437</td><td><strong>9053</strong></td><td>341</td><td><strong>11242</strong></td><td>321</td></tr><tr><td>IC-BCG</td><td>2755</td><td>576</td><td>4188</td><td>431</td><td>9521</td><td>328</td><td>11943</td><td>308</td></tr><tr><td rowspan="3">gridgena</td><td>BCG</td><td>7395</td><td>1479</td><td>8180</td><td>818</td><td>10110</td><td><strong>337</strong></td><td>10960</td><td>274</td></tr><tr><td>IB-BCG</td><td><strong>7222</strong></td><td>1493</td><td><strong>7997</strong></td><td>820</td><td><strong>9593</strong></td><td>341</td><td><strong>10266</strong></td><td>278</td></tr><tr><td>IC-BCG</td><td>7358</td><td><strong>1478</strong></td><td>8151</td><td>818</td><td>10023</td><td>338</td><td>10820</td><td>274</td></tr><tr><td rowspan="3">Kuu</td><td>BCG</td><td>1625</td><td><strong>325</strong></td><td>2290</td><td><strong>229</strong></td><td>3960</td><td>132</td><td>4480</td><td><strong>112</strong></td></tr><tr><td>IB-BCG</td><td><strong>1517</strong></td><td>335</td><td><strong>2136</strong></td><td>234</td><td><strong>3622</strong></td><td>137</td><td><strong>4115</strong></td><td>117</td></tr><tr><td>IC-BCG</td><td>1614</td><td>326</td><td>2262</td><td>230</td><td>3860</td><td>132</td><td>4363</td><td>113</td></tr><tr><td rowspan="3">s2rmt3m1</td><td>BCG</td><td>3955</td><td><strong>791</strong></td><td>4750</td><td><strong>475</strong></td><td>6360</td><td><strong>212</strong></td><td>6960</td><td><strong>174</strong></td></tr><tr><td>IB-BCG</td><td><strong>3772</strong></td><td>821</td><td><strong>4434</strong></td><td>511</td><td><strong>5696</strong></td><td>241</td><td><strong>6092</strong></td><td>201</td></tr><tr><td>IC-BCG</td><td>3909</td><td>793</td><td>4693</td><td>477</td><td>6259</td><td>213</td><td>6837</td><td>177</td></tr><tr><td rowspan="3">ted_B</td><td>BCG</td><td>135</td><td>27</td><td>270</td><td>27</td><td>870</td><td>29</td><td>1160</td><td>29</td></tr><tr><td>IB-BCG</td><td>135</td><td>27</td><td>270</td><td>27</td><td>864</td><td>33</td><td>1134</td><td>30</td></tr><tr><td>IC-BCG</td><td>135</td><td>27</td><td>270</td><td>27</td><td><strong>844</strong></td><td>29</td><td><strong>1128</strong></td><td>29</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>Figure 5: Convergence history for solving linear systems built by Kuu ($B = \text{randn}(n, p)$, $p = 20$ and $\text{maxMvps} = 5000 \times p$) with different values of the convergence threshold.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 6: Comparison convergence history of the largest/smallest backward errors $\eta_{b(i)}$ of the BCR variants with different number of right-hand sides $p$ by solving matrix apache1 with $B = \text{randn}(n, p)$, $m_d = 5000 \times p$ and $\varepsilon = 10^{-8}$.</figcaption></figure>

for the first $p/2$ right-hand sides and $\varepsilon = 10^{-8}$ for the last $p/2$ ones. It can be seen that numerically it works and enables some saving in terms of #mvps. More numerical results in terms of #mvps and block #iter
are summarized in Table 11-12, which exhibits a moderate positive benefit of this VA variant that was very effective in the block GCRO context presented in [5, Section 5.3].

<figure><img src="image_11.png" /><figcaption>Figure 7: Histories of the largest/smallest backward errors $\eta_{b(i)}$ at each #mvps for matrix crankseg_1 with convergence threshold equals to $10^{-4}$ for the first p/2 right-hand sides and $10^{-8}$ for the last p/2 ones (here p = 20).</figcaption></figure>

<table><caption>Table 11: Numerical results of BCR variants for Section 3.4 in terms of #mvps and #iter with $B = \text{randn}(n, p)$, $p = 20$ and $\text{maxMvps} = 5000 \times p$.</caption><thead><tr><th rowspan="2">Matrix</th><th colspan="2">#mvps</th><th colspan="2">#iter</th></tr><tr><th>BCR / IB-BCR / IB-BCR-VA</th><th>BCR / IB-BCR / IB-BCR-VA</th><th>BCR / IB-BCR / IB-BCR-VA</th><th>BCR / IB-BCR / IB-BCR-VA</th></tr></thead><tbody><tr><td>cbuckle</td><td>14460 / <b>14018</b> / 14049</td><td></td><td><b>723</b> / 749 / 1069</td><td></td></tr><tr><td>crankseg_1</td><td>6200 / 5554 / <b>5344</b></td><td></td><td><b>310</b> / 333 / 482</td><td></td></tr><tr><td>crankseg_2</td><td>6980 / <b>6499</b> / 6943</td><td></td><td><b>349</b> / 369 / 845</td><td></td></tr><tr><td>Kuu</td><td>3200 / <b>2963</b> / 3609</td><td></td><td><b>160</b> / 166 / 423</td><td></td></tr><tr><td>shallow_water2</td><td>560 / 560 / <b>430</b></td><td></td><td>28 / 28 / 28</td><td></td></tr><tr><td>ted_B_unscaled</td><td>540* / 540* / <b>415</b>*</td><td></td><td>27 / 27 / 27</td><td></td></tr></tbody></table>

<table><caption>Table 12: Numerical results of BCG variants for Section 3.4 in terms of #mvps and #iter with $B = \text{randn}(n, p)$, $p = 20$ and $\text{maxMvps} = 5000 \times p$.</caption><thead><tr><th rowspan="2">Matrix</th><th colspan="2">#mvps</th><th colspan="2">#iter</th></tr><tr><th>BCG / IB-BCG / IB-BCG-VA</th><th>BCG / IB-BCG / IB-BCG-VA</th><th>BCG / IB-BCG / IB-BCG-VA</th><th>BCG / IB-BCG / IB-BCG-VA</th></tr></thead><tbody><tr><td>cbuckle</td><td>15080 / 14681 / <b>14145</b></td><td></td><td><b>754</b> / 785 / 980</td><td></td></tr><tr><td>crankseg_1</td><td>6380 / 5714 / <b>5124</b></td><td></td><td><b>319</b> / 337 / 482</td><td></td></tr><tr><td>crankseg_2</td><td>7200 / 6696 / <b>6183</b></td><td></td><td><b>360</b> / 377 / 548</td><td></td></tr><tr><td>Kuu</td><td>3240 / <b>2987</b> / 3100</td><td></td><td><b>162</b> / 166 / 326</td><td></td></tr><tr><td>shallow_water2</td><td>560 / 560 / <b>420</b></td><td></td><td>28 / 28 / 28</td><td></td></tr><tr><td>ted_B_unscaled</td><td>560 / 543 / <b>411</b></td><td></td><td>28 / 28 / 28</td><td></td></tr></tbody></table>

## 3.5 Experiments with symmetric matrices

In this section, we consider testing three symmetric but not positive definite matrices described in Table 2, with linearly independent right-hand sides defined as $B = \text{randn}(n, p)$ with $p = 20$, $m_d = 5000 \times p$, $\epsilon = 10^{-8}$, and no preconditioner is applied. The corresponding convergence history and numerical results
are respectively reported in Figure 8 and Table 13. The observations are very similar to what we have seen for symmetric positive definite matrices in the previous sections, that is IB-BCR is often the best in minimizing the number of \#mvps at a possible extra cost of a few more block \#iter.

<figure><img src="image_12.png" /><figcaption>Figure 8: Comparison convergence history of the largest/smallest backward errors $\eta_{b(i)}$ at each \#mvps of the BCR variants by solving symmetry but not positive definite matrices listed in Table 2 with the right-hand sides $B = \text{randn}(n, p)$, $p = 20$, $\text{maxMvps} = 5000 \times p$ and $\varepsilon = 10^{-8}$.</figcaption></figure>

<table><caption>Table 13: Numerical results of the BCR variants in terms of \#mvps and \#iter for matrices listed in Table 2 with $B = \text{randn}(n, p)$, $p = 20$ and $\text{maxMvps} = 5000 \times p$ for Section 3.5.</caption><thead><tr><th>Matrix</th><th>Method</th><th>#mvps</th><th>#iter</th></tr></thead><tbody><tr><td rowspan="3">benzene</td><td>BCR</td><td>1380</td><td>69</td></tr><tr><td>IB-BCR</td><td><b>1231</b></td><td>70</td></tr><tr><td>IC-BCR</td><td>1325</td><td>69</td></tr><tr><td rowspan="3">rail_5177</td><td>BCR</td><td>2040</td><td>102</td></tr><tr><td>IB-BCR</td><td><b>1815</b></td><td>107</td></tr><tr><td>IC-BCR</td><td>1988</td><td>102</td></tr><tr><td rowspan="3">saylr4</td><td>BCR</td><td>380</td><td>19</td></tr><tr><td>IB-BCR</td><td><b>325</b></td><td>19</td></tr><tr><td>IC-BCR</td><td>371</td><td>19</td></tr></tbody></table>
# 4 Concluding remarks

In this work, we propose new variants of the block conjugate residual (BCR) method, incorporating mechanisms to address the partial convergence issues. These variants are denoted IB-BCR, and a cost-free alternative, IC-BCR. The IB-BCR and IC-BCR methods ensure the robustness by incorporating mechanisms to detect and handle partial convergence effectively. These enhancements improve the computational efficiency and reliability of the BCR variants in handling symmetric linear systems with multiple right-hand sides. We also explore the possibility of applying partial convergence detection, originally developed in the minimum residual norm context, to the block conjugate gradient (BCG) algorithm, in order to develop a heuristic solver called IB-BCG for symmetric positive definite linear systems. This numerical mechanism is considered heuristic because the resulting policy for expanding the search space still relies on residual norm minimization, while ideally, it should be based on minimizing the A-norm error. However, A-norm error directions are not a natural by-product of BCG, unlike the residuals in the case of BCR. The general observed trend is that the IB-BCR or IB-BCG variants perform the best in terms of the number of matrix-vector products, while the standard BCR or BCG methods excel in terms of block iterations. The best variant in terms of time depends on the trade-off between the final matrix-vector products, block iterations, and the additional cost associated with implementing the breakdown and partial convergence detecting mechanisms.

For long-term recurrence algorithms based on the Arnoldi basis, as discussed in [5], a sufficiently accurate orthonormal basis can be generated, ensuring that the true residual norm closely aligns with the least squares residual norm. This alignment enables monitoring of the true convergence behavior using the iterated residual norm. In contrast, for short-term recurrence algorithms such as BCR and BCG, residual gaps can occur. These gaps sometimes hinder the convergence of the true backward error, based on the true residual norm, even when the iterated residual norm indicates convergence. Thus, the partial convergence detecting mechanism applied to the iterated residual may become less reliable in short-term recurrence cases, especially when residual gaps occur. A potential remedy to address this issue is extending the rounding error analysis from [20] to the block case. Such an extension could provide a framework for estimating the block residual gap. Furthermore, it might facilitate the development of corresponding replacement techniques, as explored in [3] and [2], to improve the robustness of these short-term recurrence algorithms. Those topics might be the core of future research.

# References

[1] E. Agullo, L. Giraud, and Y.-F. Jing. "Block GMRES Method with Inexact Breakdowns and De- flated Restarting". In: *SIAM J. Matrix Anal. Appl.* 35.4 (2014), pp. 1625–1651. DOI: 10.1137/ 140961912.

[2] E. Carson and J. Demmel. "A Residual Replacement Strategy for Improving the Maximum Attainable Accuracy of s-Step Krylov Subspace Methods". In: *SIAM J. Matrix Anal. Appl.* 35.1 (2014), pp. 22–43. DOI: 10.1137/120893057.

[3] S. Cools, E. F. Yetkin, E. Agullo, L. Giraud, and W. Vanroose. "Analyzing the Effect of Local Round- ing Error Propagation on the Maximal Attainable Accuracy of the Pipelined Conjugate Gradient Method". In: *SIAM J. Matrix Anal. Appl.* 39.1 (2018), pp. 426–450. DOI: 10.1137/17M1117872.

[4] T. A. Davis and Y. Hu. "The University of Florida sparse matrix collection". In: *ACM Trans. Math. Softw.* 38 (2011), p. 1.

[5] L. Giraud, Y.-F. Jing, and Y.-F. Xiang. "A Block Minimum Residual Norm Subspace Solver with Partial Convergence Management for Sequences of Linear Systems". In: *SIAM J. Matrix Anal. Appl.* 43.2 (2022), pp. 710–739. DOI: 10.1137/21M1401127.

[6] M. Hestenes. "The conjugate gradient method for solving linear systems". In: *Numerical Analysis. Ed. by J. H. Curtiss*. Providence Rhode Island: American Mathematical Society, 1956, pp. 83–102.

[7] M. Hestenes and E. Stiefel. "Methods of Conjugate Gradients for Solving Linear Systems". In: *Journal of Research of the National Bureau of Standards* 49.6 (1952), pp. 409–436.

[8] H. Ji and Y.-H. Li. "A breakdown-free block conjugate gradient method". In: *BIT* 57 (2017), pp. 379–403. DOI: 10.1007/s10543-016-0631-z.
[9] J. Langou. *Iterative Methods for Solving Linear Systems with Multiple Right-Hand Sides*. Ph.D. dissertation TH/PA/03/24. CERFACS, Toulouse, France, 2003.

[10] M. Levonyak, C. Pacher, and W. N. Gansterer. "Scalable Resilience Against Node Failures for Communication-Hiding Preconditioned Conjugate Gradient and Conjugate Residual Methods". In: *Proceedings of the 2020 SIAM Conference on Parallel Processing for Scientific Computing*. 2020, pp. 81–92. DOI: 10.1137/1.9781611976137.

[11] D. Luenberger. "The Conjugate Residual Method for Constrained Minimization Problems". In: *SIAM J. Numer. Anal.* 7.3 (1970), pp. 390–398. DOI: 10.1137/070703.

[12] D. P. O'Leary. "The block conjugate gradient algorithm and related methods". In: *Linear Algebra Appl.* 29 (1980), pp. 293–322. DOI: 10.1016/0024-3795(80)90247-5.

[13] L. F. Pavarino. "Preconditioned conjugate residual methods for mixed spectral discretizations of elasticity and Stokes problems". In: *Comput. Methods Appl. Mech. Engrg.* 146.1-2 (1997), pp. 19–30. DOI: 10.1016/S0045-7825(96)01224-8.

[14] M. Robbé and M. Sadkane. "Exact and inexact breakdowns in the block GMRES method". In: *Linear Algebra Appl.* 419.1 (2006), pp. 265–285. DOI: 10.1016/j.laa.2006.04.018.

[15] Y. Saad. *Iterative Methods for Sparse Linear Systems*, 2nd ed. Philadelphia: SIAM, 2003.

[16] Y. Saad and M. H. Schultz. "GMRES: A Generalized Minimal Residual Algorithm for Solving Non-symmetric Linear Systems". In: *SIAM J. Sci. Stat. Comput.* 7.3 (1986), pp. 856–869. DOI: 10.1137/0907058.

[17] T. Sogabe. *Extensions of the Conjugate Residual Method*. Ph.D. dissertation. http://www.ist.aichi-pu.ac.jp/person/sogabe/thesis.pdf. University of Tokyo, Japan, 2003.

[18] T. Sogabe, M. Sugihara, and S. L. Zhang. "An extension of the conjugate residual method to non-symmetric linear systems". In: *Journal of Computational and Applied Mathematics* 226.1 (2009), pp. 103–113. DOI: 10.1016/j.cam.2008.05.018.

[19] E. de Sturler. "Nested Krylov methods based on GCR". In: *J. Comput. Appl. Math.* 67.1 (1996), pp. 15–41. DOI: 10.1016/0377-0427(94)00123-5.

[20] H. A. van der Vorst and Q. Ye. "Residual Replacement Strategies for Krylov Subspace Iterative Methods for the Convergence of True Residuals". In: *SIAM J. Sci. Comput.* 22.3 (2000), pp. 835–852. DOI: 10.1137/S1064827599353865.
<figure><img src="image_13.png" /></figure>

RESEARCH CENTRE
BORDEAUX - SUD-OUEST

200 avenue de la Vieille Tour
33405 Talence Cedex

Publisher
Inria
Domaine de Voluceau - Rocquencourt
BP 105 - 78153 Le Chesnay Cedex
inria.fr

ISSN 0249-6399