A flexible short recurrence Krylov subspace method for matrices arising in the time integration of port Hamiltonian systems and ODEs/DAEs with a dissipative Hamiltonian

For several classes of mathematical models that yield linear systems, the splitting of the matrix into its Hermitian and skew Hermitian parts is naturally related to properties of the underlying model. This is particularly so for discretizations of dissipative Hamiltonian ODEs, DAEs and port Hamiltonian systems where, in addition, the Hermitian part is positive definite or semi-definite. It is then possible to develop short recurrence optimal Krylov subspace methods in which the Hermitian part is used as a preconditioner. In this paper we develop new, right preconditioned variants of this approach which as their crucial new feature allow the systems with the Hermitian part to be solved only approximately in each iteration while keeping the short recurrences. This new class of methods is particularly efficient as it allows, for example, to use few steps of a multigrid solver or a (preconditioned) CG method for the Hermitian part in each iteration. We illustrate this with several numerical experiments for large scale systems.

Here, * denotes the adjoint with respect to the standard inner product, i.e., for any matrix B the matrix B * is obtained by transposing B and taking complex conjugates for each entry.
In several important applications, and in particular in linear systems arising from discretizations of dissipative Hamiltonian ODEs and certain port Hamiltonian systems, the hermitian part H is positive definite and "dominates" S in the sense that H −1 S or S is small.Then, H is an efficient preconditioner for A, and as we will discuss in Section 2, there exist optimal short recurrence Krylov subspace methods using this preconditioner.
Each iterative step of these methods requires one exact solve with the matrix H.The purpose of this paper is to go one step further and investigate practically important short recurrence methods where the solves with H are performed only approximately.For example, we might want to use an incomplete Cholesky factorization in cases where an exact factorization is too expensive, or we might want to perform just some steps of (plain or preconditioned) CG or a multigrid solver for H.We believe that this is an aspect of primordial importance: As long as systems are small enough such that direct factorizations with A are feasible, there is no real need to consider an iterative method which, after all, requires again a direct factorization of a matrix, H, of the same size.Thus we expect that variants which allow for inexact solves of systems with H to have practical impact.
The paper is organized as follows: In section 2 we review the existing methods with exact solves for H and relate them to a known method for systems where the matrix is Hermitian plus a purely imaginary multiple of the identity.In section 3 we then develop right preconditioned variants of the methods existing in the literature and use those to develop flexible variants.We describe these algorithmically and analyze some of its most relevant properties.In section 5 we present results of numerical experiments for four examples demonstrating the efficiency gains of the flexible methods over those relying on exact solves.
2. Review of methods with exact solves for H. From now on we assume that the Hermitian part H of A is positive definite.The method proposed by Concus and Golub [3] and, independently, by Widlund [27] as well as the method of Rapoport [18] are, mathematically, equivalent to the full orthogonalization method FOM (see [20], e.g.) and the minimal residual method MINRES [17], respectively, for the preconditioned system with the standard inner product replaced by the H-inner product x, y H := y * Hx.
In this section, we present the essentials of both these methods.The starting point is to observe is that the matrix H −1 S in the preconditioned matrix H −1 A = I +H −1 S is H-anti-selfadjoint according to the following definition and subsequent proposition (whose proof is straightforward).
Because I + H −1 S is a shifted H-anti-selfadjoint matrix, the Arnoldi process for the H-inner product produces a three term recurrence.Introducing the notation K m (A, b) for the Krylov subspace K m (A, b) = span(b, Ab, . . ., A m−1 b) of order m generated by the matrix A and the vector b, this short recurrence property can be formulated as where the columns and T m+1,m ∈ R (m+1)×m is tridiagonal.Denoting T m the matrix obtained from T m+1,m by removing the last row, we also have that T m,m − I is skew-symmetric.The iterates of the method of Concus and Golub / Widlund are now variationally characterized by a Galerkin condition which we denote GAL to stress the left preconditioning; the iterates of the method of Rapoport are characterized my a minimal residual condition denoted MR: Note that in the formula for the m-th iterate of the MR variant we find the standard 2-norm in the least squares problem for ζ m .
The short recurrences for computing the iterates of GAL and MR with the H-inner product arise using an LU or a QR factorization of T m or T m+1,m which is updated from one iteration to the next.All eigenvalues of T m are of the form 1 + iµ with µ ∈ R so that T m as well as all its principal minors cannot become singular, which guarantees that a non-pivoted LU-factorization always exists.For the MR variant one typically uses a QR factorization of T m+1,m with an implicit representation of Q as a sequence of Givens rotations and R having three non-zero diagonals (the diagonal and two superdiagonals).Overall this process is analogous to the implementations proposed for SYMMLQ and MINRES in [17], to the implementation of the GAL, MR and ME methods in [5], to the implementation of GAL in [3], [27] and of MR in [18] and to the implementation of QMR as described in [6].We therefore spare ourselves from reproducing technical details, the principles of which can also be found, e.g., in the textbooks [12], [20] and in [23].In here iS =: B is now Hermitian, and preconditioning with H −1 yields a shifted H-selfadjoint matrix In his 1990 paper, Freund [5] considers CG type methods for matrices of the form There, three short recurrence Krylov subspace methods, all within the framework of the standard 2 inner product are developed: GAL, MR and ME.GAL is variationally characterized by a Galerkin condition, MR by a minimal residual condition, and ME ("minimal error") minimizes the error in 2-norm over the subspace x 0 + (iσI + W ) * K m (iσI + W, r 0 ).The ME approach dates back to [7] for the case σ = 0 and is less commonly used today.In our numerical experiments it behaved quite similarly as the other methods, but with slightly larger errors and residuals than GAL, so we do not consider this variant any further in this paper.The methods in [5] can be easily adapted to work with the H-inner product rather than with the standard inner product.The requirement is that W is then an H-selfadjoint matrix, which is exactly the situation that arises when we precondition (3.1) with H −1 , yielding W = H −1 B (and σ = 1).For any vector r, the Krylov subspaces K(iσI + H −1 B, ir) and K(σI + H −1 S, r) are identical.This means that for a given initial guess x 0 , the defining Galerkin conditions produce identical iterates, independently of whether we consider (σI + H −1 S)x = b with initial residual r = b − (σI + H −1 S)x 0 or (iσI + B)x = ib with initial residual ib − (iσI + B)x 0 = ir.The same holds for the minimal residual based methods.We emphasize this relation in the following proposition for the case σ = 1.Proposition 3.1.Freund's GAL method with the H-inner product for (iI + H −1 B)x = ib is identical to the Concus and Golub / Widlund method for (I + H −1 S)x = b in the sense that it produces the same iterates if we have the same initial guess.Freund's MR with the H-inner product for iI + H −1 B is, in the same sense, identical to Rapoport's method for I + H −1 S.
There such equivalent to Freund's ME seems to have been published so far.As with left preconditioning, there is a direct and simple correspondence with the "Hermitian" formulation on which we will focus in what follows, i.e.
Herein, W = BH −1 is H −1 -selfadjoint, and we again obtain short recurrence methods GAL, MR, ME by using the H −1 inner product in the methods and algorithms of [5].
Proposition 3.2.There exist short recurrence GAL, MR and ME methods for the right preconditioned system iσI + BH −1 .These give rise to short recurrence methods for iσH + B when using the H −1 inner product.
Although it might not be obvious at first sight, the implementation of right preconditioning can be done by investing one multiplication with B, one with H and just one solve with H per iteration.The idea is to carry an additional vector ṽ which stores H −1 v for the current Lanczos vector v. Details are given in Algorithm 3.1 in which, when solving the linear system (3.1),we take the initial vector w as the initial residual ib − (iσH + B)x 0 .We note that an implementation with similar cost is possible for left preconditioning.Also note that in the algorithm we have which is why we can use Algorithm 3.1 m steps of Lanczos for iσI + BH −1 with the H −1 inner product.
•, • denotes the standard inner product One can work with σI + SH −1 rather than with iσI + BH −1 in Algorithm 3.1, adapting the matrix vector multiplication in line 4 and changing the last summand This observation represents the general transition rule from formulations for the systems with iσI + BH −1 , iσI + H −1 B or iσH + B to systems with σI + SH −1 , σI + H −1 S or σH + S. Our presentation will continue to be based on the systems using B rather than S.
Right preconditioning is the key to develop flexible variants of the GAL, MR and ME methods.The term "flexible" includes "inexact" methods in which w = H −1 w is replaced by w = H −1 w with a fixed matrix H.For example, H may result from an incomplete Cholesky factorization."Flexible" also includes the more general possibility to approximate H −1 w via a non-stationary iteration such as (possibly preconditioned) CG which, formally, means that H changes from one approximate solve to the next.

