upgrade the Cuba library from 3.2 -> 4.2

configure.ac

@@ -1,7 +1,7 @@
 AC_REVISION([m4_esyscmd_s([git describe --always])])
 
 AC_PREREQ(2.63)
-AC_INIT([musrfit],[1.2.0],[andreas.suter@psi.ch])
+AC_INIT([musrfit],[1.2.1],[andreas.suter@psi.ch])
 AC_CONFIG_AUX_DIR(admin)
 AC_CANONICAL_HOST
 #AC_MSG_RESULT([${host} ${host_cpu} ${host_vendor} ${host_os}])
@@ -36,7 +36,7 @@ dnl -----------------------------------------------
 #release versioning
 MUSR_MAJOR_VERSION=1
 MUSR_MINOR_VERSION=2
-MUSR_MICRO_VERSION=0
+MUSR_MICRO_VERSION=1
 
 #release versioning
 MUSR_ROOT_MAJOR_VERSION=1
@@ -69,7 +69,7 @@ PLUGIN_MINOR_VERSION=0
 PLUGIN_MICRO_VERSION=0
 
 #release versioning
-CUBA_MAJOR_VERSION=3
+CUBA_MAJOR_VERSION=4
 CUBA_MINOR_VERSION=2
 CUBA_MICRO_VERSION=0
 

src/external/BMWtools/BMWIntegrator.cpp

@@ -31,6 +31,7 @@
 #include "cuba.h"
 
 #define USERDATA NULL
+#define SPIN NULL
 #define SEED 0
 #define STATEFILE NULL
 
@@ -45,6 +46,7 @@ std::vector<double> TDWaveGapIntegralCuhre::fPar;
 double TDWaveGapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-4);
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
@@ -57,9 +59,9 @@ double TDWaveGapIntegralCuhre::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Cuhre(fNDim, NCOMP, Integrand, USERDATA,
+  Cuhre(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
-    KEY, STATEFILE,
+    KEY, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];
@@ -96,6 +98,7 @@ std::vector<double> TCosSqDWaveGapIntegralCuhre::fPar;
 double TCosSqDWaveGapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-8);
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
@@ -108,9 +111,9 @@ double TCosSqDWaveGapIntegralCuhre::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Cuhre(fNDim, NCOMP, Integrand, USERDATA,
+  Cuhre(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
-    KEY, STATEFILE,
+    KEY, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];
@@ -147,6 +150,7 @@ std::vector<double> TSinSqDWaveGapIntegralCuhre::fPar;
 double TSinSqDWaveGapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-8);
   const double EPSABS (1e-10);
   const unsigned int VERBOSE (0);
@@ -159,9 +163,9 @@ double TSinSqDWaveGapIntegralCuhre::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Cuhre(fNDim, NCOMP, Integrand, USERDATA,
+  Cuhre(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
-    KEY, STATEFILE,
+    KEY, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];
@@ -198,6 +202,7 @@ std::vector<double> TAnSWaveGapIntegralCuhre::fPar;
 double TAnSWaveGapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-4);
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
@@ -210,9 +215,9 @@ double TAnSWaveGapIntegralCuhre::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Cuhre(fNDim, NCOMP, Integrand, USERDATA,
+  Cuhre(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
-    KEY, STATEFILE,
+    KEY, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];
@@ -249,6 +254,7 @@ std::vector<double> TAnSWaveGapIntegralDivonne::fPar;
 double TAnSWaveGapIntegralDivonne::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-4);
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
@@ -268,10 +274,10 @@ double TAnSWaveGapIntegralDivonne::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Divonne(fNDim, NCOMP, Integrand, USERDATA,
+  Divonne(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE, SEED, MINEVAL, MAXEVAL,
     KEY1, KEY2, KEY3, MAXPASS, BORDER, MAXCHISQ, MINDEVIATION,
-    NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL, STATEFILE,
+    NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];
@@ -308,6 +314,7 @@ std::vector<double> TAnSWaveGapIntegralSuave::fPar;
 double TAnSWaveGapIntegralSuave::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-4);
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
@@ -316,14 +323,15 @@ double TAnSWaveGapIntegralSuave::IntegrateFunc()
   const unsigned int MAXEVAL (1000000);
 
   const unsigned int NNEW (1000);
