Linear solvers

struct FactorizeLinSolver <: LinSolver

This represents the linear solver associated with julia factorize(). See LinSolver and FactorizeLinSolverCreator for examples.

FactorizeLinSolverCreator(;umfpack_refinements,max_factorizations,nep,precomp_values)

FactorizeLinSolverCreator-objects can instantiate FactorizeLinSolver objects via the create_linsolver function.

The FactorizeLinSolver is based on factorize-calls. The time point of the call to factorize can be controlled by parameters to FactorizeLinSolverCreator:

• By default, the factorize call is carried out by the instantiation of the FactorizeLinSolver, i.e., when the NEP-solver calls create_linsolver.

• You can also precompute the factorization, at the time point when you instantiate FactorizeLinSolverCreator. If you set precomp_values::Vector{Number} to a non-empty vector, and set nep kwarg, the factorization (of all λ-values in the precomp_values) will be computed when the FactorizeLinSolverCreator is instantiated. If the NEP-solver calls a create_linsolver with a λ-value from that vector, the factorization will be used (otherwise it will be computed).

Further recycling is possible. If the variable max_factorizations is set to a positive value, the object will store that many factorizations for possible reuse. Every lin_solve-call then computes a factorization, unless a lin_solve-call for that λ has been computed earlier. This procedure can at most store max_factorization (which can be set Inf).

See also: FactorizeLinSolver, create_linsolver

struct BackslashLinSolver <: LinSolver

This represents a linear solver corresponding to the backslash operator (no pre-factorization).

See also: LinSolver and BackslashLinSolverCreator

struct BackslashLinSolverCreator <: LinSolverCreator

Creator to for the BackslashLinSolver, i.e., usage of backslash to make linear solves. Specify objects of this type if you want the solver to use backslash.

See also: BackslashLinSolver, create_linsolver

struct GMRESLinSolver <: LinSolver

This represents a solver done with the julia GMRES implementation.

See also: LinSolver, GMRESLinSolverCreator

GMRESLinSolverCreator(;kwargs...)

This is the creator for the GMRES-method. Instantiate this object if you want to use GMRES as your linear system solver. The kwargs are stored and used as keyword arguments in the call to gmres. See list of keyword in the IterativeSolvers.jl manual.

DefaultLinSolverCreator

This is the default linear solver if no other is specified (for most methods). It is a FactorizeLinSolverCreator.

See also: LinSolver, create_linsolver, lin_solve, FactorizeLinSolverCreator, FactorizeLinSolver

abstract type LinSolver

Structs inheriting from this type are able to solve linear systems associated with a NEP, for a specific λ-value. The most common are direct solvers such as FactorizeLinSolver, BackslashLinSolver and iterative solvers such as GMRESLinSolver.

The LinSolver objects are usually created by the NEP-algorithms through creator functions, which are passed as parameters.

Example

The most common usecase is that you want to pass a linsolvercreator-function as a parameter to the NEP-algorithm. This example shows how you can use solvers based on backslash or factorize(). In the example, BackslashLinSolver does not exploit that the system matrix remains the same throughout the algorithm and is therefore slower.

julia> nep=nep_gallery("qdep0");
julia> using BenchmarkTools
julia> v0=ones(size(nep,1));
julia> @btime λ,v=quasinewton(nep,λ=-1,v=v0, linsolvercreator=DefaultLinSolverCreator());
  181.017 ms (3779 allocations: 58.61 MiB)
julia> @btime λ,v=quasinewton(nep,λ=-1,v=v0, linsolvercreator=BackslashLinSolverCreator());
  2.040 s (4510 allocations: 553.24 MiB)

Example

The LinSolvers are constructed for extendability. This example creates our own LinSolver which uses an explicit formula for the inverse if the NEP has dimension 2x2.

Create the types and a creator.

julia> using LinearAlgebra
julia> struct MyLinSolver <: LinSolver
   M::Matrix{ComplexF64}
end
julia> function my_linsolvercreator(nep,λ)
   M=compute_Mder(nep,λ);
   return MyLinSolver(M);
end

Explicit import lin_solve to show how to solve a linear system.

julia> import NonlinearEigenproblems.LinSolvers.lin_solve;
julia> function lin_solve(solver::MyLinSolver,b::AbstractVecOrMat;tol=0)
   M=solver.M;
   invM=(1/(det(M)))*[M[2,2] -M[1,2];-M[2,1] M[1,1]]
   return invM*b
end
julia> nep=SPMF_NEP([[1.0 3.0; 4.0 5.0], [2.0 1.0; -1 2.0]], [S->S^2,S->exp(S)])
julia> λ,v=quasinewton(nep,λ=-1,v=[1;1],linsolvercreator=my_linsolvercreator);

See also: lin_solve, FactorizeLinSolver, FactorizeLinSolver, DefaultLinSolverCreator, BackslashLinSolver, BackslashLinSolverCreator, GMRESLinSolver, GMRESLinSolverCreator

lin_solve(solver::LinSolver, b::AbstractVecOrMat; tol=0)

This function solves the linear system represented in solver::LinSolver with a right-hand side b. The tol kwarg is controlling how accurate the linear system needs to be solved. A NEP-algorithm will call this solver every time a linear system associated with M(λ) needs to be solved.

This function must be overloaded if a user wants to define their own way of solving linear systems. See LinSolver for examples.

create_linsolver(creator::LinSovlerCreator,nep,λ)

Creates a LinSolver instance for the nep corresponding which is evaluated in λ. The type of the output is decided by dispatch and the type of the LinSolverCreator.

See also: LinSolver, FactorizeLinSolverCreator, BackslashLinSolverCreator, DefaultLinSolverCreator, GMRESLinSolverCreator.

struct DefaultEigSolver <: EigSolver

A linear eigenvalueproblem solver that calls checks for sparsity and accordingly assigns an appropriate solver.

See also: EigSolver, eig_solve, EigenEigSolver, ArnoldiEigSolver

struct EigenEigSolver <: EigSolver

A linear eigenvalueproblem solver that calls Julia's in-built eigen()

Constructed as EigenEigSolver(A, [B,]), and solves the problem

The paramter B is optional an default is indentity, for which a standard linear eigenproblem is solved.

See also: EigSolver and eig_solve

struct ArnoldiEigSolver <: EigSolver

A linear eigenproblem solver for large and sparse problems that calls the Arnoldi method implemented in the Julia package ArnoldiMethod.jl.

Constructed as ArnoldiEigSolver(A, [B,]), and solves the problem

The paramter B is optional an default is indentity, for which a standard linear eigenproblem is solved.

See also: EigSolver and eig_solve

abstract type EigSolver

Structs inheriting from this type are able to solve linear eigenvalue problems arising in certain methods, such as, e.g., mslp, sgiter, and polyeig.

The EigSolver objects are passed as types to the NEP-algorithms, which uses it to dispatch the correct version of the function eig_solve.

Example

The most common usecase is that you do not want to specify anything in particular, since the DefaultEigSolver will use a dense or a sparse method depending on you problem. However, this example shows how you can force mslp to use the sparse solver.

julia> nep=nep_gallery("qdep0");
julia> λ,v = mslp(nep, eigsolvertype=ArnoldiEigSolver);
julia> norm(compute_Mlincomb(nep,λ,v))
9.323110647141726e-16

Example

The EigSolvers are constructed for extendability. As an illustartion this example creates a naive EigSolver which casts the problem to a standard linear eigenproblem and calls the built-in function to solve it.

Create the types and a creator.

julia> struct MyEigSolver <: EigSolver
   A
   E
   function MyEigSolver(A,E)
      return new(A,E)
   end
end

julia> import NonlinearEigenproblems.LinSolvers.eig_solve;
julia> function eig_solve(solver::MyEigSolver;nev = 1, target = 0)
   M = solver.E \ solver.A
   eig = eigen(M)
   i = argmin(abs.(eig.values))
   return eig.values[i], eig.vectors[:,i]
end
julia> nep=nep_gallery("dep0", 50);
julia> λ,v = mslp(nep, eigsolvertype=MyEigSolver, tol=1e-5);
julia> norm(compute_Mlincomb(nep,λ,v))
3.0777795031319117e-10

See also: eig_solve, DefaultEigSolver, EigenEigSolver, ArnoldiEigSolver, eig_solve

eig_solve(solver::EigSolver; [nev,] [target,])

This function solves the linear eigenvalue problem represented in solver::EigSolver. The nev kwarg is controlling the number of eigenvalues aimed for, and target specifies around which point the eigenvalues are computed. The former has a defalut value equalt to the seize of the problem, and the latter has a defalut value 0.

Return values are of the form (Vector, Matrix) where the former contains the eigenvalues and the latter the eigenvectors.

This function must be overloaded if a user wants to define their own way of solving linear eigenvalue problems. See EigSolver for examples.