We implement the action of the matrices and matrix-free, using LinearMaps.jl:
using LinearMaps
function myA!(y, x)
for i = 1 : length(x)
@inbounds y[i] = sqrt(i) * x[i]
end
end
function myB!(y, x)
for i = 1 : length(x)
@inbounds y[i] = x[i] / sqrt(i)
end
end
n = 100_000
A = LinearMap{Complex128}(myA!, n; ismutating = true)
B = LinearMap{Complex128}(myB!, n; ismutating = true)
The order of the matrices is 100_000
. It turns out that if we target eigenvalues in the interior of the spectrum, iterative solvers used internally in Jacobi-Davidson might have trouble solving very indefinite systems.
In that case we should use a preconditioner for , where is the target. We will just use the exact inverse, which is a diagonal matrix with entries . It can be implemented matrix-free and in-place:
import Base.LinAlg.A_ldiv_B!
struct SuperPreconditioner{numT <: Number}
target::numT
end
function A_ldiv_B!(p::SuperPreconditioner, x)
for i = 1 : length(x)
@inbounds x[i] *= sqrt(i) / (i - p.target)
end
end
Now we call Jacobi-Davidson with the Near
target and pass the preconditioner. We use GMRES as the iterative solver, but we could have used BiCGStabl(l) as well.
using JacobiDavidson
τ = 50_000.1 + 0im
target = Near(τ)
P = SuperPreconditioner(τ)
schur, residuals = jdqz(A, B,
gmres_solver(n, iterations = 10),
preconditioner = P,
target = target,
pairs = 5,
ɛ = 1e-9,
min_dimension = 5,
max_dimension = 10,
max_iter = 200,
verbose = true
)
It converges to the eigenvalues 49999, 50000, 50001, 50002 and 50004:
50004.00000000014 + 3.5749921718300463e-12im
49999.999999986496 - 7.348301591250897e-12im
50001.00000000359 - 1.9761169705101647e-11im
49998.99999999998 - 1.0866253642291695e-10im
50002.00000000171 - 2.3559720511618024e-11im
It does not yet detect 50003, but that might happen when pairs
is increased a bit. As a result of our preconditioner, Jacobi-Davidson converges very quickly:
It's not easy to construct a preconditioner this good for any given problem, but usually people tend to know what works well in specific classes of problems. If no specific preconditioner is availabe, you can always try a general one such as ILU. The next section illustrates that.
As an example of how ILU can be used we generate a non-symmetric, banded matrix having a structure that typically arises in finite differences schemes of three-dimensional problems:
n = 64
N = n^3
A = spdiagm((fill(-1.0, n - 1), fill(3.0, n), fill(-2.0, n - 1)), (-1, 0, 1))
Id = speye(n)
A = kron(A, Id) + kron(Id, A)
A = kron(A, Id) + kron(Id, A)
x = ones(N)
b = A * x
The matrix has size . We want to solve the problem using for instance BiCGStab(2), but it turns out that convergence can get slow when the size of the problem grows. A quick benchmark shows it takes about 2.0 seconds to solve the problem to a reasonable tolerance:
> using BenchmarkTools, IterativeSolvers
> my_x = @btime bicgstabl($A, $b, 2, max_mv_products = 2000);
2.051 s
> norm(b - A * my_x) / norm(b)
1.6967043606691152e-9
Now let's construct the ILU factorization:
> using ILU
> LU = crout_ilu(A, τ = 0.1)
> nnz(LU) / nnz(A)
2.1180353639352374
Using the above drop tolerance , our ILU factorization stores only about twice as many entries as the original matrix, which is reasonable. Let's see what happens when we benchmark the solver again, now with ILU as a preconditioner:
> my_x = @btime bicgstabl($A, $b, 2, Pl = $LU, max_mv_products = 2000);
692.187 ms
> norm(b - A * my_x) / norm(b)
2.133397068536056e-9
It solves the problem 66% faster to the same tolerance. There is of course a caveat, as constructing the preconditioner itself takes time as well:
> LU = @btime crout_ilu($A, τ = 0.1);
611.019 ms
So all in all the problem is solved about 36% faster. However, if we have multiple right-hand sides for the same matrix, we can construct the preconditioner only once and use it multiple times. Even when the matrix changes slightly you could reuse the ILU factorization. The latter is exactly what happens in Jacobi-Davidson.
I would really want to thank my GSoC mentor Andreas Noack for the many discussions we had in chat and video calls. Also I would like to thank the Julia community in general for giving me a warm welcome, both online and at JuliaCon 2017.