monte_carlo.cpp 2.4 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273
  1. //---------------------------------------------------------------------------//
  2. // Copyright (c) 2013 Kyle Lutz <kyle.r.lutz@gmail.com>
  3. //
  4. // Distributed under the Boost Software License, Version 1.0
  5. // See accompanying file LICENSE_1_0.txt or copy at
  6. // http://www.boost.org/LICENSE_1_0.txt
  7. //
  8. // See http://boostorg.github.com/compute for more information.
  9. //---------------------------------------------------------------------------//
  10. #include <iostream>
  11. #include <boost/compute/function.hpp>
  12. #include <boost/compute/system.hpp>
  13. #include <boost/compute/algorithm/count_if.hpp>
  14. #include <boost/compute/container/vector.hpp>
  15. #include <boost/compute/iterator/buffer_iterator.hpp>
  16. #include <boost/compute/random/default_random_engine.hpp>
  17. #include <boost/compute/types/fundamental.hpp>
  18. namespace compute = boost::compute;
  19. int main()
  20. {
  21. // get default device and setup context
  22. compute::device gpu = compute::system::default_device();
  23. compute::context context(gpu);
  24. compute::command_queue queue(context, gpu);
  25. std::cout << "device: " << gpu.name() << std::endl;
  26. using compute::uint_;
  27. using compute::uint2_;
  28. #ifdef CI_BUILD // lower number of points for CI builds
  29. size_t n = 2000;
  30. #else
  31. // ten million random points
  32. size_t n = 10000000;
  33. #endif
  34. // generate random numbers
  35. compute::default_random_engine rng(queue);
  36. compute::vector<uint_> vector(n * 2, context);
  37. rng.generate(vector.begin(), vector.end(), queue);
  38. // function returing true if the point is within the unit circle
  39. BOOST_COMPUTE_FUNCTION(bool, is_in_unit_circle, (const uint2_ point),
  40. {
  41. const float x = point.x / (float) UINT_MAX - 1;
  42. const float y = point.y / (float) UINT_MAX - 1;
  43. return (x*x + y*y) < 1.0f;
  44. });
  45. // iterate over vector<uint> as vector<uint2>
  46. compute::buffer_iterator<uint2_> start =
  47. compute::make_buffer_iterator<uint2_>(vector.get_buffer(), 0);
  48. compute::buffer_iterator<uint2_> end =
  49. compute::make_buffer_iterator<uint2_>(vector.get_buffer(), vector.size() / 2);
  50. // count number of random points within the unit circle
  51. size_t count = compute::count_if(start, end, is_in_unit_circle, queue);
  52. // print out values
  53. float count_f = static_cast<float>(count);
  54. std::cout << "count: " << count << " / " << n << std::endl;
  55. std::cout << "ratio: " << count_f / float(n) << std::endl;
  56. std::cout << "pi = " << (count_f / float(n)) * 4.0f << std::endl;
  57. return 0;
  58. }