<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
* 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
#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>
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());
}
DEAL::1
DEAL::-12
DEAL::125
-DEAL::-15.3253 -67.9624 74.3805
+DEAL::6.39848 -6.55406 -60.8373