commit 2a7b25faa69b6261e2a04b78ec8870f731974a87
parent c19561bc509f9730aee1017b19f7f302c5fd0ece
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date: Thu, 31 Mar 2016 10:11:25 +0200
Move from handmade openMP parallelization to Star-MC
Stop using rsys/logger too
Diffstat:
| M | cmake/CMakeLists.txt | | | 20 | +++++++------------- |
| A | src/realization.c | | | 140 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | src/realization.h | | | 50 | ++++++++++++++++++++++++++++++++++++++++++++++++++ |
| M | src/s4vs.c | | | 235 | ++++++++++++------------------------------------------------------------------- |
4 files changed, 233 insertions(+), 212 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -34,16 +34,12 @@ set(S4VS_SOURCE_DIR ${PROJECT_SOURCE_DIR}/../src/)
################################################################################
# Check dependencies
################################################################################
-find_package(OpenMP)
find_package(RCMake 0.2 REQUIRED)
find_package(RSys 0.2.1 REQUIRED)
find_package(Star3D 0.3 REQUIRED)
find_package(Star3DAW 0.1.2 REQUIRED)
-find_package(StarSP 0.2 REQUIRED)
-
-if(NOT OPENMP_FOUND)
- message(STATUS "No OpenMP support: muti-threading is disabled")
-endif()
+find_package(StarSP 0.3 REQUIRED)
+find_package(StarMC 0.3 REQUIRED)
include_directories(
${RSys_INCLUDE_DIR}
@@ -62,24 +58,22 @@ set(VERSION_MINOR 1)
set(VERSION_PATCH 0)
set(VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH})
-set(S4VS_FILES_SRC s4vs.c)
+set(S4VS_FILES_SRC s4vs.c realization.c)
+set(S4VS_FILES_INC realization.h)
set(S4VS_FILES_DOC COPYING.fr COPYING.en README.md)
# Prepend each file in the `S4VS_FILES_<SRC|INC>' list by the absolute
# path of the directory in which they lie
rcmake_prepend_path(S4VS_FILES_SRC ${S4VS_SOURCE_DIR})
+rcmake_prepend_path(S4VS_FILES_INC ${S4VS_SOURCE_DIR})
rcmake_prepend_path(S4VS_FILES_DOC ${PROJECT_SOURCE_DIR}/../)
-if(OPENMP_FOUND)
- set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
-endif()
-
if(CMAKE_COMPILER_IS_GNUCC)
set(MATH_LIB m)
endif()
-add_executable(s4vs ${S4VS_FILES_SRC})
-target_link_libraries(s4vs RSys Star3D Star3DAW StarSP ${MATH_LIB})
+add_executable(s4vs ${S4VS_FILES_SRC} ${S4VS_FILES_INC})
+target_link_libraries(s4vs RSys Star3D Star3DAW StarSP StarMC ${MATH_LIB})
set_target_properties(s4vs PROPERTIES
VERSION ${VERSION}
SOVERSION ${VERSION_MAJOR})
diff --git a/src/realization.c b/src/realization.c
@@ -0,0 +1,140 @@
+/* Copyright (C) |Meso|Star> 2015-2016 (contact@meso-star.com)
+ *
+ * This software is governed by the CeCILL license under French law and
+ * abiding by the rules of distribution of free software. You can use,
+ * modify and/or redistribute the software under the terms of the CeCILL
+ * license as circulated by CEA, CNRS and INRIA at the following URL
+ * "http://www.cecill.info".
+ *
+ * As a counterpart to the access to the source code and rights to copy,
+ * modify and redistribute granted by the license, users are provided only
+ * with a limited warranty and the software's author, the holder of the
+ * economic rights, and the successive licensors have only limited
+ * liability.
+ *
+ * In this respect, the user's attention is drawn to the risks associated
+ * with loading, using, modifying and/or developing or reproducing the
+ * software by the user in light of its specific status of free software,
+ * that may mean that it is complicated to manipulate, and that also
+ * therefore means that it is reserved for developers and experienced
+ * professionals having in-depth computer knowledge. Users are therefore
+ * encouraged to load and test the software's suitability as regards their
+ * requirements in conditions enabling the security of their systems and/or
+ * data to be ensured and, more generally, to use and operate it in the
+ * same conditions as regards security.
+ *
+ * The fact that you are presently reading this means that you have had
+ * knowledge of the CeCILL license and that you accept its terms. */
+
+#include <rsys/float3.h>
+
+#include <star/s3d.h>
+#include <star/ssp.h>
+#include <star/smc.h>
+
+#include "realization.h"
+
+/* Maximum number of failures before an error occurs */
+#define MAX_FAILURES 10
+
+void
+realization(void* length, struct ssp_rng* rng, void* context)
+{
+ struct context* ctx = (struct context*)context;
+ struct s3d_attrib attrib;
+ struct s3d_primitive prim;
+ double lambda;
+ float normal[3], direction[4], origin[3], range[2], st[2];
+ struct s3d_hit hit;
+ int nfailures = 0;
+
+ SMC_DOUBLE(length) = 0;
+
+ /* sample a starting point on the scene surface to start the random walk from */
+ for(;;) {
+ /* sample the scene */
+ float u, v, w;
+ do { u = (float)ssp_rng_canonical(rng); } while(u >= 1);
+ do { v = (float)ssp_rng_canonical(rng); } while(v >= 1);
+ do { w = (float)ssp_rng_canonical(rng); } while(w >= 1);
+ /* get a location: primitive ID and parametric coordinates */
+ S3D(scene_sample(ctx->scene, u, v, w, &prim, st));
+
+ /* retrieve the sampled geometric normal and position */
+ S3D(primitive_get_attrib(&prim, S3D_GEOMETRY_NORMAL, st, &attrib));
+ f3_normalize(normal, attrib.value);
+ S3D(primitive_get_attrib(&prim, S3D_POSITION, st, &attrib));
+ f3_set(origin, attrib.value);
+
+ /* cosine weighted sampling of the hemisphere around the sampled normal */
+ ssp_ran_hemisphere_cos(rng, normal, direction);
+
+ /* trace a ray from the sampled location */
+ range[0] = 1.e-6f; range[1] = FLT_MAX;
+ for(;;) {
+ /* find first intersection with geometry to determine its distance */
+ S3D(scene_trace_ray(ctx->scene, origin, direction, range, &hit));
+ /* handle auto-intersection: hitting the same primitive means auto-intersection */
+ if(S3D_PRIMITIVE_EQ(&hit.prim, &prim)) {
+ /* just increase the minimum distance and retry */
+ range[0] += 1.e-6f;
+ }
+ else break; /* standard case: OK */
+ }
+
+ f3_normalize(hit.normal, hit.normal);
+
+ /* handle precision issues (no intersection means geometry leakage: error) */
+ if(S3D_HIT_NONE(&hit)) {
+ if(++nfailures >= MAX_FAILURES) {
+ /* too many errors: stop */
+ ctx->exit_failure = 1;
+ break;
+ }
+ }
+ else break; /* standard case: OK */
+ }
+ if(ctx->exit_failure)
+ return;
+
+ /* here the starting point and a propagation direction have been sampled */
+ /* and we have determined the distance of the geometry in this direction */
+
+ /* start diffuse random walk */
+ for(;;) {
+ /* sample a free path according to ks */
+ lambda = ssp_ran_exp(rng, ctx->ks);
+ if(lambda >= hit.distance) {
+ /* lambda exceeds the geometry distance: the random walk ends on geometry */
+ SMC_DOUBLE(length) += hit.distance;
+ break; /* the standard way of leaving the loop */
+ } else {
+ /* lambda doesn't exceed the geometry distance */
+ /* the random walk can continue after a diffusion */
+ SMC_DOUBLE(length) += lambda;
+ /* new location after lambda */
+ f3_add(origin, origin, f3_mulf(direction, direction, (float)lambda));
+
+ do {
+ /* sample a new direction */
+ ssp_ran_sphere_uniform(rng, direction);
+ range[0] = 0;
+ /* find the first intersection with the geometry */
+ /* to refresh the geometry's distance */
+ S3D(scene_trace_ray(ctx->scene, origin, direction, range, &hit));
+
+ /* handle precision issues (no intersection means geometry leakage: error) */
+ if(S3D_HIT_NONE(&hit)) {
+ if(++nfailures >= MAX_FAILURES) {
+ /* too many errors: stop */
+ ctx->exit_failure = 1;
+ break;
+ }
+ }
+ /* should be an infinite loop, or there is a geometry leakage */
+ } while(S3D_HIT_NONE(&hit));
+ }
+ }
+
+ ATOMIC_ADD(&ctx->nfailures, nfailures);
+}
diff --git a/src/realization.h b/src/realization.h
@@ -0,0 +1,50 @@
+/* Copyright (C) |Meso|Star> 2015-2016 (contact@meso-star.com)
+ *
+ * This software is governed by the CeCILL license under French law and
+ * abiding by the rules of distribution of free software. You can use,
+ * modify and/or redistribute the software under the terms of the CeCILL
+ * license as circulated by CEA, CNRS and INRIA at the following URL
+ * "http://www.cecill.info".
+ *
+ * As a counterpart to the access to the source code and rights to copy,
+ * modify and redistribute granted by the license, users are provided only
+ * with a limited warranty and the software's author, the holder of the
+ * economic rights, and the successive licensors have only limited
+ * liability.
+ *
+ * In this respect, the user's attention is drawn to the risks associated
+ * with loading, using, modifying and/or developing or reproducing the
+ * software by the user in light of its specific status of free software,
+ * that may mean that it is complicated to manipulate, and that also
+ * therefore means that it is reserved for developers and experienced
+ * professionals having in-depth computer knowledge. Users are therefore
+ * encouraged to load and test the software's suitability as regards their
+ * requirements in conditions enabling the security of their systems and/or
+ * data to be ensured and, more generally, to use and operate it in the
+ * same conditions as regards security.
+ *
+ * The fact that you are presently reading this means that you have had
+ * knowledge of the CeCILL license and that you accept its terms. */
+
+#ifndef REALIZATION_H_
+#define REALIZATION_H_
+
+/* forward definition */
+struct ssp_rng;
+
+struct context {
+ struct s3d_sampler* sampler;
+ struct s3d_scene* scene;
+ double ks;
+ int nfailures;
+ char exit_failure;
+};
+
+/*******************************************************************************
+ * MC realization function
+ ******************************************************************************/
+
+void
+realization(void* length, struct ssp_rng* rng, void* context);
+
+#endif
diff --git a/src/s4vs.c b/src/s4vs.c
@@ -28,171 +28,42 @@
#include <rsys/clock_time.h>
#include <rsys/float3.h>
-#include <rsys/logger.h>
+#include <rsys/mem_allocator.h>
#include <star/s3d.h>
#include <star/s3daw.h>
-#include <star/ssp.h>
+#include <star/smc.h>
-#include <omp.h>
-
-/* Maximum number of failures before an error occurs */
-#define MAX_FAILURES 10
-
-/*******************************************************************************
- * Helper function
- ******************************************************************************/
-struct context {
- struct s3d_sampler* sampler;
- struct s3d_scene* scene;
- double ks;
- double sum;
- double sum2;
- int nfailures;
- char exit_failure;
-};
-
-static void
-realization(struct ssp_rng* rng, struct context* ctx)
-{
- struct s3d_attrib attrib;
- struct s3d_primitive prim;
- double lambda, sum = 0;
- float normal[3], direction[4], origin[3], range[2], st[2];
- struct s3d_hit hit;
- int nfailures = 0;
-
- /* sample a starting point on the scene surface to start the random walk from */
- for(;;) {
- /* sample the scene */
- float u, v, w;
- do { u = (float)ssp_rng_canonical(rng); } while(u >= 1);
- do { v = (float)ssp_rng_canonical(rng); } while(v >= 1);
- do { w = (float)ssp_rng_canonical(rng); } while(w >= 1);
- /* get a location: primitive ID and parametric coordinates */
- S3D(scene_sample(ctx->scene, u, v, w, &prim, st));
-
- /* retrieve the sampled geometric normal and position */
- S3D(primitive_get_attrib(&prim, S3D_GEOMETRY_NORMAL, st, &attrib));
- f3_normalize(normal, attrib.value);
- S3D(primitive_get_attrib(&prim, S3D_POSITION, st, &attrib));
- f3_set(origin, attrib.value);
-
- /* cosine weighted sampling of the hemisphere around the sampled normal */
- ssp_ran_hemisphere_cos(rng, normal, direction);
-
- /* trace a ray from the sampled location */
- range[0] = 1.e-6f; range[1] = FLT_MAX;
- for(;;) {
- /* find first intersection with geometry to determine its distance */
- S3D(scene_trace_ray(ctx->scene, origin, direction, range, &hit));
- /* handle auto-intersection: hitting the same primitive means auto-intersection */
- if(S3D_PRIMITIVE_EQ(&hit.prim, &prim)) {
- /* just increase the minimum distance and retry */
- range[0] += 1.e-6f;
- }
- else break; /* standard case: OK */
- }
-
- f3_normalize(hit.normal, hit.normal);
-
- /* handle precision issues (no intersection means geometry leakage: error) */
- if(S3D_HIT_NONE(&hit)) {
- if(++nfailures >= MAX_FAILURES) {
- /* too many errors: stop */
- ctx->exit_failure = 1;
- break;
- }
- }
- else break; /* standard case: OK */
- }
- if(ctx->exit_failure)
- return;
-
- /* here the starting point and a propagation direction have been sampled */
- /* and we have determined the distance of the geometry in this direction */
-
- /* start diffuse random walk */
- for(;;) {
- /* sample a free path according to ks */
- lambda = ssp_ran_exp(rng, ctx->ks);
- if(lambda >= hit.distance) {
- /* lambda exceeds the geometry distance: the random walk ends on geometry */
- sum += hit.distance;
- break;
- } else {
- /* lambda doesn't exceed the geometry distance */
- /* the random walk can continue after a diffusion */
- sum += lambda;
- /* new location after lambda */
- f3_add(origin, origin, f3_mulf(direction, direction, (float)lambda));
-
- do {
- /* sample a new direction */
- ssp_ran_sphere_uniform(rng, direction);
- range[0] = 0;
- /* find the first intersection with the geometry */
- /* to refresh the geometry's distance */
- S3D(scene_trace_ray(ctx->scene, origin, direction, range, &hit));
-
- /* handle precision issues (no intersection means geometry leakage: error) */
- if(S3D_HIT_NONE(&hit)) {
- if(++nfailures >= MAX_FAILURES) {
- /* too many errors: stop */
- ctx->exit_failure = 1;
- break;
- }
- }
- } while(S3D_HIT_NONE(&hit));
- }
- }
- /* accumulate results */
-#ifdef _OPENMP
-#pragma omp atomic update
- ctx->sum += sum;
-#pragma omp atomic update
- ctx->sum2 += sum * sum;
-#pragma omp atomic update
- ctx->nfailures += nfailures;
-#else
- ctx->sum += sum;
- ctx->sum2 += sum * sum;
- ctx->nfailures += nfailures;
-#endif
-}
+#include "realization.h"
static res_T
compute_4v_s(struct s3d_scene* scene, const size_t max_steps, const double ks)
{
char buf[512];
struct time begin, end, elapsed;
- struct s3d_shape* shape;
struct context ctx;
- double length, var, sig;
float S, V, expected;
- struct ssp_rng* rng;
-#ifdef _OPENMP
- struct ssp_rng_proxy* proxy;
-#endif
- unsigned i;
+ struct smc_device* smc;
+ struct smc_integrator integrator;
+ struct smc_estimator* estimator = NULL;
+ struct smc_estimator_status estimator_status;
+ double length;
+ double sig_length;
res_T res = RES_OK;
ASSERT(max_steps > 0 && ks > 0);
- S3D(scene_instantiate(scene, &shape));
S3D(scene_begin_session(scene, S3D_SAMPLE|S3D_TRACE));
/* compute the expected result using a mesh-based method */
S3D(scene_compute_area(scene, &S));
S3D(scene_compute_volume(scene, &V));
if(eq_epsf(S, 0, 1.e-6f) || S < 0) {
- logger_print(LOGGER_DEFAULT, LOG_ERROR,
- "No surface to sample. Is the scene empty?\n");
+ printf("No surface to sample. Is the scene empty?\n");
goto error;
}
if(eq_epsf(V, 0, 1.e-6f) || V < 0) {
- logger_print(LOGGER_DEFAULT, LOG_ERROR,
- "Invalid volume \"%.2f\". The scene might not match the prerequisites:\n"
+ printf("Invalid volume \"%.2f\". The scene might not match the prerequisites:\n"
"it must be closed and its normals must point *into* the volume.\n", V);
goto error;
}
@@ -203,42 +74,16 @@ compute_4v_s(struct s3d_scene* scene, const size_t max_steps, const double ks)
ctx.nfailures = 0;
ctx.exit_failure = 0;
ctx.ks = ks;
- ctx.sum = 0;
- ctx.sum2 = 0;
-
- time_current(&begin);
-#ifdef _OPENMP
-#pragma omp parallel private(rng)
- {
- size_t me = (size_t)omp_get_thread_num();
- /* create a RNG proxy and a derived RNG by thread */
-#pragma omp single
- {
- size_t nb_t = (size_t)omp_get_num_threads();
- SSP(rng_proxy_create(NULL, &ssp_rng_threefry, nb_t, &proxy));
- }
- SSP(rng_proxy_create_rng(proxy, me, &rng));
-
-#pragma omp for
- for (i = 0; i < max_steps; i++) {
- realization(rng, &ctx);
- }
- SSP(rng_ref_put(rng));
-#pragma omp single
- {
- SSP(rng_proxy_ref_put(proxy));
- }
- }
-#else
- /* single-threaded execution: create a single RNG */
- SSP(rng_create(NULL, &ssp_rng_threefry, &rng));
- for (i = 0; i < max_steps; i++) {
- realization(rng, &ctx);
- }
- SSP(rng_ref_put(rng));
-#endif
+ /* setup Star-MC */
+ integrator.integrand = &realization; /* realization function */
+ integrator.type = &smc_double; /* realization type */
+ integrator.max_steps = max_steps; /* realization count */
+ time_current(&begin);
+ /* solve */
+ SMC(device_create(NULL, NULL, 8, NULL, &smc));
+ SMC(solve(smc, &integrator, &ctx, &estimator));
time_current(&end);
time_sub(&elapsed, &end, &begin);
time_dump(&elapsed, TIME_MIN|TIME_SEC|TIME_MSEC, NULL, buf, sizeof(buf));
@@ -246,36 +91,29 @@ compute_4v_s(struct s3d_scene* scene, const size_t max_steps, const double ks)
S3D(scene_end_session(scene));
if(ctx.exit_failure) {
- logger_print(LOGGER_DEFAULT, LOG_ERROR,
- "Too many failures. The scene might not match the prerequisites: it must be\n"
+ printf("Too many failures. The scene might not match the prerequisites: it must be\n"
"closed and its normals must point *into* the volume.\n");
goto error;
}
- /* compute mean length and standard error from MC results */
- length = ctx.sum / (double)max_steps;
- var = ctx.sum2 / (double)max_steps - length * length;
- /* numerical errors can trick us: if variance is negative pretend standard error is 0 */
- sig = var > 0 ? sqrt(var / (double)max_steps) : 0;
-
- logger_print(LOGGER_DEFAULT, LOG_OUTPUT,
- "\n4V/S = %g ~ %g +/- %g\n# failures: %d/%lu\nElapsed time: %s\n",
- expected, length, sig, ctx.nfailures, max_steps, buf);
+ /* get simulation results */
+ SMC(estimator_get_status(estimator, &estimator_status));
+ length = SMC_DOUBLE(estimator_status.E);
+ sig_length = SMC_DOUBLE(estimator_status.SE);
-#ifdef _OPENMP
- logger_print(LOGGER_DEFAULT, LOG_OUTPUT,
- "On %d threads\n", omp_get_num_threads());
-#endif
+ printf("\n4V/S = %g ~ %g +/- %g\n# failures: %d/%lu\nElapsed time: %s\n",
+ expected, length, sig_length, ctx.nfailures, max_steps, buf);
exit:
/* clean data */
+ SMC(device_ref_put(smc));
+ SMC(estimator_ref_put(estimator));
if(scene) {
int mask;
S3D(scene_get_session_mask(scene, &mask));
if(mask) S3D(scene_end_session(scene));
S3D(scene_ref_put(scene));
}
- if(shape) S3D(shape_ref_put(shape));
return res;
error:
@@ -328,18 +166,17 @@ int main(int argc, char *argv[]) {
/* check command arguments */
if(argc < 2 || argc > 4) {
if(argc < 2)
- logger_print(LOGGER_DEFAULT, LOG_ERROR, "Argument required\n");
+ printf("Argument required\n");
else
- logger_print(LOGGER_DEFAULT, LOG_ERROR, "Too many arguments\n");
- logger_print(LOGGER_DEFAULT, LOG_OUTPUT,
- "usage: %s OBJ_FILE [SAMPLES_COUNT [K_SCATTERING]]\n", argv[0]);
+ printf("Too many arguments\n");
+ printf("usage: %s OBJ_FILE [SAMPLES_COUNT [K_SCATTERING]]\n", argv[0]);
return RES_BAD_ARG;
}
/* import file's content in the scene */
res = import_obj(argv[1], &scene);
if(res != RES_OK) {
- logger_print(LOGGER_DEFAULT, LOG_ERROR, "Invalid file '%s'\n", argv[1]);
+ printf("Invalid file '%s'\n", argv[1]);
return res;
}
@@ -350,7 +187,7 @@ int main(int argc, char *argv[]) {
else {
long ns = atol(argv[2]);
if(ns <= 0) {
- logger_print(LOGGER_DEFAULT, LOG_ERROR, "Invalid number of steps `%s'\n", argv[2]);
+ printf("Invalid number of steps `%s'\n", argv[2]);
return RES_BAD_ARG;
}
nsteps = (size_t)ns;
@@ -363,13 +200,13 @@ int main(int argc, char *argv[]) {
else {
ks = atof(argv[3]);
if(ks <= 0) {
- logger_print(LOGGER_DEFAULT, LOG_ERROR,
- "Invalid k-scattering value `%s'\n", argv[3]);
+ printf("Invalid k-scattering value `%s'\n", argv[3]);
return RES_BAD_ARG;
}
}
- compute_4v_s(scene, nsteps, ks);
+ /* compute */
+ CHECK(compute_4v_s(scene, nsteps, ks), RES_OK);
/* check memory leaks */
CHECK(mem_allocated_size(), 0);
