add DeflatedNEPLinSolver again

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
+
+
 
 ##############################################################################
 """