# Some Generalizations and Modifications of Iterative Methods for Solving Large Sparse Symmetric Indefinite Linear Systems

Authors: Yu-Chien Li; Jen-Yuan Chen; David R. Kincaid
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We discuss a variety of iterative methods that are based on the Arnoldi process for solving large sparse symmetric indefinite linear systems. We describe the SYMMLQ and SYMMQR methods, as well as generalizations and modifications of them. Then, we cover the Lanczos/MSYMMLQ and Lanczos/MSYMMQR methods, which arise from a double linear system. We present pseudocodes for these algorithms.


---

## Body

## 1. Introduction

Frequently, when computing numerical solutions of partial differential equations, one needs to solve systems of very large sparse linear algebraic equations of the form(1)
A
x
=
b
,
where A is an n
×
n matrix, b is an n
×
1 vector, and one seeks a numerical solution vector x or a good approximation of it. Particularly for large linear systems arising from partial differential equations in three dimensions, well-known direct methods, such as Gaussian elimination, may become prohibitively expensive in terms of both computer storage and computer time. On the other hand, a variety of iterative methods may avoid these difficulties.For linear systems involving symmetric positive definite (SPD) matrices, the conjugate gradient (CG) method (and variations of it) may work well. On the other hand, when solving linear systems, where the coefficient matrixA is symmetric indefinite, the choice of a suitable iterative method isnot at all clear. On the other hand, the SYMMLQ and MINRES methods have been shown to be useful in certain situations (see Paige and Saunders [1]). For nonsymmetric systems, Saad and Schultz [2] generalized the MINRES method to obtain the GMRES method.In Section2, we review the Arnoldi process. In Sections 3 and 4, we describe the SYMMLQ and SYMMQR methods. Then we can generalize them, in Section 5, and we outline the modified SYMMLQ method, in Section 6. Next, in Section 7, we discuss applying the MSYMMLQ and MSYMMQR methods applied to a double linear system. Finally, we present pseudocodes in Sections 8–11.


## 2. Arnoldi Process

We begin with a review of the Arnoldi process.Theorem 1.
Suppose thatA is an n
×
n symmetric matrix. One can generate orthonormal vectors w
(
0
)
,
w
(
1
)
,
…
,
w
(
n
-
2
)
,
w
(
n
-
1
) using this short-term recurrence
(2)
w
~
(
j
+
1
)
≡
A
w
(
j
)
-
α
j
w
(
j
)
-
β
j
w
(
j
-
1
)
(
0
≤
j
≤
n
-
2
)
w
(
j
+
1
)
=
(
1
σ
j
+
1
)
w
~
(
j
+
1
)
,
w
h
e
r
e
σ
j
+
1
=
〈
w
~
(
j
+
1
)
,
w
~
(
j
+
1
)
〉
,
where
(3)
α
j
=
〈
A
w
(
j
)
,
w
(
j
)
〉
,
β
j
=
〈
A
w
(
j
)
,
w
(
j
-
1
)
〉
.
Here, one assumes that w
(
-
1
)
=
0 and σ
j
≠
0, for all j. Then the following properties hold, for (0
≤
i, j
≤
n
-
1):
(4)
〈
w
(
i
)
,
w
(
j
)
〉
=
δ
i
j
,
β
j
=
σ
j
.Proof.
If we letw
~
(
0
)
≡
r
(
0
)
=
b
-
A
u
(
0
), then the subspace
(5)
Span
⁡
{
w
(
0
)
,
w
(
1
)
,
…
,
w
(
n
-
2
)
,
w
(
n
-
1
)
}
is equivalent to the Krylov subspace
(6)
𝒦
n
(
r
(
0
)
,
A
)
=
Span
⁡
{
r
(
0
)
,
A
r
(
0
)
,
…
,
A
n
-
2
r
(
0
)
,
A
n
-
1
r
(
0
)
}
.
We obtain(7)
β
j
=
〈
A
w
(
j
)
,
w
(
j
-
1
)
〉
(
by
(
3
)
)
=
〈
w
(
j
)
,
A
w
(
j
-
1
)
〉
(
A
T
=
A
)
=
〈
w
(
j
)
,
(
σ
j
w
(
j
)
+
α
j
-
1
w
(
j
-
1
)
+
β
j
-
1
w
(
j
-
2
)
)
〉
(
by
(
2
)
)
=
σ
j
since 〈
w
(
i
)
,
w
(
j
)
〉
=
δ
i
j.From Theorem1, in matrix form, it follows that
(8)
A
W
n
-
1
=
W
n
-
1
T
n
+
σ
n
w
(
n
)
e
n
T
=
W
n
T
~
n
+
1
,
where
(9)
W
n
-
1
=
[
w
(
0
)
,
w
(
1
)
,
…
,
w
(
n
-
1
)
,
w
(
n
-
2
)
,
w
(
n
-
1
)
]
n
×
n
e
n
T
=
[
0,0
,
0
,
…
,
0,0
,
1
]
1
×
n
,
T
n
≡
[
α
0
β
1
σ
1
α
1
β
2
σ
2
α
2
β
3
⋱
⋱
⋱
σ
n
-
1
α
n
-
3
β
n
-
2
σ
n
-
2
α
n
-
2
β
n
-
1
σ
n
-
1
α
n
-
1
]
n
×
n
,
T
~
n
+
1
≡
[
α
0
β
1
σ
1
α
1
β
2
σ
2
α
2
β
3
⋱
⋱
⋱
σ
n
-
3
α
n
-
3
β
n
-
2
σ
n
-
2
α
n
-
2
β
n
-
1
σ
n
-
1
α
n
-
1
σ
n
]
(
n
+
1
)
×
n
.Example 2.
We illustrate Theorem1 for the case n
=
3.From (2) and (3), we have
(10)
σ
1
w
(
1
)
=
w
~
(
1
)
=
A
w
(
0
)
-
α
0
w
(
0
)
-
β
0
w
(
-
1
)
,
σ
2
w
(
2
)
=
w
~
(
2
)
=
A
w
(
1
)
-
α
1
w
(
1
)
-
β
1
w
(
0
)
,
σ
3
w
(
3
)
=
w
~
(
3
)
=
A
w
(
2
)
-
α
2
w
(
2
)
-
β
2
w
(
1
)
.
Consequently, we obtain, since w
(
-
1
)
=
0,
(11)
A
W
2
=
A
[
w
(
0
)
w
(
1
)
w
(
2
)
]
=
[
A
w
(
0
)
A
w
(
1
)
A
w
(
2
)
]
=
[
β
0
w
(
-
1
)
+
α
0
w
(
0
)
+
σ
1
w
(
1
)
β
1
w
(
0
)
+
α
1
w
(
1
)
+
σ
2
w
(
2
)
β
2
w
(
1
)
+
α
2
w
(
2
)
+
σ
3
w
(
3
)
]
T
=
[
w
(
0
)
w
(
1
)
w
(
2
)
w
(
3
)
]
×
[
α
0
β
1
0
σ
1
α
1
β
2
0
σ
2
α
2
0
0
σ
3
]
4
×
3
=
[
W
2
w
(
3
)
]
[
T
3
0
0
σ
3
]
.
So we obtain
(12)
A
W
2
=
W
2
T
3
+
σ
3
w
(
3
)
e
3
T
=
W
3
T
~
4


