Skip to content

Commit 22c5745

Browse files
author
Pawel Latawiec
committed
LALQMR for complex numbers
1 parent a7bfb0d commit 22c5745

File tree

5 files changed

+300
-1
lines changed

5 files changed

+300
-1
lines changed

‎src/IterativeSolvers.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@ include("gmres.jl")
2727
include("chebyshev.jl")
2828
include("idrs.jl")
2929
include("qmr.jl")
30+
include("lalqmr.jl")
3031

3132
# Eigensolvers
3233
include("simple.jl")

‎src/lalqmr.jl‎

Lines changed: 203 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,203 @@
1+
import Base: iterate
2+
3+
export lalqmr, lalqmr!
4+
5+
mutable struct LALQMRIterable{T, xT, rT, lanczosT}
6+
x::xT
7+
8+
lanczos::lanczosT
9+
resnorm::rT
10+
tol::rT
11+
maxiter::Int
12+
13+
t::Vector{T}
14+
R::UpperTriangular{T, Matrix{T}}
15+
16+
G::Vector{Givens{T}}
17+
18+
d::Matrix{T}
19+
end
20+
21+
function lalqmr_iterable!(
22+
x, A, b;
23+
abstol::Real = zero(real(eltype(b))),
24+
reltol::Real = sqrt(eps(real(eltype(b)))),
25+
maxiter::Int = size(A, 2),
26+
initially_zero::Bool = false,
27+
kwargs...
28+
)
29+
T = eltype(x)
30+
r = copy(b)
31+
if !initially_zero
32+
r -= A*x
33+
end
34+
resnorm = norm(r)
35+
v = r/resnorm
36+
w = copy(v)
37+
38+
lanczos = LookAheadLanczosDecomp(
39+
A, v, w;
40+
vw_normalized = true,
41+
kwargs...
42+
)
43+
44+
R = UpperTriangular(Matrix{T}(undef, 0, 0))
45+
# Givens rotations
46+
G = Vector{Givens{T}}()
47+
48+
# Projection operators
49+
D = similar(x, size(x, 1), 0)
50+
51+
t = Vector{T}(undef, 1)
52+
t[1] = resnorm
53+
54+
tolerance = max(reltol * resnorm, abstol)
55+
56+
return LALQMRIterable(
57+
x,
58+
lanczos, resnorm, tolerance,
59+
maxiter,
60+
t, R,
61+
G, D
62+
)
63+
end
64+
65+
converged(q::LALQMRIterable) = q.resnorm ≤ q.tol
66+
start(::LALQMRIterable) = 1
67+
done(q::LALQMRIterable, iteration::Int) = iteration > q.maxiter || converged(q)
68+
69+
function iterate(q::LALQMRIterable, n::Int=start(q))
70+
# Check for termination first
71+
if done(q, n)
72+
return nothing
73+
end
74+
75+
iterate(q.lanczos, n)
76+
77+
# Eq. 6.2, update QR factorization of L
78+
Llastcol = q.lanczos.L[:, end]
79+
for g in q.G
80+
Llastcol = g*Llastcol
81+
end
82+
gend, r = givens(Llastcol, n, n+1)
83+
push!(q.G, gend)
84+
Llastcol[end-1] = r
85+
Llastcol[end] = 0
86+
q.R = UpperTriangular([[q.R; fill(0, 1, n-1)] Llastcol[1:end-1]])
87+
@assert q.R[:, end] ≈ Llastcol[1:end-1]
88+
89+
90+
# Eq. 6.2, update t
91+
push!(q.t, 0)
92+
q.t = gend * q.t
93+
94+
# Eq. 6.3, calculate projection
95+
# Dn = [Vn Un^-1 Rn^-1]
96+
# => Dn Rn Un = Vn
97+
RU = q.R*q.lanczos.U
98+
d = q.lanczos.V[:, end-1]
99+
for i in 1:size(RU, 1)-1
100+
axpy!(-RU[i, end], q.d[:, i], d)
101+
end
102+
d = d / RU[end, end]
103+
q.d = [q.d d]
104+
105+
# iterate x_n = x_n-1 + d_n τ_n
106+
axpy!(q.t[end-1], d, q.x)
107+
108+
# Eq. 4.12, Freund 1990
109+
q.resnorm = q.resnorm * abs(gend.s) * sqrt(n+1)/sqrt(n)
110+
111+
return q.resnorm, n + 1
112+
end
113+
114+
"""
115+
lalqmr(A, b; kwargs...) -> x, [history]
116+
117+
Same as [`lalqmr!`](@ref), but allocates a solution vector `x` initialized with zeros.
118+
"""
119+
lalqmr(A, b; kwargs...) = lalqmr!(zerox(A, b), A, b; initially_zero = true, kwargs...)
120+
121+
"""
122+
lalqmr!(x, A, b; kwargs...) -> x, [history]
123+
124+
Solves the problem ``Ax = b`` with the Quasi-Minimal Residual (QMR) method with Look-ahead. See [`LookAheadLanczosDecomp`](@ref).
125+
126+
# Arguments
127+
- `x`: Initial guess, will be updated in-place;
128+
- `A`: linear operator;
129+
- `b`: right-hand side.
130+
131+
## Keywords
132+
133+
- `initally_zero::Bool`: If `true` assumes that `iszero(x)` so that one
134+
matrix-vector product can be saved when computing the initial residual
135+
vector;
136+
- `maxiter::Int = size(A, 2)`: maximum number of iterations;
137+
- `abstol::Real = zero(real(eltype(b)))`,
138+
`reltol::Real = sqrt(eps(real(eltype(b))))`: absolute and relative
139+
tolerance for the stopping condition
140+
`|r_k| ≤ max(reltol * |r_0|, abstol)`, where `r_k = A * x_k - b`
141+
- `log::Bool`: keep track of the residual norm in each iteration;
142+
- `verbose::Bool`: print convergence information during the iteration.
143+
- `max_block_size`: maximum block size during look-ahead process
144+
- `max_memory`: maximum allowed memory for look-ahead process
145+
146+
# Return values
147+
148+
**if `log` is `false`**
149+
150+
- `x`: approximate solution.
151+
152+
**if `log` is `true`**
153+
154+
- `x`: approximate solution;
155+
156+
- `history`: convergence history.
157+
158+
[^Freund1990]:
159+
Freund, W. R., & Nachtigal, N. M. (1990). QMR : for a Quasi-Minimal
160+
Residual Linear Method Systems. (December).
161+
"""
162+
function lalqmr!(
163+
x, A, b;
164+
abstol::Real = zero(real(eltype(b))),
165+
reltol::Real = sqrt(eps(real(eltype(b)))),
166+
maxiter::Int = size(A, 2),
167+
log::Bool = false,
168+
initially_zero::Bool = false,
169+
verbose::Bool = false,
170+
kwargs...
171+
)
172+
history = ConvergenceHistory(partial = !log)
173+
history[:abstol] = abstol
174+
history[:reltol] = reltol
175+
log && reserve!(history, :resnorm, maxiter)
176+
177+
iterable = lalqmr_iterable!(
178+
x, A, b;
179+
abstol=abstol,
180+
reltol=reltol,
181+
maxiter=maxiter,
182+
initially_zero=initially_zero,
183+
log=log,
184+
verbose=verbose
185+
)
186+
187+
verbose && @printf("=== qmr ===\n%4s\t%7s\n","iter","resnorm")
188+
189+
for (iteration, residual) = enumerate(iterable)
190+
if log
191+
nextiter!(history)
192+
push!(history, :resnorm, residual)
193+
end
194+
195+
verbose && @printf("%3d\t%1.2e\n", iteration, residual)
196+
end
197+
198+
verbose && println()
199+
log && setconv(history, converged(iterable))
200+
log && shrink!(history)
201+
202+
log ? (x, history) : x
203+
end

