Bug: AssertionError: LDLT.factor != nothing

[tmigot]

I tried to solve the following problem modeled with PDENLPModels.jl

using Gridap, PDENLPModels

  # Definition of the domain
  n = 100
  domain = (-1, 1, -1, 1)
  partition = (n, n)
  model = CartesianDiscreteModel(domain, partition)

  # Definition of the spaces:
  valuetype = Float64
  reffe = ReferenceFE(lagrangian, valuetype, 2)
  Xpde = TestFESpace(model, reffe; conformity = :H1, dirichlet_tags = "boundary")
  y0(x) = 0.0
  Ypde = TrialFESpace(Xpde, y0)

  reffe_con = ReferenceFE(lagrangian, valuetype, 1)
  Xcon = TestFESpace(model, reffe_con; conformity = :H1)
  Ycon = TrialFESpace(Xcon)
  Y = MultiFieldFESpace([Ypde, Ycon])

  # Integration machinery
  trian = Triangulation(model)
  degree = 1
  dΩ = Measure(trian, degree)

  # Objective function:
  yd(x) = -x[1]^2
  α = 1e-2
  function f(y, u)
    ∫(0.5 * (yd - y) * (yd - y) + 0.5 * α * u * u) * dΩ
  end

  # Definition of the constraint operator
  ω = π - 1 / 8
  h(x) = -sin(ω * x[1]) * sin(ω * x[2])
  function res(y, u, v)
    ∫(∇(v) ⊙ ∇(y) - v * u - v * h) * dΩ
  end
  op = FEOperator(res, Y, Xpde)
  npde = Gridap.FESpaces.num_free_dofs(Ypde)
  ncon = Gridap.FESpaces.num_free_dofs(Ycon)
  x0 = zeros(npde + ncon)
  nlp = GridapPDENLPModel(x0, f, trian, Ypde, Ycon, Xpde, Xcon, op, name = "Control elastic membrane")

using CaNNOLeS, NLPModelsModifiers
nls = FeasibilityResidual(nlp)
stats_cannoles = cannoles(nls)

However, I get the following mysterious error:

┌ Warning: linsolve ma57 not available. Using :ldlfactorizations instead
└ @ CaNNOLeS .julia\packages\CaNNOLeS\CGoHr\src\CaNNOLeS.jl:60
[ Info:      I      #F        fx        Δt      ‖∇L‖  ‖Fx - r‖    ‖c(x)‖         α         η         ρ         δ   in_it     nbk  
[ Info:      0       2   2.2e-04   0.0e+00   1.6e-07   0.0e+00   0.0e+00
ERROR: LoadError: AssertionError: LDLT.factor != nothing
Stacktrace:
 [1] newton_system(x::Vector{Float64}, r::Vector{Float64}, λ::Vector{Float64}, Fx::Vector{Float64}, rhs::Vector{Float64}, LDLT::CaNNOLeS.LDLFactStruct, ρold::Float64, params::Dict{Symbol, Float64}, method::Symbol, linsolve::Symbol)
   @ CaNNOLeS .julia\packages\CaNNOLeS\CGoHr\src\CaNNOLeS.jl:616
 [2] cannoles(nls::FeasibilityResidual{Float64, Vector{Float64}}; x::Vector{Float64}, λ::Vector{Float64}, method::Symbol, merit::Symbol, linsolve::Symbol, max_f::Int64, max_time::Float64, max_inner::Int64, ϵtol::Float64, check_small_residual::Bool, always_accept_extrapolation::Bool, ϵkchoice::Symbol, δdec::Float64)
   @ CaNNOLeS .julia\packages\CaNNOLeS\CGoHr\src\CaNNOLeS.jl:315
 [3] cannoles(nls::FeasibilityResidual{Float64, Vector{Float64}})
   @ CaNNOLeS .julia\packages\CaNNOLeS\CGoHr\src\CaNNOLeS.jl:50

Any idea?

My environment is the following:

  [5a1c9e79] CaNNOLeS v0.5.3
  [56d4f2e9] Gridap v0.15.5
  [40e66cde] LDLFactorizations v0.8.1
  [a4795742] NLPModels v0.18.1
  [e01155f1] NLPModelsModifiers v0.5.1
  [80da258d] PDENLPModels v0.3.0

[abelsiqueira]

This assertion. The only way to reach here is if the regularization becomes infinite (rho_max is prev_float(inf)).
This try is a catch-all to declare failure in the regularization attempt. If it is failing infinitely, it could mean that something else in that clause is throwing errors. I.e., something unrelated to LDL. Looking at it again, it can be

CaNNOLeS.jl/src/CaNNOLeS.jl

Line 561 in 23116b0

A = Matrix(Symmetric(A, :L))

which is the conversion to Matrix that ldl needed for some reason. Out of memory error, maybe?
Do you have HSL to test?

[tmigot]

So, I move forward with this. The same error appears for a couple of problems from NLSProblems as well. For instance:

using CaNNOLeS, NLSProblems
nlp = hs06()
cannoles(nlp, verbose = 1)

The error doesn't come here from the ldlt factorization, but from here

CaNNOLeS.jl/src/CaNNOLeS.jl

Line 644 in a19be83

A = sparse(LDLT.rows, LDLT.cols, LDLT.vals, N, N)

(LDLT.rows, LDLT.cols, LDLT.vals) = ([1, 1, 3, 4, 4, 3, 4, 1, 2], [1, 5108070049, 1, 1, 2, 3, 4, 1, 2], [0.0, -1.5621301775147933, -1.0, 24.0, 10.0, -1.0, -0.1, 0.0, 0.0])
ERROR: ArgumentError: column indices J[k] must satisfy 1 <= J[k] <= n

[tmigot]

Oh, I think the problem comes from the model actually. I did the following (a bit nervous) test:

julia> using NLPModels

julia> hess_structure(nlp)
([1, 267043088], [1, 167116808])

julia> hess_structure(nlp)
([1, 331019536], [1, 269359968])

julia> hess_structure(nlp)
([1, 331019536], [1, 167149984])

julia> hess_structure(nlp)
([1, 331019536], [1, 167150048])

julia> hess_structure(nlp)
([1, 5124065472], [1, 5124065472])