Changing `dotPALC` is not reflected in `BorderedPred()`?

Hi Romain,

Let me preface this by saying a huge thank you for your invaluable work on this package. At some point this year I am hoping to release my own amateurish package about numerical continuation methods in atomistic simulation of materials and it will be based around BifurcationKit.jl. I have also cited the paper about this package in a recent preprint.

Onto to the issue - title says it all really. It looks to me like any adjustement in dotPALC is not reflected in the computation of the predictor with BorderedPred(), which in my case leads to the breakdown of the continuation routine.

If in my continuation routine I use theta = 0.5 and

tangentAlgo = BorderedPred(), dotPALC = (x, y) -> dot(x, y)/length(x)

then everything runs well, but to achieve the desired step size control, I need to pass something dsmin = 0.00001, dsmax = 0.001, ds= 0.0004 (which would have to be decreased even further as if I was to increase the size of my computational domain - the number of atoms I simulate).

As noted in the documentation, I would like to instead pass

tangentAlgo = BorderedPred(), dotPALC = (x, y) -> dot(x, y)

but when I do that and run the routine, I get an error PosDefException: matrix is not positive definite; Cholesky factorization failed. If on the other hand I pass

tangentAlgo = SecantPred(), dotPALC = (x, y) -> dot(x, y)

then everything is well again and this time dsmin, dsmax, ds have to be adjusted to more reasonable values, which crucially are (more or less) domain-size independent.

I would really like to use BorderedPred() though and having browsed the source code, it looks like this might have something to do with the function `DotTheta()' in Predictor.jl. I think this should be easily fixable or perhaps I am simply misunderstanding something, but I thought I would draw your attention to it just in case. Thanks!

Hello,

Thank you for the preface, it is always nice to hear 😁

I need to pass something dsmin = 0.00001, dsmax = 0.001, ds= 0.0004 (which would have to be decreased even further as

You can mitigate this using theta. For example, theta = 0.9 would favour the parameter in the constraint equation.

but when I do that and run the routine, I get an error PosDefException: matrix is not positive definite; Cholesky factorization failed. If on the other hand I pass

I believe the normalisation is correct, see: https://github.com/rveltz/BifurcationKit.jl/blob/master/src/Predictor.jl#L135-L141

PosDefException: matrix is not positive definite; Cholesky factorization failed

I think I understand the error. You must be using a linear solver in NewtonPar.linsolver that works for positive definite matrices. However, for the bordered predictor, you must invert the matrix https://rveltz.github.io/BifurcationKit.jl/dev/Predictors/#Bordered-predictor which is not necessarily positive definite even if Fx is. You can get around this using the option linearAlgo = BorderingBLS(). See https://rveltz.github.io/BifurcationKit.jl/dev/borderedlinearsolver/ for more info on. the bordered linear solvers

Wait are you saying that?

tangentAlgo = BorderedPred(), dotPALC = (x, y) -> dot(x, y)/length(x) everything is OK. Do you pass this yourself or is it the default?
tangentAlgo = BorderedPred(), dotPALC = (x, y) -> dot(x, y) returns error PosDefException: matrix is not positive definite; Cholesky factorization failed.

Do you have your own dot defined?
what is the linear solver in NewtonPar?
what is linearAlgo?

With regards to just adjusting theta: I do get minor improvement by setting theta = 0.99 (or even closer to 1), but dividing by the length of x still forces me to choose tiny values fordsmin, dsmax, ds. Ideally I would like to shed a few zeros from them. And since the standard continuation constraint is

$N(x,p) = \frac{\theta}{L_x}\langle x-x_0, dx_0 \rangle + (1-\theta)(p - p_0)dp_0 - ds$

where L_x is the length of x, then getting rid of L_x is the only way to do it. Of course a workaround is to just divide by L_x when defining dsmin,dsmax,ds, but not settling for that led to finding this potential problem.

Here is the original setup I use which works:

optcont1 = ContinuationPar(dsmin = 0.000001, dsmax = 0.0003, ds= 0.00004, pMax = 4.1, maxSteps=100,
                   theta=0.99, newtonOptions = NewtonPar(tol = 1e-8),saveSolEveryStep = 0,doArcLengthScaling = false);

iter1 = BK.ContIterable(FF_m1, JJ_m1, x0, par, (@lens _.K0), optcont1; plot = false, verbosity = 0,
           tangentAlgo = BorderedPred(),dotPALC = (x, y) -> dot(x, y)/length(x))

If instead I pass tangentAlgo = BorderedPred(), dotPALC = (x, y) -> dot(x, y) and theta=0.5 (or something else more reasonable than 0.99) then I get error PosDefException: matrix is not positive definite; Cholesky factorization failed.. I have just checked though that if I decrease dsmin,dsmax,ds further, then sometimes this works for a few steps, but eventually I get the same error.

If however I pass tangentAlgo = SecantPred(),dotPALC = (x, y) -> dot(x, y) and theta=0.5 then everything works as intended - dsmin = 0.000001, dsmax = 0.0003, ds= 0.00004 now correspond to very tiny steps along the bifurcation diagram. The only problem is that shedding even a single zero from dsmin,dsmax,ds leads to the breakdown of the routine after several steps with the error again error PosDefException: matrix is not positive definite; Cholesky factorization failed.. This is to be expected in my kind of problems, hence the interest in using BorderedPred().

To answer your questions:

No, I use the standard dot product provided by Julia, but actually your questions prompted me to realise that I should really be using a slightly different notion of the dot product, thanks! (more on that below)
I just pass NewtonPar(tol = 1e-8), so it uses the default one. However, I have just changed that to NewtonPar(tol = 1e-8,linsolver = GMRESIterativeSolvers()) and then to `NewtonPar(tol = 1e-8,linsolver = GMRESKrylovKit())' and there is absolutely no change in the behaviour.
Again, just the default one, which judging by the documation is BorderingBLS()? I think so, because if I explicitly pass linearAlgo = BorderingBLS() then nothing changes.

