From: Wolfgang Bangerth Date: Sun, 19 Apr 2015 21:56:10 +0000 (-0500) Subject: Replace use of rand_r by equivalent BOOST functionality. X-Git-Tag: v8.3.0-rc1~239^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=28a60bc984956d02dea34f2e6676d79ff205d014;p=dealii.git Replace use of rand_r by equivalent BOOST functionality. --- diff --git a/doc/news/changes.h b/doc/news/changes.h index ce4f8efc1a..12e0a202a9 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -352,6 +352,14 @@ inconvenience this causes.
    +
  1. Fixed: Utilities::generate_normal_random_number() will now + produce a deterministic sequence of numbers on every thread. + Furthermore, it will also work on systems that do not have the + rand_r library function, such as Cygwin. +
    + (Wolfgang Bangerth, 2015/04/19) +
  2. +
  3. New: A LinearOperator class that stores the abstract concept of a linear operator. The class is fully compatible with the solver and preconditioner interfaces. The primary purpose of this class is to diff --git a/include/deal.II/base/utilities.h b/include/deal.II/base/utilities.h index 4200224e60..739fb61d8a 100644 --- a/include/deal.II/base/utilities.h +++ b/include/deal.II/base/utilities.h @@ -194,15 +194,17 @@ namespace Utilities * Generate a random number from a normalized Gaussian probability * distribution centered around @p a and with standard deviation @p sigma. * - * This function is reentrant, i.e., it can safely be called from multiple - * threads at the same time. However, if so done, then there is no guarantee - * that each thread will get the same sequence of numbers every time. - * Rather, the produced sequence of random numbers will be apportioned to - * the different threads in non-deterministic ways. If this is a problem, - * for example for exactly reproducibility, then you need to use separate - * random number facilities for separate threads, rather than this global - * function. For example, the C++11 standard offers such objects, as does - * BOOST. + * This function is reentrant, i.e., it can safely be called from + * multiple threads at the same time. In addition, each thread will + * get the same sequence of numbers every time. On the other hand, + * if you run Threads::Task objects via the Threading Building + * Blocks, then tasks will be assigned to mostly random threads, and + * may get a different sequence of random numbers in different runs + * of the program, since a previous task may already have consumed + * the first few random numbers generated for the thread you're + * on. If this is a problem, you need to create your own random + * number generator objects every time you want to start from a + * defined point. * * @note Like the system function rand(), this function produces the same * sequence of random numbers every time a program is started. This is an diff --git a/source/base/utilities.cc b/source/base/utilities.cc index e3d872b5b6..d776a0ab62 100644 --- a/source/base/utilities.cc +++ b/source/base/utilities.cc @@ -16,16 +16,17 @@ #include #include +#include #include #include +#include #include #include #include #include #include -#include #include #include #include @@ -390,41 +391,13 @@ namespace Utilities if (sigma == 0) return a; - // we want to use rand(), but that function is not reentrant in a thread - // context. thus, use rand_r. this does not produce reproducible results - // between threads either, but at least it is reentrant. if you need - // an exactly reproducible sequence even in multithreaded contexts, - // then this is probably not the function to use. -#ifdef DEAL_II_HAVE_RAND_R - static unsigned int seed = 0xabcd1234; - const double y = 1.0*rand_r(&seed)/RAND_MAX; -#else - const double y = 1.0*rand()/RAND_MAX; -#endif - - // find x such that y=erf(x). do so - // using a Newton method to find - // the zero of F(x)=erf(x)-y. start - // at x=0 - double x = 0; - unsigned int iteration = 0; - while (true) - { - const double residual = 0.5+boost::math::erf(x/std::sqrt(2.)/sigma)/2-y; - - if (std::fabs(residual) < 1e-7) - break; - - const double F_prime = 1./std::sqrt(2*3.1415926536)/sigma * - std::exp(-x*x/sigma/sigma/2); - x += -residual / F_prime; - - // make sure that we don't - // recurse endlessly - ++iteration; - Assert (iteration < 20, ExcInternalError()); - }; - return x+a; + // we would want to use rand(), but that function is not reentrant + // in a thread context. one could use rand_r, but this does not + // produce reproducible results between threads either (though at + // least it is reentrant). these two approaches being + // non-workable, use a thread-local random number generator here + static Threads::ThreadLocalStorage random_number_generator; + return boost::normal_distribution<>(a,sigma)(random_number_generator.get()); } diff --git a/tests/base/utilities_01.output b/tests/base/utilities_01.output index f8f85d1c68..aa659424e0 100644 --- a/tests/base/utilities_01.output +++ b/tests/base/utilities_01.output @@ -6,4 +6,4 @@ DEAL::413 DEAL::1 DEAL::-12 DEAL::125 -DEAL::-15.3253 -67.9624 74.3805 +DEAL::6.39848 -6.55406 -60.8373