add update_minaug_every_step option in codim 2 PO

src/Problems.jl

@@ -282,7 +282,7 @@ for (op, at) in (
                     Finp = _F
                 else
                     F = _F
-                    Finp = (o, x, p) -> o .= _F(x, p)
+                    Finp = (o, x, p) -> copyto!(o, _F(x, p))
                 end
 
                 J = isnothing(J) ? (x, p) -> ForwardDiff.jacobian(z -> F(z, p), x) : J

src/codim2/MinAugFold.jl

@@ -184,6 +184,7 @@ end
 @inline getdelta(foldpb::FoldMAProblem) = getdelta(foldpb.prob)
 @inline is_symmetric(foldpb::FoldMAProblem) = is_symmetric(foldpb.prob)
 residual(foldpb::FoldMAProblem, x, p) = foldpb.prob(x, p)
+residual!(foldpb::FoldMAProblem, out, x, p) = (copyto!(out, foldpb.prob(x, p)); out)
 jad(foldpb::FoldMAProblem, args...) = jad(foldpb.prob, args...)
 save_solution(::FoldMAProblem, x, p) = x
 

src/codim2/MinAugHopf.jl

@@ -199,6 +199,7 @@ end
 @inline is_symmetric(hopfpb::HopfMAProblem) = is_symmetric(hopfpb.prob)
 @inline getdelta(hopfpb::HopfMAProblem) = getdelta(hopfpb.prob)
 residual(hopfpb::HopfMAProblem, x, p) = hopfpb.prob(x, p)
+residual!(hopfpb::HopfMAProblem, out, x, p) = (copyto!(out, hopfpb.prob(x, p)); out)
 save_solution(::HopfMAProblem, x ,p) = x
 
 # jacobian(hopfpb::HopfMAProblem, x, p) = hopfpb.jacobian(x, p)

src/continuation/Contbase.jl

@@ -64,7 +64,7 @@ function _step_size_control!(state, contparams::ContinuationPar, verbosity)
     ds = state.ds
     if converged(state) == false
         if  abs(ds) <= contparams.dsmin
-            @error "Failure to converge with given tolerance = $(contparams.newton_options.tol).\n You can decrease the tolerance or pass a different norm using the argument `normC`.\n We reached the smallest value [dsmin] valid for ds, namely $(contparams.dsmin).\n Stopping continuation at continuation step $(state.step)."
+            @error "Failure to converge with given tolerance = $(contparams.newton_options.tol).\nStep = $(state.step)\n You can decrease the tolerance or pass a different norm using the argument `normC`.\n We reached the smallest value [dsmin] valid for ds, namely $(contparams.dsmin).\n Stopping continuation at continuation step $(state.step)."
             # we stop the continuation
             state.stopcontinuation = true
             return

src/periodicorbit/PeriodicOrbitCollocation.jl

@@ -304,7 +304,7 @@ function Base.show(io::IO, pb::PeriodicOrbitOCollProblem)
     if pb.meshadapt
         println(io, "├───── K              : ", pb.K)
     end
-    println(io, "└─ # unknowns without phase condition) : ", pb.N * (1 + m * Ntst))
+    println(io, "└─ # unknowns (without phase condition) : ", pb.N * (1 + m * Ntst))
 end
 
 function get_matrix_phase_condition(coll::PeriodicOrbitOCollProblem)
@@ -629,7 +629,7 @@ Compute the jacobian of the problem defining the periodic orbits by orthogonal c
 
             for l2 in 1:m+1
                 J[_rgX, rgNy .+ (l2-1)*n ] .= @. (-α * L[l2, l] * ρF) * J0 +
-                                                        (ρD * ∂L[l2, l] - α * L[l2, l] * ρI) * In
+                                                 (ρD * ∂L[l2, l] - α * L[l2, l] * ρI) * In
             end
             # add derivative w.r.t. the period
             J[rgNx .+ (l-1)*n, end] .= residual(VF, pj[:,l], pars) .* (-dt)
@@ -645,15 +645,13 @@ Compute the jacobian of the problem defining the periodic orbits by orthogonal c
         for k₁ = 1:m+1
             for k₂ = 1:m+1
                 # J[end, rg] .+= Ω[k₁, k₂] .* ϕc[:, (j-1)*m + k₂]
-                axpby!(Ω[k₁, k₂], ϕc[:, (j-1)*m + k₂], 1, J[end, rg])
+                axpby!(Ω[k₁, k₂] / period, ϕc[:, (j-1)*m + k₂], 1, J[end, rg])
             end
             if k₁ < m + 1
                 rg = rg .+ n
             end
         end
     end
-    J[end, 1:end-1] ./= period
-
     vj = get_tmp(coll.cache.vj, u)
     phase = _phase_condition(coll, uc, (L, ∂L), (pj, uj, ϕj, vj), period)
     J[end, end] = -phase / period

src/periodicorbit/PeriodicOrbitTrapeze.jl

@@ -344,12 +344,12 @@ function Jc(pb::PeriodicOrbitTrapProblem, outc::AbstractMatrix, u0::AbstractVect
 
     h = T * get_time_step(pb, 1)
     @views potrap_scheme!(pb, outc[:, 1], u0c[:, 1], u0c[:, M-1],
-                                         duc[:, 1], duc[:, M-1], par, h/2, tmp, true; applyf = false)
+                                          duc[:, 1], duc[:, M-1], par, h/2, tmp, true; applyf = false)
 
     for ii in 2:M-1
         h = T * get_time_step(pb, ii)
         @views potrap_scheme!(pb, outc[:, ii], u0c[:, ii], u0c[:, ii-1],
-                                              duc[:, ii], duc[:, ii-1], par, h/2, tmp, true; applyf = false)
+                                               duc[:, ii], duc[:, ii-1], par, h/2, tmp, true; applyf = false)
     end
 
     # we also return a Vector version of outc

