,

A, …, Z  matrices
a, …, z  vectors
α, β, …, ω  scalars
A^{T}  matrix transpose
A^{H}  conjugate transpose (Hermitian) of A
A^{−1}  matrix inverse
A^{−T}  the inverse of A^{T}
a_{i,j}  matrix element
a_{.,j}  jth matrix column
A_{i,j}  matrix subblock
a_{i}  vector element
u_{x}, u_{xx}  first, second derivative with respect to x
(x, y), x^{T} y  vector dot product (inner product)
x_{j}^{(i)}  jth component of vector x in the ith iteration
diag(A)  diagonal of matrix A
diag(α, β …)  diagonal matrix constructed from scalars α, β …
span(a, b …)  spanning space of vectors a, b …
R  set of real numbers
R^{n}  real n-space
||x||_{2}  2-norm
||x||_{p}  p-norm
||x||_{A}  the “A-norm”, defined as (Ax, x)^{1/2}
λ_{max}(A), λ_{min}(A) eigenvalues of A with maximum (resp. minimum) modulus
σ_{max}(A), σ_{min}(A) largest and smallest singular values of A
κ_{2}(A)  spectral condition number of matrix A
L  linear operator
\overline{α}  complex conjugate of the scalar α
max{S}  maximum value in set S
min{S}  minimum value in set S
∑  summation
O(·)  “big-oh” asymptotic bound

Conventions Used in this Book

D  diagonal matrix
L  lower triangular matrix
U  upper triangular matrix
Q  orthogonal matrix
M  preconditioner
I, I^{n×n}  n × n identity matrix
x̂  typically, the exact solution to Ax = b
h  discretization mesh width

A y = A · x z = A^{T} · x y z x Ax = b Ay − b ω (AA^{T})y = b x = A^{T} y (A^{T} A)x = b̃ b̃ = A^{T} b AA^{T} A^{T} A A^{T} x^{(k)} = Bx^{(k−1)} + c (2.1)

B c k n i

∑_{j=1}^{n} a_{i,j}x_{j} = b_{i},

x_{i}

x_{i} = (b_{i} − ∑_{j≠i} a_{i,j}x_{j})/a_{i,i}. (2.2)

x_{i}^{(k)} = (b_{i} − ∑_{j≠i} a_{i,j}x_{j}^{(k−1)})/a_{i,i}, (2.3)

x^{(k)} = D^{−1}(L + U)x^{(k−1)} + D^{−1}b, (2.4)

D −L −U

¯ x

x^{(k−1)} x^{(k)} Choose an initial guess x^{(0)} to the solution x. for k = 1, 2, …

for i = 1, 2, …, n x̄_{i} = 0 for j = 1, 2, …, i − 1, i + 1, …, n

x̄_{i} = x̄_{i} + a_{i,j}x_{j}^{(k−1)}

end

x̄_{i} = (b_{i} − x̄_{i})/a_{i,i} end x^{(k)} = x̄ check convergence; continue if necessary

end

Figure 2.1: The Jacobi Method

Lu = −u_{xx} = f  on (0, 1),  u(0) = u_{0}, u(1) = u_{1},

Lu(x_{i}) = 2u(x_{i}) − u(x_{i−1}) − u(x_{i+1}) = f(x_{i}) for x_{i} = i/N, i = 1 … N − 1.

L L n = 1 … N − 1 u_{n}(x) = sin nπx λ = 2 sin^{2} nπ/(2N) λ(B) = 1 − 1/2λ(L) = 1 − sin^{2} nπ/(2N) u_{n} 1 ≈ 1 − 5/N^{2} u(x) = sin πx max{|λ(A)|} 1 − O(h^{2}) h h = N^{−d} N d

x_{i}^{(k)} = (b_{i} − ∑_{j<i} a_{i,j}x_{j}^{(k)} − ∑_{j>i} a_{i,j}x_{j}^{(k−1)})/a_{i,i}. (2.5)

x^{(k)} x^{(k)} = (D − L)^{−1}(Ux^{(k−1)} + b). (2.6)

Choose an initial guess x^{(0)} to the solution x. for k = 1, 2, …

for i = 1, 2, …, n σ = 0 for j = 1, 2, …, i − 1

σ = σ + a_{i,j}x_{j}^{(k)} end for j = i + 1, …, n

σ = σ + a_{i,j}x_{j}^{(k−1)} end x_{i}^{(k)} = (b_{i} − σ)/a_{i,i}

end check convergence; continue if necessary

end

Figure 2.2: The Gauss-Seidel Method x_{i}^{(k)} = ωx̄_{i}^{(k)} + (1 − ω)x_{i}^{(k−1)}

¯ x

x^{(k)} = (D − ωL)^{−1}(ωU + (1 − ω)D)x^{(k−1)} + ω(D − ωL)^{−1}b. (2.7)

Choose an initial guess x^{(0)} to the solution x. for k = 1, 2, …

for i = 1, 2, …, n σ = 0 for j = 1, 2, …, i − 1

σ = σ + a_{i,j}x_{j}^{(k)} end for j = i + 1, …, n

σ = σ + a_{i,j}x_{j}^{(k−1)} end σ = (b_{i} − σ)/a_{i,i} x_{i}^{(k)} = x_{i}^{(k−1)} + ω(σ − x_{i}^{(k−1)})

end check convergence; continue if necessary

end

Figure 2.3: The SOR Method ω = 1 (0, 2)

0 < ω < 1 ω ∈ (0, 2) ω = 2 − O(h)

M a_{i,j} = a_{j,i} j x^{T} Ax > 0 a_{ii} > ∑_{j≠i} |a_{ij}| a_{i,j} ≤ 0 i ≠ j (A^{−1})_{i,j} ≥ 0 ρ

ω_{opt} = \frac{2}{1 + \sqrt{1 − ρ^{2}}} . (2.8)

x^{(k)} = B_{1}B_{2}x^{(k−1)} + ω(2 − ω)(D − ωU)^{−1}D(D − ωL)^{−1}b, (2.9)

B_{1} = (D − ωU)^{−1}(ωL + (1 − ω)D), B_{2} = (D − ωL)^{−1}(ωU + (1 − ω)D).

B_{2} B_{1} U Choose an initial guess x^{(0)} to the solution x. for k = 1, 2, …

for i = 1, 2, …, n σ = 0 for j = 1, 2, …, i − 1

σ = σ + a_{i,j}x_{j}^{(k− \frac{1}{2} )} end for j = i + 1, …, n

σ = σ + a_{i,j}x_{j}^{(k−1)} end σ = (b_{i} − σ)/a_{i,i} x_{i}^{(k− \frac{1}{2} )} = x_{i}^{(k−1)} + ω(σ − x_{i}^{(k−1)})

end for i = n, n − 1, …, 1

σ = 0 for j = 1, 2, …, i − 1

σ = σ + a_{i,j}x_{j}^{(k− \frac{1}{2} )} end for j = i + 1, …, n

σ = σ + a_{i,j}x_{j}^{(k)} end x_{i}^{(k)} = x_{i}^{(k− \frac{1}{2} )} + ω(σ − x_{i}^{(k− \frac{1}{2} )})

end check convergence; continue if necessary

end

Figure 2.4: The SSOR Method x^{(i)} (α_{i}) p^{(i)} x^{(i)} = x^{(i−1)} + α_{i}p^{(i)}. r^{(i)} = b − Ax^{(i)}

r^{(i)} = r^{(i−1)} − αq^{(i)}  where  q^{(i)} = Ap^{(i)}.  (2.10)

