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+ rng = Random. MersenneTwister (123 )
14+ n = 10
15+
16+ @testset " Matrix{$T }" for T in (Float32, Float64, ComplexF32, ComplexF64)
17+ A = rand (rng, T, n, n)
18+ b = rand (rng, T, n)
19+ F = lu (A)
20+ reltol = √ eps (real (T))
21+
22+ x, history = lalqmr (A, b, log= true , maxiter= 10 , reltol= reltol);
23+ @test isa (history, ConvergenceHistory)
24+
25+ # Left exact preconditioner
26+ # x, history = lalqmr(A, b, Pl=F, maxiter=1, reltol=reltol, log=true)
27+ # @test history.isconverged
28+ # @test norm(F \ (A * x - b)) / norm(b) ≤ reltol
29+
30+ # Right exact preconditioner
31+ # x, history = lalqmr(A, b, Pl=Identity(), Pr=F, maxiter=1, reltol=reltol, log=true)
32+ # @test history.isconverged
33+ # @test norm(A * x - b) / norm(b) ≤ reltol
34+ end
35+
36+ @testset " SparseMatrixCSC{$T , $Ti }" for T in (Float64, ComplexF64), Ti in (Int64, Int32)
37+ A = sprand (rng, T, n, n, 0.5 ) + n * I
38+ b = rand (rng, T, n)
39+ F = lu (A)
40+ reltol = √ eps (real (T))
41+
42+ x, history = lalqmr (A, b, log = true , maxiter = 10 );
43+ @test norm (A * x - b) / norm (b) ≤ reltol
44+
45+ # Left exact preconditioner
46+ # x, history = lalqmr(A, b, Pl=F, maxiter=1, log=true)
47+ # @test history.isconverged
48+ # @test norm(F \ (A * x - b)) / norm(b) ≤ reltol
49+
50+ # Right exact preconditioner
51+ # x, history = lalqmr(A, b, Pl = Identity(), Pr=F, maxiter=1, log=true)
52+ # @test history.isconverged
53+ # @test norm(A * x - b) / norm(b) ≤ reltol
54+ end
55+
56+ @testset " Block Creation {$T }" for T in (Float32, Float64, ComplexF32, ComplexF64)
57+ # Guaranteed to create blocks during Lanczos process
58+ # This satisfies the condition that in the V-W sequence, the first
59+ # iterates are orthogonal: <Av - v<A, v>, Atv - v<At, v>> under transpose inner product
60+ # These do _not_ work for regular qmr
61+ dl = fill (one (T), n- 1 )
62+ du = fill (one (T), n- 1 )
63+ d = fill (one (T), n)
64+ dl[1 ] = - 1
65+ A = Tridiagonal (dl, d, du)
66+ b = fill (zero (T), n)
67+ b[2 ] = 1.0
68+
69+ reltol = √ eps (real (T))
70+
71+ x, history = lalqmr (A, b, log = true )
72+ @test norm (A * x - b) / norm (b) ≤ reltol
73+ end
74+
75+ @testset " Linear operator defined as a function" begin
76+ A = LinearMap (identity, identity, 100 ; ismutating= false )
77+ b = rand (100 )
78+ reltol = 1e-6
79+
80+ x = lalqmr (A, b; reltol= reltol, maxiter= 2000 )
81+ @test norm (A * x - b) / norm (b) ≤ reltol
82+ end
83+
84+ @testset " Termination criterion" begin
85+ for T in (Float32, Float64, ComplexF32, ComplexF64)
86+ A = T[ 2 - 1 0
87+ - 1 2 - 1
88+ 0 - 1 2 ]
89+ n = size (A, 2 )
90+ b = ones (T, n)
91+ x0 = A \ b
92+ perturbation = 10 * sqrt (eps (real (T))) * T[(- 1 )^ i for i in 1 : n]
93+
94+ # If the initial residual is small and a small relative tolerance is used,
95+ # many iterations are necessary
96+ x = x0 + perturbation
97+ initial_residual = norm (A * x - b)
98+ x, ch = lalqmr! (x, A, b, log= true )
99+ @test 2 ≤ niters (ch) ≤ n
100+
101+ # If the initial residual is small and a large absolute tolerance is used,
102+ # no iterations are necessary
103+ x = x0 + perturbation
104+ initial_residual = norm (A * x - b)
105+ x, ch = lalqmr! (x, A, b, abstol= 2 * initial_residual, reltol= zero (real (T)), log= true )
106+ @test niters (ch) == 0
107+ end
108+ end
109+ end
110+
111+ end # module
0 commit comments