Flexible methods.
Using exact multiplications with H −1 , the Lanczos algorithm given in Algorithm 3.1 for the right preconditioned matrix iσI + BH −1 produces the relation , where now V m has H −1 -orthonormal columns.The transformation of an iterate for the right preconditioned system back to the iterate of the original system is most easily done via a short recurrence update using the preconditioned Lanczos vectors v m = H −1 v m which we compute anyways.The relation (3.2) can then be stated as We can now divise a flexible method in a manner similar to what has been done for FGMRES, the flexible GMRES method [26] and, with regard to short recurrences, in flexible QMR [24], e.g.Assume that z k is only an approximation to H −1 v k and proceed as follows: This means that the inner products w, v k H −1 and w, v k−1 H −1 are obtained as standard inner products with w where w is an approximation to H −1 w.
3. Approximately H −1 -normalize the vector resulting from the approximate H −1 orthogonalization.Algorithm 3.2 gives the details for one iteration of the resulting right preconditioned flexible Lanczos process.Algorithm 3.2 m steps of right preconditioned flexible Lanczos for iσH + B with approximate solves for H. •, • denotes the standard inner product w ≈ H −1 w 9: , and we avoided to compute w, v m−1 explicitly.In the flexible context the above equality does not hold any more, since w, v m−1 is computed with the approximation z m−1 for H −1 v m−1 , whereas β m is computed with the approximation w for H −1 w.
With the flexible Lanczos process we have the recurrence relation were T m+1,m is again tridiagonal, When used to solve the linear system (3.1), the initial vector w in Algorithm 3.2 is the residual ib − (iσH + B)x 0 , i.e., the first Lanczos vector v 1 is We emphasize that V m is now only approximately H −1 -orthonormal and Z m is only approximately equal to H −1 V m .We define the flexible GAL and MR iterates in analogy to (2.3) and (2.4) as FMR: where now V m , T m+1,m and T m obey relation (3.4).Computationally, we obtain short recurrences for the iterates x m by updating factorizations of T m,m and T m+1,m just as in the non-flexible case; see the discussion at the end of section 2. Commented Matlab implementations of FMR and FGAL are available from the authors on request.
Note that the FGAL iterates do not fulfill a strict H −1 -orthogonality relation for the residuals, nor do the FMR iterates minimize the H −1 -norm of the residual exactly.But, of course, if we solve for H exactly in the flexible Lanczos process, FGAL and FMR reduce to right preconditioned counterparts of FOM and MR from (2.3) and (2.4), respectively.
A flexibly preconditioned CG method (FCG) for Hermitian positive definite matrices using a preconditioner which may change from one iteration to the next has been proposed in [8] and further analyzed in [16].The FGAL method can be used in this situation, too, with σ = 0. FCG and FGAL are similar in spirit but not identical: FCG constructs search directions by applying the variable preconditioner to the current residual and then A-orthogonalizing against a fixed number s of the previous search directions.In contrast, FGAL orthogonalizes the next preconditioned Lanczos vector in the 2 -inner product against the two previous ones.This is different, even when s = 2.