α = α_{i} = r^{(i)^{T}} r^{(i)}/p^{(i)^{T}} Ap^{(i)} r^{(i)^{T}} A^{−1}r^{(i)} α p^{(i)} = r^{(i)} + β_{i−1}p^{(i−1)}, (2.11) β_{i} = r^{(i)^{T}} r^{(i)}/r^{(i−1)^{T}} r^{(i−1)} p^{(i)} Ap^{(i−1)} r^{(i)} r^{(i−1)} β_{i} p^{(i)} r^{(i)} Ap^{(j)} r^{(j)} M = I z^{(i−1)} r^{(i−1)} z^{(0)} r^{(0)} Compute r^{(0)} = b − Ax^{(0)} for some initial guess x^{(0)} for i = 1, 2, …

solve Mz^{(i−1)} = r^{(i−1)} ρ_{i−1} = r^{(i−1)^{T}} z^{(i−1)} if i = 1

p^{(1)} = z^{(0)} else

β_{i−1} = ρ_{i−1}/ρ_{i−2}

p^{(i)} = z^{(i−1)} + β_{i−1}p^{(i−1)} endif q^{(i)} = Ap^{(i)} α_{i} = ρ_{i−1}/p^{(i)^{T}} q^{(i)} x^{(i)} = x^{(i−1)} + α_{i}p^{(i)} r^{(i)} = r^{(i−1)} − α_{i}q^{(i)} check convergence; continue if necessary

end

Figure 2.5: The Preconditioned Conjugate Gradient Method x^{(i)} x^{(0)} + span{r^{(0)}, …, A^{i−1}r^{(0)}} (x^{(i)} − x̂)^{T} A(x^{(i)} − x̂)

ˆ x

r^{(i)} = b − Ax^{(i)} M^{−1} r^{(i)^{T}} M^{−1}r^{(j)} = 0 span{r^{(0)}, …, A^{i−1}r^{(0)}} {A^{i}x}_{i≥0} κ_{2} M^{−1}A λ_{max} λ_{min} κ_{2}(B) = λ_{max}(B)/λ_{min}(B))

ˆ x

∥x^{(i)} − x̂∥_{A} ≤ 2α^{i}∥x^{(0)} − x̂∥_{A} (2.12) α = (\sqrt{κ_{2}} − 1)/(\sqrt{κ_{2}} + 1)

∥y∥_{A}^{2} ≡ (y, Ay) ϵ

\sqrt{κ_{2}}

∥A∥/∥A^{−1}∥ = λ_{max}^{1/2}(A^{T} A)/λ_{min}^{1/2}(A^{T} A) κ_{2}(A) = O(h^{−2}) h^{−1} M^{−1}A LU 2 Ar^{(i)} = r^{(i+1)}t_{i+1,i} + r^{(i)}t_{i,i} + r^{(i−1)}t_{i−1,i}, AR_{i} = R_{i+1}T̄_{i},

¯ T_{i}

(i + 1) × i

T̄_{i} = \begin{matrix}← & i & → \\ ( \begin{matrix}... & ... \\ ... & ... & ... \\ & ... & ... & ... \\ & & ... & ... \\ . & . & . \end{matrix} ) . & \begin{matrix}↑ \\ i + 1 \\ ↓ \end{matrix} \end{matrix}

(·, ·)_{A} x^{(i)} ∈ {r^{(0)}, Ar^{(0)}, …, A^{i−1}r^{(0)}}, x^{(i)} = R_{i}ȳ ∥Ax^{(i)} − b∥_{2} = ∥AR_{i}ȳ − b∥_{2} = ∥R_{i+1}T̄_{i}y − b∥_{2}.

D_{i+1} ≡ diag(∥r^{(0)}∥_{2}, ∥r^{(1)}∥_{2}, ..., ∥r^{(i)}∥_{2})

R_{i+1}D_{i+1}^{−1} ∥Ax^{(i)} − b∥_{2} = ∥D_{i+1}T̄_{i}y − ∥r^{(0)}∥_{2}e^{(1)}∥_{2} (i + 1, i)

¯ T_{i}

T_{i}y = ∥r^{(0)}∥_{2}e^{(1)} T_{i} i × i

¯ T_{i}

LQ A^{T} Ax = A^{T} b

( \begin{matrix}I & A \\ A^{T} & 0 \end{matrix} ) ( \begin{matrix}r \\ x \end{matrix} ) = ( \begin{matrix}b \\ 0 \end{matrix} ) .

(p, A^{T} Ap)

(Ap, Ap) span{r^{(0)}, Ar^{(0)}, A^{2}r^{(0)}, …} w^{(i)} = Av^{(i)} for k = 1, …, i

w^{(i)} = w^{(i)} − (w^{(i)}, v^{(k)})v^{(k)} end v^{(i+1)} = w^{(i)}/∥w^{(i)}∥ {A^{k}r^{(0)}} (w^{(i)}, v^{(k)})

∥w^{(i)}∥ x^{(i)} = x^{(0)} + y_{1}v^{(1)} + ⋯ + y_{i}v^{(i)}, y_{k} ∥b − Ax^{(i)}∥ m (m) x^{(0)} is an initial guess for j = 1, 2, ....

Solve r from Mr = b − Ax^{(0)}

v^{(1)} = r/∥r∥_{2}
s := ∥r∥_{2}e_{1}
for i = 1, 2, ..., m

Solve w from Mw = Av^{(i)}

for k = 1, ..., i

h_{k,i} = (w, v^{(k)})
w = w − h_{k,i}v^{(k)}

end

h_{i+1,i} = ∥w∥_{2} v^{(i+1)} = w/h_{i+1,i} apply J_{1}, ..., J_{i−1} on (h_{1,i}, ..., h_{i+1,i}) construct J_{i}, acting on ith and (i + 1)st component of h_{.,i}, such that (i + 1)st component of J_{i}h_{.,i} is 0 s := J_{i}s if s(i + 1) is small enough then (UPDATE(x̃, i) and quit)

end

UPDATE(x̃, m) end

In this scheme UPDATE(x̃, i) replaces the following computations:

Compute y as the solution of Hy = s̃, in which the upper i × i triangular part of H has h_{i,j} as its elements (in least squares sense if H is singular), s̃ represents the first i components of s x̃ = x^{(0)} + y^{(1)}v^{(1)} + y^{(2)}v^{(2)} + ... + y^{(i)}v^{(i)} s^{(i+1)} = ∥b − Ax̃∥_{2} if x̃ is an accurate enough approximation then quit else x^{(0)} = x̃

Figure 2.6: The Preconditioned GMRES(m) Method r^{(i)} = r^{(i−1)} − α_{i}Ap^{(i)}, r̃^{(i)} = r̃^{(i−1)} − α_{i}A^{T} p̃^{(i)}, p^{(i)} = r^{(i−1)} + β_{i−1}p^{(i−1)}, p̃^{(i)} = r̃^{(i−1)} + β_{i−1}p̃^{(i−1)}.

α_{i} = \frac{r̃^{(i−1)^{T}} r^{(i−1)}}{p̃^{(i)^{T}} Ap^{(i)}} , β_{i} = \frac{r̃^{(i)^{T}} r^{(i)}}{r̃^{(i−1)^{T}} r^{(i−1)}}

r̃^{(i)^{T}} r^{(j)} = p̃^{(i)^{T}} Ap^{(j)} = 0 if i ≠ j.

Compute r^{(0)} = b − Ax^{(0)} for some initial guess x^{(0)}. Choose r̃^{(0)} (for example, r̃^{(0)} = r^{(0)}). for i = 1, 2, …