Final note is that it could well be the issue is with the marginally incorrect notion of the dot product that I use (I should use a discrete equivalent of the (Sobolev space) H^1 inner product, as opposed to the discrete equivalent of L^2) and perhaps dividing by L_x is what is stabilising the whole thing. I am not sure though, because the setup with tangentAlgo = BorderedPred(),dotPALC = (x, y) -> dot(x, y)/length(x) is extremely stable - it can make huge jumps forward along the bifurcation diagram and pass through fold points in a very nonchalant fashion and everything still converges. Nonetheless, I will check what happens when I pass dotPALC = (x,y) -> doth1(x,y), where doth1 will be the manually defined function implementing the correct notion of the inner product.

Happy to hear your thoughts, perhaps there is some easy way out of this. If you are reasonably assurred that there are no bugs in the code, please do not waste too much time on this :) Thanks!

If you are reasonably assurred that there are no bugs in the code, please do not waste too much time on this :) Thanks!

That's the thing, issues are a useful way to catch them. Can you please paste the full error stack here?

If your code is not too sensitive, you can also try to send it to me directly.

Here is the full error stack:

PosDefException: matrix is not positive definite; Cholesky factorization failed.

Stacktrace:
 [1] #ldlt!#10(::Float64, ::Bool, ::typeof(ldlt!), ::SuiteSparse.CHOLMOD.Factor{Float64}, ::SuiteSparse.CHOLMOD.Sparse{Float64}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/SuiteSparse/src/cholmod.jl:1448
 [2] #ldlt!#11 at ./none:0 [inlined]
 [3] ldlt! at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/SuiteSparse/src/cholmod.jl:1473 [inlined]
 [4] factorize(::Hermitian{Float64,SparseMatrixCSC{Float64,Int64}}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/SparseArrays/src/linalg.jl:1497
 [5] factorize(::SparseMatrixCSC{Float64,Int64}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/SparseArrays/src/linalg.jl:1475
 [6] #_#6(::Int64, ::Int64, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::DefaultLS, ::SparseMatrixCSC{Float64,Int64}, ::Array{Float64,1}, ::Array{Float64,1}) at /Users/maciejbuze/.julia/packages/BifurcationKit/Yvq7p/src/LinearSolver.jl:83
 [7] DefaultLS at /Users/maciejbuze/.julia/packages/BifurcationKit/Yvq7p/src/LinearSolver.jl:82 [inlined]
 [8] #_#44(::Nothing, ::BorderingBLS{DefaultLS,Float64}, ::SparseMatrixCSC{Float64,Int64}, ::Array{Float64,1}, ::Array{Float64,1}, ::Float64, ::Array{Float64,1}, ::Float64, ::Float64, ::Float64) at /Users/maciejbuze/.julia/packages/BifurcationKit/Yvq7p/src/LinearBorderSolver.jl:41
 [9] Any at ./none:0 [inlined]
 [10] #_#39 at /Users/maciejbuze/.julia/packages/BifurcationKit/Yvq7p/src/LinearBorderSolver.jl:8 [inlined]
 [11] (::BorderingBLS{DefaultLS,Float64})(::SparseMatrixCSC{Float64,Int64}, ::Array{Float64,1}, ::BorderedArray{Array{Float64,1},Float64}, ::Array{Float64,1}, ::Float64, ::Float64) at /Users/maciejbuze/.julia/packages/BifurcationKit/Yvq7p/src/LinearBorderSolver.jl:8
 [12] #newtonPALC#182(::BorderingBLS{DefaultLS,Float64}, ::typeof(norm), ::Function, ::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol},NamedTuple{(:iterationC, :z0),Tuple{Int64,BorderedArray{Array{Float64,1},Float64}}}}, ::typeof(newtonPALC), ::var"#FF_m1#221"{NCFlex.NCF,Int64}, ::var"#JJ_m1#222"{NCFlex.NCF,Int64}, ::NamedTuple{(:K0,),Tuple{Float64}}, ::Setfield.PropertyLens{:K0}, ::BorderedArray{Array{Float64,1},Float64}, ::BorderedArray{Array{Float64,1},Float64}, ::BorderedArray{Array{Float64,1},Float64}, ::Float64, ::Float64, ::ContinuationPar{Float64,GMRESIterativeSolvers{Float64,IterativeSolvers.Identity,IterativeSolvers.Identity},DefaultEig{typeof(real)}}, ::BifurcationKit.DotTheta{var"#227#228"}) at /Users/maciejbuze/.julia/packages/BifurcationKit/Yvq7p/src/Predictor.jl:495
 [13] #iterate#139(::Int64, ::typeof(iterate), ::ContIterable{var"#FF_m1#221"{NCFlex.NCF,Int64},var"#JJ_m1#222"{NCFlex.NCF,Int64},Array{Float64,1},NamedTuple{(:K0,),Tuple{Float64}},Setfield.PropertyLens{:K0},Float64,GMRESIterativeSolvers{Float64,IterativeSolvers.Identity,IterativeSolvers.Identity},DefaultEig{typeof(real)},BorderedPred,BorderingBLS{DefaultLS,Float64},BifurcationKit.var"#118#128"{BifurcationKit.var"#118#119#129"},BifurcationKit.var"#120#130"{BifurcationKit.var"#120#121#131"},typeof(norm),BifurcationKit.DotTheta{var"#227#228"},typeof(BifurcationKit.finaliseDefault),typeof(BifurcationKit.cbDefault),Nothing,String}, ::ContState{BorderedArray{Array{Float64,1},Float64},Float64,Nothing,Nothing,Tuple{Nothing,Nothing}}) at ./none:0
 [14] iterate(::ContIterable{var"#FF_m1#221"{NCFlex.NCF,Int64},var"#JJ_m1#222"{NCFlex.NCF,Int64},Array{Float64,1},NamedTuple{(:K0,),Tuple{Float64}},Setfield.PropertyLens{:K0},Float64,GMRESIterativeSolvers{Float64,IterativeSolvers.Identity,IterativeSolvers.Identity},DefaultEig{typeof(real)},BorderedPred,BorderingBLS{DefaultLS,Float64},BifurcationKit.var"#118#128"{BifurcationKit.var"#118#119#129"},BifurcationKit.var"#120#130"{BifurcationKit.var"#120#121#131"},typeof(norm),BifurcationKit.DotTheta{var"#227#228"},typeof(BifurcationKit.finaliseDefault),typeof(BifurcationKit.cbDefault),Nothing,String}, ::ContState{BorderedArray{Array{Float64,1},Float64},Float64,Nothing,Nothing,Tuple{Nothing,Nothing}}) at /Users/maciejbuze/.julia/packages/BifurcationKit/Yvq7p/src/Continuation.jl:311
 [15] #simple_continuation2#218(::Float64, ::Float64, ::Float64, ::Float64, ::Int64, ::Float64, ::BorderedPred, ::Function, ::typeof(simple_continuation2), ::NCFlex.NCF) at ./In[168]:95
 [16] (::var"#kw##simple_continuation2")(::NamedTuple{(:theta, :maxSteps, :dsmin, :dsmax, :ds, :tangentAlgo, :dotPALC),Tuple{Float64,Int64,Float64,Float64,Float64,BorderedPred,var"#227#228"}}, ::typeof(simple_continuation2), ::NCFlex.NCF) at ./none:0
 [17] top-level scope at In[170]:1

This is obtained with the following commands passed:

optcont1 = ContinuationPar(theta=0.5, maxSteps=20, dsmin = 0.00005, dsmax = 0.0005, ds= 0.00005, pMax = 4.1, 
                   newtonOptions = NewtonPar(tol = 1e-8,linsolver = GMRESIterativeSolvers()),
                   saveSolEveryStep = 0,doArcLengthScaling = false);
iter1 = BK.ContIterable(FF_m1, JJ_m1, x0, par, (@lens _.K0), optcont1; plot = false, verbosity = 0,
           tangentAlgo = BorderedPred(),dotPALC = (x, y) -> dot(x, y), linearAlgo = MatrixBLS(GMRESIterativeSolvers()))

Not sure if MatrixBLS(GMRESIterativeSolvers()) is the correct way to go about it (checked other possibilities too).

I should also perhaps add that I am actually using a customised routine through the Iterator interface. In this particular case, thanks to the tiny choices for ds,dsmin,dsmax, the routine works fine for a few tiny steps and then I get the error. If I instead pass tangentAlgo = SecantPred() then it works just fine.

If needed, I can share the full routine too.

Just to be sure:

are you on master?

This is strange because the linearsolver is :: GMRESIterativeSolvers for newton, BorderingBLS for newtonPALC. It seems the predictor for the tangent went fine. It is stuck in the corrector and perhaps it is a bug in SparseArrays. You check this by passing the jacobian like this (x,p)->lu(JJ_m1(x,p)).

Can you please try linearAlgo = MatrixBLS() and pass newtonOptions = NewtonPar(tol = 1e-8) to use the CHOLDMOD sparse solver.

If needed, I can share the full routine too.

Yes, please do. It'll be faster to debug if your code is not too big.

The code is a mess at the moment, I will try to create a MWE tomorrow. In the meantime, I updated to Julia 1.6 and to the master branch of BifurcationKit. The problem still persists, but the error message is different:

ZeroPivotException: factorization encountered one or more zero pivots. Consider switching to a pivoted LU factorization.

Stacktrace:
  [1] ldlt!(F::SuiteSparse.CHOLMOD.Factor{Float64}, A::SuiteSparse.CHOLMOD.Sparse{Float64}; shift::Float64, check::Bool)
    @ SuiteSparse.CHOLMOD /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SuiteSparse/src/cholmod.jl:1472
  [2] #ldlt!#11
    @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SuiteSparse/src/cholmod.jl:1497 [inlined]
  [3] ldlt!
    @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SuiteSparse/src/cholmod.jl:1497 [inlined]
  [4] factorize(A::Hermitian{Float64, SparseMatrixCSC{Float64, Int64}})
    @ SparseArrays /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:1634
  [5] factorize(A::SparseMatrixCSC{Float64, Int64})
    @ SparseArrays /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:1612
  [6] (::DefaultLS)(J::SparseMatrixCSC{Float64, Int64}, rhs1::Vector{Float64}, rhs2::Vector{Float64}; a₀::Int64, a₁::Int64, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/oRkQZ/src/LinearSolver.jl:83
  [7] DefaultLS
    @ ~/.julia/packages/BifurcationKit/oRkQZ/src/LinearSolver.jl:82 [inlined]
  [8] (::BorderingBLS{DefaultLS, Float64})(J::SparseMatrixCSC{Float64, Int64}, dR::Vector{Float64}, dzu::Vector{Float64}, dzp::Float64, R::Vector{Float64}, n::Float64, xiu::Float64, xip::Float64; shift::Nothing)
    @ BifurcationKit ~/.julia/packages/BifurcationKit/oRkQZ/src/LinearBorderSolver.jl:41
  [9] #_#39
    @ ~/.julia/packages/BifurcationKit/oRkQZ/src/LinearBorderSolver.jl:8 [inlined]
 [10] (::BorderingBLS{DefaultLS, Float64})(J::SparseMatrixCSC{Float64, Int64}, dR::Vector{Float64}, dz::BorderedArray{Vector{Float64}, Float64}, R::Vector{Float64}, n::Float64, theta::Float64)
    @ BifurcationKit ~/.julia/packages/BifurcationKit/oRkQZ/src/LinearBorderSolver.jl:8
 [11] newtonPALC(F::var"#FF_m1#27"{NCFlex.NCF, Int64}, Jh::var"#JJ_m1#28"{NCFlex.NCF, Int64}, par::NamedTuple{(:K0,), Tuple{Float64}}, paramlens::Setfield.PropertyLens{:K0}, z0::BorderedArray{Vector{Float64}, Float64}, τ0::BorderedArray{Vector{Float64}, Float64}, z_pred::BorderedArray{Vector{Float64}, Float64}, ds::Float64, θ::Float64, contparams::ContinuationPar{Float64, DefaultLS, DefaultEig{typeof(real)}}, dottheta::BifurcationKit.DotTheta{var"#35#36"}; linearbdalgo::BorderingBLS{DefaultLS, Float64}, normN::typeof(norm), callback::Function, kwargs::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:iterationC, :z0), Tuple{Int64, BorderedArray{Vector{Float64}, Float64}}}})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/oRkQZ/src/Predictor.jl:495
 [12] #newtonPALC#189
    @ ~/.julia/packages/BifurcationKit/oRkQZ/src/Predictor.jl:547 [inlined]
 [13] #corrector#168
    @ ~/.julia/packages/BifurcationKit/oRkQZ/src/Predictor.jl:68 [inlined]
 [14] iterate(it::ContIterable{var"#FF_m1#27"{NCFlex.NCF, Int64}, var"#JJ_m1#28"{NCFlex.NCF, Int64}, Vector{Float64}, NamedTuple{(:K0,), Tuple{Float64}}, Setfield.PropertyLens{:K0}, Float64, DefaultLS, DefaultEig{typeof(real)}, BorderedPred, BorderingBLS{DefaultLS, Float64}, BifurcationKit.var"#120#130"{BifurcationKit.var"#120#121#131"}, BifurcationKit.var"#122#132"{BifurcationKit.var"#122#123#133"}, typeof(norm), BifurcationKit.DotTheta{var"#35#36"}, typeof(BifurcationKit.finaliseDefault), typeof(BifurcationKit.cbDefault), Nothing, String}, state::ContState{BorderedArray{Vector{Float64}, Float64}, Float64, Nothing, Nothing, Tuple{Nothing, Nothing}}; _verbosity::Int64)
    @ BifurcationKit ~/.julia/packages/BifurcationKit/oRkQZ/src/Continuation.jl:328
 [15] iterate(it::ContIterable{var"#FF_m1#27"{NCFlex.NCF, Int64}, var"#JJ_m1#28"{NCFlex.NCF, Int64}, Vector{Float64}, NamedTuple{(:K0,), Tuple{Float64}}, Setfield.PropertyLens{:K0}, Float64, DefaultLS, DefaultEig{typeof(real)}, BorderedPred, BorderingBLS{DefaultLS, Float64}, BifurcationKit.var"#120#130"{BifurcationKit.var"#120#121#131"}, BifurcationKit.var"#122#132"{BifurcationKit.var"#122#123#133"}, typeof(norm), BifurcationKit.DotTheta{var"#35#36"}, typeof(BifurcationKit.finaliseDefault), typeof(BifurcationKit.cbDefault), Nothing, String}, state::ContState{BorderedArray{Vector{Float64}, Float64}, Float64, Nothing, Nothing, Tuple{Nothing, Nothing}})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/oRkQZ/src/Continuation.jl:311
 [16] simple_continuation2(ncf::NCFlex.NCF; dsmin::Float64, dsmax::Float64, ds::Float64, pMax::Float64, maxSteps::Int64, theta::Float64, tangentAlgo::BorderedPred, dotPALC::Function, linsolver::DefaultLS, linearAlgo::MatrixBLS{Nothing})
    @ Main ./In[18]:96
 [17] top-level scope
    @ In[23]:1
 [18] eval
    @ ./boot.jl:360 [inlined]
 [19] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1094

