Advanced LinSolver usage for deflated NEPs

[jarlebring]

A common use-case: We have a linear solver (eg gmres) for a nep, and want to carry out deflation. The original linear solver will not directly work for the deflated nep due to the way it is modified. However, it is possible to define your own linear solver for the deflated nep. Since this is a common use-case, I suggest we incorporate this into the package in some way.

Illustration (only works for one deflated pair):

using NonlinearEigenproblems
using LinearAlgebra
import NonlinearEigenproblems.create_linsolver
import NonlinearEigenproblems.lin_solve
mutable struct MyLinSolverCreator <: LinSolverCreator
    orglinsolvercreator
    dnep
end

mutable struct MyLinSolver <: LinSolver
    orglinsolver
    dnep
    λ
end


function create_linsolver(creator::MyLinSolverCreator,dnep,λ)
    orglinsolver=create_linsolver(creator.orglinsolvercreator,
                                  dnep.orgnep,λ);
    return MyLinSolver(orglinsolver,dnep,λ);
end

function lin_solve(solver::MyLinSolver, b::AbstractVecOrMat;tol=0)

    n0=size(solver.dnep.orgnep,1);
    b1=b[1:n0];
    b2=b[(n0+1):end];
    z1=lin_solve(solver.orglinsolver,b1);
    z2=b2;

    # Now use Schur complement, i.e., that 
    # inv([I x1 ; x2' 0])=
    # (1/α)*[(α*I+x1*x2') -x1 ; -x2' 1]
    # where x1=x/(λ-s)
    #       x2=x

    s=solver.dnep.S0[1,1];
    x=solver.dnep.V0[:,1];

    x1=x/(solver.λ-s)
    x2=x;
    α=(-x2'*x1)[1];
    α=(-x2'*x1)[1];
    q1=(α*z1+x1*(x2'*z1)[1]-x1*z2[1])/α
    q2=(-x2'*z1+z2[1])/α
    return [q1;q2];

end


nep=nep_gallery("dep0");
(λ,v)=newton(nep,v=ones(size(nep,1)));
dnep=deflate_eigpair(nep,λ,v)

# The underlying linsolver:
orglinsolver=BackslashLinSolverCreator();
creator=MyLinSolverCreator(orglinsolver,dnep);
(λ2,v2)=augnewton(dnep,
                  v=ones(size(dnep,1)),
                  linsolvercreator=creator,
                  logger=1);  # this converges to different eigval

Note that the MyLinSolver is a linear solver for the dnep:

julia> zz=0.3;
julia> linsolver=create_linsolver(creator,dnep,zz);
julia> b=randn(size(dnep,1));
julia> @show z_a=lin_solve(linsolver,b)
z_a = lin_solve(linsolver, b) = Complex{Float64}[-0.4101211405863793 - 0.0im, 1.5263830696452803 - 0.0im, 1.6157867977440057 - 0.0im, -1.1376264798061828 - 0.0im, 0.5353383364146362 - 0.0im, -0.5880736961493719 - 0.0im]
6-element Array{Complex{Float64},1}:
 -0.4101211405863793 - 0.0im
  1.5263830696452803 - 0.0im
  1.6157867977440057 - 0.0im
 -1.1376264798061828 - 0.0im
  0.5353383364146362 - 0.0im
 -0.5880736961493719 - 0.0im
julia> @show z_b=compute_Mder(dnep,zz)\b
z_b = compute_Mder(dnep, zz) \ b = Complex{Float64}[-0.4101211405863793 - 0.0im, 1.5263830696452805 + 0.0im, 1.6157867977440055 + 0.0im, -1.1376264798061826 - 0.0im, 0.5353383364146361 + 0.0im, -0.5880736961493715 - 0.0im]
6-element Array{Complex{Float64},1}:
 -0.4101211405863793 - 0.0im
  1.5263830696452805 + 0.0im
  1.6157867977440055 + 0.0im
 -1.1376264798061826 - 0.0im
  0.5353383364146361 + 0.0im
 -0.5880736961493715 - 0.0im

[eringh]

We might want to do it from scratch to get a good design. However, I think part of what we want is available in the function jd_inner_effenberger_linear_solver!, located at the end of method_jd.jl (here). It can handle a few deflated eigenvalues.