Jacobi-Davidson

An implementation of Jacobi-Davidson in Julia. This method can be used as an alternative to Julia's eigs to find a few eigenvalues and eigenvectors of a large sparse matrix.

Example

We generate two random complex matrices A and B and use JDQZ to find the eigenvalues λ of the generalized eigenvalue problem Ax = λBx.

using JacobiDavidson, LinearAlgebra, SparseArrays, Plots

function run(n = 1000)
  A = 2I + sprand(ComplexF64, n, n, 1 / n)
  B = 2I + sprand(ComplexF64, n, n, 1 / n)

  # Find all eigenvalues with a direct method
  values = eigvals(Matrix(A), Matrix(B))

  target = Near(1.5 - 0.7im)

  pschur, residuals = jdqz(
    A, B,
    solver = GMRES(n, iterations = 7),
    target = target,
    pairs = 7,
    subspace_dimensions = 10:15,
    max_iter = 300,
    verbosity = 1
  )

  # The eigenvalues found by Jacobi-Davidson
  found = pschur.alphas ./ pschur.betas

  # 
  p1 = scatter(real(values), imag(values), label = "eig")
  scatter!(real(found), imag(found), marker = :+, label = "jdqz")
  scatter!([real(target.τ)], [imag(target.τ)], marker = :star, label = "Target")

  p2 = plot(residuals, yscale = :log10, marker = :auto, label = "Residual norm")

  p1, p2
end

The first plot shows the full spectrum, together with the target we have set and the seven converged eigenvalues:

The second plot shows the residual norm of Ax - λBx during the iterations: