no method matching ( ::DeflationProblem, SArray ... ) for static immutable arrays

[gszep]

Deflated newton works with SizedArray and MArray after patching LinearAlgebra.factorize however DeflationProblem fails to initialise when using SArray. SArray still works with the normal newton without deflation though. Here a minimal working example

using BifurcationKit,Parameters,StaticArrays,ForwardDiff
using BifurcationKit: DeflationOperator
using LinearAlgebra: dot

########################################################### setup
function F(x,params) where T<:Number
	@unpack k = params
	f = first(x)*first(k)
	F = similar(x,typeof(f))

	F[1] = k[2] + x[1] - x[1]^(k[1]+1)/(k[1]+1)
	F[2] = x[1] - x[2]
	return F
end

J(x,p) = ForwardDiff.jacobian( x->F(x,p), x )
p = ( k=[1.0,1.0], )

using LinearAlgebra: lu
import LinearAlgebra: factorize
factorize(A::StaticMatrix) = lu(A) # default factorization for StaticArrays

########################################################### normal arrays
x = [-1.0,42.0]
deflated_root = [Inf,Inf]
deflation = DeflationOperator(1.0, dot, 1.0, [ deflated_root ] )

u, = newton( F, J, x, p, NewtonPar() )                      # success
u, = newton( F, J, x, p, NewtonPar(), deflation )      # success

########################################################### static immutable arrays
x = @SVector[-1.0,42.0]
deflated_root = @SVector[Inf,Inf]
deflation = DeflationOperator(1.0, dot, 1.0, [ deflated_root ] )

u, = newton( F, J, x, p, NewtonPar() )                      # success
u, = newton( F, J, x, p, NewtonPar(), deflation )      # no method matching ( ::DeflationProblem, ::SArray, ... )

########################################################### static mutable arrays
x = @MVector[-1.0,42.0]
deflated_root = @MVector[Inf,Inf]
deflation = DeflationOperator(1.0, dot, 1.0, [ deflated_root ] )

u, = newton( F, J, x, p, NewtonPar() )                      # success
u, = newton( F, J, x, p, NewtonPar(), deflation )      # success

########################################################### sized mutable arrays
x = SizedVector{2}([-1.0,42.0])
deflated_root = SizedVector{2}([Inf,Inf])
deflation = DeflationOperator(1.0, dot, 1.0, [ deflated_root ] )

u, = newton( F, J, x, p, NewtonPar() )                      # success
u, = newton( F, J, x, p, NewtonPar(), deflation )      # success

[rveltz]

Hi,

Note that you may be interested in using DefaultLS(false) to remove the need to define factorize.

One thing I noticed (maybe unimportant) is that F returns a MArray when given a SArray.

I cannot even make

########################################################### normal arrays
x = @SVector[-1.0,42.0]
deflated_root = @SVector[Inf,Inf]
deflation = DeflationOperator(1.0, dot, 1.0, [ deflated_root ] )

u, = newton( F, J, x, p, NewtonPar() )                      # success

work. Basically, the methods of KrylovKit have to be defined (see here). How did you solve this?

ERROR: setindex!(::SArray{Tuple{2},Float64,1,2}, value, ::Int) is not defined.
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] setindex!(::SArray{Tuple{2},Float64,1,2}, ::Float64, ::Int64) at /Users/rveltz/.julia/packages/StaticArrays/l7lu2/src/indexing.jl:3
 [3] axpy!(::Int64, ::MArray{Tuple{2},Float64,1,2}, ::SArray{Tuple{2},Float64,1,2}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/generic.jl:1390
 [4] minus!(::SArray{Tuple{2},Float64,1,2}, ::MArray{Tuple{2},Float64,1,2}) at /Users/rveltz/work/prog_gd/julia/dev/dev1/BifurcationKit/src/BorderedArrays.jl:99
 [5] newton(::typeof(F), ::typeof(J), ::SArray{Tuple{2},Float64,1,2}, ::NamedTuple{(:k,),Tuple{Array{Float64,1}}}, ::NewtonPar{Float64,DefaultLS,DefaultEig{typeof(real)}}; normN::typeof(LinearAlgebra.norm), callback::Function, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/rveltz/work/prog_gd/julia/dev/dev1/BifurcationKit/src/Newton.jl:116
 [6] newton(::Function, ::Function, ::SArray{Tuple{2},Float64,1,2}, ::NamedTuple{(:k,),Tuple{Array{Float64,1}}}, ::NewtonPar{Float64,DefaultLS,DefaultEig{typeof(real)}}) at /Users/rveltz/work/prog_gd/julia/dev/dev1/BifurcationKit/src/Newton.jl:92
 [7] top-level scope at REPL[47]:1

[gszep]

One thing I noticed (maybe unimportant) is that F returns a MArray when given a SArray.

this is because similar(::SArray) returns an MArray which is mutable. Otherwise you can't assign values F[1] = .... setindex!(::SArray...) is not defined on purpose to prevent mutation. Since netwon has in-place methods similar(::SArray) should be used to convert to mutable type before other methods are called.

I cannot even make

########################################################### normal arrays
x = @SVector[-1.0,42.0]
deflated_root = @SVector[Inf,Inf]
deflation = DeflationOperator(1.0, dot, 1.0, [ deflated_root ] )

u, = newton( F, J, x, p, NewtonPar() )                      # success

work. Basically, the methods of KrylovKit have to be defined (see here). How did you solve this?

I am one subversion of BifurcationKit behind you I think pkg status reads

  [0f109fa4] BifurcationKit v0.1.0
  [f6369f11] ForwardDiff v0.10.12
  [d96e819e] Parameters v0.12.1
  [90137ffa] StaticArrays v0.12.4

[gszep]

Note that you may be interested in using DefaultLS(false) to remove the need to define factorize.

this is very useful thank you! Although it is the BorderedLinearSolver that eventually calls factorize which is sadly not accessible through ContinuationPar. Perhaps it should be? The way I set it at the moment is with the last positional argument linearAlgo of

Cont/PALCIterable(Fhandle, Jhandle, x0, par, lens::Lens, contParams::ContinuationPar, linearAlgo::AbstractBorderedLinearSolver)

[rveltz]

It should be, you pass MatrixBLS(DefaultLS(false)).

[gszep]

how? currently to pass DefaultLS I have to do this. I can't see a kwarg for to place a MatrixBLS struct

ContinuationPar(newtonOptions=NewtonPar(linsolver=DefaultLS(false)))

[rveltz]

no, you pass this to continuation as linearAlgo. There are 2 continuation, chose this one

If you pass linearAlgo = MatrixBLS(DefaultLS(false)) to the usual (simple) continuation, it might be overwritten here.

[gszep]

I'm using the iterator interface so I need it initialise ContIterable rather than continuation what I'm suggesting is that rather than having to pass ContinuationPar and linearAlgo in separate positional arguments, linearAlgo should be a field in the ContinuationPar struct and let continuation extract it from there. This way the user only needs to modify ContinuationPar to configure the method

[rveltz]

Ill think about this.