nep-pack / nonlineareigenproblems.jl Goto Github PK
View Code? Open in Web Editor NEWNonlinear eigenvalue problems in Julia: Iterative methods and benchmarks
Home Page: https://nep-pack.github.io/NonlinearEigenproblems.jl
License: MIT License
Nonlinear eigenvalue problems in Julia: Iterative methods and benchmarks
Home Page: https://nep-pack.github.io/NonlinearEigenproblems.jl
License: MIT License
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?
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:
function f{T}(x::T)
, has been changed to function f(x::T) where {T}
.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)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))
eye()
has been deprecated. Instead use Matrix(...)
(not sure what exactly the recommended way is.)workspace()
is removed (don't know how to do instead). Partial replaced by revise
?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.expm
is replaced by exp
. Same for sqrtm
-> sqrt
. I fear this will cause very hidden bugs.type
need to be replaced with struct
or mutable struct
full(A)
should be replaced by Matrix(A)
Points 1-3,8,9 can be done without breaking backward compatibility with 0.6, so they are recommended already now.
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
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.
NoConvergenceException
if it does not find the specified number of eigenvaluesaugnewton
resinv
newton
rfi
quasinewton
newtonqr
implicitdet
mslp
By proper function comments I mean:
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
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.
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.
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.
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 defineerrmeasure
as a type and we need to have a functioncompute_error
(or any other name) that takes as inputerrmeasure,s,v
where(s,v)
is an eigenpair approximation andcompute_error
may have different values such asdefault
(residual) and the functioncompute_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 typeerrmeasure
and overload the functioncompute_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
.
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
and we should and inner solver helper functions corresponding to a couple of other methods, at least
(I don't think we need to include all methods we have as inner solves.)
Migrated from previous repo
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.
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}
.
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?
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%
──────────────────────────────────────────────────────────────────────────
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.
The current way of using nlevp-examples is to
GalleryNLEVP.NLEVP_NEP
by calling GalleryNLEVP.nep_gallery
GalleryNLEVP.nlevp_make_native
which creates a SPMF-objectThis 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:
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.
I will submit several issues to avoid conflict with issue tracking at gitr @ kth.
Migrated from previous repo
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:
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?
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.
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
NonlinearEigenproblems
?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 howDifferentialEquations.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 methodsNEPIAR.jl
Methods similar to infinite arnoldi methodNEPJD.jl
Jacobi-Davidson.NEPWEP.jl
optimized waveguidePkg.clone("git://github.com/nep-pack/NonlinearEigenproblems.jl")
)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.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.