## 3. SYMMLQ Method

We chooseu
(
n
), such that u
(
n
)
-
u
(
0
)
∈
𝒦
n
(
r
(
0
)
,
A
). Hence, we have
(13)
u
(
n
)
=
u
(
0
)
+
W
n
-
1
y
(
n
)
,
r
(
n
)
=
r
(
0
)
-
A
W
n
-
1
y
(
n
)
.Imposing the Galerkin conditionr
(
n
)
⊥
𝒦
n
(
r
(
0
)
,
A
), we obtain
(14)
(
r
(
n
)
)
T
W
n
-
1
=
0
,
W
n
-
1
T
r
(
n
)
=
0
,
W
n
-
1
T
r
(
0
)
=
W
n
-
1
T
A
W
n
-
1
y
(
n
)
(
by
(
13
)
)
.
We obtain
(15)
T
n
y
(
n
)
=
σ
0
e
1
because
(16)
W
n
-
1
T
A
W
n
-
1
=
T
n
,
W
n
-
1
T
r
(
0
)
=
σ
0
e
1
,
where
σ
0
=
∥
r
(
0
)
∥
.
Instead of solving for y
(
n
) directly from the triangular linear system (15), Paige and Saunders [1] factorize the matrix T
n into a lower triangular matrix with bandwidth three (resulting in the SYMMLQ method). Also, we have
(17)
T
n
Q
n
-
1
=
L
^
n
=
[
γ
0
δ
1
γ
1
ε
2
δ
2
γ
2
⋱
⋱
⋱
ε
n
-
3
δ
n
-
3
γ
n
-
3
ε
n
-
2
δ
n
-
2
γ
n
-
2
ε
n
-
1
δ
n
-
1
γ
^
n
-
1
]
n
×
n
,
where Q
n
-
1
≡
Q
1,2
Q
2,3
⋯
Q
n
-
1
,
n is an orthogonal matrix, and
(18)
Q
i
,
i
+
1
=
[
1
1
⋱
1
c
i
s
i
s
i
-
c
i
1
1
⋱
1
]
n
×
n
,
where c
i
2
+
s
i
2
=
1. Since T
n
Q
n
-
1
=
L
^
n, we have
(19)
T
n
Q
n
-
1
Q
n
-
1
-
1
y
(
n
)
=
L
^
n
Q
n
-
1
-
1
y
(
n
)
=
σ
0
e
1
.
Letting
(20)
z
^
(
n
)
=
Q
n
-
1
-
1
y
(
n
)
,
then
(21)
L
^
n
z
^
(
n
)
=
σ
0
e
1
.
Next letting
(22)
x
^
(
n
)
=
u
(
0
)
+
W
n
-
1
y
(
n
)
,
we have
(23)
x
^
(
n
)
=
u
(
0
)
+
W
n
-
1
Q
n
-
1
Q
n
-
1
-
1
y
(
n
)
=
u
(
0
)
+
W
n
-
1
Q
n
-
1
z
^
(
n
)
.
Defining
(24)
V
^
n
=
W
n
-
1
Q
n
-
1
z
^
n
=
[
ζ
0
,
ζ
1
,
…
,
ζ
n
-
2
,
ζ
^
n
-
1
]
1
×
n
T
,
we have
(25)
x
^
(
n
)
=
x
(
0
)
+
V
^
n
z
^
n
,
where
(26)
V
^
n
=
[
v
(
0
)
,
v
(
1
)
,
…
,
v
(
n
-
2
)
,
v
^
(
n
-
1
)
]
n
×
n
.
We let
(27)
z
n
=
[
ζ
0
,
ζ
1
,
…
,
ζ
n
-
2
,
ζ
n
-
1
]
1
×
n
T
,
V
n
=
[
v
(
0
)
,
v
(
1
)
,
…
,
v
(
n
-
2
)
,
v
(
n
-
1
)
]
n
×
n
,
L
n
z
(
n
)
=
σ
0
e
1
,
where
(28)
L
n
=
[
γ
0
δ
1
γ
1
ε
2
δ
2
γ
2
⋱
⋱
⋱
ε
n
-
3
δ
n
-
3
γ
n
-
3
ε
n
-
2
δ
n
-
2
γ
n
-
2
ε
n
-
1
δ
n
-
1
γ
n
-
1
]
n
×
n
,
γ
n
-
1
=
γ
^
n
-
1
2
+
β
n
-
1
2
,
c
n
=
γ
^
n
-
1
γ
n
-
1
,
s
n
=
β
n
γ
n
-
1
.
From (21) and (28), we have ζ
n
-
1
=
(
ζ
^
n
-
1
/
ζ
n
-
1
)
ζ
^
n
-
1
=
c
n
ζ
^
n
-
1. Since
(29)
V
^
n
+
1
=
[
v
(
0
)
,
v
(
1
)
,
…
,
v
(
n
-
1
)
,
v
^
(
n
)
]
=
W
n
Q
n
=
[
W
n
-
1
,
w
(
n
)
]
Q
1,2
,
Q
2,3
,
…
,
Q
n
-
1
,
n
Q
n
,
n
+
1
=
[
V
^
n
,
w
(
n
)
]
Q
n
,
n
+
1
=
[
v
(
0
)
,
v
(
1
)
,
…
,
v
^
(
n
-
1
)
,
w
(
n
)
]
Q
n
,
n
+
1
,
we have
(30)
v
(
n
-
1
)
=
c
n
v
^
(
n
-
1
)
+
s
n
w
(
n
)
,
v
^
(
n
)
=
s
n
v
^
(
n
-
1
)
-
c
n
w
(
n
)
.
If β
n
≠
0, then L
n is nonsingular. We can find z
(
n
) by solving
(31)
L
n
z
(
n
)
=
σ
0
e
1
,
u
(
n
)
=
u
(
0
)
+
V
n
z
(
n
)
=
u
(
n
-
1
)
+
ζ
n
-
1
v
(
n
-
1
)
.


## 4. SYMMQR Method