It looks an issue of similar kind have already been reported, e.g. here: https://discourse.julialang.org/t/zero-pivot-error-in-ldlt-factorization-for-large-matrices/1966

Passing the jacobian as lu(JJ_m1(x,p)) results in a different error which seems to be a consquence of the fact that one cannot combine factorize and lu. I tried passing DefaultLS(useFactorization=false) to stop it from attempting to factorize it, but this does not seem to work. Here is the error message:

MethodError: no method matching factorize(::SuiteSparse.UMFPACK.UmfpackLU{Float64, Int64})
Closest candidates are:
  factorize(::StridedMatrix{T}) where T at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/dense.jl:1241
  factorize(::Union{Hermitian{ComplexF64, var"#s832"}, Hermitian{Float64, var"#s832"}, Symmetric{Float64, var"#s832"}} where var"#s832"<:SparseArrays.AbstractSparseMatrixCSC) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:1629
  factorize(::Union{Hermitian{T, S}, Symmetric{T, S}} where {T, S}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/symmetric.jl:637
  ...

Stacktrace:
  [1] #_#6
    @ ~/.julia/packages/BifurcationKit/oRkQZ/src/LinearSolver.jl:83 [inlined]
  [2] DefaultLS
    @ ~/.julia/packages/BifurcationKit/oRkQZ/src/LinearSolver.jl:82 [inlined]
  [3] (::BorderingBLS{DefaultLS, Float64})(J::SuiteSparse.UMFPACK.UmfpackLU{Float64, Int64}, dR::Vector{Float64}, dzu::Vector{Float64}, dzp::Float64, R::Vector{Float64}, n::Float64, xiu::Float64, xip::Float64; shift::Nothing)
    @ BifurcationKit ~/.julia/packages/BifurcationKit/oRkQZ/src/LinearBorderSolver.jl:41
  [4] #_#39
    @ ~/.julia/packages/BifurcationKit/oRkQZ/src/LinearBorderSolver.jl:8 [inlined]
  [5] (::BorderingBLS{DefaultLS, Float64})(J::SuiteSparse.UMFPACK.UmfpackLU{Float64, Int64}, dR::Vector{Float64}, dz::BorderedArray{Vector{Float64}, Float64}, R::Vector{Float64}, n::Float64, theta::Float64)
    @ BifurcationKit ~/.julia/packages/BifurcationKit/oRkQZ/src/LinearBorderSolver.jl:8
  [6] newtonPALC(F::var"#FF_m1#40"{NCFlex.NCF, Int64}, Jh::var"#JJ_m1#41"{NCFlex.NCF, Int64}, par::NamedTuple{(:K0,), Tuple{Float64}}, paramlens::Setfield.PropertyLens{:K0}, z0::BorderedArray{Vector{Float64}, Float64}, τ0::BorderedArray{Vector{Float64}, Float64}, z_pred::BorderedArray{Vector{Float64}, Float64}, ds::Float64, θ::Float64, contparams::ContinuationPar{Float64, DefaultLS, DefaultEig{typeof(real)}}, dottheta::BifurcationKit.DotTheta{var"#44#45"}; linearbdalgo::BorderingBLS{DefaultLS, Float64}, normN::typeof(norm), callback::Function, kwargs::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:iterationC, :z0), Tuple{Int64, BorderedArray{Vector{Float64}, Float64}}}})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/oRkQZ/src/Predictor.jl:495
  [7] #newtonPALC#189
    @ ~/.julia/packages/BifurcationKit/oRkQZ/src/Predictor.jl:547 [inlined]
  [8] #corrector#168
    @ ~/.julia/packages/BifurcationKit/oRkQZ/src/Predictor.jl:68 [inlined]
  [9] iterate(it::ContIterable{var"#FF_m1#40"{NCFlex.NCF, Int64}, var"#JJ_m1#41"{NCFlex.NCF, Int64}, Vector{Float64}, NamedTuple{(:K0,), Tuple{Float64}}, Setfield.PropertyLens{:K0}, Float64, DefaultLS, DefaultEig{typeof(real)}, BorderedPred, BorderingBLS{DefaultLS, Float64}, BifurcationKit.var"#120#130"{BifurcationKit.var"#120#121#131"}, BifurcationKit.var"#122#132"{BifurcationKit.var"#122#123#133"}, typeof(norm), BifurcationKit.DotTheta{var"#44#45"}, typeof(BifurcationKit.finaliseDefault), typeof(BifurcationKit.cbDefault), Nothing, String}, state::ContState{BorderedArray{Vector{Float64}, Float64}, Float64, Nothing, Nothing, Tuple{Nothing, Nothing}}; _verbosity::Int64)
    @ BifurcationKit ~/.julia/packages/BifurcationKit/oRkQZ/src/Continuation.jl:328
 [10] iterate(it::ContIterable{var"#FF_m1#40"{NCFlex.NCF, Int64}, var"#JJ_m1#41"{NCFlex.NCF, Int64}, Vector{Float64}, NamedTuple{(:K0,), Tuple{Float64}}, Setfield.PropertyLens{:K0}, Float64, DefaultLS, DefaultEig{typeof(real)}, BorderedPred, BorderingBLS{DefaultLS, Float64}, BifurcationKit.var"#120#130"{BifurcationKit.var"#120#121#131"}, BifurcationKit.var"#122#132"{BifurcationKit.var"#122#123#133"}, typeof(norm), BifurcationKit.DotTheta{var"#44#45"}, typeof(BifurcationKit.finaliseDefault), typeof(BifurcationKit.cbDefault), Nothing, String}, state::ContState{BorderedArray{Vector{Float64}, Float64}, Float64, Nothing, Nothing, Tuple{Nothing, Nothing}})
    @ BifurcationKit ~/.julia/packages/BifurcationKit/oRkQZ/src/Continuation.jl:311
 [11] simple_continuation2(ncf::NCFlex.NCF; dsmin::Float64, dsmax::Float64, ds::Float64, pMax::Float64, maxSteps::Int64, theta::Float64, tangentAlgo::BorderedPred, dotPALC::Function, linsolver::DefaultLS, linearAlgo::MatrixBLS{Nothing})
    @ Main ./In[24]:43
 [12] top-level scope
    @ In[25]:1
 [13] eval
    @ ./boot.jl:360 [inlined]
 [14] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1094

