Comments (14)
rveltz
commented on May 29, 2024
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
from bifurcationkit.jl.
rveltz
commented on May 29, 2024
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)returnserror PosDefException: matrix is not positive definite; Cholesky factorization failed.
- Do you have your own
dotdefined? - what is the linear solver in
NewtonPar? - what is
linearAlgo?
from bifurcationkit.jl.
mbuze
commented on May 29, 2024
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
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 toNewtonPar(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 passlinearAlgo = 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!
from bifurcationkit.jl.
rveltz
commented on May 29, 2024
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.
from bifurcationkit.jl.
mbuze
commented on May 29, 2024
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.
from bifurcationkit.jl.
rveltz
commented on May 29, 2024
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.
from bifurcationkit.jl.
mbuze
commented on May 29, 2024
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
from bifurcationkit.jl.
rveltz
commented on May 29, 2024
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)
endfrom bifurcationkit.jl.
mbuze
commented on May 29, 2024
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.
from bifurcationkit.jl.
rveltz
commented on May 29, 2024
no problem!
from bifurcationkit.jl.
rveltz
commented on May 29, 2024
Hi,
Did you progress on this? Do you need a hand?
from bifurcationkit.jl.
rveltz
commented on May 29, 2024
hi,
Did you manage to solve it? Do you need a hand?
from bifurcationkit.jl.
rveltz
commented on May 29, 2024
Ping!
from bifurcationkit.jl.
rveltz
commented on May 29, 2024
it should be fixed with this commit
from bifurcationkit.jl.
Related Issues (20)
- Continuation of LPC from Bautin/GH point HOT 2
- Out of memory for 1D periodic non-linear system HOT 8
- Plotting with Makie HOT 3
- Getting print statements despite `verbosity=0` HOT 6
- BifurcationKit not working with ArrayPartitions ? HOT 1
- getperiod vs. getPeriod HOT 3
- _isinplace causes incorrect dispatch when both out of place and in place versions of a function is defined HOT 9
- Add `BifurcationProblem` dispatch for ModelingToolkit problems (possibly in an extension?) HOT 8
- BifurcationKit incompatible with latest ModelingToolkit version HOT 8
- record from solution title looks werid in docs HOT 2
- Incorrect calls to `EventSpecialPoint` HOT 2
- Problems with event detection HOT 5
- using `ContinuationPar` directly as input to `bifurcationdiagram` (rather than `(args...) -> opts_br`) HOT 2
- Custom BorderedLinearSolver interface requirements for PALC are missing HOT 4
- Generated function error HOT 2
- record_from_solution shoud have acces to the entire set of parameters HOT 2
- BK question HOT 8
- Problem with plotting a periodic orbit HOT 2
- (potentially long-term) Complete integration with ModelingTookit's nonlinear problem HOT 9
- Issue with passing the `record_from_solution` to the mimally augmented codim 2 continuation HOT 1
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from bifurcationkit.jl.