stardis-solver

commit 1b96de6cc62acbd0fdd845c74f1e3a8c7f6b0c01
parent 05511a3d31077d3b127da1dcf2295415ded9e41d
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Tue, 13 Feb 2018 10:55:15 +0100

Add multithreaded parallelism

Diffstat:
M
src/sdis_solve_probe.c
 | 
71
+++++++++++++++++++++++++++++++++++++++++++++++++++++------------------

1 file changed, 53 insertions(+), 18 deletions(-)
diff --git a/src/sdis_solve_probe.c b/src/sdis_solve_probe.c
@@ -27,6 +27,7 @@
 #include "sdis_solve_probe_Xd.h"
 
 #include <star/ssp.h>
+#include <omp.h>
 
 res_T
 sdis_solve_probe
@@ -39,11 +40,14 @@ sdis_solve_probe
 {
   const struct sdis_medium* medium = NULL;
   struct sdis_estimator* estimator = NULL;
-  struct ssp_rng* rng = NULL;
+  struct ssp_rng_proxy* rng_proxy = NULL;
+  struct ssp_rng** rngs = NULL;
   double weight = 0;
   double sqr_weight = 0;
   size_t irealisation = 0;
-  res_T res = RES_OK;
+  size_t N = 0; /* #realisations that do not fail */
+  size_t i;
+  ATOMIC res = RES_OK;
 
   if(!scn || !nrealisations || !position || time < 0 || fp_to_meter <= 0
   || !out_estimator) {
@@ -51,48 +55,79 @@ sdis_solve_probe
     goto error;
   }
 
-  res = scene_get_medium(scn, position, &medium);
+  /* Create the proxy RNG */
+  res = ssp_rng_proxy_create(scn->dev->allocator, &ssp_rng_mt19937_64,
+    scn->dev->nthreads, &rng_proxy);
   if(res != RES_OK) goto error;
 
-  res = ssp_rng_create(scn->dev->allocator, &ssp_rng_mt19937_64, &rng);
-  if(res != RES_OK) goto error;
+  /* Create the per thread RNG */
+  rngs = MEM_CALLOC
+    (scn->dev->allocator, scn->dev->nthreads, sizeof(struct ssp_rng*));
+  if(!rngs) {
+    res = RES_MEM_ERR;
+    goto error;
+  }
+  FOR_EACH(i, 0, scn->dev->nthreads) {
+    res = ssp_rng_proxy_create_rng(rng_proxy, i, rngs+i);
+    if(res != RES_OK) goto error;
+  }
+
+  /* Create the estimator */
   res = estimator_create(scn->dev, &estimator);
   if(res != RES_OK) goto error;
 
-  FOR_EACH(irealisation, 0, nrealisations) {
+  /* Retrieve the medium in which the submitted position lies */
+  res = scene_get_medium(scn, position, &medium);
+  if(res != RES_OK) goto error;
+
+  /* Here we go! Launch the Monte Carlo estimation */
+  #pragma omp parallel for schedule(static) reduction(+:weight,sqr_weight,N)
+  for(irealisation = 0; irealisation < nrealisations; ++irealisation) {
+    res_T res_local;
     double w;
+    const int ithread = omp_get_thread_num();
+    struct ssp_rng* rng = rngs[ithread];
+
+    if(ATOMIC_GET(&res) != RES_OK) continue; /* An error occured */
 
     if(scene_is_2d(scn)) {
-      res = probe_realisation_2d
+      res_local = probe_realisation_2d
         (scn, rng, medium, position, time, fp_to_meter, &w);
     } else {
-      res = probe_realisation_3d
+      res_local = probe_realisation_3d
         (scn, rng, medium, position, time, fp_to_meter, &w);
     }
-    if(res != RES_OK) {
-      if(res == RES_BAD_OP) {
+    if(res_local != RES_OK) {
+      if(res_local == RES_BAD_OP) {
         ++estimator->nfailures;
       } else {
-        goto error;
+        ATOMIC_SET(&res, res_local);
+        continue;
       }
     } else {
       weight += w;
       sqr_weight += w*w;
-      ++estimator->nrealisations;
+      ++N;
     }
   }
 
-  estimator->temperature.E = weight / (double)estimator->nrealisations;
+  estimator->nrealisations = N;
+  estimator->temperature.E = weight / (double)N;
   estimator->temperature.V =
-    sqr_weight / (double)estimator->nrealisations
+    sqr_weight / (double)N
   - estimator->temperature.E * estimator->temperature.E;
-  estimator->temperature.SE =
-    sqrt(estimator->temperature.V / (double)estimator->nrealisations);
+  estimator->temperature.SE = sqrt(estimator->temperature.V / (double)N);
 
 exit:
-  if(rng) SSP(rng_ref_put(rng));
+  if(rngs) {
+    FOR_EACH(i, 0, scn->dev->nthreads)  {
+      if(rngs[i]) SSP(rng_ref_put(rngs[i]));
+    }
+    MEM_RM(scn->dev->allocator, rngs);
+  }
+  if(rng_proxy) SSP(rng_proxy_ref_put(rng_proxy));
   if(out_estimator) *out_estimator = estimator;
-  return res;
+  return (res_T)res;
 error:
   if(estimator) {
     SDIS(estimator_ref_put(estimator));