These errors are difficult to understand without the call signature, so if you have a MWE it will be simpler (I guess). There is another possibility to the error PosDefException: matrix is not positive definite which is that the current state in newton iterations is not valid. This happens for example when residual are huge, like 1e9. Did you monitor the residuals?

As for the LU factorisation, I dont get the issue. One workaround is to keep your original jacobian and pass to NewtonPar the linear solver, like NewtonPar(linsolver = LinSolveLU())

struct LinSolveLU <: BifurcationKit.AbstractLinearSolver; end

function (l::LinSolveLU)(J, rhs; a₀ = 0, a₁ = 1, kwargs...)
    return lu(_axpy(J, a₀, a₁)) \ rhs, true, 1
end

function (l::LinSolveLU)(J, rhs1, rhs2; a₀ = 0, a₁ = 1, kwargs...)
    Jfact = lu(_axpy(J, a₀, a₁))
    return Jfact \ rhs1, Jfact \ rhs2, true, (1, 1)
end

Hi Romain, I'm really sorry about the delay. I'm very busy this week. I'll return to this next week if that's okay.

Hi,

Did you progress on this? Do you need a hand?

hi,

Did you manage to solve it? Do you need a hand?

it should be fixed with this commit

bifurcationkit / BifurcationKit.jl

Changing `dotPALC` is not reflected in `BorderedPred()`? #38