Skip to content

Commit 6051bf0

Browse files
author
Pawel Latawiec
committed
Look-ahead Lanczos construction
1 parent 1c5c21b commit 6051bf0

File tree

4 files changed

+348
-0
lines changed

4 files changed

+348
-0
lines changed

‎src/IterativeSolvers.jl

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,9 @@ include("history.jl")
1313
# Factorizations
1414
include("hessenberg.jl")
1515

16+
# Krylov subspace methods
17+
include("lal.jl")
18+
1619
# Linear solvers
1720
include("stationary.jl")
1821
include("stationary_sparse.jl")

‎src/lal.jl

Lines changed: 212 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,212 @@
1+
"""
2+
LookAheadLanczosDecompOptions
3+
4+
Options for [`LookAheadLanczosDecomp`](@ref).
5+
6+
# Fields
7+
- `max_iter`: Maximum number of iterations
8+
- `max_block_size`: Maximum block size allowed to be constructed
9+
- `log`: Flag determining logging
10+
- `verbose`: Flag determining verbosity
11+
"""
12+
struct LookAheadLanczosDecompOptions
13+
max_iter::Int
14+
max_block_size::Int
15+
log::Bool
16+
verbose::Bool
17+
end
18+
19+
"""
20+
LookAheadLanczosDecompLog
21+
22+
Log for [`LookAheadLanczosDecomp`](@ref). In particular, logs the sizes of blocks constructed in the P-Q and V-W sequences.
23+
"""
24+
struct LookAheadLanczosDecompLog
25+
PQ_block_count::Dict{Int, Int}
26+
VW_block_count::Dict{Int, Int}
27+
end
28+
LookAheadLanczosDecompLog() = LookAheadLanczosDecompLog(Dict{Int, Int}(), Dict{Int, Int}())
29+
30+
mutable struct LookAheadLanczosDecomp{OpT, OptT, VecT, MatT, ElT, ElRT}
31+
# Operator
32+
A::OpT
33+
At::OptT
34+
35+
# P-Q sequence
36+
p::VecT
37+
q::VecT
38+
p̂::VecT
39+
q̂::VecT
40+
P::MatT
41+
Q::MatT
42+
43+
# V-W sequence
44+
v::VecT
45+
w::VecT
46+
ṽ::VecT
47+
w̃::VecT
48+
V::MatT
49+
W::MatT
50+
51+
# matrix-vector products
52+
Ap::VecT
53+
Atq::VecT
54+
55+
# dot products - note we take tranpose(w)*v, not adjoint(w)*v
56+
qtAp::ElT
57+
w̃tṽ::ElT
58+
wtv::ElT
59+
60+
# norms
61+
normp::ElRT
62+
normq::ElRT
63+
ρ::ElRT
64+
ξ::ElRT
65+
66+
γ::Vector{ElRT}
67+
68+
D::Matrix{ElT}
69+
E::Matrix{ElT}
70+
F::Matrix{ElT}
71+
F̃lastcol::Vector{ElT}
72+
G::Matrix{ElT}
73+
H::Vector{ElT}
74+
75+
U::UpperTriangular{ElT, Matrix{ElT}}
76+
L::Matrix{ElT}
77+
78+
# Indices tracking location in block and sequence
79+
n::Int
80+
k::Int
81+
l::Int
82+
kstar::Int
83+
lstar::Int
84+
mk::Vector{Int}
85+
nl::Vector{Int}
86+
87+
# Flag determining if we are in inner block, see [^Freund1994] Alg. 5.2
88+
innerp::Bool
89+
innerv::Bool
90+
91+
# Estimate of norm(A), see [^Freund1993]
92+
nA::ElRT
93+
nA_recompute::ElRT
94+
95+
# Logs and options
96+
log::LookAheadLanczosDecompLog
97+
opts::LookAheadLanczosDecompOptions
98+
end
99+
100+
"""
101+
LookAheadLanczosDecomp(A, v, w; kwargs...)
102+
103+
Provides an iterable which constructs basis vectors for a Krylov subspace generated by `A` given by two initial vectors `v` and `w`. This implementation follows [^Freund1994], where a coupled two-term recurrence is used to generate both a V-W and a P-Q sequence. Following the reference, the Lanczos sequence is generated by `A` and `transpose(A)`.
104+
105+
# Arguments
106+
- `A`: Operator used to construct Krylov subspace
107+
- `v`: Initial vector for Krylov subspace generated by `A`
108+
- `w`: Initial vector for Krylov subspace generated by `transpose(A)`
109+
110+
# Keywords
111+
- `max_iter=size(A, 2)`: Maximum number of iterations
112+
- `max_block_size=2`: Maximum look-ahead block size to construct. Following [^Freund1994], it is rare for blocks to go beyond size 3. This pre-allocates the block storage for the computation to size `(max_block_size, length(v))`. If a block would be built that exceeds this size, the estimate of `norm(A)` is adjusted to allow the block to close.
113+
- `log=false`: Flag determining whether to log history in a [`LookAheadLanczosDecompLog`](@ref)
114+
- `verbose=false`: Flag determining verbosity of output during iteration
115+
116+
# References
117+
[^Freund1993]:
118+
Freund, R. W., Gutknecht, M. H., & Nachtigal, N. M. (1993). An Implementation of the Look-Ahead Lanczos Algorithm for Non-Hermitian Matrices. SIAM Journal on Scientific Computing, 14(1), 137–158. https://doi.org/10.1137/0914009
119+
[^Freund1994]:
120+
Freund, R. W., & Nachtigal, N. M. (1994). An Implementation of the QMR Method Based on Coupled Two-Term Recurrences. SIAM Journal on Scientific Computing, 15(2), 313–337. https://doi.org/10.1137/0915022
121+
"""
122+
function LookAheadLanczosDecomp(
123+
A, v, w;
124+
max_iter=size(A, 2),
125+
max_block_size=8,
126+
log=false,
127+
verbose=false
128+
)
129+
elT = eltype(v)
130+
131+
p = similar(v)
132+
q = similar(v)
133+
p̂ = similar(v)
134+
q̂ = similar(v)
135+
P = similar(v, size(v, 1), 0)
136+
Q = similar(v, size(v, 1), 0)
137+
Ap = similar(v)
138+
Atq = similar(v)
139+
qtAp = zero(elT)
140+
normp = zero(real(elT))
141+
normq = zero(real(elT))
142+
143+
ṽ = similar(v)
144+
w̃ = similar(v)
145+
V = reshape(copy(v), size(v, 1), 1)
146+
W = reshape(copy(w), size(v, 1), 1)
147+
w̃tṽ = zero(elT)
148+
wtv = transpose(w) * v
149+
ρ = zero(normp)
150+
ξ = zero(normp)
151+
152+
γ = Vector{real(elT)}(undef, 1)
153+
γ[1] = 1.0
154+
155+
D = Matrix{elT}(undef, 0, 0)
156+
E = Matrix{elT}(undef, 0, 0)
157+
G = Matrix{elT}(undef, 0, 0)
158+
H = Vector{elT}()
159+
160+
F = Matrix{elT}(undef, 0, 0)
161+
F̃lastcol = Vector{elT}()
162+
163+
U = UpperTriangular(Matrix{elT}(undef, 0, 0))
164+
L = Matrix{elT}(undef, 0, 0)
165+
166+
n = 1
167+
k = 1
168+
l = 1
169+
kstar = 1
170+
lstar = 1
171+
mk = [1]
172+
nl = [1]
173+
174+
# Equation following 4.5 of [^Freund1993]
175+
# initial estimate of norm(A)
176+
nA = max(
177+
norm(mul!(Ap, A, v)),
178+
norm(mul!(Atq, transpose(A), w))
179+
)
180+
if log
181+
ld_log = LookAheadLanczosDecompLog(
182+
Dict(i => 0 for i=1:max_block_size),
183+
Dict(i => 0 for i=1:max_block_size)
184+
)
185+
else
186+
ld_log = LookAheadLanczosDecompLog()
187+
end
188+
189+
return LookAheadLanczosDecomp(
190+
A, transpose(A),
191+
p, q, p̂, q̂, P, Q,
192+
v, w, ṽ, w̃, V, W,
193+
Ap, Atq,
194+
qtAp, w̃tṽ, wtv,
195+
normp, normq, ρ, ξ,
196+
γ,
197+
D, E, F, F̃lastcol, G, H,
198+
U, L,
199+
n, k, l, kstar, lstar, mk, nl,
200+
false, false, nA, nA,
201+
ld_log,
202+
LookAheadLanczosDecompOptions(
203+
max_iter,
204+
max_block_size,
205+
log,
206+
verbose
207+
)
208+
)
209+
end
210+
211+
isverbose(ld::LookAheadLanczosDecomp) = ld.opts.verbose
212+
islogging(ld::LookAheadLanczosDecomp) = ld.opts.log