We chooseu
(
0
) such that u
(
n
)
-
u
(
0
)
∈
𝒦
n
(
r
(
0
)
,
A
). Hence, we have
(32)
u
(
n
)
=
u
(
0
)
+
W
n
-
1
y
(
n
)
,
r
(
n
)
=
r
(
0
)
-
A
W
n
-
1
y
(
n
)
.
Imposing the Galerkin condition r
(
n
)
⊥
𝒦
n
(
r
(
0
)
,
A
), as before, we obtain
(33)
(
r
(
n
)
)
T
W
n
-
1
=
0
,
W
n
-
1
T
r
(
n
)
=
0
,
W
n
-
1
T
r
(
0
)
=
W
n
-
1
T
A
W
n
-
1
y
(
n
)
(
by
(
32
)
)
.
Since
(34)
W
n
-
1
T
A
W
n
-
1
=
T
n
W
n
-
1
T
r
(
0
)
=
σ
0
e
1
,
where
σ
0
=
∥
r
(
0
)
∥
,
we have
(35)
T
n
y
(
n
)
=
σ
0
e
1
.
Instead of solving for y
(
n
) directly from the triangular system (35), Paige and Saunders [1] factorized the matrix T
n into a lower triangular matrix with bandwidth three.We can use a different factorization ofT
n to obtain a slightly different method, which is called the SYMMQR method. We multiply the matrix T
n by an orthogonal matrix on the left-hand side instead of the right-hand side. We have
(36)
Q
n
-
1
T
n
y
^
n
=
Q
n
-
1
σ
0
e
1
,
where
(37)
Q
n
-
1
T
n
=
R
^
=
[
γ
0
δ
1
ε
2
γ
1
δ
2
ε
3
γ
2
δ
3
ε
4
⋱
⋱
⋱
γ
n
-
3
δ
n
-
2
ε
n
-
1
γ
n
-
2
δ
n
-
1
γ
^
n
-
1
]
n
×
n
.
We obtain the matrix Q
n
-
1
≡
Q
n
-
1
,
n
Q
n
-
2
,
n
-
1
⋯
Q
1,2, where
(38)
Q
i
,
i
+
1
=
[
1
1
⋱
1
c
i
-
s
i
s
i
c
i
1
1
⋱
1
]
n
×
n
,
with c
i
2
+
s
i
2
=
1 being the Givens rotation. Letting y
^
n be the solution of
(39)
R
^
y
^
n
=
Q
n
-
1
σ
0
e
1
,
then we have
(40)
x
^
(
n
)
=
u
(
0
)
+
W
n
-
1
y
^
n
=
u
(
0
)
+
W
n
-
1
R
^
n
-
1
Q
n
-
1
σ
0
e
1
which satisfies the Galerkin condition r
^
n
⊥
𝒦
n
(
r
(
0
)
,
A
), where r
^
n
=
b
-
A
x
^
(
n
). We note that γ
^
n
-
1 isnot always nonzero and, thus, R
^
n might be singular. We assume that R
^
n is nonsingular and then we define
(41)
P
^
n
-
1
=
W
n
-
1
R
^
n
-
1
,
where
(42)
P
^
n
-
1
=
[
p
(
0
)
,
p
(
1
)
,
…
,
p
(
n
-
2
)
,
p
^
(
n
-
1
)
]
n
×
n
,
z
^
(
n
)
=
Q
n
-
1
σ
0
e
1
=
[
ζ
0
,
ζ
1
,
…
,
ζ
n
-
2
,
ζ
^
n
-
1
]
T
.
We have
(43)
x
^
(
n
)
=
u
(
0
)
+
P
^
n
-
1
z
^
(
n
)
.For the next iteratex
^
(
n
+
1
), we need to solve
(44)
T
n
+
1
y
^
(
n
+
1
)
=
σ
0
e
1
,
where
(45)
T
n
+
1
=
[
α
0
β
1
β
1
α
1
β
2
β
2
α
2
β
3
⋱
⋱
⋱
β
n
-
2
α
n
-
2
β
n
-
1
β
n
-
1
α
n
-
1
β
n
β
n
α
n
]
n
+
1
×
n
.
Applying the Givens rotation Q
n
-
1 to both sides of (45), we have
(46)
Q
n
-
1
T
n
+
1
=
[
γ
0
δ
1
ε
2
γ
1
δ
2
ε
3
γ
2
δ
3
ε
4
⋱
⋱
⋱
γ
n
-
3
δ
n
-
2
ε
n
-
1
γ
n
-
2
δ
n
-
1
ε
n
γ
^
n
-
1
ψ
β
n
α
n
]
n
+
1
×
n
,
where [
0,0
,
0
,
…
,
0,0
,
ε
n
,
ψ
,
α
n
]
T
=
Q
n
-
1
,
n
[
0,0
,
0
,
…
,
0,0
,
β
n
,
α
n
]
T.To eliminateβ
n, we compute the nth Given rotation Q
n
,
n
+
1 by
(47)
γ
n
-
1
=
γ
^
n
-
1
2
+
β
n
2
,
c
n
=
γ
^
n
-
1
γ
n
-
1
,
s
n
=
-
β
n
γ
n
-
1
.
By multiplying Q
n
,
n
+
1 times Q
n
-
1
T
n
-
1 and times Q
n
-
1
σ
0
e
1, we have
(48)
Q
n
,
n
+
1
(
Q
n
-
1
T
n
+
1
)
=
R
^
n
+
1
=
[
γ
0
δ
1
ε
2
γ
1
δ
2
ε
3
γ
2
δ
3
ε
4
⋱
⋱
⋱
γ
n
-
3
δ
n
-
2
ε
n
-
1
γ
n
-
2
δ
n
-
1
ε
n
γ
n
-
1
δ
n
γ
^
n
]
n
+
1
×
n
,
Q
n
,
n
+
1
(
Q
n
-
1
σ
0
e
1
)
=
z
^
(
n
+
1
)
=
[
ζ
0
,
ζ
1
,
…
,
ζ
n
-
1
,
ζ
^
n
]
T
,
where
(49)
[
0,0
,
0
,
…
,
0,0
,
ε
n
,
δ
n
,
γ
^
n
]
T
=
Q
n
,
n
+
1
[
0,0
,
0
,
…
,
0,0
,
ε
n
,
ψ
,
α
n
]
T
,
ζ
n
-
1
=
c
n
ζ
^
n
-
1
,
ζ
^
n
=
s
n
ζ
^
n
-
1
.
Let
(50)
R
n
=
[
γ
0
δ
1
ε
2
γ
1
δ
2
ε
3
γ
2
δ
3
ε
4
⋱
⋱
⋱
γ
n
-
3
δ
n
-
2
ε
n
-
1
γ
n
-
2
δ
n
-
1
γ
n
-
1
]
n
×
n
.
We define z
(
n
)
=
[
ζ
0
,
ζ
1
,
…
,
ζ
n
-
2
,
ζ
n
-
1
]
T. Since β
n
=
|
|
w
(
n
)
|
|
≠
0, then γ
n
≠
0 and R
n is nonsingular. We can solve for y
(
n
) from
(51)
R
n
y
(
n
)
=
z
(
n
)
.
We discuss the case β
n
=
∥
w
(
n
)
∥ later.Consider solving the least square problem involvingy
(
n
) minimizing ∥
T
~
n
+
1
y
(
n
)
-
σ
0
e
1
∥, where
(52)
T
~
n
+
1
=
[
α
0
β
1
β
1
α
1
β
2
β
2
α
2
β
3
⋱
⋱
⋱
β
n
-
2
α
n
-
2
β
n
-
1
β
n
-
1
α
n
-
1
β
n
]
n
+
1
×
n
.
We have
(53)
∥
Q
n
,
n
+
1
Q
n
-
1
(
T
~
n
+
1
y
(
n
)
-
σ
0
e
1
)
∥
=
∥
[
R
n
0,0
,
…
,
0,0
]
n
+
1
×
n
y
(
n
)
-
[
ζ
0
ζ
1
ζ
1
⋮
ζ
n
-
1
ζ
^
n
]
∥
.
Hence, the solution y
(
n
) from R
n
y
(
n
)
=
z
(
n
) minimizes |
|
T
~
n
+
1
y
(
n
)
-
σ
0
e
1
|
| and ζ
~
n
=
|
|
T
~
n
+
1
y
(
n
)
-
σ
0
e
1
|
|.Let(54)
u
(
n
)
=
u
(
0
)
+
W
n
-
1
y
(
n
)
=
u
(
0
)
+
W
n
-
1
R
n
-
1
z
(
n
)
,
where
(55)
P
n
-
1
R
n
=
W
n
-
1
,
P
n
-
1
=
[
p
(
0
)
,
p
(
1
)
,
…
,
p
(
n
-
2
)
,
p
(
n
-
1
)
]
.
We have
(56)
p
(
n
-
1
)
=
(
1
γ
n
-
1
)
[
w
(
n
-
1
)
-
ε
n
-
1
p
(
n
-
3
)
-
δ
n
-
1
p
(
n
-
2
)
]
,
u
(
n
)
=
u
(
0
)
+
P
n
-
1
z
(
n
)
=
u
(
n
-
1
)
+
γ
n
-
1
p
(
n
-
1
)
.
Since
(57)
x
^
(
n
+
1
)
=
u
(
0
)
+
W
n
y
^
(
n
+
1
)
,
we obtain
(58)
x
^
(
n
+
1
)
=
u
(
0
)
+
W
n
R
^
n
+
1
-
1
z
^
(
n
+
1
)
=
u
(
0
)
+
P
^
n
z
^
(
n
+
1
)
=
u
(
0
)
+
P
n
-
1
z
(
n
)
+
ζ
^
n
p
^
(
n
)
=
u
(
n
)
+
ζ
^
n
p
^
(
n
)
.
We note that x
^
(
n
) is the estimated solution vector satisfying the Galerkin condition, while
(59)
u
(
n
)
=
u
(
0
)
+
W
n
-
1
y
(
n
)
with y
(
n
) minimizing |
|
T
~
n
+
1
y
(
n
)
-
σ
0
e
1
|
|.


