nep-pack/nonlineareigenproblems.jl

Nonlinear eigenvalue problems in Julia: Iterative methods and benchmarks

NEP-PACK

A nonlinear eigenvalue problem is the problem to determine a scalar λ and a vector v such that

M(λ)v=0

where M is an nxn-matrix depending on a parameter. This package aims to provide state-of-the-art algorithms to solve this problem, as well as a framework to formulate applications and easy access to benchmark problems. This currently includes (but is not restricted to) Newton-type methods, Subspace methods, Krylov methods, contour integral methods, block methods, companion matrix approaches. Problem transformation techniques such as scaling, shifting, deflating are also natively supported by the package.

How to use it?

On Julia 1.X, install it as a registered package by typing ] add ... at the REPL-prompt:

julia> ]
(v1.0) pkg> add NonlinearEigenproblems

After that, check out "Getting started" in

NEP-PACK online user's guide

or read the preprint: https://arxiv.org/abs/1811.09592

GIT Version

If you want the cutting edge development version and not the latest release, install it with the URL:

julia> ]
(v1.0) pkg> add git://github.com/nep-pack/NonlinearEigenproblems.jl.git

NEP solvers

Features and solvers (see documentation https://nep-pack.github.io/NonlinearEigenproblems.jl/methods/ for further information and references):

• Arnoldi/Krylov type
  • NLEIGS
  • Infinite Arnoldi method: (iar)
  • Tensor infinite Arnoldi method (tiar)
  • Infinite bi-Lanczos (infbilanczos)
  • Infinite Lanczos (ilan)
  • AAA CORK (AAAeigs)
• Projection methods
  • Jacobi-Davidson (jd_effenberger)
  • Jacobi-Davidson (jd_betcke)
  • Nonlinear Arnoldi method (nlar)
  • Common Rayleigh-Ritz projection interface
• Contour integral methods
  • Beyn's contour integral method
  • Block SS (Higher moments) contour integral method of Asakura & Sakurai
  • Common quadrature interface for parallelization
• Newton & Rayleigh type:
  • Classical Newton-Raphson
  • Augmented Newton
  • Residual inverse iteration
  • Quasi-Newton
  • Block Newton
  • Rayleigh functional iteration (RFI a, b)
  • Newton-QR
  • Implicit determinant method
  • Broyden's method
• Companion matrices
  • First companion form
  • Companion form for Chebyshev polynomials
• Other
  • Chebyshev interpolation
  • Transformations
  • Rayleigh-Ritz (ProjNEP and inner_solve)
  • Problem gallery (including access to the NLEVP-gallery)
  • Deflation (Effenberger style)

Development

Core developers (alphabetical): Max Bennedich, Elias Jarlebring (www.math.kth.se/~eliasj), Giampaolo Mele (www.math.kth.se/~gmele), Emil Ringh (www.math.kth.se/~eringh), Parikshit Upadhyaya (https://www.kth.se/profile/pup/). Thanks to A Koskela for involvement in initial version of the software.

How to cite

If you find this software useful please cite

@Misc{,
  author = 	 {E. Jarlebring and M. Bennedich and G. Mele and E. Ringh and P. Upadhyaya},
  title = 	 {{NEP-PACK}: A {Julia} package for nonlinear eigenproblems},
  year = 	 {2018},
  note = 	 {https://github.com/nep-pack},
  eprint = {arXiv:1811.09592},
}

If you use a specific method, please also give credit to the algorithm researcher. Reference to a corresponding algorithm paper can be found by in, e.g., by writing ?resinv.

Some links below are developer info on KTH. We will migrate them soon:

• NEP-page style "guide": https://github.com/nep-pack/NonlinearEigenproblems.jl/wiki/Style-guidelines-and-notes

• GIT-workflow: https://github.com/nep-pack/NonlinearEigenproblems.jl/wiki/Git-workflow

• GIT-usage @ KTH: https://gitr.sys.kth.se/nep-pack/nep-pack-alpha/wiki

• NEP-methods @ KTH: https://gitr.sys.kth.se/nep-pack/nep-pack-alpha/wiki/NEP-methods

• NEP-applications @ KTH: https://gitr.sys.kth.se/nep-pack/nep-pack-alpha/wiki/Applications

nonlineareigenproblems.jl's Issues

Checklist for migration of unit testing

Functionality which needs unit test (migrated from kth github # 28)

run-scripts which should be migrated:

• run_test_comp.jl @upadpup
• run_test_eigenvector_extraction*.jl @gmele
• run_test_TIAR.jl @gmele
• run_test_scaling.jl @eliasj
• run_test_sg.jl @eringh (?)
• run_test_jd_quad.jl @eringh (?)

dummy

dummy

Type stability and arithmetic type

After PR #44 simple algorithms such as augnewton are almost type stable if you also provide the arithmetic type. However, if you don't provide the arithmetic type (e.g. just augnewton(nep) instead of augnewton(Complex128,nep)) it is not type stable. It becomes type-stable if we add the types to the kwargs. It seems to work to instead of

    function augnewton(::Type{T},
                       nep::NEP;
                       errmeasure::Function = default_errmeasure(nep::NEP),
                       tol=eps(real(T))*100,
                       maxit=30,
                       λ=zero(T),
                       v=randn(real(T),size(nep,1)),
                       c=v,
                       displaylevel=0,
                       linsolvercreator::Function=backslash_linsolvercreator,
                       armijo_factor=one(real(T)),
                       armijo_max=5) where {T<:Number}

do

   function augnewton(::Type{T},
                       nep::NEP;
                       errmeasure::Function = default_errmeasure(nep::NEP),
                       tol::Real=eps(real(T))*100,
                       maxit::Integer=30,
                       λ::Number=zero(T),
                       v::Vector=randn(real(T),size(nep,1)),
                       c::Vector=v,
                       displaylevel::Int=0,
                       linsolvercreator::Function=backslash_linsolvercreator,
                       armijo_factor::Real=one(real(T)),
                       armijo_max::Integer=5) where {T<:Number}

Commit for augnewton coming soon.

We should add types to all kwargs in all NEP-solvers. In general and in particular while doing this change, we can try to remember to use Vector or Matrix instead of Array{,1} and Array{,2}.

dummy

EigSolver refactorization for consistency

Migrated from gitr @ kth # 29:

Currently a LinSolver is passed to the method by passing a function linsolvercreator::Function=default_linsolvercreator, whereas the specification of EigSolver is passed as a Type eigsolvertype::DataType=DefaultEigSolver. This is indeed inconsistent as @eringh and @pup wrote:

@eringh in #26:

I think we potentially should redo the EigSolver in the same way as the LinSolver. As it is now a LinSolver comes with a creator function that creates a LinSolver. However, the EigSolver is still in the previous thinking where a DataType is provided, and directly used for creating the solver through the constructor. Maybe the variety of EigSolvers is not that big as for LinSolvers (where for example a preconditioner is much easier to use with the "new" structure). However, is seems a bit confusing to have two different ways of working, especially since they come from the same module. Should we add this the the check list?

@pup wrote in #26:

As @eringh mentioned in one of the previous comments, there are some differences between how we implement the LinSolver and EigSolver types. EigSolver and and its derivatives are not parametric types and the there is no creator function involved. Although this doesn't make much difference to the method user, who will just type in the name of the solver/creator function as an arguement, it is good to have consistency, as Emil mentioned. The only downside I can think of is that a lot of our code assumes the current implementation of EigSolver, and could break, but that of course is fixable. @eliasj , Should I go ahead with changing the EigSolver implementation?

I agree that they are inconsistent, and this seems unnatural. I'm not yet sure what the best solution is.

Note that a LinSolver is conceptually slightly different from a EigSolver since a linear system has a right-hand side. A LinSolver object can be used for different right-hand sides, whereas an EigSolver does not a have a corresponding input (except for arguably the target and the number of eigenvalues, but still quite different).

I think Julia is designed for extensive use of types (and not function pointers), e.g., look at the modern way of doing orthogonalization in IterativeMethods, where one passes the type IterativeSolvers.OrthogonalizationMethod, where the corresponding function is called by

        H[k+1,k] = orthogonalize_and_normalize!(VV, vv, view(H,1:k,k), orthmethod)

This seems like a third approach. I think we need to think about if we can come up with a similar construction for our case, before we do any actual code refactorization.

dummy

Julia packaging of NEP-PACK

Eventually we want the NEP-PACK to become available in the julia packaging ecosystem, e.g., Pkg.add("NEPPACK") or something. Here is a list (to be updated) of things that needs to be done to achieve that. This is mostly based on what is expected here https://docs.julialang.org/en/release-0.6/manual/packages/

• Require and dependencies. Specify julia requirements in REQUIRE
• Naming the repository
  • Name should end with .jl
  • Name should "Avoid jargon". NEP might be jargon, so package name NonlinearEigenproblems?
• Introduce versioning
• Create a main/supermodule. The file src/NonlinearEigenproblems.jl is a special file which is loaded when we do Pkg.add() (I think). It should create a module named NonlinearEigenproblems, which contains necessary info. Not sure how
• Separate the code into several packages. We can use DifferentialEquations.jl https://github.com/JuliaDiffEq/DifferentialEquations.jl as a model. They have a DifferentialEquations.jl-package, which is very small and only does require on many (all?) packages . Most of the common code for differential equations are in the separate package DiffEqBase.jl. Analogous to this, we would have:
  • NonlinearEigenproblems.jl
  • NEPBase.jl essentially what is currently in NEPCore and NEPTypes.
  • NEPGallery.jl
  • NEPNewton.jl Newton-like methods
  • NEPIAR.jl Methods similar to infinite arnoldi method
  • NEPJD.jl Jacobi-Davidson.
  • NEPWEP.jl optimized waveguide
• Register the packages such that it becomes available on https://pkg.julialang.org/. (For later. No rush. We can first make it available through Pkg.clone("git://github.com/nep-pack/NonlinearEigenproblems.jl"))
• Continuous integration #36
• Algorithm points in #26

Note that the julia packaging system will change in a new version of julia. My interpretation of Stefan Karpinski's lecture at JuliaCON '17 (https://www.youtube.com/watch?v=-yUiLCGegJs) about Pkg3 is that a transition to the Pkg3 will probably not be very problematic.

dummy

Speed up unit tests

Unit tests are quite slow, particularly on the build server (Travis), which often causes the build to time out, and won't scale as we add more tests. I've added a performance report at the end of each test run where we can see how long each test takes on the build server, see an example here, or the report here:

               Test Performace                   Time         Allocations  
                                            ──────────────   ──────────────
               Total measured:                   2732s          50.0GiB    
 Section                            ncalls     time   %tot     alloc   %tot
 ──────────────────────────────────────────────────────────────────────────
 infbilanczos                            1     632s  23.1%   12.7GiB  25.5%
 jd                                      1     392s  14.3%   3.16GiB  6.32%
 tiar                                    1     361s  13.2%   2.78GiB  5.56%
 dep_distributed                         1     225s  8.25%   1.06GiB  2.13%
 compute_eigvec_from_eigval_lu           1     168s  6.14%   2.54GiB  5.08%
 iar                                     1     134s  4.89%   1.00GiB  2.01%
 wep_small                               1     119s  4.34%   11.1GiB  22.2%
 beyn                                    1     115s  4.22%    793MiB  1.55%
 transf                                  1     101s  3.69%   2.43GiB  4.86%
 compan                                  1    76.2s  2.79%    936MiB  1.83%
 compute_eigvec_from_eigval_lopcg        1    48.1s  1.76%   1.81GiB  3.63%
 deflation                               1    45.8s  1.68%    676MiB  1.32%
 proj                                    1    45.4s  1.66%   1.12GiB  2.23%
 spmf                                    1    36.0s  1.32%    825MiB  1.61%
 inner_solves                            1    34.5s  1.26%   0.94GiB  1.87%
 sgiter                                  1    33.1s  1.21%    877MiB  1.71%
 rfi                                     1    31.6s  1.16%    800MiB  1.56%
 bigfloats                               1    29.8s  1.09%    826MiB  1.61%
 newton_test                             1    25.2s  0.92%    756MiB  1.48%
 mslp                                    1    23.9s  0.88%    720MiB  1.41%
 periodicdde                             1    22.6s  0.83%   1.35GiB  2.70%
 core                                    1    17.9s  0.66%    411MiB  0.80%
 blocknewton                             1    15.8s  0.58%    522MiB  1.02%
 ──────────────────────────────────────────────────────────────────────────

Test Coverage Reports

I've added test coverage reporting through Coveralls and Codecov. You can click the badges at the top of README.md to see the reports, which are automatically generated on each successful build. From there, you can dig into the coverage for each file, and even see line-by-line what's covered by tests and what's not.

For files and folders that don't need to be unit tested, we can add them to the list of excluded files in this script, so that they don't dilute the coverage percentage. I added a few already, but feel free to add more.

dummy

nep-pack transition to 0.7 (and 1.0)

Here is an issue to collect points that need to be addressed when we switch to julia version 0.7. There is an extensive list here: https://github.com/JuliaLang/julia/blob/master/NEWS.md but we can collect the ones relevant for us here:

1. The syntax for parametric methods, function f{T}(x::T), has been changed to function f(x::T) where {T}.
2. ones(3)+1 will no longer work. Instead one should use ones(3) .+ 1. Note it is recommended to use space before the dot, to avoid ambiguity with the dot operator (as in 3 below)
3. MATLAB-style componentwise version of basic operations is deprecated, sin(ones(3)) will not work. Use sin.(ones(3)) instead. (The dot is actually an operation, so you can use it on functions you defined yourself: f=x-> (x+1)^2; f.(ones(3))
4. eye() has been deprecated. Instead use Matrix(...) (not sure what exactly the recommended way is.)
5. workspace() is removed (don't know how to do instead). Partial replaced by revise?
6. qr methods now return decomposition objects such as QR, QRPivoted, and QRCompactWY rather than tuples of arrays. We need to do a search and replace on all calls to qr.
7. expm is replaced by exp. Same for sqrtm -> sqrt. I fear this will cause very hidden bugs.
8. All type need to be replaced with struct or mutable struct
9. full(A) should be replaced by Matrix(A)
10. "global" definition is necessary to make (some variables) visible inside e.g. a for-loop. See comment below.

Points 1-3,8,9 can be done without breaking backward compatibility with 0.6, so they are recommended already now.

IAR chebyshev bug - scaling of chebyshev polys

The test in 5432caf illustrates a bug in IAR chebyshev, when specifying the interval for the scaling of chebyshev poly's. It has something to do with domain of acos. @meleg has agreed to look into it. The iar_chebyshev test is temporary disabled in runtests.jl for this reason..

dummy

dummy

dummy

dummy

Broydens method in julia

The manuscript "Broyden's method for nonlinear eigenproblems" https://arxiv.org/abs/1802.07322 (by myself) has julia code available online, but not yet incorporated into the nep-pack framework. The code is however written with nep-pack in mind, so incorporation is expected to be modestly complicated.

julia eigs buggy for generalized eigenvalue problems

Migrated from previous repo

dummy

Style guideline - trailing whitespace

The default behavior of Juno/Atom is to remove trailing whitespace in files on save. This can be turned off to preserve existing whitespace. Unless everyone uses the same behavior, we will end up with unnecessarily large diffs and merge conflicts. See for example this changeset.

I have a slight preference for removing trailing whitespace, which is also the guideline used within the Julia Language project. Not doing this, I think that we will 1) accumulate random whitespace over time which will make editing a bit less efficient, and 2) have external contributors remove trailing whitespace in entire files and create large changesets for tiny changes, since that's the default behavior in Juno and Julia.

Opinions?

Problems with Compute_Mlincomb

The pull request #53 seemed ok when we merged it. However, it has given some problems for us now. Typically we have observed that we get stuck in infinite loops, which gives Stack Overflow (I get it on my JD-branch now, or try the state given by the commit ea266fe but without typing a as a Vector). Or, because of the strict typing we get ambiguous calls (see cba5b03).

Perhaps we should not implement the transition functions but require the type developer to do that in order to use the relevant methods.

dummy

Remove MAT-package dependency

dummy

dummy

compute_Mlincomb_from_MM gives Overflow for large k

Migrated from previous repo

dummy

Sum of NEPs

Several NEPs in applications are naturally formed as sums of other NEPs, potentially of different NEP-types. (Quite close to @maxbennedich's recent work in NEPTypes for NLEIGS, who could use a sum of PEP and SPMF)

It can be achieved by

struct SUMNEP{NEP1<:NEP,NEP2<:NEP}  <: NEP 
   nep1::NEP1
   nep2::NEP2
end

which allows us to to delegate interface functions, e.g.,

  compute_Mlincomb(nep::SUMNEP, λ, V) = (compute_Mlincomb(nep1, λ, V)+compute_Mlincomb(nep2, λ, V))

For the sum of two NEPs, this seems natural and seems to be doable in a type-stable way. (Note that the types of nep1 and nep2 is stored in the parameter of SUMNEP.) However, the fact that this is only two NEPs is hard-coded. Sometimes one would like to construct the sum of several NEPs. I don't know how to do that in a type stable way.

Attempt 1: The struct

struct SUMNEP  <: NEP 
   nepv::Vector{<:NEP}
end

allows storing a sum of more than two NEPs, but does not contain the type information in the parameter of SUMNEP (so not good for type stability).
Attempt 2: The definition

struct SUMNEP{COMMONNEP<:NEP}  <: NEP 
   nepv::Vector{COMMONNEP}
end

would store type information, but it would force the NEPs to be of the same type.

Any other ideas how to handle the sum of more than two NEPs in a type stable way?

dummy

dummy

dummy

dummy

Type stability for DEPs

At the meeting 2018-07-17 @eringh proposed to look into type stability. Rather than trying to make NEP-PACK type stable right now (which seems like a huge task) we can try it out for a part of the code, and see how difficult it is and if we get perforance improvement. I think we can first try to make the functions associated with one NEP-type type stable. I suggest the DEP.

The command @code_warntype is useful when investigating type stability. For instance, if you run it on a DEP to check type stability of the size command:

julia> nep=nep_gallery("dep0");
julia> @code_warntype(size(nep,1))
Variables:
  #self# <optimized out>
  nep::NEPTypes.DEP
  dim <optimized out>

Body:
  begin 
      return (Core.getfield)(nep::NEPTypes.DEP, :n)::Int64
  end::Int64

At the end you see that the output is Int64, so julia could determine the output type based on input type => type stable.

In contrast, if we run it on issparse we get

julia> @code_warntype(issparse(nep))
Variables:
  #self# <optimized out>
  nep::NEPTypes.DEP

Body:
  begin 
      return (NEPTypes.issparse)((Base.arrayref)((Core.getfield)(nep::NEPTypes.DEP, :A)::Array{#s739,N} where N where #s739<:(AbstractArray{T,2} where T), 1)::AbstractArray{T,2} where T)::Any
  end::Any

At the end we see Any, which means julia could not determine the output type based on the input types. I think we should try to make issparse and compute_Mder() type stable first (and then look at compute_Mlincomb and maybe some other functions).

Type stable for DEP:

• size()
• issparse()
• compute_Mder
• compute_Mlincomb

This julia discourse discussion will probably be helpful. The size-function was made type-stable in ba09cd7 and 248f227.

dummy

From passing functions to passing types

Many functions that we have implemented have as arguments other functions. Quoting @jarlebring in the issue #32

I think Julia is designed for extensive use of types (and not function pointers)
Therefore we should try to avoid every function-passing and replace it with a type-passing.

For example, quoting myself from the issue #32

We should also change errmeasure that at the moment is a function. We can define errmeasure as a type and we need to have a function compute_error (or any other name) that takes as input errmeasure,s,v where (s,v) is an eigenpair approximation and compute_error may have different values such as default (residual) and the function compute_error may be dispatched to compute differently the error for the different formats of the NEP (DEP,PEP,SPMF,...), and the user can define his own type errmeasure and overload the function compute_error to measure the error in the way that he wants.

I think that this should also be the spirit for passing a linsolver (currently done with linsolvercreator). We should not pass a function but a type. Clearly a pre-computation may be needed in this case. This can be done in the same fashion of the function iar_chebyshev for the argument compute_y0_method that requires a pre-computation.

Therefore I here attach a check list for these changes, namely the following arguments should be types

• errmeasure
• linsolvercreator
• EigSolver

At the moment we agree that errmeasure should be changed. What about the others? Refer to the discussion in the issue #32 for EigSolver.

Native support of nlevp-methods

The current way of using nlevp-examples is to

1. create a GalleryNLEVP.NLEVP_NEP by calling GalleryNLEVP.nep_gallery
2. (optional) call GalleryNLEVP.nlevp_make_native which creates a SPMF-object

This construction is nice since it is general, but it requires an installation of MATLAB. For some (or all) examples in NLEVP we may want to create native support in Gallery.nep_gallery.

This can be done with by saving the files with our Serializer util: src/utils/Serializer.jl. For usage see test/serialization.jl. (Edit removed usage MAT)

I have already done it for "gun"-problem by adding an example "nlevp_native_gun" in Gallery.nep_gallery.

Suggested things to do:

• Nativize "fiber"
• Nativize "hadeler" (convert to DEP?)
• Nativize "pdde_stability" to PEP
• Nativize "railtrack"
• Nativize "loaded_string" (convert to REP?)
• Nativize "cd_player" (PEP)
• Make the NEP-PACK find installed NLEVP files in a general way.

dummy

nleigs as part of NEPSolver like other solvers

At the moment nleigs is treated in a slightly different way than other methods. It is quite a bit of software, so it probably deserves a bit separate treatment, e.g. with its own directory like now. It would however be nice if it would be accessible like the other methods, e.g., by appropriate include-commands in NEPSolver.jl.

This is how I would want to use it considering other solvers:

using Gallery
using NEPSolver
nep = nep_gallery("dep0")
S=[-1.0-1im,-1+1im,+1+1im,1-1im];
nleigs(nep,S)

This does not work at the moment. Even after loading files like in the test, I did not manage to make it work, since it throws an error I do not understand:

julia> nleigs(nep,S)
ERROR: UndefVarError: NleigsNEP not defined
Stacktrace:
 [1] get_nleigs_nep(::Type{Float64}, ::NEPTypes.DEP{Array{Float64,2}}) at /home/jarl/.julia/v0.6/NonlinearEigenproblems/src/nleigs/method_nleigs.jl:386
 [2] #nleigs#99(::Array{Float64,1}, ::Int64, ::Int64, ::Int64, ::Int64, ::Float64, ::Float64, ::Array{Complex{Float64},1}, ::NEPCore.##9#10{NEPTypes.DEP{Array{Float64,2}}}, ::Bool, ::Bool, ::Int64, ::Array{Complex{Float64},1}, ::Int64, ::Int64, ::Bool, ::Int64, ::#nleigs, ::Type{Float64}, ::NEPTypes.DEP{Array{Float64,2}}, ::Array{Complex{Float64},1}) at /home/jarl/.julia/v0.6/NonlinearEigenproblems/src/nleigs/method_nleigs.jl:77
 [3] nleigs(::Type{Float64}, ::NEPTypes.DEP{Array{Float64,2}}, ::Array{Complex{Float64},1}) at /home/jarl/.julia/v0.6/NonlinearEigenproblems/src/nleigs/method_nleigs.jl:71
 [4] #nleigs#98 at /home/jarl/.julia/v0.6/NonlinearEigenproblems/src/nleigs/method_nleigs.jl:47 [inlined]
 [5] nleigs(::NEPTypes.DEP{Array{Float64,2}}, ::Array{Complex{Float64},1}) at /home/jarl/.julia/v0.6/NonlinearEigenproblems/src/nleigs/method_nleigs.jl:47

dummy

dummy

I will submit several issues to avoid conflict with issue tracking at gitr @ kth.

dummy

dummy

First publicly announced version

This issue is checklist/wishlist of things we should plan to do before we feel we can start announcing the existance of the packade to colleagues (maybe version 0.1). Feel free to add if you think it should be included.

• NLArnoldi
  • Modernization of orthogonalization
  • Stay real when problem is real
  • One version of restarting
  • Proper function comments
  • Unit testing
  • Throw NoConvergenceException if it does not find the specified number of eigenvalues
  • Selection around disk
• IAR etc
  • One version of restart
  • Chebyshev version (as different function maybe)
  • Extraction by projection (instead of Hessenberg)
  • Proper function comments
  • Proper function comments infbilanczos
• TIAR
  • Extraction by projection
  • Modernize orthogonalization
  • Proper function comments
• SG
  • Reimplementation of such that we can compute kth eigenvalue (for minmax problems)
  • Proper function comments
• JD
  • Reimplementation with new ProjectedNEP
  • Reimplementation with modern orthogonalization
  • Reimplement the linear solve to get rid of explicit matrix access (computed via Mlincomb)
  • Implement so that it can be used to compute more than one eigenvalue
  • Proper function comments
  • Deflation a.l.a effenberger(?)
• Newton methods
  • Proper function comments augnewton
  • Proper function comments resinv
  • Proper function comments newton
  • Proper function comments rfi
  • Proper function comments quasinewton
  • Proper function comments newtonqr
  • Proper function comments implicitdet
  • Proper function comments mslp
• Block Newton method
  • Basic implementation
  • Proper function comments
• WEP
  • Direct LinSolver on Schur complement
• Contour integral methods
  • Beyns version
    • trapezoidal rule
    • Proper function comments
    • Eigenvector extraction
  • SLEPc version
    • Basic implementation
    • Proper function comments
• Transformations
  • Effenberger deflation
• Core and types
  • REP as SPMFs
  • Isolate Matlab dependence from main package (LinSolver)
  • Proper function comments for
    • compute_Mder, compute_Mlincomb, compute_MM
    • compute_rf
    • SPMF,PEP, REPs
    • Projectable NEPs
    • nep transformations
    • linear system solvers
    • linear eigenvalue problem solvers
    • companion
• Practical
  • Introduce a copyright/copying licence
  • Proper first page
  • A short/technical manuscript which we can put on arxiv such that it can be cited.

By proper function comments I mean:

• Function definition with optional parameters
• Description of function (reference newton() for the common parameters)
• A self-contained example.
• References

I have already done "proper function comments" for newton() and resinv(). We should do it using Markdown, i.e., the rules here https://docs.julialang.org/en/stable/manual/documentation/

Migrated from https://gitr.sys.kth.se/nep-pack/nep-pack-alpha/issues/26

Inner solves and projected problems

Several of the methods are inner-outer iteration methods. At some point in the algorithm you need to solve a smaller problem (typically the Galerkin projection of the problem). In commit 6554628 and the following commits, I have made an attempt to make this systematic. The inner solves are controlled by the helper functions in inner_solver.jl. I have also already included it in iar, as an alternative way to extract eigenpairs complement to the computation of eigenvalues from the hessenberg matrix.

Consistent with IAR we now need to incorporate this system into

• Jacobi-Davidson
• NLArnoldi
• TIAR

and we should and inner solver helper functions corresponding to a couple of other methods, at least

• sgiter
• iar (yes, thereby you can do inner and outer iterations with IAR)
• augnewton
• polyeig
• beyncontour

(I don't think we need to include all methods we have as inner solves.)

dummy

Unit testing in github + continuous integration

According to the julia docs
https://docs.julialang.org/en/release-0.6/manual/packages/#Generating-the-package-1
it should be possible to make the unit-tests run on a continuous integration system (Travis and/or AppVeyor) reasonably easy, i.e. on a separate server which carries out tests after each commit. I think we should start using this if we can get it for free. Fixing this issue requires more understanding of continuous integration systems and github than actual programming.

dummy

Julia and eigs in NEP-PACK

As we have observed several times e.g. #1, julias eigs is broken for a large class of generalized eigenvalue problems. As pointed out here https://github.com/JuliaLang/julia/issues/ 24668 this seems to be restricted to indefinite B-matrices, without a clear solution in sight. Our current workaround with MATLAB is not so nice for many reasons. We should try to come up with better work-arounds.