Imported Upstream version 3.0+20111124upstream/3.0+20111124upstream

1 files changed, 352 insertions, 0 deletions

diff --git a/src/common/DoSample.c b/src/common/DoSample.c
new file mode 100644
index 0000000..a637674
--- /dev/null
+++ b/src/common/DoSample.c
@@ -0,0 +1,352 @@
+/*
+	DoSample.c
+		the actual sampling routine, serial and parallel,
+		for the C versions of the Cuba routines
+		by Thomas Hahn
+		last modified 24 Nov 11 th
+*/
+
+#define MINSLICE 10
+#define MINCORES 1
+//#define MINCORES 2
+
+#if defined(VEGAS) || defined(SUAVE)
+#define VEG_ONLY(...) __VA_ARGS__
+#else
+#define VEG_ONLY(...)
+#endif
+
+#ifdef DIVONNE
+#define DIV_ONLY(...) __VA_ARGS__
+#define LDX(ldx) ldx
+#else
+#define DIV_ONLY(...)
+#define LDX(ldx) t->ndim
+#endif
+
+typedef struct {
+  real *f;
+  number n;
+  VEG_ONLY(count iter;)
+  DIV_ONLY(number neval_opt, neval_cut;
+           count ldx, phase, iregion;)
+#define NREGIONS ldx
+#define NEVAL n
+#define RETVAL phase
+} Slice;
+
+/*********************************************************************/
+
+#ifndef MSG_WAITALL
+/* Windows */
+#define MSG_WAITALL 0
+#endif
+
+static inline int readsock(int fd, void *data, size_t n)
+{
+  ssize_t got;
+  size_t remain = n;
+  do got = recv(fd, data, remain, MSG_WAITALL);
+  while( got > 0 && (data += got, remain -= got) > 0 );
+  return got;
+}
+
+static inline int writesock(int fd, const void *data, size_t n)
+{
+  ssize_t got;
+  size_t remain = n;
+  do got = send(fd, data, remain, MSG_WAITALL);
+  while( got > 0 && (data += got, remain -= got) > 0 );
+  return got;
+}
+
+/*********************************************************************/
+
+static inline bool SampleSerial(cThis *t, number n, creal *x, real *f
+  VEG_ONLY(, creal *w, ccount iter)
+  DIV_ONLY(, ccount ldx))
+{
+  while( n-- ) {
+    if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata
+          VEG_ONLY(, w++, &iter)
+          DIV_ONLY(, &t->phase)) == ABORT ) return true;
+    x += LDX(ldx);
+    f += t->ncomp;
+  }
+  return false;
+}
+
+/*********************************************************************/
+
+static void DoSample(This *t, number n, creal *x, real *f
+  VEG_ONLY(, creal *w, ccount iter)
+  DIV_ONLY(, ccount ldx))
+{
+  char s[128];
+  Slice slice;
+  int ncores;
+
+  t->neval += n;
+
+  ncores = IMin(t->ncores, n/MINSLICE);
+
+  if( ncores < MINCORES ) {
+    if( VERBOSE > 2 ) {
+      sprintf(s, "sampling " NUMBER " points serially", n);
+      Print(s);
+    }
+
+    if( SampleSerial(t, n, x, f
+          VEG_ONLY(, w, iter)
+          DIV_ONLY(, ldx)) ) longjmp(t->abort, -99);
+  }
+  else {
+    int core, abort;
+
+    slice.n = (n + ncores - 1)/ncores;
+
+    if( VERBOSE > 2 ) {
+      sprintf(s, "sampling " NUMBER " points each on %d cores",
+        slice.n, ncores);
+      Print(s);
+    }
+
+    slice.f = f;
+    VEG_ONLY(slice.iter = iter;)
+    DIV_ONLY(slice.ldx = ldx;)
+    DIV_ONLY(slice.phase = t->phase;)
+
+    for( core = 0; core < ncores; ++core ) {
+      cint fd = t->child[core];
+      writesock(fd, &slice, sizeof slice);
+      VEG_ONLY(writesock(fd, w, slice.n*sizeof *w);)
+      writesock(fd, x, slice.n*LDX(ldx)*sizeof *x);
+
+      VEG_ONLY(w += n;)
+      x += slice.n*LDX(ldx);
+      slice.f += slice.n*t->ncomp;
+      n -= slice.n;
+      slice.n = IMin(slice.n, n);
+    }
+
+    abort = 0;
+    for( core = ncores; --core >= 0; ) {
+      cint fd = t->child[core];
+      readsock(fd, &slice, sizeof slice);
+      if( slice.n == 0 ) abort = 1;
+      else readsock(fd, slice.f, slice.n*t->ncomp*sizeof *f);
+    }
+    if( abort ) longjmp(t->abort, -99);
+  }
+}
+
+/*********************************************************************/
+
+#ifdef DIVONNE
+static inline int ReadyCore(cThis *t)
+{
+  int core;
+  fd_set ready;
+
+  memcpy(&ready, &t->children, sizeof ready);
+  select(t->nchildren, &ready, NULL, NULL, NULL);
+
+  for( core = 0; core < t->ncores; ++core )
+    if( FD_ISSET(t->child[core], &ready) ) break;
+
+  return core;
+}
+
+/*********************************************************************/
+
+static int ExploreParent(This *t, cint iregion)
+{
+  TYPEDEFREGION;
+  Region *region;
+  Slice slice;
+  int ireg = iregion, core = t->running;
+
+  if( t->ncores < MINCORES ) return Explore(t, iregion);
+
+  if( t->running >= ((iregion < 0) ? 1 : t->ncores) ) {
+    Totals totals[t->ncomp];
+    count comp, succ;
+    cint fd = t->child[core = ReadyCore(t)];
+
+    --t->running;
+    readsock(fd, &slice, sizeof slice);
+//DEBSLICE("parent read", fd, slice);
+    ireg = slice.iregion;
+    region = RegionPtr(ireg);
+    succ = ireg + region->next;
+    readsock(fd, region, sizeof(Region));
+    if( --slice.NREGIONS > 0 ) {
+      region->next = t->nregions - ireg;
+      EnlargeRegions(t, slice.NREGIONS);
+      readsock(fd, RegionPtr(t->nregions), slice.NREGIONS*sizeof(Region));
+      t->nregions += slice.NREGIONS;
+      RegionPtr(t->nregions-1)->next = succ - t->nregions + 1;
+    }
+
+    readsock(fd, totals, sizeof totals);
+    for( comp = 0; comp < t->ncomp; ++comp )
+      t->totals[comp].secondspread =
+        Max(t->totals[comp].secondspread, totals[comp].secondspread);
+
+    t->neval += slice.NEVAL;
+    t->neval_opt += slice.neval_opt;
+    t->neval_cut += slice.neval_cut;
+
+    if( slice.RETVAL == -1 ) return -1;
+  }
+
+  if( iregion >= 0 ) {
+    region = RegionPtr(iregion);
+    cint fd = t->child[core];
+    slice.n = 0;
+    slice.phase = t->phase;
+    slice.iregion = iregion;
+//DEBSLICE("  parent write", fd, slice);
+    writesock(fd, &slice, sizeof slice);
+    writesock(fd, &t->samples[region->isamples], sizeof(Samples));
+    writesock(fd, region, sizeof *region);
+    writesock(fd, t->totals, sizeof *t->totals);
+    region->depth = 0;
+    ++t->running;
+  }
+
+  return ireg;
+}
+#endif
+
+/*********************************************************************/
+
+static inline void DoChild(This *t, cint fd)
+{
+  Slice slice;
+
+#ifdef DIVONNE
+  TYPEDEFREGION;
+  Totals totals[t->ncomp];
+
+  t->totals = totals;
+  t->ncores = 0;	/* no recursive forks */
+  AllocRegions(t);
+  SamplesIni(&t->samples[0]);
+  t->samples[0].n = 0;
+  SamplesIni(&t->samples[1]);
+  t->samples[1].n = 0;
+  SamplesIni(&t->samples[2]);
+  t->samples[2].n = 0;
+#endif
+
+  while( readsock(fd, &slice, sizeof slice) ) {
+    number n = slice.n;
+    DIV_ONLY(t->phase = slice.phase;)
+//DEBSLICE("  child read", fd, slice);
+    if( n > 0 ) {
+      VEG_ONLY(real w[n];)
+      real x[n*LDX(slice.ldx)];
+      real f[n*t->ncomp];
+
+      VEG_ONLY(readsock(fd, w, sizeof w);)
+      readsock(fd, x, sizeof x);
+
+      if( SampleSerial(t, n, x, f
+            VEG_ONLY(, w, slice.iter)
+            DIV_ONLY(, slice.ldx)) ) slice.n = 0;
+      writesock(fd, &slice, sizeof slice);
+      if( slice.n ) writesock(fd, f, sizeof f);
+    }
+#ifdef DIVONNE
+    else {
+      Samples *samples, psamples;
+
+      readsock(fd, &psamples, sizeof psamples);
+      readsock(fd, RegionPtr(0), sizeof(Region));
+      readsock(fd, totals, sizeof totals);
+      t->nregions = 1;
+      t->neval = t->neval_opt = t->neval_cut = 0;
+
+      samples = &t->samples[RegionPtr(0)->isamples];
+      if( psamples.n != samples->n ) {
+        SamplesFree(samples);
+        *samples = psamples;
+        SamplesAlloc(t, samples);
+      }
+
+      slice.RETVAL = Explore(t, 0);
+      slice.NREGIONS = t->nregions;
+      slice.NEVAL = t->neval;
+      slice.neval_opt = t->neval_opt;
+      slice.neval_cut = t->neval_cut;
+//DEBSLICE("child write", fd, slice);
+      writesock(fd, &slice, sizeof slice);
+      writesock(fd, RegionPtr(0), t->nregions*sizeof(Region));
+      writesock(fd, totals, sizeof totals);
+    }
+#endif
+  }
+
+  exit(0);
+}
+
+/*********************************************************************/
+
+static inline void ForkCores(This *t)
+{
+  int core;
+  cchar *env = getenv("CUBACORES");
+
+  t->ncores = env ? atoi(env) : sysconf(_SC_NPROCESSORS_ONLN);
+#ifdef HAVE_GETLOADAVG
+  if( env == NULL || t->ncores < 0 ) {
+    double load = 0;
+    getloadavg(&load, 1);
+    t->ncores = abs(t->ncores) - floor(load);
+  }
+#endif
+
+#ifdef DIVONNE
+  t->nchildren = t->running = 0;
+#endif
+
+  if( t->ncores < MINCORES ) return;
+  if( VERBOSE ) printf("using %d cores\n", t->ncores);
+  fflush(stdout);
+
+  Alloc(t->child, t->ncores);
+  for( core = 0; core < t->ncores; ++core ) {
+    int fd[2];
+    pid_t pid;
+    assert(
+      socketpair(AF_LOCAL, SOCK_STREAM, 0, fd) != -1 &&
+      (pid = fork()) != -1 );
+    if( pid == 0 ) {
+      close(fd[0]);
+      DoChild(t, fd[1]);
+    }
+    close(fd[1]);
+    t->child[core] = fd[0];
+#ifdef DIVONNE
+    FD_SET(fd[0], &t->children);
+    t->nchildren = IMax(t->nchildren, fd[0] + 1);
+#endif
+  }
+}
+
+/*********************************************************************/
+
+static inline void WaitCores(cThis *t)
+{
+  if( t->ncores >= MINCORES ) {
+    int core;
+    pid_t pid;
+    for( core = 0; core < t->ncores; ++core )
+      close(t->child[core]);
+    free(t->child);
+    for( core = 0; core < t->ncores; ++core )
+      wait(&pid);
+  }
+}
+