Mapping with openmp segmentation faults sometimes

[ernestds]

import casadi as ca
import numpy as np
import time

num_samples = 10
num_positions = 2
map_size = 2

module = ca.MX
weights_matrix_sym = module.sym('A',1,num_samples) # gets tranposed later
support_vectors_matrix_sym = module.sym('S',num_samples,num_positions)


configuration = module.sym('q',num_positions)


fk_support_vectors_matrix_sym = module.sym('fk_S',12,num_samples)
support_vector =  module.sym('q',num_positions)
case = 1
# Case 1: seg faults
# Case 2: doesn't seg fault
# Case 3: doesn't seg fault
# Case 4: doesn't seg fault
# Case 5: sometimes seg faults
# Case 6: doesn't seg fault
polyharmonic_kernel_function = ca.Function('func_name',[support_vector,fk_support_vectors_matrix_sym],[ca.norm_2((ca.veccat(support_vector,fk_support_vectors_matrix_sym)))])
if case == 1:
  score = ca.sumsqr(polyharmonic_kernel_function.call([configuration,fk_support_vectors_matrix_sym],False)[0])
  score_function = ca.Function('score',[configuration,fk_support_vectors_matrix_sym,weights_matrix_sym],[(score)],['q_sample','fk_S','A'],['g'],{'jit':False,}) 
if case == 2:
  configuration = support_vector
  score = ca.sumsqr(polyharmonic_kernel_function.call([configuration,fk_support_vectors_matrix_sym],False)[0])
  score_function = ca.Function('score',[configuration,fk_support_vectors_matrix_sym,weights_matrix_sym],[(score)],['q_sample','fk_S','A'],['g'],{'jit':False,}) 
if case == 3:
  
  score = ca.sumsqr(polyharmonic_kernel_function.call([configuration,fk_support_vectors_matrix_sym],True)[0])
  score_function = ca.Function('score',[configuration,fk_support_vectors_matrix_sym,weights_matrix_sym],[(score)],['q_sample','fk_S','A'],['g'],{'jit':False,}) 
if case == 4:
  score = ca.sumsqr(polyharmonic_kernel_function.call([configuration,fk_support_vectors_matrix_sym],False)[0])
  score_function = ca.Function('score',[configuration,fk_support_vectors_matrix_sym,weights_matrix_sym],[(score)],['q_sample','fk_S','A'],['g'],{'jit':True,'jit_options':{'flags':['-Ofast']}}) 
if case == 5:
  score = ca.sumsqr(polyharmonic_kernel_function.call([configuration,fk_support_vectors_matrix_sym],False)[0])
  score_function = ca.Function('score',[configuration,fk_support_vectors_matrix_sym,weights_matrix_sym],[(score)],['q_sample','fk_S','A'],['g'],{'jit':False,}) 
if case == 6:
  score = ca.sumsqr(polyharmonic_kernel_function.call([configuration,fk_support_vectors_matrix_sym],False)[0])
  score_function = ca.Function('score',[configuration,fk_support_vectors_matrix_sym,weights_matrix_sym],[(score)],['q_sample','fk_S','A'],['g'],{'jit':False,}) 
if case == 6:
  func = score_function.map(12,)
else:
  func = score_function.map(12,'openmp')

[buf,f_eval] = func.buffer()
res = []
for i in range(func.n_out()):
    res.append(np.zeros(func.sparsity_out(i).shape))
    buf.set_res(i, memoryview(res[-1]))
args = []
for i in range(func.n_in()):
    args.append((np.random.randn(*func.sparsity_in(i).shape)))
    buf.set_arg(i, memoryview(args[-1]))
times = []
r = 100
n = 100
for it in range(r):
    t0 = time.time()
    for i in range(n):
        if case == 5: time.sleep(0.01)
        f_eval()
        # print('aaa')
        # func.call(args)
    print(it)
    times.append(time.time()-t0)
times = np.array(times)
print(f'Mean: {times.mean()/n}; standard deviation: {times.std()/n}')

Reproducible in my system as above. Ubuntu 22.04 in Docker over WSL with Casadi 3.6.5 compiled with:
-DWITH_HIGHS=OFF;-DWITH_BUILD_MUMPS=ON;-DWITH_BUILD_HIGHS=OFF;-DWITH_BUILD_SPRAL=OFF;-DWITH_SPRAL=OFF;-DWITH_PROXQP=OFF;-DWITH_BUILD_PROXQP=OFF;-DWITH_BUILD_EIGEN3=ON;-DWITH_BUILD_SIMDE=OFF;-DWITH_SLEQP=OFF;-DWITH_BUILD_SLEQP=OFF;-DWITH_BUILD_ALPAQA=OFF;-DWITH_ALPAQA=OFF;-DWITH_BQPD=OFF;-DWITH_UNO=OFF;-DWITH_OPENMP=ON;-DWITH_IPOPT=ON;-DWITH_BUILD_IPOPT=ON;-DWITH_MUMPS=ON;-DWITH_BUILD_LAPACK=ON;-DWITH_LAPACK=ON;-DWITH_BUILD_BLASFEO=ON;-DWITH_OSQP=ON;-DWITH_HSL=ON;-DWITH_MOCKUP_HSL=ON;-DWITH_BUILD_METIS=ON;-DWITH_BUILD_OSQP=ON

Hopefully the code is understandable. Let me know if I'm actually doing something dumb.