commit c97cb074257507fb1fc001790b6dd9cee31d7270
parent 8e8d90d630d74ec8aab18acaf0c50ca7b0395b5e
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Mon, 7 Mar 2016 11:12:36 +0100
Make parallel the computation of the phase_functions
The per wavelength phase function can be computed in parallel.
Diffstat:
2 files changed, 87 insertions(+), 43 deletions(-)
diff --git a/src/sschiff_estimator.c b/src/sschiff_estimator.c
@@ -105,6 +105,12 @@ struct integrator_context {
struct sschiff_estimator* estimator;
};
+/* Store per thread data of the post integration process */
+struct solver_context {
+ gsl_multiroot_fdfsolver* solver;
+ gsl_vector* init_val;
+};
+
struct sschiff_estimator {
double* wavelengths; /* List of wavelengths in vacuum */
double* angles; /* List of scattering angles */
@@ -285,12 +291,8 @@ compute_path_length
S3D(scene_trace_ray(scn, ray_org, ray_dir, range, &hit));
} while(hit.distance < range[0] || self_hit(&hit_from, &hit, 1.e-6f));
- if(S3D_HIT_NONE(&hit)) {
-/* log_error(dev,
-"Error in computing the radiative path length. This may be due to numerical\n"
-"imprecisions or to a geometry that does not define a closed volume.\n");*/
+ if(S3D_HIT_NONE(&hit))
return RES_BAD_OP;
- }
len += hit.distance - dst;
dst = range[0] = hit.distance;
@@ -697,6 +699,7 @@ function_fdf(const gsl_vector* vec, void* params, gsl_vector* res, gsl_matrix* J
static res_T
solve_wide_angles_phase_function_model_parameters
(struct sschiff_device* dev,
+ struct solver_context* ctx,
const double angle,
const double sin_half_angle,
const double cos_angle,
@@ -707,9 +710,7 @@ solve_wide_angles_phase_function_model_parameters
double* n,
double* r)
{
- gsl_multiroot_fdfsolver* solver = NULL;
gsl_multiroot_function_fdf func;
- gsl_vector* init_val = NULL;
gsl_vector* root = NULL;
struct function_arg arg;
int status;
@@ -721,25 +722,7 @@ solve_wide_angles_phase_function_model_parameters
ASSERT(eq_eps(cos(angle), cos_angle, 1.e-6));
ASSERT(eq_eps(cos(2.0*angle), cos_2_angle, 1.e-6));
- /* Create a Newton-Raphson nonlinear solver without derivatives */
- solver = gsl_multiroot_fdfsolver_alloc(gsl_multiroot_fdfsolver_newton, 1);
- if(!solver) {
- log_error(dev,
- "Not enough memory. Couldn't allocate the GSL Newton solver.\n");
- res = RES_MEM_ERR;
- goto error;
- }
-
- /* Allocate the result vector. */
- init_val = gsl_vector_alloc(1);
- if(!init_val) {
- log_error(dev,
- "Not enough memory. Couldn't allocate the GSL init vector.\n");
- res = RES_MEM_ERR;
- goto error;
- }
-
- /* Fill system arguments */
+ /* Fill system arguments */
arg.angle = angle;
arg.scattering_cross_section = scattering_cross_section;
arg.differential_cross_section = differential_cross_section;
@@ -762,23 +745,21 @@ solve_wide_angles_phase_function_model_parameters
/* Initialise the solver. TODO perform this in an initialisation step of a
* phase_function context */
- gsl_vector_set(init_val, 0, 3);
- gsl_multiroot_fdfsolver_set(solver, &func, init_val);
+ gsl_multiroot_fdfsolver_set(ctx->solver, &func, ctx->init_val);
gsl_set_error_handler_off();
/* Launch the Newton estimation. Iterate up to `max_iter'or until the
* estimation is "good enough". */
iter = 0;
-
do {
- status = gsl_multiroot_fdfsolver_iterate(solver);
+ status = gsl_multiroot_fdfsolver_iterate(ctx->solver);
if(status == GSL_SUCCESS) {
- status = gsl_multiroot_test_residual(solver->f, 1.e-6);
+ status = gsl_multiroot_test_residual(ctx->solver->f, 1.e-6);
}
} while(status == GSL_CONTINUE && ++iter < max_iter);
/* Retrieve the estimated "n" value */
- root = gsl_multiroot_fdfsolver_root(solver);
+ root = gsl_multiroot_fdfsolver_root(ctx->solver);
*n = gsl_vector_get(root, 0);
if(status != GSL_SUCCESS) {
@@ -804,8 +785,6 @@ solve_wide_angles_phase_function_model_parameters
/ (1 + arg.sqr_cos_angle);
exit:
- if(solver) gsl_multiroot_fdfsolver_free(solver);
- if(init_val) gsl_vector_free(init_val);
return res;
error:
goto exit;
@@ -829,6 +808,7 @@ compute_differential_cross_section_cumulative_term
static res_T
compute_phase_function
(struct sschiff_device* dev,
+ struct solver_context* ctx,
const size_t iwlen,
const size_t nangles,
struct sschiff_estimator* estimator)
@@ -847,7 +827,7 @@ compute_phase_function
size_t ilimit_angle; /* Index of the limit angle of the current wavelength */
size_t iangle;
res_T res = RES_OK;
- ASSERT(dev && estimator);
+ ASSERT(dev && estimator && iwlen < sa_size(estimator->wavelengths) && ctx);
/* Fetch the limit angle and precompute some values */
ilimit_angle = estimator->limit_angles[iwlen]-1;
@@ -856,7 +836,6 @@ compute_phase_function
cos_limit_angle = cos(limit_angle);
cos_2_limit_angle = cos(2.0*limit_angle);
- /* Find the connector values between small and wide angles */
get_mc_value
(&estimator->accums[iwlen].scattering_cross_section,
estimator->nrealisations,
@@ -870,8 +849,10 @@ compute_phase_function
estimator->nrealisations,
&limit_differential_cross_section_cumulative);
+ /* Find the connector values between small and wide angles */
res = solve_wide_angles_phase_function_model_parameters
(dev,
+ ctx,
limit_angle,
sin_half_limit_angle,
cos_limit_angle,
@@ -1143,6 +1124,46 @@ error:
goto exit;
}
+static void
+solver_context_release(struct solver_context* ctx)
+{
+ ASSERT(ctx);
+ if(ctx->solver) gsl_multiroot_fdfsolver_free(ctx->solver);
+ if(ctx->init_val) gsl_vector_free(ctx->init_val);
+}
+
+static res_T
+solver_context_init(struct sschiff_device* dev, struct solver_context* ctx)
+{
+ res_T res = RES_OK;
+ ASSERT(ctx);
+ memset(ctx, 0, sizeof(*ctx));
+
+ /* Create a Newton-Raphson nonlinear solver with derivatives */
+ ctx->solver = gsl_multiroot_fdfsolver_alloc(gsl_multiroot_fdfsolver_newton, 1);
+ if(!ctx->solver) {
+ log_error(dev,
+ "Not enough memory. Couldn't allocate the GSL Newton solver.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+ /* Allocate the initial value vector */
+ ctx->init_val = gsl_vector_alloc(1);
+ if(!ctx->init_val) {
+ log_error(dev,
+ "Not enough memory. Couldn't allocate the GSL init vector.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ solver_context_release(ctx);
+ goto exit;
+}
+
static res_T
begin_realisation
(struct sschiff_device* dev,
@@ -1418,11 +1439,12 @@ sschiff_integrate
struct sschiff_estimator** out_estimator)
{
struct integrator_context* ctxs = NULL;
+ struct solver_context* solver_ctxs = NULL;
struct sschiff_estimator* estimator = NULL;
struct ssp_rng** rngs = NULL;
struct ssp_rng_proxy* rng_proxy = NULL;
size_t i;
- size_t iwlen;
+ int iwlen;
int igeom;
ATOMIC res = (ATOMIC)RES_OK;
@@ -1452,12 +1474,18 @@ sschiff_integrate
res = RES_MEM_ERR;
goto error;
}
+ if(!sa_add(solver_ctxs, dev->nthreads)) {
+ log_error(dev, "Couldn'nt allocate the solver contexts.\n");
+ res = RES_MEM_ERR;
+ goto error;
+ }
if(!sa_add(rngs, dev->nthreads)) {
log_error(dev, "Couldn't allocate the list of RNGs.\n");
res = RES_MEM_ERR;
goto error;
}
memset(ctxs, 0, sizeof(ctxs[0])*dev->nthreads);
+ memset(solver_ctxs, 0, sizeof(solver_ctxs[0])*dev->nthreads);
memset(rngs, 0, sizeof(rngs[0])*dev->nthreads);
/* Init the per thread data structures */
@@ -1476,6 +1504,13 @@ sschiff_integrate
(unsigned long)i);
goto error;
}
+ res = solver_context_init(dev, solver_ctxs + i);
+ if(res != RES_OK) {
+ log_error(dev,
+ "Couldn't initialise the solver context of the thread \"%lu\".\n",
+ (unsigned long)i);
+ goto error;
+ }
}
/* Paralell Schiff integration */
@@ -1503,6 +1538,7 @@ sschiff_integrate
/* Merge the per thread integration ressults */
FOR_EACH(i, 0, dev->nthreads) {
+ size_t iwlen;
FOR_EACH(iwlen, 0, nwavelengths) {
size_t iangle;
#define MC_ACCUM(Data) { \
@@ -1526,11 +1562,15 @@ sschiff_integrate
}
}
- /* TODO use parallelism */
- FOR_EACH(i, 0, nwavelengths) {
- if(estimator->limit_angles[i] == INVALID_LIMIT_ANGLE) continue;
+ /* Static scheduling is not necessary to ensure the reproductability;
+ * computations are purely analytics */
+ #pragma omp parallel for
+ for(iwlen=0; iwlen < (int)nwavelengths; ++iwlen) {
+ const int ithread = omp_get_thread_num();
+ if(estimator->limit_angles[iwlen] == INVALID_LIMIT_ANGLE) continue;
/* Do not handle phase function errors */
- compute_phase_function(dev, i, nangles, estimator);
+ compute_phase_function
+ (dev, &solver_ctxs[ithread], (size_t)iwlen, nangles, estimator);
}
exit:
@@ -1539,6 +1579,10 @@ exit:
FOR_EACH(i, 0, dev->nthreads) integrator_context_release(ctxs+i);
sa_release(ctxs);
}
+ if(solver_ctxs) {
+ FOR_EACH(i, 0, dev->nthreads) solver_context_release(solver_ctxs+i);
+ sa_release(solver_ctxs);
+ }
if(rngs) {
FOR_EACH(i, 0, dev->nthreads) if(rngs[i]) SSP(rng_ref_put(rngs[i]));
sa_release(rngs);
diff --git a/src/test_sschiff_estimator_rhodo.c b/src/test_sschiff_estimator_rhodo.c
@@ -602,8 +602,8 @@ main(int argc, char** argv)
const double* angles = NULL;
const double wavelength = 0.4; /* In micron */
const size_t nscatt_angles = 1000;
- const size_t ngeoms = 10000000;
- const size_t ndirs = 10;
+ const size_t ngeoms = 10000;
+ const size_t ndirs = 100;
size_t count;
size_t i;
(void)argc, (void)argv;