## 5. Generalized SYMMLQ and SYMMQR Methods

Now, we generalize the SYMMLQ and SYMMQR methods.Theorem 3.
Suppose thatE is an n
×
n symmetric positive definite (SPD) matrix and E
A is an n
×
n symmetric matrix. One can generate orthonormal vectors w
(
0
)
,
w
(
1
)
,
…
w
(
n
-
2
)
,
w
(
n
-
1
) using this short-term recurrence
(60)
w
~
(
j
+
1
)
≡
A
w
(
j
)
-
α
j
w
(
j
)
-
β
j
w
(
j
-
1
)
(
0
≤
j
≤
n
-
2
)
,
w
~
(
j
+
1
)

=

(
1
γ
j
+
1
)
w
~
(
j
+
1
)
,
w
h
e
r
e
γ
j
+
1

=

|
〈
E
w
~
(
j
+
1
)
,
w
~
(
j
+
1
)
〉
|
,
where
(61)
α
j
=
〈
E
A
w
(
j
)
,
w
(
j
)
〉
β
j
=
〈
E
A
w
(
j
)
,
w
(
j
-
1
)
〉
.
Then the following properties hold, for (0
≤
i, j
≤
n
-
1):
(62)
〈
w
(
j
)
,
E
w
(
j
)
〉
=
δ
i
j
,
β
j
=
γ
j
.Proof.
We obtain(63)
β
j
=
〈
E
A
w
(
j
)
,
w
(
j
-
1
)
〉
(
by
(
61
)
)
=
〈
w
(
j
)
,
E
A
w
(
j
-
1
)
〉
(
(
E
A
)
T
=
E
A
)
=
〈
w
(
j
)
,
E
(
γ
j
w
(
j
)
+
α
j
-
1
w
(
j
-
1
)
+
β
j
-
1
w
(
j
-
2
)
)
〉
nnnnnnnnnnnnnnnnnnnnnnnnn
(
by
(
60
)
)
=
〈
w
(
j
)
,
γ
j
E
w
(
j
)
〉
=
γ
j
.
Since 〈
w
(
i
)
,
E
w
(
j
)
〉
=
δ
i
j.As before, we let(64)
W
n
=
[
w
(
1
)
,
w
(
2
)
,
…
,
w
(
n
-
1
)
,
w
(
n
)
]
n
×
n
,
e
n
T
=
[
0,0
,
0
,
…
,
0,0
,
1
]
1
×
n
.
Moreover, we have
(65)
A
W
n
=
W
n
T
n
+
γ
n
+
1
w
(
n
+
1
)
e
n
T
,
where
(66)
T
n
=
[
α
1
β
2
σ
2
α
2
β
3
σ
3
α
3
β
4
⋱
⋱
⋱
σ
n
-
2
α
n
-
2
β
n
-
1
σ
n
-
1
α
n
-
1
β
n
σ
n
α
n
]
n
×
n
.As before, we let(67)
x
(
n
)
=
x
(
0
)
+
W
n
y
(
n
)
,
r
(
n
)
=
r
(
0
)
-
A
W
n
y
(
n
)
.
Imposing the Galerkin condition again, we have
(68)
(
E
r
(
n
)
)
T
W
n
=
0
,
(
E
W
n
)
T
r
(
n
)
=
0
,
(
E
W
n
)
T
r
(
0
)
=
(
E
W
n
)
T
A
W
n
y
(
n
)
(
by
(
67
)
)
.
We obtain
(69)
T
n
y
(
n
)
=
β
1
e
1
where
β
1
=
|
〈
E
r
(
0
)
,
r
(
0
)
〉
|
because
(70)
(
E
W
n
)
T
A
W
n
=
T
n
,
(
E
W
n
)
T
W
n
=
I
,
(
E
W
n
)
T
r
(
0
)
=
β
1
e
1
.
Since T
n is symmetric, we can apply the same techniques as in the SYMMLQ method. Also, if E
=
I, the method reduces to the SYMMQR method.


