➨  ᅇඛ㐍ࢫ࣮ࣃ࣮ࢥࣥࣆ࣮ࣗࢸ࢕ࣥࢢ⎔ቃ◊✲఍㸦$6( ◊✲఍㸧Ⓨ⾲㈨ᩱ
ᮾி኱Ꮫ᝟ሗᇶ┙ࢭࣥࢱ࣮ ≉௵෸ᩍᤵ ∦᱒Ꮥὒ  ᖺ  ᭶  ᪥㸦㔠㸧 ᫬  ศ㹼 ᫬  ศࠊᮾி኱Ꮫ᝟ሗᇶ┙ࢭࣥࢱ࣮㸲㝵㐲㝸ㅮ⩏ᐊ࡟࡚ࠊ➨㸴ᅇඛ㐍ࢫ࣮ࣃ࣮ࢥࣥࣆ࣮ࣗࢸ࢕ࣥࢢ⎔ቃ◊✲఍㸦$6( ◊✲఍㸧ࡀ㛤ദࡉࢀࡲࡋࡓࠋᮏྕ࡛ࡣࠊࢻ࢖ࢶࡢ+DPEXUJ+DUEXUJᕤ⛉኱Ꮫࡢ -HQV3HWHU=HPNH༤ኈ࡟ࡼࡿᣍ
ᚅㅮ₇ࠕ,'5 V DQG,'5 V (LJLQ3DUDOOHO&RPSXWLQJࠖࡢෆᐜ࡟㛵ࡋࠊㅮ₇ᚋ࡟ =HPNH
༤ኈࡀᡃࠎࡢࡓࡵ࡟≉ู࡟ࡈ⏝ព࠸ࡓࡔ࠸ࡓㄽᩥࢆ㈨ᩱ࡜ࡋ࡚ᥖ㍕ࡉࡏ࡚࠸ࡓࡔࡁࡲࡍࠋ
 $6( ◊✲఍࡛ࡣࠊඛ㐍ࢫ࣮ࣃ࣮ࢥࣥࣆ࣮ࣗࢸ࢕ࣥࢢ⎔ቃࢆᨭ᥼ࡍࡿ◊✲ヰ㢟ࡢᥦ౪ࢆ┠ⓗ࡟ᖺᩘᅇࡢάືࢆィ⏬ࡋ࡚࠾ࡾࡲࡍࠋࡇࢀ࠿ࡽࡶࡈᨭ᥼ࡢ࡯࡝ࢆࠊࡼࢁࡋࡃ࠾㢪࠸⏦ࡋୖ
ࡆࡲࡍࠋ

௨ୖ

IDR(s) AND IDR(s)Eig IN PARALLEL COMPUTING∗JENS-PETER M. ZEMKE†

Abstract. In this short summary of a talk held at the University of Tokyo on the 12th February2010 we consider a certain class of Krylov subspace methods: the IDR methods by Peter Sonneveld,
and their recent generalizations to incorporate more shadow vectors, the IDR(s) methods. We
sketch the main ideas behind IDR, IDR(s), and IDR(s)Eig, and indicate possible approaches fortheir parallelization.
Key words. IDR; PIA; IDR(s); IDR(s)Eig; Krylov subspace methods; Eigenvalues; Linear
Systems; Parallelization
AMS subject classiﬁcations. 65F10 (primary), 65F15, 65Y05, 65Y20

1. IDR and IDR(s). IDR [16, 17, 24] and IDR(s) [15, 18, 23, 19, 13, 14, 20, 11]
are Krylov subspace methods which have been invented by numerical analyst Peter
Sonneveld at the TU Delft in 1976 and 2006, respectively, the latter with co-author
Martin van Gijzen. Recently there has been an increased interest in these methods,
as they appear to be more stable than other transpose-free Krylov subspace methods
based on short recurrences in certain situations. An important predecessor of the
new family IDR(s) is the less well known method ML(k)BiCGStab by Man-Chung
Yeung and Tony Chan [26].
1.1. Krylov subspace methods. IDR and IDR(s) are Krylov subspace methods. The mth Krylov subspace Km is deﬁned for a given square matrix A ∈ Cn×n
and a starting vector q ∈ Cn as follows,
Km (A, q) := span {q, Aq, . . . , Am−1 q}.

(1.1)

As long as the Krylov subspace Km has full dimension dim(Km ) = m, there is anatural isomorphismv ∈ Km⇔v = ν(A)q

between vectors v in a Krylov subspace and polynomials ν ∈ Pm−1 . Obviously there
is a last index m  n with dim(Km ) = m, this index is known as the grade of thevector q with respect to the matrix A. The Krylov matrices


(1.2)
Km := q, Aq, A2 q, . . . , Am−1 qsatisfy the matrix recurrenceq, AKm = Km+1 .

(1.3)∗ Based in part on the forthcoming paper [5]. The author acknowledges the invitation by Prof.
Seiji Fujino of the Research Institute for Information Technology of Kyushu University to Japan in
February 2010 where part of this work was written.
† Institut für Numerische Simulation, Technische Universität Hamburg-Harburg, D-21073 Hamburg, Germany (zemke@tu-harburg.de).

1

The mth Krylov matrix spans a basis of the mth Krylov space Km if and only if m isless or equal to the grade of q. We assume here that this is always the case. Suppose
now that we choose upper triangular basis transformations Km =: Qm Rm ,q, AQm Rm = Qm+1 Rm+1 .
(1.4)
The Krylov matrix recurrence (1.3) can now be reformulated,
1q, AQm = Qm+1 Rm+1
o

oT
Rm

−1
.

(1.5)

These upper triangular basis transformations are natural in the context of Krylov
subspace methods, as in step k only the ﬁrst k columns are available. When we strip
oﬀ the ﬁrst column on both sides of Eqn. (1.5), we exhibit the intimate connection of
Krylov subspace methods to so-called Hessenberg decompositions.
Let the matrix Cm ∈ C(m+1)×m be deﬁned by



o


Cm


:= Rm+1

1
o

oT
Rm

−1
.

(1.6)

By the group structure of regular upper triangular matrices, this matrix is unreduced
extended Hessenberg.
In every simple Krylov subspace method, i.e., excluding block variants, we end
up with a so-called Hessenberg decomposition1
AQm = Qm+1 Cm =: Qm Cm + qm+1 cm+1,m eT
m,

(1.7)

where Cm is unreduced Hessenberg and measures the “ratio” of the basis transformations.
Krylov subspace methods can be classiﬁed into distinct approaches. We ﬁrst givethe matrix based classiﬁcation. In the following, we always assume that in the context
of the approximate solution of linear systems, the starting vector q1 has been chosenby r0 2 q1 = r0 . There are three well-known approaches based on such Hessenberg
decompositions, namely,
QOR: approximate x = A−1 r0 by xm := Qm C−1
m e1 r0 2 ,QMR: approximate x = A−1 r0 by xm := Qm C†m e1 r0 2 ,Ritz-Galërkin: approximate part of J = V−1 AV by Jm := S−1
m Cm Smand part of V by Vm := Qm Sm , where Cm Sm = Sm Jm .
To every method from one class corresponds a method of the other. This fact is usedin the forthcoming publication [5] to compute eigenvalues using IDR.
It turns out to be helpful to look at the corresponding polynomial description.
We already mentioned that Krylov subspace methods compute elements from the
Krylov subspace Km which can also be described as polynomials. The classiﬁcation
1 In [5] we name these relations in honor of Karl Hessenberg. He was to our knowledge theﬁrst who considered relations of the type AQn = Qn+1 Hn with a special unreduced extended
Hessenberg matrix Hn , see [7]. Usually the names of Lanczos [8, 9] or Arnoldi [1] are associated withsuch relations.

2

of Krylov methods given above can be rephrased using the language of polynomials.
The three classes of methods can be described using certain polynomials and polynomial interpolation, cf. [27, 28] for the general perturbed case. The Ritz values are
the eigenvalues of the Hessenberg matrix Cm and the harmonic Ritz values are theeigenvalues of the matrix (C†m Im )−1 , where Im denotes an extended identity matrix
which has an additional row of zeros attached to the bottom. The three approaches
implicitly compute polynomials with the following characteristics.
QOR: rm = Rm (A)r0 , where Rm (z) := det (Im − zC−1m ),xm = Lm [z −1 ](A)r0 , where Lm [z −1 ](z) interpolates
the function z −1 at the Ritz values,
QMR: rm = Rm (A)r0 , where Rm (z) := det (Im − zC†m Im ),xm = Lm [z −1 ](A)r0 , where Lm [z −1 ](z) interpolates
the function z −1 at the harmonic Ritz values,(m)Ritz-Galërkin: Unscaled Ritz vectors are given by vj = Am (θj , A)q1 ,where Am (θ, z) := (χm (θ) − χm (z))(θ − z)−1 , θ = z,
Cm sj = sj θj and χm (z) := det (zIm − Cm ).
These polynomial descriptions play a prominent rôle in determining the rate of convergence and understanding the behavior of Krylov subspace methods in ﬁnite precisionor subject to more general perturbations, e.g., those, which are introduced by an
inexact matrix vector product.
1.2. 1976–1980: IDR. In 1976 Sonneveld experimentally observed that for
B ∈ Cn×n and a given starting vector f0 ∈ Cn and f1 := Bf0 the following three-termrecurrence
fk+1 := B(fk − γk (fk − fk−1 )),

γk :=

pH fk
k − fk−1 )pH (falmost always terminates after 2n steps with the zero vector f2n = on [16, 17].
Analyzing this startling behavior, he discovered that the two consecutive vectors
f2j , f2j+1 constructed in this manner live in spaces Gj of shrinking dimensions. These
spaces are nowadays known as “Sonneveld spaces”. He called this property “Induced
Dimension Reduction” (IDR), and algorithms like the given three-term recurrence“IDR Algorithms”. The theorem he proved was the ﬁrst “IDR Theorem”.
Sonneveld ﬁrst made numerical experiments and then gave a rigorous proof for
his ﬁrst IDR Theorem. It is easy to see that apart from the ﬁrst two arbitrarily chosen
residuals the constructed residuals are in the B image of the space S := p⊥ . The
same argument proves that in general (observe that the ﬁrst two residuals f0 , f1 areusually not in S) for k  1
k
 k
⊥ 
⊥

B j (S) = + B −j H {p}= Kk (B −H , B −H p) .
f2k , f2k+1 ∈ Gk :=
j=1

j=1

Sonneveld proved that the dimensions of the spaces constructed are shrinking. This is
the essence of the ﬁrst IDR Theorem. He did not use the description as an orthogonal
complement of a Krylov subspace as it is done here and elsewhere [13, 11]. We remarkthat generically dim(Kn (B −H , B −H p)) = n. Using the Krylov subspace point of viewand the explicit orthogonalization against p before multiplication with B, we see thatindeed f2n = Bon = on .
3

This IDR Theorem can be turned into many IDR Algorithms. The three-termrecurrence
fk+1 = B(fk − γk (fk−1 − fk )),

where

γk =

pH fk
,
pH (fk−1 − fk )is an “implementation” of the Induced Dimension Reduction (IDR) Theorem. The
vectors constructed live in spaces of shrinking dimensions. Methods like this are called“IDR Algorithms”.
Another implementation by Sonneveld can be used to solve “genuine” linear systems. The idea is to rewrite the linear system to Richardson’s iteration form,
Ax = b⇒x = (I − A)x + b =: B x + b.

The classical Richardson’s iteration with a starting guess x0 is then given by
xk+1 = (I − A)xk + b.
With r0 := b − Ax0 , the Richardson’s iteration is carried out as follows:
xk+1 = xk + rk ,rk+1 = (I − A)rk .

In a Richardson-type IDR Algorithm, the second equation is replaced by the updaterk+1 = (I − A)(rk + γk (rk − rk−1 )),

γk =

pH rk
.
pH (rk−1 − rk )The update of the iterates has to be modiﬁed accordingly,−A(xk+1 − xk ) = rk+1 − rk = (I − A)(rk + γk (rk − rk−1 ) − rk
= (I − A)(rk − γk A(xk − xk−1 ) − rk
= −A(rk + γk (I − A)(xk − xk−1 ))⇔xk+1 − xk = rk + γk (I − A)(xk − xk−1 )= rk + γk (xk − xk−1 + rk − rk−1 ).

Sonneveld termed the outcome the “Primitive IDR Algorithm” [16]. For simplicity wewill refer to it as PIA. A very simple implementation of PIA is given in Algorithm 1.
The algorithm can be restated with less memory requirements when we allow tooverwrite the vectors and scalars. This variant can be found in Algorithm 2. In
Figure 1.1 we compare Richardson’s iteration and PIA. These pictures gives some
impressions of the “ﬁnite termination” property at step 2n and a feeling for theacceleration of Richardson’s iteration by PIA in ﬁnite precision.
Sonneveld never did use PIA, as he considered it to be too unstable, instead hewent on with a corresponding acceleration of the Gauß-Seidel method. In [17] heterms this method “Accelerated Gauß-Seidel” (AGS) and refers to it as “[t]he veryﬁrst IDR-algorithm [..]”, see page 6, Ibid. This part of the story took place “in the
background” in the year 1976. In September 1979 Sonneveld did attend the IUTAM
Symposium on Approximation Methods for Navier-Stokes Problems in Paderborn,
Germany. At this symposium he presented a new variant of IDR based on a variable
splitting I − ωj A, where ωj is ﬁxed for two steps and otherwise could be chosen
4input : A, b, x0 , p
output: rk , xk , γk
1 r0 = b − Ax0
2 x1 = x0 + r0
3 r1 = r0 − Ar0
4 for k = 1, 2 . . . do
5γk = pT rk /pT (rk−1 − rk )6
sk = rk + γk (rk − rk−1 )7
xk+1 = xk + γk (xk − xk−1 ) + sk
rk+1 = sk − Ask8
9 end
Algorithm 1: The most simple variant of PIA
input : A, b, x0 , p
output: r, x, γ
1 xold ← x0
2 rold ← b − Axold
3 xnew ← xold + rold
4 rnew ← rold − Arold
5 while not converged do
6γ ← pT rnew /pT (rold − rnew )7s ← rnew + γ(rnew − rold )8
xtmp ← xnew + γ(xnew − xold ) + s
9rtmp ← s − Asxold ← xnew , xnew ← xtmp
10
rold ← rnew , rnew ← rtmp
11
12 endAlgorithm 2: A better implementation of PIA
freely, but non-zero. This algorithm with a minimization of every second residual todetermine the new ωj is included in the proceedings from 1980 [24]. The connection
to Krylov methods, e.g., to BiCG and Lanczos, is also given there. We refer to it
as “Classical IDR” or simply “IDR”. The pseudo-code for the algorithm of IDR isgiven in Algorithm 3.
Algorithm 3 is the original IDR Algorithm from page 551 of [24]. It uses OrthoRes(1) in the ﬁrst step and a residual minimization every second step. The residuals
are the vectors −f2j . The ﬁnite termination property follows from a generalizationof the IDR Theorem based on the commutativity of the linear polynomials I − ωj A.
A numerical comparison of Richardson’s iteration, original IDR, and PIA is given in
Figure 1.2.
As is obvious from the derivation of PIA, AGS, and IDR, there are many brothersof classical and primitive IDR. In 1976 Sonneveld already considered the accelerationof Gauß-Seidel (AGS). Similarly, the IDR philosophy can be used as an accelerator forthe other classical splitting methods like Jacobi, SOR, SSOR, or Chebyshev, or evenmore general semi-iterative methods. Thus, IDR is a whole family of methods insteadof a single instant. We remark that the well known method BiCGStab [22, 21] is
5

0

5
10
matrix−vector multiplies
PIA for n = 5 and scaling

0

10

−10

10

0

5
10
matrix−vector multiplies

20

10

10

10

0

10

0

20
40
60
matrix−vector multiplies
PIA for n = 20 and scaling

0

10

−10

10

0

20
40
60
matrix−vector multiplies

PIA for n = 100 and no scalingtrue and updated residuals

−10

10

30

10

true and updated residuals

0

10

PIA for n = 20 and no scalingtrue and updated residuals

PIA for n = 5 and no scalingtrue and updated residuals

true and updated residuals
true and updated residuals

10

10

200

10

100

10

0

10

0

100
200
matrix−vector multiplies
PIA for n = 100 and scaling

0

10

−10

10

0

100
200
matrix−vector multiplies

Figure 1.1. A numerical comparison of PIA and Richardson’s iteration. The convergencecurves of PIA are in all cases, after some possible initial erratic behavior, below the convergencecurves of Richardson’s iteration. For small n, PIA manages to converge when Richardson’s iteration
actually diverges.

mathematically merely a reformulation of IDR and thus also an IDR algorithm.
Some of the methods based on IDR acceleration of classical splitting and semiiterative methods have been considered much more recently by Seiji Fujino et al. under
the names ISOR, IJacobi, IGS. The generalization of these methods to incorporate
more than one shadow vector seems very promising on distributed memory computers.
One has to look careful at the diﬀerence between preconditioning and using a variable
splitting before applying IDR acceleration or afterwards. The numerical behavior is
not very promising. But this picture changes, when we use more shadow vectors . . .
1.3. 2006–2010: IDR(s). In 2006 the rebirth of IDR took place. Almost
exactly 30 years after inventing IDR, Peter Sonneveld together with Martin van Gijzen
reconsidered IDR and introduced more shadow vectors: IDR(s)was born. They cameup with a variant called IDR(s) that used orthogonalization against a larger space,
where s denotes the dimension of that space. Algorithmically, the transition from
IDR to IDR(s) corresponds to replacing the single vector p ∈ Cn with a matrix orblock vector P ∈ Cn×s , 1  s  n.
To analyze IDR and IDR(s), we have to consider generalized Hessenberg decompositions (also referred to as rational Hessenberg decompositions) and to generalize
QOR, QMR and Ritz-Galërkin. We have to prove that the expressions for the iterates
and residuals based on polynomials are still valid. But: All these approaches extend
easily to generalized Hessenberg decompositions. The prototype IDR(s) algorithm
without the recurrences for xm and thus already slightly rewritten is given in Algorithm 4. This variant forms the basis for eigenvalue computations based on IDR in[5].
6input : A, b, x0 , p
output: fk , xk , γk , ωk
1 γ0 = 0
2 f0 = Ax0 − b
3 Δg0 = on
4 Δy0 = on
5 for k = 1, . . . do
6sk = fk−1 + γk−1 Δgk−1
tk = Ask
7
8
if k = 1 or k is even thenH9ωk = (tHk sk )/(tk tk )10else
ωk = ωk−1
11
12
end
13
Δxk = γk−1 Δyk−1 − ωk sk
Δ f k = γk−1 Δgk−1 − ωk tk
14
xk = xk−1 + Δxk
15
f k = f k−1 + Δfk
16
17
if k is even then18
Δyk = Δyk−1
19
Δgk = Δgk−1
20else
Δyk = Δxk
21
22
Δgk = Δfk23
end24
γk = −(pH fk )/(pH Δgk )25 endAlgorithm 3: Classical IDRWith mention a few remarks: We can start with any (simple) Krylov subspacemethod. The steps in the s-loop only diﬀer from the ﬁrst block in that no new ωjis computed. IDR(s)ORes is based on oblique projections and s + 1 consecutive
multiplications with the same linear factor I − ωj A. To understand IDR we haveto ﬁnd the underlying Hessenberg decomposition. We already noted that essential
features of Krylov subspace methods can be described by a Hessenberg decomposition
AQm = Qm+1 Hm = Qm Hm + qm+1 hm+1,m eT
m.

(1.8)

Here, Hm denotes an unreduced Hessenberg matrix. In the perturbed case, e.g., in ﬁnite precision and/or based on inexact matrix-vector multiplies, we obtain a perturbed
Hessenberg decomposition
AQm + Fm = Qm+1 Hm = Qm Hm + qm+1 hm+1,m eT
m.

(1.9)

The matrix Hm of the perturbed variant will, in general, still be unreduced. It turns
out that in case of IDR we have to consider generalized Hessenberg decompositions.
AQm Um = Qm+1 Hm = Qm Hm + qm+1 hm+1,m eT
m7(1.10)

−10

10

0

5
10
matrix−vector multiplies
RIP for n = 5 and scaling

0

10

−10

10

0

5
10
matrix−vector multiplies

0

10

0

20
40
60
matrix−vector multiplies
RIP for n = 20 and scaling

0

10

−10

10

0

20
40
60
matrix−vector multiplies

true and updated residuals

10

20

10

RIP for n = 100 and no scalingtrue and updated residuals

0

RIP for n = 20 and no scalingtrue and updated residuals

RIP for n = 5 and no scalingtrue and updated residuals

true and updated residuals
true and updated residuals

10

10

200

10

100

10

0

10

0

100
200
matrix−vector multiplies
RIP for n = 100 and scaling

0

10

−10

10

0

100
200
matrix−vector multiplies

Figure 1.2.
A numerical comparison of Richardson’s iteration, classical IDR, and PIA.
Classical IDR performs best for these small-scale test problems and converges in all cases to areasonable level of accuracy, even when PIA and Richardson’s iteration diverge.

and perturbed generalized Hessenberg decompositions
AQm Um + Fm = Qm+1 Hm = Qm Hm + qm+1 hm+1,m eT
m(1.11)

with upper triangular (possibly even singular) matrices Um . Generalized Hessenberg
decompositions correspond to an oblique projection of the pencil (A, I) to the pencil(Hm , Um ) as long as Qm+1 has full rank, H (A, I)Qm Um = Q
 H (AQm Um , Qm Um )Qmm
 H (Qm+1 H , Qm Um ) = (IT H , Um ) = (Hm , Um ),
=Q
m
m m
m

(1.12)

 H := IT Q† .
here Qmm m+1
In understanding IDR, we have to understand the concept of OrthoRes-type
methods. The entries of the Hessenberg matrices of Hessenberg decompositions are deﬁned in diﬀerent variations. Three well-known ways for implementing the QOR/QMR
approach are commonly denoted as OrthoRes, OrthoMin, and OrthoDir. The
prototype IDR(s) belongs to the class OrthoRes and uses short recurrences, therefore we refer to it as IDR(s)ORes. OrthoRes-type methods have a generalized
Hessenberg decomposition
ARm Um = Rm+1 H◦m = Rm H◦m + rm+1 h◦m+1,m eT
m,T
where eT H◦m = oTm , e = (1, . . . , 1), and the matrix
r0 2rm 2,...,
Rm+1 = r0 , . . . , rm = Qm+1 diag
q1 2
qm+1 28

(1.13)(1.14)

input : A, b, x0 , s, P
output: Rm+1 , cs+1 , cs+2 , . . ., ω1 , ω2 , . . .
1 r0 = b − Ax0
2 compute Rs+1 = R0:s = r0 , . . . , rs using, e.g., OrthoRes
3 ∇R1:s = ∇r1 , . . . , ∇rs = r1 − r0 , . . . , rs − rs−1
4 m ← s + 1, j ← 1
5 while not converged do
6
cm = (PH ∇Rm−s:m−1 )−1 PH rm−1
7
vm−1 = rm−1 − ∇Rm−s:m−1 cmcompute ωj8∇rm = −∇Rm−s:m−1 cm − ωj Avm−1
9rm = rm−1 + ∇rm
10
11
m←m+1

12
∇Rm−s:m−1 = ∇rm−s , . . . , ∇rm−1
13
for k = 1, . . . , s do
14
cm = (PH ∇Rm−s:m−1 )−1 PH rm−1
15
vm−1 = rm−1 − ∇Rm−s:m−1 cm∇rm = −∇Rm−s:m−1 cm − ωj Avm−1
16
rm = rm−1 + ∇rm
17
18
m←m+1

19
∇Rm−s:m−1 = ∇rm−s , . . . , ∇rm−1
20
end
21
j ←j+122 endAlgorithm 4: IDR(s)ORes
is diagonally scaled to be the matrix of residual vectors.
Next we depict the underlying generalized Hessenberg decomposition. We alreadynoted in the remarks following Algorithm 4 that the IDR recurrences of IDR(s)ORes
can be summarized by
vm−1 := rm−1 − ∇Rm−s:m−1 cm = Rm−s−1:m−1 ym(m)= (1 − γs) rm−1 +(m)(m)
s−1 (m)rm−s−1 ,
=1 (γs−+1 − γs− ) rn−−1 + γ1

1 · rm := (I − ωj A) vm−1 .
Here, m > s, and the index of the scalar ωj is deﬁned by
j :=

m
,
s+1

compare with the so-called “index functions” in [25]. Removing vm−1 from the recurrence we obtain the generalized Hessenberg decomposition= Rm+1 Y◦m ,ARm Ym D(m)ω

(1.15)

where non-zero elements in the upper triangular part of the columns of the upper
triangular matrix Ym and the unreduced extended upper Hessenberg matrix Y◦m9(m)are given by the elements in the vectors ym , the diagonal matrix Dω has s + 1
repeated copies of ωj along the diagonal and the lower diagonal of the unreduced
extended upper Hessenberg matrix Y◦m is given by minus ones, which makes it an
OrthoRes-type matrix, i.e., the columns sum to zero.
2. IDR(s)Eig. The IDR(s)Eig approach in [5] is based on the upper unreduced Hessenberg/upper triangular pencil of the generalized Hessenberg decomposition (1.15). For details we refer the reader to the forthcoming technical report [5]. We
simply sketch the mathematical transformations with the aid of small pictures basedon IDR(s)ORes with s = 3 for 12 steps.
2.1. Sonneveld pencil. In this small section we consider the Sonneveld pencil and Sonneveld matrix. The IDR(s)ORes pencil, the so-called Sonneveld pencil
(m)
◦
(Ym
, Ym Dω ), can be depicted by⎛×××× ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ⎞ ⎛×××× ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ⎞◦ ×××× ◦ ◦ ◦ ◦ ◦ ◦ ◦+×××× ◦ ◦ ◦ ◦ ◦ ◦ ◦⎜ ◦ +×××× ◦ ◦ ◦ ◦ ◦ ◦ ⎟ ⎜ ◦ ◦ ×××× ◦ ◦ ◦ ◦ ◦ ◦ ⎟⎜ ◦ ◦ +×××× ◦ ◦ ◦ ◦ ◦ ⎟ ⎜ ◦ ◦ ◦ ×××× ◦ ◦ ◦ ◦ ◦ ⎟⎟⎜⎟ ⎜⎜ ◦ ◦ ◦ +×××× ◦ ◦ ◦ ◦ ⎟ ⎜ ◦ ◦ ◦ ◦ ×××× ◦ ◦ ◦ ◦ ⎟⎜ ◦ ◦ ◦ ◦ +×××× ◦ ◦ ◦ ⎟ ⎜ ◦ ◦ ◦ ◦ ◦ ×××× ◦ ◦ ◦ ⎟⎜ ◦ ◦ ◦ ◦ ◦ +×××× ◦ ◦ ⎟ , ⎜ ◦ ◦ ◦ ◦ ◦ ◦ ×××× ◦ ◦ ⎟ .
⎜ ◦ ◦ ◦ ◦ ◦ ◦ +×××× ◦ ⎟ ⎜ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ×××× ◦ ⎟⎟⎜⎟ ⎜⎜ ◦ ◦ ◦ ◦ ◦ ◦ ◦ +××××⎟ ⎜ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ××××⎟⎝ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ +×××⎠ ⎝ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ×××⎠
◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ××
◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ +××◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦×
◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ +×
We observe a nice banded structure of this pencil. The eigenvalues can easily be
computed using the QZ algorithm. The data dependencies can be reduced to obtain
the approximate eigenvalues in a more stable manner. The upper triangular matrix(m)
Ym Dω could be inverted, which results in the Sonneveld matrix, a full unreduced
Hessenberg matrix. The QR algorithm on the Sonneveld matrix could be used toobtain the same approximations to eigenvalues.
2.2. Puriﬁed pencil. We know the eigenvalues ≈ roots of kernel polynomials
1/ωj . We are only interested in the other eigenvalues. The puriﬁed IDR(s)ORes(m)◦pencil (Ym, Um Dω ), that has only the remaining eigenvalues and some inﬁnite
ones as eigenvalues, can be depicted by⎛×××× ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ⎞ ⎛⎞××× ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦+×××× ◦ ◦ ◦ ◦ ◦ ◦ ◦
◦
××
◦
◦
◦
◦
◦
◦
◦
◦
◦
⎜ ◦ +×××× ◦ ◦ ◦ ◦ ◦ ◦ ⎟ ⎜ ◦ ◦ × ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ⎟⎜ ◦ ◦ +×××× ◦ ◦ ◦ ◦ ◦ ⎟ ⎜ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ⎟⎜⎟⎟ ⎜⎜ ◦ ◦ ◦ +×××× ◦ ◦ ◦ ◦ ⎟ ⎜ ◦ ◦ ◦ ◦ ××× ◦ ◦ ◦ ◦ ◦ ⎟⎜ ◦ ◦ ◦ ◦ +×××× ◦ ◦ ◦ ⎟ ⎜ ◦ ◦ ◦ ◦ ◦ ×× ◦ ◦ ◦ ◦ ◦ ⎟⎜ ◦ ◦ ◦ ◦ ◦ +×××× ◦ ◦ ⎟ , ⎜ ◦ ◦ ◦ ◦ ◦ ◦ × ◦ ◦ ◦ ◦ ◦ ⎟ .
⎜ ◦ ◦ ◦ ◦ ◦ ◦ +×××× ◦ ⎟ ⎜⎟⎜⎟ ⎜◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦⎟⎜ ◦ ◦ ◦ ◦ ◦ ◦ ◦ +××××⎟ ⎜ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ××× ◦ ⎟⎝ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ +×××⎠ ⎝ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ×× ◦ ⎠◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ +××◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦×◦◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ +×◦◦◦◦◦◦◦◦◦◦◦◦
We get rid of the inﬁnite eigenvalues using a change of basis (block Gauß elimination/Schur complement). The resulting deﬂated puriﬁed IDR(s)ORes pencil,(m)◦after the elimination step (YmGm , Um Dω ), where Gm denotes the block Gauß10

eliminator, can be depicted by⎛××××××× ◦ ◦ ◦ ◦ ◦ ⎞ ⎛⎞××× ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦+×××××× ◦ ◦ ◦ ◦ ◦
◦
××
◦
◦
◦
◦
◦
◦
◦
◦
◦
⎜ ◦ +××××× ◦ ◦ ◦ ◦ ◦ ⎟ ⎜ ◦ ◦ × ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ⎟⎜◦ ◦ ◦+◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦⎟ ⎜◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦⎟⎜⎟ ⎜⎟⎜ ◦ ◦ ++××××××× ◦ ⎟ ⎜ ◦ ◦ ◦ ◦ ××× ◦ ◦ ◦ ◦ ◦ ⎟⎜ ◦ ◦ ◦ ◦ +×××××× ◦ ⎟ ⎜ ◦ ◦ ◦ ◦ ◦ ×× ◦ ◦ ◦ ◦ ◦ ⎟⎜ ◦ ◦ ◦ ◦ ◦ +××××× ◦ ⎟ , ⎜ ◦ ◦ ◦ ◦ ◦ ◦ × ◦ ◦ ◦ ◦ ◦ ⎟ .
⎜◦ ◦ ◦ ◦ ◦ ◦ ◦+◦ ◦ ◦ ◦⎟ ⎜⎟⎜⎟ ⎜◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦⎟⎜ ◦ ◦ ◦ ◦ ◦ ◦ ++××××⎟ ⎜ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ××× ◦ ⎟⎝ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ +×××⎠ ⎝ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ×× ◦ ⎠◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ +××◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦×◦◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦+◦◦◦◦◦◦◦◦◦◦◦◦(m)◦
Gm we can get rid ofUsing Laplace expansion of the determinant of zUm Dω − Ym
the trivial constant factors corresponding to inﬁnite eigenvalues. This amounts to adeﬂation.

2.3. Deﬂated pencil. Let D denote an deﬂation operator that removes everys + 1th column and row from the matrix the operator is applied to. The deﬂated(m)◦Gm ), D(Um Dω )), canpuriﬁed IDR(s)ORes pencil, after the deﬂation step (D(Ymbe depicted by⎛×××××× ◦ ◦ ◦ ⎞ ⎛××× ◦ ◦ ◦ ◦ ◦ ◦ ⎞+××××× ◦ ◦ ◦◦ ×× ◦ ◦ ◦ ◦ ◦ ◦⎜ ◦ +×××× ◦ ◦ ◦ ⎟ ⎜ ◦ ◦ × ◦ ◦ ◦ ◦ ◦ ◦ ⎟⎜ ◦ ◦ +××××××⎟ ⎜ ◦ ◦ ◦ ××× ◦ ◦ ◦ ⎟⎜⎟⎟ ⎜⎜ ◦ ◦ ◦ +×××××⎟ , ⎜ ◦ ◦ ◦ ◦ ×× ◦ ◦ ◦ ⎟ .
⎜ ◦ ◦ ◦ ◦ +××××⎟ ⎜ ◦ ◦ ◦ ◦ ◦ × ◦ ◦ ◦ ⎟⎝ ◦ ◦ ◦ ◦ ◦ +×××⎠ ⎝ ◦ ◦ ◦ ◦ ◦ ◦ ×××⎠
◦ ◦ ◦ ◦ ◦ ◦ +××◦ ◦ ◦ ◦ ◦ ◦ ◦ ××
◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦×
◦ ◦ ◦ ◦ ◦ ◦ ◦ +×(m)

The block-diagonal matrix D(Um Dω ) has invertible upper triangular blocks andcan be inverted to expose the underlying Lanczos process, a Lanczos process with s
left-hand and one right-hand side. We refer to this process as BiORes(s,1), as this is
a process of OrthoRes-type. In the next section we show that IDR(s) is a Lanczos
process with multiple left-hand sides.
(m)2.4. BiORes(s,1). Inverting the block-diagonal matrix D(Um Dω )) gives analgebraic eigenvalue problem with a block-tridiagonal unreduced upper Hessenberg
matrix⎛×××××× ◦ ◦ ◦ ⎞+××××× ◦ ◦ ◦⎜ ◦ +×××× ◦ ◦ ◦ ⎟⎜ ◦ ◦ +××××××⎟⎜⎟◦−1
Gm ) · D(Um D(m)= ⎜ ◦ ◦ ◦ +×××××⎟ .
Lm := D(Ymω ))⎜ ◦ ◦ ◦ ◦ +××××⎟⎝ ◦ ◦ ◦ ◦ ◦ +×××⎠◦ ◦ ◦ ◦ ◦ ◦ +××◦ ◦ ◦ ◦ ◦ ◦ ◦ +×This is the matrix of the underlying BiORes(s, 1) process. The extended matrix
version Lm satisﬁes
AQm = Qm+1 Lm ,where the reduced residuals q js+k , k = 0, . . . , s − 1, j = 0, 1, . . ., with Ω0 (z) ≡ 1 andjΩj (z) = k=1 (1 − ωk z) are given by Ωj (A)q js+k = rj(s+1)+k . The reduced residuals
11are deﬁned by
Ωj (A)q js+k = rj(s+1)+k = (I − ωj A)vj(s+1)+k−1
and every vj(s+1)+k−1 is orthogonal to P. Thus, q js+k ⊥ Ωj−1 (AH )P. Using induction [13] one can prove that q js+k ⊥ Kj (AH , P); thus, this is a two-sided Lanczos
process with s left and one right starting vectors. This can more easily be provenusing the representations (S := P⊥ )G0 = K(A, r0 ), where K(A, r0 ) denotes the full Krylov subspace,
j−1
 j−1
⊥

Gj =
Ωk (A)−1 Ωj (A)(S) = + Ωj (A) −H Ωk (A) H {P}
k=0



k=0

⊥

⊥
= Ωj (A) −H Kj (AH , P)
= Ωj (A) Kj (AH , P)of the Sonneveld spaces. This has to be compared with Theorem 4.2 in [13] andwith Theorem 4.1 in [11] (here a similar result is obtained; the authors use a slightly
diﬀerent method of proof).
The ﬁrst equality
Gj =

j−1


Ωk (A)−1 Ωj (A)(S) =

k=0

j


(I − ωj A) · · · (I − ωk A)(S)

k=1

follows from the observations that the ﬁrst s + 1 residuals obviously are in G0 :=K(A, r0 ), the next s + 1 residuals (or any other vectors in G1 ) are in the I − ω1 A
image of S = P⊥ , the last s + 1 residuals are in the I − ωj A image of S = P⊥ , and,since they are computed as images of linear combinations of previous information, also
images of linear combinations of previously obtained images (I − ωj−1 A) · · · (I − ωk A)
of S = P⊥ .
The second equality
j−1

k=0

Ωk (A)−1 Ωj (A)(S) =

 j−1
⊥
+ Ωj (A) −H Ωk (A) H {P}
k=0

is based onBP⊥ = (B −H P)⊥andU ⊥ ∩ V ⊥ = (U ∪ V)⊥ = (U + V)⊥ .
The second relations are basic linear algebra. The ﬁrst relation follows from theobservation thatP ⊥ = v ∈ C n | PH v = o s⇒ BP⊥ = Bv ∈ Cn | PH v = os ,since, for invertible B,y ∈ BP⊥ ⇔ y = Bv ∧ PH v = os ⇔ PH v = PH B−1 y = (B−H P)H y = os .
12

The third and fourth equality
⊥ 
⊥
 j−1
+ Ωj (A) −H Ωk (A) H {P}= Ωj (A) −H Kj (AH , P)
k=0


⊥
= Ωj (A) Kj (AH , P)

are satisﬁed• since the polynomials Ωk (A), 0  k < j form a basis of the space of polynomials of degree less j, and• by the property proved above, respectively.
2.5. Generalizations of IDR(s). The residuals computed last in a complete
cycle are uniquely deﬁned. Based on the analysis of a possible breakdown of IDR(s),Sonneveld and van Gijzen came up with their new implementation IDR(s)BiO [23].
In this implementation of the IDR Theorem they use basis vectors g−1 , . . . , g−s ∈ Gj ,which are not simply residual diﬀerences, but linear combinations of these.
The new vectors vm and rm+1 are in this general setting given by the updates
vm = rm −s

gm−i γi =: rm − Gm cm ,i=1

rm+1 = (I − ωA)vm = rm − ωAvm −s


and thus,
gm−i γi ,i=1where cm is determined such that PH vm = o.
Recently, the relations between IDR(s) and BiCGStab() and combinations ofboth methods have been investigated.
• In [13] the authors derive diﬀerent implementations of ML(k)BiCGStab-likealgorithms.
• In [14] the authors combine the IDR philosophy with higher degree stabilization polynomials. The resulting method is named IDR(s)Stab(). The
approach is comparable to the one resulting in BiCGStab().
• In [20] the authors derive the algorithm GBiCGStab(s,L), which is similar to IDR(s)Stab(). In their own words: “Our algorithm is to theirs
what the Gauss-Seidel iteration is to the Jacobi iteration.” A predecessor of
GBiCGStab(s,L) seems to be the method called GIDR(s,L) in [19].
• In [12] the ideas behind BiCGStab2 [4] and GPBiCG [29] are considered.
The relation of IDR(s) to Petrov-Galërkin with a rational Krylov space motived themethod IDR-Ritz [11]. Another, simpler motivation is that the residual polynomials
should be designed to dampen the spectrum. Using the residual polynomial representation of IDR(s) we could choose the 1/ωj close but not equal to eigenvalues, atleast we should choose them in the ﬁeld of values of A.
The minimization used in IDR(s)ORes and IDR(s)BiO results in values ωjwhich are in the ﬁeld of values of A−H , thus Simoncini and Szyld suggest to use a
few steps of the Arnoldi method to compute some Ritz values, which are then used insome ordering as 1/ωj values. For real non-symmetric matrices this typically results
in an algorithm based on complex arithmetic in place of real arithmetic.
13

Last but not least: Certain old ideas have been reactivated. Sonneveld presented
the hitherto unpublished Accelerated Gauß-Seidel (AGS) method at the Kyoto Forumon Krylov Subspace Methods in 2008. Based on the algorithm in the proceedings,
Seiji Fujino et al. considered the acceleration of the classical splitting methods (Jacobi,
Gauß-Seidel and SOR). The resulting methods are called• IDR(s)-Jacobi (w/o adaptive tuning),• IDR(s)-GS,• IDR(s)-SOR.
These approaches result in a “tight packing” of preconditioning and Krylov subspace methods, compare with PIA. In most of these methods the ωj are ﬁxed by the
splitting chosen.
3. Parallelization of IDR(s) and IDR(s)Eig. To understand aspects of a
possible parallelization, we take a look at the structure of IDR(s). IDR(s) is atypical Krylov subspace method. Most known Krylov subspace methods are based on• matrix-vector products (“black-box” or with a sparse matrix),• some kind of (bi-)orthogonalization (frequently Gram-Schmidt),• solution of small linear systems,• dots (line search minimization),• axpys or gemvs (updates of vectors).
3.1. . . . a short introduction to IDR(s) parallelization. Here, we simply
sketch the IDR(s) variant IDR(s)BiO described in the technical report [23]. Its
Matlab source code can be downloaded fromhttp://ta.twi.tudelft.nl/nw/users/gijzen/IDR.html.
We remark that this version of IDR(s) is contained in release 3.0 of the IFISS package.
The structure of an unpreconditioned IDR(s)BiO is sketched as a pseudo-code in
Algorithm 5. A parallelization is easily possible, we could use, e.g.: {Sca,P}LAPACK,CUBLAS/CUDA, cloud computing, MPI & OpenMP, . . .
A naı̈ve CUBLAS implementation for instance would be based on the parallel evaluation of small sized problems and several calls to cublasSaxpy(), cublasSgemv(),
cublasSdot(), and cublasScopy(). But there is room for improvement. The implementation of the matrix-vector multiplication with A should be adjusted to the typeof matrix given, e.g.,• (fully optimized parallel) “black-box”,• reordering of a given sparse matrix using some heuristics (e.g., reverse CuthillMcKee, multi-color schemes, . . . ) for a better block-distribution to the nodes
and to reduce communication,• parallel implementation of an H-matrix format.
We could also allow for a (parallel) preconditioner.
Most important is the load balancing, especially in a non-homogeneous environment, since we have several synchronization points in the algorithm, namely the dotsfrom orthogonalization and (possibly) the computation of ωj and (possibly) the computation of the norm of the residual or backward error.
The naı̈ve CUBLAS implementation can be enhanced by using some other IDR(s)variant. The triangular bi-orthogonalisation scheme could be replaced:14• Instead of modiﬁed Gram-Schmidt we could use (iterated) classical GramSchmidt.
• Furthermore, we could adopt delayed re-orthogonalization [6] to the case of
iterated Gram-Schmidt-like bi-orthogonalisation.
• We could use other triangular basis transformations. This would result in fulls × s-systems, but the main computational eﬀort takes place elsewhere.
Similarly to the approaches used in parallel implementations of CG [10, 3] we could
try to minimize the number of synchronization barriers given by the orthogonalisation
against P, any other occurring (bi-)orthogonalization and the computation of ω byutilization of algebraic rewritings. The recent parallel variant in [2] is based on aminimization of synchronization barriers.
Next we give a few remarks about the general structure of IDR(s). We can quite
easily get rid of the synchronization points caused by the computation of the non-zeroscalars ωj . One idea is to precompute a certain amount of Ritz values like in [11] or
to use the CPU (or some nodes of the GPU) to compute rough approximations toeigenvalues based on the Sonneveld pencil [5].
If we use a method like IDR(s)-Jacobi, we typically have a lower convergence
rate, but have removed all synchronization points due to (bi)-orthogonalization of thebasis vectors in the Sonneveld spaces or their pre-images. This may give a speedupthat covers the price paid due to a slower convergence.
But: We can never get rid of the orthogonalization against P. This has to be
carried out in every IDR(s) method. This part of the algorithm should be optimized
using code adopted to the architecture we are working on.
We give a few remarks on the structure of P and the resulting cost of orthogonalization. We could adopt the generic choice for P, namely, using randomly generated
orthonormal columns to the type of problem. If the problem is the discrete version of
a 2D or 3D physical problem, we could use as test vectors pj the discretization of test
functions with compact support, e.g., we could use some (non-)overlapping Schwarzlike vectors. These vectors can be stored more compactly and could be distributed
to all nodes in a distributed memory architecture. To save memory, we could use
randomly generated vectors with only ±1 and 0. The orthogonalization against P isin this case based on the computation of two sums of subsets of the vectors and ﬁnally
one subtraction. Additional structure gives additional gain in eﬃciency. Both ideas
can be combined. One has to careful balance the beneﬁts and the risks of resulting
instabilities.
Last we give a few remarks about the structure of IDR(s)Eig and its possible
parallelization. As IDR(s)Eig is based on (a given variant) of IDR(s), the samecomments apply. We only have to store the vectors deﬁning the orthogonalization
against P (in every step one vector of length s), any triangular basis transformations(in every sweep of s + 1 steps a few s × s triangular matrices) and the ωj used fors + 1 consecutive steps. The computation of the eigenvalues should be performed onsome of the pencils using an adopted QZ algorithm working near the original band of
the banded pencil, the shift strategy should be chosen as to minimize communication
between diagonal blocks while retaining favorable convergence properties. It is notknown by now, which IDR(s)Eig algorithm is the one most stable. Thus, up to now,
nobody did consider to come up with a stable eigenvalue solver designed for the special
structure of pencils stemming from IDR(s) algorithms. Once a good candidate for an
15

IDR(s) algorithm suitable for stable eigenvalue computations is known, one can comeup with a parallel variant.
4. Conclusions. We gave a short introduction to Krylov subspace methods
with special emphasis on IDR. We sketched the IDR and IDR(s) families and indicated how to compute eigenvalues using one particular instance of IDR(s), namely,
IDR(s)ORes. We sketched the relation between the algorithm IDR(s)ORes for theapproximate solution of linear systems, the eigenvalue routine IDR(s)Eig based on it,
and a two-sided Lanczos process BiORes(s,1). We brieﬂy indicated how to parallelize
one member of the IDR(s) family, namely the more recent variant IDR(s)BiO. We
gave only a rather philosophical treatment of a parallel implementation of IDR(s)Eig.
Much work, both theoretical and practical, remains to be done.
Acknowledgments. I am indebted to my co-author Martin Gutknecht. It was(and still is) a pleasure to write the joint report [5]. Without Peter Sonneveld answering my e-mail dating to 2006 and sending me an e-mail informing me about the DCSE
symposium on IDR, June 3, 2009, nothing of this would have happened. I thank Martin van Gijzen for our interesting encounter in Harrachov 2007 and many interesting
discussions since then. I have to thank our whole IDR family. I am deeply gratefulfor the support by Prof. Seiji Fujino from the Research Institute for Information
Technology of Kyushu University who enabled my visit to Japan and thus my talks
at Kyushu University and at the University of Tokyo. Last but not least I would liketo thank Prof. Kengo Nakajima and Prof. Takahiro Katagiri at the Supercomputing
Division of the Information Technology Center of the University of Tokyo for making
my talk and this extended abstract possible.
Postscript. Just after the slides for the talk in Tokyo had been ﬁnished (5thFebruary 2010), Tijmen Collignon and Martin van Gijzen published a technical report[2] on parallelization of IDR(s) (9th February 2010). In this report some of the ideas
presented in the talk have been considered, namely, the use of specially structured P
and some algebraic rewritings to minimize the number of synchronization points.
REFERENCES[1] W. E. Arnoldi. The principle of minimized iterations in the solution of the matrix eigenvalue
problem. Quart. Appl. Math., 9:17–29, 1951.
[2] Tijmen P. Collignon and Martin B. van Gijzen. Fast solution of nonsymmetric linear systems
on Grid computers using parallel variants of IDR(s). Reports of the Department of Applied
Mathematical Analysis Report 10-05, Delft University of Technology, 2010.
[3] Eduardo F. D’Azevedo, Victor Eijkhout, and Charles H. Romine. A matrix framework for
conjugate gradient methods and some variants of CG with less synchronization overhead.
In PPSC, pages 644–646, 1993.
[4] Martin H. Gutknecht. Variants of BICGSTAB for matrices with complex spectrum. SIAM J.
Sci. Comput., 14(5):1020–1033, September 1993.
[5] Martin H. Gutknecht and Jens-Peter M. Zemke. Eigenvalue computations based on IDR, 2010.
Technical Report (to appear 2010).
[6] V. Hernández, J. E. Román, and A Tomás. A parallel variant of the Gram-Schmidt process with reorthogonalization. In G. R. Joubert, W. E. Nagel, F. J. Peters, O. G. Plata,
P. Tirado, and E. L. Zapata, editors, Proceedings of the International Conference on Parallel Computing (ParCo 2005), volume 33, pages 221–228. Central Institute for Applied
Mathematics, Jülich, Germany, 2006.
[7] Karl Hessenberg. Behandlung linearer Eigenwertaufgaben mit Hilfe der Hamilton-Cayleyschen
Gleichung. Numerische Verfahren, Bericht 1, Institut für Praktische Mathematik (IPM),
16

Technische Hochschule Darmstadt, July 1940. Scanned report and biographical sketch of
Karl Hessenberg’s life online available at http://www.hessenberg.de/karl1.html.
[8] C. Lanczos. An iteration method for the solution of the eigenvalue problem of linear diﬀerential
and integral operators. J. Res. Nat. Bureau Standards, 45:255–281, 1950.
[9] C. Lanczos. Solution of systems of linear equations by minimized iterations. J. Res. Nat.
Bureau Standards, 49:33–53, 1952.
[10] Gérard Meurant. Multitasking the conjugate gradient method on the CRAY X-MP/48. Parallel
Comput., 5(3):267–280, 1987.
[11] Valeria Simoncini and Daniel Szyld. Interpreting IDR as a Petrov-Galerkin method. Report
09-10-22, Dipartimento di Matematica, Università di Bologna and Department of Mathematics, Temple University, Philadelphia, 2009.
[12] Gerard L. G. Sleijpen and Kuniyoshi Abe, 2010. Publication in preparation (January 2010).
[13] Gerard L. G. Sleijpen, Peter Sonneveld, and Martin B. van Gijzen. Bi-CGSTAB as an induced
dimension reduction method. Reports of the Department of Applied Mathematical Analysis
Report 08-07, Delft University of Technology, 2008. ISSN 1389-6520.
[14] Gerard L. G. Sleijpen and Martin B. van Gijzen. Exploiting BiCGstab() strategies to induce
dimension reduction. Reports of the Department of Applied Mathematical Analysis Report
09-02, Delft University of Technology, 2009. ISSN 1389-6520.
[15] P. Sonneveld and M. B. van Gijzen. IDR(s): a family of simple and fast algorithms for solving
large nonsymmetric systems of linear equations. Report 07-07, Department of Applied
Mathematical Analysis, Delft University of Technology, 2007.
[16] Peter Sonneveld. History of IDR: an example of serendipity. PDF ﬁle sent by Peter Sonneveld
on Monday, 24th of July 2006, July 2006. 8 pages; evolved into [17].
[17] Peter Sonneveld. AGS-IDR-CGS-BiCGSTAB-IDR(s): The circle closed. A case of serendipity.
In Proceedings of the International Kyoto Forum 2008 on Krylov subspace methods, pages
1–14, September 2008.
[18] Peter Sonneveld and Martin B. van Gijzen. IDR(s): A family of simple and fast algorithms
for solving large nonsymmetric systems of linear equations. SIAM Journal on Scientiﬁc
Computing, 31(2):1035–1062, 2008. Received Mar. 20, 2007.
[19] Masaaki Tanio and Masaaki Sugihara. GIDR(s,L): generalized IDR(s). In The 2008 annual
conference of the Japan Society for Industrial and Applied Mathematic, pages 411–412,
Chiba, Japan, September 2008. (In Japanese).
[20] Masaaki Tanio and Masaaki Sugihara. GBi-CGSTAB(s, L): IDR(s) with higher-order stabilization polynomials. Technical Report METR 2009-16, Department of Mathematical
Informatics, Graduate School of information Science and Technology, University of Tokio,
April 2009.
[21] H. A. van der Vorst. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for thesolution of nonsymmetric linear systems. SIAM J. Sci. Statist. Comput., 13:631–644, 1992.
Received May 21, 1990.
[22] H. A. van der Vorst and P. Sonneveld. CGSTAB, a more smoothly converging variant of CG-S.
Report 90-50, Department of Mathematics and Informatics, Delft University of Technology,
1990.
[23] Martin B. van Gijzen and Peter Sonneveld. An elegant IDR(s) variant that eﬃciently exploits
bi-orthogonality properties. Reports of the Department of Applied Mathematical Analysis
Report 08-21, Delft University of Technology, 2008. ISSN 1389-6520.
[24] P. Wesseling and P. Sonneveld. Numerical experiments with a multiple grid and a preconditioned Lanczos type method. In Approximation Methods for Navier-Stokes Problems,volume 771 of Lecture Notes in Mathematics, pages 543–562. Springer, 1980.
[25] Man-Chung Yeung and Daniel Boley. Transpose-free multiple Lanczos and its application in
Padé approximation. J. Comput. Appl. Math., 177(1):101–127, 2005.
[26] Man-Chung Yeung and Tony F. Chan. ML(k)BiCGSTAB: a BiCGSTAB variant based on
multiple Lanczos starting vectors. SIAM J. Sci. Comput., 21(4):1263–1290, 1999. Received
May 16, 1997, electr. publ. Dec. 15, 1999.
[27] Jens-Peter M. Zemke. Hessenberg eigenvalue-eigenmatrix relations. Linear Algebra and its
Applications, 414(2–3):589–606, 2006.
[28] Jens-Peter M. Zemke. Abstract perturbed Krylov methods. Linear Algebra and its Applications,424(2–3):405–434, 2007.
[29] Shao-Liang Zhang. GPBi-CG: generalized product-type methods based on Bi-CG for solving
nonsymmetric linear systems. SIAM Journal on Scientiﬁc Computing, 18(2):537–551,
1997.

17

input : A, b, x0 , s, P
output: Rm+1 , cs+1 , cs+2 , . . ., ω1 , ω2 , . . .
1 x = x0
2 r = b − Ax
3 ω =1T
4 PT = P
5 G = zeros(n, s)
6 U = zeros(n, s)
7 M = eye(s)8 while not converged do
9
PT r = PT · r
10
for k = 1, . . . , s do
11
% Solve small system and make v orthogonal to P:
12
c = M(k : s, k : s)−1 PT r(k : s)
13v = r − G(:, k : s)c
14
U(:, k) = U(:, k : s)c + ωv
15
% Compute G(:,k) = AU(:,k)
16G(:, k) = AU(:, k)17% Bi-Orthogonalize the new basis vectors:
18
for j = 1, . . . , k − 1 do
19
α = (PT(j, :)G(:, k))/M(j, j)
20
G(:, k) = G(:, k) − αG(:, j)21
U(:, k) = U(:, k) − αU(:, j)
22
end
23
% New column of M = P’*G (first k-1 entries are zero)24
M(k : s, k) = PT(k : s, :)G(:, k)25% Make r orthogonal to p j, j = 1,...,k26
β = PT r(k)/M(k, k)
27
r = r − βG(:, k)28x = x + βU(:, k)29% New PTr = P’*r (first k components are zero)30
if k < s then31
PT r(k + 1 : s) = PT r(k + 1 : s) − βM(k + 1 : s, k)
32
end
33
end
34
% Note: r is already perpendicular to P so v = r
35
v=r
36
t = Av
37
select or compute ω38
r = r − ωt39x = x + ωv
40 endAlgorithm 5: IDR(s)BiO

18