ddemidov / vexcl
1
#define BOOST_TEST_MODULE RandomNumbers
2
#include <boost/test/unit_test.hpp>
3
#include <vexcl/vector.hpp>
4
#include <vexcl/element_index.hpp>
5
#include <vexcl/random.hpp>
6
#include <vexcl/reductor.hpp>
7
#include <vexcl/tagged_terminal.hpp>
8
#include <vexcl/temporary.hpp>
9
#include <vexcl/function.hpp>
10
#include <boost/math/constants/constants.hpp>
11
#include "context_setup.hpp"
12

13 6
BOOST_AUTO_TEST_CASE(random_numbers)
14
{
15
    const size_t N = 1 << 20;
16

17 6
    vex::Reductor<size_t, vex::SUM> sumi(ctx);
18 6
    vex::Reductor<double, vex::SUM> sumd(ctx);
19

20 6
    vex::Random<cl_int> rand0;
21 6
    vex::vector<cl_uint> x0(ctx, N);
22 6
    x0 = rand0(vex::element_index(), std::rand());
23

24 6
    vex::Random<cl_float4> rand1;
25 6
    vex::vector<cl_float4> x1(ctx, N);
26 6
    x1 = rand1(vex::element_index(), std::rand());
27

28 6
    vex::Random<cl_double4> rand2;
29 6
    vex::vector<cl_double4> x2(ctx, N);
30 6
    x2 = rand2(vex::element_index(), std::rand());
31

32 6
    vex::Random<cl_double> rand3;
33 6
    vex::vector<cl_double> x3(ctx, N);
34 6
    x3 = rand3(vex::element_index(), std::rand());
35

36
    // X in [0,1]
37 6
    BOOST_CHECK(sumi(x3 > 1) == 0);
38 6
    BOOST_CHECK(sumi(x3 < 0) == 0);
39

40
    // mean = 0.5
41 6
    BOOST_CHECK(std::abs((sumd(x3) / N) - 0.5) < 1e-2);
42

43 6
    vex::RandomNormal<cl_double> rand4;
44 6
    vex::vector<cl_double> x4(ctx, N);
45 6
    x4 = rand4(vex::element_index(), std::rand());
46

47
    // E(X ~ N(0,s)) = 0
48 6
    BOOST_CHECK(std::abs(sumd(x4)/N) < 1e-2);
49

50
    // E(abs(X) ~ N(0,s)) = sqrt(2/M_PI) * s
51 6
    BOOST_CHECK(std::abs(sumd(fabs(x4))/N - std::sqrt(2 / boost::math::constants::pi<double>())) < 1e-2);
52

53 6
    vex::Random<cl_double, vex::random::threefry> rand5;
54 6
    vex::vector<cl_double> x5(ctx, N);
55 6
    x5 = rand5(vex::element_index(), std::rand());
56

57 6
    BOOST_CHECK(std::abs(sumd(x5)/N - 0.5) < 1e-2);
58
}
59

60 6
BOOST_AUTO_TEST_CASE(monte_carlo_pi)
61
{
62 6
    vex::Random<double, vex::random::threefry> rnd;
63

64 6
    vex::Reductor<size_t, vex::SUM> sum(ctx);
65

66
    const size_t n = 1 << 20;
67

68 6
    auto i = vex::tag<0>(vex::element_index(0, n));
69

70 6
    auto x = vex::make_temp<1>(rnd(i, std::rand()));
71 6
    auto y = vex::make_temp<2>(rnd(i, std::rand()));
72

73 6
    double pi = 4.0 * sum( (x * x + y * y) < 1 ) / n;
74

75 6
    BOOST_CHECK_CLOSE(pi, boost::math::constants::pi<double>(), 0.5);
76
}
77

78 6
BOOST_AUTO_TEST_SUITE_END()
79

Read our documentation on viewing source code .

Loading