+  const unsigned int NMIN (2);
   const double FLATNESS (25.);
 
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Suave(fNDim, NCOMP, Integrand, USERDATA,
+  Suave(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE | LAST, SEED, MINEVAL, MAXEVAL,
-    NNEW, FLATNESS, STATEFILE,
+    NNEW, NMIN, FLATNESS, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];
@@ -360,6 +368,7 @@ std::vector<double> TNonMonDWave1GapIntegralCuhre::fPar;
 double TNonMonDWave1GapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-4);
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
@@ -372,9 +381,9 @@ double TNonMonDWave1GapIntegralCuhre::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Cuhre(fNDim, NCOMP, Integrand, USERDATA,
+  Cuhre(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
-    KEY, STATEFILE,
+    KEY, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];
@@ -411,6 +420,7 @@ std::vector<double> TNonMonDWave2GapIntegralCuhre::fPar;
 double TNonMonDWave2GapIntegralCuhre::IntegrateFunc()
 {
   const unsigned int NCOMP(1);
+  const unsigned int NVEC(1);
   const double EPSREL (1e-4);
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
@@ -423,9 +433,9 @@ double TNonMonDWave2GapIntegralCuhre::IntegrateFunc()
   int nregions, neval, fail;
   double integral[NCOMP], error[NCOMP], prob[NCOMP];
 
-  Cuhre(fNDim, NCOMP, Integrand, USERDATA,
+  Cuhre(fNDim, NCOMP, Integrand, USERDATA, NVEC,
     EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
-    KEY, STATEFILE,
+    KEY, STATEFILE, SPIN,
     &nregions, &neval, &fail, integral, error, prob);
 
   return integral[0];

src/external/libCuba/src/Makefile.am

@@ -11,7 +11,7 @@ lib_LTLIBRARIES = libcuba.la
 
 libcuba_la_SOURCES =
 
-libcuba_la_LIBADD = cuhre/libcuhre.la divonne/libdivonne.la suave/libsuave.la vegas/libvegas.la -lm
+libcuba_la_LIBADD = common/libcommon.la cuhre/libcuhre.la divonne/libdivonne.la suave/libsuave.la vegas/libvegas.la -lm
 libcuba_la_LDFLAGS = -version-info $(CUBA_LIBRARY_VERSION) -release $(CUBA_RELEASE) $(AM_LDFLAGS)
 
 pkgconfigdir = $(libdir)/pkgconfig

src/external/libCuba/src/common/CSample.c

@@ -3,19 +3,27 @@
 		the serial sampling routine
 		for the C versions of the Cuba routines
 		by Thomas Hahn
-		last modified 19 Dec 11 th
+		last modified 9 Oct 14 th
 */
 
 
-static inline number SampleRaw(cThis *t, number n, creal *x, real *f
+coreinit cubafun_;
-  VES_ONLY(, creal *w, ccount iter))
+extern int cubaverb_;
+extern corespec cubaworkers_;
+
+
+static inline number SampleRaw(This *t, number n, creal *x, real *f,
+  cint core VES_ONLY(, creal *w, ccount iter))
 {
-  for( ; n; --n ) {
+  number nvec;
-    if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata
+  for( nvec = t->nvec; n > 0; n -= nvec ) {
-          VES_ONLY(, w++, &iter)
+    nvec = IMin(n, nvec);
+    if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata, &nvec, &core
+          VES_ONLY(, w, &iter)
           DIV_ONLY(, &t->phase)) == ABORT ) return -1;
-    x += t->ndim;
+    VES_ONLY(w += nvec;)
-    f += t->ncomp;
+    x += nvec*t->ndim;
+    f += nvec*t->ncomp;
   }
   return 0;
 }
@@ -25,8 +33,9 @@ static inline number SampleRaw(cThis *t, number n, creal *x, real *f
 static inline void DoSampleSerial(This *t, cnumber n, creal *x, real *f
   VES_ONLY(, creal *w, ccount iter))
 {
+  MasterInit();
   t->neval += n;
-  if( SampleRaw(t, n, x, f VES_ONLY(, w, iter)) ) 
+  if( SampleRaw(t, n, x, f, -1 VES_ONLY(, w, iter)) ) 
     longjmp(t->abort, -99);
 }
 
@@ -43,7 +52,13 @@ DIV_ONLY(static int Explore(This *t, cint iregion);)
 #define DoSample DoSampleSerial
 #define Explore ExploreSerial
 #define ForkCores(t)
-#define WaitCores(t)
+
+static inline void WaitCores(This *t, Spin **pspin)
+{
+  if( Invalid(pspin) ) MasterExit();
+}
+
+#define WaitCores(t, pspin)
 
 #endif
 
@@ -51,7 +66,7 @@ DIV_ONLY(static int Explore(This *t, cint iregion);)
 static inline count SampleExtra(This *t, cBounds *b)
 {
   number n = t->nextra;
-  t->peakfinder(&t->ndim, b, &n, t->xextra);
+  t->peakfinder(&t->ndim, b, &n, t->xextra, t->userdata);
   DoSample(t, n, t->xextra, t->fextra);
   return n;
 }
@@ -60,7 +75,7 @@ static inline count SampleExtra(This *t, cBounds *b)
 #include "common.c"
 
 #ifdef HAVE_FORK
-#include "Fork.c"
+#include "Parallel.c"
 #endif
 
 #include "Integrate.c"

src/external/libCuba/src/common/ChiSquare.c

@@ -3,7 +3,7 @@
 		the chi-square cdf
 		after W.J. Kennedy and J.E. Gentle,
 		Statistical computing, p. 116
-		last modified 9 Feb 05 th
+		last modified 12 Mar 15 th
 */
 
 #ifdef HAVE_ERF
@@ -31,7 +31,7 @@ static real ChiSquare(creal x, cint df)
   if( df > 1000 ) {
     if( x < 2 ) return 0;
     y = 2./(9*df);
-    y = (pow(x/df, 1/3.) - (1 - y))/sqrt(y);
+    y = (powx(x/df, 1/3.) - (1 - y))/sqrtx(y);
     if( y > 5 ) return 1;
     if( y < -18.8055 ) return 0;
     return Normal(y);
@@ -40,13 +40,13 @@ static real ChiSquare(creal x, cint df)
   y = .5*x;
 
   if( df & 1 ) {
-    creal sqrty = sqrt(y);
+    creal sqrty = sqrtx(y);
     real h = Erf(sqrty);
     count i;
 
     if( df == 1 ) return h;
 
-    y = sqrty*exp(-y)/.8862269254527579825931;
+    y = sqrty*expx(-y)/.8862269254527579825931;
     for( i = 3; i < df; i += 2 ) {
       h -= y;
       y *= x/i;
@@ -54,7 +54,7 @@ static real ChiSquare(creal x, cint df)
     y = h - y;
   }
   else {
-    real term = exp(-y), sum = term;
+    real term = expx(-y), sum = term;
     count i;
 
     for( i = 1; i < df/2; ++i )

src/external/libCuba/src/common/Erf.c

@@ -4,7 +4,7 @@
 		= 2/Sqrt[Pi] Integrate[Exp[-t^2], {t, 0, x}]
 		Code from Takuya Ooura's gamerf2a.f
 		http://www.kurims.kyoto-u.ac.jp/~ooura/gamerf.html
-		last modified 8 Feb 05 th
+		last modified 12 Mar 15 th
 */
 
 
@@ -21,12 +21,12 @@ static real Erfc(creal x)
     1.66642447174307753e-07, 1.48455557345597957e+01,
     6.10399733098688199e+00, 1.26974899965115684e+01 };
   real y = x*x;
-  y = exp(-y)*x*(
+  y = expx(-y)*x*(
     c[0]/(y + c[1]) + c[2]/(y + c[3]) +
     c[4]/(y + c[5]) + c[6]/(y + c[7]) +
     c[8]/(y + c[9]) + c[10]/(y + c[11]) +
     c[12]/(y + c[13]) + c[14]/(y + c[15]) );
-  if( x < c[16] ) y += 2/(exp(c[17]*x) + 1);
+  if( x < c[16] ) y += 2/(expx(c[17]*x) + 1);
   return y;
 }
 
@@ -40,7 +40,7 @@ static real Erf(creal x)
    -2.68661698447642378e-02,
     5.22387877685618101e-03,
    -8.49202435186918470e-04 };
-  real y = fabs(x);
+  real y = fabsx(x);
   if( y > .125 ) {
     y = 1 - Erfc(y);
     return (x > 0) ? y : -y;

src/external/libCuba/src/common/Fork.c

@@ -1,391 +1,98 @@
 /*
 	Fork.c
-		the parallel sampling routine
+		fork the cores for parallel sampling
-		for the C versions of the Cuba routines
+		(C version only)
 		by Thomas Hahn
-		last modified 25 Sep 13 th
+		last modified 23 Apr 15 th
 */
 
-#include <sys/select.h>
 
-#define MINSLICE 10
+#define ROUTINE "cubafork"
+#include "stddecl.h"
+
+#ifdef HAVE_FORK
+
+#include "sock.h"
+
 #define MINCORES 1
-/*#define MINCORES 2*/
 
-typedef struct {
+coreinit cubafun_;
-  number n, m, i;
+extern int cubaverb_;
-  VES_ONLY(count iter;)
+extern corespec cubaworkers_;
-  DIV_ONLY(int phase SHM_ONLY(, shmid);)
-} Slice;
 
-workerini cubaini;
+/*********************************************************************/
 
-#if defined HAVE_SHMGET && (defined SUAVE || defined DIVONNE)
+static inline void Child(cint fd, cint core)
-#define FRAMECOPY
+{
-#endif
+  dispatch d;
 
-#ifdef DEBUG
+  while( readsock(fd, &d, sizeof d) == sizeof d ) {
-#define TERM_RED "\e[31m"
+    if( d.thissize ) {
-#define TERM_BLUE "\e[34m"
+      MemAlloc(d.thisptr, d.thissize);
-#define TERM_RESET "\e[0m\n"
+      WORKER("reading This (%lu)", d.thissize);
-#define MASTER(s, ...) \
+      readsock(fd, d.thisptr, d.thissize);
-fprintf(stderr, TERM_RED ROUTINE " master %d(%d): " s TERM_RESET, core, getpid(), ##__VA_ARGS__)
+    }
-#define WORKER(s, ...) \
+    WORKER("running %p on fd %d", d.thisptr, fd);
-fprintf(stderr, TERM_BLUE ROUTINE " worker %d(%d): " s TERM_RESET, core, getpid(), ##__VA_ARGS__)
+    d.worker(d.thisptr, d.thissize, core, fd);
+    if( d.thissize ) free(d.thisptr);
+  }
+}
+
+/*********************************************************************/
+
+Extern void SUFFIX(cubafork)(Spin **pspin)
+{
+  char out[128];
+  int cores, core;
+  fdpid *pfp;
+  Spin *spin;
+
+  VerboseInit();
+
+  EnvInit(cubaworkers_.paccel, "CUBAACCELMAX", 1000);
+  EnvInit(cubaworkers_.pcores, "CUBACORESMAX", 10000);
+  EnvInit(cubaworkers_.naccel, "CUBAACCEL", 0);
+  EnvInit(cubaworkers_.ncores, "CUBACORES", -sysconf(_SC_NPROCESSORS_ONLN));
+
+#ifdef HAVE_GETLOADAVG
+  if( cubaworkers_.ncores < 0 ) {
+    static int load = uninitialized;
+    if( load == uninitialized ) {
+      double loadavg;
+      getloadavg(&loadavg, 1);
+      load = floor(loadavg);
+    }
+    cubaworkers_.ncores = IMax(-cubaworkers_.ncores - load, 0);
+  }
 #else
-#define MASTER(s, ...)
+  cubaworkers_.ncores = abs(cubaworkers_.ncores);
-#define WORKER(s, ...)
 #endif
 
-/*********************************************************************/
+  cores = cubaworkers_.naccel + cubaworkers_.ncores;
-
+  if( cores < MINCORES ) {
-#ifndef MSG_WAITALL
+    *pspin = NULL;
-/* Windows */
+    return;
-#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 void DoSample(This *t, number n, creal *x, real *f
-  VES_ONLY(, creal *w, ccount iter))
-{
-  cint ncores = IMin(t->ncores, n/MINSLICE);
-
-  if( ncores < MINCORES ) DoSampleSerial(t, n, x, f VES_ONLY(, w, iter));
-  else {
-    Slice slice;
-    int core, abort;
-    number nx;
-    char s[128];
-
-    t->neval += n;
-
-    nx = n % ncores;
-    slice.m = slice.n = n/ncores + 1;
-    if( VERBOSE > 2 ) {
-      sprintf(s, "sampling " NUMBER " points each on %d cores",
-        slice.n, ncores);
-      Print(s);
-    }
-
-    slice.i = 0;
-    VES_ONLY(slice.iter = iter;)
-    DIV_ONLY(slice.phase = t->phase;)
-
-#ifdef DIVONNE
-    if( n > t->nframe ) {
-      FrameFree(t, ShmRm(t));
-      t->nframe = n;
-      FrameAlloc(t);
-    }
-    SHM_ONLY(slice.shmid = t->shmid;)
-#endif
-
-    SHM_ONLY(if( t->shmid != -1 ) {
-      slice.m = n;
-#ifdef FRAMECOPY
-      VES_ONLY(Copy(t->frame, w, n);)
-      Copy(t->frame + n*NW, x, n*t->ndim);
-#endif
-    })
-
-    for( core = 0; core < ncores; ++core ) {
-      cint fd = t->child[core];
-      slice.n -= (core == nx);
-      MASTER("sending samples (sli:%lu[+" VES_ONLY(NUMBER "w:%lu+")
-        NUMBER "x:%lu]) to fd %d",
-        sizeof slice, VES_ONLY(slice.n, sizeof *w,)
-        slice.n, t->ndim*sizeof *x, fd);
-      writesock(fd, &slice, sizeof slice);
-      SHM_ONLY(if( t->shmid == -1 )) {
-        VES_ONLY(writesock(fd, w, slice.n*sizeof *w);
-                 w += slice.n;)
-        writesock(fd, x, slice.n*t->ndim*sizeof *x);
-        x += slice.n*t->ndim;
-      }
-      slice.i += slice.n;
-    }
-
-    abort = 0;
-    for( core = 0; core < ncores; ++core ) {
-      cint fd = t->child[core];
-      readsock(fd, &slice, sizeof slice);
-      MASTER("reading samples (sli:%lu[+" NUMBER "f:%lu]) from fd %d",
-        sizeof slice, slice.n, t->ncomp*sizeof *f, fd);
-      if( slice.n == -1 ) abort = 1;
-      else SHM_ONLY(if( t->shmid == -1 )) readsock(fd,
-        f + slice.i*t->ncomp, slice.n*t->ncomp*sizeof *f);
-    }
-    if( abort ) longjmp(t->abort, -99);
-
-#ifdef FRAMECOPY
-    if( t->shmid != -1 )
-      Copy(f, t->frame + slice.m*(NW + t->ndim), slice.m*t->ncomp);
-#endif
-  }
-}
-
-/*********************************************************************/
-
-#ifdef DIVONNE
-
-static inline int ReadyCore(cThis *t)
-{
-  fd_set ready;
-  int core, n = 0;
-
-  FD_ZERO(&ready);
-  for( core = 0; core < t->ncores; ++core ) {
-    FD_SET(t->child[core], &ready);
-    n = IMax(n, t->child[core]);
-  }
-  select(n + 1, &ready, NULL, NULL, NULL);
-
-  for( core = 0; core < t->ncores; ++core )
-    if( FD_ISSET(t->child[core], &ready) ) break;
-
-  return core;
-}
-
-/*********************************************************************/
-
-typedef struct {
-  number neval, neval_opt, neval_cut;
-  count nregions, iregion, retval;
-} ExploreResult;
-
-static int Explore(This *t, cint iregion)
-{
-  csize_t regionsize = RegionSize;
-  Region *region;
-  int ireg = iregion, core = t->running;
-  Vector(Totals, totals, NCOMP);
-
-  if( t->ncores < MINCORES ) return ExploreSerial(t, iregion);
-
-  if( t->running >= ((iregion < 0) ? 1 : t->ncores) ) {
-    cint fd = t->child[core = ReadyCore(t)];
-    ExploreResult res;
-    count comp, succ;
-
-    --t->running;
-    MASTER("reading res + region (res:%lu+reg:%lu) from fd %d",
-      sizeof res, regionsize, fd);
-    readsock(fd, &res, sizeof res);
-    ireg = res.iregion;
-    region = RegionPtr(ireg);
-    succ = ireg + region->next;
-    readsock(fd, region, regionsize);
-    if( --res.nregions > 0 ) {
-      region->next = t->nregions - ireg;
-      EnlargeRegions(t, res.nregions);
-      MASTER("reading regions (%dreg:%lu) from fd %d",
-        res.nregions, regionsize, fd);
-      readsock(fd, RegionPtr(t->nregions), res.nregions*regionsize);
-      t->nregions += res.nregions;
-
-      RegionPtr(t->nregions-1)->next = succ - t->nregions + 1;
-    }
-
-    MASTER("reading totals (tot:%lu) from fd %d",
-      t->ncomp*sizeof(Totals), fd);
-    readsock(fd, totals, t->ncomp*sizeof(Totals));
-    for( comp = 0; comp < t->ncomp; ++comp )
-      t->totals[comp].secondspread =
-        Max(t->totals[comp].secondspread, totals[comp].secondspread);
-
-    t->neval += res.neval;
-    t->neval_opt += res.neval_opt;
-    t->neval_cut += res.neval_cut;
-
-    if( res.retval == -1 ) return -1;
   }
 
-  if( iregion >= 0 ) {
+  if( cubaverb_ ) {
-    Slice slice;
+    sprintf(out, "using %d cores %d accelerators via "
-    cint fd = t->child[core];
-    slice.n = 0;
-    slice.i = iregion;
-    slice.phase = t->phase;
-    region = RegionPtr(iregion);
-    MASTER("writing region (sli:%lu+sam:%lu+reg:%lu+tot:%lu) to fd %d",
-      sizeof slice, sizeof(Samples), regionsize,
-      t->ncomp*sizeof(Totals), fd);
-    writesock(fd, &slice, sizeof slice);
-    writesock(fd, &t->samples[region->isamples], sizeof(Samples));
-    writesock(fd, region, regionsize);
-    writesock(fd, t->totals, t->ncomp*sizeof(Totals));
-    region->depth = 0;
-    ++t->running;
-  }
-
-  return ireg;
-}
-#endif
-
-/*********************************************************************/
-
-static void DoChild(This *t, cint core, cint fd)
-{
-  Slice slice;
-
-#ifdef DIVONNE
-  csize_t regionsize = RegionSize;
-  Vector(Totals, totals, NCOMP);
-  ExploreResult res;
-
-  t->totals = totals;
-  t->ncores = 0;	/* no recursive forks */
-  t->size = 2*t->ndim + 2;
-  AllocRegions(t);
-#endif
-
-#ifdef SUAVE
-  SHM_ONLY(if( t->shmid == -1 ))
-    MemAlloc(t->frame, t->nframe*SAMPLESIZE);
-#endif
-
-  if( cubaini.initfun ) cubaini.initfun(cubaini.initarg);
-
-  while( readsock(fd, &slice, sizeof slice) ) {
-    number n = slice.n;
-    DIV_ONLY(t->phase = slice.phase;)
-    if( n > 0 ) {
-      real VES_ONLY(*w,) *x, *f;
-      WORKER("reading samples (sli:%lu[+" VES_ONLY(NUMBER "w:%lu+")
-        NUMBER "x:%lu]) from fd %d",
-        sizeof slice, VES_ONLY(n, sizeof *w,) n, t->ndim*sizeof *x, fd);
-
-#ifdef DIVONNE
-      if( slice.m > t->nframe ) {
-        FrameFree(t);
-        t->nframe = slice.m;
-        SHM_ONLY(t->shmid = slice.shmid; ShmMap(t) else)
-        MemAlloc(t->frame, t->nframe*SAMPLESIZE);
-      }
-#endif
-
-      VES_ONLY(w = t->frame;)
-      x = t->frame + slice.m*NW;
-      f = x + slice.m*t->ndim;
-
-      SHM_ONLY(if( t->shmid != -1 ) {
-        VES_ONLY(w += slice.i;)
-        x += slice.i*t->ndim;
-        f += slice.i*t->ncomp;
-      }
-      else) {
-        VES_ONLY(readsock(fd, w, n*sizeof *w);)
-        readsock(fd, x, n*t->ndim*sizeof *x);
-      }
-
-      slice.n |= SampleRaw(t, n, x, f VES_ONLY(, w, slice.iter));
-      WORKER("writing samples (sli:%lu[+" NUMBER "f:%lu]) to fd %d",
-        sizeof slice, slice.n, t->ncomp*sizeof *f, fd);
-      writesock(fd, &slice, sizeof slice);
-      if( SHM_ONLY(t->shmid == -1 &&) slice.n != -1 )
-        writesock(fd, f, slice.n*t->ncomp*sizeof *f);
-    }
-#ifdef DIVONNE
-    else {
-      Samples *samples, psamples;
-
-      WORKER("reading region (sli:%lu+sam:%lu+reg:%lu+tot:%lu) from fd %d",
-        sizeof slice, sizeof(Samples), regionsize,
-        t->ncomp*sizeof(Totals), fd);
-      readsock(fd, &psamples, sizeof(Samples));
-      readsock(fd, t->region, regionsize);
-      readsock(fd, totals, t->ncomp*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);
-      }
-
-      res.retval = ExploreSerial(t, 0);
-      res.neval = t->neval;
-      res.neval_opt = t->neval_opt;
-      res.neval_cut = t->neval_cut;
-      res.nregions = t->nregions;
-      res.iregion = slice.i;
-      WORKER("writing regions (res:%lu+%dreg:%lu+tot:%lu) to fd %d",
-        sizeof res, t->nregions, regionsize,
-        t->ncomp*sizeof(Totals), fd);
-      writesock(fd, &res, sizeof res);
-      writesock(fd, t->region, t->nregions*regionsize);
-      writesock(fd, totals, t->ncomp*sizeof(Totals));
-    }
-#endif
-  }
-
-  WORKER("wrapping up");
-  if( cubaini.exitfun ) cubaini.exitfun(cubaini.exitarg);
-
-  exit(0);
-}
-
-/*********************************************************************/
-
-#ifdef HAVE_GETLOADAVG
-double cubaloadavg_;
-#endif
-
-static inline void ForkCores(This *t)
-{
-  int core;
-  char s[128];
-  cchar *env = getenv("CUBACORES");
-
-  t->ncores = env ? atoi(env) : sysconf(_SC_NPROCESSORS_ONLN);
-#ifdef HAVE_GETLOADAVG
-  if( env == NULL || t->ncores < 0 ) {
-    if( cubaloadavg_ < 0 ) getloadavg(&cubaloadavg_, 1);
-    t->ncores = abs(t->ncores) - floor(cubaloadavg_);
-  }
-#endif
-
-  DIV_ONLY(t->running = 0;)
-
-  if( t->ncores < MINCORES ) return;
-  if( VERBOSE ) {
-    sprintf(s, "using %d cores via "
 #ifdef HAVE_SHMGET
       "shared memory",
 #else
       "pipes",
 #endif
-      t->ncores);
+      cubaworkers_.ncores, cubaworkers_.naccel);
-    Print(s);
+    Print(out);
   }
 
   fflush(NULL);		/* make sure all buffers are flushed,
 			   or else buffered content will be written
 			   out multiply, at each child's exit(0) */
 
-  Alloc(t->child, t->ncores);
+  MemAlloc(spin, sizeof *spin + cores*sizeof *spin->fp);
-  for( core = 0; core < t->ncores; ++core ) {
+  spin->spec = cubaworkers_;
+  pfp = spin->fp;
+  for( core = -spin->spec.naccel; core < spin->spec.ncores; ++core ) {
     int fd[2];
     pid_t pid;
     assert(
@@ -393,32 +100,65 @@ static inline void ForkCores(This *t)
       (pid = fork()) != -1 );
     if( pid == 0 ) {
       close(fd[0]);
-      DoChild(t, core, fd[1]);
+      free(spin);
+      Child(fd[1], core);
+      exit(0);
     }
     MASTER("forked pid %d pipe %d(master) -> %d(worker)",
       pid, fd[0], fd[1]);
     close(fd[1]);
-    t->child[core] = fd[0];
+    pfp->fd = fd[0];
+    pfp->pid = pid;
+    ++pfp;
   }
+
+  *pspin = spin;
 }
 
 /*********************************************************************/
 
-static inline void WaitCores(cThis *t)
+Extern void SUFFIX(cubawait)(Spin **pspin)
 {
-  if( t->ncores >= MINCORES ) {
+  int cores, core, status;
-    int core;
+  Spin *spin;
-    pid_t pid;
+
-    for( core = 0; core < t->ncores; ++core ) {
+  MasterExit();
-      MASTER("closing fd %d", t->child[core]);
+
-      close(t->child[core]);
+  if( Invalid(pspin) || (spin = *pspin) == NULL ) return;
-    }
+
-    free(t->child);
+  cores = spin->spec.naccel + spin->spec.ncores;
-    for( core = 0; core < t->ncores; ++core ) {
+
-      MASTER("waiting for child");
+  for( core = 0; core < cores; ++core ) {
-      wait(&pid);
+    MASTER("closing fd %d", spin->fp[core].fd);
-      MASTER("pid %d terminated", pid);
+    close(spin->fp[core].fd);
-    }
   }
+
+#ifdef KILL_WORKERS
+  for( core = 0; core < cores; ++core ) {
+    MASTER("killing pid %d", spin->fp[core].pid);
+    kill(spin->fp[core].pid, SIGKILL);
+  }
+#endif
+
+  for( core = 0; core < cores; ++core ) {
+    DEB_ONLY(pid_t pid;)
+    MASTER("waiting for child");
+    DEB_ONLY(pid =) wait(&status);
+    MASTER("pid %d terminated with exit code %d", pid, status);
+  }
+
+  free(spin);
+  *pspin = NULL;
 }
 
+#else
+
+Extern void SUFFIX(cubafork)(Spin **pspin) {}
+
+Extern void SUFFIX(cubawait)(Spin **pspin)
+{
+  MasterExit();
+}
+
+#endif
+

src/external/libCuba/src/common/MSample.c

@@ -3,7 +3,7 @@
 		the sampling routine for the
 		Mathematica versions of the Cuba routines
 		by Thomas Hahn
-		last modified 19 Mar 12 th
+		last modified 13 Mar 15 th
 */
 
 
@@ -11,20 +11,20 @@ static void DoSample(This *t, cnumber n, real *x, real *f
   VES_ONLY(, real *w, ccount iter))
 {
   real *mma_f;
-  long mma_n;
+  int mma_n;
 
   if( MLAbort ) longjmp(t->abort, -99);
 
   MLPutFunction(stdlink, "EvaluatePacket", 1);
   MLPutFunction(stdlink, "Cuba`" ROUTINE "`sample", 1 VES_ONLY(+2) DIV_ONLY(+1));
-  MLPutRealList(stdlink, x, n*t->ndim);
+  MLPutRealxList(stdlink, x, n*t->ndim);
-  VES_ONLY(MLPutRealList(stdlink, w, n);
+  VES_ONLY(MLPutRealxList(stdlink, w, n);
            MLPutInteger(stdlink, iter);)
   DIV_ONLY(MLPutInteger(stdlink, t->phase);)
   MLEndPacket(stdlink);
 
   MLNextPacket(stdlink);
-  if( !MLGetRealList(stdlink, &mma_f, &mma_n) ) {
+  if( !MLGetRealxList(stdlink, &mma_f, &mma_n) ) {
     MLClearError(stdlink);
     MLNewPacket(stdlink);
     longjmp(t->abort, -99);
@@ -33,12 +33,12 @@ static void DoSample(This *t, cnumber n, real *x, real *f
   t->neval += mma_n;
 
   if( mma_n != n*t->ncomp ) {
-    MLDisownRealList(stdlink, mma_f, mma_n);
+    MLReleaseRealxList(stdlink, mma_f, mma_n);
     longjmp(t->abort, -3);
   }
  
   Copy(f, mma_f, n*t->ncomp);
-  MLDisownRealList(stdlink, mma_f, mma_n);
+  MLReleaseRealxList(stdlink, mma_f, mma_n);
 }
 
 /*********************************************************************/
@@ -50,16 +50,16 @@ static count SampleExtra(This *t, cBounds *b)
 {
   count n, nget;
   real *mma_f;
-  long mma_n;
+  int mma_n;
 
   MLPutFunction(stdlink, "EvaluatePacket", 1);
   MLPutFunction(stdlink, "Cuba`Divonne`findpeak", 2);
-  MLPutRealList(stdlink, (real *)b, 2*t->ndim);
+  MLPutRealxList(stdlink, (real *)b, 2*t->ndim);
   MLPutInteger(stdlink, t->phase);
   MLEndPacket(stdlink);
 
   MLNextPacket(stdlink);
-  if( !MLGetRealList(stdlink, &mma_f, &mma_n) ) {
+  if( !MLGetRealxList(stdlink, &mma_f, &mma_n) ) {
     MLClearError(stdlink);
     MLNewPacket(stdlink);
     longjmp(t->abort, -99);
@@ -73,7 +73,7 @@ static count SampleExtra(This *t, cBounds *b)
     Copy(t->fextra, mma_f + nget*t->ndim, n*t->ncomp);
   }
 
-  MLDisownRealList(stdlink, mma_f, mma_n);
+  MLReleaseRealxList(stdlink, mma_f, mma_n);
 
   return n;
 }

src/external/libCuba/src/common/Makefile.am

@@ -1,12 +1,14 @@
 ## Process this file with automake to create Makefile.in
 
-c_sources = WorkerIni.c
+c_sources = \
+  Global.c \
+  Data.c
 
 AM_CPPFLAGS = -I. -I.. -I../common -DNOUNDERSCORE
 AM_CFLAGS = $(LOCAL_CUBA_LIB_CFLAGS)
 AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS)
 
-noinst_LTLIBRARIES = libworkerini.la
+noinst_LTLIBRARIES = libcommon.la
 
-libworkerini_la_SOURCES = $(c_sources)
+libcommon_la_SOURCES = $(c_sources)
-libworkerini_la_LDFLAGS = $(AM_LDFLAGS)
+libcommon_la_LDFLAGS = $(AM_LDFLAGS)

src/external/libCuba/src/common/Random.c

@@ -1,7 +1,7 @@
 /*
 	Random.c
 		quasi- and pseudo-random-number generation
-		last modified 7 Aug 13 th
+		last modified 18 Mar 14 th
 */
 
 
@@ -88,14 +88,13 @@ static inline void SobolIni(This *t)
     299,   1,   3,   3,   9,   9,  25, 107,  39 };
 
   count dim, bit, nbits;
-  number max, *pini = ini;
+  number *pini = ini, max;
-  cnumber nmax = 2*t->maxeval;
 
-  for( nbits = 0, max = 1; max <= nmax; max <<= 1 ) ++nbits;
+  for( nbits = 0, max = t->maxeval; max; max >>= 1 ) ++nbits;
-  t->rng.sobol.norm = 1./max;
+  t->rng.sobol.norm = ldexp(.5, -nbits);
 
-  for( bit = 0; bit < nbits; ++bit )
+  for( bit = 0; bit <= nbits; ++bit )
-    t->rng.sobol.v[0][bit] = (max >>= 1);
+    t->rng.sobol.v[0][bit] = (number)1 << (nbits - bit);
 
   for( dim = 1; dim < t->ndim; ++dim ) {
     number *pv = t->rng.sobol.v[dim], *pvv = pv;
@@ -103,10 +102,10 @@ static inline void SobolIni(This *t)
     int inibits = -1, bit;
     for( j = powers; j; j >>= 1 ) ++inibits;
 
-    memcpy(pv, pini, inibits*sizeof(*pini));
+    memcpy(pv, pini, inibits*sizeof *pini);
     pini += 8;
 
-    for( bit = inibits; bit < nbits; ++bit ) {
+    for( bit = inibits; bit <= nbits; ++bit ) {
       number newv = *pvv, j = powers;
       int b;
       for( b = 0; b < inibits; ++b ) {
@@ -117,8 +116,8 @@ static inline void SobolIni(This *t)
       ++pvv;
     }
 
-    for( bit = 0; bit < nbits - 1; ++bit )
+    for( bit = 0; bit < nbits; ++bit )
-      pv[bit] <<= nbits - bit - 1;
+      pv[bit] <<= nbits - bit;
   }
 
   t->rng.sobol.seq = 0;

src/external/libCuba/src/common/WorkerIni.c

@@ -1,37 +0,0 @@
-/*
-	WorkerIni.c
-		set/run the init/exit functions for worker processes
-		by Thomas Hahn
-		last modified 6 Sep 12 th
-*/
-
-
-#include "stddecl.h"
-
-extern workerini cubaini;
-
-Extern void SUFFIX(cubasetinit)(subroutine f, void *arg)
-{
-  cubaini.initfun = f;
-  cubaini.initarg = arg;
-}
-
-
-Extern void SUFFIX(cubasetexit)(subroutine f, void *arg)
-{
-  cubaini.exitfun = f;
-  cubaini.exitarg = arg;
-}
-
-
-Extern void SUFFIX(cubaruninit)()
-{
-  if( cubaini.initfun ) cubaini.initfun(cubaini.initarg);
-}
-
-
-Extern void SUFFIX(cubarunexit)()
-{
-  if( cubaini.exitfun ) cubaini.exitfun(cubaini.exitarg);
-}
-

src/external/libCuba/src/common/stddecl.h

@@ -1,7 +1,7 @@
 /*
 	stddecl.h
 		declarations common to all Cuba routines
-		last modified 17 Sep 13 th
+		last modified 23 Apr 15 th
 */
 
 
@@ -13,7 +13,6 @@
 #endif
 
 #define _DEFAULT_SOURCE
-#define _XOPEN_SOURCE
 
 #include <stdio.h>
 #include <stdlib.h>
@@ -30,6 +29,7 @@
 #ifdef HAVE_FORK
 #include <sys/wait.h>
 #include <sys/socket.h>
+#include <signal.h>
 #ifdef HAVE_SHMGET
 #include <sys/ipc.h>
 #include <sys/shm.h>
@@ -83,10 +83,31 @@ void *alloca (size_t);
 
 #define SAMPLESIZE (NW + t->ndim + t->ncomp)*sizeof(real)
 
+
+enum { uninitialized = 0x61627563 };
+
+#define EnvInit(var, name, default) \
+  if( var == uninitialized ) { \
+    cchar *env = getenv(name); \
+    if( env == NULL ) var = default; \
+    else { \
+      var = atoi(env); \
+      if( cubaverb_ ) { \
+        char out[64]; \
+        sprintf(out, "env " name " = %d", (int)var); \
+        Print(out); \
+      } \
+    } \
+  }
+
+#define VerboseInit() EnvInit(cubaverb_, "CUBAVERBOSE", 0)
+#define MaxVerbose(flags) (flags + IDim(IMin(cubaverb_, 3) - ((flags) & 3)))
+
 #define VERBOSE (t->flags & 3)
 #define LAST (t->flags & 4)
 #define SHARPEDGES (t->flags & 8)
 #define KEEPFILE (t->flags & 16)
+#define ZAPSTATE (t->flags & 32)
 #define REGIONS (t->flags & 128)
 #define RNG (t->flags >> 8)
 
@@ -120,7 +141,7 @@ void *alloca (size_t);
 
 #define Zap(d) memset(d, 0, sizeof(d))
 
-#define MaxErr(avg) Max(t->epsrel*fabs(avg), t->epsabs)
+#define MaxErr(avg) Max(t->epsrel*fabsx(avg), t->epsabs)
 
 #ifdef __cplusplus
 #define mallocset(p, n) (*(void **)&p = malloc(n))
@@ -157,8 +178,21 @@ void *alloca (size_t);
 
 #ifdef MLVERSION
 #define ML_ONLY(...) __VA_ARGS__
+#define ML_NOT(...)
 #else
 #define ML_ONLY(...)
+#define ML_NOT(...) __VA_ARGS__
+
+#define CORE_MASTER (int []){32768}
+#define MasterInit() do if( !cubafun_.init ) { \
+  cubafun_.init = true; \
+  if( cubafun_.initfun ) cubafun_.initfun(cubafun_.initarg, CORE_MASTER); \
+} while( 0 )
+#define MasterExit() do if( cubafun_.init ) { \
+  cubafun_.init = false; \
+  if( cubafun_.exitfun ) cubafun_.exitfun(cubafun_.exitarg, CORE_MASTER); \
+} while( 0 )
+#define Invalid(s) ((s) == NULL || *(int *)(s) == -1)
 
 #ifdef HAVE_FORK
 #undef FORK_ONLY
@@ -168,37 +202,40 @@ void *alloca (size_t);
 #undef SHM_ONLY
 #define SHM_ONLY(...) __VA_ARGS__
 
-#define ShmMap(t, ...) if( t->shmid != -1 ) { \
+#define MasterAlloc(t) \
-  t->frame = shmat(t->shmid, NULL, 0); \
+  t->shmid = shmget(IPC_PRIVATE, t->nframe*SAMPLESIZE, IPC_CREAT | 0600)
-  if( t->frame == (void *)-1 ) Abort("shmat"); \
+#define MasterFree(t) shmctl(t->shmid, IPC_RMID, NULL)
-  __VA_ARGS__ \
+#define WorkerAlloc(t)
-}
+#define WorkerFree(r)
-
-#define ShmRm(t) shmctl(t->shmid, IPC_RMID, NULL);
 
 #undef ShmAlloc
-#define ShmAlloc(t, ...) \
+#define ShmAlloc(t, who) \
-  t->shmid = shmget(IPC_PRIVATE, t->nframe*SAMPLESIZE, IPC_CREAT | 0600); \
+  who##Alloc(t); \
-  ShmMap(t, __VA_ARGS__)
+  if( t->shmid != -1 ) { \
+    t->frame = shmat(t->shmid, NULL, 0); \
+    if( t->frame == (void *)-1 ) Abort("shmat"); \
+  }
 
 #undef ShmFree
-#define ShmFree(t, ...) if( t->shmid != -1 ) { \
+#define ShmFree(t, who) \
-  shmdt(t->frame); \
+  if( t->shmid != -1 ) { \
-  __VA_ARGS__ \
+    shmdt(t->frame); \
-}
+    who##Free(t); \
+  }
 
 #endif
 #endif
 #endif
   
-#define FrameAlloc(t, ...) \
+#define FrameAlloc(t, who) \
-  SHM_ONLY(ShmAlloc(t, __VA_ARGS__) else) \
+  SHM_ONLY(ShmAlloc(t, who) else) \
   MemAlloc(t->frame, t->nframe*SAMPLESIZE);
 
-#define FrameFree(t, ...) DIV_ONLY(if( t->nframe )) { \
+#define FrameFree(t, who) \
-  SHM_ONLY(ShmFree(t, __VA_ARGS__) else) \
+  DIV_ONLY(if( t->nframe )) { \
-  free(t->frame); \
+    SHM_ONLY(ShmFree(t, who) else) \
-}
+    free(t->frame); \
+  }
 
 
 #define StateDecl \
@@ -219,7 +256,9 @@ struct stat st
 
 typedef long long int signature_t;
 
-#define StateSignature(t, i) (0x41425543 + \
+enum { signature = 0x41425543 };
+
+#define StateSignature(t, i) (signature + \
   ((signature_t)(i) << 60) + \
   ((signature_t)(t)->ncomp << 48) + \
   ((signature_t)(t)->ndim << 32))
@@ -306,39 +345,100 @@ typedef const count ccount;
 #define PREFIX(s) ll##s
 #define NUMBER "%lld"
 #define NUMBER7 "%7lld"
+#define NUMBER_MAX LLONG_MAX
 typedef long long int number;
 #else
 #define PREFIX(s) s
 #define NUMBER "%d"
 #define NUMBER7 "%7d"
+#define NUMBER_MAX INT_MAX
 typedef int number;
 #endif
 typedef const number cnumber;
 
 #define REAL "%g"
 #define REALF "%f"
-typedef /*long*/ double real;
+#define SHOW(r) (double)(r)
-	/* Switching to long double is not as trivial as it
+	/* floating-point numbers are printed with SHOW */
-	   might seem here.  sqrt, erf, exp, pow need to be
+
-	   replaced by their long double versions (sqrtl, ...),
+#if REALSIZE == 16
-	   printf formats need to be updated similarly, and
+#include <quadmath.h>
-	   ferrying long doubles to Mathematica is of course
+typedef __float128 real;
-	   quite another matter, too. */
+#define RC(x) x##Q
+#define sqrtx sqrtq
+#define expx expq
+#define powx powq
+#define erfx erfq
+#define fabsx fabsq
+#define ldexpx ldexpq
+#define REAL_MAX_EXP FLT128_MAX_EXP
+#define REAL_MAX FLT128_MAX
+#elif REALSIZE == 10
+typedef long double real;
+#define RC(x) x##L
+#define sqrtx sqrtl
+#define expx expl
+#define powx powl
+#define erfx erfl
+#define fabsx fabsl
+#define ldexpx ldexpl
+#define REAL_MAX_EXP LDBL_MAX_EXP
+#define REAL_MAX LDBL_MAX
+#define MLPutRealxList MLPutReal128List
+#define MLGetRealxList MLGetReal128List
+#define MLReleaseRealxList MLReleaseReal128List
+#else
+typedef double real;
+#define RC(x) x
+#define sqrtx sqrt
+#define expx exp
+#define powx pow
+#define erfx erf
+#define fabsx fabs
+#define ldexpx ldexp
+#define REAL_MAX_EXP DBL_MAX_EXP
+#define REAL_MAX DBL_MAX
+#define MLPutRealxList MLPutReal64List
+#define MLGetRealxList MLGetReal64List
+#define MLReleaseRealxList MLReleaseReal64List
+#endif
 
 typedef const real creal;
 
-typedef void (*subroutine)();
+typedef void (*subroutine)(void *, cint *);
 
 typedef struct {
   subroutine initfun;
   void *initarg;
   subroutine exitfun;
   void *exitarg;
-} workerini;
+  bool init;
+} coreinit;
+
+typedef struct {
+  int ncores, naccel;
+  int pcores, paccel;
+} corespec;
+
+typedef struct {
+  int fd, pid;
+} fdpid;
+
+typedef struct {
+  corespec spec;
+  fdpid fp[];
+} Spin;
 
 
 struct _this;
 
+typedef struct {
+  void (*worker)(struct _this *, csize_t, cint, cint);
+  struct _this *thisptr;
+  size_t thissize;
+} dispatch;
+
+
 typedef unsigned int state_t;
 
 #define SOBOL_MINDIM 1
@@ -407,7 +507,7 @@ static inline real Max(creal a, creal b) {
 }
 
 static inline real Weight(creal sum, creal sqsum, cnumber n) {
-  creal w = sqrt(sqsum*n);
+  creal w = sqrtx(sqsum*n);
   return (n - 1)/Max((w + sum)*(w - sum), NOTZERO);
 }
 
@@ -452,7 +552,7 @@ static inline void Print(MLCONST char *s)
 
 #else
 
-#define Print(s) puts(s); fflush(stdout)
+#define Print(s) { puts(s); fflush(stdout); }
 
 #endif
 

src/external/libCuba/src/cuba.h

@@ -2,116 +2,122 @@
 	cuba.h
 		Prototypes for the Cuba library
 		this file is part of Cuba
-		last modified 30 Apr 13 th
+		last modified 13 Mar 15 th
 */
 
+typedef double cubareal;
+
+	/* integrand_t is intentionally a minimalistic integrand type.
+	   It includes neither the nvec and core arguments nor the
+	   extra arguments passed by Vegas/Suave (weight, iter) and
+	   Divonne (phase).
+	   In most cases, integrand_t is just what you want, otherwise
+	   simply use an explicit typecast to integrand_t in the Cuba
+	   invocation. */
+typedef int (*integrand_t)(const int *ndim, const cubareal x[],
+  const int *ncomp, cubareal f[], void *userdata);
+
+typedef void (*peakfinder_t)(const int *ndim, const cubareal b[],
+  int *n, cubareal x[], void *userdata);
+
 #ifdef __cplusplus
 extern "C" {
 #endif
 
-	/* NB: Divonne actually passes a fifth argument, a const int *
-	   which points to the integration phase.  This is used only
-	   rarely and most users are confused by the warnings the
-	   compiler emits if the `correct' prototype is used.  Thus,
-	   if you need to access this argument, use an explicit cast
-	   to integrand_t when invoking Divonne. */
-typedef int (*integrand_t)(const int *ndim, const double x[],
-  const int *ncomp, double f[], void *userdata);
-
-typedef void (*peakfinder_t)(const int *ndim, const double b[],
-  int *n, double x[]);
-
 void Vegas(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
+  integrand_t integrand, void *userdata, const int nvec,
-  const double epsrel, const double epsabs,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags, const int seed,
   const int mineval, const int maxeval,
   const int nstart, const int nincrease, const int nbatch,
-  const int gridno, const char *statefile,
+  const int gridno, const char *statefile, void *spin,
   int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
 void llVegas(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
+  integrand_t integrand, void *userdata, const long long int nvec,
-  const double epsrel, const double epsabs,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags, const int seed,
   const long long int mineval, const long long int maxeval,
   const long long int nstart, const long long int nincrease,
   const long long int nbatch,
-  const int gridno, const char *statefile,
+  const int gridno, const char *statefile, void *spin,
   long long int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
 void Suave(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
+  integrand_t integrand, void *userdata, const int nvec,
-  const double epsrel, const double epsabs,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags, const int seed,
   const int mineval, const int maxeval,
-  const int nnew, const double flatness,
+  const int nnew, const int nmin,
-  const char *statefile,
+  const cubareal flatness, const char *statefile, void *spin,
   int *nregions, int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
 void llSuave(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
+  integrand_t integrand, void *userdata, const long long int nvec,
-  const double epsrel, const double epsabs,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags, const int seed,
   const long long int mineval, const long long int maxeval,
-  const long long int nnew, const double flatness,
+  const long long int nnew, const long long int nmin,
-  const char *statefile,
+  const cubareal flatness, const char *statefile, void *spin,
   int *nregions, long long int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
 void Divonne(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
+  integrand_t integrand, void *userdata, const int nvec,
-  const double epsrel, const double epsabs,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags, const int seed,
   const int mineval, const int maxeval,
   const int key1, const int key2, const int key3, const int maxpass,
-  const double border, const double maxchisq, const double mindeviation,
+  const cubareal border, const cubareal maxchisq, const cubareal mindeviation,
-  const int ngiven, const int ldxgiven, double xgiven[],
+  const int ngiven, const int ldxgiven, cubareal xgiven[],
   const int nextra, peakfinder_t peakfinder,
-  const char *statefile,
+  const char *statefile, void *spin,
   int *nregions, int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
 void llDivonne(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
+  integrand_t integrand, void *userdata, const long long int nvec,
-  const double epsrel, const double epsabs,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags, const int seed,
   const long long int mineval, const long long int maxeval,
   const int key1, const int key2, const int key3, const int maxpass,
-  const double border, const double maxchisq, const double mindeviation,
+  const cubareal border, const cubareal maxchisq, const cubareal mindeviation,
-  const long long int ngiven, const int ldxgiven, double xgiven[],
+  const long long int ngiven, const int ldxgiven, cubareal xgiven[],
-  const long long int nextra,
+  const long long int nextra, peakfinder_t peakfinder,
-  void (*peakfinder)(const int *, const double [], int *, double []),
+  const char *statefile, void *spin,
-  const char *statefile,
   int *nregions, long long int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
 void Cuhre(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
+  integrand_t integrand, void *userdata, const int nvec,
-  const double epsrel, const double epsabs,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags, const int mineval, const int maxeval,
   const int key,
-  const char *statefile,
+  const char *statefile, void *spin,
   int *nregions, int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
 void llCuhre(const int ndim, const int ncomp,
-  integrand_t integrand, void *userdata,
+  integrand_t integrand, void *userdata, const long long int nvec,
-  const double epsrel, const double epsabs,
+  const cubareal epsrel, const cubareal epsabs,
   const int flags,
   const long long int mineval, const long long int maxeval,
   const int key,
-  const char *statefile,
+  const char *statefile, void *spin,
   int *nregions, long long int *neval, int *fail,
-  double integral[], double error[], double prob[]);
+  cubareal integral[], cubareal error[], cubareal prob[]);
 
-void cubasetinit(void (*)(), void *);
+void cubafork(void *pspin);
-void cubasetexit(void (*)(), void *);
+void cubawait(void *pspin);
-void cubaruninit(void);
+
-void cubaruninit(void);
+void cubacores(const int n, const int p);
+void cubaaccel(const int n, const int p);
+
+void cubainit(void (*f)(), void *arg);
+void cubaexit(void (*f)(), void *arg);
 
 #ifdef __cplusplus
 }

src/external/libCuba/src/cuhre/Cuhre.c

@@ -2,7 +2,7 @@
 	Cuhre.c
 		Adaptive integration using cubature rules
 		by Thomas Hahn
-		last modified 17 Sep 13 th
+		last modified 22 Jul 14 th
 */
 
 
@@ -15,55 +15,70 @@
 /*********************************************************************/
 
 Extern void EXPORT(Cuhre)(ccount ndim, ccount ncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber nvec,
   creal epsrel, creal epsabs,
   cint flags, cnumber mineval, cnumber maxeval,
-  ccount key, cchar *statefile,
+  ccount key, cchar *statefile, Spin **pspin,
   count *pnregions, number *pneval, int *pfail,
   real *integral, real *error, real *prob)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = ndim;
   t.ncomp = ncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = nvec;
   t.epsrel = epsrel;
   t.epsabs = epsabs;
-  t.flags = flags;
+  t.flags = MaxVerbose(flags);
   t.mineval = mineval;
   t.maxeval = maxeval;
   t.key = key;
   t.statefile = statefile;
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
  
   *pfail = Integrate(&t, integral, error, prob);
   *pnregions = t.nregions;
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
 
 /*********************************************************************/
 
 Extern void EXPORT(cuhre)(ccount *pndim, ccount *pncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber *pnvec,
   creal *pepsrel, creal *pepsabs,
   cint *pflags, cnumber *pmineval, cnumber *pmaxeval,
-  ccount *pkey, cchar *statefile,
+  ccount *pkey, cchar *statefile, Spin **pspin,
   count *pnregions, number *pneval, int *pfail,
   real *integral, real *error, real *prob, cint statefilelen)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = *pndim;
   t.ncomp = *pncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = *pnvec;
   t.epsrel = *pepsrel;
   t.epsabs = *pepsabs;
-  t.flags = *pflags;
+  t.flags = MaxVerbose(*pflags);
   t.mineval = *pmineval;
   t.maxeval = *pmaxeval;
   t.key = *pkey;
   CString(t.statefile, statefile, statefilelen);
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
 
   *pfail = Integrate(&t, integral, error, prob);
   *pnregions = t.nregions;
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
+

src/external/libCuba/src/cuhre/Integrate.c

@@ -3,7 +3,7 @@
 		integrate over the unit hypercube
 		this file is part of Cuhre
 		checkpointing by B. Chokoufe
-		last modified 17 Sep 13 th
+		last modified 14 Mar 15 th
 */
 
 
@@ -11,6 +11,9 @@
 
 typedef struct pool {
   struct pool *next;
+#if REALSIZE > 8
+  void *dummy; /* for alignment */
+#endif
   char region[];
 } Pool;
 
@@ -42,12 +45,14 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   if( VERBOSE > 1 ) {
     sprintf(out, "Cuhre input parameters:\n"
       "  ndim " COUNT "\n  ncomp " COUNT "\n"
+      ML_NOT("  nvec " NUMBER "\n")
       "  epsrel " REAL "\n  epsabs " REAL "\n"
       "  flags %d\n  mineval " NUMBER "\n  maxeval " NUMBER "\n"
       "  key " COUNT "\n"
       "  statefile \"%s\"",
       t->ndim, t->ncomp,
-      t->epsrel, t->epsabs,
+      ML_NOT(t->nvec,)
+      SHOW(t->epsrel), SHOW(t->epsabs),
       t->flags, t->mineval, t->maxeval,
       t->key,
       t->statefile);
@@ -61,7 +66,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
 
   RuleAlloc(t);
   t->mineval = IMax(t->mineval, t->rule.n + 1);
-  FrameAlloc(t, ShmRm(t));
+  FrameAlloc(t, Master);
   ForkCores(t);
 
   if( (fail = setjmp(t->abort)) ) goto abort;
@@ -125,7 +130,8 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       for( tot = state->totals, comp = 0; tot < Tot; ++tot )
         oe += sprintf(oe, "\n[" COUNT "] "
           REAL " +- " REAL "  \tchisq " REAL " (" COUNT " df)",
-          ++comp, tot->avg, tot->err, tot->chisq, t->nregions - 1);
+          ++comp, SHOW(tot->avg), SHOW(tot->err),
+          SHOW(tot->chisq), t->nregions - 1);
       Print(out);
     }
 
@@ -188,7 +194,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
 
       tot->lastavg += diff = resL->avg + resR->avg - res->avg;
 
-      diff = fabs(.25*diff);
+      diff = fabsx(.25*diff);
       err = resL->err + resR->err;
       if( err > 0 ) {
         creal c = 1 + 2*diff/err;
@@ -213,7 +219,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       }
       else {
         tot->avg = avg;
-        tot->err = sqrt(sigsq);
+        tot->err = sqrtx(sigsq);
       }
     }
     ++t->nregions;
@@ -249,13 +255,13 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
         Result *Res;
 
         MLPutFunction(stdlink, "Cuba`Cuhre`region", 2);
-        MLPutRealList(stdlink, (real *)region->bounds, 2*t->ndim);
+        MLPutRealxList(stdlink, (real *)region->bounds, 2*t->ndim);
 
         MLPutFunction(stdlink, "List", t->ncomp);
         for( Res = (res = RegionResult(region)) + t->ncomp;
              res < Res; ++res ) {
           real r[] = {res->avg, res->err};
-          MLPutRealList(stdlink, r, Elements(r));
+          MLPutRealxList(stdlink, r, Elements(r));
         }
       }
   }
@@ -266,9 +272,7 @@ abort:
     cur = cur->next;
     free(pool);
   }
-
+  FrameFree(t, Master);
-  WaitCores(t);
-  FrameFree(t);
   RuleFree(t);
 
   StateRemove(t);

src/external/libCuba/src/cuhre/Rule.c

@@ -4,11 +4,12 @@
 		code lifted with minor modifications from DCUHRE
 		by J. Berntsen, T. Espelid, and A. Genz
 		this file is part of Cuhre
-		last modified 5 Aug 13 th
+		last modified 7 May 15 th
 */
 
 
 #define NextSet(p) p = (Set *)((char *)p + setsize)
+#define IndexSet(p, n) ((Set *)((char *)p + n*setsize))
 
 /*********************************************************************/
 
@@ -153,7 +154,7 @@ static void Rule13Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -298,7 +299,7 @@ static void Rule11Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -309,32 +310,55 @@ static void Rule11Alloc(This *t)
 static void Rule9Alloc(This *t)
 {
   static creal w[] = {
-    -.0023611709677855117884,   .11415390023857325268,
+    RC(-.002361170967785511788400941242259231309691),
-    -.63833920076702389094,     .74849988504685208004,
+    RC(.1141539002385732526821323741697655347686),
-    -.0014324017033399125142,   .057471507864489725949,
+    RC(-.6383392007670238909386026193674701393074),
-    -.14225104571434243234,    -.062875028738286979989,
+    RC(.7484998850468520800423030047583803945205),
-     .254591133248959089,     -1.207328566678236261,
+    RC(-.001432401703339912514196154599769007103671),
-     .89567365764160676508,    -.36479356986049146661,
+    RC(.05747150786448972594860897296200006759892),
-     .0035417564516782676826,  -.072609367395893679605,
+    RC(-.1422510457143424323449521620935950679394),
-     .10557491625218991012,     .0021486025550098687713,
+    RC(-.06287502873828697998942424881040490136987),
-    -.032268563892953949998,    .010636783990231217481,
+    RC(.2545911332489590890011611142429070613156),
-     .014689102496143490175,    .51134708346467591431,
+    RC(-1.207328566678236261002219995185143356737),
-     .45976448120806344646,     .18239678493024573331,
+    RC(.8956736576416067650809467826488567200939),
-    -.04508628929435784076,     .21415883524352793401,
+    RC(-.3647935698604914666100134551377381205297),
-    -.027351546526545644722,    .054941067048711234101,
+    RC(.003541756451678267682601411863388846964536),
-     .11937596202570775297,     .65089519391920250593,
+    RC(-.07260936739589367960492815865074633743652),
-     .14744939829434460168,     .057693384490973483573,
+    RC(.1055749162521899101218622863269817454540),
-     .034999626602143583822,  -1.3868627719278281436,
+    RC(.002148602555009868771294231899653510655506),
-    -.2386668732575008879,      .015532417276607053264,
+    RC(-.03226856389295394999786630399875134318006),
-     .0035328099607090870236,   .09231719987444221619,
+    RC(.01063678399023121748083624225818915724455),
-     .02254314464717892038,     .013675773263272822361,
+    RC(.01468910249614349017540783437728097691502),
-    -.32544759695960125297,     .0017708782258391338413,
+    RC(.5113470834646759143109387357149329909126),
-     .0010743012775049343856,   .25150011495314791996 };
+    RC(.4597644812080634464633352781605214342691),
+    RC(.1823967849302457333050067275688690602649),
+    RC(-.04508628929435784075980562738240804429658),
+    RC(.2141588352435279340097929526588394300172),
+    RC(-.02735154652654564472203690086290223507436),
+    RC(.05494106704871123410060080562462135546101),
+    RC(.1193759620257077529708962121565290178730),
+    RC(.6508951939192025059314756320878023215278),
+    RC(.1474493982943446016775696826942585013243),
+    RC(.05769338449097348357291272840392627722165),
+    RC(.03499962660214358382244159694487155861542),
+    RC(-1.386862771927828143599782668709014266770),
+    RC(-.2386668732575008878964134721962088068396),
+    RC(.01553241727660705326386197156586357005224),
+    RC(.003532809960709087023561817517751309380604),
+    RC(.09231719987444221619017126187763868745587),
+    RC(.02254314464717892037990281369120402214829),
+    RC(.01367577326327282236101845043145111753718),
+    RC(-.3254475969596012529657378160439011607639),
+    RC(.001770878225839133841300705931694423482268),
+    RC(.001074301277504934385647115949826755327753),
+    RC(.2515001149531479199576969952416196054795) };
 
   static creal g[] = {
-     .47795365790226950619,     .20302858736911986780,
+    RC(.4779536579022695061928604197171830064732),
-     .44762735462617812882,     .125,
+    RC(.2030285873691198677998034402373279133258),
-     .34303789878087814570 };
+    RC(.4476273546261781288207704806530998539285),
+    RC(.125),
+    RC(.3430378987808781457001426145164678603407) };
 
   enum { nsets = 9 };
 
@@ -440,7 +464,7 @@ static void Rule9Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -451,21 +475,33 @@ static void Rule9Alloc(This *t)
 static void Rule7Alloc(This *t)
 {
   static creal w[] = {
-     .019417866674748388428,   -.40385257701150182546,
+    RC(.01941786667474838842844534313920462333850),
-     .64485668767465982223,     .01177982690775806141,
+    RC(-.4038525770115018254611834753723880293161),
-    -.18041318740733609012,    -.088785828081335044443,
+    RC(.6448566876746598222277360730193089551024),
-     .056328645808285941374,   -.0097089333373741942142,
+    RC(.01177982690775806141012214458820955067854),
-    -.99129176779582358138,    -.17757165616267008889,
+    RC(-.1804131874073360901182293138710989490609),
-     .12359398032043233572,     .074978148702033690681,
+    RC(-.08878582808133504444306598174517276122439),
-     .55489147051423559776,     .088041241522692771226,
+    RC(.05632864580828594137378124255408286479947),
-     .021118358455513385083,   -.0099302203239653333087,
+    RC(-.009708933337374194214222671569602311669249),
-    -.064100053285010904179,    .030381729038221007659,
+    RC(-.9912917677958235813775106862002319060386),
-     .0058899134538790307051,  -.0048544666686870971071,
+    RC(-.1775716561626700888861319634903455224488),
-     .35514331232534017777 };
+    RC(.1235939803204323357183625846672135876752),
+    RC(.07497814870203369068087999555157339703666),
+    RC(.5548914705142355977605994477355651401434),
+    RC(.08804124152269277122645182458858273865209),
+    RC(.02111835845551338508329573367808085283304),
+    RC(-.009930220323965333308685820460105538586058),
+    RC(-.06410005328501090417895544042025034295870),
+    RC(.03038172903822100765927778829870429682489),
+    RC(.005889913453879030705061072294104775339268),
+    RC(-.004854466668687097107111335784801155834624),
+    RC(.3551433123253401777722639269806910448976) };
 
   static creal g[] = {
-     .47795365790226950619,     .20302858736911986780,
+    RC(.4779536579022695061928604197171830064732),
-     .375,                      .34303789878087814570 };
+    RC(.2030285873691198677998034402373279133258),
+    RC(.375),
+    RC(.3430378987808781457001426145164678603407) };
 
   enum { nsets = 6 };
 
@@ -541,7 +577,7 @@ static void Rule7Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -630,7 +666,8 @@ static void Sample(This *t, Region *region)
   Bounds *b, *B = region->bounds + t->ndim;
   Result *result = RegionResult(region), *res, *Res = result + t->ncomp;
   creal *errcoeff = t->rule.errcoeff;
-  creal ratio = Sq(first[2].gen[0]/first[1].gen[0]);
+  creal ratio = Sq(IndexSet(first,2)->gen[0]/
+                   IndexSet(first,1)->gen[0]);
 
   ccount offset = 2*t->ndim*t->ncomp;
   count dim, rul, n, maxdim = 0;
@@ -659,7 +696,7 @@ static void Sample(This *t, Region *region)
     for( dim = 0; dim < t->ndim; ++dim ) {
       creal *fp = f1 + t->ncomp;
       creal *fm = fp + t->ncomp;
-      creal fourthdiff = fabs(base +
+      creal fourthdiff = fabsx(base +
         ratio*(fp[0] + fm[0]) - (fp[offset] + fm[offset]));
       f1 = fm;
       if( fourthdiff > maxdiff ) {
@@ -688,7 +725,7 @@ static void Sample(This *t, Region *region)
       real maxerr = 0;
       for( s = first; s <= last; NextSet(s) )
         maxerr = Max(maxerr,
-          fabs(sum[rul + 1] + s->scale[rul]*sum[rul])*s->norm[rul]);
+          fabsx(sum[rul + 1] + s->scale[rul]*sum[rul])*s->norm[rul]);
       sum[rul] = maxerr;
     }
 
@@ -712,7 +749,7 @@ static void Sample(This *t, Region *region)
 
     for( res = result, comp = 0; res < Res; ++res )
       oe += sprintf(oe, "\n[" COUNT "] "
-        REAL " +- " REAL, ++comp, res->avg, res->err);
+        REAL " +- " REAL, ++comp, SHOW(res->avg), SHOW(res->err));
 
     Print(out);
   }

src/external/libCuba/src/cuhre/decl.h

@@ -2,7 +2,7 @@
 	decl.h
 		Type declarations
 		this file is part of Cuhre
-		last modified 26 Jul 13 th
+		last modified 21 Jul 14 th
 */
 
 
@@ -47,16 +47,18 @@ typedef struct {
 
 typedef const Rule cRule;
 
-typedef int (*Integrand)(ccount *, creal *, ccount *, real *, void *);
+typedef int (*Integrand)(ccount *, creal *, ccount *, real *,
+  void *, cnumber *, cint *);
 
 typedef struct _this {
   count ndim, ncomp;
 #ifndef MLVERSION
   Integrand integrand;
   void *userdata;
+  number nvec;
 #ifdef HAVE_FORK
-  int ncores, *child;
   SHM_ONLY(int shmid;)
+  Spin *spin;
 #endif
 #endif
   real *frame;

src/external/libCuba/src/divonne/Divonne.c

@@ -4,7 +4,7 @@
 		originally by J.H. Friedman and M.H. Wright
 		(CERNLIB subroutine D151)
 		this version by Thomas Hahn
-		last modified 17 Sep 13 th
+		last modified 22 Jul 14 th
 */
 
 #define DIVONNE
@@ -16,7 +16,7 @@
 /*********************************************************************/
 
 Extern void EXPORT(Divonne)(ccount ndim, ccount ncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber nvec,
   creal epsrel, creal epsabs,
   cint flags, cint seed,
   cnumber mineval, cnumber maxeval,
@@ -24,18 +24,22 @@ Extern void EXPORT(Divonne)(ccount ndim, ccount ncomp,
   creal border, creal maxchisq, creal mindeviation,
   cnumber ngiven, ccount ldxgiven, real *xgiven,
   cnumber nextra, PeakFinder peakfinder,
-  cchar *statefile,
+  cchar *statefile, Spin **pspin,
   int *pnregions, number *pneval, int *pfail,
   real *integral, real *error, real *prob)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = ndim;
   t.ncomp = ncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = nvec;
   t.epsrel = epsrel;
   t.epsabs = epsabs;
-  t.flags = flags;
+  t.flags = MaxVerbose(flags);
   t.seed = seed;
   t.mineval = mineval;
   t.maxeval = maxeval;
@@ -52,16 +56,19 @@ Extern void EXPORT(Divonne)(ccount ndim, ccount ncomp,
   t.nextra = nextra;
   t.peakfinder = peakfinder;
   t.statefile = statefile;
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
 
   *pfail = Integrate(&t, integral, error, prob);
   *pnregions = t.nregions;
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
 
 /*********************************************************************/
 
 Extern void EXPORT(divonne)(ccount *pndim, ccount *pncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber *pnvec,
   creal *pepsrel, creal *pepsabs,
   cint *pflags, cint *pseed,
   cnumber *pmineval, cnumber *pmaxeval,
@@ -69,18 +76,22 @@ Extern void EXPORT(divonne)(ccount *pndim, ccount *pncomp,
   creal *pborder, creal *pmaxchisq, creal *pmindeviation,
   cnumber *pngiven, ccount *pldxgiven, real *xgiven,
   cnumber *pnextra, PeakFinder peakfinder,
-  cchar *statefile,
+  cchar *statefile, Spin **pspin,
   int *pnregions, number *pneval, int *pfail,
   real *integral, real *error, real *prob, cint statefilelen)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = *pndim;
   t.ncomp = *pncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = *pnvec;
   t.epsrel = *pepsrel;
   t.epsabs = *pepsabs;
-  t.flags = *pflags;
+  t.flags = MaxVerbose(*pflags);
   t.seed = *pseed;
   t.mineval = *pmineval;
   t.maxeval = *pmaxeval;
@@ -97,9 +108,12 @@ Extern void EXPORT(divonne)(ccount *pndim, ccount *pncomp,
   t.nextra = *pnextra;
   t.peakfinder = peakfinder;
   CString(t.statefile, statefile, statefilelen);
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
 
   *pfail = Integrate(&t, integral, error, prob);
   *pnregions = t.nregions;
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
 

src/external/libCuba/src/divonne/Explore.c

@@ -2,7 +2,7 @@
 	Explore.c
 		sample region, determine min and max, split if necessary
 		this file is part of Divonne
-		last modified 2 Aug 13 th
+		last modified 12 Mar 15 th
 */
 
 
@@ -19,10 +19,11 @@ static int ExploreSerial(This *t, ccount iregion)
   Region *region = RegionPtr(iregion);
   cBounds *bounds = region->bounds;
   Result *result = RegionResult(region);
+  real *minmax = RegionMinMax(region);
 
   Vector(Extrema, extrema, NCOMP);
   Vector(real, xtmp, NDIM);
-  Result *r, *r0;
+  Result *r;
   creal *x;
   real *f;
   real halfvol, maxerr;
@@ -102,7 +103,7 @@ skip:
       ftmp = FindMinimum(t, bounds, xtmp, e->fmin);
       if( ftmp < r->fmin ) {
         r->fmin = ftmp;
-        XCopy(&r->xminmax[0], xtmp);
+        XCopy(&minmax[2*comp*t->ndim], xtmp);
       }
 
       t->selectedcomp = Tag(comp);
@@ -110,12 +111,12 @@ skip:
       ftmp = -FindMinimum(t, bounds, xtmp, -e->fmax);
       if( ftmp > r->fmax ) {
         r->fmax = ftmp;
-        XCopy(&r->xminmax[t->ndim], xtmp);
+        XCopy(&minmax[(2*comp + 1)*t->ndim], xtmp);
       }
     }
 
     r->spread = halfvol*(r->fmax - r->fmin);
-    err = r->spread/Max(fabs(r->avg), NOTZERO);
+    err = r->spread/Max(fabsx(r->avg), NOTZERO);
     if( err > maxerr ) {
       maxerr = err;
       maxcomp = comp;
@@ -130,22 +131,22 @@ skip:
   }
 
   region->cutcomp = maxcomp;
-  r0 = RegionResult(region);
+  r = RegionResult(region) + maxcomp;
-  r = r0 + maxcomp;
   if( halfvol*(r->fmin + r->fmax) > r->avg ) {
     region->fminor = r->fmin;
     region->fmajor = r->fmax;
-    region->xmajor = &r->xminmax[t->ndim] - (real *)r0;
+    region->xmajor = (2*maxcomp + 1)*t->ndim;
   }
   else {
     region->fminor = r->fmax;
     region->fmajor = r->fmin;
-    region->xmajor = &r->xminmax[0] - (real *)r0;
+    region->xmajor = 2*maxcomp*t->ndim;
   }
 
   if( region->isamples == 0 ) {
     if( (region->depth < INIDEPTH && r->spread < samples->neff*r->err) ||
-        r->spread < t->totals[maxcomp].secondspread ) region->depth = 0;
+        r->spread < t->totals[maxcomp].secondspread )
+      region->depth = 0;
     if( region->depth == 0 )
       for( comp = 0; comp < t->ncomp; ++comp )
         t->totals[comp].secondspread =

src/external/libCuba/src/divonne/FindMinimum.c

@@ -2,7 +2,7 @@
 	FindMinimum.c
 		find minimum (maximum) of hyperrectangular region
 		this file is part of Divonne
-		last modified 7 Aug 13 th
+		last modified 12 Mar 15 th
 */
 
 
@@ -42,7 +42,7 @@ static inline real Dot(ccount n, creal *a, creal *b)
 
 static inline real Length(ccount n, creal *vec)
 {
-  return sqrt(Dot(n, vec, vec));
+  return sqrtx(Dot(n, vec, vec));
 }
 
 /*********************************************************************/
@@ -120,7 +120,7 @@ static void UpdateCholesky(cThis *t, ccount n, real *hessian,
     p[i] = dir;
     gamma += Sq(dir)/Hessian(i, i);
   }
-  gamma = Max(fabs(1 - gamma), EPS);
+  gamma = Max(fabsx(1 - gamma), EPS);
 
   while( --i >= 0 ) {
     creal dir = z[i] = p[i];
@@ -154,7 +154,7 @@ static inline void BFGS(cThis *t, ccount n, real *hessian,
 
   c = Dot(n, g, p);
   if( c >= 0 ) return;
-  c = 1/sqrt(-c);
+  c = 1/sqrtx(-c);
   for( i = 0; i < n; ++i )
     y[i] = c*g[i];
   UpdateCholesky(t, n, hessian, y, p);
@@ -199,7 +199,7 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
      c) the gradient is positive, i.e. we'd move uphill */
 
   if( step > 0 && range > tol2 && grad <= 0 ) {
-    creal eps = RTEPS*fabs(range) + ftol;
+    creal eps = RTEPS*fabsx(range) + ftol;
     creal mingrad = -1e-4*grad, maxgrad = -gtol*grad;
 
     real end = range + eps;
@@ -225,7 +225,7 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
         maxstep = maxstep*(1 + .75*RTEPS) + .75*tol;
       }
 
-      cur.dx = (fabs(step) >= tol) ? step : (step > 0) ? tol : -tol;
+      cur.dx = (fabsx(step) >= tol) ? step : (step > 0) ? tol : -tol;
       dist = distmin + cur.dx;
       for( i = 0; i < nfree; ++i ) {
         ccount dim = ifree[i];
@@ -248,7 +248,7 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
         if( cur.dx < 0 ) b = w;
         else a = w;
 
-        tol = RTEPS*fabs(distmin) + ftol;
+        tol = RTEPS*fabsx(distmin) + ftol;
         tol2 = tol + tol;
       }
       else {
@@ -260,14 +260,14 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
 
       if( distmin + b.dx <= xtol ) break;
       if( min.f < fini &&
-          a.f - min.f <= fabs(a.dx)*maxgrad &&
+          a.f - min.f <= fabsx(a.dx)*maxgrad &&
-          (fabs(distmin - range) > tol || maxstep < b.dx) ) break;
+          (fabsx(distmin - range) > tol || maxstep < b.dx) ) break;
 
       mid = .5*(a.dx + b.dx);
-      if( fabs(mid) <= tol2 - .5*(b.dx - a.dx) ) break;
+      if( fabsx(mid) <= tol2 - .5*(b.dx - a.dx) ) break;
 
       r = q = s = 0;
-      if( fabs(end) > tol ) {
+      if( fabsx(end) > tol ) {
         if( first ) {
           creal s1 = w.dx*grad;
           creal s2 = w.f - min.f;
@@ -281,7 +281,7 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
           q = 2*(s2 - s1);
         }
         if( q > 0 ) s = -s;
-        q = fabs(q);
+        q = fabsx(q);
         r = end;
         if( step != b1 || b.dx <= maxstep ) end = step;
       }
@@ -290,15 +290,15 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
       else if( b.dx > maxstep ) step = (step < b.dx) ? -4*a.dx : maxstep;
       else {
         real num = a.dx, den = b.dx;
-        if( fabs(b.dx) <= tol || (w.dx > 0 && fabs(a.dx) > tol) )
+        if( fabsx(b.dx) <= tol || (w.dx > 0 && fabsx(a.dx) > tol) )
           num = b.dx, den = a.dx;
         num /= -den;
-        step = (num < 1) ? .5*den*sqrt(num) : 5/11.*den*(.1 + 1/num);
+        step = (num < 1) ? .5*den*sqrtx(num) : 5/11.*den*(.1 + 1/num);
       }
 
       if( step > 0 ) a1 = a.dx, b1 = step;
       else a1 = step, b1 = b.dx;
-      if( fabs(s) < fabs(.5*q*r) && s > q*a1 && s < q*b1 ) {
+      if( fabsx(s) < fabsx(.5*q*r) && s > q*a1 && s < q*b1 ) {
         step = s/q;
         if( step - a.dx < tol2 || b.dx - step < tol2 )
           step = (mid > 0) ? tol : -tol;
@@ -307,7 +307,7 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree,
     }
 
     first = true;
-    if( fabs(distmin - range) < tol ) {
+    if( fabsx(distmin - range) < tol ) {
       distmin = range;
       if( maxstep > b.dx ) first = false;
     }
@@ -372,7 +372,7 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree,
      or we come close to a border. */
 
   XCopy(y, x);
-  ftest = SUFTOL*(1 + fabs(fx));
+  ftest = SUFTOL*(1 + fabsx(fx));
   delta = RTDELTA/5;
   do {
     delta = Min(5*delta, smax);
@@ -381,7 +381,7 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree,
       y[dim] = x[dim] + delta*p[i];
     }
     fy = Sample(t, y);
-    if( fabs(fy - fx) > ftest ) break;
+    if( fabsx(fy - fx) > ftest ) break;
   } while( delta != smax );
 
   /* Construct a second direction p' orthogonal to p, i.e. p.p' = 0.
@@ -425,7 +425,7 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree,
      or we come close to a border. */
 
   XCopy(z, y);
-  ftest = SUFTOL*(1 + fabs(fy));
+  ftest = SUFTOL*(1 + fabsx(fy));
   delta = RTDELTA/5;
   do {
     delta = Min(5*delta, smax);
@@ -434,7 +434,7 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree,
       z[dim] = y[dim] + delta*p[i];
     }
     fz = Sample(t, z);
-    if( fabs(fz - fy) > ftest ) break;
+    if( fabsx(fz - fy) > ftest ) break;
   } while( delta != smax );
 
   if( fy != fz ) {
@@ -541,12 +541,12 @@ static real FindMinimum(This *t, cBounds *b, real *xmin, real fmin)
     bool resample = false;
     nfree = nfix = 0;
     for( dim = 0; dim < t->ndim; ++dim ) {
-      if( xmin[dim] < b[dim].lower + (1 + fabs(b[dim].lower))*QEPS ) {
+      if( xmin[dim] < b[dim].lower + (1 + fabsx(b[dim].lower))*QEPS ) {
         xmin[dim] = b[dim].lower;
         ifix[nfix++] = dim;
         resample = true;
       }
-      else if( xmin[dim] > b[dim].upper - (1 + fabs(b[dim].upper))*QEPS ) {
+      else if( xmin[dim] > b[dim].upper - (1 + fabsx(b[dim].upper))*QEPS ) {
         xmin[dim] = b[dim].upper;
         ifix[nfix++] = Tag(dim);
         resample = true;
@@ -562,7 +562,7 @@ static real FindMinimum(This *t, cBounds *b, real *xmin, real fmin)
     if( local || Length(nfree, gfree) > GTOL ) break;
 
     ftmp = LocalSearch(t, nfree, ifree, b, xmin, fmin, tmp);
-    if( ftmp > fmin - (1 + fabs(fmin))*RTEPS )
+    if( ftmp > fmin - (1 + fabsx(fmin))*RTEPS )
       goto releasebounds;
     fmin = ftmp;
     XCopy(xmin, tmp);
@@ -586,7 +586,7 @@ static real FindMinimum(This *t, cBounds *b, real *xmin, real fmin)
       minstep = INFTY;
       for( i = 0; i < nfree; ++i ) {
         count dim = Untag(ifree[i]);
-        if( fabs(p[i]) > EPS ) {
+        if( fabsx(p[i]) > EPS ) {
           real step;
           count fix;
           if( p[i] < 0 ) {
@@ -642,11 +642,11 @@ fixbound:
         BFGS(t, nfree, hessian, tmp, gfree, p, low.dx);
         XCopy(gfree, tmp);
 
-        if( fabs(low.dx - minstep) < QEPS*minstep ) goto fixbound;
+        if( fabsx(low.dx - minstep) < QEPS*minstep ) goto fixbound;
 
         fdiff = fini - fmin;
         fini = fmin;
-        if( fdiff > (1 + fabs(fmin))*FTOL ||
+        if( fdiff > (1 + fabsx(fmin))*FTOL ||
             low.dx*plen > (1 + Length(t->ndim, xmin))*FTOL ) continue;
       }
     }

src/external/libCuba/src/divonne/Integrate.c

@@ -5,7 +5,7 @@
 		then do a main integration over all regions
 		this file is part of Divonne
 		checkpointing by B. Chokoufe
-		last modified 5 Aug 13 th
+		last modified 13 Mar 15 th
 */
 
 
@@ -29,8 +29,8 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   Totals *tot, *Tot = state->totals + t->ncomp;
   Bounds *b, *B;
   Result *res;
-  count comp, err, iregion;
+  count comp, iregion;
-  number nwant;
+  number nwant, err;
   real nneed;
   ML_ONLY(number neff;)
   int fail;
@@ -38,6 +38,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   if( VERBOSE > 1 ) {
     sprintf(out, "Divonne input parameters:\n"
       "  ndim " COUNT "\n  ncomp " COUNT "\n"
+      ML_NOT("  nvec " NUMBER "\n")
       "  epsrel " REAL "\n  epsabs " REAL "\n"
       "  flags %d\n  seed %d\n"
       "  mineval " NUMBER "\n  maxeval " NUMBER "\n"
@@ -46,11 +47,12 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       "  ngiven " NUMBER "\n  nextra " NUMBER "\n"
       "  statefile \"%s\"",
       t->ndim, t->ncomp,
-      t->epsrel, t->epsabs,
+      ML_NOT(t->nvec,)
+      SHOW(t->epsrel), SHOW(t->epsabs),
       t->flags, t->seed,
       t->mineval, t->maxeval,
       t->key1, t->key2, t->key3, t->maxpass,
-      t->border.lower, t->maxchisq, t->mindeviation,
+      SHOW(t->border.lower), SHOW(t->maxchisq), SHOW(t->mindeviation),
       t->ngiven, t->nextra,
       t->statefile);
     Print(out);
@@ -132,7 +134,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   /* Step 1: partition the integration region */
 
   if( t->phase == 1 ) {
-    if( VERBOSE ) Print("Partitioning phase:");
+    if( VERBOSE ) Print("\nPartitioning phase:");
 
     if( ini ) Iterate(t, 0, INIDEPTH, 0, NULL);
 
@@ -164,7 +166,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
         integral[comp] = tot->avg;
         valid += tot->avg == tot->avg;
         if( tot->spreadsq > maxtot->spreadsq ) maxtot = tot;
-        tot->spread = sqrt(tot->spreadsq);
+        tot->spread = sqrtx(tot->spreadsq);
         error[comp] = tot->spread/t->samples[0].neff;
       }
 
@@ -197,7 +199,7 @@ if( StateWriteTest(t) ) { \
         for( comp = 0; comp < t->ncomp; ++comp )
           oe += sprintf(oe, "\n[" COUNT "] "
             REAL " +- " REAL,
-            comp + 1, integral[comp], error[comp]);
+            comp + 1, SHOW(integral[comp]), SHOW(error[comp]));
         Print(out);
       }
 
@@ -230,7 +232,7 @@ if( StateWriteTest(t) ) { \
     tot->maxerrsq = Sq(maxerr);
     tot->mindevsq = tot->maxerrsq*Sq(t->mindeviation);
   }
-  nwant = (number)Min(ceil(nneed), MARKMASK/40.);
+  nwant = (number)Min(ceil(nneed), NWANTMAX/40.);
 
   err = SamplesLookup(t, &t->samples[1], t->key2, nwant,
     (t->maxeval - t->neval)/t->nregions + 1, t->samples[0].n + 1);
@@ -308,11 +310,11 @@ refine:
                 Iterate(t, state->iregion, POSTDEPTH, 1, state->totals);
 
                 if( can_adjust ) {
-                  cnumber nnew = (tot->spreadsq/Sq(MARKMASK) > tot->maxerrsq) ?
+                  cnumber nnew = (tot->spreadsq/Sq(NWANTMAX) > tot->maxerrsq) ?
-                    MARKMASK :
+                    NWANTMAX :
-                    (number)ceil(sqrt(tot->spreadsq/tot->maxerrsq));
+                    (number)ceil(sqrtx(tot->spreadsq/tot->maxerrsq));
                   if( nnew > nwant + nwant/64 ) {
-                    ccount err = SamplesLookup(t, &t->samples[1], t->key2, nnew,
+                    cnumber err = SamplesLookup(t, &t->samples[1], t->key2, nnew,
                       (t->maxeval - t->neval)/t->nregions + 1, t->samples[1].n);
                     fail += Unmark(err)*t->nregions;
                     nwant = nnew;
@@ -404,14 +406,14 @@ refine:
         if( chisq > EPS ) chisq /= Max(chiden, NOTZERO);
 
         if( VERBOSE > 2 ) {
-#define Out2(f, r) (r)->avg, res->spread/t->samples[f].neff, (r)->err
+#define Out2(f, r) SHOW((r)->avg), SHOW(res->spread/t->samples[f].neff), SHOW((r)->err)
 #define Out(f) Out2(f, &tot->phase[f])
           oe += sprintf(oe, "\n[" COUNT "] "
             REAL " +- " REAL "(" REAL ")\n    "
             REAL " +- " REAL "(" REAL ")", ++comp, Out(0), Out(1));
           if( todo == 3 ) oe += sprintf(oe, "\n    "
             REAL " +- " REAL "(" REAL ")", Out2(2, res));
-          oe += sprintf(oe, "  \tchisq " REAL, chisq);
+          oe += sprintf(oe, "  \tchisq " REAL, SHOW(chisq));
         }
 
         tot->integral += avg;
@@ -419,7 +421,7 @@ refine:
         tot->chisq += chisq;
 
         res->avg = avg;
-        res->spread = sqrt(sigsq);
+        res->spread = sqrtx(sigsq);
         res->chisq = chisq;
       }
 
@@ -433,7 +435,7 @@ refine:
 
     for( tot = state->totals, comp = 0; tot < Tot; ++tot, ++comp ) {
       integral[comp] = tot->integral;
-      error[comp] = sqrt(tot->sigsq);
+      error[comp] = sqrtx(tot->sigsq);
       prob[comp] = ChiSquare(tot->chisq, df);
     }
 
@@ -442,7 +444,8 @@ refine:
       for( tot = state->totals, comp = 0; tot < Tot; ++tot, ++comp )
         oe += sprintf(oe, "\n[" COUNT "] "
           REAL " +- " REAL "  \tchisq " REAL " (" COUNT " df)",
-          comp + 1, integral[comp], error[comp], tot->chisq, df);
+          comp + 1, SHOW(integral[comp]), SHOW(error[comp]),
+          SHOW(tot->chisq), df);
       Print(out);
     }
 
@@ -468,12 +471,12 @@ refine:
 
       MLPutFunction(stdlink, "Cuba`Divonne`region", 4);
 
-      MLPutRealList(stdlink, bounds, 2*t->ndim);
+      MLPutRealxList(stdlink, bounds, 2*t->ndim);
 
       MLPutFunction(stdlink, "List", t->ncomp);
       for( Res = (res = RegionResult(region)) + t->ncomp; res < Res; ++res ) {
         real r[] = {res->avg, res->spread/neff, res->chisq};
-        MLPutRealList(stdlink, r, Elements(r));
+        MLPutRealxList(stdlink, r, Elements(r));
       }
 
       MLPutInteger(stdlink, region->depth + 1);  /* misused for df */
@@ -482,8 +485,7 @@ refine:
 #endif
 
 abort:
-  WaitCores(t);
+  FORK_ONLY(FrameFree(t, Master);)
-  FORK_ONLY(FrameFree(t, ShmRm(t));)
 
   RuleFree(t);
   SamplesFree(&t->samples[2]);

src/external/libCuba/src/divonne/Iterate.c

@@ -2,7 +2,7 @@
 	Iterate.c
 		recursion over regions
 		this file is part of Divonne
-		last modified 2 Aug 13 th
+		last modified 12 Mar 15 th
 */
 
 
@@ -85,10 +85,10 @@ FORK_ONLY(more:)
 
       norm = 1./nsplit--;
       for( res = RegionResult(parent), c = corr; c < C; ++res, ++c ) {
-        creal diff = fabs(res->avg - c->avg)*norm;
+        creal diff = fabsx(res->avg - c->avg)*norm;
         c->avg = diff*norm*nsplit;
         c->err = (c->err == 0) ? 1 : 1 + diff/c->err;
-        c->spread = (c->spread == 0) ? 1 : 1 + diff/sqrt(c->spread);
+        c->spread = (c->spread == 0) ? 1 : 1 + diff/sqrtx(c->spread);
       }
 
       for( td = todo; td < tdmax; ++td )

src/external/libCuba/src/divonne/Rule.c

@@ -4,14 +4,10 @@
 		code lifted with minor modifications from DCUHRE
 		by J. Berntsen, T. Espelid, and A. Genz
 		this file is part of Divonne
-		last modified 26 Jul 13 th
+		last modified 7 May 15 th
 */
 
 
-#define NextSet(p) p = (Set *)((char *)p + setsize)
-
-/*********************************************************************/
-
 static void Rule13Alloc(This *t)
 {
   static creal w[][nrules] = {
@@ -153,7 +149,7 @@ static void Rule13Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -298,7 +294,7 @@ static void Rule11Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -309,32 +305,55 @@ static void Rule11Alloc(This *t)
 static void Rule9Alloc(This *t)
 {
   static creal w[] = {
-    -.0023611709677855117884,   .11415390023857325268,
+    RC(-.002361170967785511788400941242259231309691),
-    -.63833920076702389094,     .74849988504685208004,
+    RC(.1141539002385732526821323741697655347686),
-    -.0014324017033399125142,   .057471507864489725949,
+    RC(-.6383392007670238909386026193674701393074),
-    -.14225104571434243234,    -.062875028738286979989,
+    RC(.7484998850468520800423030047583803945205),
-     .254591133248959089,     -1.207328566678236261,
+    RC(-.001432401703339912514196154599769007103671),
-     .89567365764160676508,    -.36479356986049146661,
+    RC(.05747150786448972594860897296200006759892),
-     .0035417564516782676826,  -.072609367395893679605,
+    RC(-.1422510457143424323449521620935950679394),
-     .10557491625218991012,     .0021486025550098687713,
+    RC(-.06287502873828697998942424881040490136987),
-    -.032268563892953949998,    .010636783990231217481,
+    RC(.2545911332489590890011611142429070613156),
-     .014689102496143490175,    .51134708346467591431,
+    RC(-1.207328566678236261002219995185143356737),
-     .45976448120806344646,     .18239678493024573331,
+    RC(.8956736576416067650809467826488567200939),
-    -.04508628929435784076,     .21415883524352793401,
+    RC(-.3647935698604914666100134551377381205297),
-    -.027351546526545644722,    .054941067048711234101,
+    RC(.003541756451678267682601411863388846964536),
-     .11937596202570775297,     .65089519391920250593,
+    RC(-.07260936739589367960492815865074633743652),
-     .14744939829434460168,     .057693384490973483573,
+    RC(.1055749162521899101218622863269817454540),
-     .034999626602143583822,  -1.3868627719278281436,
+    RC(.002148602555009868771294231899653510655506),
-    -.2386668732575008879,      .015532417276607053264,
+    RC(-.03226856389295394999786630399875134318006),
-     .0035328099607090870236,   .09231719987444221619,
+    RC(.01063678399023121748083624225818915724455),
-     .02254314464717892038,     .013675773263272822361,
+    RC(.01468910249614349017540783437728097691502),
-    -.32544759695960125297,     .0017708782258391338413,
+    RC(.5113470834646759143109387357149329909126),
-     .0010743012775049343856,   .25150011495314791996 };
+    RC(.4597644812080634464633352781605214342691),
+    RC(.1823967849302457333050067275688690602649),
+    RC(-.04508628929435784075980562738240804429658),
+    RC(.2141588352435279340097929526588394300172),
+    RC(-.02735154652654564472203690086290223507436),
+    RC(.05494106704871123410060080562462135546101),
+    RC(.1193759620257077529708962121565290178730),
+    RC(.6508951939192025059314756320878023215278),
+    RC(.1474493982943446016775696826942585013243),
+    RC(.05769338449097348357291272840392627722165),
+    RC(.03499962660214358382244159694487155861542),
+    RC(-1.386862771927828143599782668709014266770),
+    RC(-.2386668732575008878964134721962088068396),
+    RC(.01553241727660705326386197156586357005224),
+    RC(.003532809960709087023561817517751309380604),
+    RC(.09231719987444221619017126187763868745587),
+    RC(.02254314464717892037990281369120402214829),
+    RC(.01367577326327282236101845043145111753718),
+    RC(-.3254475969596012529657378160439011607639),
+    RC(.001770878225839133841300705931694423482268),
+    RC(.001074301277504934385647115949826755327753),
+    RC(.2515001149531479199576969952416196054795) };
 
   static creal g[] = {
-     .47795365790226950619,     .20302858736911986780,
+    RC(.4779536579022695061928604197171830064732),
-     .44762735462617812882,     .125,
+    RC(.2030285873691198677998034402373279133258),
-     .34303789878087814570 };
+    RC(.4476273546261781288207704806530998539285),
+    RC(.125),
+    RC(.3430378987808781457001426145164678603407) };
 
   enum { nsets = 9 };
 
@@ -440,7 +459,7 @@ static void Rule9Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -451,21 +470,33 @@ static void Rule9Alloc(This *t)
 static void Rule7Alloc(This *t)
 {
   static creal w[] = {
-     .019417866674748388428,   -.40385257701150182546,
+    RC(.01941786667474838842844534313920462333850),
-     .64485668767465982223,     .01177982690775806141,
+    RC(-.4038525770115018254611834753723880293161),
-    -.18041318740733609012,    -.088785828081335044443,
+    RC(.6448566876746598222277360730193089551024),
-     .056328645808285941374,   -.0097089333373741942142,
+    RC(.01177982690775806141012214458820955067854),
-    -.99129176779582358138,    -.17757165616267008889,
+    RC(-.1804131874073360901182293138710989490609),
-     .12359398032043233572,     .074978148702033690681,
+    RC(-.08878582808133504444306598174517276122439),
-     .55489147051423559776,     .088041241522692771226,
+    RC(.05632864580828594137378124255408286479947),
-     .021118358455513385083,   -.0099302203239653333087,
+    RC(-.009708933337374194214222671569602311669249),
-    -.064100053285010904179,    .030381729038221007659,
+    RC(-.9912917677958235813775106862002319060386),
-     .0058899134538790307051,  -.0048544666686870971071,
+    RC(-.1775716561626700888861319634903455224488),
-     .35514331232534017777 };
+    RC(.1235939803204323357183625846672135876752),
+    RC(.07497814870203369068087999555157339703666),
+    RC(.5548914705142355977605994477355651401434),
+    RC(.08804124152269277122645182458858273865209),
+    RC(.02111835845551338508329573367808085283304),
+    RC(-.009930220323965333308685820460105538586058),
+    RC(-.06410005328501090417895544042025034295870),
+    RC(.03038172903822100765927778829870429682489),
+    RC(.005889913453879030705061072294104775339268),
+    RC(-.004854466668687097107111335784801155834624),
+    RC(.3551433123253401777722639269806910448976) };
 
   static creal g[] = {
-     .47795365790226950619,     .20302858736911986780,
+    RC(.4779536579022695061928604197171830064732),
-     .375,                      .34303789878087814570 };
+    RC(.2030285873691198677998034402373279133258),
+    RC(.375),
+    RC(.3430378987808781457001426145164678603407) };
 
   enum { nsets = 6 };
 
@@ -541,7 +572,7 @@ static void Rule7Alloc(This *t)
         -s->weight[r + 1]/s->weight[r];
       real sum = 0;
       for( x = first; x <= last; NextSet(x) )
-        sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]);
+        sum += x->n*fabsx(x->weight[r + 1] + scale*x->weight[r]);
       s->scale[r] = scale;
       s->norm[r] = 1/sum;
     }
@@ -666,7 +697,7 @@ static void SampleRule(This *t, ccount iregion)
       real maxerr = 0;
       for( s = first; s <= last; NextSet(s) )
         maxerr = Max(maxerr,
-          fabs(sum[rul + 1] + s->scale[rul]*sum[rul])*s->norm[rul]);
+          fabsx(sum[rul + 1] + s->scale[rul]*sum[rul])*s->norm[rul]);
       sum[rul] = maxerr;
     }
 

src/external/libCuba/src/divonne/Sample.c

@@ -2,14 +2,16 @@
 	Sample.c
 		most of what is related to sampling
 		this file is part of Divonne
-		last modified 30 Aug 13 th
+		last modified 12 Mar 15 th
 */
 
 
-#define MARKMASK 0xfffffff
+#define MARKMASK NUMBER_MAX
 #define Marked(x) ((x) & ~MARKMASK)
 #define Unmark(x) ((x) & MARKMASK)
 
+#define NWANTMAX NUMBER_MAX
+
 #define EXTRAPOLATE_EPS (.25*t->border.lower)
 /*#define EXTRAPOLATE_EPS 0x1p-26*/
 
@@ -96,7 +98,7 @@ static void SampleKorobov(This *t, ccount iregion)
   }
 
   if( dist > 0 ) {
-    dist = sqrt(dist)/EXTRAPOLATE_EPS;
+    dist = sqrtx(dist)/EXTRAPOLATE_EPS;
     for( dim = 0; dim < t->ndim; ++dim ) {
       real x2 = x[dim], dx = x2 - t->border.upper;
       if( dx > 0 ) {
@@ -148,7 +150,7 @@ static void SampleKorobov(This *t, ccount iregion)
          1..39 = multiplicator,        Korobov numbers,
        40..inf = absolute # of points, Korobov numbers.  */
 
-static count SamplesLookup(This *t, Samples *samples, cint key,
+static number SamplesLookup(This *t, Samples *samples, cint key,
   cnumber nwant, cnumber nmax, number nmin)
 {
   number n;
@@ -191,7 +193,8 @@ static count SamplesLookup(This *t, Samples *samples, cint key,
 static void SamplesAlloc(cThis *t, Samples *samples)
 {
 #define FIRST -INT_MAX
-#define MarkLast(x) (x | Marked(INT_MAX))
+#define MarkLast(x) ((x) | 0x40000000)
+#define UnmarkLast(x) ((x) & 0x3fffffff)
 
 #include "KorobovCoeff.c"
 
@@ -205,12 +208,12 @@ static void SamplesAlloc(cThis *t, Samples *samples)
 
     while( i = IMin(IDim(i), max),
            n > (p = prime[i + 1]) || n <= prime[i] ) {
-      cint d = (n - Unmark(p)) >> ++shift;
+      cint d = (n - UnmarkLast(p)) >> ++shift;
       i += Min1(d);
     }
 
     samples->coeff = coeff[i][t->ndim - KOROBOV_MINDIM];
-    samples->neff = p = Unmark(p);
+    samples->neff = p = UnmarkLast(p);
     samples->n = p/2 + 1;
   }
 
@@ -240,7 +243,7 @@ static real Sample(This *t, creal *x0)
   }
 
   if( dist > 0 ) {
-    dist = sqrt(dist)/EXTRAPOLATE_EPS;
+    dist = sqrtx(dist)/EXTRAPOLATE_EPS;
     for( dim = 0; dim < t->ndim; ++dim ) {
       real x2 = xtmp[dim], dx, b;
       if( (dx = x2 - (b = t->border.lower)) < 0 ||

src/external/libCuba/src/divonne/Split.c

@@ -2,7 +2,7 @@
 	Split.c
 		determine optimal cuts for splitting a region
 		this file is part of Divonne
-		last modified 31 Aug 13 th
+		last modified 12 Mar 15 th
 */
 
 
@@ -31,7 +31,7 @@ typedef struct {
 
 static inline real Div(creal a, creal b)
 {
-  return (b != 0 /*&& fabs(a) > SMALL*fabs(b)*/) ? a/b : a;
+  return (b != 0 /*&& fabsx(a) > SMALL*fabsx(b)*/) ? a/b : a;
 }
 
 /*********************************************************************/
@@ -58,7 +58,7 @@ static void SomeCut(This *t, Cut *cut, Bounds *b)
     yupper = Sample(t, xmid);
     xmid[dim] = x;
 
-    dev = fabs(ymid - .5*(ylower + yupper));
+    dev = fabsx(ymid - .5*(ylower + yupper));
     if( dev >= maxdev ) {
       maxdev = dev;
       maxdim = dim;
@@ -178,7 +178,7 @@ static count FindCuts(This *t, Cut *cut, Bounds *bounds, creal vol,
 
     for( icut = 0; icut < ncuts; ++icut ) {
       Cut *c = &cut[icut];
-      creal diff = fabs(fmajor - c->f);
+      creal diff = fabsx(fmajor - c->f);
       if( diff <= mindiff ) {
         mindiff = diff;
         mincut = c;
@@ -230,18 +230,18 @@ repeat:
     if( lhssqnew <= lhssq ) {
       real fmax;
 
-      if( fabs(gammanew - gamma) < GAMMATOL*gamma ) break;
+      if( fabsx(gammanew - gamma) < GAMMATOL*gamma ) break;
       gamma = gammanew;
 
-      fmax = fabs(fgamma);
+      fmax = fabsx(fgamma);
       for( icut = 0; icut < ncuts; ++icut ) {
         Cut *c = &cut[icut];
         creal dfmin = SINGTOL*c->df;
         creal sol = c->sol/div;
         real df = c->f - c->fold;
-        df = (fabs(df) > SMALL*fabs(sol)) ? df/sol : 1;
+        df = (fabsx(df) > SMALL*fabsx(sol)) ? df/sol : 1;
-        c->df = (fabs(df) < fabs(dfmin)) ? dfmin : df;
+        c->df = (fabsx(df) < fabsx(dfmin)) ? dfmin : df;
-        fmax = Max(fmax, fabs(c->f));
+        fmax = Max(fmax, fabsx(c->f));
         c->fold = c->f;
       }
 
@@ -276,7 +276,7 @@ static void Split(This *t, ccount iregion)
   t->selectedcomp = region->cutcomp;
   t->neval_cut -= t->neval;
   ncuts = FindCuts(t, cut, region->bounds, region->vol,
-    (real *)RegionResult(region) + region->xmajor, region->fmajor,
+    RegionMinMax(region) + region->xmajor, region->fmajor,
     region->fmajor - region->fminor);
   t->neval_cut += t->neval;
 

src/external/libCuba/src/divonne/decl.h

@@ -2,7 +2,7 @@
 	decl.h
 		Type declarations
 		this file is part of Divonne
-		last modified 26 Jul 13 th
+		last modified 9 Oct 14 th
 */
 
 
@@ -44,6 +44,8 @@ typedef struct {
 
 #define SetSize (sizeof(Set) + t->ndim*sizeof(real))
 
+#define NextSet(p) p = (Set *)((char *)p + setsize)
+
 typedef struct {
   Set *first, *last;
   real errcoeff[3];
@@ -71,12 +73,11 @@ typedef const Errors cErrors;
 typedef struct {
   real avg, err, spread, chisq;
   real fmin, fmax;
-  real xminmax[];
 } Result;
 
 typedef const Result cResult;
 
-#define ResultSize (sizeof(Result) + t->ndim*2*sizeof(real))
+#define MinMaxSize (t->ncomp*t->ndim*2*sizeof(real))
 
 typedef struct region {
   int depth, next;
@@ -85,29 +86,34 @@ typedef struct region {
   Bounds bounds[];
 } Region;
 
-#define RegionSize (sizeof(Region) + t->ndim*sizeof(Bounds) + t->ncomp*ResultSize)
+#define RegionSize (sizeof(Region) + t->ndim*sizeof(Bounds) + t->ncomp*sizeof(Result) + MinMaxSize)
 
 #define RegionResult(r) ((Result *)(r->bounds + t->ndim))
 
+#define RegionMinMax(r) ((real *)(RegionResult(r) + t->ncomp))
+
 #define RegionPtr(n) ((Region *)((char *)t->region + (n)*regionsize))
 
 
-typedef int (*Integrand)(ccount *, creal *, ccount *, real *, void *, cint *);
+typedef int (*Integrand)(ccount *, creal *, ccount *, real *,
+  void *, cnumber *, cint *, cint *);
 
-typedef void (*PeakFinder)(ccount *, cBounds *, number *, real *);
+typedef void (*PeakFinder)(ccount *, cBounds *, number *, real *, void *);
 
 typedef struct _this {
   count ndim, ncomp;
 #ifndef MLVERSION
   Integrand integrand;
   void *userdata;
-  PeakFinder peakfinder;
+  number nvec;
 #ifdef HAVE_FORK
-  int ncores, running, *child;
+  SHM_ONLY(int shmid;)
+  Spin *spin;
   real *frame;
   number nframe;
-  SHM_ONLY(int shmid;)
+  int running;
 #endif
+  PeakFinder peakfinder;
 #endif
   real epsrel, epsabs;
   int flags, seed;

src/external/libCuba/src/suave/Fluct.c

@@ -2,35 +2,43 @@
 	Fluct.c
 		compute the fluctuation in the left and right half
 		this file is part of Suave
-		last modified 29 Jul 13 th
+		last modified 14 Mar 15 th
 */
 
 
-#if defined(HAVE_LONG_DOUBLE) && defined(HAVE_POWL)
+#if defined(HAVE_LONG_DOUBLE) && defined(HAVE_POWL) && REALSIZE <= 10
 
-typedef long double realx;
+typedef long double realL;
-#define XDBL_MAX_EXP LDBL_MAX_EXP
+#define REALL_MAX_EXP LDBL_MAX_EXP
-#define XDBL_MAX LDBL_MAX
+#define REALL_MAX LDBL_MAX
-#define powx powl
+#define powL powl
-#define ldexpx ldexpl
+#define ldexpL ldexpl
 
 #else
 
-typedef double realx;
+typedef real realL;
-#define XDBL_MAX_EXP DBL_MAX_EXP
+#define REALL_MAX_EXP REAL_MAX_EXP
-#define XDBL_MAX DBL_MAX
+#define REALL_MAX REAL_MAX
-#define powx pow
+#define powL powx
-#define ldexpx ldexp
+#define ldexpL ldexpx
 
 #endif
 
-typedef const realx crealx;
+typedef const realL crealL;
 
 typedef struct {
-  realx fluct;
+  realL fluct;
   number n;
 } Var;
 
+static inline realL MinL(crealL a, crealL b) {
+  return (a < b) ? a : b;
+}
+
+static inline realL MaxL(crealL a, crealL b) {
+  return (a > b) ? a : b;
+}
+
 /*********************************************************************/
 
 static void Fluct(cThis *t, Var *var,
@@ -41,27 +49,27 @@ static void Fluct(cThis *t, Var *var,
   count nvar = 2*t->ndim;
   creal norm = 1/(err*Max(fabs(avg), err));
   creal flat = 2/3./t->flatness;
-  crealx max = ldexpx(1., (int)((XDBL_MAX_EXP - 2)/t->flatness));
+  crealL max = ldexpL(1., (int)((REALL_MAX_EXP - 2)/t->flatness));
 
   Clear(var, nvar);
 
   while( n-- ) {
     count dim;
-    crealx arg = 1 + fabs(*w++)*Sq(*f - avg)*norm;
+    crealL arg = 1 + fabs(*w++)*Sq(*f - avg)*norm;
-    crealx ft = powx(arg < max ? arg : max, t->flatness);
+    crealL ft = powL(MinL(arg, max), t->flatness);
 
     f += t->ncomp;
 
     for( dim = 0; dim < t->ndim; ++dim ) {
       Var *v = &var[2*dim + (*x++ >= .5*(b[dim].lower + b[dim].upper))];
-      crealx f = v->fluct + ft;
+      crealL f = v->fluct + ft;
-      v->fluct = (f > XDBL_MAX/2) ? XDBL_MAX/2 : f;
+      v->fluct = MaxL(f, REALL_MAX/2);
       ++v->n;
     }
   }
 
   while( nvar-- ) {
-    var->fluct = powx(var->fluct, flat);
+    var->fluct = powL(var->fluct, flat);
     ++var;
   }
 }

src/external/libCuba/src/suave/Grid.c

@@ -2,7 +2,7 @@
 	Grid.c
 		utility functions for the Vegas grid
 		this file is part of Suave
-		last modified 7 Aug 13 th
+		last modified 12 Mar 15 th
 */
 
 
@@ -33,7 +33,7 @@ static void RefineGrid(cThis *t, Grid grid, Grid margsum)
     real impfun = 0;
     if( margsum[bin] > 0 ) {
       creal r = margsum[bin]*norm;
-      avgperbin += impfun = pow((r - 1)/log(r), 1.5);
+      avgperbin += impfun = powx((r - 1)/log(r), 1.5);
     }
     imp[bin] = impfun;
   }

src/external/libCuba/src/suave/Integrate.c

@@ -3,7 +3,7 @@
 		integrate over the unit hypercube
 		this file is part of Suave
 		checkpointing by B. Chokoufe
-		last modified 5 Aug 13 th
+		last modified 13 Mar 15 th
 */
 
 
@@ -26,22 +26,27 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   Result *tot, *Tot = state->totals + t->ncomp;
   Result *res, *resL, *resR;
   Bounds *b, *B;
+  cnumber minsamples = IMax(t->nmin, MINSAMPLES);
   count dim, comp;
   int fail;
 
   if( VERBOSE > 1 ) {
     sprintf(out, "Suave input parameters:\n"
       "  ndim " COUNT "\n  ncomp " COUNT "\n"
+      ML_NOT("  nvec " NUMBER "\n")
       "  epsrel " REAL "\n  epsabs " REAL "\n"
       "  flags %d\n  seed %d\n"
       "  mineval " NUMBER "\n  maxeval " NUMBER "\n"
-      "  nnew " NUMBER "\n  flatness " REAL "\n"
+      "  nnew " NUMBER "\n  nmin " NUMBER "\n"
-      "  statefile \"%s\"\n",
+      "  flatness " REAL "\n"
+      "  statefile \"%s\"",
       t->ndim, t->ncomp,
-      t->epsrel, t->epsabs,
+      ML_NOT(t->nvec,)
+      SHOW(t->epsrel), SHOW(t->epsabs),
       t->flags, t->seed,
       t->mineval, t->maxeval,
-      t->nnew, t->flatness,
+      t->nnew, t->nmin,
+      SHOW(t->flatness),
       t->statefile);
     Print(out);
   }
@@ -49,7 +54,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   if( BadComponent(t) ) return -2;
   if( BadDimension(t) ) return -1;
 
-  ShmAlloc(t, ShmRm(t));
+  ShmAlloc(t, Master);
   ForkCores(t);
 
   if( (fail = setjmp(t->abort)) ) goto abort;
@@ -132,7 +137,8 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       for( tot = state->totals, comp = 0; tot < Tot; ++tot )
         oe += sprintf(oe, "\n[" COUNT "] "
           REAL " +- " REAL "  \tchisq " REAL " (" COUNT " df)",
-          ++comp, tot->avg, tot->err, tot->chisq, state->df);
+          ++comp, SHOW(tot->avg), SHOW(tot->err),
+          SHOW(tot->chisq), state->df);
       Print(out);
     }
 
@@ -173,7 +179,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       region->result[maxcomp].avg, Max(maxerr, t->epsabs));
 
     bias = (t->epsrel < 1e-50) ? 2 :
-      Max(pow(2., -(real)region->div/t->ndim)/t->epsrel, 2.);
+      Max(powx(2., -(real)region->div/t->ndim)/t->epsrel, 2.);
     minfluct = INFTY;
     bisectdim = 0;
     for( dim = 0; dim < t->ndim; ++dim ) {
@@ -190,9 +196,9 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
     minfluct = vLR[0].fluct + vLR[1].fluct;
     nnewL = IMax(
       (minfluct == 0) ? t->nnew/2 : (count)(vLR[0].fluct/minfluct*t->nnew),
-      MINSAMPLES );
+      minsamples );
     nL = vLR[0].n + nnewL;
-    nnewR = IMax(t->nnew - nnewL, MINSAMPLES);
+    nnewR = IMax(t->nnew - nnewL, minsamples);
     nR = vLR[1].n + nnewR;
 
     regionL = RegionAlloc(t, nL, nnewL);
@@ -211,7 +217,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
     while( n-- ) {
       cbool final = (*w < 0);
       if( x[bisectdim] < mid ) {
-        if( final && wR > RegionW(regionR) ) wR[-1] = -fabs(wR[-1]);
+        if( final && wR > RegionW(regionR) ) wR[-1] = -fabsx(wR[-1]);
         *wL++ = *w++;
         XCopy(xL, x);
         xL += t->ndim;
@@ -219,7 +225,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
         fL += t->ncomp;
       }
       else {
-        if( final && wL > RegionW(regionL) ) wL[-1] = -fabs(wL[-1]);
+        if( final && wL > RegionW(regionL) ) wL[-1] = -fabsx(wL[-1]);
         *wR++ = *w++;
         XCopy(xR, x);
         xR += t->ndim;
@@ -259,15 +265,15 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       diff = Sq(.25*diff);
       sigsq = resL->sigsq + resR->sigsq;
       if( sigsq > 0 ) {
-        creal c = Sq(1 + sqrt(diff/sigsq));
+        creal c = Sq(1 + sqrtx(diff/sigsq));
         resL->sigsq *= c;
         resR->sigsq *= c;
       }
-      resL->err = sqrt(resL->sigsq += diff);
+      resL->err = sqrtx(resL->sigsq += diff);
-      resR->err = sqrt(resR->sigsq += diff);
+      resR->err = sqrtx(resR->sigsq += diff);
 
       tot->sigsq += resL->sigsq + resR->sigsq - res->sigsq;
-      tot->err = sqrt(tot->sigsq);
+      tot->err = sqrtx(tot->sigsq);
 
       tot->chisq += resL->chisq + resR->chisq - res->chisq;
     }
@@ -319,12 +325,12 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
 
       MLPutFunction(stdlink, "Cuba`Suave`region", 3);
 
-      MLPutRealList(stdlink, bounds, 2*t->ndim);
+      MLPutRealxList(stdlink, bounds, 2*t->ndim);
 
       MLPutFunction(stdlink, "List", t->ncomp);
       for( Res = (res = region->result) + t->ncomp; res < Res; ++res ) {
         real r[] = {res->avg, res->err, res->chisq};
-        MLPutRealList(stdlink, r, Elements(r));
+        MLPutRealxList(stdlink, r, Elements(r));
       }
 
       MLPutInteger(stdlink, region->df);
@@ -338,8 +344,7 @@ abort:
     anchor = anchor->next;
     free(region);
   }
-  WaitCores(t);
+  ShmFree(t, Master);
-  ShmFree(t);
 
   StateRemove(t);
 

src/external/libCuba/src/suave/Sample.c

@@ -2,7 +2,7 @@
 	Sample.c
 		the sampling step of Suave
 		this file is part of Suave
-		last modified 30 Jul 13 th
+		last modified 13 Mar 15 th
 */
 
 
@@ -71,7 +71,7 @@ static void Sample(This *t, cnumber nnew, Region *region,
 
   while( w < lastw ) {
     cbool final = (*w < 0);
-    creal weight = fabs(*w++);
+    creal weight = fabsx(*w++);
     ++n;
 
     for( c = cumul, comp = 0; c < C; ++c ) {
@@ -80,20 +80,20 @@ static void Sample(This *t, cnumber nnew, Region *region,
       c->sqsum += Sq(wfun);
 
       if( final ) {
-        if( n > 1 ) {
+        if( n >= t->nmin ) {
           real w = Weight(c->sum, c->sqsum, n);
           c->weightsum += c->weight = w;
           c->avgsum += c->avg = w*c->sum;
 
           if( VERBOSE > 2 ) {
-            creal sig = sqrt(1/w);
+            creal sig = sqrtx(1/w);
             ss[comp] += (df == 0) ?
               sprintf(ss[comp], "\n[" COUNT "] "
                 REAL " +- " REAL " (" NUMBER ")", comp + 1,
-                c->sum, sig, n) :
+                SHOW(c->sum), SHOW(sig), n) :
               sprintf(ss[comp], "\n    "
                 REAL " +- " REAL " (" NUMBER ")",
-                c->sum, sig, n);
+                SHOW(c->sum), SHOW(sig), n);
           }
 
           if( df == 0 ) c->guess = c->sum;
@@ -107,7 +107,7 @@ static void Sample(This *t, cnumber nnew, Region *region,
       }
     }
 
-    if( final ) ++df, n = 0;
+    if( final ) df -= NegQ(t->nmin - n - 1), n = 0;
   }
 
   region->df = --df;
@@ -124,7 +124,7 @@ static void Sample(This *t, cnumber nnew, Region *region,
       res->sigsq = sigsq;
       res->avg = avg;
     }
-    res->err = sqrt(res->sigsq);
+    res->err = sqrtx(res->sigsq);
 
     res->chisq = (sigsq < .9*NOTZERO) ? 0 : c->chisqsum - avg*c->chisum;
       /* This catches the special case where the integrand is constant
@@ -153,7 +153,7 @@ static void Sample(This *t, cnumber nnew, Region *region,
     for( comp = 0, res = region->result;
          comp < t->ncomp; ++comp, ++res ) {
       p += sprintf(p, "%s  \tchisq " REAL " (" COUNT " df)",
-        p0, res->chisq, df);
+        p0, SHOW(res->chisq), df);
       p0 += chars;
     }
 

src/external/libCuba/src/suave/Suave.c

@@ -1,8 +1,8 @@
 /*
 	Suave.c
-		Subregion-adaptive Vegas Monte-Carlo integration
+		Subregion-adaptive Vegas Monte Carlo integration
 		by Thomas Hahn
-		last modified 17 Sep 13 th
+		last modified 28 Nov 14 th
 */
 
 
@@ -15,64 +15,80 @@
 /*********************************************************************/
 
 Extern void EXPORT(Suave)(ccount ndim, ccount ncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber nvec,
   creal epsrel, creal epsabs,
   cint flags, cint seed,
   cnumber mineval, cnumber maxeval,
-  cnumber nnew, creal flatness,
+  cnumber nnew, cnumber nmin, creal flatness,
-  cchar *statefile,
+  cchar *statefile, Spin **pspin,
   count *pnregions, number *pneval, int *pfail,
   real *integral, real *error, real *prob)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = ndim;
   t.ncomp = ncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = nvec;
   t.epsrel = epsrel;
   t.epsabs = epsabs;
-  t.flags = flags;
+  t.flags = MaxVerbose(flags);
   t.seed = seed;
   t.mineval = mineval;
   t.maxeval = maxeval;
   t.nnew = nnew;
+  t.nmin = IMax(nmin, 2);
   t.flatness = flatness;
   t.statefile = statefile;
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
 
   *pfail = Integrate(&t, integral, error, prob);
   *pnregions = t.nregions;
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
 
 /*********************************************************************/
 
 Extern void EXPORT(suave)(ccount *pndim, ccount *pncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber *pnvec,
   creal *pepsrel, creal *pepsabs,
   cint *pflags, cint *pseed,
   cnumber *pmineval, cnumber *pmaxeval,
-  cnumber *pnnew, creal *pflatness,
+  cnumber *pnnew, cnumber *pnmin, creal *pflatness,
-  cchar *statefile,
+  cchar *statefile, Spin **pspin,
   count *pnregions, number *pneval, int *pfail,
   real *integral, real *error, real *prob, cint statefilelen)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = *pndim;
   t.ncomp = *pncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = *pnvec;
   t.epsrel = *pepsrel;
   t.epsabs = *pepsabs;
-  t.flags = *pflags;
+  t.flags = MaxVerbose(*pflags);
   t.seed = *pseed;
   t.mineval = *pmineval;
   t.maxeval = *pmaxeval;
   t.nnew = *pnnew;
+  t.nmin = IMax(*pnmin, 2);
   t.flatness = *pflatness;
   CString(t.statefile, statefile, statefilelen);
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
 
   *pfail = Integrate(&t, integral, error, prob);
   *pnregions = t.nregions;
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
 

src/external/libCuba/src/suave/decl.h

@@ -2,7 +2,7 @@
 	decl.h
 		Type declarations
 		this file is part of Suave
-		last modified 29 Jul 13 th
+		last modified 25 Nov 14 th
 */
 
 
@@ -35,23 +35,24 @@ typedef struct {
 typedef const Bounds cBounds;
 
 typedef int (*Integrand)(ccount *, creal *, ccount *, real *,
-  void *, creal *, cint *);
+  void *, cnumber *, cint *, creal *, cint *);
 
 typedef struct _this {
   count ndim, ncomp;
 #ifndef MLVERSION
   Integrand integrand;
   void *userdata;
+  number nvec;
 #ifdef HAVE_FORK
-  int ncores, *child;
-  real *frame;
   SHM_ONLY(int shmid;)
+  Spin *spin;
+  real *frame;
 #endif
 #endif
   real epsrel, epsabs;
   int flags, seed;
   number mineval, maxeval;
-  number nnew;
+  number nnew, nmin;
   real flatness;
   cchar *statefile;
   count nregions;

src/external/libCuba/src/vegas/Grid.c

@@ -2,7 +2,7 @@
 	Grid.c
 		utility functions for the Vegas grid
 		this file is part of Vegas
-		last modified 13 Dec 11 th
+		last modified 12 Mar 15 th
 */
 
 
@@ -70,7 +70,7 @@ static void RefineGrid(cThis *t, Grid grid, Grid margsum)
     real impfun = 0;
     if( margsum[bin] > 0 ) {
       creal r = margsum[bin]*norm;
-      avgperbin += impfun = pow((r - 1)/log(r), 1.5);
+      avgperbin += impfun = powx((r - 1)/log(r), 1.5);
     }
     imp[bin] = impfun;
   }

src/external/libCuba/src/vegas/Integrate.c

@@ -2,7 +2,7 @@
 	Integrate.c
 		integrate over the unit hypercube
 		this file is part of Vegas
-		last modified 8 Aug 13 th
+		last modified 13 Mar 15 th
 */
 
 
@@ -31,6 +31,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   if( VERBOSE > 1 ) {
     sprintf(out, "Vegas input parameters:\n"
       "  ndim " COUNT "\n  ncomp " COUNT "\n"
+      ML_NOT("  nvec " NUMBER "\n")
       "  epsrel " REAL "\n  epsabs " REAL "\n"
       "  flags %d\n  seed %d\n"
       "  mineval " NUMBER "\n  maxeval " NUMBER "\n"
@@ -38,7 +39,8 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       "  nbatch " NUMBER "\n  gridno %d\n"
       "  statefile \"%s\"",
       t->ndim, t->ncomp,
-      t->epsrel, t->epsabs,
+      ML_NOT(t->nvec,)
+      SHOW(t->epsrel), SHOW(t->epsabs),
       t->flags, t->seed,
       t->mineval, t->maxeval,
       t->nstart, t->nincrease, t->nbatch,
@@ -49,7 +51,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
   if( BadComponent(t) ) return -2;
   if( BadDimension(t) ) return -1;
 
-  FrameAlloc(t, ShmRm(t));
+  FrameAlloc(t, Master);
   ForkCores(t);
   Alloc(bins, t->nbatch*t->ndim);
 
@@ -68,12 +70,12 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
     t->rng.skiprandom(t, t->neval);
   }
 
-  if( ini ) {
+  if( ini | ZAPSTATE ) {
+    t->neval = 0;
     state->niter = 0;
     state->nsamples = t->nstart;
     FClear(state->cumul);
-    GetGrid(t, state_grid);
+    if( ini ) GetGrid(t, state_grid);
-    t->neval = 0;
   }
 
   /* main iteration loop */
@@ -143,7 +145,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       real avg = sigsq*(c->avgsum += w*c->sum);
 
       c->avg = LAST ? (sigsq = 1/w, c->sum) : avg;
-      c->err = sqrt(sigsq);
+      c->err = sqrtx(sigsq);
       fail |= (c->err > MaxErr(c->avg));
 
       if( state->niter == 0 ) c->guess = c->sum;
@@ -163,7 +165,8 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
       for( c = state->cumul, comp = 0; c < C; ++c )
         oe += sprintf(oe, "\n[" COUNT "] "
           REAL " +- " REAL "  \tchisq " REAL " (" COUNT " df)",
-          ++comp, c->avg, c->err, c->chisq, state->niter);
+          ++comp, SHOW(c->avg), SHOW(c->err),
+          SHOW(c->chisq), state->niter);
       Print(out);
     }
 
@@ -215,8 +218,7 @@ static int Integrate(This *t, real *integral, real *error, real *prob)
 abort:
   PutGrid(t, state_grid);
   free(bins);
-  WaitCores(t);
+  FrameFree(t, Master);
-  FrameFree(t);
 
   StateRemove(t);
 

src/external/libCuba/src/vegas/Vegas.c

@@ -1,8 +1,8 @@
 /*
 	Vegas.c
-		Vegas Monte-Carlo integration
+		Vegas Monte Carlo integration
 		by Thomas Hahn
-		last modified 17 Sep 13 th
+		last modified 25 Nov 14 th
 */
 
 
@@ -15,22 +15,27 @@
 /*********************************************************************/
 
 Extern void EXPORT(Vegas)(ccount ndim, ccount ncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber nvec,
   creal epsrel, creal epsabs, cint flags, cint seed,
   cnumber mineval, cnumber maxeval,
-  cnumber nstart, cnumber nincrease, cnumber nbatch,
+  cnumber nstart, cnumber nincrease,
-  cint gridno, cchar *statefile,
+  cnumber nbatch, cint gridno,
+  cchar *statefile, Spin **pspin,
   number *pneval, int *pfail,
   real *integral, real *error, real *prob)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = ndim;
   t.ncomp = ncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = nvec;
   t.epsrel = epsrel;
   t.epsabs = epsabs;
-  t.flags = flags;
+  t.flags = MaxVerbose(flags);
   t.seed = seed;
   t.mineval = mineval;
   t.maxeval = maxeval;
@@ -39,30 +44,38 @@ Extern void EXPORT(Vegas)(ccount ndim, ccount ncomp,
   t.nbatch = nbatch;
   t.gridno = gridno;
   t.statefile = statefile;
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
 
   *pfail = Integrate(&t, integral, error, prob);
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
 
 /*********************************************************************/
 
 Extern void EXPORT(vegas)(ccount *pndim, ccount *pncomp,
-  Integrand integrand, void *userdata,
+  Integrand integrand, void *userdata, cnumber *pnvec,
   creal *pepsrel, creal *pepsabs, cint *pflags, cint *pseed,
   cnumber *pmineval, cnumber *pmaxeval,
   cnumber *pnstart, cnumber *pnincrease, 
-  cnumber *pnbatch, cint *pgridno, cchar *statefile,
+  cnumber *pnbatch, cint *pgridno,
+  cchar *statefile, Spin **pspin,
   number *pneval, int *pfail,
   real *integral, real *error, real *prob, cint statefilelen)
 {
   This t;
+
+  VerboseInit();
+
   t.ndim = *pndim;
   t.ncomp = *pncomp;
   t.integrand = integrand;
   t.userdata = userdata;
+  t.nvec = *pnvec;
   t.epsrel = *pepsrel;
   t.epsabs = *pepsabs;
-  t.flags = *pflags;
+  t.flags = MaxVerbose(*pflags);
   t.seed = *pseed;
   t.mineval = *pmineval;
   t.maxeval = *pmaxeval;
@@ -71,8 +84,11 @@ Extern void EXPORT(vegas)(ccount *pndim, ccount *pncomp,
   t.nbatch = *pnbatch;
   t.gridno = *pgridno;
   CString(t.statefile, statefile, statefilelen);
+  FORK_ONLY(t.spin = Invalid(pspin) ? NULL : *pspin;)
 
   *pfail = Integrate(&t, integral, error, prob);
   *pneval = t.neval;
+
+  WaitCores(&t, pspin);
 }
 

src/external/libCuba/src/vegas/decl.h

@@ -2,7 +2,7 @@
 	decl.h
 		Type declarations
 		this file is part of Vegas
-		last modified 21 Dec 11 th
+		last modified 21 Jul 14 th
 */
 
 
@@ -29,16 +29,17 @@ typedef struct {
 typedef const Cumulants cCumulants;
 
 typedef int (*Integrand)(ccount *, creal *, ccount *, real *,
-  void *, creal *, cint *);
+  void *, cnumber *, cint *, creal *, cint *);
 
 typedef struct _this {
   count ndim, ncomp;
 #ifndef MLVERSION
   Integrand integrand;
   void *userdata;
+  number nvec;
 #ifdef HAVE_FORK
-  int ncores, *child;
   SHM_ONLY(int shmid;)
+  Spin *spin;
 #endif
 #endif
   real *frame;