Replace Marsaglia polar method with Box-muller to generate a normally distributed random number

[Shihab-Shahriar]

On GPU, Box-Muller is expected to significantly outperform current Marsaglia polar. The reasoning can be found here:

The polar method (Press et. al. 1992) is simple and relatively efficient, but the probability of looping per thread is 14 percent. This leads to an expected 1.6 iterations per generated sample turning into an expected 3.1 iterations when warp effects are taken into account.

On the CPU side, there doesn't seem to be any difference.

Before

N	CPU	CPU	CUDA	CUDA
	64	1024	64	1024
1e4	0.061071	0.0757334	0.0150474	0.016794
1e5	0.600972	0.761684	0.120827	0.28777
2e5	1.19428	1.50356	0.236278	0.586949

After

N	CPU	CPU	CUDA	CUDA
	64	1024	64	1024
1e4	0.0610803	0.0734002	0.0119627	0.0119717
1e5	0.60062	0.73712	0.0957707	0.0968437
2e5	1.17503	1.47986	0.195854	0.19791

(The timings are taken using one of the kokkos examples here- generate a 2D array with N rows and 1000 columns. All values in seconds. CPU is AMD 3700x, GPU is Nvidia 1650 super, 16GB RAM)

I'm not quite sure why there is such a performance gap in Polar method on GPU between the 64 and 1024 bit versions. Maybe this can be tuned out by trying different kernel launch configurations (i.e. #blocks, #threads/block). I wanted to present the results before diving deeper.

[dalg24-jenkins]

Can one of the admins verify this patch?

[fnrizzi]

ok to test

[fnrizzi]

can you please also include the code you used to generate all the timings?

[fnrizzi]

format is wrong , you need to run the /scripts/apply-clang-format before pushing the repo

[dalg24]

format is wrong , you need to run the /scripts/apply-clang-format before pushing the repo

diff --git a/algorithms/src/Kokkos_Random.hpp b/algorithms/src/Kokkos_Random.hpp
index 0057c79cb..0ee96098c 100644
--- a/algorithms/src/Kokkos_Random.hpp
+++ b/algorithms/src/Kokkos_Random.hpp
@@ -855,9 +855,9 @@ class Random_XorShift64 {
   double normal() {
     constexpr double M_PI2 = 2.0 * Kokkos::numbers::pi_v<double>;
 
-    double u = drand();
-    double v = drand();
-    double r = Kokkos::sqrt(-2.0 * Kokkos::log(u));
+    double u     = drand();
+    double v     = drand();
+    double r     = Kokkos::sqrt(-2.0 * Kokkos::log(u));
     double theta = v * M_PI2;
     return r * Kokkos::cos(theta);
   }
@@ -1099,9 +1099,9 @@ class Random_XorShift1024 {
   double normal() {
     constexpr double M_PI2 = 2.0 * Kokkos::numbers::pi_v<double>;
 
-    double u = drand();
-    double v = drand();
-    double r = Kokkos::sqrt(-2.0 * Kokkos::log(u));
+    double u     = drand();
+    double v     = drand();
+    double r     = Kokkos::sqrt(-2.0 * Kokkos::log(u));
     double theta = v * M_PI2;
     return r * Kokkos::cos(theta);
   }

[Shihab-Shahriar]

@fnrizzi , here is the code I used to benchmark. Please note that the random example isn't automatically built by kokkos even with Kokkos_ENABLE_EXAMPLES=ON, so I had to edit some CMakeLists.txt files. The CMake options I used are in the build-job.sh file at the root.

I just pushed the formatting-related changes.

[fnrizzi]

thank you for working on this!

[Shihab-Shahriar]

Thanks for your feedback and help.