solve Mz^{(i−1)} = r^{(i−1)} solve M^{T} z̃^{(i−1)} = r̃^{(i−1)} ρ_{i−1} = z^{(i−1)^{T}} r̃^{(i−1)} if ρ_{i−1} = 0, method fails if i = 1

p^{(i)} = z^{(i−1)}

p̃^{(i)} = z̃^{(i−1)}

else

β_{i−1} = ρ_{i−1}/ρ_{i−2}

p^{(i)} = z^{(i−1)} + β_{i−1}p^{(i−1)} p̃^{(i)} = z̃^{(i−1)} + β_{i−1}p̃^{(i−1)}

endif q^{(i)} = Ap^{(i)} q̃^{(i)} = A^{T} p̃^{(i)} α_{i} = ρ_{i−1}/p̃^{(i)^{T}} q^{(i)} x^{(i)} = x^{(i−1)} + α_{i}p^{(i)} r^{(i)} = r^{(i−1)} − α_{i}q^{(i)} r̃^{(i)} = r̃^{(i−1)} − α_{i}q̃^{(i)} check convergence; continue if necessary

end

Figure 2.7: The Preconditioned BiConjugate Gradient Method z^{(i−1)^{T}} r̃^{(i−1)} ≈ 0 p̃^{(i)^{T}} q^{(i)} ≈ 0 Ap^{(k)} A^{T} p̃^{(k)} Compute r^{(0)} = b − Ax^{(0)} for some initial guess x^{(0)} ṽ^{(1)} = r^{(0)}; solve M_{1}y = ṽ^{(1)}; ρ_{1} = ∥y∥_{2} Choose w̃^{(1)} , for example w̃^{(1)} = r^{(0)} solve M_{2}z = w̃^{(1)}; ξ_{1} = ∥z∥_{2} γ_{0} = 1; η_{0} = −1 for i = 1, 2, …

if ρ_{i} = 0 or ξ_{i} = 0 method fails v^{(i)} = ṽ^{(i)}/ρ_{i}; y = y/ρ_{i} w^{(i)} = w̃^{(i)}/ξ_{i}; z = z/ξ_{i} δ_{i} = z^{T} y; if δ_{i} = 0 method fails solve M_{2}ỹ = y solve M_{1}^{T} z̃ = z if i = 1

p^{(1)} = ỹ; q^{(1)} = z̃

else

p^{(i)} = ỹ − (ξ_{i}δ_{i}/ϵ_{i−1})p^{(i−1)} q^{(i)} = z̃ − (ρ_{i}δ_{i}/ϵ_{i−1})q^{(i−1)}

endif p̃ = Ap^{(i)} ϵ_{i} = q^{(i)^{T}} p̃; if ϵ_{i} = 0 method fails β_{i} = ϵ_{i}/δ_{i}; if β_{i} = 0 method fails ṽ^{(i+1)} = p̃ − β_{i}v^{(i)} solve M_{1}y = ṽ^{(i+1)} ρ_{i+1} = ∥y∥_{2} w̃^{(i+1)} = A^{T} q^{(i)} − β_{i}w^{(i)} solve M_{2}^{T} z = w̃^{(i+1)} ξ_{i+1} = ∥z∥_{2} θ_{i} = ρ_{i+1}/(γ_{i−1}|β_{i}|); γ_{i} = 1/\sqrt{1 + θ_{i}^{2}}; if γ_{i} = 0 method fails η_{i} = −η_{i−1}ρ_{i}γ_{i}^{2}/(β_{i}γ_{i−1}^{2}) if i = 1

d^{(1)} = η_{1}p^{(1)}; s^{(1)} = η_{1}p̃ else

d^{(i)} = η_{i}p^{(i)} + (θ_{i−1}γ_{i})^{2}d^{(i−1)} s^{(i)} = η_{i}p̃ + (θ_{i−1}γ_{i})^{2}s^{(i−1)}

endif x^{(i)} = x^{(i−1)} + d^{(i)} r^{(i)} = r^{(i−1)} − s^{(i)} check convergence; continue if necessary

end

Figure 2.8: The Preconditioned Quasi Minimal Residual Method without Look-ahead

ˆ x

M = M_{1}M_{2} M_{1} = I ∥r^{(i)}∥ r^{(i)}

\sqrt{i + 1}

r^{(i)} r^{(0)} r^{(i)} = P_{i}(A)r^{(0)}. (2.13) r̃^{(i)} = P_{i}(A)r̃^{(0)} ρ_{i} = (r̃^{(i)}, r^{(i)}) = (P_{i}(A^{T})r̃^{(0)}, P_{i}(A)r^{(0)}) = (r̃^{(i)}, P_{i}^{2}(A)r^{(0)}). (2.14)

P_{i}(A) r^{(0)} r^{(i)} P_{i}^{2}(A)r^{(0)} Compute r^{(0)} = b − Ax^{(0)} for some initial guess x^{(0)} Choose r̃ (for example, r̃ = r^{(0)}) for i = 1, 2, …

ρ_{i−1} = r̃^{T} r^{(i−1)} if ρ_{i−1} = 0 method fails if i = 1

u^{(1)} = r^{(0)}

p^{(1)} = u^{(1)} else

β_{i−1} = ρ_{i−1}/ρ_{i−2}

u^{(i)} = r^{(i−1)} + β_{i−1}q^{(i−1)}

p^{(i)} = u^{(i)} + β_{i−1}(q^{(i−1)} + β_{i−1}p^{(i−1)}) endif solve Mp̂ = p^{(i)} v̂ = Ap̂ α_{i} = ρ_{i−1}/r̃^{T} v̂ q^{(i)} = u^{(i)} − α_{i}v̂ solve Mû = u^{(i)} + q^{(i)} x^{(i)} = x^{(i−1)} + α_{i}û q̂ = Aû r^{(i)} = r^{(i−1)} − α_{i}q̂ check convergence; continue if necessary

end

Figure 2.9: The Preconditioned Conjugate Gradient Squared Method r^{(0)} r^{(k)} = P_{k}(A)r^{(0)} i ↦ P_{i}^{2}(A)r^{(0)} i ↦ Q_{i}(A)P_{i}(A)r^{(0)} Q_{i} Compute r^{(0)} = b − Ax^{(0)} for some initial guess x^{(0)} Choose r̃ (for example, r̃ = r^{(0)}) for i = 1, 2, …

ρ_{i−1} = r̃^{T} r^{(i−1)} if ρ_{i−1} = 0 method fails if i = 1

p^{(i)} = r^{(i−1)} else

β_{i−1} = (ρ_{i−1}/ρ_{i−2})(α_{i−1}/ω_{i−1})

p^{(i)} = r^{(i−1)} + β_{i−1}(p^{(i−1)} − ω_{i−1}v^{(i−1)}) endif solve Mp̂ = p^{(i)} v^{(i)} = Ap̂ α_{i} = ρ_{i−1}/r̃^{T} v^{(i)} s = r^{(i−1)} − α_{i}v^{(i)} check norm of s; if small enough: set x^{(i)} = x^{(i−1)} + α_{i}p̂ and stop solve Mŝ = s t = Aŝ ω_{i} = t^{T} s/t^{T} t x^{(i)} = x^{(i−1)} + α_{i}p̂^{(i)} + ω_{i}ŝ r^{(i)} = s − ω_{i}t check convergence; continue if necessary for continuation it is necessary that ω_{i} ≠ 0

end

Figure 2.10: The Preconditioned BiConjugate Gradient Stabilized Method ω_{i} d > 0 d + c d − c r

r = \frac{a + \sqrt{a^{2} − c^{2}}}{d + \sqrt{d^{2} − c^{2}}} , (2.15)