‎src/qmr.jl‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -266,7 +266,6 @@ function qmr!(x, A, b;
266266
abstol::Real = zero(real(eltype(b))),
267267
reltol::Real = sqrt(eps(real(eltype(b)))),
268268
maxiter::Int = size(A, 2),
269-
lookahead::Bool = false,
270269
log::Bool = false,
271270
initially_zero::Bool = false,
272271
verbose::Bool = false)

‎test/lalqmr.jl‎

Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
module TestLALQMR
2+
3+
using IterativeSolvers
4+
using Test
5+
using LinearMaps
6+
using LinearAlgebra
7+
using Random
8+
using SparseArrays
9+
10+
#LALQMR
11+
@testset "LALQMR" begin
12+
13+
Random.seed!(1234321)
14+
n = 10
15+
16+
@testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
17+
A = rand(T, n, n)
18+
b = rand(T, n)
19+
F = lu(A)
20+
reltol = √eps(real(T))
21+
22+
# Test optimality condition: residual should be non-increasing
23+
x, history = lalqmr(A, b, log=true, maxiter=10, reltol=reltol);
24+
@test isa(history, ConvergenceHistory)
25+
@test all(diff(history[:resnorm]) .<= 0.0)
26+
27+
# Left exact preconditioner
28+
#x, history = lalqmr(A, b, Pl=F, maxiter=1, reltol=reltol, log=true)
29+
#@test history.isconverged
30+
#@test norm(F \ (A * x - b)) / norm(b) ≤ reltol
31+
32+
# Right exact preconditioner
33+
#x, history = lalqmr(A, b, Pl=Identity(), Pr=F, maxiter=1, reltol=reltol, log=true)
34+
#@test history.isconverged
35+
#@test norm(A * x - b) / norm(b) ≤ reltol
36+
end
37+
38+
@testset "SparseMatrixCSC{$T, $Ti}" for T in (Float64, ComplexF64), Ti in (Int64, Int32)
39+
A = sprand(T, n, n, 0.5) + I
40+
b = rand(T, n)
41+
F = lu(A)
42+
reltol = √eps(real(T))
43+
44+
# Test optimality condition: residual should be non-increasing
45+
x, history = lalqmr(A, b, log = true, maxiter = 10);
46+
@test all(diff(history[:resnorm]) .<= 0.0)
47+
48+
# Left exact preconditioner
49+
#x, history = lalqmr(A, b, Pl=F, maxiter=1, log=true)
50+
#@test history.isconverged
51+
#@test norm(F \ (A * x - b)) / norm(b) ≤ reltol
52+
53+
# Right exact preconditioner
54+
#x, history = lalqmr(A, b, Pl = Identity(), Pr=F, maxiter=1, log=true)
55+
#@test history.isconverged
56+
#@test norm(A * x - b) / norm(b) ≤ reltol
57+
end
58+
59+
@testset "Linear operator defined as a function" begin
60+
A = LinearMap(cumsum!, 100; ismutating=true)
61+
b = rand(100)
62+
reltol = 1e-5
63+
64+
x = lalqmr(A, b; reltol=reltol, maxiter=2000)
65+
@test norm(A * x - b) / norm(b) ≤ reltol
66+
end
67+
68+
@testset "Termination criterion" begin
69+
for T in (Float32, Float64, ComplexF32, ComplexF64)
70+
A = T[ 2 -1 0
71+
-1 2 -1
72+
0 -1 2]
73+
n = size(A, 2)
74+
b = ones(T, n)
75+
x0 = A \ b
76+
perturbation = 10 * sqrt(eps(real(T))) * T[(-1)^i for i in 1:n]
77+
78+
# If the initial residual is small and a small relative tolerance is used,
79+
# many iterations are necessary
80+
x = x0 + perturbation
81+
initial_residual = norm(A * x - b)
82+
x, ch = lalqmr!(x, A, b, log=true)
83+
@test 2 ≤ niters(ch) ≤ n
84+
85+
# If the initial residual is small and a large absolute tolerance is used,
86+
# no iterations are necessary
87+
x = x0 + perturbation
88+
initial_residual = norm(A * x - b)
89+
x, ch = lalqmr!(x, A, b, abstol=2*initial_residual, reltol=zero(real(T)), log=true)
90+
@test niters(ch) == 0
91+
end
92+
end
93+
end
94+
95+
end # module

‎test/runtests.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@ include("qmr.jl")
3030

3131
# Look-Ahead Lanczos
3232
include("lal.jl")
33+
include("lalqmr.jl")
3334

3435
#Chebyshev
3536
include("chebyshev.jl")

0 commit comments

Comments
 (0)