From c323e2df7a4e6c7300e62e5282895c2440b61c5c Mon Sep 17 00:00:00 2001 From: Luochenghuang Date: Fri, 15 May 2026 22:47:18 -0700 Subject: [PATCH] Add --enable-openmp to exercise ctlgeom thread-safety in test suite Adds m4/ax_openmp.m4 (autoconf-archive AX_OPENMP) and a --enable-openmp configure flag that appends OPENMP_CFLAGS to CFLAGS/LDFLAGS. geomtst.c and test-prism.c are wrapped with #pragma omp parallel for on the random-point/random-segment test loops so the ctlgeom entry points they exercise are hammered from multiple threads concurrently. random_object_and_lattice() in geomtst.c now takes a perturb_lattice flag; the parallel loops call it with 0 so they don't race on the global geometry_lattice (perturbed once in serial up-front instead). test-prism's file writes are guarded with #pragma omp critical, and counters use reduction(+:...). Default ./configure leaves OpenMP off and make check is unchanged. With --enable-openmp, make check will surface races in ctlgeom (notably reinit_prism via geom_fix_object_ptr). --- configure.ac | 12 +++++ m4/ax_openmp.m4 | 123 +++++++++++++++++++++++++++++++++++++++++++++ utils/geomtst.c | 51 +++++++++++++------ utils/test-prism.c | 36 +++++++++---- 4 files changed, 197 insertions(+), 25 deletions(-) create mode 100644 m4/ax_openmp.m4 diff --git a/configure.ac b/configure.ac index aabdf73..614f0fb 100644 --- a/configure.ac +++ b/configure.ac @@ -52,6 +52,18 @@ fi ########################################################################### +# Optional OpenMP support, used to exercise thread-safety of ctlgeom +# functions in the test suite via #pragma omp parallel loops. +AC_ARG_ENABLE(openmp, [ --enable-openmp enable OpenMP for tests], + enable_openmp=$enableval, enable_openmp=no) +if test "x$enable_openmp" = xyes; then + AX_OPENMP([CFLAGS="$CFLAGS $OPENMP_CFLAGS" + LDFLAGS="$LDFLAGS $OPENMP_CFLAGS"], + [AC_MSG_ERROR([--enable-openmp specified but no OpenMP support detected])]) +fi + +########################################################################### + # Find Guile library, flags, etcetera: AC_ARG_WITH(guile, [AS_HELP_STRING([--without-guile],[disable use of Guile])], diff --git a/m4/ax_openmp.m4 b/m4/ax_openmp.m4 new file mode 100644 index 0000000..9868fd0 --- /dev/null +++ b/m4/ax_openmp.m4 @@ -0,0 +1,123 @@ +# =========================================================================== +# https://www.gnu.org/software/autoconf-archive/ax_openmp.html +# =========================================================================== +# +# SYNOPSIS +# +# AX_OPENMP([ACTION-IF-FOUND[, ACTION-IF-NOT-FOUND]]) +# +# DESCRIPTION +# +# This macro tries to find out how to compile programs that use OpenMP a +# standard API and set of compiler directives for parallel programming +# (see http://www-unix.mcs/) +# +# On success, it sets the OPENMP_CFLAGS/OPENMP_CXXFLAGS/OPENMP_F77FLAGS +# output variable to the flag (e.g. -omp) used both to compile *and* link +# OpenMP programs in the current language. +# +# NOTE: You are assumed to not only compile your program with these flags, +# but also link it with them as well. +# +# If you want to compile everything with OpenMP, you should set: +# +# CFLAGS="$CFLAGS $OPENMP_CFLAGS" +# #OR# CXXFLAGS="$CXXFLAGS $OPENMP_CXXFLAGS" +# #OR# FFLAGS="$FFLAGS $OPENMP_FFLAGS" +# +# (depending on the selected language). +# +# The user can override the default choice by setting the corresponding +# environment variable (e.g. OPENMP_CFLAGS). +# +# ACTION-IF-FOUND is a list of shell commands to run if an OpenMP flag is +# found, and ACTION-IF-NOT-FOUND is a list of commands to run it if it is +# not found. If ACTION-IF-FOUND is not specified, the default action will +# define HAVE_OPENMP. +# +# LICENSE +# +# Copyright (c) 2008 Steven G. Johnson +# Copyright (c) 2015 John W. Peterson +# Copyright (c) 2016 Nick R. Papior +# +# This program is free software: you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation, either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program. If not, see . +# +# As a special exception, the respective Autoconf Macro's copyright owner +# gives unlimited permission to copy, distribute and modify the configure +# scripts that are the output of Autoconf when processing the Macro. You +# need not follow the terms of the GNU General Public License when using +# or distributing such scripts, even though portions of the text of the +# Macro appear in them. The GNU General Public License (GPL) does govern +# all other use of the material that constitutes the Autoconf Macro. +# +# This special exception to the GPL applies to versions of the Autoconf +# Macro released by the Autoconf Archive. When you make and distribute a +# modified version of the Autoconf Macro, you may extend this special +# exception to the GPL to apply to your modified version as well. + +#serial 13 + +AC_DEFUN([AX_OPENMP], [ +AC_PREREQ([2.69]) dnl for _AC_LANG_PREFIX + +AC_CACHE_CHECK([for OpenMP flag of _AC_LANG compiler], ax_cv_[]_AC_LANG_ABBREV[]_openmp, [save[]_AC_LANG_PREFIX[]FLAGS=$[]_AC_LANG_PREFIX[]FLAGS +ax_cv_[]_AC_LANG_ABBREV[]_openmp=unknown +# Flags to try: -fopenmp (gcc), -mp (SGI & PGI), +# -qopenmp (icc>=15), -openmp (icc), +# -xopenmp (Sun), -omp (Tru64), +# -qsmp=omp (AIX), +# none +ax_openmp_flags="-fopenmp -openmp -qopenmp -mp -xopenmp -omp -qsmp=omp none" +if test "x$OPENMP_[]_AC_LANG_PREFIX[]FLAGS" != x; then + ax_openmp_flags="$OPENMP_[]_AC_LANG_PREFIX[]FLAGS $ax_openmp_flags" +fi +for ax_openmp_flag in $ax_openmp_flags; do + case $ax_openmp_flag in + none) []_AC_LANG_PREFIX[]FLAGS=$save[]_AC_LANG_PREFIX[] ;; + *) []_AC_LANG_PREFIX[]FLAGS="$save[]_AC_LANG_PREFIX[]FLAGS $ax_openmp_flag" ;; + esac + AC_LINK_IFELSE([AC_LANG_SOURCE([[ +@%:@include + +static void +parallel_fill(int * data, int n) +{ + int i; +@%:@pragma omp parallel for + for (i = 0; i < n; ++i) + data[i] = i; +} + +int +main(void) +{ + int arr[100000]; + omp_set_num_threads(2); + parallel_fill(arr, 100000); + return 0; +} +]])],[ax_cv_[]_AC_LANG_ABBREV[]_openmp=$ax_openmp_flag; break],[]) +done +[]_AC_LANG_PREFIX[]FLAGS=$save[]_AC_LANG_PREFIX[]FLAGS +]) +if test "x$ax_cv_[]_AC_LANG_ABBREV[]_openmp" = "xunknown"; then + m4_default([$2],:) +else + if test "x$ax_cv_[]_AC_LANG_ABBREV[]_openmp" != "xnone"; then + OPENMP_[]_AC_LANG_PREFIX[]FLAGS=$ax_cv_[]_AC_LANG_ABBREV[]_openmp + fi + m4_default([$1], [AC_DEFINE(HAVE_OPENMP,1,[Define if OpenMP is enabled])]) +fi +]) diff --git a/utils/geomtst.c b/utils/geomtst.c index da24aed..e6bdc33 100644 --- a/utils/geomtst.c +++ b/utils/geomtst.c @@ -3,6 +3,10 @@ #include #include +#ifdef _OPENMP +#include +#endif + #include "ctlgeom.h" /************************************************************************/ @@ -145,15 +149,18 @@ static double simple_ellip_overlap(geom_box b, geometric_object o, double tol) { return olap0; } -geometric_object random_object_and_lattice(void) { +/* If perturb_lattice is nonzero, this function mutates the global + geometry_lattice variable. When called from a parallel (e.g. OpenMP) + context, pass perturb_lattice=0 to avoid the data race. */ +geometric_object random_object_and_lattice(int perturb_lattice) { geometric_object o = random_object(); -#if 1 - geometry_lattice.basis1 = random_unit_vector3(); - geometry_lattice.basis2 = random_unit_vector3(); - geometry_lattice.basis3 = random_unit_vector3(); - geom_fix_lattice(); + if (perturb_lattice) { + geometry_lattice.basis1 = random_unit_vector3(); + geometry_lattice.basis2 = random_unit_vector3(); + geometry_lattice.basis3 = random_unit_vector3(); + geom_fix_lattice(); + } geom_fix_object_ptr(&o); -#endif return o; } @@ -195,11 +202,11 @@ void check_overlap(double tol, double olap0, double olap, int dim, geometric_obj printf("Got %s %dd overlap %g vs. %g with tol = %e\n", object_name(o), dim, olap, olap0, tol); } -static void test_overlap(double tol, +static void test_overlap(double tol, int perturb_lattice, number (*box_overlap_with_object)(geom_box b, geometric_object o, number tol, integer maxeval), double (*simple_overlap)(geom_box b, geometric_object o, double tol)) { - geometric_object o = random_object_and_lattice(); + geometric_object o = random_object_and_lattice(perturb_lattice); vector3 dir = random_unit_vector3(); geom_box b; double d, olap, olap0; @@ -221,8 +228,8 @@ static void test_overlap(double tol, geometric_object_destroy(o); } -static void test_volume(double tol) { - geometric_object o = random_object_and_lattice(); +static void test_volume(double tol, int perturb_lattice) { + geometric_object o = random_object_and_lattice(perturb_lattice); geom_box b; double olap1, olap2; @@ -240,19 +247,35 @@ int main(void) { const int ntest = 100; const double tol = 1e-2; int i; + int perturb_lattice = 1; srand(time(NULL)); geom_initialize(); +#ifdef _OPENMP + /* Perturb the global geometry_lattice once in serial; the parallel loops + below must not race on it, so they pass perturb_lattice=0 to the test + helpers. */ + geometry_lattice.basis1 = random_unit_vector3(); + geometry_lattice.basis2 = random_unit_vector3(); + geometry_lattice.basis3 = random_unit_vector3(); + geom_fix_lattice(); + perturb_lattice = 0; + printf("**** whole box overlap (OpenMP, %d threads): ****\n", omp_get_max_threads()); +#else printf("**** whole box overlap: ****\n"); +#endif + +#pragma omp parallel for schedule(dynamic) for (i = 0; i < ntest; ++i) - test_volume(tol); + test_volume(tol, perturb_lattice); +#pragma omp parallel for schedule(dynamic) for (i = 0; i < ntest; ++i) { printf("**** box overlap: ****\n"); - test_overlap(tol, box_overlap_with_object, simple_overlap); + test_overlap(tol, perturb_lattice, box_overlap_with_object, simple_overlap); printf("**** ellipsoid overlap: ****\n"); - test_overlap(tol, ellipsoid_overlap_with_object, simple_ellip_overlap); + test_overlap(tol, perturb_lattice, ellipsoid_overlap_with_object, simple_ellip_overlap); } return 0; diff --git a/utils/test-prism.c b/utils/test-prism.c index f88692f..af89c16 100644 --- a/utils/test-prism.c +++ b/utils/test-prism.c @@ -29,6 +29,10 @@ #include #include +#ifdef _OPENMP +#include +#endif + #include "ctlgeom.h" vector3 normal_to_plane(vector3 o, vector3 v1, vector3 v2, vector3 p, double *min_distance); @@ -274,6 +278,9 @@ int test_point_inclusion(geometric_object the_block, geometric_object the_prism, vector3 max_corner = vector3_scale(+1.0, size); FILE *f = write_log ? fopen("/tmp/test-prism.points", "w") : 0; int num_failed = 0, num_adjusted = 0, n; + char *exclude_env = getenv("LIBCTL_EXCLUDE_BOUNDARIES"); + boolean libctl_include_boundaries = !(exclude_env && exclude_env[0] == '1'); +#pragma omp parallel for schedule(dynamic) reduction(+ : num_failed, num_adjusted) for (n = 0; n < num_tests; n++) { vector3 p = random_point_in_box(min_corner, max_corner); boolean in_block = point_in_objectp(p, the_block); @@ -281,16 +288,16 @@ int test_point_inclusion(geometric_object the_block, geometric_object the_prism, if (in_block != in_prism) { // retry with boundary exclusion/inclusion reversed - boolean libctl_include_boundaries = 1; - char *s = getenv("LIBCTL_EXCLUDE_BOUNDARIES"); - if (s && s[0] == '1') libctl_include_boundaries = 0; in_prism = point_in_or_on_prism(the_prism.subclass.prism_data, p, 1 - libctl_include_boundaries); if (in_block == in_prism) num_adjusted++; } if (in_block != in_prism) num_failed++; - if (f) fprintf(f, "%i %i %e %e %e \n", in_block, in_prism, p.x, p.y, p.z); + if (f) { +#pragma omp critical(test_prism_log) + fprintf(f, "%i %i %e %e %e \n", in_block, in_prism, p.x, p.y, p.z); + } } if (f) fclose(f); @@ -314,6 +321,7 @@ int test_normal_to_object(geometric_object the_block, geometric_object the_prism double tolerance = 1.0e-6; int n; +#pragma omp parallel for schedule(dynamic) reduction(+ : num_failed) for (n = 0; n < num_tests; n++) { // with probability PFACE, generate random base point lying on one // of the 6 faces of the prism. @@ -326,10 +334,12 @@ int test_normal_to_object(geometric_object the_block, geometric_object the_prism vector3 nhat_prism = standardize(normal_to_object(p, the_prism)); if (!vector3_nearly_equal(nhat_block, nhat_prism, tolerance)) num_failed++; - if (f) + if (f) { +#pragma omp critical(test_prism_log) fprintf(f, "%e %e %e %e %e %e %e %e %e %i\n\n\n", p.x, p.y, p.z, nhat_block.x, nhat_block.y, nhat_block.z, nhat_prism.x, nhat_prism.y, nhat_prism.z, vector3_nearly_equal(nhat_block, nhat_prism, tolerance)); + } } if (f) fclose(f); @@ -349,6 +359,7 @@ int test_line_segment_intersection(geometric_object the_block, geometric_object int num_failed = 0; int n; +#pragma omp parallel for schedule(dynamic) reduction(+ : num_failed) for (n = 0; n < num_tests; n++) { // randomly generated base point within enlarged bounding box vector3 p = random_point_in_box(min_corner, max_corner); @@ -362,13 +373,16 @@ int test_line_segment_intersection(geometric_object the_block, geometric_object if (f) { int success = fabs(sblock - sprism) <= 1.0e-6 * fmax(fabs(sblock), fabs(sprism)); - fprintf(f, " %e %e %s\n", sblock, sprism, success ? "success" : "fail"); - if (success == 0) { - fprintf(f, "#%e %e %e %e %e %e %e %e\n", p.x, p.y, p.z, d.x, d.y, d.z, a, b); - fprintf(f, "%e %e %e\n%e %e %e\n%e %e %e\n", p.x, p.y, p.z, p.x + a * d.x, p.y + a * d.y, - p.z + a * d.z, p.x + b * d.x, p.y + b * d.y, p.z + b * d.z); +#pragma omp critical(test_prism_log) + { + fprintf(f, " %e %e %s\n", sblock, sprism, success ? "success" : "fail"); + if (success == 0) { + fprintf(f, "#%e %e %e %e %e %e %e %e\n", p.x, p.y, p.z, d.x, d.y, d.z, a, b); + fprintf(f, "%e %e %e\n%e %e %e\n%e %e %e\n", p.x, p.y, p.z, p.x + a * d.x, p.y + a * d.y, + p.z + a * d.z, p.x + b * d.x, p.y + b * d.y, p.z + b * d.z); + } + fprintf(f, "\n"); } - fprintf(f, "\n"); } } if (f) fclose(f);