NLEIGS TOAR

[jarlebring]

Here is a reimplementation of NLEIGS TOAR. This is inspired / follows the implementation in SLEPc: https://slepc.upv.es/slepc-master/src/nep/impls/nleigs/. In contrast to our other implementation of NLEIGS, this incorporates a compact representation of the basis, ie the same basic idea as we already have in tiar. I have kept it separate from our other nleigs implementation so we can compare with SLEPc.

Our implementation is separated into an approximate function (leja bagby approximation) NLEIGS_NEPand the actual rational Krylov solver nleigs_toar. The approximation function returns a new NEP-object which we send to the nleigs_toar-function. The nleigs_toar can only solve NEPs of the type NLEIGS_NEP, consistent with the fact that this is a fixed approximation method.

Try it out like this:

julia> import Base.exp
julia> using LinearAlgebra
julia> exp(A::LowerTriangular)=exp(Matrix(A)); # Only needed for the delay example. Matrix exponential for triangular matrix required and not available in native julia: see https://github.com/JuliaLang/julia/issues/32899#issuecomment-521432236
julia> using NonlinearEigenproblems
julia> using Random
julia> Random.seed!(0);
julia> n=50;
julia> A0=randn(n,n)
julia> A1=randn(n,n);
julia> tau=0.5
julia> nep=DEP([A0,A1],[0,tau]);
julia> pshifts=[-1.3+0.2im ] # shifts
julia> prg=[-20,5,-0.001,0.001]; # region
julia> ptarget=-1.1+0.00001im
julia> nleigs_nep=NLEIGS_NEP(nep,prg,shifts=pshifts);  ## Create an approximate NEP
julia> neigs=1
julia> (λv,VV)=nleigs_toar(nleigs_nep,prg,neigs=neigs,target=ptarget,logger=0); ## Call the solver
julia>  minimum(svdvals(compute_Mder(nep,λv[1])))
4.7247867700526746e-7

It is not ready to be merged. I put it here for comments, and progress documentation.

TODO:

• Extract eigenvectors
• Improve domain discretization
• Move RQ-restart to separate function
• Improve breakdown in toar_extendbasis and nleigs_toar_run
• More efficient CheapKrylovErrmeasure
• unit tests
• multiple shifts in rkcontinutation