## 6. Modified SYMMLQ Method

Next, we outline the modified SYMMLQ method.Theorem 4.
Suppose thatE is an n
×
n symmetric (not necessary positive definite) matrix and E
A is an n
×
n symmetric matrix. One can generate orthonormal vectors w
(
0
)
,
w
(
1
)
,
…
w
(
n
-
2
)
,
w
(
n
-
1
) using this short-term recurrence
(71)
w
~
(
j
+
1
)
≡
A
w
(
j
)
-
α
j
w
(
j
)
-
β
j
w
(
j
-
1
)
,
(
0
≤
j
≤
n
-
2
)
,
w
(
j
+
1
)

=

(
1
γ
j
+
1
)
w
~
(
j
+
1
)
,
w
h
e
r
e
γ
i
+
1

=

|
E
w
~
(
j
+
1
)
,
w
~
(
j
+
1
)
〉
|
,
where
(72)
α
j
=
〈
E
A
w
(
j
)
,
w
(
j
)
〉
〈
E
w
(
j
)
,
w
(
j
)
〉
,
β
j
=
〈
E
A
w
(
j
)
,
w
(
j
-
1
)
〉
〈
E
w
(
j
-
1
)
,
w
(
j
-
1
)
〉
.
Then the following properties hold, for (
0
≤
j
≤
n
-
2
):
(73)
d
j
+
1
=
〈
E
w
(
j
+
1
)
,
w
(
j
+
1
)
〉
=
{
1
,
i
f
〈
E
w
~
(
j
+
1
)
,
w
~
(
j
+
1
)
〉
>
0
,
-
1
,
i
f
〈
E
w
~
(
j
+
1
)
,
w
~
(
j
+
1
)
〉
<
0
,
and, for (0
≤
i, j
≤
n
-
1),
(74)
〈
E
w
(
i
)
,
w
(
j
)
〉
=
δ
i
j
.From Theorem4, in matrix form, we obtain
(75)
A
W
n
=
W
n
T
n
+
γ
n
+
1
w
(
n
+
1
)
e
n
T
,
where
(76)
W
n
=
[
w
(
1
)
,
w
(
2
)
,
…
,
w
(
n
-
2
)
,
w
(
n
-
1
)
,
w
(
n
)
]
n
×
n
,
e
n
T
=
[
0,0
,
0
,
…
,
0,0
,
1
]
1
×
n
,
T
n
=
[
α
1
β
2
γ
2
α
2
β
3
γ
3
α
3
β
4
⋱
⋱
⋱
γ
n
-
2
α
n
β
n
-
1
γ
n
-
1
α
n
-
1
β
n
γ
n
α
n
]
n
×
n
.
Moreover, from Theorem 4, we obtain
(77)
(
E
W
n
)
T
W
n
=
[
d
1
d
2
d
3
⋱
d
n
-
2
d
n
-
1
d
n
]
n
×
n
≡
D
n
.
Then, we have
(78)
(
E
W
n
)
T
A
W
n
=
(
E
W
n
)
T
W
n
T
n
+
(
E
W
n
)
T
(
γ
n
+
1
w
(
n
+
1
)
e
n
T
)
=
D
n
T
n
+
0
.
Here the second term on the right-hand side is the zero matrix!In addition, we have(79)
x
(
n
)
=
x
(
0
)
+
W
n
y
(
n
)
,
r
(
n
)
=
r
(
0
)
-
A
W
n
y
(
n
)
.
Imposing the Galerkin condition, r
(
n
)
⊥
𝒦
n
(
r
(
0
)
,
A
), as we did before, we obtain
(80)
(
E
r
(
n
)
)
T
W
n
=
0
,
(
E
W
n
)
T
r
(
n
)
=
0
,
(
E
W
n
)
T
r
(
0
)
=
(
E
W
n
)
T
A
W
n
y
(
n
)
(
by
(
79
)
)
,
(
E
W
n
)
T
r
(
0
)
=
D
n
T
n
y
(
n
)
(
by
(
78
)
)
.
In other words, we use
(81)
(
E
W
n
)
T
A
W
n
=
D
n
T
n
,
(
E
W
n
)
T
W
n
=
D
n
.
We obtain
(82)
(
E
W
n
)
T
r
(
0
)
=
β
1
e
1
because
(83)
D
n
T
n
y
(
n
)
=
β
1
e
1
,
where
β
1
=
|
〈
E
r
(
0
)
,
r
(
0
)
〉
|
.
Here(84)
D
n
T
n
=
[
d
1
α
1
d
1
β
2
d
2
γ
2
d
2
α
2
d
2
β
3
d
3
γ
3
d
3
α
3
d
3
α
4
⋱
⋱
⋱
d
n
-
2
γ
n
-
2
d
n
-
2
α
n
-
2
d
n
-
2
β
n
-
1
d
n
-
1
γ
n
-
1
d
n
-
1
α
n
-
1
d
n
-
1
β
n
d
n
γ
n
d
n
α
n
]
n
×
n
.We note that D
n
T
n is symmetric, for (
1
≤
i
≤
n
-
1
):
(85)
d
i
β
i
+
1
=
d
i
〈
E
A
w
(
i
+
1
)
,
w
(
i
)
〉
〈
E
w
(
i
)
,
w
(
i
)
〉
(
by
(
72
)
)
=
d
i
〈
w
(
i
+
1
)
,
E
A
w
(
i
)
〉
〈
E
w
(
i
)
,
w
(
i
)
〉
(
(
E
A
)
T
=
E
A
)
=
d
i
〈
w
(
i
+
1
)
,
E
(
γ
i
+
1
w
(
i
+
1
)
+
α
i
w
(
i
)
+
β
i
w
(
i
-
1
)
〉
〈
E
w
(
i
)
,
w
(
i
)
〉
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
l
(
by
(
71
)
)
=
d
i
+
1
γ
j
+
1
.


## 7. Lanczos/MSYMMLQ Method

Next, we consider this2
n
×
2
n double linear system:
(86)
[
A
0
0
A
T
]
[
x
x
^
]
=
[
b
b
^
]
.
We obtain the block symmetric matrices 𝒜, ℰ, and ℰ
𝒜, where
(87)
𝒜
=
[
A
0
0
A
T
]
,
ℰ
=
[
0
I
I
0
]
,
ℰ
𝒜
=
[
0
A
T
A
0
]
.For example, the modified SYMMLQ method and the modified SYMMQR method can be applied to the double linear system (86). This leads us to the LAN/MSYMMLQ method and the LAN/MSYMMQR method. The pseudocodes for these methods are given in the following sections. For additional details, see Li [3]. See the books by Golub and Van Loan [4] and Saad [5], as well as the papers by Lanczos [6] and Kincaid et al. [7], among others.


## 8. MSYMMLQ Pseudocode

(88)
r
(
0
)
=
b
-
A
x
(
0
)
,
β
1
=
|
E
r
(
0
)
,
r
(
0
)
|
,
w
(
1
)
=
(
1
β
1
)
r
(
0
)
,
d
1
=
{
1
,
if
〈
E
r
(
0
)
,
E
r
(
0
)
〉
>
0
,
-
1
,
if
〈
E
r
(
0
)
,
E
r
(
0
)
〉
<
0
,
ε
1
=
ε
2
=
0
,
s
0
=
0
,
c
0
=
-
1
,
α
1
=
(
1
d
1
)
〈
E
A
w
(
1
)
,
w
(
1
)
〉
,
w
~
(
2
)
=
A
w
(
1
)
-
α
1
w
(
1
)
,
β
2
=
|
〈
E
w
~
(
2
)
,
w
~
(
2
)
〉
|
,
w
(
2
)
=
(
1
β
2
)
w
~
(
2
)
,
d
2
=
{
1
,
if
〈
E
w
~
(
2
)
,
w
~
(
2
)
〉
>
0
,
-
1
,
if
〈
E
w
~
(
2
)
,
w
~
(
2
)
〉
<
0
,
ζ
^
1
=
β
1
α
1
,
v
^
(
1
)
=
w
(
1
)
,
x
^
(
1
)
=
x
(
0
)
+
ζ
^
1
v
(
1
)
,
c
1
=
d
1
α
1
α
1
2
+
β
2
2
,
s
1
=
d
1
β
2
α
1
2
+
β
2
2
,
γ
1
=
(
d
1
α
1
)
c
1
+
(
d
1
β
2
)
s
1
,
ζ
1
=
(
β
1
γ
1
)
,
v
(
1
)
=
c
i
v
^
(
1
)
+
s
1
w
(
2
)
,
x
(
1
)
=
x
(
0
)
+
ζ
1
v
(
1
)
,
(89)
f
o
r
i
=
2,3
,
…
,
N
,
α
i
=
(
1
d
i
)
〈
E
A
w
(
i
)
,
w
(
i
)
〉
,
β
i
=
(
1
d
i
-
1
)
〈
E
A
w
(
i
)
,
w
(
i
-
1
)
〉
,
w
~
(
i
+
1
)
=
A
w
(
i
)
-
α
i
w
(
i
)
-
β
i
w
(
i
-
1
)
,
β
i
+
1
=
|
〈
E
w
~
(
i
+
1
)
,
w
~
(
i
+
1
)
〉
|
,
w
(
i
+
1
)
=
(
1
β
i
+
1
)
w
~
(
i
+
1
)
,
d
i
+
1
=
{
1
,
if
〈
E
w
~
(
i
+
1
)
,
w
~
(
i
+
1
)
〉
>
0
,
-
1
,
if
〈
E
w
~
(
i
+
1
)
,
w
~
(
i
+
1
)
〉
<
0
,
ε
i
=
(
d
i
-
1
β
i
)
s
i
-
2
,
δ
i
=
-
(
d
i
-
1
β
i
)
c
i
-
2
,
h
h
=
δ
i
,
δ
i
=
(
h
h
)
c
i
-
1
+
(
d
i
α
i
)
s
i
-
1
,
γ
^
i
=
(
h
h
)
s
i
-
1
-
(
d
i
α
i
)
c
i
-
1
,
ζ
^
i
=
(
1
γ
^
i
)
(
-
ε
i
ζ
i
-
2
-
δ
i
ζ
i
-
1
)
,
v
^
(
i
)
=
s
i
-
1
v
^
(
i
-
1
)
-
c
i
-
1
w
(
i
)
,
x
^
(
i
)
=
x
(
i
-
1
)
+
ζ
^
i
v
^
(
i
)
,
γ
i
=
α
i
2
+
β
i
+
1
2
,
c
i
=
d
i
(
α
i
γ
i
)
,
s
i
=
d
i
(
β
i
+
1
γ
i
)
,
ζ
i
=
(
1
γ
1
)
(
-
ε
i
ζ
i
-
2
-
δ
i
ζ
i
-
1
)
,
v
(
i
)
=
c
1
v
^
(
i
)
+
s
i
w
(
i
+
1
)
,
x
(
i
)
=
x
(
i
-
1
)
+
ζ
i
v
(
i
)
e
n
d
f
o
r


## 9. MSYMMQR Pseudocode

(90)
r
(
0
)
=
b
-
A
x
(
0
)
,
β
1
=
|
E
r
(
0
)
,
r
(
0
)
|
,
w
(
1
)
=
(
1
β
1
)
r
(
0
)
,
d
1
=
{
1
,
if
〈
E
r
(
0
)
,
E
r
(
0
)
〉
>
0
,
-
1
,
if
〈
E
r
(
0
)
,
E
r
(
0
)
〉
<
0
,
ε
1
=
ε
2
=
0
,
s
0
=
0
,
c
0
=
1
,
α
1
=
(
1
d
1
)
〈
E
A
w
(
1
)
,
w
(
1
)
〉
,
w
~
(
2
)
=
A
w
(
1
)
-
α
1
w
(
1
)
,
β
2
=
|
〈
E
w
~
(
2
)
,
w
~
(
2
)
〉
|
,
w
(
2
)
=
(
1
β
2
)
w
~
(
2
)
,
d
2
=
{
1
,
if
〈
E
w
~
(
2
)
,
w
~
(
2
)
〉
>
0
,
-
1
,
if
〈
E
w
~
(
2
)
,
w
~
(
2
)
〉
<
0
,
ζ
^
1
=
β
1
d
1
,
p
~
(
1
)
=
(
1
α
1
d
1
)
w
(
1
)
,
x
^
(
1
)
=
x
(
0
)
+
ζ
^
1
p
^
(
1
)
,
c
1
=
d
1
α
1
α
1
2
+
β
2
2
,
s
1
=
-
d
1
β
2
α
1
2
+
β
2
2
,
γ
1
=
(
d
1
α
1
)
c
1
-
(
d
1
β
2
)
s
1
,
ζ
1
=
c
1
ζ
^
1
,
p
(
1
)
=
(
1
γ
1
)
w
(
1
)
,
x
(
1
)
=
x
(
0
)
+
ζ
1
p
(
1
)
,
(91)
f
o
r
i
=
2,3
,
…
,
N
,
α
i
=
(
1
d
i
)
〈
E
A
w
(
i
)
,
w
(
i
)
〉
,
β
i
=
(
1
d
i
-
1
)
〈
E
A
w
(
i
)
,
w
(
i
-
1
)
〉
,
w
~
(
i
+
1
)
=
A
w
(
i
)
-
α
i
w
(
i
)
-
β
i
w
(
i
-
1
)
,
β
i
+
1
=
|
〈
E
w
~
(
i
+
1
)
,
w
~
(
i
+
1
)
〉
|
,
w
(
i
+
1
)
=
(
1
β
i
+
1
)
w
~
(
i
+
1
)
,
d
i
+
1
=
{
1
,
if
〈
E
w
~
(
i
+
1
)
,
w
~
(
i
+
1
)
〉
>
0
,
-
1
,
if
〈
E
w
~
(
i
+
1
)
,
w
~
i
+
1
〉
<
0
,
ε
i
=
-
(
d
i
-
1
β
i
)
s
i
-
2
,
δ
i
=
-
(
d
i
-
1
β
i
)
c
i
-
2
,
h
h
=
δ
i
,
δ
i
=
(
h
h
)
c
i
-
1
-
(
d
i
α
i
)
s
i
-
1
,
γ
^
i
=
(
h
h
)
s
i
-
1
+
(
d
i
α
i
)
c
i
-
1
,
ζ
^
i
=
s
i
-
1
ζ
^
i
-
1
,
p
^
(
i
)
=
(
1
γ
~
i
)
[
w
(
i
)
-
ε
i
p
(
i
-
2
)
-
δ
i
p
(
i
-
1
)
]
,
x
^
(
i
)
=
x
(
i
-
1
)
+
ζ
^
i
p
^
(
i
)
,
γ
i
=
ζ
^
i
2
+
β
i
+
1
2
,
c
i
=
γ
^
i
γ
i
,
s
i
=
-
(
1
γ
i
)
(
d
i
β
i
+
1
)
,
ζ
i
=
c
i
-
1
ζ
^
i
,
p
(
i
)
=
(
1
γ
i
)
[
w
(
i
)
-
ε
i
p
(
i
-
2
)
-
δ
i
p
(
i
-
1
)
]
,
x
(
i
)
=
x
(
i
-
1
)
+
ζ
i
p
(
i
)
e
n
d
f
o
r


## 10. LAN/MSYMMLQ Pseudocode

(92)
r
(
0
)
=
b
-
A
x
(
0
)
,
r
^
(
0
)
=
b
^
-
A
T
x
(
0
)
,
β
1
=
|
2

〈
r
^
(
0
)
,
r
(
0
)
〉
|
,
w
(
1
)
=
(
1
β
1
)
r
(
0
)
,
w
^
(
1
)
=
(
1
β
1
)
r
^
(
0
)
,
d
1
=
{
1
,
if
2

〈
r
^
(
0
)
,
r
(
0
)
〉

>

0
,
-
1
,
if
2

〈
r
^
(
0
)
,
r
(
0
)
〉

<

0
,
ε
1
=
ε
2
=
0
,
s
0
=
0
,
c
0
=
-
1
,
α
1
=
(
2
d
1
)
〈
w
^
(
1
)
,
w
(
1
)
〉
,
w
~
(
2
)
=
A
w
(
1
)
-
α
1
w
(
1
)
,
w
^
~
(
2
)
=
A
T
w
(
1
)
-
α
1
w
^
(
1
)
,
β
2
=
|
2

〈
w
^
~
(
2
)
,
w
~
(
2
)
〉
|
,
w
2
=
(
1
β
2
)
w
~
(
2
)
,
w
^
(
2
)
=
(
1
β
2
)
w
^
~
(
2
)
,
d
2
=
{
1
,
if
2

〈
w
^
~
(
2
)
,
w
~
(
2
)
〉
>
0
,
-
1
,
if
2

〈
w
^
~
(
2
)
,
w
~
(
2
)
〉
<
0
,
ζ
^
1
=
β
1
α
1
,
v
^
(
1
)
=
w
(
1
)
,
x
^
(
1
)
=
x
(
0
)
+
ζ
^
1
v
^
(
1
)
,
c
1
=
d
1
α
1
α
1
2
+
β
2
2
,
s
1
=
d
1
β
2
α
1
2
+
β
2
2
,
γ
1
=
(
d
1
α
1
)
c
1
+
(
d
1
β
2
)
s
1
,
ζ
1
=
β
1
γ
1
,
v
(
1
)
=
c
1
v
^
(
1
)
+
s
1
w
(
2
)
,
x
(
1
)
=
x
(
0
)
+
ζ
1
v
(
1
)
,
(93)
f
o
r
i
=
2,3
,
…
,
N
,
α
i
=
(
2
d
i
)
〈
A
w
(
i
)
,
w
^
(
i
)
〉
,
β
i
=
(
1
d
i
-
1
)
[
〈
A
w
(
i
)
,
w
^
(
i
-
1
)
〉
+
〈
w
^
(
i
)
,
A
w
(
i
-
1
)
〉
]
,
w
~
(
i
+
1
)
=
A
w
(
i
)
-
α
i
w
(
i
)
-
β
i
w
(
i
-
1
)
,
w
^
~
(
i
+
1
)
=
A
T
w
^
(
i
)
-
α
i
w
^
(
i
)
-
β
i
w
^
(
i
-
1
)
,
β
i
+
1
=
|
2

〈
w
~
(
i
+
1
)
,
w
~
(
i
+
1
)
〉
|
,
w
(
i
+
1
)
=
(
1
β
i
+
1
)
w
~
(
i
+
1
)
,
w
^
(
i
+
1
)
=
(
1
β
i
+
1
)
w
^
~
(
i
+
1
)
,
d
i
+
1
=
{
1
,
if
2

〈
w
~
(
i
+
1
)
,
w
~
(
i
+
1
)
〉

>

0
,
-
1
,
if
2

〈
w
~
(
i
+
1
)
,
w
~
i
+
1
〉

<

0
,
ε
i
=
d
i
-
1
β
i
s
i
-
2
,
δ
i
=
-
(
d
i
-
1
β
i
)
c
i
-
2
,
h
h
=
δ
i
,
δ
i
=
(
h
h
)
c
i
-
1
+
(
d
i
α
i
)
s
i
-
1
,
γ
^
i
=
(
h
h
)
s
i
-
1
-
(
d
i
α
i
)
c
i
-
1
,
ζ
^
i
=
(
1
γ
^
i
)
(
-
ε
i
ζ
i
-
2
-
δ
i
ζ
i
-
1
)
,
v
^
(
i
)
=
s
i
-
1
v
^
(
i
-
1
)
-
c
i
-
1
w
(
i
)
,
x
^
(
i
)
=
x
(
i
-
1
)
+
ζ
^
i
v
^
(
i
)
,
γ
i
=
α
i
2
+
β
i
+
1
2
,
c
i
=
d
i
(
α
i
γ
i
)
,
s
i
=
d
i
(
β
i
+
1
γ
i
)
,
ζ
i
=
(
1
γ
i
)
(
-
ε
i
ζ
i
-
2
-
δ
i
ζ
i
-
1
)
,
v
(
i
)
=
c
i
v
^
(
i
)
+
s
i
w
(
i
+
1
)
,
x
(
i
)
=
x
(
i
-
1
)
+
ζ
i
v
(
i
)
e
n
d
f
o
r


## 11. LAN/MSYMMQR Pseudocode

(94)
r
(
0
)
=
b
-
A
x
(
0
)
,
r
^
(
0
)
=
b
^
-
A
T
x
(
0
)
,
β
1
=
|
2
〈
r
^
(
0
)
,
r
(
0
)
〉
|
w
(
1
)
=
(
1
β
1
)
r
(
0
)
,
w
^
(
1
)
=
(
1
β
1
)
r
^
(
0
)
,
d
1
=
{
1
,
if
2

〈
r
^
(
0
)
,
r
(
0
)
〉

>

0
,
-
1
,
if
2

〈
r
^
(
0
)
,
r
(
0
)
〉

<

0
,
ε
1
=
ε
2
=
0
,
s
0
=
0
,
c
0
=
1
,
α
1
=
(
2
d
1
)

〈
w
^
(
1
)
,
w
(
1
)
〉
,
w
~
(
2
)
=
A
w
(
1
)
-
α
1
w
(
1
)
,
w
^
~
(
2
)
=
A
T
w
(
1
)
-
α
1
w
^
(
2
)
,
β
2
=
|
2
〈
w
^
~
(
2
)
,
w
~
(
2
)
〉
|
,
w
(
2
)
=
(
1
β
2
)
w
~
(
2
)
,
w
^
(
2
)
=
(
1
β
2
)
w
^
~
(
2
)
,
d
2
=
{
1
,
if
2

〈
w
^
~
(
2
)
,
w
~
(
2
)
〉
>
0
,
-
1
,
if
2

〈
w
^
~
(
2
)
,
w
~
(
2
)
〉
<
0
,
p
^
(
1
)
=
(
1
α
i
d
i
)
w
(
1
)
,
ζ
^
1
=
β
1
α
1
,
x
^
(
1
)
=
x
(
0
)
+
ζ
^
1
p
^
(
1
)
,
c
1
=
d
1
α
1
α
1
2
+
β
2
2
,
s
1
=
-
d
1
β
2
α
1
2
+
β
2
2
,
γ
1
=
(
d
1
α
1
)
c
1
-
(
d
1
β
2
)
s
1
,
ζ
1
=
β
1
γ
^
1
,
p
(
1
)
=
(
1
γ
1
)
w
(
1
)
,
x
(
1
)
=
x
(
0
)
+
ζ
1
p
(
1
)
,
(95)
f
o
r
i
=
2,3
,
…
,
N
,
α
i
=
(
2
d
i
)
〈
A
w
(
i
)
,
w
^
(
i
)
〉
,
β
i
=
(
1
d
i
-
1
)
[
〈
A
w
(
i
)
,
w
^
(
i
-
1
)
〉
+
〈
w
^
(
i
)
,
A
w
(
i
-
1
)
〉
]
,
w
~
(
i
+
1
)
=
A
w
(
i
)
-
α
i
w
(
i
)
-
β
i
w
(
i
-
1
)
,
w
^
~
(
i
+
1
)
=
A
T
w
^
(
i
)
-
α
i
w
^
(
i
)
-
β
i
w
^
(
i
-
1
)
,
β
i
+
1
=
|
2

〈
w
~
(
i
+
1
)
,
w
~
(
i
+
1
)
〉
|
,
w
(
i
+
1
)
=
(
1
β
i
+
1
)
w
~
(
i
+
1
)
,
w
^
(
i
+
1
)
=
(
1
β
i
+
1
)
w
^
~
(
i
+
1
)
,
d
i
+
1
=
{
1
,
if
2

〈
w
~
(
i
+
1
)
,
w
~
(
i
+
1
)
〉
>
0
,
-
1
,
if
2

〈
w
~
(
i
+
1
)
,
w
~
(
i
+
1
)
〉

<

0
,
ε
i
=
(
d
i
-
1
β
i
)
s
i
-
2
,
δ
i
=
(
d
i
-
1
β
i
)
c
i
-
2
,
h
h
=
δ
i
,
δ
i
=
(
h
h
)
c
i
-
1
-
d
i
(
α
i
s
i
-
1
)
,
γ
^
i
=
(
h
h
)
s
c
i
-
1
+
d
i
(
α
i
c
i
-
1
)
,
ζ
^
i
=
ζ
~
^
i
-
1
s
i
-
1
,
p
^
(
i
)
=
(
1
γ
^
i
)
[
w
(
i
)
-
ε
i
p
(
i
)
-
δ
i
p
(
i
-
1
)
]
,
x
^
(
i
)
=
x
(
i
-
1
)
+
ζ
^
i
p
^
(
i
)
,
γ
i
=
γ
^
2
+
β
i
+
1
2
,
c
i
=
(
γ
^
i
γ
i
)
,
s
i
=
-
d
i
(
β
i
+
1
γ
i
)
,
ζ
i
=
c
i
ζ
^
i
,
p
(
i
)
=
(
1
γ
i
)
[
w
(
i
)
-
ε
i
p
(
i
-
2
)
-
δ
i
p
(
i
-
1
)
]
,
x
(
i
)
=
x
(
i
-
1
)
+
ζ
i
p
(
i
)
e
n
d
f
o
r



---

*Source: 237808-2014-04-03.xml*