diff --git a/src/LinSolvers.jl b/src/LinSolvers.jl index 97936cc52..1a42eb690 100644 --- a/src/LinSolvers.jl +++ b/src/LinSolvers.jl @@ -7,11 +7,15 @@ module LinSolvers using LinearMaps using ArnoldiMethod + using ..NEPTypes:DeflatedNEP + using ..NEPTypes:deflated_nep_compute_Q + # Linear system of equation solvers export LinSolver export FactorizeLinSolver export BackslashLinSolver export GMRESLinSolver + export DeflatedNEPLinSolver export lin_solve # Eigenvalue solvers @@ -185,6 +189,69 @@ See also: [`LinSolver`](@ref), [`GMRESLinSolverCreator`](@ref) +############################################################################## +""" + struct DeflatedNEPLinSolver <: LinSolver +The `DeflatedNEPLinSolver` solves the linear system connected with a deflated NEP. +The extended linear system is solveed with a Schur complement strategy, recycling +the linear solver of the original NEP. such that +```math +[M, U; X^T, 0][v1; v2] = [b1; b2] +``` +NB1: The implementation assumes minimality index = 1. +NB2: The Schur complement is explicitly formed. Hence it is only efficient for a +few deflated eigenvalues. +See also: [`LinSolver`](@ref), [`DeflatedNEPLinSolverCreator`](@ref), +[`deflate_eigpair`](@ref) +# References +* C. Effenberger, Robust successive computation of eigenpairs for nonlinear eigenvalue problems. SIAM J. Matrix Anal. Appl. 34, 3 (2013), pp. 1231-1256. +""" + struct DeflatedNEPLinSolver{T_NEP, T_num, T_LinSolver} <: LinSolver + deflated_nep::T_NEP + λ::T_num + orglinsolver::T_LinSolver + end + + + function DeflatedNEPLinSolver(nep::DeflatedNEP, λ, orglinsolver::LinSolver) + DeflatedNEPLinSolver{typeof(nep), typeof(λ), typeof(orglinsolver)}(nep, λ, orglinsolver) + end + + + function lin_solve(solver::DeflatedNEPLinSolver, b; tol=0) + deflated_nep = solver.deflated_nep + λ = solver.λ + orglinsolver = solver.orglinsolver + orgnep = deflated_nep.orgnep + + X = deflated_nep.V0 + Λ = deflated_nep.S0 + + n = size(orgnep,1) + m = size(Λ,1) + T = eltype(b) + v = zeros(T,n+m) + + v1 = view(v, 1:n) + v2 = view(v, (n+1):(n+m)) + b1 = b[1:n] + b2 = b[(n+1):(n+m)] + U = deflated_nep_compute_Q(deflated_nep, λ, 0) + + # Precompute some reused entities + b1tilde = lin_solve(orglinsolver, b1, tol=tol) # b1tilde = M^{-1}b1 + Z::Matrix{T} = zeros(T,n,m) + for i = 1:m + Z[:,i] = lin_solve(orglinsolver, vec(U[:,i]), tol=tol) # Z = M^{-1}U + end + S = -X'*Z #Schur complement + v2[:] = S\(b2 - X'*b1tilde) + v1[:] = b1tilde - Z*v2 + + return v + end + + ############################################################################## """