]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Replace use of rand_r by equivalent BOOST functionality.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 19 Apr 2015 21:56:10 +0000 (16:56 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 19 Apr 2015 22:00:25 +0000 (17:00 -0500)
doc/news/changes.h
include/deal.II/base/utilities.h
source/base/utilities.cc
tests/base/utilities_01.output

index ce4f8efc1a929f894ca5554f02a6fbdb39f58b6c..12e0a202a9e719964790990095545b1850ebf08d 100644 (file)
@@ -352,6 +352,14 @@ inconvenience this causes.
 
 
 <ol>
+  <li> 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
+  <code>rand_r</code> library function, such as Cygwin.
+  <br>
+  (Wolfgang Bangerth, 2015/04/19)
+  </li>
+  
   <li> 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
index 4200224e60647ef3b95f4cf574538e23d81da251..739fb61d8a79cf74eeaaec62b82e780d133e09c8 100644 (file)
@@ -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
index e3d872b5b6e8e04fe0c9bb1e2c76b90ba64843b1..d776a0ab623457327bf40ba93ba2d4f4f228053d 100644 (file)
 
 #include <deal.II/base/utilities.h>
 #include <deal.II/base/exceptions.h>
+#include <deal.II/base/thread_local_storage.h>
 
 #include <boost/math/special_functions/erf.hpp>
 #include <boost/lexical_cast.hpp>
+#include <boost/random.hpp>
 
 #include <algorithm>
 #include <cerrno>
 #include <cmath>
 #include <cstddef>
 #include <cstdio>
-#include <cstdlib>
 #include <ctime>
 #include <fstream>
 #include <iomanip>
@@ -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<boost::mt19937> random_number_generator;
+    return boost::normal_distribution<>(a,sigma)(random_number_generator.get());
   }
 
 
index f8f85d1c683cc4c27f15e5fe8b998d518803eb6b..aa659424e0a2e0d93f3bc8fe353d7a6a20509c03 100644 (file)
@@ -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 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.