‎test/lal.jl

Lines changed: 130 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,130 @@
1+
using IterativeSolvers; const IS = IterativeSolvers
2+
using LinearAlgebra
3+
using SparseArrays
4+
using Test
5+
6+
# Equation references and identities from:
7+
# Freund, R. W., & Nachtigal, N. M. (1994). An Implementation of the QMR Method Based on Coupled Two-Term Recurrences. SIAM Journal on Scientific Computing, 15(2), 313–337. https://doi.org/10.1137/0915022
8+
9+
function test_intermediate_identities(ld)
10+
# 2.7, 3.23
11+
@test transpose(ld.W[:, 1:end-1]) * ld.V[:, 1:end-1] ≈ ld.D
12+
# 3.14, 3.26
13+
@test transpose(ld.Q) * ld.A * ld.P ≈ ld.E
14+
# 3.10
15+
@test ld.A * ld.V[:, 1:end-1] ≈ ld.V * ld.L * ld.U
16+
# 3.7
17+
@test ld.V[:, 1:end-1] ≈ ld.P * ld.U
18+
# 3.7
19+
@test ld.A * ld.P ≈ ld.V * ld.L
20+
# 5.1
21+
@test ld.G ≈ ld.U * ld.L[1:end-1, 1:end-1]
22+
# 3.11
23+
@test ld.H ≈ ld.L * ld.U
24+
25+
Γ = Diagonal(ld.γ[1:end-1])
26+
# 3.15
27+
@test ld.D * Γ ≈ transpose(ld.D * Γ)
28+
# 3.16
29+
@test ld.E * Γ ≈ transpose(ld.E * Γ)
30+
31+
# 3.8
32+
@test ld.W[:, 1:end-1] ≈ ld.Q * ((Γ \ ld.U) * Γ)
33+
# 3.8
34+
@test ld.At * ld.Q[:, 1:end-1] ≈ ld.W[:, 1:end-1] * (Γ \ ld.L[1:end-1, 1:end-1]) * Γ[1:end-1, 1:end-1]
35+
# 4.1
36+
@test ld.A * ld.P[:, 1:end-1] ≈ ld.P * ld.U * ld.L[1:end-1, 1:end-1]
37+
@test ld.At * ld.Q[:, 1:end-1] ≈ ld.Q * (Γ \ ld.U) * ld.L[1:end-1, 1:end-1] * Γ[1:end-1, 1:end-1]
38+
39+
# 3.9
40+
@test all(diag(ld.U) .≈ 1)
41+
42+
# Definition, by 5.1
43+
F = transpose(ld.W[:, 1:end-1]) * ld.A * ld.P
44+
F̃ = transpose(ld.Q) * (A * ld.V[:, 1:end-1])
45+
# Lemma 5.1
46+
@test F * Γ ≈ transpose(F̃ * Γ)
47+
@test ld.F[:, 1:end-1] ≈ F[:, 1:end-1]
48+
49+
@test all([norm(ld.V[:, i]) for i in axes(ld.V, 2)] .≈ 1)
50+
@test all([norm(ld.W[:, i]) for i in axes(ld.W, 2)] .≈ 1)
51+
end
52+
53+
function test_finished_identities(ld)
54+
# 2.7, 3.23
55+
@test transpose(ld.W) * ld.V ≈ ld.D
56+
# 3.14, 3.26
57+
@test transpose(ld.Q[:, 1:end-1]) * ld.A * ld.P[:, 1:end-1] ≈ ld.E
58+
# 3.10
59+
@test ld.A * ld.V[:, 1:end-1] ≈ ld.V * ld.L * ld.U[1:end-1, 1:end-1]
60+
# 3.7
61+
@test ld.V ≈ ld.P * ld.U
62+
# 3.7
63+
@test ld.A * ld.P[:, 1:end-1] ≈ ld.V * ld.L
64+
# 5.1
65+
@test ld.G ≈ ld.U * ld.L[1:end-1, 1:end-1]
66+
# 3.11
67+
@test ld.H ≈ ld.L * ld.U
68+
69+
Γ = Diagonal(ld.γ)
70+
# 3.15
71+
@test ld.D * Γ ≈ transpose(ld.D * Γ)
72+
# 3.16
73+
@test ld.E * Γ[1:end-1, 1:end-1] ≈ transpose(ld.E * Γ[1:end-1, 1:end-1])
74+
75+
# 3.8
76+
@test ld.W ≈ ld.Q * ((Γ \ ld.U) * Γ)
77+
# 3.8
78+
@test ld.At * ld.Q[:, 1:end-1] ≈ ld.W * (Γ \ ld.L) * Γ[1:end-1, 1:end-1]
79+
# 4.1
80+
@test ld.A * ld.P[:, 1:end-1] ≈ ld.P * ld.U * ld.L
81+
@test ld.At * ld.Q[:, 1:end-1] ≈ ld.Q * (Γ \ ld.U) * ld.L * Γ[1:end-1, 1:end-1]
82+
83+
# 3.9
84+
@test all(diag(ld.U) .≈ 1)
85+
86+
# Definition, by 5.1
87+
F = transpose(ld.W) * ld.A * ld.P
88+
F̃ = transpose(ld.Q) * (A * ld.V)
89+
# Lemma 5.1
90+
@test F * Γ ≈ transpose(F̃ * Γ)
91+
@test ld.F[:, 1:end-1] ≈ F[:, 1:end-1]
92+
@test ld.F̃lastcol ≈ F̃[1:end-1, end]
93+
94+
# 3.13
95+
@test ld.L ≈ (ld.D \ Γ) * transpose(ld.U) * (Γ \ transpose(ld.Q)) * A * ld.P[:, end]
96+
97+
# 3.35
98+
@test ld.U[1:end-1, 1:end-1] ≈ (ld.E \ transpose(ld.Q[1:end-1])) * (ld.A * ld.V[:, 1:end-1])
99+
# 3.42
100+
# Likewise, similar to 3.13 above
101+
@test ld.L ≈ (ld.D \ transpose(ld.W)) * (ld.A * ld.P[:, 1:end-1])
102+
103+
@test all([norm(ld.V[:, i]) for i in axes(ld.V, 2)] .≈ 1)
104+
@test all([norm(ld.W[:, i]) for i in axes(ld.W, 2)] .≈ 1)
105+
end
106+
107+
function test_regular_identities(ld)
108+
# If all regular vectors, expect:
109+
# 3.22
110+
@test abs(sum(triu(ld.U, 2))) < 1e-8
111+
# 3.20
112+
@test abs(sum(triu(ld.L, 1))) < 1e-8
113+
@test abs(sum(triu(ld.D, 1)) + sum(tril(ld.D, -1))) < 1e-8
114+
@test abs(sum(triu(ld.E, 1)) + sum(tril(ld.E, -1))) < 1e-8
115+
@test ld.log.PQ_block_count[1] == ld.n
116+
@test ld.log.VW_block_count[1] == ld.n - 1
117+
end
118+
119+
@testset "A = I" begin
120+
# A = I terminates immediately (because p1 = v1 -> v2 = Ap1 - v1 = 0)
121+
A = Diagonal(fill(1.0, 5))
122+
v = rand(5)
123+
v ./= norm(v)
124+
w = rand(5)
125+
w ./= norm(w)
126+
ld = IS.LookAheadLanczosDecomp(A, v, w)
127+
for l in ld end
128+
129+
@test ld.n == 1
130+
end

‎test/runtests.jl

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,9 @@ include("idrs.jl")
2727
# QMR
2828
include("qmr.jl")
2929

30+
# Look-Ahead Lanczos
31+
include("lal.jl")
32+
3033
#Chebyshev
3134
include("chebyshev.jl")
3235

0 commit comments

Comments
 (0)