From 6d7b889d76da70d973e6b559db228e61d6b41ca5 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Wed, 5 Feb 2025 10:10:02 -0800 Subject: [PATCH 1/6] fix recursion-related nvlink warning in RandomGamma --- Src/Base/AMReX_Random.H | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/Src/Base/AMReX_Random.H b/Src/Base/AMReX_Random.H index a8bda8b036d..f378cb77bf8 100644 --- a/Src/Base/AMReX_Random.H +++ b/Src/Base/AMReX_Random.H @@ -132,6 +132,7 @@ namespace amrex */ Real RandomGamma (Real alpha, Real beta); + template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real RandomGamma (Real alpha, Real beta, RandomEngine const& random_engine) { @@ -141,8 +142,12 @@ namespace amrex AMREX_IF_ON_DEVICE(( if (alpha < 1) { - Real u = amrex::Random(random_engine); - return RandomGamma(1.0_rt + alpha, beta, random_engine) * std::pow(u, 1.0_rt / alpha); + if constexpr (depth == 0) + { + // note that alpha is assumed to be > 0, so this function will recurse at most once. + Real u = amrex::Random(random_engine); + return RandomGamma<1>(1.0_rt + alpha, beta, random_engine) * std::pow(u, 1.0_rt / alpha); + } } { From 19e0de431ba52516cd3ec7b8f0e3fdc121382af6 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Thu, 6 Feb 2025 11:33:54 -0800 Subject: [PATCH 2/6] remove recursion altogether --- Src/Base/AMReX_Random.H | 74 ++++++++++++++++++++++------------------- 1 file changed, 39 insertions(+), 35 deletions(-) diff --git a/Src/Base/AMReX_Random.H b/Src/Base/AMReX_Random.H index f378cb77bf8..70e6a65c8c1 100644 --- a/Src/Base/AMReX_Random.H +++ b/Src/Base/AMReX_Random.H @@ -119,6 +119,39 @@ namespace amrex #endif } + namespace { + + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + Real RandomGamma_alpha_ge_1 (Real alpha, Real beta, RandomEngine const& random_engine) + { + AMREX_ASSERT(alpha >= 1); + AMREX_ASSERT(beta > 0); + + Real x, v, u; + Real d = alpha - 1.0_rt / 3.0_rt; + Real c = (1.0_rt / 3.0_rt) / std::sqrt(d); + + while (1) { + do { + x = amrex::RandomNormal(0.0_rt, 1.0_rt, random_engine); + v = 1.0_rt + c * x; + } while (v <= 0.0_rt); + + v = v * v * v; + u = amrex::Random(random_engine); + + if (u < 1.0_rt - 0.0331_rt * x * x * x * x) { + break; + } + + if (std::log(u) < 0.5_rt * x * x + d * (1.0_rt - v + std::log(v))) { + break; + } + } + return beta * d * v; + } + } + /** * \brief Generate a psuedo-random floating point number from the Gamma distribution * @@ -132,7 +165,6 @@ namespace amrex */ Real RandomGamma (Real alpha, Real beta); - template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real RandomGamma (Real alpha, Real beta, RandomEngine const& random_engine) { @@ -142,45 +174,17 @@ namespace amrex AMREX_IF_ON_DEVICE(( if (alpha < 1) { - if constexpr (depth == 0) - { - // note that alpha is assumed to be > 0, so this function will recurse at most once. - Real u = amrex::Random(random_engine); - return RandomGamma<1>(1.0_rt + alpha, beta, random_engine) * std::pow(u, 1.0_rt / alpha); - } - } - - { - Real x, v, u; - Real d = alpha - 1.0_rt / 3.0_rt; - Real c = (1.0_rt / 3.0_rt) / std::sqrt(d); - - while (1) { - do { - x = amrex::RandomNormal(0.0_rt, 1.0_rt, random_engine); - v = 1.0_rt + c * x; - } while (v <= 0.0_rt); - - v = v * v * v; - u = amrex::Random(random_engine); - - if (u < 1.0_rt - 0.0331_rt * x * x * x * x) { - break; - } - - if (std::log(u) < 0.5_rt * x * x + d * (1.0_rt - v + std::log(v))) { - break; - } - } - return beta * d * v; + Real u = amrex::Random(random_engine); + return RandomGamma_alpha_ge_1(1.0_rt + alpha, beta, random_engine) * std::pow(u, 1.0_rt / alpha); + } else { + RandomGamma_alpha_ge_1(alpha, beta, random_engine); } )) AMREX_IF_ON_HOST(( - amrex::ignore_unused(random_engine); - return RandomGamma(alpha, beta); + amrex::ignore_unused(random_engine); + return RandomGamma(alpha, beta); )) - } /** From 895cda56f7ba242f9c64353828e32587ca465327 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Thu, 6 Feb 2025 11:42:33 -0800 Subject: [PATCH 3/6] fix missing return statement --- Src/Base/AMReX_Random.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Src/Base/AMReX_Random.H b/Src/Base/AMReX_Random.H index 70e6a65c8c1..2d98e1ca3cb 100644 --- a/Src/Base/AMReX_Random.H +++ b/Src/Base/AMReX_Random.H @@ -177,7 +177,7 @@ namespace amrex Real u = amrex::Random(random_engine); return RandomGamma_alpha_ge_1(1.0_rt + alpha, beta, random_engine) * std::pow(u, 1.0_rt / alpha); } else { - RandomGamma_alpha_ge_1(alpha, beta, random_engine); + return RandomGamma_alpha_ge_1(alpha, beta, random_engine); } )) @@ -185,7 +185,7 @@ namespace amrex amrex::ignore_unused(random_engine); return RandomGamma(alpha, beta); )) - } + /** * \brief Generates one pseudorandom unsigned integer which is From 33fd30efa6c819ad62d76987cb6b3cce9a5df31f Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Thu, 6 Feb 2025 11:46:46 -0800 Subject: [PATCH 4/6] fix syntax error --- Src/Base/AMReX_Random.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Src/Base/AMReX_Random.H b/Src/Base/AMReX_Random.H index 2d98e1ca3cb..396675500de 100644 --- a/Src/Base/AMReX_Random.H +++ b/Src/Base/AMReX_Random.H @@ -185,7 +185,7 @@ namespace amrex amrex::ignore_unused(random_engine); return RandomGamma(alpha, beta); )) - + } /** * \brief Generates one pseudorandom unsigned integer which is From eba0105f22e2eb268f68c53fc311880000dd58e5 Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Thu, 6 Feb 2025 12:13:00 -0800 Subject: [PATCH 5/6] remove anon namespace --- Src/Base/AMReX_Random.H | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Src/Base/AMReX_Random.H b/Src/Base/AMReX_Random.H index 396675500de..ffa03a33af4 100644 --- a/Src/Base/AMReX_Random.H +++ b/Src/Base/AMReX_Random.H @@ -119,7 +119,7 @@ namespace amrex #endif } - namespace { + namespace random_util { AMREX_GPU_DEVICE AMREX_FORCE_INLINE Real RandomGamma_alpha_ge_1 (Real alpha, Real beta, RandomEngine const& random_engine) @@ -175,9 +175,9 @@ namespace amrex if (alpha < 1) { Real u = amrex::Random(random_engine); - return RandomGamma_alpha_ge_1(1.0_rt + alpha, beta, random_engine) * std::pow(u, 1.0_rt / alpha); + return amrex::random_util::RandomGamma_alpha_ge_1(1.0_rt + alpha, beta, random_engine) * std::pow(u, 1.0_rt / alpha); } else { - return RandomGamma_alpha_ge_1(alpha, beta, random_engine); + return amrex::random_util::RandomGamma_alpha_ge_1(alpha, beta, random_engine); } )) From eecb62764f27d2c4ca9878e3d12adbe57c2a3d1b Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Thu, 6 Feb 2025 12:44:39 -0800 Subject: [PATCH 6/6] while (true) to pacify clang --- Src/Base/AMReX_Random.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Src/Base/AMReX_Random.H b/Src/Base/AMReX_Random.H index ffa03a33af4..09b2c18c696 100644 --- a/Src/Base/AMReX_Random.H +++ b/Src/Base/AMReX_Random.H @@ -131,7 +131,7 @@ namespace amrex Real d = alpha - 1.0_rt / 3.0_rt; Real c = (1.0_rt / 3.0_rt) / std::sqrt(d); - while (1) { + while (true) { do { x = amrex::RandomNormal(0.0_rt, 1.0_rt, random_engine); v = 1.0_rt + c * x;