GLPK not thread safe?

[rambunctiousapple]

I want to multithread set computations for speed. This example gives the gist:

pupil(diameter) = overapproximate(Ball2([0.0,0.0],diameter),.05)
regular_polygon(numsides::Int) =[[sin(2π*i/numsides),cos(2π*i/numsides)] for i in 0:numsides-1]

function test(numtimes, eye)
    display = VPolygon(regular_polygon(6))
    result = Vector(undef,0)

    Threads.@threads for i in 1:numtimes
        push!(result,area(LazySets.translate(display,[i,0])∩eye))
    end
    return result
end

This code works fine if I don't multithread the for loop but it is very slow: 14ms on my core i7 11800H 8 core laptop. But the multithreaded code crashes:

julia> eye = pupil(4.0); test(6,eye)
glp_free: memory allocation error
Error detected in file env/alloc.c at line 72
error in running finalizer: GLPK.GLPKFatalError(msg="GLPK call failed. All GLPK objects you defined so far are now invalidated.")
_err_hook at C:\Users\seatt\.julia\packages\GLPK\oTTtu\src\GLPK.jl:223
unknown function (ip: 0000000000f43634)
errfunc at /workspace/srcdir/glpk-4.64/src\env\error.c:54
dma at /workspace/srcdir/glpk-4.64/src\env\alloc.c:72
_glp_lufint_delete at /workspace/srcdir/glpk-4.64/src\bflib\lufint.c:162
_glp_fhvint_delete at /workspace/srcdir/glpk-4.64/src\bflib\fhvint.c:157
_glp_bfd_delete_it at /workspace/srcdir/glpk-4.64/src\draft\bfd.c:528
glp_delete_prob at /workspace/srcdir/glpk-4.64/src\api\prob1.c:1583
delete_prob at C:\Users\seatt\.julia\packages\GLPK\oTTtu\src\GLPK.jl:229
jl_apply at /cygdrive/c/buildbot/worker/package_win64/build/src\julia.h:1788 [inlined]
run_finalizer at /cygdrive/c/buildbot/worker/package_win64/build/src\gc.c:278
jl_gc_run_finalizers_in_list at /cygdrive/c/buildbot/worker/package_win64/build/src\gc.c:365
run_finalizers at /cygdrive/c/buildbot/worker/package_win64/build/src\gc.c:394 [inlined]
run_finalizers at /cygdrive/c/buildbot/worker/package_win64/build/src\gc.c:372
jl_gc_run_pending_finalizers at /cygdrive/c/buildbot/worker/package_win64/build/src\gc.c:405
jl_mutex_unlock at /cygdrive/c/buildbot/worker/package_win64/build/src\julia_locks.h:131 [inlined]
jl_generate_fptr at /cygdrive/c/buildbot/worker/package_win64/build/src\jitlayers.cpp:359
jl_compile_method_internal at /cygdrive/c/buildbot/worker/package_win64/build/src\gf.c:1980
jl_compile_method_internal at /cygdrive/c/buildbot/worker/package_win64/build/src\gf.c:1934 [inlined]
_jl_invoke at /cygdrive/c/buildbot/worker/package_win64/build/src\gf.c:2239 [inlined]
jl_apply_generic at /cygdrive/c/buildbot/worker/package_win64/build/src\gf.c:2429
display at .\multimedia.jl:328
jl_apply at /cygdrive/c/buildbot/worker/package_win64/build/src\julia.h:1788 [inlined]
jl_f__call_latest at /cygdrive/c/buildbot/worker/package_win64/build/src\builtins.c:757
#invokelatest#2 at .\essentials.jl:716 [inlined]
invokelatest at .\essentials.jl:714 [inlined]
print_response at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.7\REPL\src\REPL.jl:291
#45 at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.7\REPL\src\REPL.jl:275
unknown function (ip: 0000000068e5cfd6)
with_repl_linfo at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.7\REPL\src\REPL.jl:508
print_response at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.7\REPL\src\REPL.jl:273
do_respond at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.7\REPL\src\REPL.jl:844
unknown function (ip: 0000000068e46070)
jl_apply at /cygdrive/c/buildbot/worker/package_win64/build/src\julia.h:1788 [inlined]
jl_f__call_latest at /cygdrive/c/buildbot/worker/package_win64/build/src\builtins.c:757
#invokelatest#2 at .\essentials.jl:716 [inlined]
invokelatest at .\essentials.jl:714 [inlined]
run_interface at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.7\REPL\src\LineEdit.jl:2493
unknown function (ip: 0000000060d64d6f)
run_frontend at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.7\REPL\src\REPL.jl:1230
#49 at .\task.jl:423
unknown function (ip: 0000000060d5ff93)
jl_apply at /cygdrive/c/buildbot/worker/package_win64/build/src\julia.h:1788 [inlined]
start_task at /cygdrive/c/buildbot/worker/package_win64/build/src\task.c:877

Is the GLPK library not thread safe? Is it possible to multithread operations in LazySets?

[schillic]

We have not really¹ considered multithreading so far, so I am not surprised that there are problems.

Here is what happens. The area function computes the concrete intersection via intersection. That function dispatches on the set types, which here are VPolygon and HPolygon, and thus chooses a method for general polytopes (maybe this should be changed, but not the point here). The problem is then that we remove redundant constraints (by default) with an LP solver, for which we use GLPK (by default). And as you observed, that fails.

As a remedy you can use the intersection method for polygons, which is chosen automatically if the display set is an HPolygon. You can convert it via tohrep. So a thread-safe version could be:

function test(numtimes, eye)
    display = tohrep(VPolygon(regular_polygon(6)))
    result = Vector(undef,0)

    Threads.@threads for i in 1:numtimes
        push!(result,area(LazySets.translate(display,[i,0])∩eye))
    end
    return result
end

Of course that does not fix the problem with GLPK (which you need if you want to work with higher-dimensional sets). For that it really seems that the solver is not thread-safe. If you need a solution for that case, I think it is better to ask the GLPK developers (see their issue #157) and then we can think about a solution.

¹Not totally true, but anyway.

[rambunctiousapple]

Thanks, I'll try using tohrep and look at the GLPK issue #157.

[schillic]

I close this because there is probably nothing we can do. #2935 added a method that does not use GLPK for polygons anymore.

[mforets]

I want to multithread set computations for speed.

For speed, my suggestion would be that you consider using statically sized arrays, instead of normal Julia arrays (see StaticArrays.jl).