From 2818739ec656c4f6d7b7664830843dba523f5fca Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Fri, 20 Dec 2013 09:58:23 +0000 Subject: [PATCH] upgrade of cuba to version 3.2. Merge in from BMW. --- ChangeLog | 1 + configure.ac | 26 +- src/external/BMWtools/BMWIntegrator.cpp | 17 +- src/external/libCuba/src/Makefile.am | 15 +- src/external/libCuba/src/common/CSample.c | 67 +++ src/external/libCuba/src/common/ChiSquare.c | 20 - src/external/libCuba/src/common/Erf.c | 20 - src/external/libCuba/src/common/Fork.c | 422 ++++++++++++++ src/external/libCuba/src/common/MSample.c | 90 +++ src/external/libCuba/src/common/Makefile.am | 13 + src/external/libCuba/src/common/Random.c | 34 +- src/external/libCuba/src/common/WorkerIni.c | 37 ++ src/external/libCuba/src/common/debug.c | 105 ---- src/external/libCuba/src/common/stddecl.h | 319 ++++++++-- src/external/libCuba/src/cuba.h | 33 +- src/external/libCuba/src/cuhre/Cuhre.c | 56 +- src/external/libCuba/src/cuhre/Integrate.c | 251 ++++---- src/external/libCuba/src/cuhre/Makefile.am | 13 + src/external/libCuba/src/cuhre/Rule.c | 210 +++---- src/external/libCuba/src/cuhre/common.c | 28 +- src/external/libCuba/src/cuhre/decl.h | 59 +- src/external/libCuba/src/cuhre/util.c | 37 -- src/external/libCuba/src/divonne/Divonne.c | 100 +--- src/external/libCuba/src/divonne/Explore.c | 96 +--- .../libCuba/src/divonne/FindMinimum.c | 77 +-- src/external/libCuba/src/divonne/Integrate.c | 543 +++++++++--------- src/external/libCuba/src/divonne/Iterate.c | 132 +++++ .../libCuba/src/divonne/KorobovCoeff.c | 20 - src/external/libCuba/src/divonne/Makefile.am | 13 + src/external/libCuba/src/divonne/Rule.c | 207 +++---- src/external/libCuba/src/divonne/Sample.c | 102 ++-- src/external/libCuba/src/divonne/Split.c | 179 ++---- src/external/libCuba/src/divonne/common.c | 76 ++- src/external/libCuba/src/divonne/decl.h | 125 ++-- src/external/libCuba/src/divonne/util.c | 51 -- src/external/libCuba/src/suave/Fluct.c | 46 +- src/external/libCuba/src/suave/Grid.c | 29 +- src/external/libCuba/src/suave/Integrate.c | 306 +++++----- src/external/libCuba/src/suave/Makefile.am | 13 + src/external/libCuba/src/suave/Sample.c | 81 +-- src/external/libCuba/src/suave/Suave.c | 53 +- src/external/libCuba/src/suave/common.c | 44 +- src/external/libCuba/src/suave/decl.h | 52 +- src/external/libCuba/src/suave/util.c | 43 -- src/external/libCuba/src/vegas/Grid.c | 26 +- src/external/libCuba/src/vegas/Integrate.c | 201 +++---- src/external/libCuba/src/vegas/Makefile.am | 13 + src/external/libCuba/src/vegas/Vegas.c | 62 +- src/external/libCuba/src/vegas/common.c | 32 +- src/external/libCuba/src/vegas/decl.h | 29 +- src/external/libCuba/src/vegas/util.c | 46 -- 51 files changed, 2454 insertions(+), 2216 deletions(-) create mode 100644 src/external/libCuba/src/common/CSample.c create mode 100644 src/external/libCuba/src/common/Fork.c create mode 100644 src/external/libCuba/src/common/MSample.c create mode 100644 src/external/libCuba/src/common/Makefile.am create mode 100644 src/external/libCuba/src/common/WorkerIni.c delete mode 100644 src/external/libCuba/src/common/debug.c create mode 100644 src/external/libCuba/src/cuhre/Makefile.am delete mode 100644 src/external/libCuba/src/cuhre/util.c create mode 100644 src/external/libCuba/src/divonne/Iterate.c create mode 100644 src/external/libCuba/src/divonne/Makefile.am delete mode 100644 src/external/libCuba/src/divonne/util.c create mode 100644 src/external/libCuba/src/suave/Makefile.am delete mode 100644 src/external/libCuba/src/suave/util.c create mode 100644 src/external/libCuba/src/vegas/Makefile.am delete mode 100644 src/external/libCuba/src/vegas/util.c diff --git a/ChangeLog b/ChangeLog index 62d57de8..b3e4837c 100644 --- a/ChangeLog +++ b/ChangeLog @@ -60,6 +60,7 @@ FIXED 2012-09-23 fixed wrong chisq output in musrview if expected chisq is present. FIXED 2012-05-30 fixed RRF bug in single histo plotting. FIXED 2012-05-18 fixed wrong forward/backward tag for ROOT-PPC (MUSR-215) +CHANGED 2013-12-20 upgrade of cuba to version 3.2. Merge in from BMW CHANGED 2013-12-16 prettyfied the Noakes-Kalvius formulae. Furthermore added a sub-folder with cross checks for these formulae. CHANGED 2013-11-12 changed normalization of log max likelihood according to S. diff --git a/configure.ac b/configure.ac index 3620058f..5bff03eb 100644 --- a/configure.ac +++ b/configure.ac @@ -69,8 +69,8 @@ PLUGIN_MINOR_VERSION=0 PLUGIN_MICRO_VERSION=0 #release versioning -CUBA_MAJOR_VERSION=2 -CUBA_MINOR_VERSION=1 +CUBA_MAJOR_VERSION=3 +CUBA_MINOR_VERSION=2 CUBA_MICRO_VERSION=0 #API version @@ -200,9 +200,9 @@ dnl ----------------------------------------------- dnl Automake initialization and program checks dnl ----------------------------------------------- -AM_INIT_AUTOMAKE(@PACKAGE_NAME@, @PACKAGE_VERSION@) -m4_ifdef([AM_SILENT_RULES], - [AM_SILENT_RULES([yes])]) +AM_INIT_AUTOMAKE +# m4_ifdef([AM_SILENT_RULES], +# [AM_SILENT_RULES([yes])]) AC_CONFIG_HEADER([config.h]) AC_LANG([C++]) AC_PROG_LN_S @@ -657,6 +657,17 @@ AC_ARG_ENABLE([BMWlibs], [AC_HELP_STRING([--enable-BMWlibs],[build optional BMW AC_CHECK_FUNCS([powl]) AC_CHECK_FUNCS([erf]) +# AC_FUNC_FORK + AC_FUNC_ALLOCA + + AC_DEFUN([chk_shmget], [dnl + AC_REQUIRE([AC_CANONICAL_HOST]) + AS_CASE([$host_os], + [*cygwin*], [], + [AC_CHECK_FUNCS([shmget])]) + ]) + chk_shmget + AC_CHECK_FUNCS([getloadavg]) MAXDIM=${MAXDIM:-16} AC_ARG_WITH(maxdim, @@ -1112,6 +1123,11 @@ AC_CONFIG_FILES([Makefile \ src/external/libCuba/Makefile \ src/external/libCuba/src/Makefile \ src/external/libCuba/src/cuba.pc \ + src/external/libCuba/src/cuhre/Makefile \ + src/external/libCuba/src/divonne/Makefile \ + src/external/libCuba/src/suave/Makefile \ + src/external/libCuba/src/vegas/Makefile \ + src/external/libCuba/src/common/Makefile \ src/external/BMWtools/Makefile \ src/external/libFitPofB/Makefile \ src/external/libFitPofB/classes/Makefile \ diff --git a/src/external/BMWtools/BMWIntegrator.cpp b/src/external/BMWtools/BMWIntegrator.cpp index 72c2d523..78af1928 100644 --- a/src/external/BMWtools/BMWIntegrator.cpp +++ b/src/external/BMWtools/BMWIntegrator.cpp @@ -34,6 +34,7 @@ #define USERDATA NULL #define SEED 0 +#define STATEFILE NULL std::vector TDWaveGapIntegralCuhre::fPar; @@ -60,7 +61,7 @@ double TDWaveGapIntegralCuhre::IntegrateFunc() Cuhre(fNDim, NCOMP, Integrand, USERDATA, EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, - KEY, + KEY, STATEFILE, &nregions, &neval, &fail, integral, error, prob); return integral[0]; @@ -111,7 +112,7 @@ double TCosSqDWaveGapIntegralCuhre::IntegrateFunc() Cuhre(fNDim, NCOMP, Integrand, USERDATA, EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, - KEY, + KEY, STATEFILE, &nregions, &neval, &fail, integral, error, prob); return integral[0]; @@ -162,7 +163,7 @@ double TSinSqDWaveGapIntegralCuhre::IntegrateFunc() Cuhre(fNDim, NCOMP, Integrand, USERDATA, EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, - KEY, + KEY, STATEFILE, &nregions, &neval, &fail, integral, error, prob); return integral[0]; @@ -213,7 +214,7 @@ double TAnSWaveGapIntegralCuhre::IntegrateFunc() Cuhre(fNDim, NCOMP, Integrand, USERDATA, EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, - KEY, + KEY, STATEFILE, &nregions, &neval, &fail, integral, error, prob); return integral[0]; @@ -272,7 +273,7 @@ double TAnSWaveGapIntegralDivonne::IntegrateFunc() Divonne(fNDim, NCOMP, Integrand, USERDATA, EPSREL, EPSABS, VERBOSE, SEED, MINEVAL, MAXEVAL, KEY1, KEY2, KEY3, MAXPASS, BORDER, MAXCHISQ, MINDEVIATION, - NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL, + NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL, STATEFILE, &nregions, &neval, &fail, integral, error, prob); return integral[0]; @@ -324,7 +325,7 @@ double TAnSWaveGapIntegralSuave::IntegrateFunc() Suave(fNDim, NCOMP, Integrand, USERDATA, EPSREL, EPSABS, VERBOSE | LAST, SEED, MINEVAL, MAXEVAL, - NNEW, FLATNESS, + NNEW, FLATNESS, STATEFILE, &nregions, &neval, &fail, integral, error, prob); return integral[0]; @@ -375,7 +376,7 @@ double TNonMonDWave1GapIntegralCuhre::IntegrateFunc() Cuhre(fNDim, NCOMP, Integrand, USERDATA, EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, - KEY, + KEY, STATEFILE, &nregions, &neval, &fail, integral, error, prob); return integral[0]; @@ -426,7 +427,7 @@ double TNonMonDWave2GapIntegralCuhre::IntegrateFunc() Cuhre(fNDim, NCOMP, Integrand, USERDATA, EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, - KEY, + KEY, STATEFILE, &nregions, &neval, &fail, integral, error, prob); return integral[0]; diff --git a/src/external/libCuba/src/Makefile.am b/src/external/libCuba/src/Makefile.am index 2df617cd..afff0736 100644 --- a/src/external/libCuba/src/Makefile.am +++ b/src/external/libCuba/src/Makefile.am @@ -1,25 +1,18 @@ ## Process this file with automake to create Makefile.in ## $Id$ -h_sources = cuba.h - -c_sources = cuhre/Cuhre.c \ - divonne/Divonne.c \ - suave/Suave.c \ - vegas/Vegas.c +SUBDIRS = cuhre divonne suave vegas common include_HEADERS = cuba.h -AM_CPPFLAGS = -I. -Icommon -AM_CFLAGS = $(LOCAL_CUBA_LIB_CFLAGS) - AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS) CLEANFILES = common/*~ cuhre/*~ divonne/*~ suave/*~ vegas/*~ *~ core lib_LTLIBRARIES = libcuba.la -libcuba_la_SOURCES = $(h_sources) $(c_sources) -libcuba_la_LIBADD = -lm +libcuba_la_SOURCES = + +libcuba_la_LIBADD = 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 diff --git a/src/external/libCuba/src/common/CSample.c b/src/external/libCuba/src/common/CSample.c new file mode 100644 index 00000000..5705be06 --- /dev/null +++ b/src/external/libCuba/src/common/CSample.c @@ -0,0 +1,67 @@ +/* + CSample.c + the serial sampling routine + for the C versions of the Cuba routines + by Thomas Hahn + last modified 19 Dec 11 th +*/ + + +static inline number SampleRaw(cThis *t, number n, creal *x, real *f + VES_ONLY(, creal *w, ccount iter)) +{ + for( ; n; --n ) { + if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata + VES_ONLY(, w++, &iter) + DIV_ONLY(, &t->phase)) == ABORT ) return -1; + x += t->ndim; + f += t->ncomp; + } + return 0; +} + +/*********************************************************************/ + +static inline void DoSampleSerial(This *t, cnumber n, creal *x, real *f + VES_ONLY(, creal *w, ccount iter)) +{ + t->neval += n; + if( SampleRaw(t, n, x, f VES_ONLY(, w, iter)) ) + longjmp(t->abort, -99); +} + +/*********************************************************************/ + +#ifdef HAVE_FORK + +static void DoSample(This *t, number n, creal *x, real *f + VES_ONLY(, creal *w, ccount iter)); +DIV_ONLY(static int Explore(This *t, cint iregion);) + +#else + +#define DoSample DoSampleSerial +#define Explore ExploreSerial +#define ForkCores(t) +#define WaitCores(t) + +#endif + +#ifdef DIVONNE +static inline count SampleExtra(This *t, cBounds *b) +{ + number n = t->nextra; + t->peakfinder(&t->ndim, b, &n, t->xextra); + DoSample(t, n, t->xextra, t->fextra); + return n; +} +#endif + +#include "common.c" + +#ifdef HAVE_FORK +#include "Fork.c" +#endif + +#include "Integrate.c" + diff --git a/src/external/libCuba/src/common/ChiSquare.c b/src/external/libCuba/src/common/ChiSquare.c index 2d8f5d65..baee4131 100644 --- a/src/external/libCuba/src/common/ChiSquare.c +++ b/src/external/libCuba/src/common/ChiSquare.c @@ -6,26 +6,6 @@ last modified 9 Feb 05 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - #ifdef HAVE_ERF #define Erf erf #else diff --git a/src/external/libCuba/src/common/Erf.c b/src/external/libCuba/src/common/Erf.c index 0bfe1269..635b6e5c 100644 --- a/src/external/libCuba/src/common/Erf.c +++ b/src/external/libCuba/src/common/Erf.c @@ -7,26 +7,6 @@ last modified 8 Feb 05 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - static real Erfc(creal x) { diff --git a/src/external/libCuba/src/common/Fork.c b/src/external/libCuba/src/common/Fork.c new file mode 100644 index 00000000..ac7d3047 --- /dev/null +++ b/src/external/libCuba/src/common/Fork.c @@ -0,0 +1,422 @@ +/* + Fork.c + the parallel sampling routine + for the C versions of the Cuba routines + by Thomas Hahn + last modified 25 Sep 13 th +*/ + +#define MINSLICE 10 +#define MINCORES 1 +/*#define MINCORES 2*/ + +typedef struct { + number n, m, i; + VES_ONLY(count iter;) + DIV_ONLY(int phase SHM_ONLY(, shmid);) +} Slice; + +workerini cubaini; + +#if defined HAVE_SHMGET && (defined SUAVE || defined DIVONNE) +#define FRAMECOPY +#endif + +#ifdef DEBUG +#define TERM_RED "\e[31m" +#define TERM_BLUE "\e[34m" +#define TERM_RESET "\e[0m\n" +#define MASTER(s, ...) \ +fprintf(stderr, TERM_RED ROUTINE " master %d(%d): " s TERM_RESET, core, getpid(), ##__VA_ARGS__) +#define WORKER(s, ...) \ +fprintf(stderr, TERM_BLUE ROUTINE " worker %d(%d): " s TERM_RESET, core, getpid(), ##__VA_ARGS__) +#else +#define MASTER(s, ...) +#define WORKER(s, ...) +#endif + +/*********************************************************************/ + +#ifndef MSG_WAITALL +/* Windows */ +#define MSG_WAITALL 0 +#endif + +static inline int readsock(int fd, void *data, size_t n) +{ + ssize_t got; + size_t remain = n; + do got = recv(fd, data, remain, MSG_WAITALL); + while( got > 0 && (data += got, remain -= got) > 0 ); + return got; +} + +static inline int writesock(int fd, const void *data, size_t n) +{ + ssize_t got; + size_t remain = n; + do got = send(fd, data, remain, MSG_WAITALL); + while( got > 0 && (data += got, remain -= got) > 0 ); + return got; +} + +/*********************************************************************/ + +static 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 ) { + Slice slice; + 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); + Print(s); + } + + 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); + for( core = 0; core < t->ncores; ++core ) { + int fd[2]; + pid_t pid; + assert( + socketpair(AF_LOCAL, SOCK_STREAM, 0, fd) != -1 && + (pid = fork()) != -1 ); + if( pid == 0 ) { + close(fd[0]); + DoChild(t, core, fd[1]); + } + MASTER("forked pid %d pipe %d(master) -> %d(worker)", + pid, fd[0], fd[1]); + close(fd[1]); + t->child[core] = fd[0]; + } +} + +/*********************************************************************/ + +static inline void WaitCores(cThis *t) +{ + if( t->ncores >= MINCORES ) { + int core; + pid_t pid; + for( core = 0; core < t->ncores; ++core ) { + MASTER("closing fd %d", t->child[core]); + close(t->child[core]); + } + free(t->child); + for( core = 0; core < t->ncores; ++core ) { + MASTER("waiting for child"); + wait(&pid); + MASTER("pid %d terminated", pid); + } + } +} + diff --git a/src/external/libCuba/src/common/MSample.c b/src/external/libCuba/src/common/MSample.c new file mode 100644 index 00000000..4c146346 --- /dev/null +++ b/src/external/libCuba/src/common/MSample.c @@ -0,0 +1,90 @@ +/* + MSample.c + the sampling routine for the + Mathematica versions of the Cuba routines + by Thomas Hahn + last modified 19 Mar 12 th +*/ + + +static void DoSample(This *t, cnumber n, real *x, real *f + VES_ONLY(, real *w, ccount iter)) +{ + real *mma_f; + long 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); + VES_ONLY(MLPutRealList(stdlink, w, n); + MLPutInteger(stdlink, iter);) + DIV_ONLY(MLPutInteger(stdlink, t->phase);) + MLEndPacket(stdlink); + + MLNextPacket(stdlink); + if( !MLGetRealList(stdlink, &mma_f, &mma_n) ) { + MLClearError(stdlink); + MLNewPacket(stdlink); + longjmp(t->abort, -99); + } + + t->neval += mma_n; + + if( mma_n != n*t->ncomp ) { + MLDisownRealList(stdlink, mma_f, mma_n); + longjmp(t->abort, -3); + } + + Copy(f, mma_f, n*t->ncomp); + MLDisownRealList(stdlink, mma_f, mma_n); +} + +/*********************************************************************/ + +#ifdef DIVONNE +#define Explore ExploreSerial + +static count SampleExtra(This *t, cBounds *b) +{ + count n, nget; + real *mma_f; + long mma_n; + + MLPutFunction(stdlink, "EvaluatePacket", 1); + MLPutFunction(stdlink, "Cuba`Divonne`findpeak", 2); + MLPutRealList(stdlink, (real *)b, 2*t->ndim); + MLPutInteger(stdlink, t->phase); + MLEndPacket(stdlink); + + MLNextPacket(stdlink); + if( !MLGetRealList(stdlink, &mma_f, &mma_n) ) { + MLClearError(stdlink); + MLNewPacket(stdlink); + longjmp(t->abort, -99); + } + + t->neval += nget = mma_n/(t->ndim + t->ncomp); + + n = IMin(nget, t->nextra); + if( n ) { + Copy(t->xextra, mma_f, n*t->ndim); + Copy(t->fextra, mma_f + nget*t->ndim, n*t->ncomp); + } + + MLDisownRealList(stdlink, mma_f, mma_n); + + return n; +} +#endif + +/*********************************************************************/ + +#include "common.c" + +#define ForkCores(t) +#define WaitCores(t) + +#include "Integrate.c" + diff --git a/src/external/libCuba/src/common/Makefile.am b/src/external/libCuba/src/common/Makefile.am new file mode 100644 index 00000000..e21290e0 --- /dev/null +++ b/src/external/libCuba/src/common/Makefile.am @@ -0,0 +1,13 @@ +## Process this file with automake to create Makefile.in +## $Id: Makefile.am 4808 2011-03-21 14:33:34Z l_wojek $ + +c_sources = WorkerIni.c + +AM_CPPFLAGS = -I. -I.. -I../common -DNOUNDERSCORE +AM_CFLAGS = $(LOCAL_CUBA_LIB_CFLAGS) +AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS) + +noinst_LTLIBRARIES = libworkerini.la + +libworkerini_la_SOURCES = $(c_sources) +libworkerini_la_LDFLAGS = $(AM_LDFLAGS) \ No newline at end of file diff --git a/src/external/libCuba/src/common/Random.c b/src/external/libCuba/src/common/Random.c index 72abf187..4c2777c2 100644 --- a/src/external/libCuba/src/common/Random.c +++ b/src/external/libCuba/src/common/Random.c @@ -1,29 +1,9 @@ /* Random.c quasi- and pseudo-random-number generation - last modified 8 Jun 10 th + last modified 7 Aug 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - /* PART 1: Sobol quasi-random-number generator @@ -142,7 +122,7 @@ static inline void SobolIni(This *t) } t->rng.sobol.seq = 0; - VecClear(t->rng.sobol.prev); + XClear(t->rng.sobol.prev); t->rng.getrandom = SobolGet; t->rng.skiprandom = SobolSkip; @@ -300,10 +280,10 @@ static void RanluxGet(This *t, real *x) cint nskip = (--t->rng.ranlux.n24 >= 0) ? 0 : (t->rng.ranlux.n24 = 24, t->rng.ranlux.nskip); cint s = RanluxInt(t, 1 + nskip); - x[dim] = s*0x1p-24; + x[dim] = ldexp(s, -24); /* small numbers (with less than 12 significant bits) are "padded" */ if( s < (1 << 12) ) - x[dim] += t->rng.ranlux.state[t->rng.ranlux.j24]*0x1p-48; + x[dim] += ldexp(t->rng.ranlux.state[t->rng.ranlux.j24], -48); } } @@ -320,11 +300,11 @@ static inline void RanluxIni(This *t) cint skip[] = {24, 48, 97, 223, 389, 223, 223, 223, 223, 223, 223, 223, 223, 223, 223, 223, 223, 223, 223, 223, 223, 223, 223, 223, 223}; - state_t seed = t->seed; - state_t level = RNG; + int seed = t->seed; + int level = RNG; count i; - if( level < sizeof skip ) level = skip[level]; + if( level < Elements(skip) ) level = skip[level]; t->rng.ranlux.nskip = level - 24; t->rng.ranlux.i24 = 23; diff --git a/src/external/libCuba/src/common/WorkerIni.c b/src/external/libCuba/src/common/WorkerIni.c new file mode 100644 index 00000000..52c9978a --- /dev/null +++ b/src/external/libCuba/src/common/WorkerIni.c @@ -0,0 +1,37 @@ +/* + 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); +} + diff --git a/src/external/libCuba/src/common/debug.c b/src/external/libCuba/src/common/debug.c deleted file mode 100644 index e1169ac8..00000000 --- a/src/external/libCuba/src/common/debug.c +++ /dev/null @@ -1,105 +0,0 @@ -/*************************************************************************** - * Copyright (C) 2004-2009 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - -#define DEBFile stdout - -#define DEBF "%9.5f" -//#define DEBF "%21.17f" -//#define DEBF "%a" - - -#ifdef DEBFile - -#define DEBOpen() -#define DEBClose() - -#else - -FILE *DEBFile; - -static inline void DEBOpen() -{ -#ifdef MLVERSION - DEBFile = fopen("log-mma", "w"); -#else - DEBFile = fopen("log-c", "w"); -#endif -} - -static inline void DEBClose() -{ - fclose(DEBFile); -} - -#endif - - -#define DEB(...) fprintf(DEBFile, __VA_ARGS__); fflush(DEBFile) - - -static inline void DEBVec(const char *s, creal *d) -{ - char space[strlen(s) + 2]; - count dim; - - memset(space, ' ', sizeof(space)); - space[sizeof(space) - 1] = 0; - - DEB("%s=" DEBF "\n", s, d[0]); - for( dim = 1; dim < ndim_; ++dim ) - DEB("%s" DEBF "\n", space, d[dim]); -} - - -/* -static inline void DEBRegion(const char *s, cBounds *b) -{ - char space[strlen(s) + 3]; - count dim; - - memset(space, ' ', sizeof(space)); - space[sizeof(space) - 1] = 0; - - DEB("%s: " DEBF " - " DEBF "\n", s, b[0].lower, b[0].upper); - for( dim = 1; dim < ndim_; ++dim ) - DEB("%s" DEBF " - " DEBF "\n", space, b[dim].lower, b[dim].upper); -} -*/ - - -static inline void DEBMem(const char *s) -{ - int kbytes = -1; - FILE *f; - char procfile[128]; - - sprintf(procfile, "/proc/%d/status", getpid()); - f = fopen(procfile, "r"); - while( !feof(f) ) { - char s[128]; - *s = 0; - fgets(s, sizeof(s), f); - if( sscanf(s, "VmSize: %d", &kbytes) == 1 ) break; - } - fclose(f); - - DEB("MEM %s: %d kbytes\n", s, kbytes); -} - diff --git a/src/external/libCuba/src/common/stddecl.h b/src/external/libCuba/src/common/stddecl.h index e43ba2b8..fed006f3 100644 --- a/src/external/libCuba/src/common/stddecl.h +++ b/src/external/libCuba/src/common/stddecl.h @@ -1,37 +1,20 @@ /* stddecl.h - Type declarations common to all Cuba routines - last modified 16 Jun 10 th + declarations common to all Cuba routines + last modified 17 Sep 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - -#ifndef __stddecl_h__ -#define __stddecl_h__ +#ifndef _stddecl_h_ +#define _stddecl_h_ #ifdef HAVE_CONFIG_H #include "config.h" #endif +#define _BSD_SOURCE +#define _XOPEN_SOURCE + #include #include #include @@ -39,28 +22,83 @@ #include #include #include +#include #include #include #include +#include +#ifdef HAVE_FORK +#include +#include +#ifdef HAVE_SHMGET +#include +#include +#endif +#endif +#ifdef HAVE_ALLOCA_H +#include +#elif defined __GNUC__ +#define alloca __builtin_alloca +#elif defined _AIX +#define alloca __alloca +#elif defined _MSC_VER +#include +#define alloca _alloca +#else +#include +#ifdef __cplusplus +extern "C" +#endif +void *alloca (size_t); +#endif #ifndef NDIM #define NDIM t->ndim -#endif -#ifndef NCOMP -#define NCOMP t->ncomp +#define MAXDIM 1024 +#else +#define MAXDIM NDIM #endif +#ifndef NCOMP +#define NCOMP t->ncomp +#define MAXCOMP 1024 +#else +#define MAXCOMP NCOMP +#endif + +#if defined(VEGAS) || defined(SUAVE) +#define VES_ONLY(...) __VA_ARGS__ +#define NW 1 +#else +#define VES_ONLY(...) +#define NW 0 +#endif + +#ifdef DIVONNE +#define DIV_ONLY(...) __VA_ARGS__ +#else +#define DIV_ONLY(...) +#endif + +#define SAMPLESIZE (NW + t->ndim + t->ncomp)*sizeof(real) #define VERBOSE (t->flags & 3) #define LAST (t->flags & 4) #define SHARPEDGES (t->flags & 8) +#define KEEPFILE (t->flags & 16) #define REGIONS (t->flags & 128) #define RNG (t->flags >> 8) #define INFTY DBL_MAX -#define NOTZERO 0x1p-104 +#if __STDC_VERSION__ >= 199901L +#define POW2(n) 0x1p-##n +#else +#define POW2(n) ldexp(1., -n) +#endif + +#define NOTZERO POW2(104) #define ABORT -999 @@ -68,15 +106,17 @@ #define Copy(d, s, n) memcpy(d, s, (n)*sizeof(*(d))) -#define VecCopy(d, s) Copy(d, s, t->ndim) +#define Move(d, s, n) memmove(d, s, (n)*sizeof(*(d))) -#define ResCopy(d, s) Copy(d, s, t->ncomp) +#define XCopy(d, s) Copy(d, s, t->ndim) + +#define FCopy(d, s) Copy(d, s, t->ncomp) #define Clear(d, n) memset(d, 0, (n)*sizeof(*(d))) -#define VecClear(d) Clear(d, t->ndim) +#define XClear(d) Clear(d, t->ndim) -#define ResClear(d) Clear(d, t->ncomp) +#define FClear(d) Clear(d, t->ncomp) #define Zap(d) memset(d, 0, sizeof(d)) @@ -90,14 +130,155 @@ #define reallocset(p, n) (p = realloc(p, n)) #endif -#define ChkAlloc(r) if( r == NULL ) { \ - fprintf(stderr, "Out of memory in " __FILE__ " line %d.\n", __LINE__); \ - exit(1); \ +#define Abort(s) abort1(s, __LINE__) +#define abort1(s, line) abort2(s, line) +#define abort2(s, line) { perror(s " " __FILE__ "(" #line ")"); exit(1); } + +#define Die(p) if( (p) == NULL ) Abort("malloc") + +#define MemAlloc(p, n) Die(mallocset(p, n)) +#define ReAlloc(p, n) Die(reallocset(p, n)) +#define Alloc(p, n) MemAlloc(p, (n)*sizeof(*p)) + +#if __STDC_VERSION__ >= 199901L +#define Sized(type, var, size) char var##_[size]; type *var = (type *)var##_ +#define Vector(type, var, n1) type var[n1] +#define Array(type, var, n1, n2) type var[n1][n2] +#else +#define Sized(type, var, size) type *var = alloca(size) +#define Vector(type, var, n1) type *var = alloca((n1)*sizeof(type)) +#define Array(type, var, n1, n2) type (*var)[n2] = alloca((n1)*(n2)*sizeof(type)) +#endif + +#define FORK_ONLY(...) +#define SHM_ONLY(...) +#define ShmAlloc(...) +#define ShmFree(...) + +#ifdef MLVERSION +#define ML_ONLY(...) __VA_ARGS__ +#else +#define ML_ONLY(...) + +#ifdef HAVE_FORK +#undef FORK_ONLY +#define FORK_ONLY(...) __VA_ARGS__ + +#ifdef HAVE_SHMGET +#undef SHM_ONLY +#define SHM_ONLY(...) __VA_ARGS__ + +#define ShmMap(t, ...) if( t->shmid != -1 ) { \ + t->frame = shmat(t->shmid, NULL, 0); \ + if( t->frame == (void *)-1 ) Abort("shmat"); \ + __VA_ARGS__ \ } -#define Alloc(p, n) MemAlloc(p, (n)*sizeof(*p)) -#define MemAlloc(p, n) ChkAlloc(mallocset(p, n)) -#define ReAlloc(p, n) ChkAlloc(reallocset(p, n)) +#define ShmRm(t) shmctl(t->shmid, IPC_RMID, NULL); + +#undef ShmAlloc +#define ShmAlloc(t, ...) \ + t->shmid = shmget(IPC_PRIVATE, t->nframe*SAMPLESIZE, IPC_CREAT | 0600); \ + ShmMap(t, __VA_ARGS__) + +#undef ShmFree +#define ShmFree(t, ...) if( t->shmid != -1 ) { \ + shmdt(t->frame); \ + __VA_ARGS__ \ +} + +#endif +#endif +#endif + +#define FrameAlloc(t, ...) \ + SHM_ONLY(ShmAlloc(t, __VA_ARGS__) else) \ + MemAlloc(t->frame, t->nframe*SAMPLESIZE); + +#define FrameFree(t, ...) DIV_ONLY(if( t->nframe )) { \ + SHM_ONLY(ShmFree(t, __VA_ARGS__) else) \ + free(t->frame); \ +} + + +#define StateDecl \ +char *statefile_tmp = NULL, *statefile_XXXXXX = NULL; \ +int statemsg = VERBOSE; \ +ssize_t ini = 1; \ +struct stat st + +#define StateSetup(t) if( (t)->statefile ) { \ + if( *(t)->statefile == 0 ) (t)->statefile = NULL; \ + else { \ + ccount len = strlen((t)->statefile); \ + statefile_tmp = alloca(len + 8); \ + strcpy(statefile_tmp, (t)->statefile); \ + statefile_XXXXXX = statefile_tmp + len; \ + } \ +} + +typedef long long int signature_t; + +#define StateSignature(t, i) (0x41425543 + \ + ((signature_t)(i) << 60) + \ + ((signature_t)(t)->ncomp << 48) + \ + ((signature_t)(t)->ndim << 32)) + +#define StateReadTest(t) (t)->statefile && \ + stat((t)->statefile, &st) == 0 && (st.st_mode & 0400) + +#define StateReadOpen(t, fd) do { \ + int fd; \ + if( (fd = open((t)->statefile, O_RDONLY)) != -1 ) { \ + do + +#define StateRead(fd, buf, size) \ + ini += size - read(fd, buf, size) + +#define StateReadClose(t, fd) \ + while( (--ini, 0) ); \ + close(fd); \ + } \ + if( ini | statemsg ) { \ + char s[512]; \ + sprintf(s, ini ? \ + "\nError restoring state from %s, starting from scratch." : \ + "\nRestored state from %s.", (t)->statefile); \ + Print(s); \ + } \ +} while( 0 ) + + +#define StateWriteTest(t) ((t)->statefile) + +#define StateWriteOpen(t, fd) do { \ + ssize_t fail = 1; \ + int fd; \ + strcpy(statefile_XXXXXX, "-XXXXXX"); \ + if( (fd = mkstemp(statefile_tmp)) != -1 ) { \ + do + +#define StateWrite(fd, buf, size) \ + fail += size - write(fd, buf, size) + +#define StateWriteClose(t, fd) \ + while( (--fail, 0) ); \ + close(fd); \ + if( fail == 0 ) fail |= rename(statefile_tmp, (t)->statefile); \ + } \ + if( fail | statemsg ) { \ + char s[512]; \ + sprintf(s, fail ? \ + "\nError saving state to %s." : \ + "\nSaved state to %s.", (t)->statefile); \ + Print(s); \ + statemsg &= fail & -2; \ + } \ +} while( 0 ) + + +#define StateRemove(t) \ +if( fail == 0 && (t)->statefile && KEEPFILE == 0 ) unlink((t)->statefile) #ifdef __cplusplus @@ -115,6 +296,8 @@ typedef const int cint; typedef const long clong; +typedef const size_t csize_t; + #define COUNT "%d" typedef /*unsigned*/ int count; typedef const count ccount; @@ -144,6 +327,15 @@ typedef /*long*/ double real; typedef const real creal; +typedef void (*subroutine)(); + +typedef struct { + subroutine initfun; + void *initarg; + subroutine exitfun; + void *exitarg; +} workerini; + struct _this; @@ -179,33 +371,42 @@ typedef struct { } RNGState; -#ifdef UNDERSCORE -#define SUFFIX(s) s##_ -#else +#if NOUNDERSCORE #define SUFFIX(s) s +#else +#define SUFFIX(s) s##_ #endif #define EXPORT(s) EXPORT_(PREFIX(s)) #define EXPORT_(s) SUFFIX(s) -static inline real Sq(creal x) -{ +#define CString(cs, fs, len) { \ + char *_s = NULL; \ + if( fs ) { \ + int _l = len; \ + while( _l > 0 && fs[_l - 1] == ' ' ) --_l; \ + if( _l > 0 && (_s = alloca(_l + 1)) ) { \ + memcpy(_s, fs, _l); \ + _s[_l] = 0; \ + } \ + } \ + cs = _s; \ +} + +static inline real Sq(creal x) { return x*x; } -static inline real Min(creal a, creal b) -{ +static inline real Min(creal a, creal b) { return (a < b) ? a : b; } -static inline real Max(creal a, creal b) -{ +static inline real Max(creal a, creal b) { return (a > b) ? a : b; } -static inline real Weight(creal sum, creal sqsum, cnumber n) -{ +static inline real Weight(creal sum, creal sqsum, cnumber n) { creal w = sqrt(sqsum*n); return (n - 1)/Max((w + sum)*(w - sum), NOTZERO); } @@ -235,5 +436,25 @@ static inline real Weight(creal sum, creal sqsum, cnumber n) /* abs(a) + (a == 0) */ #define Abs1(a) (((a) ^ NegQ(a)) - NegQ((a) - 1)) + +#ifdef MLVERSION + +static inline void Print(MLCONST char *s) +{ + MLPutFunction(stdlink, "EvaluatePacket", 1); + MLPutFunction(stdlink, "Print", 1); + MLPutString(stdlink, s); + MLEndPacket(stdlink); + + MLNextPacket(stdlink); + MLNewPacket(stdlink); +} + +#else + +#define Print(s) puts(s); fflush(stdout) + +#endif + #endif diff --git a/src/external/libCuba/src/cuba.h b/src/external/libCuba/src/cuba.h index b9f188a9..48b1c7ac 100644 --- a/src/external/libCuba/src/cuba.h +++ b/src/external/libCuba/src/cuba.h @@ -2,29 +2,9 @@ cuba.h Prototypes for the Cuba library this file is part of Cuba - last modified 16 Jun 10 th + last modified 30 Apr 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - #ifdef __cplusplus extern "C" { #endif @@ -68,6 +48,7 @@ void Suave(const int ndim, const int ncomp, const int flags, const int seed, const int mineval, const int maxeval, const int nnew, const double flatness, + const char *statefile, int *nregions, int *neval, int *fail, double integral[], double error[], double prob[]); @@ -77,6 +58,7 @@ void llSuave(const int ndim, const int ncomp, 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 char *statefile, int *nregions, long long int *neval, int *fail, double integral[], double error[], double prob[]); @@ -89,6 +71,7 @@ void Divonne(const int ndim, const int ncomp, const double border, const double maxchisq, const double mindeviation, const int ngiven, const int ldxgiven, double xgiven[], const int nextra, peakfinder_t peakfinder, + const char *statefile, int *nregions, int *neval, int *fail, double integral[], double error[], double prob[]); @@ -102,6 +85,7 @@ void llDivonne(const int ndim, const int ncomp, const long long int ngiven, const int ldxgiven, double xgiven[], const long long int nextra, void (*peakfinder)(const int *, const double [], int *, double []), + const char *statefile, int *nregions, long long int *neval, int *fail, double integral[], double error[], double prob[]); @@ -110,6 +94,7 @@ void Cuhre(const int ndim, const int ncomp, const double epsrel, const double epsabs, const int flags, const int mineval, const int maxeval, const int key, + const char *statefile, int *nregions, int *neval, int *fail, double integral[], double error[], double prob[]); @@ -119,9 +104,15 @@ void llCuhre(const int ndim, const int ncomp, const int flags, const long long int mineval, const long long int maxeval, const int key, + const char *statefile, int *nregions, long long int *neval, int *fail, double integral[], double error[], double prob[]); +void cubasetinit(void (*)(), void *); +void cubasetexit(void (*)(), void *); +void cubaruninit(void); +void cubaruninit(void); + #ifdef __cplusplus } #endif diff --git a/src/external/libCuba/src/cuhre/Cuhre.c b/src/external/libCuba/src/cuhre/Cuhre.c index c2194de1..33e17f62 100644 --- a/src/external/libCuba/src/cuhre/Cuhre.c +++ b/src/external/libCuba/src/cuhre/Cuhre.c @@ -2,56 +2,23 @@ Cuhre.c Adaptive integration using cubature rules by Thomas Hahn - last modified 16 Jun 10 th + last modified 17 Sep 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ +#define CUHRE +#define ROUTINE "Cuhre" #include "decl.h" - -#define Print(s) puts(s); fflush(stdout) +#include "CSample.c" /*********************************************************************/ -static inline void DoSample(This *t, count n, creal *x, real *f) -{ - t->neval += n; - while( n-- ) { - if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata) == ABORT ) - longjmp(t->abort, 1); - x += t->ndim; - f += t->ncomp; - } -} - -/*********************************************************************/ - -#include "common.c" - Extern void EXPORT(Cuhre)(ccount ndim, ccount ncomp, Integrand integrand, void *userdata, creal epsrel, creal epsabs, cint flags, cnumber mineval, cnumber maxeval, - ccount key, + ccount key, cchar *statefile, count *pnregions, number *pneval, int *pfail, real *integral, real *error, real *prob) { @@ -66,8 +33,7 @@ Extern void EXPORT(Cuhre)(ccount ndim, ccount ncomp, t.mineval = mineval; t.maxeval = maxeval; t.key = key; - t.nregions = 0; - t.neval = 0; + t.statefile = statefile; *pfail = Integrate(&t, integral, error, prob); *pnregions = t.nregions; @@ -80,9 +46,9 @@ Extern void EXPORT(cuhre)(ccount *pndim, ccount *pncomp, Integrand integrand, void *userdata, creal *pepsrel, creal *pepsabs, cint *pflags, cnumber *pmineval, cnumber *pmaxeval, - ccount *pkey, + ccount *pkey, cchar *statefile, count *pnregions, number *pneval, int *pfail, - real *integral, real *error, real *prob) + real *integral, real *error, real *prob, cint statefilelen) { This t; t.ndim = *pndim; @@ -95,11 +61,9 @@ Extern void EXPORT(cuhre)(ccount *pndim, ccount *pncomp, t.mineval = *pmineval; t.maxeval = *pmaxeval; t.key = *pkey; - t.nregions = 0; - t.neval = 0; - + CString(t.statefile, statefile, statefilelen); + *pfail = Integrate(&t, integral, error, prob); *pnregions = t.nregions; *pneval = t.neval; } - diff --git a/src/external/libCuba/src/cuhre/Integrate.c b/src/external/libCuba/src/cuhre/Integrate.c index 38d66806..25c40225 100644 --- a/src/external/libCuba/src/cuhre/Integrate.c +++ b/src/external/libCuba/src/cuhre/Integrate.c @@ -2,158 +2,174 @@ Integrate.c integrate over the unit hypercube this file is part of Cuhre - last modified 7 Jun 10 th + checkpointing by B. Chokoufe + last modified 17 Sep 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - #define POOLSIZE 1024 +typedef struct pool { + struct pool *next; + char region[]; +} Pool; + +typedef struct { + signature_t signature; + count nregions, ncur; + number neval; + Totals totals[]; +} State; + static int Integrate(This *t, real *integral, real *error, real *prob) { - TYPEDEFREGION; - typedef struct pool { - struct pool *next; - Region region[POOLSIZE]; - } Pool; + StateDecl; + csize_t statesize = sizeof(State) + NCOMP*sizeof(Totals); + Sized(State, state, statesize); + csize_t regionsize = RegionSize; + csize_t poolsize = sizeof(Pool) + POOLSIZE*regionsize; + Vector(Result, result, NCOMP); + Vector(char, out, 128*NCOMP + 256); - count dim, comp, ncur, ipool, npool; - int fail = -99; - Totals totals[NCOMP]; + Totals *tot, *Tot = state->totals + t->ncomp; + Result *res, *resL, *resR; + Bounds *b, *B; Pool *cur = NULL, *pool; Region *region; + count comp, ipool, npool; + int fail; if( VERBOSE > 1 ) { - char s[256]; - sprintf(s, "Cuhre input parameters:\n" + sprintf(out, "Cuhre input parameters:\n" " ndim " COUNT "\n ncomp " COUNT "\n" " epsrel " REAL "\n epsabs " REAL "\n" " flags %d\n mineval " NUMBER "\n maxeval " NUMBER "\n" - " key " COUNT, + " key " COUNT "\n" + " statefile \"%s\"", t->ndim, t->ncomp, t->epsrel, t->epsabs, t->flags, t->mineval, t->maxeval, - t->key); - Print(s); + t->key, + t->statefile); + Print(out); } if( BadComponent(t) ) return -2; if( BadDimension(t) ) return -1; - RuleAlloc(t); t->epsabs = Max(t->epsabs, NOTZERO); + + RuleAlloc(t); t->mineval = IMax(t->mineval, t->rule.n + 1); + FrameAlloc(t, ShmRm(t)); + ForkCores(t); - if( setjmp(t->abort) ) goto abort; + if( (fail = setjmp(t->abort)) ) goto abort; - Alloc(cur, 1); - cur->next = NULL; - ncur = 1; + StateSetup(t); - region = cur->region; - region->div = 0; - for( dim = 0; dim < t->ndim; ++dim ) { - Bounds *b = ®ion->bounds[dim]; - b->lower = 0; - b->upper = 1; + if( StateReadTest(t) ) { + StateReadOpen(t, fd) { + Pool *prev = NULL; + int size; + if( read(fd, state, statesize) != statesize || + state->signature != StateSignature(t, 4) ) break; + t->neval = state->neval; + t->nregions = state->nregions; + do { + MemAlloc(cur, poolsize); + cur->next = prev; + prev = cur; + size = read(fd, cur, poolsize); + } while( size == poolsize ); + if( size != state->ncur*regionsize ) break; + } StateReadClose(t, fd); } - Sample(t, region); + if( ini ) { + MemAlloc(cur, poolsize); + cur->next = NULL; + state->ncur = t->nregions = 1; - for( comp = 0; comp < t->ncomp; ++comp ) { - Totals *tot = &totals[comp]; - Result *r = ®ion->result[comp]; - tot->avg = tot->lastavg = tot->guess = r->avg; - tot->err = tot->lasterr = r->err; - tot->weightsum = 1/Max(Sq(r->err), NOTZERO); - tot->avgsum = tot->weightsum*r->avg; - tot->chisq = tot->chisqsum = tot->chisum = 0; + region = (Region *)cur->region; + region->div = 0; + for( B = (b = region->bounds) + t->ndim; b < B; ++b ) { + b->lower = 0; + b->upper = 1; + } + + t->neval = 0; + Sample(t, region); + + for( res = RegionResult(region), tot = state->totals; + tot < Tot; ++res, ++tot ) { + tot->avg = tot->lastavg = tot->guess = res->avg; + tot->err = tot->lasterr = res->err; + tot->weightsum = 1/Max(Sq(res->err), NOTZERO); + tot->avgsum = tot->weightsum*res->avg; + tot->chisq = tot->chisqsum = tot->chisum = 0; + } } - for( t->nregions = 1; ; ++t->nregions ) { + /* main iteration loop */ + for( ; ; ) { count maxcomp, bisectdim; real maxratio, maxerr; - Result result[NCOMP]; Region *regionL, *regionR; Bounds *bL, *bR; if( VERBOSE ) { - char s[128 + 128*NCOMP], *p = s; - - p += sprintf(p, "\n" + char *oe = out + sprintf(out, "\n" "Iteration " COUNT ": " NUMBER " integrand evaluations so far", t->nregions, t->neval); - - for( comp = 0; comp < t->ncomp; ++comp ) { - cTotals *tot = &totals[comp]; - p += sprintf(p, "\n[" COUNT "] " + for( tot = state->totals, comp = 0; tot < Tot; ++tot ) + oe += sprintf(oe, "\n[" COUNT "] " REAL " +- " REAL " \tchisq " REAL " (" COUNT " df)", - comp + 1, tot->avg, tot->err, tot->chisq, t->nregions - 1); - } - - Print(s); + ++comp, tot->avg, tot->err, tot->chisq, t->nregions - 1); + Print(out); } maxratio = -INFTY; maxcomp = 0; - for( comp = 0; comp < t->ncomp; ++comp ) { - creal ratio = totals[comp].err/MaxErr(totals[comp].avg); + for( tot = state->totals, comp = 0; tot < Tot; ++tot, ++comp ) { + creal ratio = tot->err/MaxErr(tot->avg); if( ratio > maxratio ) { maxratio = ratio; maxcomp = comp; } } - if( maxratio <= 1 && t->neval >= t->mineval ) { - fail = 0; + if( maxratio <= 1 && t->neval >= t->mineval ) break; + + if( t->neval >= t->maxeval ) { + fail = 1; break; } - if( t->neval >= t->maxeval ) break; - maxerr = -INFTY; - regionL = cur->region; - npool = ncur; + regionL = (Region *)cur->region; + npool = state->ncur; for( pool = cur; pool; npool = POOLSIZE, pool = pool->next ) for( ipool = 0; ipool < npool; ++ipool ) { - Region *region = &pool->region[ipool]; - creal err = region->result[maxcomp].err; + Region *region = RegionPtr(pool, ipool); + creal err = RegionResult(region)[maxcomp].err; if( err > maxerr ) { maxerr = err; regionL = region; } } - if( ncur == POOLSIZE ) { + if( state->ncur == POOLSIZE ) { Pool *prev = cur; - Alloc(cur, 1); + MemAlloc(cur, poolsize); cur->next = prev; - ncur = 0; + state->ncur = 0; } - regionR = &cur->region[ncur++]; + regionR = RegionPtr(cur, state->ncur++); regionR->div = ++regionL->div; - ResCopy(result, regionL->result); - VecCopy(regionR->bounds, regionL->bounds); + FCopy(result, RegionResult(regionL)); + XCopy(regionR->bounds, regionL->bounds); bisectdim = result[maxcomp].bisectdim; bL = ®ionL->bounds[bisectdim]; @@ -163,25 +179,25 @@ static int Integrate(This *t, real *integral, real *error, real *prob) Sample(t, regionL); Sample(t, regionR); - for( comp = 0; comp < t->ncomp; ++comp ) { - cResult *r = &result[comp]; - Result *rL = ®ionL->result[comp]; - Result *rR = ®ionR->result[comp]; - Totals *tot = &totals[comp]; + for( res = result, + resL = RegionResult(regionL), + resR = RegionResult(regionR), + tot = state->totals; + tot < Tot; ++res, ++resL, ++resR, ++tot ) { real diff, err, w, avg, sigsq; - tot->lastavg += diff = rL->avg + rR->avg - r->avg; + tot->lastavg += diff = resL->avg + resR->avg - res->avg; diff = fabs(.25*diff); - err = rL->err + rR->err; + err = resL->err + resR->err; if( err > 0 ) { creal c = 1 + 2*diff/err; - rL->err *= c; - rR->err *= c; + resL->err *= c; + resR->err *= c; } - rL->err += diff; - rR->err += diff; - tot->lasterr += rL->err + rR->err - r->err; + resL->err += diff; + resR->err += diff; + tot->lasterr += resL->err + resR->err - res->err; tot->weightsum += w = 1/Max(Sq(tot->lasterr), NOTZERO); sigsq = 1/tot->weightsum; @@ -200,10 +216,22 @@ static int Integrate(This *t, real *integral, real *error, real *prob) tot->err = sqrt(sigsq); } } + ++t->nregions; + + if( StateWriteTest(t) ) { + StateWriteOpen(t, fd) { + Pool *prev = cur; + state->signature = StateSignature(t, 4); + state->nregions = t->nregions; + state->neval = t->neval; + StateWrite(fd, state, statesize); + while( (prev = prev->next) ) StateWrite(fd, prev, poolsize); + StateWrite(fd, cur, state->ncur*regionsize); + } StateWriteClose(t, fd); + } } - for( comp = 0; comp < t->ncomp; ++comp ) { - cTotals *tot = &totals[comp]; + for( tot = state->totals, comp = 0; tot < Tot; ++tot, ++comp ) { integral[comp] = tot->avg; error[comp] = tot->err; prob[comp] = ChiSquare(tot->chisq, t->nregions - 1); @@ -214,27 +242,20 @@ static int Integrate(This *t, real *integral, real *error, real *prob) MLPutFunction(stdlink, "List", 2); MLPutFunction(stdlink, "List", t->nregions); - npool = ncur; + npool = state->ncur; for( pool = cur; pool; npool = POOLSIZE, pool = pool->next ) for( ipool = 0; ipool < npool; ++ipool ) { - Region const *region = &pool->region[ipool]; - real lower[NDIM], upper[NDIM]; + Region const *region = RegionPtr(pool, ipool); + Result *Res; - for( dim = 0; dim < t->ndim; ++dim ) { - cBounds *b = ®ion->bounds[dim]; - lower[dim] = b->lower; - upper[dim] = b->upper; - } - - MLPutFunction(stdlink, "Cuba`Cuhre`region", 3); - MLPutRealList(stdlink, lower, t->ndim); - MLPutRealList(stdlink, upper, t->ndim); + MLPutFunction(stdlink, "Cuba`Cuhre`region", 2); + MLPutRealList(stdlink, (real *)region->bounds, 2*t->ndim); MLPutFunction(stdlink, "List", t->ncomp); - for( comp = 0; comp < t->ncomp; ++comp ) { - cResult *r = ®ion->result[comp]; - real res[] = {r->avg, r->err}; - MLPutRealList(stdlink, res, Elements(res)); + for( Res = (res = RegionResult(region)) + t->ncomp; + res < Res; ++res ) { + real r[] = {res->avg, res->err}; + MLPutRealList(stdlink, r, Elements(r)); } } } @@ -246,8 +267,12 @@ abort: free(pool); } + WaitCores(t); + FrameFree(t); RuleFree(t); + StateRemove(t); + return fail; } diff --git a/src/external/libCuba/src/cuhre/Makefile.am b/src/external/libCuba/src/cuhre/Makefile.am new file mode 100644 index 00000000..5ab6b4eb --- /dev/null +++ b/src/external/libCuba/src/cuhre/Makefile.am @@ -0,0 +1,13 @@ +## Process this file with automake to create Makefile.in +## $Id: Makefile.am 4808 2011-03-21 14:33:34Z l_wojek $ + +c_sources = Cuhre.c + +AM_CPPFLAGS = -I. -I.. -I../common -DNOUNDERSCORE +AM_CFLAGS = $(LOCAL_CUBA_LIB_CFLAGS) +AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS) + +noinst_LTLIBRARIES = libcuhre.la + +libcuhre_la_SOURCES = $(c_sources) +libcuhre_la_LDFLAGS = $(AM_LDFLAGS) \ No newline at end of file diff --git a/src/external/libCuba/src/cuhre/Rule.c b/src/external/libCuba/src/cuhre/Rule.c index fdc23c04..e3e931f3 100644 --- a/src/external/libCuba/src/cuhre/Rule.c +++ b/src/external/libCuba/src/cuhre/Rule.c @@ -3,38 +3,12 @@ integration with cubature rules code lifted with minor modifications from DCUHRE by J. Berntsen, T. Espelid, and A. Genz - this file is part of Divonne - last modified 8 Jun 10 th + this file is part of Cuhre + last modified 5 Aug 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ -enum { nrules = 5 }; - -#define TYPEDEFSET \ - typedef struct { \ - count n; \ - real weight[nrules], scale[nrules], norm[nrules]; \ - real gen[NDIM]; \ - } Set +#define NextSet(p) p = (Set *)((char *)p + setsize) /*********************************************************************/ @@ -83,86 +57,84 @@ static void Rule13Alloc(This *t) enum { nsets = 14, ndim = 2 }; - TYPEDEFSET; - count n, r; Set *first, *last, *s, *x; + csize_t setsize = SetSize; - Alloc(first, nsets); - Clear(first, nsets); + Die(first = calloc(nsets, setsize)); last = first; n = last->n = 1; Copy(last->weight, w[0], nrules); - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[1], nrules); last->gen[0] = g[0]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[2], nrules); last->gen[0] = g[1]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[3], nrules); last->gen[0] = g[2]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[4], nrules); last->gen[0] = g[3]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[5], nrules); last->gen[0] = g[4]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); Copy(last->weight, w[6], nrules); last->gen[0] = g[5]; last->gen[1] = g[5]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); Copy(last->weight, w[7], nrules); last->gen[0] = g[6]; last->gen[1] = g[6]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); Copy(last->weight, w[8], nrules); last->gen[0] = g[7]; last->gen[1] = g[7]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); Copy(last->weight, w[9], nrules); last->gen[0] = g[8]; last->gen[1] = g[8]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); Copy(last->weight, w[10], nrules); last->gen[0] = g[9]; last->gen[1] = g[9]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1); Copy(last->weight, w[11], nrules); last->gen[0] = g[10]; last->gen[1] = g[11]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1); Copy(last->weight, w[12], nrules); last->gen[0] = g[12]; last->gen[1] = g[13]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1); Copy(last->weight, w[13], nrules); last->gen[0] = g[14]; @@ -175,12 +147,12 @@ static void Rule13Alloc(This *t) t->rule.errcoeff[2] = 5; t->rule.n = n; - for( s = first; s <= last; ++s ) + for( s = first; s <= last; NextSet(s) ) for( r = 1; r < nrules - 1; ++r ) { creal scale = (s->weight[r] == 0) ? 100 : -s->weight[r + 1]/s->weight[r]; real sum = 0; - for( x = first; x <= last; ++x ) + for( x = first; x <= last; NextSet(x) ) sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]); s->scale[r] = scale; s->norm[r] = 1/sum; @@ -231,84 +203,82 @@ static void Rule11Alloc(This *t) enum { nsets = 13, ndim = 3 }; - TYPEDEFSET; - count n, r; Set *first, *last, *s, *x; + csize_t setsize = SetSize; - Alloc(first, nsets); - Clear(first, nsets); + Die(first = calloc(nsets, setsize)); last = first; n = last->n = 1; Copy(last->weight, w[0], nrules); - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[1], nrules); last->gen[0] = g[0]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[2], nrules); last->gen[0] = g[1]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[3], nrules); last->gen[0] = g[2]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[4], nrules); last->gen[0] = g[3]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[5], nrules); last->gen[0] = g[4]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); Copy(last->weight, w[6], nrules); last->gen[0] = g[5]; last->gen[1] = g[5]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); Copy(last->weight, w[7], nrules); last->gen[0] = g[6]; last->gen[1] = g[6]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3; Copy(last->weight, w[8], nrules); last->gen[0] = g[7]; last->gen[1] = g[7]; last->gen[2] = g[7]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3; Copy(last->weight, w[9], nrules); last->gen[0] = g[8]; last->gen[1] = g[8]; last->gen[2] = g[8]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3; Copy(last->weight, w[10], nrules); last->gen[0] = g[9]; last->gen[1] = g[9]; last->gen[2] = g[9]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1)*(ndim - 2); Copy(last->weight, w[11], nrules); last->gen[0] = g[10]; last->gen[1] = g[11]; last->gen[2] = g[11]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1)*(ndim - 2); Copy(last->weight, w[12], nrules); last->gen[0] = g[12]; @@ -322,12 +292,12 @@ static void Rule11Alloc(This *t) t->rule.errcoeff[2] = 3; t->rule.n = n; - for( s = first; s <= last; ++s ) + for( s = first; s <= last; NextSet(s) ) for( r = 1; r < nrules - 1; ++r ) { creal scale = (s->weight[r] == 0) ? 100 : -s->weight[r + 1]/s->weight[r]; real sum = 0; - for( x = first; x <= last; ++x ) + for( x = first; x <= last; NextSet(x) ) sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]); s->scale[r] = scale; s->norm[r] = 1/sum; @@ -368,15 +338,13 @@ static void Rule9Alloc(This *t) enum { nsets = 9 }; - TYPEDEFSET; - ccount ndim = t->ndim; ccount twondim = 1 << ndim; count dim, n, r; Set *first, *last, *s, *x; + csize_t setsize = SetSize; - Alloc(first, nsets); - Clear(first, nsets); + Die(first = calloc(nsets, setsize)); last = first; n = last->n = 1; @@ -386,7 +354,7 @@ static void Rule9Alloc(This *t) last->weight[3] = ndim*(ndim*w[9] + w[10]) - 1 + last->weight[0]; last->weight[4] = ndim*w[11] + 1 - last->weight[0]; - ++last; + NextSet(last); n += last->n = 2*ndim; last->weight[0] = ndim*(ndim*w[12] + w[13]) + w[14]; last->weight[1] = ndim*(ndim*w[15] + w[16]) + w[17]; @@ -395,7 +363,7 @@ static void Rule9Alloc(This *t) last->weight[4] = w[21] - last->weight[0]; last->gen[0] = g[0]; - ++last; + NextSet(last); n += last->n = 2*ndim; last->weight[0] = ndim*w[22] + w[23]; last->weight[1] = ndim*w[24] + w[25]; @@ -404,7 +372,7 @@ static void Rule9Alloc(This *t) last->weight[4] = -last->weight[0]; last->gen[0] = g[1]; - ++last; + NextSet(last); n += last->n = 2*ndim; last->weight[0] = w[29]; last->weight[1] = w[30]; @@ -413,12 +381,12 @@ static void Rule9Alloc(This *t) last->weight[4] = -w[29]; last->gen[0] = g[2]; - ++last; + NextSet(last); n += last->n = 2*ndim; last->weight[2] = w[32]; last->gen[0] = g[3]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); last->weight[0] = w[33] - ndim*w[12]; last->weight[1] = w[34] - ndim*w[15]; @@ -428,7 +396,7 @@ static void Rule9Alloc(This *t) last->gen[0] = g[0]; last->gen[1] = g[0]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1); last->weight[0] = w[36]; last->weight[1] = w[37]; @@ -438,7 +406,7 @@ static void Rule9Alloc(This *t) last->gen[0] = g[0]; last->gen[1] = g[1]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3; last->weight[0] = w[39]; last->weight[1] = w[40]; @@ -449,7 +417,7 @@ static void Rule9Alloc(This *t) last->gen[1] = g[0]; last->gen[2] = g[0]; - ++last; + NextSet(last); n += last->n = twondim; last->weight[0] = w[41]/twondim; last->weight[1] = w[7]/twondim; @@ -466,12 +434,12 @@ static void Rule9Alloc(This *t) t->rule.errcoeff[2] = 5; t->rule.n = n; - for( s = first; s <= last; ++s ) + for( s = first; s <= last; NextSet(s) ) for( r = 1; r < nrules - 1; ++r ) { creal scale = (s->weight[r] == 0) ? 100 : -s->weight[r + 1]/s->weight[r]; real sum = 0; - for( x = first; x <= last; ++x ) + for( x = first; x <= last; NextSet(x) ) sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]); s->scale[r] = scale; s->norm[r] = 1/sum; @@ -501,15 +469,13 @@ static void Rule7Alloc(This *t) enum { nsets = 6 }; - TYPEDEFSET; - ccount ndim = t->ndim; ccount twondim = 1 << ndim; count dim, n, r; Set *first, *last, *s, *x; + csize_t setsize = SetSize; - Alloc(first, nsets); - Clear(first, nsets); + Die(first = calloc(nsets, setsize)); last = first; n = last->n = 1; @@ -519,7 +485,7 @@ static void Rule7Alloc(This *t) last->weight[3] = ndim*(ndim*w[7] + w[8]) - w[9]; last->weight[4] = 1 - last->weight[0]; - ++last; + NextSet(last); n += last->n = 2*ndim; last->weight[0] = w[10]; last->weight[1] = w[11]; @@ -528,7 +494,7 @@ static void Rule7Alloc(This *t) last->weight[4] = -w[10]; last->gen[0] = g[1]; - ++last; + NextSet(last); n += last->n = 2*ndim; last->weight[0] = w[13] - ndim*w[0]; last->weight[1] = w[14] - ndim*w[3]; @@ -537,12 +503,12 @@ static void Rule7Alloc(This *t) last->weight[4] = -last->weight[0]; last->gen[0] = g[0]; - ++last; + NextSet(last); n += last->n = 2*ndim; last->weight[2] = w[17]; last->gen[0] = g[2]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); last->weight[0] = -w[7]; last->weight[1] = w[18]; @@ -552,7 +518,7 @@ static void Rule7Alloc(This *t) last->gen[0] = g[0]; last->gen[1] = g[0]; - ++last; + NextSet(last); n += last->n = twondim; last->weight[0] = w[20]/twondim; last->weight[1] = w[5]/twondim; @@ -569,12 +535,12 @@ static void Rule7Alloc(This *t) t->rule.errcoeff[2] = 5; t->rule.n = n; - for( s = first; s <= last; ++s ) + for( s = first; s <= last; NextSet(s) ) for( r = 1; r < nrules - 1; ++r ) { creal scale = (s->weight[r] == 0) ? 100 : -s->weight[r + 1]/s->weight[r]; real sum = 0; - for( x = first; x <= last; ++x ) + for( x = first; x <= last; NextSet(x) ) sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]); s->scale[r] = scale; s->norm[r] = 1/sum; @@ -594,15 +560,12 @@ static inline void RuleAlloc(This *t) else if( t->ndim == 3 ) Rule11Alloc(t); else Rule9Alloc(t); } - Alloc(t->rule.x, t->rule.n*(t->ndim + t->ncomp)); - t->rule.f = t->rule.x + t->rule.n*t->ndim; } /*********************************************************************/ static inline void RuleFree(cThis *t) { - free(t->rule.x); free(t->rule.first); } @@ -657,25 +620,23 @@ next: /*********************************************************************/ -static void Sample(This *t, void *voidregion) +static void Sample(This *t, Region *region) { - TYPEDEFREGION; - TYPEDEFSET; - - Region *const region = (Region *)voidregion; + csize_t setsize = SetSize; creal vol = ldexp(1., -region->div); - real *x = t->rule.x, *f = t->rule.f; - Set *first = (Set *)t->rule.first, *last = (Set *)t->rule.last, *s; + real *x = t->frame, *f = x + t->rule.n*t->ndim; + Set *first = t->rule.first, *last = t->rule.last, *s; + 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]); ccount offset = 2*t->ndim*t->ncomp; - count dim, comp, rul, n, maxdim = 0; + count dim, rul, n, maxdim = 0; real maxrange = 0; - for( dim = 0; dim < t->ndim; ++dim ) { - cBounds *b = ®ion->bounds[dim]; + for( b = region->bounds, dim = 0; b < B; ++b, ++dim ) { creal range = b->upper - b->lower; if( range > maxrange ) { maxrange = range; @@ -683,13 +644,12 @@ static void Sample(This *t, void *voidregion) } } - for( s = first; s <= last; ++s ) + for( s = first; s <= last; NextSet(s) ) if( s->n ) x = ExpandFS(t, region->bounds, s->gen, x); - DoSample(t, t->rule.n, t->rule.x, f); + DoSample(t, t->rule.n, t->frame, f); - for( comp = 0; comp < t->ncomp; ++comp ) { - Result *r = ®ion->result[comp]; + for( res = result; res < Res; ++res ) { real sum[nrules]; creal *f1 = f; creal base = *f1*2*(1 - ratio); @@ -707,11 +667,11 @@ static void Sample(This *t, void *voidregion) bisectdim = dim; } } - r->bisectdim = bisectdim; + res->bisectdim = bisectdim; f1 = f++; Zap(sum); - for( s = first; s <= last; ++s ) + for( s = first; s <= last; NextSet(s) ) for( n = s->n; n; --n ) { creal fun = *f1; f1 += t->ncomp; @@ -726,37 +686,35 @@ static void Sample(This *t, void *voidregion) for( rul = 1; rul < nrules - 1; ++rul ) { real maxerr = 0; - for( s = first; s <= last; ++s ) + for( s = first; s <= last; NextSet(s) ) maxerr = Max(maxerr, fabs(sum[rul + 1] + s->scale[rul]*sum[rul])*s->norm[rul]); sum[rul] = maxerr; } - r->avg = vol*sum[0]; - r->err = vol*( + res->avg = vol*sum[0]; + res->err = vol*( (errcoeff[0]*sum[1] <= sum[2] && errcoeff[0]*sum[2] <= sum[3]) ? errcoeff[1]*sum[1] : errcoeff[2]*Max(Max(sum[1], sum[2]), sum[3]) ); } if( VERBOSE > 2 ) { - char s[64*NDIM + 128*NCOMP], *p = s; + Vector(char, out, 64*NDIM + 128*NCOMP); + char *oe = out; + count comp; + cchar *msg = "\nRegion (" REALF ") - (" REALF ")"; - for( dim = 0; dim < t->ndim; ++dim ) { - cBounds *b = ®ion->bounds[dim]; - p += sprintf(p, - (dim == 0) ? "\nRegion (" REALF ") - (" REALF ")" : - "\n (" REALF ") - (" REALF ")", - b->lower, b->upper); + for( b = region->bounds; b < B; ++b ) { + oe += sprintf(oe, msg, b->lower, b->upper); + msg = "\n (" REALF ") - (" REALF ")"; } - for( comp = 0; comp < t->ncomp; ++comp ) { - cResult *r = ®ion->result[comp]; - p += sprintf(p, "\n[" COUNT "] " - REAL " +- " REAL, comp + 1, r->avg, r->err); - } + for( res = result, comp = 0; res < Res; ++res ) + oe += sprintf(oe, "\n[" COUNT "] " + REAL " +- " REAL, ++comp, res->avg, res->err); - Print(s); + Print(out); } } diff --git a/src/external/libCuba/src/cuhre/common.c b/src/external/libCuba/src/cuhre/common.c index 4dd8741f..fcd4403f 100644 --- a/src/external/libCuba/src/cuhre/common.c +++ b/src/external/libCuba/src/cuhre/common.c @@ -2,44 +2,22 @@ common.c includes most of the modules this file is part of Cuhre - last modified 7 Jun 10 th + last modified 2 Aug 13 11 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - #include "ChiSquare.c" #include "Rule.c" static inline bool BadDimension(cThis *t) { - if( t->ndim > NDIM ) return true; + if( t->ndim > MAXDIM ) return true; return t->ndim < 2; } static inline bool BadComponent(cThis *t) { - if( t->ncomp > NCOMP ) return true; + if( t->ncomp > MAXCOMP ) return true; return t->ncomp < 1; } -#include "Integrate.c" - diff --git a/src/external/libCuba/src/cuhre/decl.h b/src/external/libCuba/src/cuhre/decl.h index a0512f15..c846e9d2 100644 --- a/src/external/libCuba/src/cuhre/decl.h +++ b/src/external/libCuba/src/cuhre/decl.h @@ -2,27 +2,9 @@ decl.h Type declarations this file is part of Cuhre - last modified 7 Jun 10 th */ + last modified 26 Jul 13 th +*/ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ #include "stddecl.h" @@ -47,9 +29,18 @@ typedef struct { typedef const Bounds cBounds; +enum { nrules = 5 }; + typedef struct { - real *x, *f; - void *first, *last; + count n; + real weight[nrules], scale[nrules], norm[nrules]; + real gen[]; +} Set; + +#define SetSize (sizeof(Set) + t->ndim*sizeof(real)) + +typedef struct { + Set *first, *last; real errcoeff[3]; count n; } Rule; @@ -63,22 +54,34 @@ typedef struct _this { #ifndef MLVERSION Integrand integrand; void *userdata; +#ifdef HAVE_FORK + int ncores, *child; + SHM_ONLY(int shmid;) #endif +#endif + real *frame; real epsrel, epsabs; int flags; number mineval, maxeval; count key, nregions; + cchar *statefile; number neval; Rule rule; jmp_buf abort; } This; +#define nframe rule.n + typedef const This cThis; -#define TYPEDEFREGION \ - typedef struct region { \ - count div; \ - Result result[NCOMP]; \ - Bounds bounds[NDIM]; \ - } Region +typedef struct region { + count div; + Bounds bounds[]; +} Region; + +#define RegionSize (sizeof(Region) + t->ndim*sizeof(Bounds) + t->ncomp*sizeof(Result)) + +#define RegionResult(r) ((Result *)(r->bounds + t->ndim)) + +#define RegionPtr(p, n) ((Region *)((char *)p->region + (n)*regionsize)) diff --git a/src/external/libCuba/src/cuhre/util.c b/src/external/libCuba/src/cuhre/util.c deleted file mode 100644 index c3dd291d..00000000 --- a/src/external/libCuba/src/cuhre/util.c +++ /dev/null @@ -1,37 +0,0 @@ -/* - util.c - Utility functions - this file is part of Cuhre - last modified 9 Jan 05 th -*/ - -/*************************************************************************** - * Copyright (C) 2004-2009 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - -#include "decl.h" - -static count ndim_, ncomp_, nregions_; -static number neval_; - - -#ifdef DEBUG -#include "debug.c" -#endif - diff --git a/src/external/libCuba/src/divonne/Divonne.c b/src/external/libCuba/src/divonne/Divonne.c index abd5db09..a0577cc6 100644 --- a/src/external/libCuba/src/divonne/Divonne.c +++ b/src/external/libCuba/src/divonne/Divonne.c @@ -4,83 +4,17 @@ originally by J.H. Friedman and M.H. Wright (CERNLIB subroutine D151) this version by Thomas Hahn - last modified 16 Jun 10 th + last modified 17 Sep 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ +#define DIVONNE +#define ROUTINE "Divonne" #include "decl.h" - -#define Print(s) puts(s); fflush(stdout) +#include "CSample.c" /*********************************************************************/ -static inline void DoSample(This *t, number n, ccount ldx, creal *x, real *f) -{ - t->neval += n; - while( n-- ) { - if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata, &t->phase) == ABORT ) - longjmp(t->abort, 1); - x += ldx; - f += t->ncomp; - } -} - -/*********************************************************************/ - -static inline count SampleExtra(This *t, cBounds *b) -{ - number n = t->nextra; - t->peakfinder(&t->ndim, b, &n, t->xextra); - DoSample(t, n, t->ldxgiven, t->xextra, t->fextra); - return n; -} - -/*********************************************************************/ - -static inline void AllocGiven(This *t, creal *xgiven) -{ - if( t->ngiven | t->nextra ) { - cnumber nxgiven = t->ngiven*(t->ldxgiven = IMax(t->ldxgiven, t->ndim)); - cnumber nxextra = t->nextra*t->ldxgiven; - cnumber nfgiven = t->ngiven*t->ncomp; - cnumber nfextra = t->nextra*t->ncomp; - - Alloc(t->xgiven, nxgiven + nxextra + nfgiven + nfextra); - t->xextra = t->xgiven + nxgiven; - t->fgiven = t->xextra + nxextra; - t->fextra = t->fgiven + nfgiven; - - if( nxgiven ) { - t->phase = 0; - Copy(t->xgiven, xgiven, nxgiven); - DoSample(t, t->ngiven, t->ldxgiven, t->xgiven, t->fgiven); - } - } -} - -/*********************************************************************/ - -#include "common.c" - Extern void EXPORT(Divonne)(ccount ndim, ccount ncomp, Integrand integrand, void *userdata, creal epsrel, creal epsabs, @@ -88,8 +22,9 @@ Extern void EXPORT(Divonne)(ccount ndim, ccount ncomp, cnumber mineval, cnumber maxeval, cint key1, cint key2, cint key3, ccount maxpass, creal border, creal maxchisq, creal mindeviation, - cnumber ngiven, ccount ldxgiven, creal *xgiven, + cnumber ngiven, ccount ldxgiven, real *xgiven, cnumber nextra, PeakFinder peakfinder, + cchar *statefile, int *pnregions, number *pneval, int *pfail, real *integral, real *error, real *prob) { @@ -112,20 +47,15 @@ Extern void EXPORT(Divonne)(ccount ndim, ccount ncomp, t.maxchisq = maxchisq; t.mindeviation = mindeviation; t.ngiven = ngiven; - t.xgiven = NULL; + t.xgiven = xgiven; t.ldxgiven = ldxgiven; t.nextra = nextra; t.peakfinder = peakfinder; - t.nregions = 0; - t.neval = 0; - - AllocGiven(&t, xgiven); + t.statefile = statefile; *pfail = Integrate(&t, integral, error, prob); *pnregions = t.nregions; *pneval = t.neval; - - free(t.xgiven); } /*********************************************************************/ @@ -137,10 +67,11 @@ Extern void EXPORT(divonne)(ccount *pndim, ccount *pncomp, cnumber *pmineval, cnumber *pmaxeval, cint *pkey1, cint *pkey2, cint *pkey3, ccount *pmaxpass, creal *pborder, creal *pmaxchisq, creal *pmindeviation, - cnumber *pngiven, ccount *pldxgiven, creal *xgiven, + cnumber *pngiven, ccount *pldxgiven, real *xgiven, cnumber *pnextra, PeakFinder peakfinder, + cchar *statefile, int *pnregions, number *pneval, int *pfail, - real *integral, real *error, real *prob) + real *integral, real *error, real *prob, cint statefilelen) { This t; t.ndim = *pndim; @@ -161,19 +92,14 @@ Extern void EXPORT(divonne)(ccount *pndim, ccount *pncomp, t.maxchisq = *pmaxchisq; t.mindeviation = *pmindeviation; t.ngiven = *pngiven; - t.xgiven = NULL; + t.xgiven = xgiven; t.ldxgiven = *pldxgiven; t.nextra = *pnextra; t.peakfinder = peakfinder; - t.nregions = 0; - t.neval = 0; - - AllocGiven(&t, xgiven); + CString(t.statefile, statefile, statefilelen); *pfail = Integrate(&t, integral, error, prob); *pnregions = t.nregions; *pneval = t.neval; - - free(t.xgiven); } diff --git a/src/external/libCuba/src/divonne/Explore.c b/src/external/libCuba/src/divonne/Explore.c index 691e7986..4f27ac25 100644 --- a/src/external/libCuba/src/divonne/Explore.c +++ b/src/external/libCuba/src/divonne/Explore.c @@ -2,29 +2,9 @@ Explore.c sample region, determine min and max, split if necessary this file is part of Divonne - last modified 8 Jun 10 th + last modified 2 Aug 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - typedef struct { real fmin, fmax; @@ -33,38 +13,21 @@ typedef struct { /*********************************************************************/ -static bool Explore(This *t, count iregion, cSamples *samples, - cint depth, cint flags) +static int ExploreSerial(This *t, ccount iregion) { -#define SPLICE (flags & 1) -#define HAVESAMPLES (flags & 2) + csize_t regionsize = RegionSize; + Region *region = RegionPtr(iregion); + cBounds *bounds = region->bounds; + Result *result = RegionResult(region); - TYPEDEFREGION; - - count n, dim, comp, maxcomp; - Extrema extrema[NCOMP]; - Result *r; + Vector(Extrema, extrema, NCOMP); + Vector(real, xtmp, NDIM); + Result *r, *r0; creal *x; real *f; real halfvol, maxerr; - Region *region; - Bounds *bounds; - Result *result; - - /* needed as of gcc 3.3 to make gcc correctly address region #@$&! */ - sizeof(*region); - - if( SPLICE ) { - if( t->nregions == t->size ) { - t->size += CHUNKSIZE; - ReAlloc(t->voidregion, t->size*sizeof(Region)); - } - VecCopy(RegionPtr(t->nregions)->bounds, RegionPtr(iregion)->bounds); - iregion = t->nregions++; - } - region = RegionPtr(iregion); - bounds = region->bounds; - result = region->result; + count n, dim, comp, maxcomp; + cSamples *samples = &t->samples[region->isamples]; for( comp = 0; comp < t->ncomp; ++comp ) { Extrema *e = &extrema[comp]; @@ -73,7 +36,7 @@ static bool Explore(This *t, count iregion, cSamples *samples, e->xmin = e->xmax = NULL; } - if( !HAVESAMPLES ) { + if( region->isamples == 0 ) { /* others already sampled */ real vol = 1; for( dim = 0; dim < t->ndim; ++dim ) { cBounds *b = &bounds[dim]; @@ -108,7 +71,7 @@ skip: f += t->ncomp; } - samples->sampler(t, samples, bounds, vol); + samples->sampler(t, iregion); } x = samples->x; @@ -131,30 +94,27 @@ skip: for( comp = 0; comp < t->ncomp; ++comp ) { Extrema *e = &extrema[comp]; Result *r = &result[comp]; - real xtmp[NDIM], ftmp, err; + real ftmp, err; if( e->xmin ) { /* not all NaNs */ t->selectedcomp = comp; - VecCopy(xtmp, e->xmin); + XCopy(xtmp, e->xmin); ftmp = FindMinimum(t, bounds, xtmp, e->fmin); if( ftmp < r->fmin ) { r->fmin = ftmp; - VecCopy(r->xmin, xtmp); + XCopy(&r->xminmax[0], xtmp); } t->selectedcomp = Tag(comp); - VecCopy(xtmp, e->xmax); + XCopy(xtmp, e->xmax); ftmp = -FindMinimum(t, bounds, xtmp, -e->fmax); if( ftmp > r->fmax ) { r->fmax = ftmp; - VecCopy(r->xmax, xtmp); + XCopy(&r->xminmax[t->ndim], xtmp); } } - r->avg = samples->avg[comp]; - r->err = samples->err[comp]; r->spread = halfvol*(r->fmax - r->fmin); - err = r->spread/Max(fabs(r->avg), NOTZERO); if( err > maxerr ) { maxerr = err; @@ -166,26 +126,25 @@ skip: if( maxcomp == -1 ) { /* all NaNs */ region->depth = 0; - return false; + return -1; } region->cutcomp = maxcomp; - r = ®ion->result[maxcomp]; + r0 = RegionResult(region); + r = r0 + maxcomp; if( halfvol*(r->fmin + r->fmax) > r->avg ) { region->fminor = r->fmin; region->fmajor = r->fmax; - region->xmajor = r->xmax - (real *)region->result; + region->xmajor = &r->xminmax[t->ndim] - (real *)r0; } else { region->fminor = r->fmax; region->fmajor = r->fmin; - region->xmajor = r->xmin - (real *)region->result; + region->xmajor = &r->xminmax[0] - (real *)r0; } - region->depth = IDim(depth); - - if( !HAVESAMPLES ) { - if( samples->weight*r->spread < r->err || + if( region->isamples == 0 ) { + if( (region->depth < INIDEPTH && r->spread < samples->neff*r->err) || r->spread < t->totals[maxcomp].secondspread ) region->depth = 0; if( region->depth == 0 ) for( comp = 0; comp < t->ncomp; ++comp ) @@ -193,7 +152,8 @@ skip: Max(t->totals[comp].secondspread, result[comp].spread); } - if( region->depth ) Split(t, iregion, region->depth); - return true; + if( region->depth ) Split(t, iregion); + + return iregion; } diff --git a/src/external/libCuba/src/divonne/FindMinimum.c b/src/external/libCuba/src/divonne/FindMinimum.c index 968a5187..c7e00b2d 100644 --- a/src/external/libCuba/src/divonne/FindMinimum.c +++ b/src/external/libCuba/src/divonne/FindMinimum.c @@ -2,37 +2,17 @@ FindMinimum.c find minimum (maximum) of hyperrectangular region this file is part of Divonne - last modified 8 Jun 10 th + last modified 7 Aug 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ +#define EPS POW2(52) +#define RTEPS POW2(26) +#define QEPS POW2(13) -#define EPS 0x1p-52 -#define RTEPS 0x1p-26 -#define QEPS 0x1p-13 - -#define DELTA 0x1p-16 -#define RTDELTA 0x1p-8 -#define QDELTA 0x1p-4 +#define DELTA POW2(16) +#define RTDELTA POW2(8) +#define QDELTA POW2(4) /* #define DELTA 1e-5 @@ -162,7 +142,8 @@ static void UpdateCholesky(cThis *t, ccount n, real *hessian, static inline void BFGS(cThis *t, ccount n, real *hessian, creal *gnew, creal *g, real *p, creal dx) { - real y[NDIM], c; + Vector(real, y, NDIM); + real c; count i, j; for( i = 0; i < n; ++i ) @@ -210,7 +191,7 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree, real tol = ftol, tol2 = tol + tol; Point cur = {0, fini}; - VecCopy(x, xini); + XCopy(x, xini); /* don't even try if a) we'd walk backwards, @@ -360,9 +341,10 @@ static Point LineSearch(This *t, ccount nfree, ccount *ifree, static real LocalSearch(This *t, ccount nfree, ccount *ifree, cBounds *b, creal *x, creal fx, real *z) { + Vector(real, y, NDIM); + Vector(real, p, NDIM); real delta, smax, sopp, spmax, snmax; - real y[NDIM], fy, fz, ftest; - real p[NDIM]; + real fy, fz, ftest; int sign; count i; @@ -389,7 +371,7 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree, /* Move along p until the integrand changes appreciably or we come close to a border. */ - VecCopy(y, x); + XCopy(y, x); ftest = SUFTOL*(1 + fabs(fx)); delta = RTDELTA/5; do { @@ -442,7 +424,7 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree, /* Move along p' until the integrand changes appreciably or we come close to a border. */ - VecCopy(z, y); + XCopy(z, y); ftest = SUFTOL*(1 + fabs(fy)); delta = RTDELTA/5; do { @@ -468,7 +450,7 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree, grad = (fy - fz)/delta; range = sopp/.9 + delta; step = Min(delta + delta, sopp); - VecCopy(y, z); + XCopy(y, z); fy = fz; for( i = 0; i < nfree; ++i ) p[i] = -p[i]; @@ -535,15 +517,18 @@ static real LocalSearch(This *t, ccount nfree, ccount *ifree, static real FindMinimum(This *t, cBounds *b, real *xmin, real fmin) { - real hessian[NDIM*NDIM]; - real gfree[NDIM], p[NDIM]; - real tmp[NDIM], ftmp, fini = fmin; + Vector(real, hessian, NDIM*NDIM); + Vector(real, gfree, NDIM); + Vector(real, p, NDIM); + Vector(real, tmp, NDIM); + Vector(count, ifree, NDIM); + Vector(count, ifix, NDIM); + real ftmp, fini = fmin; ccount maxeval = t->neval + 50*t->ndim; count nfree, nfix; - count ifree[NDIM], ifix[NDIM]; count dim, local; - Zap(hessian); + Clear(hessian, t->ndim*t->ndim); for( dim = 0; dim < t->ndim; ++dim ) Hessian(dim, dim) = 1; @@ -580,7 +565,7 @@ static real FindMinimum(This *t, cBounds *b, real *xmin, real fmin) if( ftmp > fmin - (1 + fabs(fmin))*RTEPS ) goto releasebounds; fmin = ftmp; - VecCopy(xmin, tmp); + XCopy(xmin, tmp); } while( t->neval <= maxeval ) { @@ -632,13 +617,13 @@ fixbound: tmp[i] = Hessian(i + 1, mini); for( i = mini; i < nfree; ++i ) { - Copy(&Hessian(i, 0), &Hessian(i + 1, 0), i); + Move(&Hessian(i, 0), &Hessian(i + 1, 0), i); Hessian(i, i) = Hessian(i + 1, i + 1); } RenormalizeCholesky(t, nfree, hessian, tmp, diag); - Copy(&ifree[mini], &ifree[mini + 1], nfree - mini); - Copy(&gfree[mini], &gfree[mini + 1], nfree - mini); + Move(&ifree[mini], &ifree[mini + 1], nfree - mini); + Move(&gfree[mini], &gfree[mini + 1], nfree - mini); } continue; } @@ -651,11 +636,11 @@ fixbound: real fdiff; fmin = low.f; - VecCopy(xmin, tmp); + XCopy(xmin, tmp); Gradient(t, nfree, ifree, b, xmin, fmin, tmp); BFGS(t, nfree, hessian, tmp, gfree, p, low.dx); - VecCopy(gfree, tmp); + XCopy(gfree, tmp); if( fabs(low.dx - minstep) < QEPS*minstep ) goto fixbound; @@ -696,7 +681,7 @@ releasebounds: ++nfree; --nfix; - Copy(&ifix[mini], &ifix[mini + 1], nfix - mini); + Move(&ifix[mini], &ifix[mini + 1], nfix - mini); continue; } } diff --git a/src/external/libCuba/src/divonne/Integrate.c b/src/external/libCuba/src/divonne/Integrate.c index 72aa9f75..c40e9502 100644 --- a/src/external/libCuba/src/divonne/Integrate.c +++ b/src/external/libCuba/src/divonne/Integrate.c @@ -4,189 +4,226 @@ has approximately equal spread = 1/2 vol (max - min), then do a main integration over all regions this file is part of Divonne - last modified 8 Jun 10 th + checkpointing by B. Chokoufe + last modified 5 Aug 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - -#define INIDEPTH 3 -#define DEPTH 5 -#define POSTDEPTH 15 - -/*********************************************************************/ +typedef struct { + signature_t signature; + number neval, neval_opt, neval_cut; + number nmin, nrand; + count iregion, nregions; + count phase, iter, pass, size; + Totals totals[]; +} State; static int Integrate(This *t, real *integral, real *error, real *prob) { - TYPEDEFREGION; + StateDecl; + csize_t statesize = sizeof(State) + NCOMP*sizeof(Totals); + Sized(State, state, statesize); + csize_t regionsize = RegionSize; + Vector(char, out, 64*NDIM + 256*NCOMP); - Totals totals[NCOMP]; - real nneed, weight; - count dim, comp, iter, pass = 0, err, iregion; - number nwant, nmin = INT_MAX; - int fail = -99; + Totals *tot, *Tot = state->totals + t->ncomp; + Bounds *b, *B; + Result *res; + count comp, err, iregion; + number nwant; + real nneed; + ML_ONLY(number neff;) + int fail; if( VERBOSE > 1 ) { - char s[512]; - sprintf(s, "Divonne input parameters:\n" + sprintf(out, "Divonne input parameters:\n" " ndim " COUNT "\n ncomp " COUNT "\n" " epsrel " REAL "\n epsabs " REAL "\n" " flags %d\n seed %d\n" " mineval " NUMBER "\n maxeval " NUMBER "\n" " key1 %d\n key2 %d\n key3 %d\n maxpass " COUNT "\n" " border " REAL "\n maxchisq " REAL "\n mindeviation " REAL "\n" - " ngiven " NUMBER "\n nextra " NUMBER "\n", + " ngiven " NUMBER "\n nextra " NUMBER "\n" + " statefile \"%s\"", t->ndim, t->ncomp, t->epsrel, 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, - t->ngiven, t->nextra); - Print(s); + t->ngiven, t->nextra, + t->statefile); + Print(out); } if( BadComponent(t) ) return -2; if( BadDimension(t, t->key1) || - BadDimension(t, t->key2) || + BadDimension(t, t->key2) || ((t->key3 & -2) && BadDimension(t, t->key3)) ) return -1; - t->neval_opt = t->neval_cut = 0; + FORK_ONLY(t->nframe = 0;) - t->size = CHUNKSIZE; - MemAlloc(t->voidregion, t->size*sizeof(Region)); - for( dim = 0; dim < t->ndim; ++dim ) { - Bounds *b = &RegionPtr(0)->bounds[dim]; - b->lower = 0; - b->upper = 1; - } - - RuleIni(&t->rule7); - RuleIni(&t->rule9); - RuleIni(&t->rule11); - RuleIni(&t->rule13); + AllocGiven(t); SamplesIni(&t->samples[0]); SamplesIni(&t->samples[1]); SamplesIni(&t->samples[2]); - - if( setjmp(t->abort) ) goto abort; - - t->epsabs = Max(t->epsabs, NOTZERO); - - /* Step 1: partition the integration region */ - - if( VERBOSE ) Print("Partitioning phase:"); - - if( IsSobol(t->key1) || IsSobol(t->key2) || IsSobol(t->key3) ) + RuleAlloc(t); + if( IsSobol(t->key1) | IsSobol(t->key2) | IsSobol(t->key3) ) IniRandom(t); + t->epsabs = Max(t->epsabs, NOTZERO); + t->totals = state->totals; + + ForkCores(t); + + if( (fail = setjmp(t->abort)) ) goto abort; SamplesLookup(t, &t->samples[0], t->key1, (number)47, (number)INT_MAX, (number)0); SamplesAlloc(t, &t->samples[0]); - t->totals = totals; - Zap(totals); - t->phase = 1; + StateSetup(t); - Explore(t, 0, &t->samples[0], INIDEPTH, 1); + if( StateReadTest(t) ) { + StateReadOpen(t, fd) { + if( read(fd, state, statesize) != statesize || + state->signature != StateSignature(t, 3) ) break; + t->nregions = state->nregions; + if( st.st_size != statesize + t->nregions*regionsize ) break; + t->neval = state->neval; + t->neval_opt = state->neval_opt; + t->neval_cut = state->neval_cut; + t->nrand = state->nrand; + t->phase = state->phase; + t->size = state->size; + AllocRegions(t); + StateRead(fd, t->region, t->nregions*regionsize); + } StateReadClose(t, fd); - for( iter = 1; ; ++iter ) { - Totals *maxtot; - count valid; + if( IsSobol(t->key1) | IsSobol(t->key2) | IsSobol(t->key3) ) + t->rng.skiprandom(t, t->nrand); + } - for( comp = 0; comp < t->ncomp; ++comp ) { - Totals *tot = &totals[comp]; - tot->avg = tot->spreadsq = 0; - tot->spread = tot->secondspread = -INFTY; + if( ini ) { +#if MLVERSION + t->neval = t->ngiven; +#else + t->neval = 0; +#endif + t->neval_opt = 0; + t->neval_cut = 0; + t->nrand = 0; + + t->size = CHUNKSIZE; + AllocRegions(t); + for( B = (b = t->region->bounds) + t->ndim; b < B; ++b ) { + b->lower = 0; + b->upper = 1; } + t->nregions = 1; - for( iregion = 0; iregion < t->nregions; ++iregion ) { - Region *region = RegionPtr(iregion); - for( comp = 0; comp < t->ncomp; ++comp ) { - cResult *r = ®ion->result[comp]; - Totals *tot = &totals[comp]; - tot->avg += r->avg; - tot->spreadsq += Sq(r->spread); - if( r->spread > tot->spread ) { - tot->secondspread = tot->spread; - tot->spread = r->spread; - tot->iregion = iregion; + t->phase = 1; + state->iter = 1; + state->pass = 0; + state->nmin = INT_MAX; + state->iregion = 0; + FClear(state->totals); + } + + /* Step 1: partition the integration region */ + + if( t->phase == 1 ) { + if( VERBOSE ) Print("Partitioning phase:"); + + if( ini ) Iterate(t, 0, INIDEPTH, 0, NULL); + + for( ; ; ++state->iter ) { + Totals *maxtot; + count valid; + + for( tot = state->totals; tot < Tot; ++tot ) { + tot->avg = tot->spreadsq = 0; + tot->spread = tot->secondspread = -INFTY; + } + + for( iregion = 0; iregion < t->nregions; ++iregion ) + for( res = RegionResult(RegionPtr(iregion)), tot = state->totals; + tot < Tot; ++res, ++tot ) { + tot->avg += res->avg; + tot->spreadsq += Sq(res->spread); + if( res->spread > tot->spread ) { + tot->secondspread = tot->spread; + tot->spread = res->spread; + tot->iregion = iregion; + } + else if( res->spread > tot->secondspread ) + tot->secondspread = res->spread; } - else if( r->spread > tot->secondspread ) - tot->secondspread = r->spread; + + valid = 0; + for( maxtot = tot = state->totals, comp = 0; tot < Tot; ++tot, ++comp ) { + integral[comp] = tot->avg; + valid += tot->avg == tot->avg; + if( tot->spreadsq > maxtot->spreadsq ) maxtot = tot; + tot->spread = sqrt(tot->spreadsq); + error[comp] = tot->spread/t->samples[0].neff; } - } - maxtot = totals; - valid = 0; - for( comp = 0; comp < t->ncomp; ++comp ) { - Totals *tot = &totals[comp]; - integral[comp] = tot->avg; - valid += tot->avg == tot->avg; - if( tot->spreadsq > maxtot->spreadsq ) maxtot = tot; - tot->spread = sqrt(tot->spreadsq); - error[comp] = tot->spread*t->samples[0].weight; - } +#define WriteState(t) \ +if( StateWriteTest(t) ) { \ + StateWriteOpen(t, fd) { \ + state->signature = StateSignature(t, 3); \ + state->neval = t->neval; \ + state->neval_opt = t->neval_opt; \ + state->neval_cut = t->neval_cut; \ + state->nrand = t->nrand; \ + state->nregions = t->nregions; \ + state->phase = t->phase; \ + state->size = t->size; \ + StateWrite(fd, state, statesize); \ + StateWrite(fd, t->region, t->nregions*regionsize); \ + } StateWriteClose(t, fd); \ +} - if( VERBOSE ) { - char s[128 + 64*NCOMP], *p = s; + WriteState(t); - p += sprintf(p, "\n" - "Iteration " COUNT " (pass " COUNT "): " COUNT " regions\n" - NUMBER7 " integrand evaluations so far,\n" - NUMBER7 " in optimizing regions,\n" - NUMBER7 " in finding cuts", - iter, pass, t->nregions, t->neval, t->neval_opt, t->neval_cut); - - for( comp = 0; comp < t->ncomp; ++comp ) - p += sprintf(p, "\n[" COUNT "] " - REAL " +- " REAL, - comp + 1, integral[comp], error[comp]); - - Print(s); - } - - if( valid == 0 ) goto abort; /* all NaNs */ - - if( t->neval > t->maxeval ) break; - - nneed = maxtot->spread/MaxErr(maxtot->avg); - if( nneed < MAXPRIME ) { - cnumber n = t->neval + t->nregions*(number)ceil(nneed); - if( n < nmin ) { - nmin = n; - pass = 0; + if( VERBOSE ) { + char *oe = out + sprintf(out, "\n" + "Iteration " COUNT " (pass " COUNT "): " COUNT " regions\n" + NUMBER7 " integrand evaluations so far,\n" + NUMBER7 " in optimizing regions,\n" + NUMBER7 " in finding cuts", + state->iter, state->pass, t->nregions, + t->neval, t->neval_opt, t->neval_cut); + for( comp = 0; comp < t->ncomp; ++comp ) + oe += sprintf(oe, "\n[" COUNT "] " + REAL " +- " REAL, + comp + 1, integral[comp], error[comp]); + Print(out); } - else if( ++pass > t->maxpass && n >= t->mineval ) break; - } - Split(t, maxtot->iregion, DEPTH); + if( valid == 0 ) goto abort; /* all NaNs */ + + if( t->neval > t->maxeval ) break; + + nneed = maxtot->spread/MaxErr(maxtot->avg); + if( nneed < MAXPRIME ) { + cnumber n = t->neval + t->nregions*(number)ceil(nneed); + if( n < state->nmin ) { + state->nmin = n; + state->pass = 0; + } + else if( ++state->pass > t->maxpass && n >= t->mineval ) break; + } + + Iterate(t, maxtot->iregion, DEPTH, -1, NULL); + } } /* Step 2: do a "full" integration on each region */ /* nneed = t->samples[0].neff + 1; */ nneed = 2*t->samples[0].neff; - for( comp = 0; comp < t->ncomp; ++comp ) { - Totals *tot = &totals[comp]; + for( tot = state->totals; tot < Tot; ++tot ) { creal maxerr = MaxErr(tot->avg); tot->nneed = tot->spread/maxerr; nneed = Max(nneed, tot->nneed); @@ -205,85 +242,73 @@ static int Integrate(This *t, real *integral, real *error, real *prob) if( VERBOSE ) Print("\nNot enough samples left for main integration."); for( comp = 0; comp < t->ncomp; ++comp ) prob[comp] = -999; - weight = t->samples[0].weight; + ML_ONLY(neff = t->samples[0].neff;) } else { - bool can_adjust = (t->key3 == 1 && t->samples[1].sampler != SampleRule && + bool can_adjust = (t->key3 == 1 && + t->samples[1].sampler != SampleRule && (t->key2 < 0 || t->samples[1].neff < MAXPRIME)); count df, nlimit; SamplesAlloc(t, &t->samples[1]); if( VERBOSE ) { - char s[128]; - sprintf(s, "\nMain integration on " COUNT + sprintf(out, "\nMain integration on " COUNT " regions with " NUMBER " samples per region.", t->nregions, t->samples[1].neff); - Print(s); + Print(out); } - ResClear(integral); - ResClear(error); - ResClear(prob); - nlimit = t->maxeval - t->nregions*t->samples[1].n; df = 0; - for( iregion = 0; iregion < t->nregions; ++iregion ) { - Region *region = RegionPtr(iregion); - char s[64*NDIM + 256*NCOMP], *p = s; +#define CopyPhaseResults(f) \ + for( res = RegionResult(region), tot = state->totals; tot < Tot; ++res, ++tot ) { \ + PhaseResult *p = &tot->phase[f]; \ + p->avg = res->avg; \ + p->err = res->err; \ + } + +#define Var2(f, r) Sq((r)->err ? (r)->err : res->spread/t->samples[f].neff) +#define Var(f) Var2(f, &tot->phase[f]) + + while( state->iregion < t->nregions ) { + Region *region; + char *oe = out; int todo; refine: + region = RegionPtr(state->iregion); + CopyPhaseResults(0); t->phase = 2; - t->samples[1].sampler(t, &t->samples[1], region->bounds, region->vol); + region->isamples = 1; + t->samples[1].sampler(t, state->iregion); + CopyPhaseResults(1); if( can_adjust ) - for( comp = 0; comp < t->ncomp; ++comp ) - totals[comp].spreadsq -= Sq(region->result[comp].spread); + for( res = RegionResult(region), tot = state->totals; + tot < Tot; ++res, ++tot ) + tot->spreadsq -= Sq(res->spread); nlimit += t->samples[1].n; todo = 0; - for( comp = 0; comp < t->ncomp; ++comp ) { - cResult *r = ®ion->result[comp]; - Totals *tot = &totals[comp]; - - t->samples[0].avg[comp] = r->avg; - t->samples[0].err[comp] = r->err; - + for( res = RegionResult(region), tot = state->totals; + tot < Tot; ++tot ) { if( t->neval < nlimit ) { - creal avg2 = t->samples[1].avg[comp]; - creal err2 = t->samples[1].err[comp]; - creal diffsq = Sq(avg2 - r->avg); + creal avg2 = tot->phase[1].avg; + creal diffsq = Sq(avg2 - tot->phase[0].avg); -#define Var(s) Sq((s.err[comp] == 0) ? r->spread*s.weight : s.err[comp]) - - if( err2*tot->nneed > r->spread || - diffsq > Max(t->maxchisq*(Var(t->samples[0]) + Var(t->samples[1])), - EPS*Sq(avg2)) ) { + if( res->err*tot->nneed > res->spread || + diffsq > Max(t->maxchisq*(Var(0) + Var(1)), EPS*Sq(avg2)) ) { if( t->key3 && diffsq > tot->mindevsq ) { if( t->key3 == 1 ) { - ccount xregion = t->nregions; - if( VERBOSE > 2 ) Print("\nSplit"); - t->phase = 1; - Explore(t, iregion, &t->samples[1], POSTDEPTH, 2); + Iterate(t, state->iregion, POSTDEPTH, 1, state->totals); if( can_adjust ) { - number nnew; - count ireg, xreg; - - for( ireg = iregion, xreg = xregion; - ireg < t->nregions; ireg = xreg++ ) { - cResult *result = RegionPtr(ireg)->result; - count c; - for( c = 0; c < t->ncomp; ++c ) - totals[c].spreadsq += Sq(result[c].spread); - } - - nnew = (tot->spreadsq/Sq(MARKMASK) > tot->maxerrsq) ? + cnumber nnew = (tot->spreadsq/Sq(MARKMASK) > tot->maxerrsq) ? MARKMASK : (number)ceil(sqrt(tot->spreadsq/tot->maxerrsq)); if( nnew > nwant + nwant/64 ) { @@ -298,15 +323,13 @@ refine: can_adjust = false; if( VERBOSE > 2 ) { - char s[128]; - sprintf(s, "Sampling remaining " COUNT + sprintf(out, "Sampling remaining " COUNT " regions with " NUMBER " points per region.", t->nregions, t->samples[1].neff); - Print(s); + Print(out); } } } - goto refine; } todo |= 3; @@ -316,15 +339,15 @@ refine: } } - if( can_adjust ) { - for( comp = 0; comp < t->ncomp; ++comp ) - totals[comp].maxerrsq -= - Sq(region->result[comp].spread*t->samples[1].weight); - } + if( can_adjust ) + for( res = RegionResult(region), tot = state->totals; + tot < Tot; ++res, ++tot ) + tot->maxerrsq -= Sq(res->spread/t->samples[1].neff); switch( todo ) { case 1: /* get spread right */ - Explore(t, iregion, &t->samples[1], 0, 2); + region->isamples = 1; + ExploreSerial(t, state->iregion); break; case 3: /* sample region again with more points */ @@ -334,144 +357,142 @@ refine: SamplesAlloc(t, &t->samples[2]); } t->phase = 3; - t->samples[2].sampler(t, &t->samples[2], region->bounds, region->vol); - Explore(t, iregion, &t->samples[2], 0, 2); + region->isamples = 2; + t->samples[2].sampler(t, state->iregion); + ExploreSerial(t, state->iregion); ++region->depth; /* misused for df here */ ++df; } - ++region->depth; /* misused for df here */ - if( VERBOSE > 2 ) { - for( dim = 0; dim < t->ndim; ++dim ) { - cBounds *b = ®ion->bounds[dim]; - p += sprintf(p, - (dim == 0) ? "\nRegion (" REALF ") - (" REALF ")" : - "\n (" REALF ") - (" REALF ")", - b->lower, b->upper); + cchar *msg = "\nRegion (" REALF ") - (" REALF ")"; + for( B = (b = region->bounds) + t->ndim; b < B; ++b ) { + oe += sprintf(oe, msg, b->lower, b->upper); + msg = "\n (" REALF ") - (" REALF ")"; } } - for( comp = 0; comp < t->ncomp; ++comp ) { - Result *r = ®ion->result[comp]; - - creal x1 = t->samples[0].avg[comp]; - creal s1 = Var(t->samples[0]); - creal x2 = t->samples[1].avg[comp]; - creal s2 = Var(t->samples[1]); - creal r2 = (s1 == 0) ? Sq(t->samples[1].neff*t->samples[0].weight) : s2/s1; + for( tot = state->totals, res = RegionResult(region), comp = 0; + tot < Tot; ++tot, ++res ) { + creal x1 = tot->phase[0].avg; + creal v1 = Var(0); + creal x2 = tot->phase[1].avg; + creal v2 = Var(1); + creal r2 = v1 ? v2/v1 : + Sq(t->samples[1].neff/(real)t->samples[0].neff); real norm = 1 + r2; real avg = x2 + r2*x1; - real sigsq = s2; + real sigsq = v2; real chisq = Sq(x2 - x1); - real chiden = s1 + s2; + real chiden = v1 + v2; if( todo == 3 ) { - creal x3 = t->samples[2].avg[comp]; - creal s3 = Var(t->samples[2]); - creal r3 = (s2 == 0) ? Sq(t->samples[2].neff*t->samples[1].weight) : s3/s2; + creal x3 = res->avg; + creal v3 = Var2(2, res); + creal r3 = v2 ? v3/v2 : + Sq(t->samples[2].neff/(real)t->samples[1].neff); norm = 1 + r3*norm; avg = x3 + r3*avg; - sigsq = s3; - chisq = s1*Sq(x3 - x2) + s2*Sq(x3 - x1) + s3*chisq; - chiden = s1*s2 + s3*chiden; + sigsq = v3; + chisq = v1*Sq(x3 - x2) + v2*Sq(x3 - x1) + v3*chisq; + chiden = v1*v2 + v3*chiden; } - avg = LAST ? r->avg : (sigsq *= norm = 1/norm, avg*norm); + avg = LAST ? res->avg : (sigsq *= norm = 1/norm, avg*norm); if( chisq > EPS ) chisq /= Max(chiden, NOTZERO); -#define Out(s) s.avg[comp], r->spread*s.weight, s.err[comp] - if( VERBOSE > 2 ) { - p += sprintf(p, "\n[" COUNT "] " +#define Out2(f, r) (r)->avg, res->spread/t->samples[f].neff, (r)->err +#define Out(f) Out2(f, &tot->phase[f]) + oe += sprintf(oe, "\n[" COUNT "] " REAL " +- " REAL "(" REAL ")\n " - REAL " +- " REAL "(" REAL ")", - comp + 1, Out(t->samples[0]), Out(t->samples[1])); - if( todo == 3 ) p += sprintf(p, "\n " - REAL " +- " REAL "(" REAL ")", - Out(t->samples[2])); - p += sprintf(p, " \tchisq " REAL, chisq); + 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); } - integral[comp] += avg; - error[comp] += sigsq; - prob[comp] += chisq; + tot->integral += avg; + tot->sigsq += sigsq; + tot->chisq += chisq; - r->avg = avg; - r->spread = sqrt(sigsq); - r->chisq = chisq; + res->avg = avg; + res->spread = sqrt(sigsq); + res->chisq = chisq; } - if( VERBOSE > 2 ) Print(s); - } + if( VERBOSE > 2 ) Print(out); + ++state->iregion; - for( comp = 0; comp < t->ncomp; ++comp ) - error[comp] = sqrt(error[comp]); + WriteState(t); + } df += t->nregions; - if( VERBOSE > 2 ) { - char s[16 + 128*NCOMP], *p = s; - - p += sprintf(p, "\nTotals:"); - - for( comp = 0; comp < t->ncomp; ++comp ) - p += sprintf(p, "\n[" COUNT "] " - REAL " +- " REAL " \tchisq " REAL " (" COUNT " df)", - comp + 1, integral[comp], error[comp], prob[comp], df); - - Print(s); + for( tot = state->totals, comp = 0; tot < Tot; ++tot, ++comp ) { + integral[comp] = tot->integral; + error[comp] = sqrt(tot->sigsq); + prob[comp] = ChiSquare(tot->chisq, df); } - for( comp = 0; comp < t->ncomp; ++comp ) - prob[comp] = ChiSquare(prob[comp], df); + if( VERBOSE > 2 ) { + char *oe = out + sprintf(out, "\nTotals:"); + 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); + Print(out); + } - weight = 1; + ML_ONLY(neff = 1;) } #ifdef MLVERSION if( REGIONS ) { + Vector(real, bounds, t->ndim*2); + MLPutFunction(stdlink, "List", 2); + MLPutFunction(stdlink, "List", t->nregions); for( iregion = 0; iregion < t->nregions; ++iregion ) { Region *region = RegionPtr(iregion); - cBounds *b = region->bounds; - real lower[NDIM], upper[NDIM]; + cResult *Res; + real *d = bounds; - for( dim = 0; dim < t->ndim; ++dim ) { - lower[dim] = b[dim].lower; - upper[dim] = b[dim].upper; + for( B = (b = region->bounds) + t->ndim; b < B; ++b ) { + *d++ = b->lower; + *d++ = b->upper; } MLPutFunction(stdlink, "Cuba`Divonne`region", 4); - MLPutRealList(stdlink, lower, t->ndim); - MLPutRealList(stdlink, upper, t->ndim); + MLPutRealList(stdlink, bounds, 2*t->ndim); MLPutFunction(stdlink, "List", t->ncomp); - for( comp = 0; comp < t->ncomp; ++comp ) { - cResult *r = ®ion->result[comp]; - real res[] = {r->avg, r->spread*weight, r->chisq}; - MLPutRealList(stdlink, res, Elements(res)); + for( Res = (res = RegionResult(region)) + t->ncomp; res < Res; ++res ) { + real r[] = {res->avg, res->spread/neff, res->chisq}; + MLPutRealList(stdlink, r, Elements(r)); } - MLPutInteger(stdlink, region->depth); /* misused for df */ + MLPutInteger(stdlink, region->depth + 1); /* misused for df */ } } #endif abort: + WaitCores(t); + FORK_ONLY(FrameFree(t, ShmRm(t));) + + RuleFree(t); SamplesFree(&t->samples[2]); SamplesFree(&t->samples[1]); SamplesFree(&t->samples[0]); - RuleFree(&t->rule13); - RuleFree(&t->rule11); - RuleFree(&t->rule9); - RuleFree(&t->rule7); + free(t->region); + free(t->xgiven); - free(t->voidregion); + StateRemove(t); return fail; } diff --git a/src/external/libCuba/src/divonne/Iterate.c b/src/external/libCuba/src/divonne/Iterate.c new file mode 100644 index 00000000..3d9488de --- /dev/null +++ b/src/external/libCuba/src/divonne/Iterate.c @@ -0,0 +1,132 @@ +/* + Iterate.c + recursion over regions + this file is part of Divonne + last modified 2 Aug 13 th +*/ + + +static void Iterate(This *t, count iregion, cint depth, cint isamples, + Totals *totals) +{ + csize_t regionsize = RegionSize; + Region *parent, *region; + typedef struct { + real avg, err, spread, spreadsq; + } Corr; + Vector(Corr, corr, NCOMP); + Corr *c, *C = corr + t->ncomp; + Result *res; + count ireg, mreg = iregion; + count comp, maxsplit; + int last, idest, isrc; + + region = RegionPtr(iregion); + region->depth = depth; + region->next = -iregion - 1; + if( isamples < 0 ) Split(t, iregion); + else { + region->isamples = isamples; + ExploreSerial(t, iregion); + } + + ireg = iregion + RegionPtr(iregion)->next; + + do { + region = RegionPtr(ireg); + if( region->depth > 0 ) { + --region->depth; +FORK_ONLY(more:) + ireg = Explore(t, ireg); + if( ireg == -1 ) return; + region = RegionPtr(ireg); + } + if( region->depth < 0 ) mreg = IMax(mreg, ireg); + ireg += region->next; + } while( ireg > 0 ); + + FORK_ONLY(if( t->running ) goto more;) + + maxsplit = 1; + for( ireg = mreg; ireg >= iregion; --ireg ) { + parent = RegionPtr(ireg); + maxsplit -= NegQ(parent->depth); + if( parent->depth < 0 ) { + count xreg; + struct { + count from, to; + } todo[maxsplit], *tdmax = todo, *td; + count nsplit = 0; + real norm; + + FClear(corr); + + tdmax->from = ireg + parent->next; + tdmax->to = tdmax->from - parent->depth; + ++tdmax; + for( td = todo; td < tdmax; ++td ) { + for( xreg = td->from; xreg < td->to; ++xreg ) { + Region *region = RegionPtr(xreg); + if( region->depth < 0 ) { + tdmax->from = xreg + region->next; + tdmax->to = tdmax->from - region->depth; + ++tdmax; + } + else { + ++nsplit; + for( res = RegionResult(region), c = corr; c < C; ++res, ++c ) { + c->avg += res->avg; + c->err += res->err; + c->spread += Sq(res->spread); + } + } + } + } + + norm = 1./nsplit--; + for( res = RegionResult(parent), c = corr; c < C; ++res, ++c ) { + creal diff = fabs(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); + } + + for( td = todo; td < tdmax; ++td ) + for( xreg = td->from; xreg < td->to; ++xreg ) { + Region *region = RegionPtr(xreg); + if( region->depth >= 0 ) { + cnumber neff = t->samples[region->isamples].neff; + for( res = RegionResult(region), c = corr; c < C; ++res, ++c ) { + if( res->err > 0 ) res->err = res->err*c->err + c->avg; + res->spread = res->spread*c->spread + c->avg*neff; + c->spreadsq += Sq(res->spread); + } + } + } + } + } + + if( totals ) + for( comp = 0; comp < t->ncomp; ++comp ) + totals[comp].spreadsq += corr[comp].spreadsq; + + for( last = -1, idest = isrc = iregion; iregion <= mreg; ++iregion ) { + Region *region = RegionPtr(iregion); + cint cur = NegQ(region->depth); + switch( cur - last ) { + case -1: + memmove(RegionPtr(idest), RegionPtr(isrc), + (iregion - isrc)*regionsize); + idest += iregion - isrc; + break; + case 1: + isrc = iregion; + } + last = cur; + } + + memmove(RegionPtr(idest), RegionPtr(iregion), + (t->nregions - iregion)*regionsize); + t->nregions += idest - iregion; +} + diff --git a/src/external/libCuba/src/divonne/KorobovCoeff.c b/src/external/libCuba/src/divonne/KorobovCoeff.c index ae172002..e8ace85e 100644 --- a/src/external/libCuba/src/divonne/KorobovCoeff.c +++ b/src/external/libCuba/src/divonne/KorobovCoeff.c @@ -1,23 +1,3 @@ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - #define KOROBOV_MINDIM 2 #define KOROBOV_MAXDIM 33 #define MAXPRIME 9689 diff --git a/src/external/libCuba/src/divonne/Makefile.am b/src/external/libCuba/src/divonne/Makefile.am new file mode 100644 index 00000000..51f30f38 --- /dev/null +++ b/src/external/libCuba/src/divonne/Makefile.am @@ -0,0 +1,13 @@ +## Process this file with automake to create Makefile.in +## $Id: Makefile.am 4808 2011-03-21 14:33:34Z l_wojek $ + +c_sources = Divonne.c + +AM_CPPFLAGS = -I. -I.. -I../common -DNOUNDERSCORE +AM_CFLAGS = $(LOCAL_CUBA_LIB_CFLAGS) +AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS) + +noinst_LTLIBRARIES = libdivonne.la + +libdivonne_la_SOURCES = $(c_sources) +libdivonne_la_LDFLAGS = $(AM_LDFLAGS) \ No newline at end of file diff --git a/src/external/libCuba/src/divonne/Rule.c b/src/external/libCuba/src/divonne/Rule.c index 770bd06e..5bb6c28b 100644 --- a/src/external/libCuba/src/divonne/Rule.c +++ b/src/external/libCuba/src/divonne/Rule.c @@ -4,37 +4,11 @@ code lifted with minor modifications from DCUHRE by J. Berntsen, T. Espelid, and A. Genz this file is part of Divonne - last modified 16 Jun 10 th + last modified 26 Jul 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ -enum { nrules = 5 }; - -#define TYPEDEFSET \ - typedef struct { \ - count n; \ - real weight[nrules], scale[nrules], norm[nrules]; \ - real gen[NDIM]; \ - } Set +#define NextSet(p) p = (Set *)((char *)p + setsize) /*********************************************************************/ @@ -83,86 +57,84 @@ static void Rule13Alloc(This *t) enum { nsets = 14, ndim = 2 }; - TYPEDEFSET; - count n, r; Set *first, *last, *s, *x; + csize_t setsize = SetSize; - Alloc(first, nsets); - Clear(first, nsets); + Die(first = calloc(nsets, setsize)); last = first; n = last->n = 1; Copy(last->weight, w[0], nrules); - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[1], nrules); last->gen[0] = g[0]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[2], nrules); last->gen[0] = g[1]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[3], nrules); last->gen[0] = g[2]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[4], nrules); last->gen[0] = g[3]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[5], nrules); last->gen[0] = g[4]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); Copy(last->weight, w[6], nrules); last->gen[0] = g[5]; last->gen[1] = g[5]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); Copy(last->weight, w[7], nrules); last->gen[0] = g[6]; last->gen[1] = g[6]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); Copy(last->weight, w[8], nrules); last->gen[0] = g[7]; last->gen[1] = g[7]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); Copy(last->weight, w[9], nrules); last->gen[0] = g[8]; last->gen[1] = g[8]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); Copy(last->weight, w[10], nrules); last->gen[0] = g[9]; last->gen[1] = g[9]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1); Copy(last->weight, w[11], nrules); last->gen[0] = g[10]; last->gen[1] = g[11]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1); Copy(last->weight, w[12], nrules); last->gen[0] = g[12]; last->gen[1] = g[13]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1); Copy(last->weight, w[13], nrules); last->gen[0] = g[14]; @@ -175,12 +147,12 @@ static void Rule13Alloc(This *t) t->rule13.errcoeff[2] = 5; t->rule13.n = n; - for( s = first; s <= last; ++s ) + for( s = first; s <= last; NextSet(s) ) for( r = 1; r < nrules - 1; ++r ) { creal scale = (s->weight[r] == 0) ? 100 : -s->weight[r + 1]/s->weight[r]; real sum = 0; - for( x = first; x <= last; ++x ) + for( x = first; x <= last; NextSet(x) ) sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]); s->scale[r] = scale; s->norm[r] = 1/sum; @@ -231,84 +203,82 @@ static void Rule11Alloc(This *t) enum { nsets = 13, ndim = 3 }; - TYPEDEFSET; - count n, r; Set *first, *last, *s, *x; + csize_t setsize = SetSize; - Alloc(first, nsets); - Clear(first, nsets); + Die(first = calloc(nsets, setsize)); last = first; n = last->n = 1; Copy(last->weight, w[0], nrules); - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[1], nrules); last->gen[0] = g[0]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[2], nrules); last->gen[0] = g[1]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[3], nrules); last->gen[0] = g[2]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[4], nrules); last->gen[0] = g[3]; - ++last; + NextSet(last); n += last->n = 2*ndim; Copy(last->weight, w[5], nrules); last->gen[0] = g[4]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); Copy(last->weight, w[6], nrules); last->gen[0] = g[5]; last->gen[1] = g[5]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); Copy(last->weight, w[7], nrules); last->gen[0] = g[6]; last->gen[1] = g[6]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3; Copy(last->weight, w[8], nrules); last->gen[0] = g[7]; last->gen[1] = g[7]; last->gen[2] = g[7]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3; Copy(last->weight, w[9], nrules); last->gen[0] = g[8]; last->gen[1] = g[8]; last->gen[2] = g[8]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3; Copy(last->weight, w[10], nrules); last->gen[0] = g[9]; last->gen[1] = g[9]; last->gen[2] = g[9]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1)*(ndim - 2); Copy(last->weight, w[11], nrules); last->gen[0] = g[10]; last->gen[1] = g[11]; last->gen[2] = g[11]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1)*(ndim - 2); Copy(last->weight, w[12], nrules); last->gen[0] = g[12]; @@ -322,12 +292,12 @@ static void Rule11Alloc(This *t) t->rule11.errcoeff[2] = 3; t->rule11.n = n; - for( s = first; s <= last; ++s ) + for( s = first; s <= last; NextSet(s) ) for( r = 1; r < nrules - 1; ++r ) { creal scale = (s->weight[r] == 0) ? 100 : -s->weight[r + 1]/s->weight[r]; real sum = 0; - for( x = first; x <= last; ++x ) + for( x = first; x <= last; NextSet(x) ) sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]); s->scale[r] = scale; s->norm[r] = 1/sum; @@ -368,15 +338,13 @@ static void Rule9Alloc(This *t) enum { nsets = 9 }; - TYPEDEFSET; - ccount ndim = t->ndim; ccount twondim = 1 << ndim; count dim, n, r; Set *first, *last, *s, *x; + csize_t setsize = SetSize; - Alloc(first, nsets); - Clear(first, nsets); + Die(first = calloc(nsets, setsize)); last = first; n = last->n = 1; @@ -386,7 +354,7 @@ static void Rule9Alloc(This *t) last->weight[3] = ndim*(ndim*w[9] + w[10]) - 1 + last->weight[0]; last->weight[4] = ndim*w[11] + 1 - last->weight[0]; - ++last; + NextSet(last); n += last->n = 2*ndim; last->weight[0] = ndim*(ndim*w[12] + w[13]) + w[14]; last->weight[1] = ndim*(ndim*w[15] + w[16]) + w[17]; @@ -395,7 +363,7 @@ static void Rule9Alloc(This *t) last->weight[4] = w[21] - last->weight[0]; last->gen[0] = g[0]; - ++last; + NextSet(last); n += last->n = 2*ndim; last->weight[0] = ndim*w[22] + w[23]; last->weight[1] = ndim*w[24] + w[25]; @@ -404,7 +372,7 @@ static void Rule9Alloc(This *t) last->weight[4] = -last->weight[0]; last->gen[0] = g[1]; - ++last; + NextSet(last); n += last->n = 2*ndim; last->weight[0] = w[29]; last->weight[1] = w[30]; @@ -413,12 +381,12 @@ static void Rule9Alloc(This *t) last->weight[4] = -w[29]; last->gen[0] = g[2]; - ++last; + NextSet(last); n += last->n = 2*ndim; last->weight[2] = w[32]; last->gen[0] = g[3]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); last->weight[0] = w[33] - ndim*w[12]; last->weight[1] = w[34] - ndim*w[15]; @@ -428,7 +396,7 @@ static void Rule9Alloc(This *t) last->gen[0] = g[0]; last->gen[1] = g[0]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1); last->weight[0] = w[36]; last->weight[1] = w[37]; @@ -438,7 +406,7 @@ static void Rule9Alloc(This *t) last->gen[0] = g[0]; last->gen[1] = g[1]; - ++last; + NextSet(last); n += last->n = 4*ndim*(ndim - 1)*(ndim - 2)/3; last->weight[0] = w[39]; last->weight[1] = w[40]; @@ -449,7 +417,7 @@ static void Rule9Alloc(This *t) last->gen[1] = g[0]; last->gen[2] = g[0]; - ++last; + NextSet(last); n += last->n = twondim; last->weight[0] = w[41]/twondim; last->weight[1] = w[7]/twondim; @@ -466,12 +434,12 @@ static void Rule9Alloc(This *t) t->rule9.errcoeff[2] = 5; t->rule9.n = n; - for( s = first; s <= last; ++s ) + for( s = first; s <= last; NextSet(s) ) for( r = 1; r < nrules - 1; ++r ) { creal scale = (s->weight[r] == 0) ? 100 : -s->weight[r + 1]/s->weight[r]; real sum = 0; - for( x = first; x <= last; ++x ) + for( x = first; x <= last; NextSet(x) ) sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]); s->scale[r] = scale; s->norm[r] = 1/sum; @@ -501,15 +469,13 @@ static void Rule7Alloc(This *t) enum { nsets = 6 }; - TYPEDEFSET; - ccount ndim = t->ndim; ccount twondim = 1 << ndim; count dim, n, r; Set *first, *last, *s, *x; + csize_t setsize = SetSize; - Alloc(first, nsets); - Clear(first, nsets); + Die(first = calloc(nsets, setsize)); last = first; n = last->n = 1; @@ -519,7 +485,7 @@ static void Rule7Alloc(This *t) last->weight[3] = ndim*(ndim*w[7] + w[8]) - w[9]; last->weight[4] = 1 - last->weight[0]; - ++last; + NextSet(last); n += last->n = 2*ndim; last->weight[0] = w[10]; last->weight[1] = w[11]; @@ -528,7 +494,7 @@ static void Rule7Alloc(This *t) last->weight[4] = -w[10]; last->gen[0] = g[1]; - ++last; + NextSet(last); n += last->n = 2*ndim; last->weight[0] = w[13] - ndim*w[0]; last->weight[1] = w[14] - ndim*w[3]; @@ -537,12 +503,12 @@ static void Rule7Alloc(This *t) last->weight[4] = -last->weight[0]; last->gen[0] = g[0]; - ++last; + NextSet(last); n += last->n = 2*ndim; last->weight[2] = w[17]; last->gen[0] = g[2]; - ++last; + NextSet(last); n += last->n = 2*ndim*(ndim - 1); last->weight[0] = -w[7]; last->weight[1] = w[18]; @@ -552,7 +518,7 @@ static void Rule7Alloc(This *t) last->gen[0] = g[0]; last->gen[1] = g[0]; - ++last; + NextSet(last); n += last->n = twondim; last->weight[0] = w[20]/twondim; last->weight[1] = w[5]/twondim; @@ -569,12 +535,12 @@ static void Rule7Alloc(This *t) t->rule7.errcoeff[2] = 5; t->rule7.n = n; - for( s = first; s <= last; ++s ) + for( s = first; s <= last; NextSet(s) ) for( r = 1; r < nrules - 1; ++r ) { creal scale = (s->weight[r] == 0) ? 100 : -s->weight[r + 1]/s->weight[r]; real sum = 0; - for( x = first; x <= last; ++x ) + for( x = first; x <= last; NextSet(x) ) sum += x->n*fabs(x->weight[r + 1] + scale*x->weight[r]); s->scale[r] = scale; s->norm[r] = 1/sum; @@ -583,23 +549,33 @@ static void Rule7Alloc(This *t) /*********************************************************************/ -static inline void RuleIni(Rule *rule) +static inline void RuleAlloc(This *t) { - rule->first = NULL; + if( (t->ndim - 2) | (t->key1 - 13)*(t->key2 - 13)*(t->key3 - 13) ) + t->rule13.first = NULL; + else Rule13Alloc(t); + + if( (t->ndim - 3) | (t->key1 - 11)*(t->key2 - 11)*(t->key3 - 11) ) + t->rule11.first = NULL; + else Rule11Alloc(t); + + if( (t->key1 - 9)*(t->key2 - 9)*(t->key3 - 9) ) + t->rule9.first = NULL; + else Rule9Alloc(t); + + if( (t->key1 - 7)*(t->key2 - 7)*(t->key3 - 7) ) + t->rule7.first = NULL; + else Rule7Alloc(t); } /*********************************************************************/ -static inline bool RuleIniQ(Rule *rule) +static inline void RuleFree(cThis *t) { - return rule->first == NULL; -} - -/*********************************************************************/ - -static inline void RuleFree(Rule *rule) -{ - free(rule->first); + free(t->rule7.first); + free(t->rule9.first); + free(t->rule11.first); + free(t->rule13.first); } /*********************************************************************/ @@ -653,29 +629,28 @@ next: /*********************************************************************/ -static void SampleRule(This *t, cSamples *samples, cBounds *b, creal vol) +static void SampleRule(This *t, ccount iregion) { - TYPEDEFSET; - - real *x = samples->x, *f = samples->f; - Set *first = (Set *)samples->rule->first; - Set *last = (Set *)samples->rule->last; + SAMPLERDEFS; + Set *first = samples->rule->first; + Set *last = samples->rule->last; Set *s; creal *errcoeff = samples->rule->errcoeff; - count comp, rul, n; + count comp, rul, sn; + csize_t setsize = SetSize; - for( s = first; s <= last; ++s ) + for( s = first; s <= last; NextSet(s) ) if( s->n ) x = ExpandFS(t, b, s->gen, x); - DoSample(t, samples->n, t->ndim, samples->x, f); + DoSample(t, n, samples->x, f); for( comp = 0; comp < t->ncomp; ++comp ) { real sum[nrules]; creal *f1 = f++; Zap(sum); - for( s = first; s <= last; ++s ) - for( n = s->n; n; --n ) { + for( s = first; s <= last; NextSet(s) ) + for( sn = s->n; sn; --sn ) { creal fun = *f1; f1 += t->ncomp; for( rul = 0; rul < nrules; ++rul ) @@ -689,17 +664,17 @@ static void SampleRule(This *t, cSamples *samples, cBounds *b, creal vol) for( rul = 1; rul < nrules - 1; ++rul ) { real maxerr = 0; - for( s = first; s <= last; ++s ) + for( s = first; s <= last; NextSet(s) ) maxerr = Max(maxerr, fabs(sum[rul + 1] + s->scale[rul]*sum[rul])*s->norm[rul]); sum[rul] = maxerr; } - samples->avg[comp] = vol*sum[0]; - samples->err[comp] = vol*( + res[comp].avg = region->vol*sum[0]; + res[comp].err = region->vol*( (errcoeff[0]*sum[1] <= sum[2] && errcoeff[0]*sum[2] <= sum[3]) ? errcoeff[1]*sum[1] : - errcoeff[2]*Max(Max(sum[1], sum[2]), sum[3])); + errcoeff[2]*Max(Max(sum[1], sum[2]), sum[3]) ); } } diff --git a/src/external/libCuba/src/divonne/Sample.c b/src/external/libCuba/src/divonne/Sample.c index 86295e9b..7778f581 100644 --- a/src/external/libCuba/src/divonne/Sample.c +++ b/src/external/libCuba/src/divonne/Sample.c @@ -2,29 +2,9 @@ Sample.c most of what is related to sampling this file is part of Divonne - last modified 16 Jun 10 th + last modified 30 Aug 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - #define MARKMASK 0xfffffff #define Marked(x) ((x) & ~MARKMASK) @@ -38,6 +18,7 @@ static inline void SamplesIni(Samples *samples) { samples->x = NULL; + samples->n = 0; } /*********************************************************************/ @@ -56,11 +37,11 @@ static inline void SamplesFree(cSamples *samples) /*********************************************************************/ -static void SampleSobol(This *t, cSamples *samples, cBounds *b, creal vol) +static void SampleSobol(This *t, ccount iregion) { - creal norm = vol*samples->weight; - real *x = samples->x, *f = samples->f, *avg = samples->avg; - cnumber n = samples->n; + SAMPLERDEFS; + Vector(real, avg, NCOMP); + real norm; number i; count dim, comp; @@ -70,35 +51,40 @@ static void SampleSobol(This *t, cSamples *samples, cBounds *b, creal vol) *x = b[dim].lower + *x*(b[dim].upper - b[dim].lower); } - DoSample(t, n, t->ndim, samples->x, f); + t->nrand += n; - ResCopy(avg, f); + DoSample(t, n, samples->x, f); + + FCopy(avg, f); f += t->ncomp; - for( i = 1; i < n; ++i ) + for( i = 2; i < n; ++i ) for( comp = 0; comp < t->ncomp; ++comp ) avg[comp] += *f++; - for( comp = 0; comp < t->ncomp; ++comp ) - avg[comp] *= norm; + norm = region->vol/samples->neff; + for( comp = 0; comp < t->ncomp; ++comp ) { + res[comp].avg = norm*avg[comp]; + res[comp].err = 0; + } } /*********************************************************************/ -static void SampleKorobov(This *t, cSamples *samples, cBounds *b, creal vol) +static void SampleKorobov(This *t, ccount iregion) { - creal norm = vol*samples->weight; - real *x = samples->x, *xlast = x + t->ndim; - real *f = samples->f, *flast = f + t->ncomp; - real *avg = samples->avg; - cnumber n = samples->n, neff = samples->neff; - number nextra = n, i; + SAMPLERDEFS; + real *xlast = x + t->ndim, *flast = f + t->ncomp; + Vector(real, avg, NCOMP); + real norm; + cnumber neff = samples->neff; + number nextra = 0, i; real dist = 0; count dim, comp; for( i = 1; i < n; ++i ) { number c = i; for( dim = 0; dim < t->ndim; ++dim ) { - creal dx = abs(2*c - neff)*samples->weight; + creal dx = abs(2*c - neff)/(real)neff; *xlast++ = b[dim].lower + dx*(b[dim].upper - b[dim].lower); c = c*samples->coeff % neff; } @@ -119,26 +105,29 @@ static void SampleKorobov(This *t, cSamples *samples, cBounds *b, creal vol) } xlast[dim] = x2; } - ++nextra; + nextra = 1; } - DoSample(t, nextra, t->ndim, x, f); + DoSample(t, n + nextra, x, f); - ResCopy(avg, flast); + FCopy(avg, flast); flast += t->ncomp; for( i = 2; i < n; ++i ) for( comp = 0; comp < t->ncomp; ++comp ) avg[comp] += *flast++; - if( nextra > n ) { + if( nextra ) { for( comp = 0; comp < t->ncomp; ++comp ) f[comp] += dist*(f[comp] - flast[comp]); for( dim = 0; dim < t->ndim; ++dim ) x[dim] = b[dim].upper; } - for( comp = 0; comp < t->ncomp; ++comp ) - avg[comp] = (avg[comp] + avg[comp] + f[comp])*norm; + norm = region->vol/samples->neff; + for( comp = 0; comp < t->ncomp; ++comp ) { + res[comp].avg = norm*(avg[comp] + avg[comp] + f[comp]); + res[comp].err = 0; + } } /*********************************************************************/ @@ -165,25 +154,21 @@ static count SamplesLookup(This *t, Samples *samples, cint key, number n; if( key == 13 && t->ndim == 2 ) { - if( RuleIniQ(&t->rule13) ) Rule13Alloc(t); samples->rule = &t->rule13; samples->n = n = nmin = t->rule13.n; samples->sampler = SampleRule; } else if( key == 11 && t->ndim == 3 ) { - if( RuleIniQ(&t->rule11) ) Rule11Alloc(t); samples->rule = &t->rule11; samples->n = n = nmin = t->rule11.n; samples->sampler = SampleRule; } else if( key == 9 ) { - if( RuleIniQ(&t->rule9) ) Rule9Alloc(t); samples->rule = &t->rule9; samples->n = n = nmin = t->rule9.n; samples->sampler = SampleRule; } else if( key == 7 ) { - if( RuleIniQ(&t->rule7) ) Rule7Alloc(t); samples->rule = &t->rule7; samples->n = n = nmin = t->rule7.n; samples->sampler = SampleRule; @@ -234,21 +219,18 @@ static void SamplesAlloc(cThis *t, Samples *samples) Alloc(samples->x, nx + nf + t->ncomp + t->ncomp); samples->f = samples->x + nx; - samples->avg = samples->f + nf; - samples->err = samples->avg + t->ncomp; - ResClear(samples->err); - - samples->weight = 1./samples->neff; } /*********************************************************************/ static real Sample(This *t, creal *x0) { - real xtmp[2*NDIM], ftmp[2*NCOMP], *xlast = xtmp, f; + Vector(real, xtmp, 2*NDIM); + Vector(real, ftmp, 2*NCOMP); + real *xlast = xtmp, f; real dist = 0; count dim, comp; - number nextra = 1; + number n = 1; for( dim = 0; dim < t->ndim; ++dim ) { creal x1 = *xlast++ = Min(Max(*x0++, 0.), 1.); @@ -268,14 +250,16 @@ static real Sample(This *t, creal *x0) } *xlast++ = x2; } - nextra = 2; + n = 2; } - DoSample(t, nextra, t->ndim, xtmp, ftmp); + DoSample(t, n, xtmp, ftmp); + +#define fin(x) Min(Max(x, -.5*INFTY), .5*INFTY) comp = Untag(t->selectedcomp); - f = ftmp[comp]; - if( nextra > 1 ) f += dist*(f - ftmp[comp + t->ncomp]); + f = fin(ftmp[comp]); + if( n > 1 ) f += dist*(f - fin(ftmp[comp + t->ncomp])); return Sign(t->selectedcomp)*f; } diff --git a/src/external/libCuba/src/divonne/Split.c b/src/external/libCuba/src/divonne/Split.c index c7c86b79..570258f8 100644 --- a/src/external/libCuba/src/divonne/Split.c +++ b/src/external/libCuba/src/divonne/Split.c @@ -2,33 +2,13 @@ Split.c determine optimal cuts for splitting a region this file is part of Divonne - last modified 20 Jul 10 th + last modified 31 Aug 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - #define BNDTOL .05 #define FRACT .5 -#define BIG 1e10 +#define SMALL 1e-10 #define SINGTOL 1e-4 #define LHSTOL .1 @@ -47,26 +27,21 @@ typedef struct { real lhs, row, sol; } Cut; -typedef struct { - real diff, err, spread; -} Errors; - -typedef const Errors cErrors; - /*********************************************************************/ static inline real Div(creal a, creal b) { - return (b != 0 && fabs(b) < BIG*fabs(a)) ? a/b : a; + return (b != 0 /*&& fabs(a) > SMALL*fabs(b)*/) ? a/b : a; } /*********************************************************************/ static void SomeCut(This *t, Cut *cut, Bounds *b) { + Vector(real, xmid, NDIM); + real ymid, maxdev; count dim, maxdim; static count nextdim = 0; - real xmid[NDIM], ymid, maxdev; for( dim = 0; dim < t->ndim; ++dim ) xmid[dim] = .5*(b[dim].upper + b[dim].lower); @@ -111,10 +86,10 @@ static inline real Volume(cThis *t, creal *delta) /*********************************************************************/ -static inline real SetupEqs(Cut *cut, ccount ncut, real f) +static inline real SetupEqs(Cut *cut, ccount ncuts, real f) { real sqsum = 0; - Cut *c = &cut[ncut]; + Cut *c = &cut[ncuts]; while( --c >= cut ) { sqsum += Sq(c->lhs = f - c->f); f = c->f; @@ -124,7 +99,7 @@ static inline real SetupEqs(Cut *cut, ccount ncut, real f) /*********************************************************************/ -static inline void SolveEqs(Cut *cut, count ncut, +static inline void SolveEqs(Cut *cut, count ncuts, creal *delta, creal diff) { real last = 0; @@ -135,7 +110,7 @@ static inline void SolveEqs(Cut *cut, count ncut, ccount dim = Dim(c->i); c->row = r -= Div(diff, (delta[Lower(dim)] + delta[Upper(dim)])*c->df); - if( --ncut == 0 ) break; + if( --ncuts == 0 ) break; last += r*c->lhs; } @@ -157,8 +132,8 @@ static count FindCuts(This *t, Cut *cut, Bounds *bounds, creal vol, { cint sign = (fdiff < 0) ? -1 : 1; - count ncut = 0, icut; - real delta[2*NDIM]; + count ncuts = 0, icut; + Vector(real, delta, 2*NDIM); real gamma, fgamma, lhssq; count dim, div; @@ -167,7 +142,7 @@ static count FindCuts(This *t, Cut *cut, Bounds *bounds, creal vol, creal xsave = xmajor[dim]; real dist = b->upper - xsave; if( dist >= BNDTOL*(b->upper - b->lower) ) { - Cut *c = &cut[ncut++]; + Cut *c = &cut[ncuts++]; c->i = Upper(dim); c->save = dist; xmajor[dim] += dist *= FRACT; @@ -182,7 +157,7 @@ static count FindCuts(This *t, Cut *cut, Bounds *bounds, creal vol, creal xsave = xmajor[dim]; real dist = xsave - b->lower; if( dist >= BNDTOL*(b->upper - b->lower) ) { - Cut *c = &cut[ncut++]; + Cut *c = &cut[ncuts++]; c->i = Lower(dim); c->save = dist; xmajor[dim] -= dist *= FRACT; @@ -192,7 +167,7 @@ static count FindCuts(This *t, Cut *cut, Bounds *bounds, creal vol, delta[Lower(dim)] = dist; } - if( ncut == 0 ) { + if( ncuts == 0 ) { SomeCut(t, cut, bounds); return 1; } @@ -201,7 +176,7 @@ static count FindCuts(This *t, Cut *cut, Bounds *bounds, creal vol, real mindiff = INFTY; Cut *mincut = cut; - for( icut = 0; icut < ncut; ++icut ) { + for( icut = 0; icut < ncuts; ++icut ) { Cut *c = &cut[icut]; creal diff = fabs(fmajor - c->f); if( diff <= mindiff ) { @@ -215,30 +190,30 @@ static count FindCuts(This *t, Cut *cut, Bounds *bounds, creal vol, if( sign*(mincut->f - fgamma) < 0 ) break; - if( --ncut == 0 ) { + if( --ncuts == 0 ) { SomeCut(t, cut, bounds); return 1; } delta[mincut->i] = mincut->save; - memcpy(mincut, mincut + 1, (char *)&cut[ncut] - (char *)mincut); + memmove(mincut, mincut + 1, (char *)&cut[ncuts] - (char *)mincut); } - for( icut = 0; icut < ncut; ++icut ) { + for( icut = 0; icut < ncuts; ++icut ) { Cut *c = &cut[icut]; c->fold = c->f; c->df = (c->f - fmajor)/delta[c->i]; } - lhssq = SetupEqs(cut, ncut, fgamma); + lhssq = SetupEqs(cut, ncuts, fgamma); repeat: - SolveEqs(cut, ncut, delta, gamma*fdiff); + SolveEqs(cut, ncuts, delta, gamma*fdiff); for( div = 1; div <= 16; div *= 4 ) { real gammanew, lhssqnew; - for( icut = 0; icut < ncut; ++icut ) { + for( icut = 0; icut < ncuts; ++icut ) { Cut *c = &cut[icut]; real *x = &xmajor[Dim(c->i)]; creal xsave = *x; @@ -250,7 +225,7 @@ repeat: gammanew = Volume(t, delta)/vol; fgamma = fmajor + (gammanew - 1)*fdiff; - lhssqnew = SetupEqs(cut, ncut, fgamma); + lhssqnew = SetupEqs(cut, ncuts, fgamma); if( lhssqnew <= lhssq ) { real fmax; @@ -259,12 +234,12 @@ repeat: gamma = gammanew; fmax = fabs(fgamma); - for( icut = 0; icut < ncut; ++icut ) { + 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(sol) < BIG*fabs(df)) ? df/sol : 1; + df = (fabs(df) > SMALL*fabs(sol)) ? df/sol : 1; c->df = (fabs(df) < fabs(dfmin)) ? dfmin : df; fmax = Max(fmax, fabs(c->f)); c->fold = c->f; @@ -276,90 +251,62 @@ repeat: } } - for( icut = 0; icut < ncut; ++icut ) { + for( icut = 0; icut < ncuts; ++icut ) { Cut *c = &cut[icut]; real *b = (real *)bounds + c->i; c->save = *b; *b = xmajor[Dim(c->i)] + SignedDelta(c->i); } - return ncut; + return ncuts; } /*********************************************************************/ -static void Split(This *t, count iregion, int depth) +static void Split(This *t, ccount iregion) { - TYPEDEFREGION; + csize_t regionsize = RegionSize; + Region *region = RegionPtr(iregion); + Vector(Cut, cut, 2*NDIM); + Cut *c; + count ncuts, succ; + int depth; + real *b; - Cut cut[2*NDIM]; - Errors errors[NCOMP]; - count comp, ncut, nsplit, xregion, ireg, xreg; - real tmp; - -{ - Region *const region = RegionPtr(iregion); t->selectedcomp = region->cutcomp; t->neval_cut -= t->neval; - ncut = FindCuts(t, cut, region->bounds, region->vol, - (real *)region->result + region->xmajor, region->fmajor, + ncuts = FindCuts(t, cut, region->bounds, region->vol, + (real *)RegionResult(region) + region->xmajor, region->fmajor, region->fmajor - region->fminor); t->neval_cut += t->neval; - for( comp = 0; comp < t->ncomp; ++comp ) { - Errors *e = &errors[comp]; - e->diff = region->result[comp].avg; - e->spread = e->err = 0; - } -} - - xregion = t->nregions; - - depth -= ncut; - if( Explore(t, iregion, &t->samples[0], depth, 1) ) { - Cut *c; - for( c = cut; ncut--; ++c ) { - real *b = (real *)RegionPtr(iregion)->bounds; - ccount c0 = c->i, c1 = c0 ^ 1; - creal tmp = b[c1]; - b[c1] = b[c0]; - b[c0] = c->save; - if( !Explore(t, iregion, &t->samples[0], depth++, ncut != 0) ) break; - if( ncut ) ((real *)RegionPtr(iregion)->bounds)[c1] = tmp; - } - } - - nsplit = t->nregions - xregion + 1; - - for( ireg = iregion, xreg = xregion; ireg < t->nregions; ireg = xreg++ ) { - cResult *result = RegionPtr(ireg)->result; - for( comp = 0; comp < t->ncomp; ++comp ) { - cResult *r = &result[comp]; - Errors *e = &errors[comp]; - e->diff -= r->avg; - e->err += r->err; - e->spread += Sq(r->spread); - } - } - - tmp = 1./nsplit; - for( comp = 0; comp < t->ncomp; ++comp ) { - Errors *e = &errors[comp]; - e->diff = tmp*fabs(e->diff); - e->err = (e->err == 0) ? 1 : 1 + e->diff/e->err; - e->spread = (e->spread == 0) ? 1 : 1 + e->diff/sqrt(e->spread); - } - - tmp = 1 - tmp; - for( ireg = iregion, xreg = xregion; ireg < t->nregions; ireg = xreg++ ) { - Result *result = RegionPtr(ireg)->result; - for( comp = 0; comp < t->ncomp; ++comp ) { - Result *r = &result[comp]; - cErrors *e = &errors[comp]; - creal c = tmp*e->diff; - if( r->err > 0 ) r->err = r->err*e->err + c; - r->spread = r->spread*e->spread + c*t->samples[0].neff; - } + depth = region->depth - ncuts; + + EnlargeRegions(t, ++ncuts); + region = RegionPtr(iregion); + region->depth = -ncuts; + succ = iregion + region->next; + region->next = t->nregions - iregion; + b = (real *)region->bounds; + + region = RegionPtr(t->nregions); + XCopy(region->bounds, b); + region->depth = IDim(depth) + 1; + region->next = 1; + region->isamples = 0; + for( c = cut; --ncuts; ++c ) { + ccount ci = c->i; + creal tmp = b[ci ^ 1]; + b[ci ^ 1] = b[ci]; + b[ci] = c->save; + region = RegionPtr(++t->nregions); + XCopy(region->bounds, b); + region->depth = IDim(depth) + 1; + region->next = 1; + region->isamples = 0; + ++depth; + b[ci ^ 1] = tmp; } + region->next = succ - t->nregions++; } diff --git a/src/external/libCuba/src/divonne/common.c b/src/external/libCuba/src/divonne/common.c index 8e5b2b98..34be0c2b 100644 --- a/src/external/libCuba/src/divonne/common.c +++ b/src/external/libCuba/src/divonne/common.c @@ -2,47 +2,22 @@ common.c includes most of the modules this file is part of Divonne - last modified 8 Jun 10 th + last modified 29 Jul 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - #include "Random.c" #include "ChiSquare.c" #include "Rule.c" #include "Sample.c" #include "FindMinimum.c" - -static bool Explore(This *t, count iregion, cSamples *samples, - cint depth, cint flags); - -static void Split(This *t, count iregion, int depth); - -#include "Explore.c" #include "Split.c" +#include "Explore.c" +#include "Iterate.c" static inline bool BadDimension(cThis *t, ccount key) { - if( t->ndim > NDIM ) return true; + if( t->ndim > MAXDIM ) return true; if( IsSobol(key) ) return t->ndim < SOBOL_MINDIM || (t->seed == 0 && t->ndim > SOBOL_MAXDIM); if( IsRule(key, t->ndim) ) return t->ndim < 1; @@ -51,9 +26,48 @@ static inline bool BadDimension(cThis *t, ccount key) static inline bool BadComponent(cThis *t) { - if( t->ncomp > NCOMP ) return true; + if( t->ncomp > MAXCOMP ) return true; return t->ncomp < 1; } -#include "Integrate.c" +static inline void AllocGiven(This *t) +{ + real *xgiven = NULL, *fgiven = NULL; + + if( t->ngiven | t->nextra ) { + cnumber nxgiven = t->ngiven*t->ndim; + cnumber nxextra = t->nextra*t->ndim; + cnumber nfgiven = t->ngiven*t->ncomp; + cnumber nfextra = t->nextra*t->ncomp; + + Alloc(xgiven, nxgiven + nxextra + nfgiven + nfextra); + t->xextra = xgiven + nxgiven; + fgiven = t->xextra + nxextra; + t->fextra = fgiven + nfgiven; + + if( nxgiven ) { +#ifdef MLVERSION + Copy(xgiven, t->xgiven, nxgiven); + Copy(fgiven, t->fgiven, nfgiven); +#else + if( t->ldxgiven == t->ndim ) + Copy(xgiven, t->xgiven, nxgiven); + else { + number i; + real *sgiven = t->xgiven, *dgiven = xgiven; + for( i = 0; i < t->ngiven; ++i ) { + Copy(dgiven, sgiven, t->ndim); + sgiven += t->ldxgiven; + dgiven += t->ndim; + } + } + t->phase = 0; + DoSample(t, t->ngiven, xgiven, fgiven); +#endif + } + } + + t->xgiven = xgiven; + t->fgiven = fgiven; +} diff --git a/src/external/libCuba/src/divonne/decl.h b/src/external/libCuba/src/divonne/decl.h index 8fc714cc..cf127851 100644 --- a/src/external/libCuba/src/divonne/decl.h +++ b/src/external/libCuba/src/divonne/decl.h @@ -2,32 +2,16 @@ decl.h Type declarations this file is part of Divonne - last modified 8 Jun 10 th + last modified 26 Jul 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - #include "stddecl.h" +#define INIDEPTH 3 +#define DEPTH 5 +#define POSTDEPTH 15 + #define Tag(x) ((x) | INT_MIN) #define Untag(x) ((x) & INT_MAX) @@ -37,15 +21,31 @@ typedef struct { typedef const Bounds cBounds; +typedef struct { + real avg, err; +} PhaseResult; + typedef struct { real avg, spreadsq; real spread, secondspread; real nneed, maxerrsq, mindevsq; + real integral, sigsq, chisq; + PhaseResult phase[2]; int iregion; } Totals; +enum { nrules = 5 }; + typedef struct { - void *first, *last; + count n; + real weight[nrules], scale[nrules], norm[nrules]; + real gen[]; +} Set; + +#define SetSize (sizeof(Set) + t->ndim*sizeof(real)) + +typedef struct { + Set *first, *last; real errcoeff[3]; count n; } Rule; @@ -53,16 +53,45 @@ typedef struct { typedef const Rule cRule; typedef struct samples { - real weight; - real *x, *f, *avg, *err; - void (*sampler)(struct _this *t, const struct samples *, cBounds *, creal); + real *x, *f; + void (*sampler)(struct _this *t, ccount); cRule *rule; - count coeff; number n, neff; + count coeff; } Samples; typedef const Samples cSamples; +typedef struct { + real diff, err, spread; +} Errors; + +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)) + +typedef struct region { + int depth, next; + count isamples, cutcomp, xmajor; + real fmajor, fminor, vol; + Bounds bounds[]; +} Region; + +#define RegionSize (sizeof(Region) + t->ndim*sizeof(Bounds) + t->ncomp*ResultSize) + +#define RegionResult(r) ((Result *)(r->bounds + t->ndim)) + +#define RegionPtr(n) ((Region *)((char *)t->region + (n)*regionsize)) + + typedef int (*Integrand)(ccount *, creal *, ccount *, real *, void *, cint *); typedef void (*PeakFinder)(ccount *, cBounds *, number *, real *); @@ -73,6 +102,12 @@ typedef struct _this { Integrand integrand; void *userdata; PeakFinder peakfinder; +#ifdef HAVE_FORK + int ncores, running, *child; + real *frame; + number nframe; + SHM_ONLY(int shmid;) +#endif #endif real epsrel, epsabs; int flags, seed; @@ -82,37 +117,39 @@ typedef struct _this { Bounds border; real maxchisq, mindeviation; number ngiven, nextra; - real *xgiven, *xextra, *fgiven, *fextra; + real *xgiven, *fgiven; + real *xextra, *fextra; count ldxgiven; count nregions; - number neval, neval_opt, neval_cut; - int phase; + cchar *statefile; + number neval, neval_opt, neval_cut, nrand; + count phase; count selectedcomp, size; Samples samples[3]; Totals *totals; Rule rule7, rule9, rule11, rule13; RNGState rng; - void *voidregion; + Region *region; jmp_buf abort; } This; typedef const This cThis; + #define CHUNKSIZE 4096 -#define TYPEDEFREGION \ - typedef struct { \ - real avg, err, spread, chisq; \ - real fmin, fmax; \ - real xmin[NDIM], xmax[NDIM]; \ - } Result; \ - typedef const Result cResult; \ - typedef struct region { \ - count cutcomp, depth, xmajor; \ - real fmajor, fminor, vol; \ - Bounds bounds[NDIM]; \ - Result result[NCOMP]; \ - } Region +#define AllocRegions(t) \ + MemAlloc((t)->region, (t)->size*regionsize) -#define RegionPtr(n) (&((Region *)t->voidregion)[n]) +#define EnlargeRegions(t, n) if( (t)->nregions + n > (t)->size ) \ + ReAlloc((t)->region, ((t)->size += CHUNKSIZE)*regionsize) + +#define SAMPLERDEFS \ + csize_t regionsize = RegionSize; \ + Region *region = RegionPtr(iregion); \ + cBounds *b = region->bounds; \ + Result *res = RegionResult(region); \ + cSamples *samples = &t->samples[region->isamples]; \ + real *x = samples->x, *f = samples->f; \ + cnumber n = samples->n diff --git a/src/external/libCuba/src/divonne/util.c b/src/external/libCuba/src/divonne/util.c deleted file mode 100644 index dd470281..00000000 --- a/src/external/libCuba/src/divonne/util.c +++ /dev/null @@ -1,51 +0,0 @@ -/* - util.c - Utility functions - this file is part of Divonne - last modified 9 Apr 09 th -*/ - -/*************************************************************************** - * Copyright (C) 2004-2009 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - -#include "decl.h" - -static count ndim_, ncomp_, selectedcomp_, nregions_; -static number neval_, neval_opt_, neval_cut_; -static int sign_, phase_; - -static Bounds border_; - -static Samples samples_[3]; -static Rule rule7_, rule9_, rule11_, rule13_; -static real *xgiven_, *fgiven_, *xextra_, *fextra_; -static count ldxgiven_; -static number ngiven_, nextra_; - -static Totals *totals_; - -static void *voidregion_; -#define region_ ((Region *)voidregion_) -static count size_; - -#ifdef DEBUG -#include "debug.c" -#endif - diff --git a/src/external/libCuba/src/suave/Fluct.c b/src/external/libCuba/src/suave/Fluct.c index 37ba0574..ad1180ab 100644 --- a/src/external/libCuba/src/suave/Fluct.c +++ b/src/external/libCuba/src/suave/Fluct.c @@ -2,48 +2,32 @@ Fluct.c compute the fluctuation in the left and right half this file is part of Suave - last modified 9 Feb 05 th + last modified 29 Jul 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - #if defined(HAVE_LONG_DOUBLE) && defined(HAVE_POWL) -typedef long double xdouble; +typedef long double realx; #define XDBL_MAX_EXP LDBL_MAX_EXP #define XDBL_MAX LDBL_MAX #define powx powl +#define ldexpx ldexpl #else -typedef double xdouble; +typedef double realx; #define XDBL_MAX_EXP DBL_MAX_EXP #define XDBL_MAX DBL_MAX #define powx pow +#define ldexpx ldexp #endif +typedef const realx crealx; + typedef struct { - xdouble fluct; + realx fluct; number n; } Var; @@ -54,23 +38,23 @@ static void Fluct(cThis *t, Var *var, { creal *x = w + n; creal *f = x + n*t->ndim + comp; - creal flat = 2/3./t->flatness; - creal max = ldexp(1., (int)((XDBL_MAX_EXP - 2)/t->flatness)); - creal norm = 1/(err*Max(fabs(avg), err)); 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)); Clear(var, nvar); while( n-- ) { count dim; - const xdouble ft = - powx(Min(1 + fabs(*w++)*Sq(*f - avg)*norm, max), t->flatness); + crealx arg = 1 + fabs(*w++)*Sq(*f - avg)*norm; + crealx ft = powx(arg < max ? arg : max, t->flatness); f += t->ncomp; for( dim = 0; dim < t->ndim; ++dim ) { - Var *v = &var[2*dim + (*x++ >= b[dim].mid)]; - const xdouble f = v->fluct + ft; + Var *v = &var[2*dim + (*x++ >= .5*(b[dim].lower + b[dim].upper))]; + crealx f = v->fluct + ft; v->fluct = (f > XDBL_MAX/2) ? XDBL_MAX/2 : f; ++v->n; } diff --git a/src/external/libCuba/src/suave/Grid.c b/src/external/libCuba/src/suave/Grid.c index 344040ae..8e52fde8 100644 --- a/src/external/libCuba/src/suave/Grid.c +++ b/src/external/libCuba/src/suave/Grid.c @@ -2,28 +2,9 @@ Grid.c utility functions for the Vegas grid this file is part of Suave - last modified 1 Jun 10 th + last modified 7 Aug 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ static void RefineGrid(cThis *t, Grid grid, Grid margsum) { @@ -72,7 +53,7 @@ static void RefineGrid(cThis *t, Grid grid, Grid margsum) delta = (cur - prev)*thisbin; newgrid[newbin] = SHARPEDGES ? cur - delta/imp[bin] : - (newcur = Max(newcur + 0x1p-48, + (newcur = Max(newcur + 16*DBL_EPSILON, cur - 2*delta/(imp[bin] + imp[IDim(bin - 1)]))); } Copy(grid, newgrid, NBINS - 1); @@ -84,8 +65,8 @@ static void RefineGrid(cThis *t, Grid grid, Grid margsum) static void Reweight(cThis *t, Bounds *b, creal *w, creal *f, creal *lastf, cResult *total) { - Grid margsum[NDIM]; - real scale[NCOMP]; + Vector(Grid, margsum, NDIM); + Vector(real, scale, NCOMP); cbin_t *bin = (cbin_t *)lastf; count dim, comp; @@ -95,7 +76,7 @@ static void Reweight(cThis *t, Bounds *b, scale[comp] = (total[comp].avg == 0) ? 0 : 1/total[comp].avg; } - Zap(margsum); + XClear(margsum); while( f < lastf ) { real fsq = 0; diff --git a/src/external/libCuba/src/suave/Integrate.c b/src/external/libCuba/src/suave/Integrate.c index b0cf3027..9c06d81f 100644 --- a/src/external/libCuba/src/suave/Integrate.c +++ b/src/external/libCuba/src/suave/Integrate.c @@ -2,90 +2,122 @@ Integrate.c integrate over the unit hypercube this file is part of Suave - last modified 16 Jun 10 th + checkpointing by B. Chokoufe + last modified 5 Aug 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ +typedef struct { + signature_t signature; + count nregions, df; + number neval; + Result totals[]; +} State; static int Integrate(This *t, real *integral, real *error, real *prob) { - TYPEDEFREGION; + StateDecl; + csize_t statesize = sizeof(State) + NCOMP*sizeof(Result); + Sized(State, state, statesize); + Array(Var, var, NDIM, 2); + Vector(char, out, 128*NCOMP + 256); - count dim, comp, df; - int fail = -99; - Result totals[NCOMP]; Region *anchor = NULL, *region = NULL; + Result *tot, *Tot = state->totals + t->ncomp; + Result *res, *resL, *resR; + Bounds *b, *B; + count dim, comp; + int fail; if( VERBOSE > 1 ) { - char s[256]; - sprintf(s, "Suave input parameters:\n" + sprintf(out, "Suave input parameters:\n" " ndim " COUNT "\n ncomp " COUNT "\n" " epsrel " REAL "\n epsabs " REAL "\n" " flags %d\n seed %d\n" " mineval " NUMBER "\n maxeval " NUMBER "\n" - " nnew " NUMBER "\n flatness " REAL, + " nnew " NUMBER "\n flatness " REAL "\n" + " statefile \"%s\"\n", t->ndim, t->ncomp, t->epsrel, t->epsabs, t->flags, t->seed, t->mineval, t->maxeval, - t->nnew, t->flatness); - Print(s); + t->nnew, t->flatness, + t->statefile); + Print(out); } if( BadComponent(t) ) return -2; if( BadDimension(t) ) return -1; - if( setjmp(t->abort) ) goto abort; + ShmAlloc(t, ShmRm(t)); + ForkCores(t); + + if( (fail = setjmp(t->abort)) ) goto abort; t->epsabs = Max(t->epsabs, NOTZERO); IniRandom(t); - RegionAlloc(anchor, t->nnew, t->nnew); - anchor->next = NULL; - anchor->div = 0; + StateSetup(t); - for( dim = 0; dim < t->ndim; ++dim ) { - Bounds *b = &anchor->bounds[dim]; - b->lower = 0; - b->upper = 1; - b->mid = .5; - - if( dim == 0 ) { - count bin; - /* define the initial distribution of bins */ - for( bin = 0; bin < NBINS; ++bin ) - b->grid[bin] = (bin + 1)/(real)NBINS; - } - else Copy(b->grid, anchor->bounds[0].grid, NBINS); + if( StateReadTest(t) ) { + size_t *regsize = NULL; + StateReadOpen(t, fd) { + count iregion; + size_t totsize; + Region **last = &anchor; + if( read(fd, state, statesize) != statesize || + state->signature != StateSignature(t, 2) ) break; + t->nregions = state->nregions; + totsize = t->nregions*sizeof *regsize; + MemAlloc(regsize, totsize); + StateRead(fd, regsize, totsize); + for( iregion = 0; iregion < t->nregions; ++iregion ) + totsize += regsize[iregion]; + if( st.st_size != statesize + totsize ) break; + for( iregion = 0; iregion < t->nregions; ++iregion ) { + Region *reg; + MemAlloc(reg, regsize[iregion]); + StateRead(fd, reg, regsize[iregion]); + *last = reg; + last = ®->next; + } + *last = NULL; + } StateReadClose(t, fd); + free(regsize); + t->neval = state->neval; + t->rng.skiprandom(t, t->neval); } - Sample(t, t->nnew, anchor, anchor->w, - anchor->w + t->nnew, - anchor->w + t->nnew + t->ndim*t->nnew); - df = anchor->df; - ResCopy(totals, anchor->result); + if( ini ) { + count bin; + Bounds *b0; - for( t->nregions = 1; ; ++t->nregions ) { - Var var[NDIM][2], *vLR; + anchor = RegionAlloc(t, t->nnew, t->nnew); + anchor->next = NULL; + anchor->div = 0; + t->nregions = 1; + + b0 = RegionBounds(anchor); + b0->lower = 0; + b0->upper = 1; + /* define the initial distribution of bins */ + for( bin = 0; bin < NBINS; ++bin ) + b0->grid[bin] = (bin + 1)/(real)NBINS; + + for( b = b0 + 1, B = b0 + t->ndim; b < B; ++b ) + Copy(b, b0, 1); + + t->neval = 0; + Sample(t, t->nnew, anchor, RegionW(anchor), + RegionW(anchor) + t->nnew, + RegionW(anchor) + t->nnew + t->ndim*t->nnew); + state->df = anchor->df; + FCopy(state->totals, anchor->result); + } + + /* main iteration loop */ + for( ; ; ) { + Var *vLR; real maxratio, maxerr, minfluct, bias, mid; Region *regionL, *regionR, *reg, **parent, **par; Bounds *bounds, *boundsL, *boundsR; @@ -94,39 +126,33 @@ static int Integrate(This *t, real *integral, real *error, real *prob) real *w, *wL, *wR, *x, *xL, *xR, *f, *fL, *fR, *wlast, *flast; if( VERBOSE ) { - char s[128 + 128*NCOMP], *p = s; - - p += sprintf(p, "\n" + char *oe = out + sprintf(out, "\n" "Iteration " COUNT ": " NUMBER " integrand evaluations so far", t->nregions, t->neval); - - for( comp = 0; comp < t->ncomp; ++comp ) { - cResult *tot = &totals[comp]; - p += sprintf(p, "\n[" COUNT "] " + for( tot = state->totals, comp = 0; tot < Tot; ++tot ) + oe += sprintf(oe, "\n[" COUNT "] " REAL " +- " REAL " \tchisq " REAL " (" COUNT " df)", - comp + 1, tot->avg, tot->err, tot->chisq, df); - } - - Print(s); + ++comp, tot->avg, tot->err, tot->chisq, state->df); + Print(out); } maxratio = -INFTY; maxcomp = 0; - for( comp = 0; comp < t->ncomp; ++comp ) { - creal ratio = totals[comp].err/MaxErr(totals[comp].avg); + for( tot = state->totals, comp = 0; tot < Tot; ++tot ) { + creal ratio = tot->err/MaxErr(tot->avg); if( ratio > maxratio ) { maxratio = ratio; maxcomp = comp; } } - if( maxratio <= 1 && t->neval >= t->mineval ) { - fail = 0; + if( maxratio <= 1 && t->neval >= t->mineval ) break; + + if( t->neval >= t->maxeval ) { + fail = 1; break; } - if( t->neval >= t->maxeval ) break; - maxerr = -INFTY; parent = &anchor; region = anchor; @@ -139,8 +165,11 @@ static int Integrate(This *t, real *integral, real *error, real *prob) } } - Fluct(t, var[0], - region->bounds, region->w, region->n, maxcomp, + bounds = RegionBounds(region); + w = RegionW(region); + + /* find the bisectdim which minimizes the fluctuations */ + Fluct(t, var[0], bounds, w, region->n, maxcomp, region->result[maxcomp].avg, Max(maxerr, t->epsabs)); bias = (t->epsrel < 1e-50) ? 2 : @@ -148,15 +177,15 @@ static int Integrate(This *t, real *integral, real *error, real *prob) minfluct = INFTY; bisectdim = 0; for( dim = 0; dim < t->ndim; ++dim ) { - cBounds *b = ®ion->bounds[dim]; creal fluct = (var[dim][0].fluct + var[dim][1].fluct)* - (bias - b->upper + b->lower); + (bias - bounds[dim].upper + bounds[dim].lower); if( fluct < minfluct ) { minfluct = fluct; bisectdim = dim; } } + /* apply stratified sampling to distribute points in bisected region */ vLR = var[bisectdim]; minfluct = vLR[0].fluct + vLR[1].fluct; nnewL = IMax( @@ -166,37 +195,35 @@ static int Integrate(This *t, real *integral, real *error, real *prob) nnewR = IMax(t->nnew - nnewL, MINSAMPLES); nR = vLR[1].n + nnewR; - RegionAlloc(regionL, nL, nnewL); - RegionAlloc(regionR, nR, nnewR); + regionL = RegionAlloc(t, nL, nnewL); + regionR = RegionAlloc(t, nR, nnewR); *parent = regionL; regionL->next = regionR; regionR->next = region->next; regionL->div = regionR->div = region->div + 1; - bounds = ®ion->bounds[bisectdim]; - mid = bounds->mid; + mid = .5*(bounds[bisectdim].lower + bounds[bisectdim].upper); n = region->n; - w = wlast = region->w; x = w + n; f = flast = x + n*t->ndim; - wL = regionL->w; xL = wL + nL; fL = xL + nL*t->ndim; - wR = regionR->w; xR = wR + nR; fR = xR + nR*t->ndim; - + wlast = w; x = w + n; f = flast = x + n*t->ndim; + wL = RegionW(regionL); xL = wL + nL; fL = xL + nL*t->ndim; + wR = RegionW(regionR); xR = wR + nR; fR = xR + nR*t->ndim; while( n-- ) { cbool final = (*w < 0); if( x[bisectdim] < mid ) { - if( final && wR > regionR->w ) *(wR - 1) = -fabs(*(wR - 1)); + if( final && wR > RegionW(regionR) ) wR[-1] = -fabs(wR[-1]); *wL++ = *w++; - VecCopy(xL, x); + XCopy(xL, x); xL += t->ndim; - ResCopy(fL, f); + FCopy(fL, f); fL += t->ncomp; } else { - if( final && wL > regionL->w ) *(wL - 1) = -fabs(*(wL - 1)); + if( final && wL > RegionW(regionL) ) wL[-1] = -fabs(wL[-1]); *wR++ = *w++; - VecCopy(xR, x); + XCopy(xR, x); xR += t->ndim; - ResCopy(fR, f); + FCopy(fR, f); fR += t->ncomp; } x += t->ndim; @@ -204,80 +231,100 @@ static int Integrate(This *t, real *integral, real *error, real *prob) if( n && final ) wlast = w, flast = f; } - Reweight(t, region->bounds, wlast, flast, f, totals); - VecCopy(regionL->bounds, region->bounds); - VecCopy(regionR->bounds, region->bounds); + Reweight(t, bounds, wlast, flast, f, state->totals); + boundsL = RegionBounds(regionL); + XCopy(boundsL, bounds); + boundsR = RegionBounds(regionR); + XCopy(boundsR, bounds); - boundsL = ®ionL->bounds[bisectdim]; - boundsR = ®ionR->bounds[bisectdim]; - boundsL->mid = .5*(boundsL->lower + (boundsL->upper = mid)); - boundsR->mid = .5*((boundsR->lower = mid) + boundsR->upper); - - StretchGrid(bounds->grid, boundsL->grid, boundsR->grid); + boundsL[bisectdim].upper = mid; + boundsR[bisectdim].lower = mid; + StretchGrid(bounds[bisectdim].grid, + boundsL[bisectdim].grid, boundsR[bisectdim].grid); Sample(t, nnewL, regionL, wL, xL, fL); Sample(t, nnewR, regionR, wR, xR, fR); - df += regionL->df + regionR->df - region->df; + state->df += regionL->df + regionR->df - region->df; - for( comp = 0; comp < t->ncomp; ++comp ) { - cResult *r = ®ion->result[comp]; - Result *rL = ®ionL->result[comp]; - Result *rR = ®ionR->result[comp]; - Result *tot = &totals[comp]; + for( res = region->result, + resL = regionL->result, + resR = regionR->result, + tot = state->totals; + tot < Tot; ++res, ++resL, ++resR, ++tot ) { real diff, sigsq; - tot->avg += diff = rL->avg + rR->avg - r->avg; + tot->avg += diff = resL->avg + resR->avg - res->avg; diff = Sq(.25*diff); - sigsq = rL->sigsq + rR->sigsq; + sigsq = resL->sigsq + resR->sigsq; if( sigsq > 0 ) { creal c = Sq(1 + sqrt(diff/sigsq)); - rL->sigsq *= c; - rR->sigsq *= c; + resL->sigsq *= c; + resR->sigsq *= c; } - rL->err = sqrt(rL->sigsq += diff); - rR->err = sqrt(rR->sigsq += diff); + resL->err = sqrt(resL->sigsq += diff); + resR->err = sqrt(resR->sigsq += diff); - tot->sigsq += rL->sigsq + rR->sigsq - r->sigsq; + tot->sigsq += resL->sigsq + resR->sigsq - res->sigsq; tot->err = sqrt(tot->sigsq); - tot->chisq += rL->chisq + rR->chisq - r->chisq; + tot->chisq += resL->chisq + resR->chisq - res->chisq; } free(region); region = NULL; + ++t->nregions; + + if( StateWriteTest(t) ) { + StateWriteOpen(t, fd) { + Region *reg; + size_t totsize, *regsize, *s; + state->signature = StateSignature(t, 2); + state->neval = t->neval; + state->nregions = t->nregions; + StateWrite(fd, state, statesize); + MemAlloc(regsize, totsize = t->nregions*sizeof(size_t)); + s = regsize; + for( reg = anchor; reg; reg = reg->next ) + *s++ = reg->size; + StateWrite(fd, regsize, totsize); + free(regsize); + for( reg = anchor; reg; reg = reg->next ) + StateWrite(fd, reg, reg->size); + } StateWriteClose(t, fd); + } } - for( comp = 0; comp < t->ncomp; ++comp ) { - cResult *tot = &totals[comp]; + for( tot = state->totals, comp = 0; tot < Tot; ++tot, ++comp ) { integral[comp] = tot->avg; error[comp] = tot->err; - prob[comp] = ChiSquare(tot->chisq, df); + prob[comp] = ChiSquare(tot->chisq, state->df); } #ifdef MLVERSION if( REGIONS ) { + Vector(real, bounds, NDIM*2); + MLPutFunction(stdlink, "List", 2); + MLPutFunction(stdlink, "List", t->nregions); for( region = anchor; region; region = region->next ) { - real lower[NDIM], upper[NDIM]; - - for( dim = 0; dim < t->ndim; ++dim ) { - cBounds *b = ®ion->bounds[dim]; - lower[dim] = b->lower; - upper[dim] = b->upper; + cResult *Res; + real *d = bounds; + for( B = (b = RegionBounds(region)) + t->ndim; b < B; ++b ) { + *d++ = b->lower; + *d++ = b->upper; } - MLPutFunction(stdlink, "Cuba`Suave`region", 4); - MLPutRealList(stdlink, lower, t->ndim); - MLPutRealList(stdlink, upper, t->ndim); + MLPutFunction(stdlink, "Cuba`Suave`region", 3); + + MLPutRealList(stdlink, bounds, 2*t->ndim); MLPutFunction(stdlink, "List", t->ncomp); - for( comp = 0; comp < t->ncomp; ++comp ) { - cResult *r = ®ion->result[comp]; - real res[] = {r->avg, r->err, r->chisq}; - MLPutRealList(stdlink, res, Elements(res)); + for( Res = (res = region->result) + t->ncomp; res < Res; ++res ) { + real r[] = {res->avg, res->err, res->chisq}; + MLPutRealList(stdlink, r, Elements(r)); } MLPutInteger(stdlink, region->df); @@ -287,11 +334,14 @@ static int Integrate(This *t, real *integral, real *error, real *prob) abort: free(region); - while( (region = anchor) ) { anchor = anchor->next; free(region); } + WaitCores(t); + ShmFree(t); + + StateRemove(t); return fail; } diff --git a/src/external/libCuba/src/suave/Makefile.am b/src/external/libCuba/src/suave/Makefile.am new file mode 100644 index 00000000..22fbfaea --- /dev/null +++ b/src/external/libCuba/src/suave/Makefile.am @@ -0,0 +1,13 @@ +## Process this file with automake to create Makefile.in +## $Id: Makefile.am 4808 2011-03-21 14:33:34Z l_wojek $ + +c_sources = Suave.c + +AM_CPPFLAGS = -I. -I.. -I../common -DNOUNDERSCORE +AM_CFLAGS = $(LOCAL_CUBA_LIB_CFLAGS) +AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS) + +noinst_LTLIBRARIES = libsuave.la + +libsuave_la_SOURCES = $(c_sources) +libsuave_la_LDFLAGS = $(AM_LDFLAGS) \ No newline at end of file diff --git a/src/external/libCuba/src/suave/Sample.c b/src/external/libCuba/src/suave/Sample.c index 9e5cdaa0..624c816c 100644 --- a/src/external/libCuba/src/suave/Sample.c +++ b/src/external/libCuba/src/suave/Sample.c @@ -2,29 +2,9 @@ Sample.c the sampling step of Suave this file is part of Suave - last modified 13 Sep 10 th + last modified 30 Jul 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - typedef struct { real sum, sqsum; @@ -34,16 +14,16 @@ typedef struct { /*********************************************************************/ -static void Sample(This *t, cnumber nnew, void *voidregion, +static void Sample(This *t, cnumber nnew, Region *region, real *lastw, real *lastx, real *lastf) { - TYPEDEFREGION; - - Region *const region = (Region *)voidregion; - count comp, dim, df; + count comp, df; number n; - Cumulants cumul[NCOMP]; - char **ss, *s; + Vector(Cumulants, cumul, NCOMP); + Cumulants *C = cumul + t->ncomp, *c; + Bounds *bounds = RegionBounds(region), *B = bounds + t->ndim, *b; + Result *res; + char **ss = NULL, *s = NULL; ccount chars = 128*(region->div + 1); creal jacobian = 1/ldexp((real)nnew, region->div); @@ -55,8 +35,7 @@ static void Sample(This *t, cnumber nnew, void *voidregion, t->rng.getrandom(t, f); - for( dim = 0; dim < t->ndim; ++dim ) { - cBounds *b = ®ion->bounds[dim]; + for( b = bounds; b < B; ++b ) { creal pos = *f*NBINS; ccount ipos = (count)pos; creal prev = (ipos == 0) ? 0 : b->grid[ipos - 1]; @@ -69,11 +48,11 @@ static void Sample(This *t, cnumber nnew, void *voidregion, *w++ = weight; } - DoSample(t, nnew, lastw, lastx, lastf, region->div + 1); + DoSample(t, nnew, lastx, lastf, lastw, region->div + 1); w[-1] = -w[-1]; lastw = w; - w = region->w; + w = RegionW(region); region->n = lastw - w; if( VERBOSE > 2 ) { @@ -87,7 +66,7 @@ static void Sample(This *t, cnumber nnew, void *voidregion, } } - Zap(cumul); + FClear(cumul); df = n = 0; while( w < lastw ) { @@ -95,9 +74,7 @@ static void Sample(This *t, cnumber nnew, void *voidregion, creal weight = fabs(*w++); ++n; - for( comp = 0; comp < t->ncomp; ++comp ) { - Cumulants *c = &cumul[comp]; - + for( c = cumul, comp = 0; c < C; ++c ) { creal wfun = weight*(*f++); c->sum += wfun; c->sqsum += Sq(wfun); @@ -135,23 +112,21 @@ static void Sample(This *t, cnumber nnew, void *voidregion, region->df = --df; - for( comp = 0; comp < t->ncomp; ++comp ) { - Result *r = ®ion->result[comp]; - Cumulants *c = &cumul[comp]; + for( c = cumul, res = region->result; c < C; ++c, ++res ) { creal sigsq = 1/c->weightsum; creal avg = sigsq*c->avgsum; if( LAST ) { - r->sigsq = 1/c->weight; - r->avg = r->sigsq*c->avg; + res->sigsq = 1/c->weight; + res->avg = res->sigsq*c->avg; } else { - r->sigsq = sigsq; - r->avg = avg; + res->sigsq = sigsq; + res->avg = avg; } - r->err = sqrt(r->sigsq); + res->err = sqrt(res->sigsq); - r->chisq = (sigsq < .9*NOTZERO) ? 0 : c->chisqsum - avg*c->chisum; + res->chisq = (sigsq < .9*NOTZERO) ? 0 : c->chisqsum - avg*c->chisum; /* This catches the special case where the integrand is constant over the entire region. Unless that constant is zero, only the first set of samples will have zero variance, and hence weight @@ -168,19 +143,17 @@ static void Sample(This *t, cnumber nnew, void *voidregion, if( VERBOSE > 2 ) { char *p = s; char *p0 = p + t->ndim*64; + char *msg = "\nRegion (" REALF ") - (" REALF ")"; - for( dim = 0; dim < t->ndim; ++dim ) { - cBounds *b = ®ion->bounds[dim]; - p += sprintf(p, - (dim == 0) ? "\nRegion (" REALF ") - (" REALF ")" : - "\n (" REALF ") - (" REALF ")", - b->lower, b->upper); + for( b = bounds; b < B; ++b ) { + p += sprintf(p, msg, b->lower, b->upper); + msg = "\n (" REALF ") - (" REALF ")"; } - for( comp = 0; comp < t->ncomp; ++comp ) { - cResult *r = ®ion->result[comp]; + for( comp = 0, res = region->result; + comp < t->ncomp; ++comp, ++res ) { p += sprintf(p, "%s \tchisq " REAL " (" COUNT " df)", - p0, r->chisq, df); + p0, res->chisq, df); p0 += chars; } diff --git a/src/external/libCuba/src/suave/Suave.c b/src/external/libCuba/src/suave/Suave.c index 88d98f10..79323c14 100644 --- a/src/external/libCuba/src/suave/Suave.c +++ b/src/external/libCuba/src/suave/Suave.c @@ -2,59 +2,25 @@ Suave.c Subregion-adaptive Vegas Monte-Carlo integration by Thomas Hahn - last modified 13 Sep 10 th + last modified 17 Sep 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ +#define SUAVE +#define ROUTINE "Suave" #include "decl.h" - -#define Print(s) puts(s); fflush(stdout) +#include "CSample.c" /*********************************************************************/ -static inline void DoSample(This *t, number n, - creal *w, creal *x, real *f, cint iter) -{ - t->neval += n; - while( n-- ) { - if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata, - w++, &iter) == ABORT ) - longjmp(t->abort, 1); - x += t->ndim; - f += t->ncomp; - } -} - -/*********************************************************************/ - -#include "common.c" - Extern void EXPORT(Suave)(ccount ndim, ccount ncomp, Integrand integrand, void *userdata, creal epsrel, creal epsabs, cint flags, cint seed, cnumber mineval, cnumber maxeval, cnumber nnew, creal flatness, + cchar *statefile, count *pnregions, number *pneval, int *pfail, real *integral, real *error, real *prob) { @@ -71,8 +37,7 @@ Extern void EXPORT(Suave)(ccount ndim, ccount ncomp, t.maxeval = maxeval; t.nnew = nnew; t.flatness = flatness; - t.nregions = 0; - t.neval = 0; + t.statefile = statefile; *pfail = Integrate(&t, integral, error, prob); *pnregions = t.nregions; @@ -87,8 +52,9 @@ Extern void EXPORT(suave)(ccount *pndim, ccount *pncomp, cint *pflags, cint *pseed, cnumber *pmineval, cnumber *pmaxeval, cnumber *pnnew, creal *pflatness, + cchar *statefile, count *pnregions, number *pneval, int *pfail, - real *integral, real *error, real *prob) + real *integral, real *error, real *prob, cint statefilelen) { This t; t.ndim = *pndim; @@ -103,8 +69,7 @@ Extern void EXPORT(suave)(ccount *pndim, ccount *pncomp, t.maxeval = *pmaxeval; t.nnew = *pnnew; t.flatness = *pflatness; - t.nregions = 0; - t.neval = 0; + CString(t.statefile, statefile, statefilelen); *pfail = Integrate(&t, integral, error, prob); *pnregions = t.nregions; diff --git a/src/external/libCuba/src/suave/common.c b/src/external/libCuba/src/suave/common.c index 3e9f3210..163f1116 100644 --- a/src/external/libCuba/src/suave/common.c +++ b/src/external/libCuba/src/suave/common.c @@ -2,45 +2,34 @@ common.c includes most of the modules this file is part of Suave - last modified 2 Jun 10 th + last modified 29 Jul 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - -#define RegionAlloc(p, n, nnew) MemAlloc(p, \ - sizeof(Region) + \ - (n)*(t->ndim + t->ncomp + 1)*sizeof(real) + \ - (nnew)*t->ndim*sizeof(bin_t)) +static inline Region *RegionAlloc(cThis *t, cnumber n, cnumber nnew) +{ + csize_t size = sizeof(Region) + + t->ncomp*sizeof(Result) + + t->ndim*sizeof(Bounds) + + t->ncomp*t->ndim*2*sizeof(real) + + n*SAMPLESIZE + + nnew*t->ndim*sizeof(bin_t); + Region *p; + MemAlloc(p, size); + p->size = size; + return p; +} static inline bool BadDimension(cThis *t) { - if( t->ndim > NDIM ) return true; + if( t->ndim > MAXDIM ) return true; return t->ndim < SOBOL_MINDIM || (t->seed == 0 && t->ndim > SOBOL_MAXDIM); } static inline bool BadComponent(cThis *t) { - if( t->ncomp > NCOMP ) return true; + if( t->ncomp > MAXCOMP ) return true; return t->ncomp < 1; } @@ -49,5 +38,4 @@ static inline bool BadComponent(cThis *t) #include "Grid.c" #include "Sample.c" #include "Fluct.c" -#include "Integrate.c" diff --git a/src/external/libCuba/src/suave/decl.h b/src/external/libCuba/src/suave/decl.h index 02fec798..8944e60f 100644 --- a/src/external/libCuba/src/suave/decl.h +++ b/src/external/libCuba/src/suave/decl.h @@ -2,29 +2,9 @@ decl.h Type declarations this file is part of Suave - last modified 13 Sep 10 th + last modified 29 Jul 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - #include "stddecl.h" @@ -48,7 +28,7 @@ typedef struct { typedef const Result cResult; typedef struct { - real lower, upper, mid; + real lower, upper; Grid grid; } Bounds; @@ -62,28 +42,36 @@ typedef struct _this { #ifndef MLVERSION Integrand integrand; void *userdata; +#ifdef HAVE_FORK + int ncores, *child; + real *frame; + SHM_ONLY(int shmid;) +#endif #endif real epsrel, epsabs; int flags, seed; number mineval, maxeval; number nnew; real flatness; + cchar *statefile; count nregions; number neval; RNGState rng; jmp_buf abort; } This; +#define nframe nnew + typedef const This cThis; -#define TYPEDEFREGION \ - typedef struct region { \ - struct region *next; \ - count div, df; \ - number n; \ - Result result[NCOMP]; \ - Bounds bounds[NDIM]; \ - real fluct[NCOMP][NDIM][2]; \ - real w[]; \ - } Region +typedef struct region { + struct region *next; + size_t size; + count div, df; + number n; + Result result[]; +} Region; + +#define RegionBounds(r) ((Bounds *)(r->result + t->ncomp)) +#define RegionW(r) ((real *)(RegionBounds(r) + t->ndim)) diff --git a/src/external/libCuba/src/suave/util.c b/src/external/libCuba/src/suave/util.c deleted file mode 100644 index aa6c5c73..00000000 --- a/src/external/libCuba/src/suave/util.c +++ /dev/null @@ -1,43 +0,0 @@ -/* - util.c - Utility functions - this file is part of Suave - last modified 9 Feb 05 th -*/ - -/*************************************************************************** - * Copyright (C) 2004-2009 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - -#include "decl.h" - -static count ndim_, ncomp_, nregions_; -static number neval_; - - -#define RegionAlloc(p, n, nnew) \ - MemAlloc(p, sizeof(Region) + \ - (n)*(ndim_ + ncomp_ + 1)*sizeof(real) + \ - (nnew)*ndim_*sizeof(bin_t)) - - -#ifdef DEBUG -#include "debug.c" -#endif - diff --git a/src/external/libCuba/src/vegas/Grid.c b/src/external/libCuba/src/vegas/Grid.c index 6a90eda4..0ad8a79e 100644 --- a/src/external/libCuba/src/vegas/Grid.c +++ b/src/external/libCuba/src/vegas/Grid.c @@ -2,29 +2,9 @@ Grid.c utility functions for the Vegas grid this file is part of Vegas - last modified 28 May 10 th + last modified 13 Dec 11 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - static inline void GetGrid(cThis *t, Grid *grid) { @@ -35,7 +15,7 @@ static inline void GetGrid(cThis *t, Grid *grid) if( slot < MAXGRIDS && gridptr_[slot] ) { if( griddim_[slot] == t->ndim ) { - VecCopy(grid, gridptr_[slot]); + XCopy(grid, gridptr_[slot]); return; } free(gridptr_[slot]); @@ -57,7 +37,7 @@ static inline void PutGrid(cThis *t, Grid *grid) if( slot < MAXGRIDS ) { if( gridptr_[slot] == NULL ) Alloc(gridptr_[slot], t->ndim); griddim_[slot] = t->ndim; - VecCopy(gridptr_[slot], grid); + XCopy(gridptr_[slot], grid); } } diff --git a/src/external/libCuba/src/vegas/Integrate.c b/src/external/libCuba/src/vegas/Integrate.c index dd5f7e9c..14557871 100644 --- a/src/external/libCuba/src/vegas/Integrate.c +++ b/src/external/libCuba/src/vegas/Integrate.c @@ -2,110 +2,94 @@ Integrate.c integrate over the unit hypercube this file is part of Vegas - last modified 13 Sep 10 th + last modified 8 Aug 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ +typedef struct { + signature_t signature; + count niter; + number nsamples, neval; + Cumulants cumul[]; +} State; static int Integrate(This *t, real *integral, real *error, real *prob) { - real *sample; + bin_t *bins; count dim, comp; - int fail = -99; - struct { - count niter; - number nsamples, neval; - Cumulants cumul[NCOMP]; - Grid grid[NDIM]; - } state; - int statemsg = VERBOSE; - struct stat st; + int fail; + + StateDecl; + csize_t statesize = sizeof(State) + + NCOMP*sizeof(Cumulants) + NDIM*sizeof(Grid); + Sized(State, state, statesize); + Cumulants *c, *C = state->cumul + t->ncomp; + Grid *state_grid = (Grid *)C; + Array(Grid, margsum, NCOMP, NDIM); + Vector(char, out, 128*NCOMP + 256); if( VERBOSE > 1 ) { - char s[512]; - sprintf(s, "Vegas input parameters:\n" + sprintf(out, "Vegas input parameters:\n" " ndim " COUNT "\n ncomp " COUNT "\n" " epsrel " REAL "\n epsabs " REAL "\n" " flags %d\n seed %d\n" " mineval " NUMBER "\n maxeval " NUMBER "\n" " nstart " NUMBER "\n nincrease " NUMBER "\n" " nbatch " NUMBER "\n gridno %d\n" - " statefile \"%s\"\n", + " statefile \"%s\"", t->ndim, t->ncomp, t->epsrel, t->epsabs, t->flags, t->seed, t->mineval, t->maxeval, t->nstart, t->nincrease, t->nbatch, t->gridno, t->statefile); - Print(s); + Print(out); } if( BadComponent(t) ) return -2; if( BadDimension(t) ) return -1; - SamplesAlloc(sample); + FrameAlloc(t, ShmRm(t)); + ForkCores(t); + Alloc(bins, t->nbatch*t->ndim); - if( setjmp(t->abort) ) goto abort; + if( (fail = setjmp(t->abort)) ) goto abort; IniRandom(t); - if( t->statefile && *t->statefile && - stat(t->statefile, &st) == 0 && - st.st_size == sizeof state && (st.st_mode & 0400) ) { - cint h = open(t->statefile, O_RDONLY); - read(h, &state, sizeof state); - close(h); - t->rng.skiprandom(t, t->neval = state.neval); + StateSetup(t); - if( VERBOSE ) { - char s[256]; - sprintf(s, "\nRestoring state from %s.", t->statefile); - Print(s); - } + if( StateReadTest(t) ) { + StateReadOpen(t, fd) { + if( read(fd, state, statesize) != statesize || + state->signature != StateSignature(t, 1) ) break; + } StateReadClose(t, fd); + t->neval = state->neval; + t->rng.skiprandom(t, t->neval); } - else { - t->statefile = NULL; - state.niter = 0; - state.nsamples = t->nstart; - Zap(state.cumul); - GetGrid(t, state.grid); + + if( ini ) { + state->niter = 0; + state->nsamples = t->nstart; + FClear(state->cumul); + GetGrid(t, state_grid); + t->neval = 0; } /* main iteration loop */ - for( ; ; ) { - number nsamples = state.nsamples; + number nsamples = state->nsamples; creal jacobian = 1./nsamples; - Grid margsum[NCOMP][NDIM]; - Zap(margsum); + FClear(margsum); for( ; nsamples > 0; nsamples -= t->nbatch ) { cnumber n = IMin(t->nbatch, nsamples); - real *w = sample; + real *w = t->frame; real *x = w + n; real *f = x + n*t->ndim; real *lastf = f + n*t->ncomp; - bin_t *bin = (bin_t *)lastf; + bin_t *bin = bins; while( x < f ) { real weight = jacobian; @@ -115,8 +99,8 @@ static int Integrate(This *t, real *integral, real *error, real *prob) for( dim = 0; dim < t->ndim; ++dim ) { creal pos = *x*NBINS; ccount ipos = (count)pos; - creal prev = (ipos == 0) ? 0 : state.grid[dim][ipos - 1]; - creal diff = state.grid[dim][ipos] - prev; + creal prev = (ipos == 0) ? 0 : state_grid[dim][ipos - 1]; + creal diff = state_grid[dim][ipos] - prev; *x++ = prev + (pos - ipos)*diff; *bin++ = ipos; weight *= diff*NBINS; @@ -125,25 +109,24 @@ static int Integrate(This *t, real *integral, real *error, real *prob) *w++ = weight; } - DoSample(t, n, sample, w, f, state.niter + 1); + DoSample(t, n, w, f, t->frame, state->niter + 1); - w = sample; - bin = (bin_t *)lastf; + bin = bins; + w = t->frame; while( f < lastf ) { creal weight = *w++; + Grid *m = &margsum[0][0]; - for( comp = 0; comp < t->ncomp; ++comp ) { + for( c = state->cumul; c < C; ++c ) { real wfun = weight*(*f++); if( wfun ) { - Cumulants *c = &state.cumul[comp]; - Grid *m = margsum[comp]; - c->sum += wfun; c->sqsum += wfun *= wfun; for( dim = 0; dim < t->ndim; ++dim ) m[dim][bin[dim]] += wfun; } + m += t->ndim; } bin += t->ndim; @@ -154,19 +137,16 @@ static int Integrate(This *t, real *integral, real *error, real *prob) /* compute the integral and error values */ - for( comp = 0; comp < t->ncomp; ++comp ) { - Cumulants *c = &state.cumul[comp]; - real avg, sigsq; - real w = Weight(c->sum, c->sqsum, state.nsamples); - - sigsq = 1/(c->weightsum += w); - avg = sigsq*(c->avgsum += w*c->sum); + for( c = state->cumul; c < C; ++c ) { + real w = Weight(c->sum, c->sqsum, state->nsamples); + real sigsq = 1/(c->weightsum += w); + real avg = sigsq*(c->avgsum += w*c->sum); c->avg = LAST ? (sigsq = 1/w, c->sum) : avg; c->err = sqrt(sigsq); fail |= (c->err > MaxErr(c->avg)); - if( state.niter == 0 ) c->guess = c->sum; + if( state->niter == 0 ) c->guess = c->sum; else { c->chisum += w *= c->sum - c->guess; c->chisqsum += w*c->sum; @@ -177,38 +157,29 @@ static int Integrate(This *t, real *integral, real *error, real *prob) } if( VERBOSE ) { - char s[128 + 128*NCOMP], *p = s; - - p += sprintf(p, "\n" + char *oe = out + sprintf(out, "\n" "Iteration " COUNT ": " NUMBER " integrand evaluations so far", - state.niter + 1, t->neval); - - for( comp = 0; comp < t->ncomp; ++comp ) { - cCumulants *c = &state.cumul[comp]; - p += sprintf(p, "\n[" COUNT "] " + state->niter + 1, t->neval); + for( c = state->cumul, comp = 0; c < C; ++c ) + oe += sprintf(oe, "\n[" COUNT "] " REAL " +- " REAL " \tchisq " REAL " (" COUNT " df)", - comp + 1, c->avg, c->err, c->chisq, state.niter); - } - - Print(s); + ++comp, c->avg, c->err, c->chisq, state->niter); + Print(out); } - if( fail == 0 && t->neval >= t->mineval ) { - if( t->statefile ) unlink(t->statefile); - break; - } + if( fail == 0 && t->neval >= t->mineval ) break; - if( t->neval >= t->maxeval && t->statefile == NULL ) break; + if( t->neval >= t->maxeval && !StateWriteTest(t) ) break; if( t->ncomp == 1 ) for( dim = 0; dim < t->ndim; ++dim ) - RefineGrid(t, state.grid[dim], margsum[0][dim]); + RefineGrid(t, state_grid[dim], margsum[0][dim]); else { for( dim = 0; dim < t->ndim; ++dim ) { Grid wmargsum; Zap(wmargsum); for( comp = 0; comp < t->ncomp; ++comp ) { - real w = state.cumul[comp].avg; + real w = state->cumul[comp].avg; if( w != 0 ) { creal *m = margsum[comp][dim]; count bin; @@ -217,41 +188,37 @@ static int Integrate(This *t, real *integral, real *error, real *prob) wmargsum[bin] += w*m[bin]; } } - RefineGrid(t, state.grid[dim], wmargsum); + RefineGrid(t, state_grid[dim], wmargsum); } } - ++state.niter; - state.nsamples += t->nincrease; + ++state->niter; + state->nsamples += t->nincrease; - if( t->statefile ) { - cint h = creat(t->statefile, 0666); - if( h != -1 ) { - state.neval = t->neval; - write(h, &state, sizeof state); - close(h); - - if( statemsg ) { - char s[256]; - sprintf(s, "\nSaving state to %s.", t->statefile); - Print(s); - statemsg = false; - } - } + if( StateWriteTest(t) ) { + state->signature = StateSignature(t, 1); + state->neval = t->neval; + StateWriteOpen(t, fd) { + StateWrite(fd, state, statesize); + } StateWriteClose(t, fd); if( t->neval >= t->maxeval ) break; } } for( comp = 0; comp < t->ncomp; ++comp ) { - cCumulants *c = &state.cumul[comp]; + cCumulants *c = &state->cumul[comp]; integral[comp] = c->avg; error[comp] = c->err; - prob[comp] = ChiSquare(c->chisq, state.niter); + prob[comp] = ChiSquare(c->chisq, state->niter); } abort: - free(sample); - PutGrid(t, state.grid); + PutGrid(t, state_grid); + free(bins); + WaitCores(t); + FrameFree(t); + + StateRemove(t); return fail; } diff --git a/src/external/libCuba/src/vegas/Makefile.am b/src/external/libCuba/src/vegas/Makefile.am new file mode 100644 index 00000000..85bbb378 --- /dev/null +++ b/src/external/libCuba/src/vegas/Makefile.am @@ -0,0 +1,13 @@ +## Process this file with automake to create Makefile.in +## $Id: Makefile.am 4808 2011-03-21 14:33:34Z l_wojek $ + +c_sources = Vegas.c + +AM_CPPFLAGS = -I. -I.. -I../common -DNOUNDERSCORE +AM_CFLAGS = $(LOCAL_CUBA_LIB_CFLAGS) +AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS) + +noinst_LTLIBRARIES = libvegas.la + +libvegas_la_SOURCES = $(c_sources) +libvegas_la_LDFLAGS = $(AM_LDFLAGS) \ No newline at end of file diff --git a/src/external/libCuba/src/vegas/Vegas.c b/src/external/libCuba/src/vegas/Vegas.c index 7f193ae1..43a608f7 100644 --- a/src/external/libCuba/src/vegas/Vegas.c +++ b/src/external/libCuba/src/vegas/Vegas.c @@ -2,53 +2,18 @@ Vegas.c Vegas Monte-Carlo integration by Thomas Hahn - last modified 13 Sep 10 th + last modified 17 Sep 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ +#define VEGAS +#define ROUTINE "Vegas" #include "decl.h" - -#define Print(s) puts(s); fflush(stdout) +#include "CSample.c" /*********************************************************************/ -static inline void DoSample(This *t, number n, - creal *w, creal *x, real *f, cint iter) -{ - t->neval += n; - while( n-- ) { - if( t->integrand(&t->ndim, x, &t->ncomp, f, t->userdata, - w++, &iter) == ABORT ) - longjmp(t->abort, 1); - x += t->ndim; - f += t->ncomp; - } -} - -/*********************************************************************/ - -#include "common.c" - Extern void EXPORT(Vegas)(ccount ndim, ccount ncomp, Integrand integrand, void *userdata, creal epsrel, creal epsabs, cint flags, cint seed, @@ -74,7 +39,6 @@ Extern void EXPORT(Vegas)(ccount ndim, ccount ncomp, t.nbatch = nbatch; t.gridno = gridno; t.statefile = statefile; - t.neval = 0; *pfail = Integrate(&t, integral, error, prob); *pneval = t.neval; @@ -89,9 +53,8 @@ Extern void EXPORT(vegas)(ccount *pndim, ccount *pncomp, cnumber *pnstart, cnumber *pnincrease, cnumber *pnbatch, cint *pgridno, cchar *statefile, number *pneval, int *pfail, - real *integral, real *error, real *prob, int statefilelen) + real *integral, real *error, real *prob, cint statefilelen) { - char *s = NULL; This t; t.ndim = *pndim; t.ncomp = *pncomp; @@ -107,22 +70,9 @@ Extern void EXPORT(vegas)(ccount *pndim, ccount *pncomp, t.nincrease = *pnincrease; t.nbatch = *pnbatch; t.gridno = *pgridno; - t.neval = 0; - - if( statefile ) { - /* strip trailing spaces */ - while( statefilelen > 0 && statefile[statefilelen - 1] == ' ' ) - --statefilelen; - if( statefilelen > 0 && (s = malloc(statefilelen + 1)) ) { - memcpy(s, statefile, statefilelen); - s[statefilelen] = 0; - } - } - t.statefile = s; + CString(t.statefile, statefile, statefilelen); *pfail = Integrate(&t, integral, error, prob); *pneval = t.neval; - - free(s); } diff --git a/src/external/libCuba/src/vegas/common.c b/src/external/libCuba/src/vegas/common.c index 1676b26d..eaa59a42 100644 --- a/src/external/libCuba/src/vegas/common.c +++ b/src/external/libCuba/src/vegas/common.c @@ -2,50 +2,24 @@ common.c Code common to Vegas.c and Vegas.tm this file is part of Vegas - last modified 6 Jun 10 th + last modified 29 Jul 13 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - #include "Random.c" #include "ChiSquare.c" #include "Grid.c" -#define SamplesAlloc(p) MemAlloc(p, \ - t->nbatch*((t->ndim + t->ncomp + 1)*sizeof(real) + \ - t->ndim*sizeof(bin_t))) - static inline bool BadDimension(cThis *t) { - if( t->ndim > NDIM ) return true; + if( t->ndim > MAXDIM ) return true; return t->ndim < SOBOL_MINDIM || (t->seed == 0 && t->ndim > SOBOL_MAXDIM); } static inline bool BadComponent(cThis *t) { - if( t->ncomp > NCOMP ) return true; + if( t->ncomp > MAXCOMP ) return true; return t->ncomp < 1; } -#include "Integrate.c" - diff --git a/src/external/libCuba/src/vegas/decl.h b/src/external/libCuba/src/vegas/decl.h index 8365bbbb..0bd2614d 100644 --- a/src/external/libCuba/src/vegas/decl.h +++ b/src/external/libCuba/src/vegas/decl.h @@ -2,29 +2,9 @@ decl.h Type declarations this file is part of Vegas - last modified 13 Sep 10 th + last modified 21 Dec 11 th */ -/*************************************************************************** - * Copyright (C) 2004-2010 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - #include "stddecl.h" @@ -56,7 +36,12 @@ typedef struct _this { #ifndef MLVERSION Integrand integrand; void *userdata; +#ifdef HAVE_FORK + int ncores, *child; + SHM_ONLY(int shmid;) #endif +#endif + real *frame; real epsrel, epsabs; int flags, seed; number mineval, maxeval; @@ -68,6 +53,8 @@ typedef struct _this { jmp_buf abort; } This; +#define nframe nbatch + typedef const This cThis; static Grid *gridptr_[MAXGRIDS]; diff --git a/src/external/libCuba/src/vegas/util.c b/src/external/libCuba/src/vegas/util.c deleted file mode 100644 index ae54a7ad..00000000 --- a/src/external/libCuba/src/vegas/util.c +++ /dev/null @@ -1,46 +0,0 @@ -/* - util.c - Utility functions - this file is part of Vegas - last modified 2 Mar 06 th -*/ - -/*************************************************************************** - * Copyright (C) 2004-2009 by Thomas Hahn * - * hahn@feynarts.de * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU Lesser General Public * - * License as published by the Free Software Foundation; either * - * version 2.1 of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * Lesser General Public License for more details. * - * * - * You should have received a copy of the GNU Lesser General Public * - * License along with this library; if not, write to the * - * Free Software Foundation, Inc., * - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA * - ***************************************************************************/ - -#include "decl.h" - -static count ndim_, ncomp_; -static number neval_; -static Grid *gridptr_[MAXGRIDS]; -static count griddim_[MAXGRIDS]; -int EXPORT(vegasnbatch) = 1000; -int EXPORT(vegasgridno) = 0; -char EXPORT(vegasstate)[MAXSTATESIZE] = ""; - - -#define SamplesAlloc(p, n) \ - MemAlloc(p, (n)*((ndim_ + ncomp_ + 1)*sizeof(real) + ndim_*sizeof(bin_t))) - - -#ifdef DEBUG -#include "debug.c" -#endif -