Properties.
In this section we first state two known convergence bounds for the left preconditioned methods where systems with H are solved exactly.We then show that a result from [5] can be used to improve one of these bounds, a fact that seems to have gone unnoticed so far.We continue the section considering the flexible variants for which we establish two somewhat more basic but algorithmically important properties.Theorem 4.1.Let λ ≥ 0 be the smallest value such that the interval [−iλ, iλ] contains all the (purely imaginary) eigenvalues of the H-anti-selfadjoint matrix H −1 S and let x * denote the solution of the system Ax = b.Then (i) For the GAL iterates x m from (2.1) we have Part (i) of the theorem goes back to [4,11,25], and part (ii) to [25]; see [9] for a detailed discussion.What seems to have gone unnoticed so far is that a better bound for the MR iterates results from [5,Theorem 4], adapted to H-selfadjoint matrices.We first formulate the result for an arbitrary value for the shift σ before comparing to Theorem 4.1(ii).Note that the original theorem in [5] considers residuals for the matrix iσI + W with W Hermitian and the 2 -norm.For the residuals belonging to matrices iσI + W with W H-selfadjoint it holds in exactly the same manner, but now with the H-inner product and associated norm.When translating back to our original matrix H + S, we have σ = 1, and H-norms for the residual w.r.t.iI + H −1 B turn into H −1 -norms for the residuals w.r.t.H + S.This gives us the following theorem.Theorem 4.2.Let spec(H −1 S) ⊆ i[α, β] with α < β and let R be the unique solution of (4.1) Then the residuals of the MR iterates Compared to Theorem 4.1 this theorem allows for a more narrow interval to contain the spectrum of H −1 S as it does not need to be symmetric with respect to the real axis.Note, however, that for real matrices spec(H −1 S) is always symmetric with respect to the real axis, so that optimal bounds on the spectrum are of the form i[−λ, λ] as in Theorem 4.1.But even with symmetric bounds for the spectrum, the bound of Theorem 4.2 is always better than that of Theorem 4.1 as we will explain now.In order to not unnecessarily burden our discussion, we will omit details of somewhat longish but elementary algebraic manipulations in the paragraph to follow.
We first note that the expression ( β 2 + 1 + √ α 2 + 1)/(β − α) increases monotonically in β and decreases monotonically in α.Thus, for λ ≥ max{β, −α}, we have We now turn to the flexible methods.Conceivably, a detailed convergence analysis, depending on the "degree of inexactness" of the preconditioner could be based on prior work for flexible CG [16], flexible GMRES [19,26] and, in particular, flexible QMR [24].We anticipate this to be quite involved and also quite technical, so that we do not address this in the context of the present paper.We rather state and prove two basic properties which are essential for the algorithmic aspects.