src/periodicorbit/PoincareShooting.jl

@@ -290,6 +290,11 @@ function (psh::PoincareShootingProblem)(x_bar::AbstractVector, par; verbose = fa
     return out_bar
 end
 
+function residual!(pb::PoincareShootingProblem, out, x, p)
+    copyto!(out, pb(x, p))
+    out
+end
+
 """
 This function computes the derivative of the Poincare return map Π(x) = ϕ(t(x),x) where t(x) is the return time of x to the section.
 """

src/periodicorbit/StandardShooting.jl

@@ -289,6 +289,12 @@ end
 
 # out of place version
 (sh::ShootingProblem)(::Val{:JacobianMatrix}, x::AbstractVector, par) = sh(Val(:JacobianMatrixInplace), zeros(eltype(x), length(x), length(x)), x, par)
+
+function residual!(pb::ShootingProblem, out, x, p)
+    copyto!(out, pb(x, p))
+    out
+end
+residual(pb::ShootingProblem, x, p) = pb(x, p)
 ####################################################################################################
 
 function _get_extremum(prob::ShootingProblem, x::AbstractVector, p; ratio = 1, op = (max, maximum))

src/periodicorbit/codim2/MinAugNS.jl

@@ -260,6 +260,7 @@ function (pdls::NSLinearSolverMinAug)(Jns, rhs::BorderedArray{vectype, 𝒯}; de
 end
 ###################################################################################################
 residual(nspb::NSMAProblem, x, p) = nspb.prob(x, p)
+residual!(nspb::NSMAProblem, out, x, p) = (copyto!(out, nspb.prob(x, p)); out)
 @inline getdelta(nspb::NSMAProblem) = getdelta(nspb.prob)
 save_solution(::NSMAProblem, x ,p) = x
 

src/periodicorbit/codim2/MinAugPD.jl

@@ -218,6 +218,7 @@ end
 @inline is_symmetric(pdpb::PDMAProblem) = is_symmetric(pdpb.prob)
 @inline getdelta(pdpb::PDMAProblem) = getdelta(pdpb.prob)
 residual(pdpb::PDMAProblem, x, p) = pdpb.prob(x, p)
+residual!(pdpb::PDMAProblem, out, x, p) = (copyto!(out, pdpb.prob(x, p)); out)
 save_solution(::PDMAProblem, x, p) = x
 
 jacobian(pdpb::PDMAProblem{Tprob, Nothing, Tu0, Tp, Tl, Tplot, Trecord}, x, p) where {Tprob, Tu0, Tp, Tl <: Union{AllOpticTypes, Nothing}, Tplot, Trecord} = (x = x, params = p, prob = pdpb.prob)

src/periodicorbit/codim2/PeriodicOrbitCollocation.jl

@@ -46,19 +46,23 @@ function continuation(br::AbstractResult{Tkind, Tprob},
         return continuation_coll_fold(br,
                                     ind_bif,
                                     lens2,
-                                    _options_cont; compute_eigen_elements,
+                                    _options_cont; 
+                                    compute_eigen_elements,
+                                    update_minaug_every_step,
                                     kwargs... )
     elseif biftype == :pd
         return continuation_coll_pd(br,
                                     ind_bif,
                                     lens2,
                                     _options_cont; compute_eigen_elements,
+                                    update_minaug_every_step,
                                     kwargs... )
     elseif biftype == :ns
         return continuation_coll_ns(br,
                                     ind_bif,
                                     lens2,
                                     _options_cont; compute_eigen_elements,
+                                    update_minaug_every_step,
                                     kwargs... )
     else
         throw("We continue only Fold / PD / NS points of periodic orbits. Please report this error on the website.")

src/periodicorbit/codim2/StandardShooting.jl

@@ -155,7 +155,7 @@ function continuation(br::AbstractResult{Tkind, Tprob},
     _options_cont = detect_codim2_parameters(detect_codim2_bifurcation, options_cont; update_minaug_every_step, kwargs...)
 
     if biftype == :bp || biftype == :fold
-        return continuation_sh_fold(br, ind_bif, lens2, _options_cont; compute_eigen_elements, kwargs... )
+        return continuation_sh_fold(br, ind_bif, lens2, _options_cont; compute_eigen_elements, update_minaug_every_step, kwargs... )
     elseif biftype == :pd
         return continuation_sh_pd(br, ind_bif, lens2, _options_cont; compute_eigen_elements, update_minaug_every_step, kwargs... )
     elseif biftype == :ns

src/periodicorbit/cop.jl

@@ -189,22 +189,18 @@ Solve the linear system associated with the collocation problem for computing pe
     𝒯 = eltype(coll)
     In = I(N)
 
-    if _DEBUG
-        P = Matrix{𝒯}(LinearAlgebra.I(nⱼ))
-        Jtmp = zeros(𝒯, nbcoll + δn + 1, nbcoll)
-    end
-
     rhs = condensation_of_parameters!(cop_cache, coll, J, In, rhs0)
     Jcond = cop_cache.Jcoll
 
     if _DEBUG
+        P = Matrix{𝒯}(LinearAlgebra.I(nⱼ))
+        Jtmp = zeros(𝒯, nbcoll + δn + 1, nbcoll)
         Fₚ = lu(P); Jcond = Fₚ \ J; rhs = Fₚ \ rhs0
     end
 
     # last_row_𝐅𝐬⁻¹_analytical = zeros(𝒯, δn + 1, nⱼ) # last row of 𝐅𝐬⁻¹
     # last_row_𝐅𝐬 = zeros(𝒯, δn + 1, nⱼ) # last row of 𝐅𝐬
-    @unpack last_row_𝐅𝐬⁻¹_analytical,
-            last_row_𝐅𝐬 = cop_cache
+    @unpack last_row_𝐅𝐬⁻¹_analytical = cop_cache
 
     if dim == 0 
         d = dot(last_row_𝐅𝐬⁻¹_analytical, 
@@ -213,29 +209,31 @@ Solve the linear system associated with the collocation problem for computing pe
         rhs[end] = dot(last_row_𝐅𝐬⁻¹_analytical, 
                 rhs0[eachindex(last_row_𝐅𝐬⁻¹_analytical)]) +
                 rhs0[end]
+        Jcond[end-δn:end, end-δn:end] .= d
     else
-        d = last_row_𝐅𝐬⁻¹_analytical * 
-            J[axes(last_row_𝐅𝐬⁻¹_analytical, 2), end-δn:end] .+ 
-            J[end-δn:end, end-δn:end]
+        # d = last_row_𝐅𝐬⁻¹_analytical * 
+        #     J[axes(last_row_𝐅𝐬⁻¹_analytical, 2), end-δn:end] .+ 
+        #     J[end-δn:end, end-δn:end]
+        # Jcond[end-δn:end, end-δn:end] .= d
+
+        Jcond[end-δn:end, end-δn:end] .= J[end-δn:end, end-δn:end]
+        mul!(Jcond[end-δn:end, end-δn:end], last_row_𝐅𝐬⁻¹_analytical, J[axes(last_row_𝐅𝐬⁻¹_analytical, 2), end-δn:end], true, true)
+
         rhs[end-δn:end] .= last_row_𝐅𝐬⁻¹_analytical *
             rhs0[axes(last_row_𝐅𝐬⁻¹_analytical, 2)] .+
             rhs0[end-δn:end]
     end
-    Jcond[end-δn:end, end-δn:end] .= d
-
-    # plot(heatmap(abs.(abs.(inv(P))) .> 1e-5; yflip = true, title = "invP"), 
-        # heatmap(abs.(Jcond - Jcop) .> 1e-5; yflip = true, title = "δJ")) |> display
 
     # we build the linear system for the external variables in Jext and rhs_ext
     rhs_ext = build_external_system!(Jext, Jcond, rhs, In, Ntst, nbcoll, Npo, δn, N, m)
 
-    if !_USELU
+    if _USELU
+        F = lu(Jext)
+        sol_ext = F \ rhs_ext
+    else
         # gaussian elimination plus backward substitution to invert Jext
         _gaussian_elimination_external_pivoted!(Jext, rhs_ext, N, Ntst, δn)
         sol_ext = _backward_substitution_pivoted(Jext, rhs_ext, N, Ntst, Val(dim))
-    else
-        F = lu(Jext)
-        sol_ext = F \ rhs_ext
     end
 
     return _solve_for_internal_variables(coll, Jcond, rhs, sol_ext, Val(dim))
@@ -256,7 +254,7 @@ end
     # δn = 0 for newton
     # δn = 1 for palc
     @assert δn >= 0
-    @assert δn == dim
+    @assert δn == dim "δn = $δn == dim = $dim"
 
     𝒯 = eltype(coll)
 
@@ -402,27 +400,21 @@ end
                                        δn::Int,
                                        N::Int,
                                        m::Int) where {𝒯}
-    Aᵢ = Matrix{𝒯}(undef, N, N)
-    Bᵢ = Matrix{𝒯}(undef, N, N)
-
     r1 = 1:N
     r2 = N*(m-1)+1:(m*N)
     rN = 1:N
 
     # building the external variables
     fill!(Jext, 0)
-    Jext[end-δn-N:end-δn-1,end-δn-N:end-δn-1] .= In
-    Jext[end-δn-N:end-δn-1,1:N] .= (-1) .* In
-    Jext[end-δn:end, end-δn:end] = Jcond[end-δn:end, end-δn:end]
+    Jext[end-δn-N:end-δn-1, end-δn-N:end-δn-1] .= In
+    Jext[end-δn-N:end-δn-1, 1:N] .= (-1) .* In
+    Jext[end-δn:end, end-δn:end] .= Jcond[end-δn:end, end-δn:end]
     rhs_ext = zeros(𝒯, size(Jext, 1))
 
     # we solve for the external unknowns
     for _ in 1:Ntst
-        Aᵢ .= Jcond[r2, r1]
-        Bᵢ .= Jcond[r2, r1 .+ nbcoll]
-
-        Jext[rN, rN] .= Aᵢ
-        Jext[rN, rN .+ N] .= Bᵢ
+        Jext[rN, rN] .= Jcond[r2, r1]
+        Jext[rN, rN .+ N] .= Jcond[r2, r1 .+ nbcoll]
 
         Jext[rN, end-δn:end] .= Jcond[r2, Npo:(Npo+δn)]
 

test/stuartLandauCollocation.jl

@@ -263,7 +263,7 @@ prob_col       = BK.PeriodicOrbitOCollProblem(Ntst, m; prob_vf = prob_ana,
 _ci = BK.generate_solution(prob_col, t->cos(t) .* ones(N), 2pi);
 Jco_sp = BK.analytical_jacobian_sparse(prob_col, _ci, par_sl);
 Jco = BK.analytical_jacobian(prob_col_dense, _ci, par_sl);
-@assert norminf(Jco - Array(Jco_sp)) == 0
+@assert norminf(Jco - Array(Jco_sp)) < 1e-15
 Jco2 = copy(Jco) |> sparse;
 Jco2 .= 0
 _indx = BifurcationKit.get_blocks(prob_col, Jco2);

test/stuartLandauTrap.jl

@@ -54,15 +54,18 @@ opts_po_cont = ContinuationPar(dsmax = 0.02, ds = 0.001, p_max = 2.2, max_steps
 lsdef = DefaultLS()
 lsit = GMRESKrylovKit()
 for (ind, jacobianPO) in enumerate((:Dense, :DenseAD, :FullLU, :BorderedLU, :FullSparseInplace, :BorderedSparseInplace, :FullMatrixFree, :FullMatrixFreeAD, :BorderedMatrixFree))
-    @show jacobianPO, ind
-    _ls = ind > 5 ? lsit : lsdef
+    _ls = ind > 6 ? lsit : lsdef
+    @info jacobianPO, ind, _ls
     outpo_f = newton((@set poTrap.jacobian = jacobianPO),
         orbitguess_f, (@set optn_po.linsolver = _ls);
         normN = norminf)
     @test BK.converged(outpo_f)
 
-    br_po = continuation((@set poTrap.jacobian = jacobianPO), outpo_f.u,
-        PALC(),    (@set opts_po_cont.newton_options.linsolver = _ls);
+    br_po = continuation(
+        (@set poTrap.jacobian = jacobianPO), 
+        outpo_f.u,
+        PALC(),
+        (@set opts_po_cont.newton_options.linsolver = _ls);
         verbosity = 0, plot = false,
         linear_algo = BorderingBLS(solver = _ls, check_precision = false),
         record_from_solution = (u, p; k...) -> BK.getamplitude(poTrap, u, par_hopf; ratio = 1), normC = norminf)
@@ -76,7 +79,7 @@ let
 
     # computation of the Jacobian at out_pof
     _J1 = poTrap(Val(:JacFullSparse), outpo_f.u, par_hopf)
-    _Jfd = ForwardDiff.jacobian(z-> poTrap(z,par_hopf), outpo_f.u)
+    _Jfd = ForwardDiff.jacobian(z-> BK.residual(poTrap, z, par_hopf), outpo_f.u)
 
     # test of the jacobian against automatic differentiation
     @test norm(_Jfd - Array(_J1), Inf) < 1e-7