a [λ_{min}, λ_{max}] λ_{min} λ_{max} M^{−1}A λ_{min} λ_{max} M^{−1}A λ_{max} M^{−1}A w p Compute r^{(0)} = b − Ax^{(0)} for some initial guess x^{(0)}. d = (λ_{max} + λ_{min})/2, c = (λ_{max} − λ_{min})/2. for i = 1, 2, …

solve Mz^{(i−1)} = r^{(i)}. if i = 1

p^{(1)} = z^{(0)}

α_{1} = 2/d else

β_{i−1} = (cα_{i−1}/2)^{2}

α_{i} = 1/(d − β_{i−1})

p^{(i)} = z^{(i−1)} + β_{i−1}p^{(i−1)}. endif x^{(i)} = x^{(i−1)} + α_{i}p^{(i)}. r^{(i)} = b − Ax^{(i)} (= r^{(i−1)} − α_{i}Ap^{(i)}). check convergence; continue if necessary

end

Figure 2.11: The Preconditioned Chebyshev Method M^{−1}A [λ_{min}, λ_{max}]

Matrix-

Inner  Vector Precond
Method  Product SAXPY Product Solve
JACOBI  1^{a}
GS  1  1^{a}
SOR  1  1^{a}
CG  2  3  1  1
GMRES  i + 1  i + 1  1  1
BiCG  2  5  1/1  1/1
QMR  2  8+4^{bc}  1/1  1/1
CGS  2  6  2  2
Bi-CGSTAB  4  6  2  2
CHEBYSHEV  2  1  1

Table 2.1: Summary of Operations for Iteration i. “a/b” means “a” multiplications with the matrix and “b” with its transpose.

aThis method performs no real matrix vector product or preconditioner solve, but the number of operations is equivalent to a matrix-vector multiply.

bTrue SAXPY operations + vector scalings.

cLess for implementations that do not recursively update the residual.

Method  Storage

Reqmts

JACOBI  matrix + 3n
SOR  matrix + 2n
CG  matrix + 6n
GMRES  matrix + (i + 5)n
BiCG  matrix + 10n
CGS  matrix + 11n
Bi-CGSTAB  matrix + 10n
QMR  matrix + 16n^{c}
CHEBYSHEV  matrix + 5n

Table 2.2: Storage Requirements for the Methods in iteration i: n denotes the order of the matrix. ω > 1 T M^{−1}Ax = M^{−1}b M^{−1}A A^{−1} A → M^{−1}A M_{1}^{−1}AM_{2}^{−1}(M_{2}x) = M_{1}^{−1}b.

cLess for implementations that do not recursively update the residual.

M_{1} M_{2} b ← M_{1}^{−1}b M_{1}^{−1}AM_{2}^{−1} x = M_{2}^{−1}y M_{1}^{−1}AM_{2}^{−1} M_{1} = M_{2}^{T} M^{−1}A solve u from Mu = v,