The solution
The first property is that in typical situations we cannot encounter unlucky breakdowns in the flexible Lanczos process, Algorithm 3.2.An unlucky breakdown at step m means that we get β m = w, w = 0 although w = 0.If there is no unlucky breakdown, the FMR iterates all exist up to the step where we obtain the solution, as in that case the matrix T m+1,m of (3.5) has full rank for all relevant m.Though, its square part T m,m need not be non-singular, which means that the FGAL iterate does not necessarily exist.However, as is done in the SYMMLQ-algorithm (see [12,17]), the QR-factorization of T m,m can still be updated from that of T m−1,m−1 , which allows to efficiently update the next FGAL iterate from the previous existing one using the QR-factorization.Proposition 4.3.Assume that the inexact solves w ≈ H −1 w in Algorithm 3.2 are done in one of the following ways (i) by performing k steps of the CG method with initial guess 0, (ii) by performing k steps of the preconditioned CG method with initial guess 0 and a Hermitian positive definite preconditioner K, (iii) as w = Gw with G ∈ C n×n Hermitian and positive definite.Then, if w = 0, we have β = w, w > 0.
Proof.In case (i), the initial residual for the CG method is w, and with the columns of W ∈ C n×k being the Lanczos vectors which form an orthonormal basis of the Krylov subspace K k (H, w), we know that w = V k (V * k HV k ) −1 V * k w; see e.g.[12] or [20].Thus, , and this quantity is positive as V * k HV k is Hermitian positive definite and V * k w = w 2 e 1 , e 1 the first canonical unit vector in C k , is non-zero.
Case (ii) follows in a similar manner as (i), replacing V k by the matrix of Korthonormal Lanczos vectors for K k (K −1 H, K −1 w).
In case (iii), as G is Hermitian positive definite and w = 0, we immediately have β = w, w = w, Gw > 0.
An important example for case (iii) above is when G = L −1 L − * with L a sparse lower triangular matrix arising from an incomplete Cholesky factorization of H, H = L * L−R.Such G is also a typical preconditioner K −1 for case (ii).Another important example for (iii) are Galerkin based multigrid methods for Hermitian positive definite systems with symmetric pre-and post-smoothing and exact coarsest level solves.
The second property we want to stress deals with b − Ax m H −1 as a measure of the error of iterate Thus, if H dominates S, i.e., if H −1 B is small, the matrix AH −1 A is close to σ 2 H (and we have equality if B = 0).Consequently, when H dominates B, up to the trivial factor σ 2 , the norm b − Ax m H −1 is a good approximation to the H-norm of the error.There are multiple reasons why in the case B = 0 the H-norm of the error-and not its 2 -norm nor the 2 -norm of the residual, e.g.-should be considered the most adequate measure for the error, in particular when the linear systems arise from a discretization of an underlying continuous, infinite-dimensional equation; see [1,13], e.g.Even if B = 0, the H-norm can still be regarded as the canonical measure for the error.By the following proposition, an approximate bound for this norm is available at almost no cost for the FMR iterates.Then and, if the flexible Lanczos vectors z m approximate the "exactly inverted" vector H −1 v m to a relative accuracy of ε in the H-norm, i.e.
Proof.Denote by V m ∈ C n×m the matrix which contains the flexible Lanczos vectors v k as its k-th column.We have which is (4.2).To obtain (4.4) we now show that the 2 -norm of each column of H −1/2 V m+1 is bounded by 1/(1 − ε).The inequality (4.4) then follows as the 2norm of a matrix never exceeds its Frobenius norm. With from which we get, via the Cauchy-Schwarz inequality, With assumption (4.3) we therefore have There are two important observations regarding Proposition 4.4.First, practically, we cannot monitor H-norms of the error, so inequality (4.3) cannot be checked directly.However, a relative reduction of the ( 2 -norm of the) residual, which can be measured, typically results in a similar reduction of the error.Moreover, in the case that we use the CG method to approximate H −1 v k , it is precisely the H-norm of the error which is made as small as possible by the CG iterates.We refer to [21] for a further discussion of this and related aspects.
Secondly, the quantity m is available in the algorithm at no additional cost.The situation is similar as in GMRES [20] or QMR [6]: To solve the least sqares problem (3.8), we update the QR-factorization T m+1,m = Q m+1,m+1 R m+1,m in each step.The value m is the last, i.e. the (m + 1)st, entry of Q * m+1,m r 0 H e 1 , and this value is updated to m+1 using the Givens rotation which updates Q m+1,m to Q m+2,m+1 .
We further note that the proof of Proposition 4.4 also shows that in the case where we solve for H exactly, we find m = b − Ax m H −1 as V m H −1 = 1.

