Skip to content

Add --enable-openmp to exercise ctlgeom thread-safety in test suite #79

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
stevengj merged 1 commit into from
May 18, 2026
Merged

Add --enable-openmp to exercise ctlgeom thread-safety in test suite #79

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 12 additions & 0 deletions configure.ac
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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])],
Expand Down
123 changes: 123 additions & 0 deletions m4/ax_openmp.m4
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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 <stevenj@alum.mit.edu>
# Copyright (c) 2015 John W. Peterson <jwpeterson@gmail.com>
# Copyright (c) 2016 Nick R. Papior <nickpapior@gmail.com>
#
# 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 <https://www.gnu.org/licenses/>.
#
# 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 <omp.h>

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
])
51 changes: 37 additions & 14 deletions utils/geomtst.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,10 @@
#include <math.h>
#include <time.h>

#ifdef _OPENMP
#include <omp.h>
#endif

#include "ctlgeom.h"

/************************************************************************/
Expand Down Expand Up @@ -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;
}

Expand Down Expand Up @@ -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;
Expand All @@ -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;

Expand All @@ -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;
Expand Down
36 changes: 25 additions & 11 deletions utils/test-prism.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,10 @@
#include <time.h>
#include <string.h>

#ifdef _OPENMP
#include <omp.h>
#endif

#include "ctlgeom.h"

vector3 normal_to_plane(vector3 o, vector3 v1, vector3 v2, vector3 p, double *min_distance);
Expand Down Expand Up @@ -274,23 +278,26 @@ 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);
boolean in_prism = point_in_objectp(p, 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);

Expand All @@ -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.
Expand All @@ -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);

Expand All @@ -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);
Expand All @@ -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);
Expand Down