m_{i,j} = \{ \begin{matrix}a_{i,i} & if i = j \\ 0 & otherwise. \end{matrix}

S = {1, …, n} S = ⋃_{i} S_{i} S_{i} m_{i,j} = \{ \begin{matrix}a_{i,j} & if i and j are in the same index subset \\ 0 & otherwise. \end{matrix} P_{n}(x) = x^{n} A = D + L + L^{T} M = (D + L)D^{−1}(D + L)^{T} ,

M(ω) = \frac{1}{2 − ω} (\frac{1}{ω} D + L)(\frac{1}{ω} D)^{−1}(\frac{1}{ω} D + L)^{T} .

κ(M_{ω_{opt}}^{−1}A) = O(\sqrt{κ(A)})

M = LU A = LU Q A = LQ L^{T} A = LL^{T} A ≠ LU (i, j) i + j q a_{i,j} = 0 j < i − p j > i + q A + αI S a_{i,j} ≠ 0 ILU(0) k + 1

for each k, i, j > k: a_{i,j} ← \{ \begin{matrix}a_{i,j} − a_{i,k}a_{k,k}^{−1}a_{k,j} & if (i, j) ∈ S \\ a_{i,j} & otherwise. \end{matrix}

A = D + L + U M = (D + L)D^{−1}(D + U)

(i − 1, j)

(i, j − 1)

d_{i,i} = \{ \begin{matrix}a_{1,1} & if i = 1 \\ a_{i,i} − a_{i,i−1}d_{i−1}^{−1}a_{i−1,i} & if 1 < i ≤ n \\ a_{i,i} − a_{i,i−n}d_{i−n}^{−1}a_{i−n,i} & if i = kn + 1 with k > 1 \\ \begin{aligned}a_{i,i} & − & a_{i,i−1}d_{i−1}^{−1}a_{i−1,i} \\ & & − a_{i,i−n}d_{i−n}^{−1}a_{i−n,i} \end{aligned} & otherwise. \end{matrix}

Initially: d_{i,i} = a_{i,i} for all i for i = 1..nm do:

d_{i+1,i+1} = d_{i+1,i+1} − a_{i+1,i}d_{i,i}^{−1}a_{i,i+1} if there is no k

such that i = kn

d_{i+n,i+n} = d_{i+n,i+n} − a_{i+n,i}d_{i,i}^{−1}a_{i,i+n} if i + n ≤ nm

a_{i,k}a_{k,k}^{−1}a_{k,j} a_{i,i} O(h^{−2})

A + O(h^{2})D_{A} κ(M^{−1}A) = O(h^{−1})

0 < α < 1 n × n i + j = k i + j = k − 1 I − L

(I − L)^{−1} = \{ \begin{matrix}I + L + L^{2} + L^{3} + ⋯ \\ (I + L)(I + L^{2})(I + L^{4}) ⋯ \end{matrix}

d_{i,i} = 1/a_{i,i} (A_{i,i})

{X_{i}}_{i=1..n} X_{1} = A_{1,1} for i ≥ 1

let Y_{i} ≈ X_{i}^{−1}

X_{i+1} = A_{i+1,i+1} − A_{i+1,i}Y_{i}A_{i,i+1} {Y_{i}}_{i=1..n} A = (D + L)D^{−1}(D + U) = (D + L)(I + D^{−1}U)

(I + LD^{−1})(D + U)

X = diag(X_{i})

Y = diag(Y_{i})

C = (X + L)(I + X^{−1}U)

C = (Y^{−1} + L)(I + Y U).

X_{i} C Y_{i} L^{−1}AA^{T} L^{−T} y = b A = I − B A^{−1} = ∑_{k=0}^{∞} B^{k} A = M − N

M_{p}^{−1} = (∑_{i=0}^{p−1}(I − M^{−1}A)^{i})M^{−1}.

kp M = P_{n}(A)

P (0) = 1 ∥I − M^{−1}A∥ max_{i} ∑_{j} |A_{i,j}| (A + A^{T})/2

c_{1} ≤ \frac{x^{T} Ax}{x^{T} Bx} ≤ c_{2} for all x,

c_{1} c_{2} −(a(x, y)u_{x})_{x} − (b(x, y)u_{y})_{y} = f, −(ã(x)u_{x})_{x} − (b̃(y)u_{y})_{y} = f.

L = L_{1} + L_{2}, where L_{1} = −\frac{∂^{2}}{∂x^{2}} , L_{2} = −\frac{∂^{2}}{∂y^{2}} .

L_{1} L_{2} L_{1} L_{2} L_{1} + L_{2} = (I + L_{1})(I + L_{2}) − I − L_{1}L_{2} (1 + αL_{1})(1 + αL_{2})u^{(m+1)} = [(1 + βL_{1})(1 + βL_{2})] u^{(m)} β u^{(m)} u^{(m)} I + αL_{1} A^{H}Ax = A^{H}b (x, y) = x̄^{T} y (x, y) = x^{T} y x^{T} x = 0 x ≠ 0 δA δb x^{(i)} (A + δA)x^{(i)} = b + δb {x^{(i)}} e^{(i)} ≡ x^{(i)} − x maxit ∥b∥ ∥A∥ b r^{(i)} = Ax^{(i)} − b x^{(i)} 10^{−6} ±10^{−6}∥A∥ ±10^{−6}∥b∥ ε ε = 2^{−24} ≈ 10^{−7} ε = 2^{−53} ≈ 10^{−16} i = 0 repeat

i = i + 1 Compute the approximate solution x^{(i)}. Compute the residual r^{(i)} = Ax^{(i)} − b. Compute ∥r^{(i)}∥ and ∥x^{(i)}∥.

until i ≥ maxit or ∥r^{(i)}∥ ≤ stop tol · (∥A∥ · ∥x^{(i)}∥ + ∥b∥). x^{(i)} ∥x^{(i)}∥ until i ≥ maxit or ∥r^{(i)}∥ ≤ stop tol · ∥b∥ .

∥e^{(i)}∥ ≤ ∥A^{−1}∥ · ∥r^{(i)}∥ ∥A^{−1}∥ until i ≥ maxit or ∥r^{(i)}∥ ≤ stop tol · ∥x^{(i)}∥/∥A^{−1}∥ , ∥e^{(i)}∥/∥x^{(i)}∥ e^{(i)} = x^{(i)} − x e^{(i)} r^{(i)} = Ax^{(i)} − b e^{(i)} r^{(i)} e^{(i)} r^{(i)}

∥x∥_{∞} ≡ max_{j} |x_{j}| ,
∥x∥_{1} ≡ ∑_{j} |x_{j}| , and
∥x∥_{2} ≡ (∑_{j} |x_{j}|^{2})^{1/2} .

∥x∥_{B,α} ≡ ∥Bx∥_{α} ∞

∥A∥_{∞} ≡ max_{j} ∑_{k} |a_{j,k}| ,
∥A∥_{1} ≡ max_{k} ∑_{j} |a_{j,k}| , and
∥A∥_{E} ≡ (∑_{jk} |a_{j,k}|^{2})^{1/2} ,

∥A∥_{B,α} ≡ ∥BAB^{−1}∥_{α} ∥A∥_{2} = (λ_{max}(AA^{T}))^{1/2} λ_{max} ∥x∥ ∥x∥_{2} ∥A∥_{E} ∥A∥_{2} ∥x + y∥ ≤ ∥x∥ + ∥y∥ ∥A + B∥ ≤ ∥A∥ + ∥B∥ ∥Ax∥ ≤ ∥A∥ · ∥x∥ ∥x∥_{∞} ≤ 1

\sqrt{n}

∥x∥_{1}

\sqrt{n}

∥e^{(i)}∥ max{∥δA∥/∥A∥, ∥δb∥/∥b∥} x^{(i)} (A + δA)x^{(i)} = (b + δb) δ r^{(i)} = Ax^{(i)} − b e^{(i)} = x^{(i)} − x = A^{−1}(Ax^{(i)} − b) = A^{−1}r^{(i)}, ∥e^{(i)}∥ ≤ ∥A^{−1}∥ · ∥r^{(i)}∥ ∥r^{(i)}∥ ≤ τ ∥e^{(i)}∥ ≤ τ · ∥A^{−1}∥ ∥e^{(i)}∥ ≤ ∥ |A^{−1}| · |r^{(i)}| ∥ |X| X x^{(i)} (A + δA)x^{(i)} = b + δb ∥δA∥ ≤ ε∥A∥ ∥δb∥ ≤ ε∥b∥ ∥r^{(i)}∥ ≤ S_{1} ≡ stop tol · (∥A∥ · ∥x^{(i)}∥ + ∥b∥)

∥δA∥ ≤ stop tol · ∥A∥ ∥δb∥ ≤ stop tol · ∥b∥ ∥e^{(i)}∥ ≤ ∥A^{−1}∥ · ∥r^{(i)}∥ ≤ stop tol · ∥A^{−1}∥ · (∥A∥ · ∥x^{(i)}∥ + ∥b∥) .

∥r^{(i)}∥ ≤ S_{2} ≡ stop tol · ∥b∥ δA = 0 ∥δb∥ ≤ tol · ∥b∥ ∥A∥ · ∥x∥ ≫ ∥b∥

1 ≪ \frac{∥A∥ · ∥x∥}{∥b∥} = \frac{∥A∥ · ∥A^{−1}b∥}{∥b∥} ≤ ∥A∥ · ∥A^{−1}∥ .

∥e^{(i)}∥ ≤ ∥A^{−1}∥ · ∥r^{(i)}∥ ≤ stop tol · ∥A^{−1}∥ · ∥b∥ ∥A^{−1}∥ ∥A^{−1}∥ · ∥r^{(i)}∥ ∥r^{(i)}∥ ≤ S_{3} ≡ stop tol · ∥x^{(i)}∥/∥A^{−1}∥

\frac{∥e^{(i)}∥}{∥x^{(i)}∥} ≤ \frac{∥A^{−1}∥ · ∥r^{(i)}∥}{∥x^{(i)}∥} ≤ stop tol ,

x^{(i)} ∥δA∥ ∥δb∥ δa_{j,k} δb_{j} S_{4} ≡ max_{j}(|r^{(i)}|_{j}/(E · |x^{(i)}| + f)_{j}) ≤ stop tol E f |z| (A + δA)x^{(i)} = b + δb |δa_{j,k}| ≤ tol · e_{j,k} |δb_{j}| ≤ tol · f_{j} e_{j,k} = ∥A∥_{∞} f_{j} = ∥b∥_{∞} ∥r^{(i)}∥_{∞}/(n∥A∥_{∞}∥x^{(i)}∥_{∞} + ∥b∥_{∞}) e_{j,k} = |a_{j,k}| f_{j} = |b_{j}| x^{(i)} ∥e^{(i)}∥_{∞} ≤ ∥ |A^{−1}| · |r^{(i)}| ∥ ≤ S_{4} · ∥ |A^{−1}|(E|x^{(i)}| + f)∥_{∞} |A^{−1}| A^{−1} ∥r^{(i)}∥ ≤ S_{5} ≡ stop tol · ∥r^{(0)}∥ x^{(0)} x^{(0)} = 0 r^{(0)} = b ∥b∥ ≪ ∥A∥ · ∥x∥ x^{(0)} ∥r^{(0)}∥ S_{5} ∥e^{(i)}∥ ≤ S_{5} · ∥A^{−1}∥ r^{(i)} ∥r^{(i)}∥ r^{(i)} = Ax^{(i)} − b ∥r^{(i)}∥ x^{(i)} = M^{−1}Nx^{(i−1)} + M^{−1}b ≡ Gx^{(i−1)} + c r̂^{(i)} = x^{(i)} − Gx^{(i)} − c = M^{−1}(Ax^{(i)} − b) = M^{−1}r^{(i)}

r̂^{(i)}

M^{−1}Ax = M^{−1}b

r̂^{(i)}

∥e^{(i)}∥ = ∥A^{−1}Mr̂^{(i)}∥ ≤ ∥A^{−1}M∥ · ∥r̂^{(i)}∥ ∥A^{−1}M∥ A^{−1}M = (M − N)^{−1}M = (I − G)^{−1} ∥A^{−1}M∥ = ∥(I − G)^{−1}∥ ≤ 1/(1 − ∥G∥)

∥r^{(i)}∥_{M^{−1/2},2} = (r^{(i)T} M^{−1}r^{(i)})^{1/2} ∥r^{(i)}∥_{2} ∥r^{(i)}∥_{M^{−1/2},2}/∥r^{(0)}∥_{M^{−1/2},2} ≤ tol ∥e^{(i)}∥ ≤ ∥A^{−1}M^{1/2}∥ · ∥r^{(i)}∥_{M^{−1/2},2} ∥A^{−1}∥ ∥e^{(i)}∥ A^{−1} e^{(i)} = A^{−1}r^{(i)} A^{−1} x^{(i)} = M^{−1}Nx^{(i−1)} + M^{−1}b = Gx^{(i−1)} + c, I − G ∥G∥ ∥(I − G)^{−1}∥ ≤ 1/(1 − ∥G∥) λ_{max}(A) λ_{min}(A)

∥A^{−1}∥_{2} = λ_{min}^{−1}(A)

∥A^{−1}∥_{2} = 1/λ_{min}^{1/2}(AA^{T}) = 1/λ_{min}^{1/2}(A^{T} A)

∥A^{−1}∥_{∞} O(n^{2})

O(n^{3})

∥e^{(i)}∥_{∞} ≤ ∥ |A^{−1}| · |r^{(i)}| ∥_{∞} ∥A^{−1}∥_{∞} · ∥r^{(i)}∥_{∞} ∥ |A^{−1}| · |r^{(i)}| ∥_{∞} R |r^{(i)}| ∥ |A^{−1}| · |r^{(i)}| ∥_{∞} = ∥A^{−1}R∥_{∞} ∥A^{−1}R∥_{∞} A^{−1}R (A^{−1}R)^{T} = R^{T} A^{−T} A^{−1} A^{−T} ∥x∥_{∞} = max_{j} |x_{j}| x_{j} ±∞ ∥x∥_{2} = (∑_{j} |x_{j}|^{2})^{1/2} ∥x∥_{∞} r^{(i)} = Ax^{(i)} − b Ax^{(i)} δr^{(i)} r^{(i)} ∥δr^{(i)}∥ ≤ O(ε)(∥A∥ · ∥x^{(i)}∥ + ∥b∥)

O(ε) nε \sqrt{n}ε stop tol ≤ ε r^{(i)} e^{(i)} = A^{−1}r^{(i)} O(ε)∥A^{−1}∥ · (∥A∥ · ∥x^{(i)}∥ + ∥b∥)

(δr^{(i)})_{j} r^{(i)} |A| · |x^{(i)}| + |b| |δr^{(i)}| ≤ O(ε)(|A| · |x^{(i)}| + |b|) e^{(i)} O(ε)∥ |A^{−1}| · (|A| · |x^{(i)}| + |b|)∥ e^{(i)} 3 val(k) = a_{i,j}

col ind(k) = j

val(k) = a_{i,j}

row ptr(i) ≤ k < row ptr(i + 1)

row ptr(n + 1) = nnz + 1 nnz n^{2} 2nnz + n + 1

A = ( \begin{matrix}10 & 0 & 0 & 0 & −2 & 0 \\ 3 & 9 & 0 & 0 & 0 & 3 \\ 0 & 7 & 8 & 7 & 0 & 0 \\ 3 & 0 & 8 & 7 & 5 & 0 \\ 0 & 8 & 0 & 9 & 9 & 13 \\ 0 & 4 & 0 & 0 & 2 & −1 \end{matrix} ) . (4.1)

val 10 -2 3 9 3 7 8 7 3 ⋯ 9 13 4 2 -1
col ind 1 5 1 2 6 2 3 4 1 ⋯ 5 6 2 5 6

row ptr 1 3 6 9 13 17 20

val 10 3 3 9 7 8 4 8 8 ⋯ 9 2 3 13 -1
row ind 1 2 4 2 3 5 6 3 4 ⋯ 5 6 2 5 6

col ptr 1 4 8 10 13 17 20

n_{b} nnzb nnz = nnzb × n_{b}^{2} n_{d} n_{d} = n/n_{b} 1 : nnzb 1 : n_{b} 1, 1 1 : n_{d} + 1 A = (a_{i,j}) a_{i,j} ≠ 0 i − p ≤ j ≤ i + q

p + q

A = ( \begin{matrix}10 & −3 & 0 & 0 & 0 & 0 \\ 3 & 9 & 6 & 0 & 0 & 0 \\ 0 & 7 & 8 & 7 & 0 & 0 \\ 0 & 0 & 8 & 7 & 5 & 0 \\ 0 & 0 & 0 & 9 & 9 & 13 \\ 0 & 0 & 0 & 0 & 2 & −1 \end{matrix} ) . (4.2)

val(i, j) = a_{i,i+j}. (4.3)

val(:,-1) 0 3 7 8 9 2
val(:, 0) 10 9 8 7 9 -1
val(:,+1) -3 6 7 5 13 0

S = {(i, σ(i)); i ∈ I ⊆ I_{n}} I_{n} = {1, …, n} σ (i, σ(i))

(j, σ(j)) i < j → σ(i) < σ(j). y = Ax (i, σ_{k}(i))

S_{k} x_{σ_{k}(i)} y_{σ_{k}(i)} y_{i}

A = ( \begin{matrix}10 & −3 & 0 & 1 & 0 & 0 \\ 0 & 9 & 6 & 0 & −2 & 0 \\ 3 & 0 & 8 & 7 & 0 & 0 \\ 0 & 6 & 0 & 7 & 5 & 4 \\ 0 & 0 & 0 & 0 & 9 & 13 \\ 0 & 0 & 0 & 0 & 5 & −1 \end{matrix} ) , (4.4)

4

val(:,-1) 0 0 3 6 0 5
val(:, 0) 10 9 8 7 9 -1
val(:,+1) 0 -3 6 7 5 13
val(:,+2) 0 1 -2 0 4 0

( \begin{matrix}10 & −3 & 0 & 1 & 0 & 0 \\ 0 & 9 & 6 & 0 & −2 & 0 \\ 3 & 0 & 8 & 7 & 0 & 0 \\ 0 & 6 & 0 & 7 & 5 & 4 \\ 0 & 0 & 0 & 0 & 9 & 13 \\ 0 & 0 & 0 & 0 & 5 & −1 \end{matrix} ) ⟶ ( \begin{matrix}10 & −3 & 1 \\ 9 & 6 & −2 \\ 3 & 8 & 7 \\ 6 & 7 & 5 & 4 \\ 9 & 13 \\ 5 & −1 \end{matrix} )

val(:, 1) 10 9 3 6 9 5 val(:, 2) −3 6 8 7 13 −1 val(:, 3) 1 −2 7 5 0 0 val(:, 4) 0 0 0 4 0 0 , colind(:, 1) 1 2 1 2 5 5 colind(:, 2) 2 3 3 4 6 6 colind(:, 3) 4 5 4 5 0 0 colind(:, 4) 0 0 0 6 0 0 .

jdiag 1 3 7 8 10 2; 9 9 8 ⋯ -1; 9 6 7 5; 13
col ind 1 1 2 3 1 5; 4 2 3 ⋯ 6; 5 3 4 5; 6

perm 5 2 3 4 1 6

jd ptr 1 7 13 17

+ x x
x + x x
x x + x  x
x x + x x x x
x + x x x
x x x + x x
x x + x x x
x x x + x x x  x
x x x x + x x  x
x + x x x
x x x x + x x
x x + x  x
x x x + x x x
x x + x x
+ x x
x x x + x

+

Figure 4.1: Profile of a nonsymmetric skyline or variable-band matrix.

y = Ax  and  y = A^{T} x.

y_{i} = ∑_{j} a_{i,j}x_{j},

2n^{2} y = A^{T} x

y_{i} = ∑_{j}(A^{T})_{i,j}x_{j} = ∑_{j} a_{j,i}x_{j},

for all j, do for all i: y_{i} ← y_{i} + a_{j,i}x_{j}. j → i + j y_{i} ← y_{i} + a_{i,j}x_{j} ⇒ y_{i} ← y_{i} + a_{i,i+j}x_{i+j} . a_{i,i+j} 1 ≤ i, i + j ≤ n.

y_{i} ← y_{i} + a_{i+j,i}x_{j}
= y_{i} + a_{i+j,i+j−j}x_{i+j}

A = D_{A} + L_{A} + U_{A} (D_{A} + L_{A})D_{A}^{−1}(D_{A} + U_{A})

LUy = x  and  (LU)^{T} y = x.

val(diag ptr(i)) = a_{i,i}

a_{i,j} j > i a_{j,i} a_{i,j} a_{j,j} LUy = x Lz = x, Uy = z.

LU = (D + L_{A})D^{−1}(D + U_{A})
= (I + L_{A}D^{−1})(D + U_{A})
= (D + L_{A})(I + D^{−1}U_{A})
= (I + L_{A}D^{−1})D(I + D^{−1}U_{A})

(LU)^{T} y = x U^{T} L = (l_{∗1}, l_{∗2}, …)

Ly = x x = l_{∗1}y_{1} + l_{∗2}y_{2} + ⋯ y_{1} x ← x − l_{∗1}y_{1} Ap^{(i)} A^{T} p^{(i)} x^{(i+1)} r^{(i+1)} α_{i} p^{(i)} β_{i−1} p^{(i)^{T}} Ap^{(i)} M = LL^{T} L^{−1}r^{(i)} r^{(i)^{T}} M^{−1}r^{(i)} s^{T} s s = L^{−1}r^{(i)} L^{−T} s x^{(−1)} = x^{(0)}= initial guess; r^{(0)} = b − Ax^{(0)}; p^{(−1)} = 0; β_{−1} = 0; α_{−1} = 0; s = L^{−1}r^{(0)}; ρ_{0} = (s, s) for i = 0, 1, 2, ....

w^{(i)} = L^{−T} s; p^{(i)} = w^{(i)} + β_{i−1}p^{(i−1)}; q^{(i)} = Ap^{(i)}; γ = (p^{(i)}, q^{(i)}); x^{(i)} = x^{(i−1)} + α_{i−1}p^{(i−1)}; α_{i} = ρ_{i}/γ; r^{(i+1)} = r^{(i)} − α_{i}q^{(i)}; s = L^{−1}r^{(i+1)}; ρ_{i+1} = (s, s); if ∥r^{(i+1)}∥ small enough then

x^{(i+1)} = x^{(i)} + α_{i}p^{(i)} quit;

endif β_{i} = ρ_{i+1}/ρ_{i};

end;

Figure 4.2: A rearrangement of Conjugate Gradient for parallelism γ x^{(i)} ρ_{i+1} p^{(i)} q^{(i)} p^{(i)} q^{(i)} LDL^{T} a_{i,j} Ap_{i} p_{i} D − L i + d M^{−1}Av^{(j)} v^{(1)} v^{(2)} v^{(j)} v^{(1)} Av^{(1)} A^{m}v^{(1)} n × k R^{(k)} R_{k} = [r^{(0)}, r^{(1)}, …, r^{(k−1)}], k × k B_{k}

B_{k} = [ \begin{matrix}1 & −β_{2} & & ⋯ & 0 \\ & 1 & −β_{3} & & ⋮ \\ & ... & ... & ... \\ ⋮ & & ... & ... & −β_{k} \\ 0 & ⋯ & & & 1 \end{matrix} ] ,

{r^{(k)}} {β_{k}} p^{(j)} = r^{(j−1)} + β_{j}p^{(j−1)}, j = 2, 3, …, k , p^{(1)} = r^{(0)} R_{k} = P_{k}B_{k} P_{k} = [p^{(1)}, p^{(2)}, …, p^{(k)}].

{p^{(j)}} T̂_{k} = R_{k}^{T} AR_{k} = B_{k}^{T} Λ̂_{k}B_{k}

Λ̂_{k} = [ \begin{matrix}p^{(1)^{T}} Ap^{(1)} & 0 & & ⋯ & 0 \\ 0 & p^{(2)^{T}} Ap^{(2)} & & & ⋮ \\ & ... & ... & ... \\ ⋮ & & ... & ... & 0 \\ 0 & ⋯ & & 0 & p^{(k)^{T}} Ap^{(k)} \end{matrix} ] .

p^{(1)}, p^{(2)}, …, p^{(j)} r^{(0)}, r^{(1)}, …, r^{(j−1)} {r^{(j)}} Q_{k} = R_{k}Δ^{−1} span{b, Ab, …, A^{k−1}b} Δ ∥r^{(i)}∥_{2} Q_{k} T_{k} = Δ^{−1}B_{k}^{T} Λ_{k}B_{k}Δ^{−1} , (5.1)

Λ_{k} p^{(i)^{T}} p^{(i)} T_{k} x^{(k)} n_{b} − 1 M^{−1}Ax = M^{−1}f Ax = f Mz = r 5

Ax = [ \begin{matrix}D_{R} & C \\ C^{T} & D_{B} \end{matrix} ] [ \begin{matrix}x_{R} \\ x_{B} \end{matrix} ] = [ \begin{matrix}f_{R} \\ f_{B} \end{matrix} ] , (5.2)

D_{R} D_{B}

[ \begin{matrix}I & O \\ −C^{T} D_{R}^{−1} & I \end{matrix} ] [ \begin{matrix}D_{R} & C \\ C^{T} & D_{B} \end{matrix} ] [ \begin{matrix}x_{R} \\ x_{B} \end{matrix} ] = [ \begin{matrix}I & O \\ −C^{T} D_{R}^{−1} & I \end{matrix} ] [ \begin{matrix}f_{R} \\ f_{B} \end{matrix} ]

[ \begin{matrix}D_{R} & C \\ O & D_{B} − C^{T} D_{R}^{−1}C \end{matrix} ] [ \begin{matrix}x_{R} \\ x_{B} \end{matrix} ] = [ \begin{matrix}f_{R} \\ f_{B} − C^{T} D_{R}^{−1}f_{R} \end{matrix} ] .

D_{B}^{−1/2}

(I − H^{T} H)y = g , (5.3)

H = D_{R}^{−1/2}CD_{B}^{−1/2} y = D_{B}^{1/2}x_{B} g = D_{B}^{−1/2}(f_{B} − C^{T} D_{R}^{−1}f_{R}) n^{2}/2 (n^{2} − 1)/2 n ≥ 250

−∇ · (a(x, y)∇u) = f(x, y)  (5.4)

a(x, y)

Ω u T_{H} Ω_{i}, i = 1, ..., p n_{H} Ω_{i} T_{h} H, h Au = f S_{B}u_{B} = g_{B} Av Sw v K^{−1}v Ω′_{i} n′_{i} T′_{i} ⊂ T_{h} A′_{i}, A_{H} T′_{i} R_{i}^{T} R_{i} R_{H}^{T} R_{H} K_{as}

K_{as}^{−1}v = R_{H}^{T}A_{H}^{−1}R_{H}v + ∑_{i=1}^{p} R_{i}^{T} A′_{i}^{−1}R_{i}v.

A′_{i} A_{H} R_{i}^{T} A′_{i}^{−1}R_{i}v v_{i} = R_{i}v w_{i} A′_{i}w_{i} = v_{i} y_{i} = R_{i}^{T} w_{i} K_{as}^{−1}A H δ/H ∂Ω′_{i} ∂Ω_{i} h → 0 O(1/H^{2})

K_{as} I_{i} n′_{i} × n′_{i} n × n′_{i} u_{i} (R_{i}^{T} u_{i})_{j} = (u_{i})_{j} j ∈ I_{i} R_{i}u n × n_{H} R_{H}, R_{H}^{T} R_{i}^{T} , R_{i}, R_{H}^{T}, R_{H} K_{as} R_{i}, A′_{i} A_{H}, R_{H} I ⋃_{i} ∂Ω′_{i} u ≡ (u_{I}, u_{B})^{T} f ≡ (f_{I}, f_{B})^{T}

( \begin{matrix}A_{I} & A_{IB} \\ A_{IB}^{T} & A_{B} \end{matrix} ) ( \begin{matrix}u_{I} \\ u_{B} \end{matrix} ) = ( \begin{matrix}f_{I} \\ f_{B} \end{matrix} ) .  (5.5)

A_{I} = blockdiagonal(A_{i})

A_{i}

( \begin{matrix}I & 0 \\ A_{IB}^{T}A_{I}^{−1} & I \end{matrix} ) ( \begin{matrix}A_{I} & A_{IB} \\ 0 & S_{B} \end{matrix} ) ( \begin{matrix}u_{I} \\ u_{B} \end{matrix} ) = ( \begin{matrix}f_{I} \\ f_{B} \end{matrix} ) ,  (5.6)

S_{B} ≡ A_{B} − A_{IB}^{T}A_{I}^{−1}A_{IB} A_{B} u_{I} u_{B}

S_{B}u_{B} = g_{B} ≡ f_{B} − A_{IB}A_{I}^{−1}f_{I}.  (5.7)

S_{B} O(h^{−1})

O(h^{−2}) v_{B} S_{B}v_{B} A_{B}v_{B} − A_{IB}^{T}(A_{I}^{−1}(A_{IB}v_{B}))

A_{I}^{−1} A_{i}^{−1} g_{B}

B = E ∪ V_{H}  (5.8)

V_{H} = ⋃_{k} V_{k} V_{k} E = ⋃_{i} E_{i} E_{i} V_{H} R_{E_{i}} n_{E_{i}} × n n_{E_{i}} (R_{E_{i}}u_{B})_{j} = (u_{B})_{j} j ∈ E_{i}

S_{B} = ( \begin{matrix}S_{E} & S_{EV} \\ S_{EV}^{T} & S_{V} \end{matrix} ) .  (5.9)

S_{E} S_{E_{i}} S_{E} S_{E_{i}} E_{i} S_{E_{i}}^{−1} S_{E_{i}} S_{E_{i}} S_{E_{i}} M_{E_{i}} = α_{E_{i}}K^{1/2} K n_{E_{i}} × n_{E_{i}} −1 α_{E_{i}} M_{E_{i}}

K_{BP S}^{−1}v_{B} = R_{H}^{T}A_{H}^{−1}R_{H}v_{B} + ∑_{E_{i}} R_{E_{i}}^{T}M_{E_{i}}^{−1}R_{E_{i}}v_{B}.  (5.10)

K_{BP S}^{−1}S_{B} O(1 + log^{2}(H/h)) log^{2}(H/h)

K_{BP S} K_{BP S} X = ⋃_{k} X_{k} X_{k} S_{X_{k}} R_{X_{k}}, R_{X_{k}}^{T} S_{X_{k}} K^{1/2} M_{X_{k}} S_{X_{k}}

K_{V S}^{−1}v_{B} = R_{H}^{T}A_{H}^{−1}R_{H}v_{B} + ∑_{E_{i}} R_{E_{i}}^{T}M_{E_{i}}^{−1}R_{E_{i}}v_{B}
+ ∑_{X_{k}} R_{X_{k}}^{T}M_{X_{k}}^{−1}R_{X_{k}}v_{B}.  (5.11)

K_{V S}^{−1}S_{B} A′_{i}^{−1}, A_{i}^{−1} A_{H}^{−1} K_{as}, K_{BP S}, K_{V S} Ã′_{i}^{−1}, Ã_{i}^{−1}

˜ A_{H}^{−}

K̃_{as}^{−1}v = R_{H}^{T} Ã_{H}^{−1}R_{H}v + ∑_{i=1}^{p} R_{i}^{T} Ã′_{i}^{−1}R_{i}v.

A_{H}^{−1}

˜ A_{H}^{−1}

K_{BP S} K_{V S}

˜ K_{BP S}

˜ K_{V S}

A_{I}^{−1}

˜ A_{i}

A_{I}^{−1} A_{I}, S_{B} Ã_{I}, S̃_{B}

˜ S_{B}

˜ K_{BP S}

˜ K_{V S}

˜ A_{H}

˜ A_{i}

Ã_{H}^{−1}A_{H} Ã_{i}^{−1}A_{i} A_{IB} 2 × 2

A^{(i)} = ( \begin{matrix}A_{1,1}^{(i)} & A_{1,2}^{(i)} \\ A_{2,1}^{(i)} & A_{2,2}^{(i)} \end{matrix} )

A^{(i+1)} ≈ S^{(i)} = A_{2,2}^{(i)} − A_{2,1}^{(i)}A_{1,1}^{(i)^{−1}}A_{1,2}^{(i)}.

A^{T} = [A_{1}, …, A_{m}] P_{i}x = A_{i}(A_{i}^{T} A_{i})^{−1}A_{i}^{T} x.

shar filename  contents
sctemplates.shar Single precision C routines
dctemplates.shar Double precision C routines
sftemplates.shar Single Precision Fortran 77 routines
dftemplates.shar Double Precision Fortran 77 routines
mltemplates.shar MATLAB routines

T x = b a_{i,j} = Ax + b b − Ax̂

ˆ x

rk T x x^{(i)} (A + δA)x^{(i)} = b + δb a_{i,j} = 0 {x^{T} Ax : x^{T} x = 1} [λ_{min}(A), λ_{max}(A)] {A^{i}x}_{i≥0} a_{i,j} = 0 a_{i,j} = a_{j,i} a_{i,i} > ∑_{j≠i} |a_{i,j}| min_{i}{a_{i,i} − ∑_{j≠i} |a_{i,j}|} a_{i,j} ≤ 0 (A^{−1})_{i,j} ≥ 0 AA^{T} y = b A^{H} f : R^{n} → R f(x) ≥ 0 f(x) = 0 x = 0 f(αx) = |α|f(x) f(x + y) ≤ f(x) + f(y)

∥ · ∥ ∥Ax∥ ≤ ∥A∥ ∥x∥.

∥A∥_{2} · ∥A^{−1}∥_{2} = \frac{λ_{max}^{1/2}(A^{T} A)}{λ_{min}^{1/2}(A^{T} A)},

λ_{max} λ_{min} a_{i,j} = 0 j < i m × n

A = [ \begin{matrix}a_{1,1} & ⋯ & a_{1,n} \\ ⋮ & & ⋮ \\ a_{m,1} & ⋯ & a_{m,n} \end{matrix} ] a_{i,j} ∈ R.

x = [ x_{1} ; ⋮ ; x_{n} ] x_{i} ∈ R.

I^{n×n} (a_{i,i}) a_{i,i} x_{i}^{(k)} s ℓ