Numerical results.
In this section we report results for several examples which have been considered in the literature.The general setup for all experiments is as follows: We report residual norm plots for both, the FMR and the FGAL method.In each iteration, the systems with the hermitian positive definite part H are solved approximately using the CG method, asking for the initial 2 -norm of the residual to be reduced by a factor of ε cg .We typically vary ε cg from 10 −1 up to 10 −12 .In light of Proposition 4.4 and the discussion preceding it, the FGAL and FMR iteration is stopped once the H −1 -norm of the initial residual is reduced by ε f = 10 −12 .The case ε cg = 10 −12 in the approximate solves may thus actually be considered as a equivalent to exact solves in our context.
The residual norm plots report (relative) approximate H −1 -norms of the residual, i.e. the quantity √ r * r where r is the current residual and r is the result of CG when solving H r = r with initial guess 0 and a required reduction of the residual 2 -norm of ε cg .This number can be somewhat off the exact H −1 -norm; the advantage is that it is computed with an effort comparable with what we need to update an iterate.When the iterations stopped, we double-checked that the targeted decrease of the H −1 -norm was indeed achieved by recomputing √ r * r for the final residual with now ε cg = 10 −12 in the CG method which computes r.Finally, as an illustration, our plots also give the bound on the H −1 norm of the residuals of the FMR iterates according to Proposition 4.4.Let us stress that, as opposed to the approximate H −1 -norms, these bounds can be computed at virtually no extra cost, so that in a production setting we could base a stopping criterion for the FMR iterations merely on these bounds.We start reconsidering the example in [27].on Ω = [0, 1] × [0, 1] with Dirichlet boundary conditions.We discretize on a uniform square mesh with 127×127 interior grid points, using central finite differences in both, ∆u and u x , and choose a = 10 4 .The right hand side b was chosen to have random entries.
Figure 1 shows convergence plots for four different choices of ε cg .We see that for all ε cg the residual norms of the FGAL and the FMR iterates differ only marginally, with the FMR residuals, as is to be expected due to their (approximate) minimization property, being (slightly) smaller and decreasing in a smoother manner.We also see that the bound from Proposition 4.4 captures the actual convergence behavior quite accurately.With ε cg = 10 −1 , we already need only about twice as many iterations as with "exact" solves, ε cg = 10 −12 , while an approximate solve with ε cg = 10 −1 requires about 50 CG iterations on average as opposed to about 470 CG iterations on average for ε cg = 10 −12 .
Example 2. This is the matrix treated in section 6.2 of [9].It is based on a discretized Stokes problem obtained by using the Q 1 -Q 1 finite element discretization of the unsteady channel domain problem in IFISS [22].This example was treated for two different grid parameters in [9].We only consider the larger of the two, i.e grid parameter 10, resulting in a linear system of size n = 3, 151, 875, and we took a step size of τ /2 = 10 −3 in the midpoint Euler rule.
The results are reported in Figure 2. To solve the systems with H we use preconditioned CG with the ILU(0) preconditioner.In the solves for H, every about 12 preconditioned CG iterations reduce the norm of the residual by another factor of 10 −1 .We observe that the dependence on the CG tolerance ε cg is similar as in the previous example, with ε cg = 10 −1 already requiring just about half times more the number of iterations than with "exact" solves (ε cg = 10 −12 ).FGAL and FMR behave similarly with FGAL showing a slightly more oscillatory convergence behavior than before.This example also shows stagnation once the residual norm is reduced by a factor of approximately 10 −10 for FGAL and 10 −11 for FMR.Note that here as in all the other plots we always report the norms of "computed" residuals, i.e. residuals which were explicitly computed from the current iterate in each iteration.We attribute the observed stagnation to the relatively high condition number κ of the matrix which makes us expect that the reported residual norms, which rely on explicitly computed residuals, will not go beyond roughly κ • ε mach , with ε mach ≈ 10 −16 being the machine roundoff unit.The reported residual bounds are based on updated quantities, which is why they continue to decrease beyond the point where the computed residuals stagnate.
Interestingly, the numerical computations in [9] used only an incomplete Cholesky factorization to solve the systems with H-the incomplete Cholesky drop tolerance in Matlab's ichol was set to the relatively small value of 10 −9 .Thus, the experiments reported in [9] perform approximate solves only, while the algorithms used there are not flexible.This might be the reason why the experiments in [9], as the authors say, could not use larger values than 10 −3 for τ /2.
The remaining two examples are linear systems that arise when using the implicit midpoint Euler rule when integrating a DOE resulting from a port Hamiltonian where E is hermitian and positive definite, J is anti-selfadjoint and R is hermitian and positive semidefinite.With a step size of τ > 0, the linear systems to solve at each time step are then of the form and we have H = E + τ 2 R, positive definite, and S = − τ 2 J.Note that implicit integration methods are mandatory in this situation if one wants to preserve the dissipation of energy; see [14].
Example 3. We take the homonomically constrained spring mass system from [15], which is part of the port Hamiltonian system benchmark collection1 .We refer to [10] and [9] on how the matrices E, R and J arise from the descriptor system formulation of the system.We took the system with n = 2 • 10 6 as it was also considered in [9] and τ /2 = 0.1.The numerical results are given in Figure 3.For all choices of ε cg , the norms of the residuals of FGAL and FMR are so close that they are indistinguishable in the plots, which is why we report only the results for FMR.The bounds from Proposition 4.4 again match the convergence behavior quite exactly.Similar to the previous two examples, we see that with ε cg = 10 −1 we need less than twice as many iterations than with "exact" solves (ε cg = 10 −12 ), and this number of iterations is this time exceeded by just 1 for ε cg = 10 −2 .We need between 4 and 6 (unpreconditioned) CG iterations with H to reduce the initial residual by the factor ε cg = 10 −1 in each step of MR, and between 43 and 49 to reduce it by the factor ε cg = 10 −12 .
In the three plots for ε cg = 10 −1 , 10 −2 and 10 −12 we in addition plot the H-norms of the residuals in a variant of FMR which differs from FMR only by the fact that in the underlying preconditioned Lanczos algorithm, Algorithm 3.2, we do not compute γ k but rather put γ k = β k−1 .This method is termed "non-flexible MR".The plots show that for the larger values of ε cg this modification results in a significantly to dramatically delayed convergence.
In a similar spirit, the bottom right plot of Figure 3 show the H −1 -norm of the residual for the non-flexible MR method of Rapoport (left preconditioning), where we solve for H with CG at different accuracies.Not really surprisingly, and convincingly, the plot illustrates that this "inexact" version of MR stagnates sooner the less accurately we solve for H.A similar behavior can be observed for non-flexible inexact GAL, for which we do not report results here.In contrast, the flexible methods reach high accuracies for the overall system even when ε cg is large.
Example 4. We take the coupled electro-thermal DAE system from [2], from which we treat the thermal part after applying the decoupling described in [2].This again yields an ODE of the form (5.1).Similar to the previous example, we consider the linear systems arising when using the implicit midpoint Euler method for time integration.The system has dimension n = 7.2 × 10 5 , and we take τ /2 = 10 −2 .Since the positive definite diagonal matrix E is highly ill-conditioned, we scaled the system from left and right with the square root of the inverse of E.
The results are reported in Figure 4. To solve the systems with H we use preconditioned CG with a threshold dropping incomplete Cholesky preconditioner with drop tolerance 10 −4 .This results in a very efficient preconditioner, since on average we need just slightly less than two iterations to decrease the residual by one order of magnitude.Similarly to the previous example, FMR and FGAL perform very similarly and would be indistinguishable in the plots, which is why we only report the results for FMR.As before, we observe that FMR (and FGAL) require about twice the iterations with ε cg = 10 −1 than with "exact" solves (ε cg = 10 −12 ).Both flexible methods again achieve small final residuals even when large values of ε cg are used.
For non-flexible MR, the iterations stagnate quite early, except for ε cg = 10 −12 , in which case the non-flexible method cannot be distinguished from FMR.

1 .
Introduction.Consider the linear system Ax = b, A ∈ C n×n , and the splitting A = H + S with H = 1 2 (A + A * ) the Hermitian and S = 1 2 (A − A * ) the skew Hermitian part of A, H * = H, S * = −S.

3 .
Right preconditioned and flexible methods.If we multiply the original equation Ax = (H + S)x = b with the imaginary unit i we obtain (3.1) (iH + iS)x = ib.