diff --git a/include/boost/math/distributions/chi_squared.hpp b/include/boost/math/distributions/chi_squared.hpp index f5daddc0ad..3944569e89 100644 --- a/include/boost/math/distributions/chi_squared.hpp +++ b/include/boost/math/distributions/chi_squared.hpp @@ -9,14 +9,17 @@ #ifndef BOOST_MATH_DISTRIBUTIONS_CHI_SQUARED_HPP #define BOOST_MATH_DISTRIBUTIONS_CHI_SQUARED_HPP +#include +#include +#include +#include +#include #include #include // for incomplete beta. #include // complements #include // error checks #include -#include - namespace boost{ namespace math{ template > @@ -26,20 +29,20 @@ class chi_squared_distribution using value_type = RealType; using policy_type = Policy; - explicit chi_squared_distribution(RealType i) : m_df(i) + BOOST_MATH_GPU_ENABLED explicit chi_squared_distribution(RealType i) : m_df(i) { RealType result; detail::check_df( "boost::math::chi_squared_distribution<%1%>::chi_squared_distribution", m_df, &result, Policy()); } // chi_squared_distribution - RealType degrees_of_freedom()const + BOOST_MATH_GPU_ENABLED RealType degrees_of_freedom()const { return m_df; } // Parameter estimation: - static RealType find_degrees_of_freedom( + BOOST_MATH_GPU_ENABLED static RealType find_degrees_of_freedom( RealType difference_from_variance, RealType alpha, RealType beta, @@ -66,16 +69,16 @@ chi_squared_distribution(RealType)->chi_squared_distribution -inline std::pair range(const chi_squared_distribution& /*dist*/) +BOOST_MATH_GPU_ENABLED inline boost::math::pair range(const chi_squared_distribution& /*dist*/) { // Range of permissible values for random variable x. - if (std::numeric_limits::has_infinity) + BOOST_MATH_IF_CONSTEXPR (boost::math::numeric_limits::has_infinity) { - return std::pair(static_cast(0), std::numeric_limits::infinity()); // 0 to + infinity. + return boost::math::pair(static_cast(0), boost::math::numeric_limits::infinity()); // 0 to + infinity. } else { using boost::math::tools::max_value; - return std::pair(static_cast(0), max_value()); // 0 to + max. + return boost::math::pair(static_cast(0), max_value()); // 0 to + max. } } @@ -84,21 +87,21 @@ inline std::pair range(const chi_squared_distribution -inline std::pair support(const chi_squared_distribution& /*dist*/) +BOOST_MATH_GPU_ENABLED inline boost::math::pair support(const chi_squared_distribution& /*dist*/) { // Range of supported values for random variable x. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. - return std::pair(static_cast(0), tools::max_value()); // 0 to + infinity. + return boost::math::pair(static_cast(0), tools::max_value()); // 0 to + infinity. } template -RealType pdf(const chi_squared_distribution& dist, const RealType& chi_square) +BOOST_MATH_GPU_ENABLED RealType pdf(const chi_squared_distribution& dist, const RealType& chi_square) { BOOST_MATH_STD_USING // for ADL of std functions RealType degrees_of_freedom = dist.degrees_of_freedom(); // Error check: RealType error_result; - static const char* function = "boost::math::pdf(const chi_squared_distribution<%1%>&, %1%)"; + constexpr auto function = "boost::math::pdf(const chi_squared_distribution<%1%>&, %1%)"; if(false == detail::check_df( function, degrees_of_freedom, &error_result, Policy())) @@ -132,12 +135,12 @@ RealType pdf(const chi_squared_distribution& dist, const RealT } // pdf template -inline RealType cdf(const chi_squared_distribution& dist, const RealType& chi_square) +BOOST_MATH_GPU_ENABLED inline RealType cdf(const chi_squared_distribution& dist, const RealType& chi_square) { RealType degrees_of_freedom = dist.degrees_of_freedom(); // Error check: RealType error_result; - static const char* function = "boost::math::cdf(const chi_squared_distribution<%1%>&, %1%)"; + constexpr auto function = "boost::math::cdf(const chi_squared_distribution<%1%>&, %1%)"; if(false == detail::check_df( function, degrees_of_freedom, &error_result, Policy())) @@ -153,10 +156,10 @@ inline RealType cdf(const chi_squared_distribution& dist, cons } // cdf template -inline RealType quantile(const chi_squared_distribution& dist, const RealType& p) +BOOST_MATH_GPU_ENABLED inline RealType quantile(const chi_squared_distribution& dist, const RealType& p) { RealType degrees_of_freedom = dist.degrees_of_freedom(); - static const char* function = "boost::math::quantile(const chi_squared_distribution<%1%>&, %1%)"; + constexpr auto function = "boost::math::quantile(const chi_squared_distribution<%1%>&, %1%)"; // Error check: RealType error_result; if(false == @@ -170,11 +173,11 @@ inline RealType quantile(const chi_squared_distribution& dist, } // quantile template -inline RealType cdf(const complemented2_type, RealType>& c) +BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type, RealType>& c) { RealType const& degrees_of_freedom = c.dist.degrees_of_freedom(); RealType const& chi_square = c.param; - static const char* function = "boost::math::cdf(const chi_squared_distribution<%1%>&, %1%)"; + constexpr auto function = "boost::math::cdf(const chi_squared_distribution<%1%>&, %1%)"; // Error check: RealType error_result; if(false == detail::check_df( @@ -191,11 +194,11 @@ inline RealType cdf(const complemented2_type -inline RealType quantile(const complemented2_type, RealType>& c) +BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type, RealType>& c) { RealType const& degrees_of_freedom = c.dist.degrees_of_freedom(); RealType const& q = c.param; - static const char* function = "boost::math::quantile(const chi_squared_distribution<%1%>&, %1%)"; + constexpr auto function = "boost::math::quantile(const chi_squared_distribution<%1%>&, %1%)"; // Error check: RealType error_result; if(false == ( @@ -208,22 +211,22 @@ inline RealType quantile(const complemented2_type -inline RealType mean(const chi_squared_distribution& dist) +BOOST_MATH_GPU_ENABLED inline RealType mean(const chi_squared_distribution& dist) { // Mean of Chi-Squared distribution = v. return dist.degrees_of_freedom(); } // mean template -inline RealType variance(const chi_squared_distribution& dist) +BOOST_MATH_GPU_ENABLED inline RealType variance(const chi_squared_distribution& dist) { // Variance of Chi-Squared distribution = 2v. return 2 * dist.degrees_of_freedom(); } // variance template -inline RealType mode(const chi_squared_distribution& dist) +BOOST_MATH_GPU_ENABLED inline RealType mode(const chi_squared_distribution& dist) { RealType df = dist.degrees_of_freedom(); - static const char* function = "boost::math::mode(const chi_squared_distribution<%1%>&)"; + constexpr auto function = "boost::math::mode(const chi_squared_distribution<%1%>&)"; if(df < 2) return policies::raise_domain_error( @@ -234,7 +237,7 @@ inline RealType mode(const chi_squared_distribution& dist) } template -inline RealType skewness(const chi_squared_distribution& dist) +BOOST_MATH_GPU_ENABLED inline RealType skewness(const chi_squared_distribution& dist) { BOOST_MATH_STD_USING // For ADL RealType df = dist.degrees_of_freedom(); @@ -242,14 +245,14 @@ inline RealType skewness(const chi_squared_distribution& dist) } template -inline RealType kurtosis(const chi_squared_distribution& dist) +BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const chi_squared_distribution& dist) { RealType df = dist.degrees_of_freedom(); return 3 + 12 / df; } template -inline RealType kurtosis_excess(const chi_squared_distribution& dist) +BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const chi_squared_distribution& dist) { RealType df = dist.degrees_of_freedom(); return 12 / df; @@ -264,12 +267,12 @@ namespace detail template struct df_estimator { - df_estimator(RealType a, RealType b, RealType variance, RealType delta) + BOOST_MATH_GPU_ENABLED df_estimator(RealType a, RealType b, RealType variance, RealType delta) : alpha(a), beta(b), ratio(delta/variance) { // Constructor } - RealType operator()(const RealType& df) + BOOST_MATH_GPU_ENABLED RealType operator()(const RealType& df) { if(df <= tools::min_value()) return 1; @@ -297,14 +300,14 @@ struct df_estimator } // namespace detail template -RealType chi_squared_distribution::find_degrees_of_freedom( +BOOST_MATH_GPU_ENABLED RealType chi_squared_distribution::find_degrees_of_freedom( RealType difference_from_variance, RealType alpha, RealType beta, RealType variance, RealType hint) { - static const char* function = "boost::math::chi_squared_distribution<%1%>::find_degrees_of_freedom(%1%,%1%,%1%,%1%,%1%)"; + constexpr auto function = "boost::math::chi_squared_distribution<%1%>::find_degrees_of_freedom(%1%,%1%,%1%,%1%,%1%)"; // Check for domain errors: RealType error_result; if(false == @@ -321,8 +324,8 @@ RealType chi_squared_distribution::find_degrees_of_freedom( detail::df_estimator f(alpha, beta, variance, difference_from_variance); tools::eps_tolerance tol(policies::digits()); - std::uintmax_t max_iter = policies::get_max_root_iterations(); - std::pair r = + boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); + boost::math::pair r = tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy()); RealType result = r.first + (r.second - r.first) / 2; if(max_iter >= policies::get_max_root_iterations()) diff --git a/include/boost/math/distributions/weibull.hpp b/include/boost/math/distributions/weibull.hpp index ca4bbd7b53..eb4de106c8 100644 --- a/include/boost/math/distributions/weibull.hpp +++ b/include/boost/math/distributions/weibull.hpp @@ -1,4 +1,5 @@ // Copyright John Maddock 2006. +// Copyright Matt Borland 2024. // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) @@ -9,6 +10,10 @@ // http://www.itl.nist.gov/div898/handbook/eda/section3/eda3668.htm // http://mathworld.wolfram.com/WeibullDistribution.html +#include +#include +#include +#include #include #include #include @@ -16,14 +21,12 @@ #include #include -#include - namespace boost{ namespace math { namespace detail{ template -inline bool check_weibull_shape( +BOOST_MATH_GPU_ENABLED inline bool check_weibull_shape( const char* function, RealType shape, RealType* result, const Policy& pol) @@ -39,7 +42,7 @@ inline bool check_weibull_shape( } template -inline bool check_weibull_x( +BOOST_MATH_GPU_ENABLED inline bool check_weibull_x( const char* function, RealType const& x, RealType* result, const Policy& pol) @@ -55,7 +58,7 @@ inline bool check_weibull_x( } template -inline bool check_weibull( +BOOST_MATH_GPU_ENABLED inline bool check_weibull( const char* function, RealType scale, RealType shape, @@ -73,19 +76,19 @@ class weibull_distribution using value_type = RealType; using policy_type = Policy; - explicit weibull_distribution(RealType l_shape, RealType l_scale = 1) + BOOST_MATH_GPU_ENABLED explicit weibull_distribution(RealType l_shape, RealType l_scale = 1) : m_shape(l_shape), m_scale(l_scale) { RealType result; detail::check_weibull("boost::math::weibull_distribution<%1%>::weibull_distribution", l_scale, l_shape, &result, Policy()); } - RealType shape()const + BOOST_MATH_GPU_ENABLED RealType shape()const { return m_shape; } - RealType scale()const + BOOST_MATH_GPU_ENABLED RealType scale()const { return m_scale; } @@ -107,28 +110,28 @@ weibull_distribution(RealType,RealType)->weibull_distribution -inline std::pair range(const weibull_distribution& /*dist*/) +BOOST_MATH_GPU_ENABLED inline boost::math::pair range(const weibull_distribution& /*dist*/) { // Range of permissible values for random variable x. using boost::math::tools::max_value; - return std::pair(static_cast(0), max_value()); + return boost::math::pair(static_cast(0), max_value()); } template -inline std::pair support(const weibull_distribution& /*dist*/) +BOOST_MATH_GPU_ENABLED inline boost::math::pair support(const weibull_distribution& /*dist*/) { // Range of supported values for random variable x. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. using boost::math::tools::max_value; using boost::math::tools::min_value; - return std::pair(min_value(), max_value()); + return boost::math::pair(min_value(), max_value()); // A discontinuity at x == 0, so only support down to min_value. } template -inline RealType pdf(const weibull_distribution& dist, const RealType& x) +BOOST_MATH_GPU_ENABLED inline RealType pdf(const weibull_distribution& dist, const RealType& x) { BOOST_MATH_STD_USING // for ADL of std functions - static const char* function = "boost::math::pdf(const weibull_distribution<%1%>, %1%)"; + constexpr auto function = "boost::math::pdf(const weibull_distribution<%1%>, %1%)"; RealType shape = dist.shape(); RealType scale = dist.scale(); @@ -158,11 +161,11 @@ inline RealType pdf(const weibull_distribution& dist, const Re } template -inline RealType logpdf(const weibull_distribution& dist, const RealType& x) +BOOST_MATH_GPU_ENABLED inline RealType logpdf(const weibull_distribution& dist, const RealType& x) { BOOST_MATH_STD_USING // for ADL of std functions - static const char* function = "boost::math::logpdf(const weibull_distribution<%1%>, %1%)"; + constexpr auto function = "boost::math::logpdf(const weibull_distribution<%1%>, %1%)"; RealType shape = dist.shape(); RealType scale = dist.scale(); @@ -192,11 +195,11 @@ inline RealType logpdf(const weibull_distribution& dist, const } template -inline RealType cdf(const weibull_distribution& dist, const RealType& x) +BOOST_MATH_GPU_ENABLED inline RealType cdf(const weibull_distribution& dist, const RealType& x) { BOOST_MATH_STD_USING // for ADL of std functions - static const char* function = "boost::math::cdf(const weibull_distribution<%1%>, %1%)"; + constexpr auto function = "boost::math::cdf(const weibull_distribution<%1%>, %1%)"; RealType shape = dist.shape(); RealType scale = dist.scale(); @@ -213,11 +216,11 @@ inline RealType cdf(const weibull_distribution& dist, const Re } template -inline RealType logcdf(const weibull_distribution& dist, const RealType& x) +BOOST_MATH_GPU_ENABLED inline RealType logcdf(const weibull_distribution& dist, const RealType& x) { BOOST_MATH_STD_USING // for ADL of std functions - static const char* function = "boost::math::logcdf(const weibull_distribution<%1%>, %1%)"; + constexpr auto function = "boost::math::logcdf(const weibull_distribution<%1%>, %1%)"; RealType shape = dist.shape(); RealType scale = dist.scale(); @@ -234,11 +237,11 @@ inline RealType logcdf(const weibull_distribution& dist, const } template -inline RealType quantile(const weibull_distribution& dist, const RealType& p) +BOOST_MATH_GPU_ENABLED inline RealType quantile(const weibull_distribution& dist, const RealType& p) { BOOST_MATH_STD_USING // for ADL of std functions - static const char* function = "boost::math::quantile(const weibull_distribution<%1%>, %1%)"; + constexpr auto function = "boost::math::quantile(const weibull_distribution<%1%>, %1%)"; RealType shape = dist.shape(); RealType scale = dist.scale(); @@ -258,11 +261,11 @@ inline RealType quantile(const weibull_distribution& dist, con } template -inline RealType cdf(const complemented2_type, RealType>& c) +BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type, RealType>& c) { BOOST_MATH_STD_USING // for ADL of std functions - static const char* function = "boost::math::cdf(const weibull_distribution<%1%>, %1%)"; + constexpr auto function = "boost::math::cdf(const weibull_distribution<%1%>, %1%)"; RealType shape = c.dist.shape(); RealType scale = c.dist.scale(); @@ -279,11 +282,11 @@ inline RealType cdf(const complemented2_type -inline RealType logcdf(const complemented2_type, RealType>& c) +BOOST_MATH_GPU_ENABLED inline RealType logcdf(const complemented2_type, RealType>& c) { BOOST_MATH_STD_USING // for ADL of std functions - static const char* function = "boost::math::logcdf(const weibull_distribution<%1%>, %1%)"; + constexpr auto function = "boost::math::logcdf(const weibull_distribution<%1%>, %1%)"; RealType shape = c.dist.shape(); RealType scale = c.dist.scale(); @@ -300,11 +303,11 @@ inline RealType logcdf(const complemented2_type -inline RealType quantile(const complemented2_type, RealType>& c) +BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type, RealType>& c) { BOOST_MATH_STD_USING // for ADL of std functions - static const char* function = "boost::math::quantile(const weibull_distribution<%1%>, %1%)"; + constexpr auto function = "boost::math::quantile(const weibull_distribution<%1%>, %1%)"; RealType shape = c.dist.shape(); RealType scale = c.dist.scale(); @@ -325,11 +328,11 @@ inline RealType quantile(const complemented2_type -inline RealType mean(const weibull_distribution& dist) +BOOST_MATH_GPU_ENABLED inline RealType mean(const weibull_distribution& dist) { BOOST_MATH_STD_USING // for ADL of std functions - static const char* function = "boost::math::mean(const weibull_distribution<%1%>)"; + constexpr auto function = "boost::math::mean(const weibull_distribution<%1%>)"; RealType shape = dist.shape(); RealType scale = dist.scale(); @@ -343,12 +346,12 @@ inline RealType mean(const weibull_distribution& dist) } template -inline RealType variance(const weibull_distribution& dist) +BOOST_MATH_GPU_ENABLED inline RealType variance(const weibull_distribution& dist) { RealType shape = dist.shape(); RealType scale = dist.scale(); - static const char* function = "boost::math::variance(const weibull_distribution<%1%>)"; + constexpr auto function = "boost::math::variance(const weibull_distribution<%1%>)"; RealType result = 0; if(false == detail::check_weibull(function, scale, shape, &result, Policy())) @@ -363,11 +366,11 @@ inline RealType variance(const weibull_distribution& dist) } template -inline RealType mode(const weibull_distribution& dist) +BOOST_MATH_GPU_ENABLED inline RealType mode(const weibull_distribution& dist) { BOOST_MATH_STD_USING // for ADL of std function pow. - static const char* function = "boost::math::mode(const weibull_distribution<%1%>)"; + constexpr auto function = "boost::math::mode(const weibull_distribution<%1%>)"; RealType shape = dist.shape(); RealType scale = dist.scale(); @@ -384,11 +387,11 @@ inline RealType mode(const weibull_distribution& dist) } template -inline RealType median(const weibull_distribution& dist) +BOOST_MATH_GPU_ENABLED inline RealType median(const weibull_distribution& dist) { BOOST_MATH_STD_USING // for ADL of std function pow. - static const char* function = "boost::math::median(const weibull_distribution<%1%>)"; + constexpr auto function = "boost::math::median(const weibull_distribution<%1%>)"; RealType shape = dist.shape(); // Wikipedia k RealType scale = dist.scale(); // Wikipedia lambda @@ -404,11 +407,11 @@ inline RealType median(const weibull_distribution& dist) } template -inline RealType skewness(const weibull_distribution& dist) +BOOST_MATH_GPU_ENABLED inline RealType skewness(const weibull_distribution& dist) { BOOST_MATH_STD_USING // for ADL of std functions - static const char* function = "boost::math::skewness(const weibull_distribution<%1%>)"; + constexpr auto function = "boost::math::skewness(const weibull_distribution<%1%>)"; RealType shape = dist.shape(); RealType scale = dist.scale(); @@ -429,11 +432,11 @@ inline RealType skewness(const weibull_distribution& dist) } template -inline RealType kurtosis_excess(const weibull_distribution& dist) +BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const weibull_distribution& dist) { BOOST_MATH_STD_USING // for ADL of std functions - static const char* function = "boost::math::kurtosis_excess(const weibull_distribution<%1%>)"; + constexpr auto function = "boost::math::kurtosis_excess(const weibull_distribution<%1%>)"; RealType shape = dist.shape(); RealType scale = dist.scale(); @@ -457,15 +460,15 @@ inline RealType kurtosis_excess(const weibull_distribution& di } template -inline RealType kurtosis(const weibull_distribution& dist) +BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const weibull_distribution& dist) { return kurtosis_excess(dist) + 3; } template -inline RealType entropy(const weibull_distribution& dist) +BOOST_MATH_GPU_ENABLED inline RealType entropy(const weibull_distribution& dist) { - using std::log; + BOOST_MATH_STD_USING RealType k = dist.shape(); RealType lambda = dist.scale(); return constants::euler()*(1-1/k) + log(lambda/k) + 1; diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index f0ac6318b5..00b8e45bf2 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -136,7 +136,7 @@ BOOST_MATH_GPU_ENABLED T beta_imp(T a, T b, const Lanczos&, const Policy& pol) // Generic implementation of Beta(a,b) without Lanczos approximation support // (Caution this is slow!!!): // -#ifndef BOOST_MATH_HAS_NVRTC +#ifndef BOOST_MATH_HAS_GPU_SUPPORT template BOOST_MATH_GPU_ENABLED T beta_imp(T a, T b, const lanczos::undefined_lanczos& l, const Policy& pol) { @@ -461,7 +461,7 @@ BOOST_MATH_GPU_ENABLED T ibeta_power_terms(T a, // // This version is generic, slow, and does not use the Lanczos approximation. // -#ifndef BOOST_MATH_HAS_NVRTC +#ifndef BOOST_MATH_HAS_GPU_SUPPORT template BOOST_MATH_GPU_ENABLED T ibeta_power_terms(T a, T b, @@ -741,7 +741,7 @@ BOOST_MATH_GPU_ENABLED T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool // // Incomplete Beta series again, this time without Lanczos support: // -#ifndef BOOST_MATH_HAS_NVRTC +#ifndef BOOST_MATH_HAS_GPU_SUPPORT template BOOST_MATH_GPU_ENABLED T ibeta_series(T a, T b, T x, T s0, const boost::math::lanczos::undefined_lanczos& l, bool normalised, T* p_derivative, T y, const Policy& pol) { @@ -958,7 +958,7 @@ struct Pn_size #endif }; -#ifndef BOOST_MATH_HAS_NVRTC +#ifndef BOOST_MATH_HAS_GPU_SUPPORT template BOOST_MATH_GPU_ENABLED T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Policy& pol, bool normalised) { diff --git a/include/boost/math/special_functions/detail/gamma_inva.hpp b/include/boost/math/special_functions/detail/gamma_inva.hpp index 75ac89e433..61c64a6486 100644 --- a/include/boost/math/special_functions/detail/gamma_inva.hpp +++ b/include/boost/math/special_functions/detail/gamma_inva.hpp @@ -1,4 +1,5 @@ // (C) Copyright John Maddock 2006. +// (C) Copyright Matt Borland 2024. // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) @@ -17,7 +18,7 @@ #pragma once #endif -#include +#include #include namespace boost{ namespace math{ namespace detail{ @@ -25,8 +26,8 @@ namespace boost{ namespace math{ namespace detail{ template struct gamma_inva_t { - gamma_inva_t(T z_, T p_, bool invert_) : z(z_), p(p_), invert(invert_) {} - T operator()(T a) + BOOST_MATH_GPU_ENABLED gamma_inva_t(T z_, T p_, bool invert_) : z(z_), p(p_), invert(invert_) {} + BOOST_MATH_GPU_ENABLED T operator()(T a) { return invert ? p - boost::math::gamma_q(a, z, Policy()) : boost::math::gamma_p(a, z, Policy()) - p; } @@ -36,7 +37,7 @@ struct gamma_inva_t }; template -T inverse_poisson_cornish_fisher(T lambda, T p, T q, const Policy& pol) +BOOST_MATH_GPU_ENABLED T inverse_poisson_cornish_fisher(T lambda, T p, T q, const Policy& pol) { BOOST_MATH_STD_USING // mean: @@ -67,7 +68,7 @@ T inverse_poisson_cornish_fisher(T lambda, T p, T q, const Policy& pol) } template -T gamma_inva_imp(const T& z, const T& p, const T& q, const Policy& pol) +BOOST_MATH_GPU_ENABLED T gamma_inva_imp(const T& z, const T& p, const T& q, const Policy& pol) { BOOST_MATH_STD_USING // for ADL of std lib math functions // @@ -151,7 +152,7 @@ T gamma_inva_imp(const T& z, const T& p, const T& q, const Policy& pol) } // namespace detail template -inline typename tools::promote_args::type +BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type gamma_p_inva(T1 x, T2 p, const Policy& pol) { typedef typename tools::promote_args::type result_type; @@ -181,7 +182,7 @@ inline typename tools::promote_args::type } template -inline typename tools::promote_args::type +BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type gamma_q_inva(T1 x, T2 q, const Policy& pol) { typedef typename tools::promote_args::type result_type; @@ -211,14 +212,14 @@ inline typename tools::promote_args::type } template -inline typename tools::promote_args::type +BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type gamma_p_inva(T1 x, T2 p) { return boost::math::gamma_p_inva(x, p, policies::policy<>()); } template -inline typename tools::promote_args::type +BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type gamma_q_inva(T1 x, T2 q) { return boost::math::gamma_q_inva(x, q, policies::policy<>()); diff --git a/include/boost/math/special_functions/detail/igamma_inverse.hpp b/include/boost/math/special_functions/detail/igamma_inverse.hpp index f6bbcd72d5..4efd4f78a3 100644 --- a/include/boost/math/special_functions/detail/igamma_inverse.hpp +++ b/include/boost/math/special_functions/detail/igamma_inverse.hpp @@ -1,4 +1,5 @@ // (C) Copyright John Maddock 2006. +// (C) Copyright Matt Borland 2024. // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) @@ -10,6 +11,8 @@ #pragma once #endif +#include +#include #include #include #include @@ -21,7 +24,7 @@ namespace boost{ namespace math{ namespace detail{ template -T find_inverse_s(T p, T q) +BOOST_MATH_GPU_ENABLED T find_inverse_s(T p, T q) { // // Computation of the Incomplete Gamma Function Ratios and their Inverse @@ -41,8 +44,8 @@ T find_inverse_s(T p, T q) { t = sqrt(-2 * log(q)); } - static const double a[4] = { 3.31125922108741, 11.6616720288968, 4.28342155967104, 0.213623493715853 }; - static const double b[5] = { 1, 6.61053765625462, 6.40691597760039, 1.27364489782223, 0.3611708101884203e-1 }; + BOOST_MATH_STATIC const double a[4] = { 3.31125922108741, 11.6616720288968, 4.28342155967104, 0.213623493715853 }; + BOOST_MATH_STATIC const double b[5] = { 1, 6.61053765625462, 6.40691597760039, 1.27364489782223, 0.3611708101884203e-1 }; T s = t - tools::evaluate_polynomial(a, t) / tools::evaluate_polynomial(b, t); if(p < T(0.5)) s = -s; @@ -50,7 +53,7 @@ T find_inverse_s(T p, T q) } template -T didonato_SN(T a, T x, unsigned N, T tolerance = 0) +BOOST_MATH_GPU_ENABLED T didonato_SN(T a, T x, unsigned N, T tolerance = 0) { // // Computation of the Incomplete Gamma Function Ratios and their Inverse @@ -77,7 +80,7 @@ T didonato_SN(T a, T x, unsigned N, T tolerance = 0) } template -inline T didonato_FN(T p, T a, T x, unsigned N, T tolerance, const Policy& pol) +BOOST_MATH_GPU_ENABLED inline T didonato_FN(T p, T a, T x, unsigned N, T tolerance, const Policy& pol) { // // Computation of the Incomplete Gamma Function Ratios and their Inverse @@ -93,7 +96,7 @@ inline T didonato_FN(T p, T a, T x, unsigned N, T tolerance, const Policy& pol) } template -T find_inverse_gamma(T a, T p, T q, const Policy& pol, bool* p_has_10_digits) +BOOST_MATH_GPU_ENABLED T find_inverse_gamma(T a, T p, T q, const Policy& pol, bool* p_has_10_digits) { // // In order to understand what's going on here, you will @@ -233,7 +236,7 @@ T find_inverse_gamma(T a, T p, T q, const Policy& pol, bool* p_has_10_digits) } else { - T D = (std::max)(T(2), T(a * (a - 1))); + T D = BOOST_MATH_GPU_SAFE_MAX(T(2), T(a * (a - 1))); T lg = boost::math::lgamma(a, pol); T lb = log(q) + lg; if(lb < -D * T(2.3)) @@ -315,7 +318,7 @@ T find_inverse_gamma(T a, T p, T q, const Policy& pol, bool* p_has_10_digits) template struct gamma_p_inverse_func { - gamma_p_inverse_func(T a_, T p_, bool inv) : a(a_), p(p_), invert(inv) + BOOST_MATH_GPU_ENABLED gamma_p_inverse_func(T a_, T p_, bool inv) : a(a_), p(p_), invert(inv) { // // If p is too near 1 then P(x) - p suffers from cancellation @@ -333,7 +336,7 @@ struct gamma_p_inverse_func } } - boost::math::tuple operator()(const T& x)const + BOOST_MATH_GPU_ENABLED boost::math::tuple operator()(const T& x)const { BOOST_FPU_EXCEPTION_GUARD // @@ -395,11 +398,11 @@ struct gamma_p_inverse_func }; template -T gamma_p_inv_imp(T a, T p, const Policy& pol) +BOOST_MATH_GPU_ENABLED T gamma_p_inv_imp(T a, T p, const Policy& pol) { BOOST_MATH_STD_USING // ADL of std functions. - static const char* function = "boost::math::gamma_p_inv<%1%>(%1%, %1%)"; + constexpr auto function = "boost::math::gamma_p_inv<%1%>(%1%, %1%)"; BOOST_MATH_INSTRUMENT_VARIABLE(a); BOOST_MATH_INSTRUMENT_VARIABLE(p); @@ -442,7 +445,9 @@ T gamma_p_inv_imp(T a, T p, const Policy& pol) // // Go ahead and iterate: // - std::uintmax_t max_iter = policies::get_max_root_iterations(); + boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); + + #ifndef BOOST_MATH_HAS_GPU_SUPPORT guess = tools::halley_iterate( detail::gamma_p_inverse_func(a, p, false), guess, @@ -450,6 +455,16 @@ T gamma_p_inv_imp(T a, T p, const Policy& pol) tools::max_value(), digits, max_iter); + #else + guess = tools::newton_raphson_iterate( + detail::gamma_p_inverse_func(a, p, false), + guess, + lower, + tools::max_value(), + digits, + max_iter); + #endif + policies::check_root_iterations(function, max_iter, pol); BOOST_MATH_INSTRUMENT_VARIABLE(guess); if(guess == lower) @@ -458,11 +473,11 @@ T gamma_p_inv_imp(T a, T p, const Policy& pol) } template -T gamma_q_inv_imp(T a, T q, const Policy& pol) +BOOST_MATH_GPU_ENABLED T gamma_q_inv_imp(T a, T q, const Policy& pol) { BOOST_MATH_STD_USING // ADL of std functions. - static const char* function = "boost::math::gamma_q_inv<%1%>(%1%, %1%)"; + constexpr auto function = "boost::math::gamma_q_inv<%1%>(%1%, %1%)"; if(a <= 0) return policies::raise_domain_error(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol); @@ -501,7 +516,9 @@ T gamma_q_inv_imp(T a, T q, const Policy& pol) // // Go ahead and iterate: // - std::uintmax_t max_iter = policies::get_max_root_iterations(); + boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); + + #ifndef BOOST_MATH_HAS_GPU_SUPPORT guess = tools::halley_iterate( detail::gamma_p_inverse_func(a, q, true), guess, @@ -509,6 +526,16 @@ T gamma_q_inv_imp(T a, T q, const Policy& pol) tools::max_value(), digits, max_iter); + #else + guess = tools::newton_raphson_iterate( + detail::gamma_p_inverse_func(a, q, true), + guess, + lower, + tools::max_value(), + digits, + max_iter); + #endif + policies::check_root_iterations(function, max_iter, pol); if(guess == lower) guess = policies::raise_underflow_error(function, "Expected result known to be non-zero, but is smaller than the smallest available number.", pol); @@ -518,7 +545,7 @@ T gamma_q_inv_imp(T a, T q, const Policy& pol) } // namespace detail template -inline typename tools::promote_args::type +BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type gamma_p_inv(T1 a, T2 p, const Policy& pol) { typedef typename tools::promote_args::type result_type; @@ -528,7 +555,7 @@ inline typename tools::promote_args::type } template -inline typename tools::promote_args::type +BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type gamma_q_inv(T1 a, T2 p, const Policy& pol) { typedef typename tools::promote_args::type result_type; @@ -538,14 +565,14 @@ inline typename tools::promote_args::type } template -inline typename tools::promote_args::type +BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type gamma_p_inv(T1 a, T2 p) { return gamma_p_inv(a, p, policies::policy<>()); } template -inline typename tools::promote_args::type +BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type gamma_q_inv(T1 a, T2 p) { return gamma_q_inv(a, p, policies::policy<>()); diff --git a/include/boost/math/special_functions/detail/igamma_large.hpp b/include/boost/math/special_functions/detail/igamma_large.hpp index 1fa13c6921..8e0ad1b0dd 100644 --- a/include/boost/math/special_functions/detail/igamma_large.hpp +++ b/include/boost/math/special_functions/detail/igamma_large.hpp @@ -61,6 +61,7 @@ #endif #include +#include namespace boost{ namespace math{ namespace detail{ @@ -68,7 +69,7 @@ namespace boost{ namespace math{ namespace detail{ // when T is unsuitable to be passed to these routines: // template -inline T igamma_temme_large(T, T, const Policy& /* pol */, const std::integral_constant&) +BOOST_MATH_GPU_ENABLED inline T igamma_temme_large(T, T, const Policy& /* pol */, const boost::math::integral_constant&) { // stub function, should never actually be called BOOST_MATH_ASSERT(0); @@ -78,8 +79,11 @@ inline T igamma_temme_large(T, T, const Policy& /* pol */, const std::integral_c // This version is accurate for up to 64-bit mantissa's, // (80-bit long double, or 10^-20). // + +#ifndef BOOST_MATH_HAS_GPU_SUPPORT + template -T igamma_temme_large(T a, T x, const Policy& pol, const std::integral_constant&) +BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const boost::math::integral_constant&) { BOOST_MATH_STD_USING // ADL of std functions T sigma = (x - a) / a; @@ -91,7 +95,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, const std::integral_constant -BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const std::integral_constant&) +BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const boost::math::integral_constant&) { BOOST_MATH_STD_USING // ADL of std functions T sigma = (x - a) / a; @@ -417,7 +424,18 @@ BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const s if(x < a) result = -result; + #ifdef BOOST_MATH_HAS_NVRTC + if (boost::math::is_same_v) + { + result += ::erfcf(::sqrtf(y)) / 2; + } + else + { + result += ::erfc(::sqrt(y)) / 2; + } + #else result += boost::math::erfc(sqrt(y), pol) / 2; + #endif return result; } @@ -426,7 +444,7 @@ BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const s // (IEEE float precision, or 10^-8) // template -BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const std::integral_constant&) +BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const boost::math::integral_constant&) { BOOST_MATH_STD_USING // ADL of std functions T sigma = (x - a) / a; @@ -470,7 +488,18 @@ BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const s if(x < a) result = -result; + #ifdef BOOST_MATH_HAS_NVRTC + if (boost::math::is_same_v) + { + result += ::erfcf(::sqrtf(y)) / 2; + } + else + { + result += ::erfc(::sqrt(y)) / 2; + } + #else result += boost::math::erfc(sqrt(y), pol) / 2; + #endif return result; } @@ -481,8 +510,10 @@ BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const s // It's use for a < 200 is not recommended, that would // require many more terms in the polynomials. // +#ifndef BOOST_MATH_HAS_GPU_SUPPORT + template -T igamma_temme_large(T a, T x, const Policy& pol, const std::integral_constant&) +BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const boost::math::integral_constant&) { BOOST_MATH_STD_USING // ADL of std functions T sigma = (x - a) / a; @@ -494,7 +525,7 @@ T igamma_temme_large(T a, T x, const Policy& pol, const std::integral_constant #include +#include +#include +#include #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128) // @@ -32,13 +35,13 @@ namespace boost{ namespace math{ namespace detail{ template BOOST_MATH_GPU_ENABLED T gamma_imp(T z, const Policy& pol, const Lanczos& l); template -T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l); +BOOST_MATH_GPU_ENABLED T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l); // // lgamma for small arguments: // template -BOOST_MATH_GPU_ENABLED T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, const Policy& /* l */, const Lanczos&) +BOOST_MATH_GPU_ENABLED T lgamma_small_imp(T z, T zm1, T zm2, const boost::math::integral_constant&, const Policy& /* l */, const Lanczos&) { // This version uses rational approximations for small // values of z accurate enough for 64-bit mantissas @@ -226,8 +229,10 @@ BOOST_MATH_GPU_ENABLED T lgamma_small_imp(T z, T zm1, T zm2, const std::integral } return result; } + +#ifndef BOOST_MATH_HAS_GPU_SUPPORT template -T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, const Policy& /* l */, const Lanczos&) +T lgamma_small_imp(T z, T zm1, T zm2, const boost::math::integral_constant&, const Policy& /* l */, const Lanczos&) { // // This version uses rational approximations for small @@ -484,7 +489,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, c return result; } template -T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, const Policy& pol, const Lanczos& l) +BOOST_MATH_GPU_ENABLED T lgamma_small_imp(T z, T zm1, T zm2, const boost::math::integral_constant&, const Policy& pol, const Lanczos& l) { // // No rational approximations are available because either @@ -528,6 +533,8 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant&, con return result; } +#endif // BOOST_MATH_HAS_GPU_SUPPORT + }}} // namespaces #endif // BOOST_MATH_SPECIAL_FUNCTIONS_DETAIL_LGAMMA_SMALL diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 9268ba4156..186befb612 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -15,14 +15,13 @@ #endif #include - -#ifndef BOOST_MATH_HAS_NVRTC - #include #include #include #include #include +#include +#include #include #include #include @@ -36,12 +35,12 @@ #include #include #include + +// Only needed for types larger than double +#ifndef BOOST_MATH_HAS_GPU_SUPPORT #include #include - -#include -#include -#include +#endif #ifdef _MSC_VER # pragma warning(push) @@ -60,13 +59,13 @@ namespace boost{ namespace math{ namespace detail{ template -BOOST_MATH_GPU_ENABLED inline bool is_odd(T v, const std::true_type&) +BOOST_MATH_GPU_ENABLED inline bool is_odd(T v, const boost::math::true_type&) { int i = static_cast(v); return i&1; } template -BOOST_MATH_GPU_ENABLED inline bool is_odd(T v, const std::false_type&) +BOOST_MATH_GPU_ENABLED inline bool is_odd(T v, const boost::math::false_type&) { // Oh dear can't cast T to int! BOOST_MATH_STD_USING @@ -76,7 +75,7 @@ BOOST_MATH_GPU_ENABLED inline bool is_odd(T v, const std::false_type&) template BOOST_MATH_GPU_ENABLED inline bool is_odd(T v) { - return is_odd(v, ::std::is_convertible()); + return is_odd(v, ::boost::math::is_convertible()); } template @@ -208,7 +207,7 @@ BOOST_MATH_GPU_ENABLED BOOST_MATH_FORCEINLINE T gamma_imp(T z, const Policy& pol result = -boost::math::constants::pi() / result; if(result == 0) return policies::raise_underflow_error(function, "Result of tgamma is too small to represent.", pol); - if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL) + if((boost::math::fpclassify)(result) == (int)BOOST_MATH_FP_SUBNORMAL) return policies::raise_denorm_error(function, "Result of tgamma is denormalized.", result, pol); BOOST_MATH_INSTRUMENT_VARIABLE(result); return result; @@ -259,7 +258,7 @@ BOOST_MATH_GPU_ENABLED T lgamma_imp_final(T z, const Policy& pol, const Lanczos& else if(z < 15) { typedef typename policies::precision::type precision_type; - typedef std::integral_constant(z, T(z - 1), T(z - 2), tag_type(), pol, l); } - else if((z >= 3) && (z < 100) && (std::numeric_limits::max_exponent >= 1024)) + else if((z >= 3) && (z < 100) && (boost::math::numeric_limits::max_exponent >= 1024)) { // taking the log of tgamma reduces the error, no danger of overflow here: result = log(gamma_imp(z, pol, l)); @@ -349,7 +348,7 @@ struct upper_incomplete_gamma_fract T z, a; int k; public: - typedef std::pair result_type; + typedef boost::math::pair result_type; BOOST_MATH_GPU_ENABLED upper_incomplete_gamma_fract(T a1, T z1) : z(z1-a1+1), a(a1), k(0) @@ -399,26 +398,28 @@ BOOST_MATH_GPU_ENABLED inline T lower_gamma_series(T a, T z, const Policy& pol, // lower incomplete integral. Then divide by tgamma(a) // to get the normalised value. lower_incomplete_gamma_series s(a, z); - std::uintmax_t max_iter = policies::get_max_series_iterations(); + boost::math::uintmax_t max_iter = policies::get_max_series_iterations(); T factor = policies::get_epsilon(); T result = boost::math::tools::sum_series(s, factor, max_iter, init_value); policies::check_series_iterations("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol); return result; } +#ifndef BOOST_MATH_HAS_GPU_SUPPORT + // // Fully generic tgamma and lgamma use Stirling's approximation // with Bernoulli numbers. // template -std::size_t highest_bernoulli_index() +boost::math::size_t highest_bernoulli_index() { - const float digits10_of_type = (std::numeric_limits::is_specialized - ? static_cast(std::numeric_limits::digits10) + const float digits10_of_type = (boost::math::numeric_limits::is_specialized + ? static_cast(boost::math::numeric_limits::digits10) : static_cast(boost::math::tools::digits() * 0.301F)); // Find the high index n for Bn to produce the desired precision in Stirling's calculation. - return static_cast(18.0F + (0.6F * digits10_of_type)); + return static_cast(18.0F + (0.6F * digits10_of_type)); } template @@ -426,8 +427,8 @@ int minimum_argument_for_bernoulli_recursion() { BOOST_MATH_STD_USING - const float digits10_of_type = (std::numeric_limits::is_specialized - ? (float) std::numeric_limits::digits10 + const float digits10_of_type = (boost::math::numeric_limits::is_specialized + ? (float) boost::math::numeric_limits::digits10 : (float) (boost::math::tools::digits() * 0.301F)); int min_arg = (int) (digits10_of_type * 1.7F); @@ -449,7 +450,7 @@ int minimum_argument_for_bernoulli_recursion() const float d2_minus_one = ((digits10_of_type / 0.301F) - 1.0F); const float limit = ceil(exp((d2_minus_one * log(2.0F)) / 20.0F)); - min_arg = (int) ((std::min)(digits10_of_type * 1.7F, limit)); + min_arg = (int) (BOOST_MATH_GPU_SAFE_MIN(digits10_of_type * 1.7F, limit)); } return min_arg; @@ -468,7 +469,7 @@ T scaled_tgamma_no_lanczos(const T& z, const Policy& pol, bool islog = false) // Perform the Bernoulli series expansion of Stirling's approximation. - const std::size_t number_of_bernoullis_b2n = policies::get_max_series_iterations(); + const boost::math::size_t number_of_bernoullis_b2n = policies::get_max_series_iterations(); T one_over_x_pow_two_n_minus_one = 1 / z; const T one_over_x2 = one_over_x_pow_two_n_minus_one * one_over_x_pow_two_n_minus_one; @@ -477,11 +478,11 @@ T scaled_tgamma_no_lanczos(const T& z, const Policy& pol, bool islog = false) const T half_ln_two_pi_over_z = sqrt(boost::math::constants::two_pi() / z); T last_term = 2 * sum; - for (std::size_t n = 2U;; ++n) + for (boost::math::size_t n = 2U;; ++n) { one_over_x_pow_two_n_minus_one *= one_over_x2; - const std::size_t n2 = static_cast(n * 2U); + const boost::math::size_t n2 = static_cast(n * 2U); const T term = (boost::math::bernoulli_b2n(static_cast(n)) * one_over_x_pow_two_n_minus_one) / (n2 * (n2 - 1U)); @@ -629,7 +630,7 @@ T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&) if(gamma_value == 0) return policies::raise_underflow_error(function, "Result of tgamma is too small to represent.", pol); - if((boost::math::fpclassify)(gamma_value) == static_cast(FP_SUBNORMAL)) + if((boost::math::fpclassify)(gamma_value) == static_cast(BOOST_MATH_FP_SUBNORMAL)) return policies::raise_denorm_error(function, "Result of tgamma is denormalized.", gamma_value, pol); } @@ -775,6 +776,21 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sig return log_gamma_value; } +#endif // BOOST_MATH_HAS_GPU_SUPPORT + +// In order for tgammap1m1_imp to compile we need a forward decl of boost::math::tgamma +// The rub is that we can't just use math_fwd so we provide one here only in that circumstance +#ifdef BOOST_MATH_HAS_NVRTC +template +BOOST_MATH_GPU_ENABLED tools::promote_args_t tgamma(RT z); + +template +BOOST_MATH_GPU_ENABLED tools::promote_args_t tgamma(RT1 a, RT2 z); + +template +BOOST_MATH_GPU_ENABLED tools::promote_args_t tgamma(RT1 a, RT2 z, const Policy& pol); +#endif + // // This helper calculates tgamma(dz+1)-1 without cancellation errors, // used by the upper incomplete gamma with z < 1: @@ -786,7 +802,7 @@ BOOST_MATH_GPU_ENABLED T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& typedef typename policies::precision::type precision_type; - typedef std::integral_constant inline T tgammap1m1_imp(T z, Policy const& pol, const ::boost::math::lanczos::undefined_lanczos&) @@ -841,6 +867,8 @@ inline T tgammap1m1_imp(T z, Policy const& pol, return boost::math::expm1(boost::math::lgamma(1 + z, pol)); } +#endif // BOOST_MATH_HAS_GPU_SUPPORT + // // Series representation for upper fraction when z is small: // @@ -914,7 +942,7 @@ BOOST_MATH_GPU_ENABLED T full_igamma_prefix(T a, T z, const Policy& pol) // This error handling isn't very good: it happens after the fact // rather than before it... // - if((boost::math::fpclassify)(prefix) == (int)FP_INFINITE) + if((boost::math::fpclassify)(prefix) == (int)BOOST_MATH_FP_INFINITE) return policies::raise_overflow_error("boost::math::detail::full_igamma_prefix<%1%>(%1%, %1%)", "Result of incomplete gamma function is too large to represent.", pol); return prefix; @@ -971,16 +999,16 @@ BOOST_MATH_GPU_ENABLED T regularised_gamma_prefix(T a, T z, const Policy& pol, c // T alz = a * log(z / agh); T amz = a - z; - if(((std::min)(alz, amz) <= tools::log_min_value()) || ((std::max)(alz, amz) >= tools::log_max_value())) + if((BOOST_MATH_GPU_SAFE_MIN(alz, amz) <= tools::log_min_value()) || (BOOST_MATH_GPU_SAFE_MAX(alz, amz) >= tools::log_max_value())) { T amza = amz / a; - if(((std::min)(alz, amz)/2 > tools::log_min_value()) && ((std::max)(alz, amz)/2 < tools::log_max_value())) + if((BOOST_MATH_GPU_SAFE_MIN(alz, amz)/2 > tools::log_min_value()) && (BOOST_MATH_GPU_SAFE_MAX(alz, amz)/2 < tools::log_max_value())) { // compute square root of the result and then square it: T sq = pow(z / agh, a / 2) * exp(amz / 2); prefix = sq * sq; } - else if(((std::min)(alz, amz)/4 > tools::log_min_value()) && ((std::max)(alz, amz)/4 < tools::log_max_value()) && (z > a)) + else if((BOOST_MATH_GPU_SAFE_MIN(alz, amz)/4 > tools::log_min_value()) && (BOOST_MATH_GPU_SAFE_MAX(alz, amz)/4 < tools::log_max_value()) && (z > a)) { // compute the 4th root of the result then square it twice: T sq = pow(z / agh, a / 4) * exp(amz / 4); @@ -1004,6 +1032,9 @@ BOOST_MATH_GPU_ENABLED T regularised_gamma_prefix(T a, T z, const Policy& pol, c prefix *= sqrt(agh / boost::math::constants::e()) / Lanczos::lanczos_sum_expG_scaled(a); return prefix; } + +#ifndef BOOST_MATH_HAS_GPU_SUPPORT + // // And again, without Lanczos support: // @@ -1073,6 +1104,9 @@ T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined } } } + +#endif // BOOST_MATH_HAS_GPU_SUPPORT + // // Upper gamma fraction for very small a: // @@ -1084,7 +1118,14 @@ BOOST_MATH_GPU_ENABLED inline T tgamma_small_upper_part(T a, T x, const Policy& // Compute the full upper fraction (Q) when a is very small: // + #ifdef BOOST_MATH_HAS_NVRTC + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename lanczos::lanczos::type evaluation_type; + T result {detail::tgammap1m1_imp(static_cast(a), pol, evaluation_type())}; + #else T result { boost::math::tgamma1pm1(a, pol) }; + #endif if(pgam) *pgam = (result + 1) / a; @@ -1092,7 +1133,7 @@ BOOST_MATH_GPU_ENABLED inline T tgamma_small_upper_part(T a, T x, const Policy& result -= p; result /= a; detail::small_gamma2_series s(a, x); - std::uintmax_t max_iter = policies::get_max_series_iterations() - 10; + boost::math::uintmax_t max_iter = policies::get_max_series_iterations() - 10; p += 1; if(pderivative) *pderivative = p / (*pgam * exp(x)); @@ -1141,7 +1182,21 @@ BOOST_MATH_GPU_ENABLED T finite_half_gamma_q(T a, T x, T* p_derivative, const Po // Calculates normalised Q when a is a half-integer: // BOOST_MATH_STD_USING + + #ifdef BOOST_MATH_HAS_NVRTC + T e; + if (boost::math::is_same_v) + { + e = ::erfcf(::sqrtf(x)); + } + else + { + e = ::erfc(::sqrt(x)); + } + #else T e = boost::math::erfc(sqrt(x), pol); + #endif + if((e != 0) && (a > 1)) { T term = exp(-x) / sqrt(constants::pi() * x); @@ -1192,7 +1247,7 @@ BOOST_MATH_GPU_ENABLED T incomplete_tgamma_large_x(const T& a, const T& x, const { BOOST_MATH_STD_USING incomplete_tgamma_large_x_series s(a, x); - std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations(); + boost::math::uintmax_t max_iter = boost::math::policies::get_max_series_iterations(); T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon(), max_iter); boost::math::policies::check_series_iterations("boost::math::tgamma<%1%>(%1%,%1%)", max_iter, pol); return result; @@ -1203,7 +1258,7 @@ BOOST_MATH_GPU_ENABLED T incomplete_tgamma_large_x(const T& a, const T& x, const // Main incomplete gamma entry point, handles all four incomplete gamma's: // template -BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, +BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp_final(T a, T x, bool normalised, bool invert, const Policy& pol, T* p_derivative) { constexpr auto function = "boost::math::gamma_p<%1%>(%1%, %1%)"; @@ -1218,70 +1273,6 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool in T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used - if(a >= max_factorial::value && !normalised) - { - // - // When we're computing the non-normalized incomplete gamma - // and a is large the result is rather hard to compute unless - // we use logs. There are really two options - if x is a long - // way from a in value then we can reliably use methods 2 and 4 - // below in logarithmic form and go straight to the result. - // Otherwise we let the regularized gamma take the strain - // (the result is unlikely to underflow in the central region anyway) - // and combine with lgamma in the hopes that we get a finite result. - // - if(invert && (a * 4 < x)) - { - // This is method 4 below, done in logs: - result = a * log(x) - x; - if(p_derivative) - *p_derivative = exp(result); - result += log(upper_gamma_fraction(a, x, policies::get_epsilon())); - } - else if(!invert && (a > 4 * x)) - { - // This is method 2 below, done in logs: - result = a * log(x) - x; - if(p_derivative) - *p_derivative = exp(result); - T init_value = 0; - result += log(detail::lower_gamma_series(a, x, pol, init_value) / a); - } - else - { - result = gamma_incomplete_imp(a, x, true, invert, pol, p_derivative); - if(result == 0) - { - if(invert) - { - // Try http://functions.wolfram.com/06.06.06.0039.01 - result = 1 + 1 / (12 * a) + 1 / (288 * a * a); - result = log(result) - a + (a - 0.5f) * log(a) + log(boost::math::constants::root_two_pi()); - if(p_derivative) - *p_derivative = exp(a * log(x) - x); - } - else - { - // This is method 2 below, done in logs, we're really outside the - // range of this method, but since the result is almost certainly - // infinite, we should probably be OK: - result = a * log(x) - x; - if(p_derivative) - *p_derivative = exp(result); - T init_value = 0; - result += log(detail::lower_gamma_series(a, x, pol, init_value) / a); - } - } - else - { - result = log(result) + boost::math::lgamma(a, pol); - } - } - if(result > tools::log_max_value()) - return policies::raise_overflow_error(function, nullptr, pol); - return exp(result); - } - BOOST_MATH_ASSERT((p_derivative == nullptr) || normalised); bool is_int, is_half_int; @@ -1357,7 +1348,7 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool in // series and continued fractions are slow to converge: // bool use_temme = false; - if(normalised && std::numeric_limits::is_specialized && (a > 20)) + if(normalised && boost::math::numeric_limits::is_specialized && (a > 20)) { T sigma = fabs((x-a)/a); if((a > 200) && (policies::digits() <= 113)) @@ -1414,14 +1405,40 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool in { result = finite_gamma_q(a, x, pol, p_derivative); if(!normalised) + { + #ifdef BOOST_MATH_HAS_NVRTC + if (boost::math::is_same_v) + { + result *= ::tgammaf(a); + } + else + { + result *= ::tgamma(a); + } + #else result *= boost::math::tgamma(a, pol); + #endif + } break; } case 1: { result = finite_half_gamma_q(a, x, p_derivative, pol); if(!normalised) + { + #ifdef BOOST_MATH_HAS_NVRTC + if (boost::math::is_same_v) + { + result *= ::tgammaf(a); + } + else + { + result *= ::tgamma(a); + } + #else result *= boost::math::tgamma(a, pol); + #endif + } if(p_derivative && (*p_derivative == 0)) *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type()); break; @@ -1450,7 +1467,19 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool in bool optimised_invert = false; if(invert) { + #ifdef BOOST_MATH_HAS_NVRTC + if (boost::math::is_same_v) + { + init_value = (normalised ? 1 : ::tgammaf(a)); + } + else + { + init_value = (normalised ? 1 : ::tgamma(a)); + } + #else init_value = (normalised ? 1 : boost::math::tgamma(a, pol)); + #endif + if(normalised || (result >= 1) || (tools::max_value() * result > init_value)) { init_value /= result; @@ -1507,7 +1536,7 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool in // typedef typename policies::precision::type precision_type; - typedef std::integral_constant) + { + result = ::powf(x, a) / ::tgammaf(a + 1); + } + else + { + result = ::pow(x, a) / ::tgamma(a + 1); + } + #else result = pow(x, a) / boost::math::tgamma(a + 1, pol); + #endif #ifndef BOOST_MATH_NO_EXCEPTIONS } catch (const std::overflow_error&) @@ -1565,7 +1605,19 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool in result = 1; if(invert) { + #ifdef BOOST_MATH_HAS_NVRTC + T gam; + if (boost::math::is_same_v) + { + gam = normalised ? 1 : ::tgammaf(a); + } + else + { + gam = normalised ? 1 : ::tgamma(a); + } + #else T gam = normalised ? 1 : boost::math::tgamma(a, pol); + #endif result = gam - result; } if(p_derivative) @@ -1585,6 +1637,101 @@ BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool in return result; } +// Need to implement this dispatch to avoid recursion for device compilers +template +BOOST_MATH_GPU_ENABLED T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, + const Policy& pol, T* p_derivative) +{ + constexpr auto function = "boost::math::gamma_p<%1%>(%1%, %1%)"; + if(a <= 0) + return policies::raise_domain_error(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol); + if(x < 0) + return policies::raise_domain_error(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol); + + BOOST_MATH_STD_USING + + + T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used + + if(a >= max_factorial::value && !normalised) + { + // + // When we're computing the non-normalized incomplete gamma + // and a is large the result is rather hard to compute unless + // we use logs. There are really two options - if x is a long + // way from a in value then we can reliably use methods 2 and 4 + // below in logarithmic form and go straight to the result. + // Otherwise we let the regularized gamma take the strain + // (the result is unlikely to underflow in the central region anyway) + // and combine with lgamma in the hopes that we get a finite result. + // + if(invert && (a * 4 < x)) + { + // This is method 4 below, done in logs: + result = a * log(x) - x; + if(p_derivative) + *p_derivative = exp(result); + result += log(upper_gamma_fraction(a, x, policies::get_epsilon())); + } + else if(!invert && (a > 4 * x)) + { + // This is method 2 below, done in logs: + result = a * log(x) - x; + if(p_derivative) + *p_derivative = exp(result); + T init_value = 0; + result += log(detail::lower_gamma_series(a, x, pol, init_value) / a); + } + else + { + result = gamma_incomplete_imp_final(T(a), T(x), true, invert, pol, p_derivative); + if(result == 0) + { + if(invert) + { + // Try http://functions.wolfram.com/06.06.06.0039.01 + result = 1 + 1 / (12 * a) + 1 / (288 * a * a); + result = log(result) - a + (a - 0.5f) * log(a) + log(boost::math::constants::root_two_pi()); + if(p_derivative) + *p_derivative = exp(a * log(x) - x); + } + else + { + // This is method 2 below, done in logs, we're really outside the + // range of this method, but since the result is almost certainly + // infinite, we should probably be OK: + result = a * log(x) - x; + if(p_derivative) + *p_derivative = exp(result); + T init_value = 0; + result += log(detail::lower_gamma_series(a, x, pol, init_value) / a); + } + } + else + { + #ifdef BOOST_MATH_HAS_NVRTC + if (boost::math::is_same_v) + { + result = ::logf(result) + ::lgammaf(a); + } + else + { + result = ::log(result) + ::lgamma(a); + } + #else + result = log(result) + boost::math::lgamma(a, pol); + #endif + } + } + if(result > tools::log_max_value()) + return policies::raise_overflow_error(function, nullptr, pol); + return exp(result); + } + + // If no special handling is required then we proceeds as normal + return gamma_incomplete_imp_final(T(a), T(x), normalised, invert, pol, p_derivative); +} + // // Ratios of two gamma functions: // @@ -1612,7 +1759,18 @@ BOOST_MATH_GPU_ENABLED T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Poli } else { + #ifdef BOOST_MATH_HAS_NVRTC + if (boost::math::is_same_v) + { + return 1 / (z * ::tgammaf(z + delta)); + } + else + { + return 1 / (z * ::tgamma(z + delta)); + } + #else return 1 / (z * boost::math::tgamma(z + delta, pol)); + #endif } } T zgh = static_cast(z + T(Lanczos::g()) - constants::half()); @@ -1651,6 +1809,8 @@ BOOST_MATH_GPU_ENABLED T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Poli // // And again without Lanczos support this time: // +#ifndef BOOST_MATH_HAS_GPU_SUPPORT + template T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos& l) { @@ -1707,6 +1867,8 @@ T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos: return ratio; } +#endif + template BOOST_MATH_GPU_ENABLED T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol) { @@ -1715,7 +1877,18 @@ BOOST_MATH_GPU_ENABLED T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol) if((z <= 0) || (z + delta <= 0)) { // This isn't very sophisticated, or accurate, but it does work: + #ifdef BOOST_MATH_HAS_NVRTC + if (boost::math::is_same_v) + { + return ::tgammaf(z) / ::tgammaf(z + delta); + } + else + { + return ::tgamma(z) / ::tgamma(z + delta); + } + #else return boost::math::tgamma(z, pol) / boost::math::tgamma(z + delta, pol); + #endif } if(floor(delta) == delta) @@ -1775,17 +1948,32 @@ BOOST_MATH_GPU_ENABLED T tgamma_ratio_imp(T x, T y, const Policy& pol) if((y <= 0) || (boost::math::isinf)(y)) return policies::raise_domain_error("boost::math::tgamma_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got b=%1%).", y, pol); + // We don't need to worry about the denorm case on device + // And this has the added bonus of removing recursion + #ifndef BOOST_MATH_HAS_GPU_SUPPORT if(x <= tools::min_value()) { // Special case for denorms...Ugh. T shift = ldexp(T(1), tools::digits()); return shift * tgamma_ratio_imp(T(x * shift), y, pol); } + #endif if((x < max_factorial::value) && (y < max_factorial::value)) { // Rather than subtracting values, lets just call the gamma functions directly: + #ifdef BOOST_MATH_HAS_NVRTC + if (boost::math::is_same_v) + { + return ::tgammaf(x) / ::tgammaf(y); + } + else + { + return ::tgamma(x) / ::tgamma(y); + } + #else return boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol); + #endif } T prefix = 1; if(x < 1) @@ -1801,12 +1989,35 @@ BOOST_MATH_GPU_ENABLED T tgamma_ratio_imp(T x, T y, const Policy& pol) y -= 1; prefix /= y; } + + #ifdef BOOST_MATH_HAS_NVRTC + if (boost::math::is_same_v) + { + return prefix * ::tgammaf(x) / ::tgammaf(y); + } + else + { + return prefix * ::tgamma(x) / ::tgamma(y); + } + #else return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol); + #endif } // // result is almost certainly going to underflow to zero, try logs just in case: // + #ifdef BOOST_MATH_HAS_NVRTC + if (boost::math::is_same_v) + { + return ::expf(::lgammaf(x) - ::lgammaf(y)); + } + else + { + return ::exp(::lgamma(x) - ::lgamma(y)); + } + #else return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol)); + #endif } if(y < 1) { @@ -1821,17 +2032,44 @@ BOOST_MATH_GPU_ENABLED T tgamma_ratio_imp(T x, T y, const Policy& pol) x -= 1; prefix *= x; } + + #ifdef BOOST_MATH_HAS_NVRTC + if (boost::math::is_same_v) + { + return prefix * ::tgammaf(x) / ::tgammaf(y); + } + else + { + return prefix * ::tgamma(x) / ::tgamma(y); + } + #else return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol); + #endif } // // Result will almost certainly overflow, try logs just in case: // + #ifdef BOOST_MATH_HAS_NVRTC + if (boost::math::is_same_v) + { + return ::expf(::lgammaf(x) - ::lgammaf(y)); + } + else + { + return ::exp(::lgamma(x) - ::lgamma(y)); + } + #else return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol)); + #endif } // // Regular case, x and y both large and similar in magnitude: // + #ifdef BOOST_MATH_HAS_NVRTC + return detail::tgamma_delta_ratio_imp(x, y - x, pol); + #else return boost::math::tgamma_delta_ratio(x, y - x, pol); + #endif } template @@ -1866,7 +2104,18 @@ BOOST_MATH_GPU_ENABLED T gamma_p_derivative_imp(T a, T x, const Policy& pol) if(f1 == 0) { // Underflow in calculation, use logs instead: + #ifdef BOOST_MATH_HAS_NVRTC + if (boost::math::is_same_v) + { + f1 = a * ::logf(x) - x - ::lgammaf(a) - ::logf(x); + } + else + { + f1 = a * ::log(x) - x - ::lgamma(a) - ::log(x); + } + #else f1 = a * log(x) - x - lgamma(a, pol) - log(x); + #endif f1 = exp(f1); } else @@ -1901,7 +2150,7 @@ struct igamma_initializer { typedef typename policies::precision::type precision_type; - typedef std::integral_constant - BOOST_MATH_GPU_ENABLED static void do_init(const std::integral_constant&) + BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant&) { // If std::numeric_limits::digits is zero, we must not call // our initialization code here as the precision presumably // varies at runtime, and will not have been set yet. Plus the // code requiring initialization isn't called when digits == 0. - if (std::numeric_limits::digits) + if (boost::math::numeric_limits::digits) { boost::math::gamma_p(static_cast(400), static_cast(400), Policy()); } } - BOOST_MATH_GPU_ENABLED static void do_init(const std::integral_constant&){} - void force_instantiate()const{} + BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant&){} + BOOST_MATH_GPU_ENABLED void force_instantiate()const{} }; BOOST_MATH_STATIC const init initializer; BOOST_MATH_GPU_ENABLED static void force_instantiate() @@ -1945,7 +2194,7 @@ struct lgamma_initializer BOOST_MATH_GPU_ENABLED init() { typedef typename policies::precision::type precision_type; - typedef std::integral_constant&) + BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant&) { boost::math::lgamma(static_cast(2.5), Policy()); boost::math::lgamma(static_cast(1.25), Policy()); boost::math::lgamma(static_cast(1.75), Policy()); } - BOOST_MATH_GPU_ENABLED static void do_init(const std::integral_constant&) + BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant&) { boost::math::lgamma(static_cast(2.5), Policy()); boost::math::lgamma(static_cast(1.25), Policy()); boost::math::lgamma(static_cast(1.5), Policy()); boost::math::lgamma(static_cast(1.75), Policy()); } - BOOST_MATH_GPU_ENABLED static void do_init(const std::integral_constant&) + BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant&) { } BOOST_MATH_GPU_ENABLED void force_instantiate()const{} @@ -2016,13 +2265,6 @@ BOOST_MATH_GPU_ENABLED inline tools::promote_args_t } // namespace detail -template -BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type - tgamma(T z) -{ - return tgamma(z, policies::policy<>()); -} - template BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type lgamma(T z, int* sign, const Policy&) @@ -2079,7 +2321,7 @@ BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - return policies::checked_narrowing_cast::type, forwarding_policy>(detail::tgammap1m1_imp(static_cast(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma1pm1<%!%>(%1%)"); + return policies::checked_narrowing_cast::type, forwarding_policy>(detail::tgammap1m1_imp(static_cast(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma1pm1<%!%>(%1%)"); } template @@ -2111,6 +2353,12 @@ BOOST_MATH_GPU_ENABLED inline tools::promote_args_t using result_type = tools::promote_args_t; return static_cast(detail::tgamma(a, z, pol, boost::math::false_type())); } +template +BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type + tgamma(T z) +{ + return tgamma(z, policies::policy<>()); +} // // Full lower incomplete gamma: // @@ -2284,73 +2532,4 @@ BOOST_MATH_GPU_ENABLED inline tools::promote_args_t #include #include -#else - -#include -#include - -namespace boost { -namespace math { - -inline BOOST_MATH_GPU_ENABLED float tgamma(float x) { return ::tgammaf(x); } -inline BOOST_MATH_GPU_ENABLED double tgamma(double x) { return ::tgamma(x); } - -template -BOOST_MATH_GPU_ENABLED T tgamma(T x, const Policy&) -{ - return boost::math::tgamma(x); -} - -inline BOOST_MATH_GPU_ENABLED float lgamma(float x) { return ::lgammaf(x); } -inline BOOST_MATH_GPU_ENABLED double lgamma(double x) { return ::lgamma(x); } - -template -BOOST_MATH_GPU_ENABLED T lgamma(T x, const Policy&) -{ - return boost::math::lgamma(x); -} - -template -BOOST_MATH_GPU_ENABLED T lgamma(T x, int* sign, const Policy&) -{ - auto res = boost::math::lgamma(x); - if (sign != nullptr) - { - if (res < 0) - { - *sign = -1; - } - else - { - *sign = 1; - } - } - - return res; -} - -template -BOOST_MATH_GPU_ENABLED T tgamma1pm1(T z) -{ - using namespace boost::math; - - if (fabs(z) < T(0.55)) - { - return expm1(lgamma(z)); - } - - return expm1(lgamma(1 + z)); -} - -template -BOOST_MATH_GPU_ENABLED T tgamma1pm1(T x, const Policy&) -{ - return tgamma1pm1(x); -} - -} // namespace math -} // namespace boost - -#endif // __CUDACC_RTC__ - #endif // BOOST_MATH_SF_GAMMA_HPP diff --git a/include/boost/math/special_functions/log1p.hpp b/include/boost/math/special_functions/log1p.hpp index cdec8ee86c..758f606687 100644 --- a/include/boost/math/special_functions/log1p.hpp +++ b/include/boost/math/special_functions/log1p.hpp @@ -13,15 +13,13 @@ #endif #include - -#ifndef BOOST_MATH_HAS_NVRTC - -#include -#include -#include #include #include #include +#include +#include +#include +#include #include #include #include @@ -82,7 +80,7 @@ namespace detail // it performs no better than log(1+x): which is to say not very well at all. // template -BOOST_MATH_GPU_ENABLED T log1p_imp(T const & x, const Policy& pol, const std::integral_constant&) +BOOST_MATH_GPU_ENABLED T log1p_imp(T const & x, const Policy& pol, const boost::math::integral_constant&) { // The function returns the natural logarithm of 1 + x. typedef typename tools::promote_args::type result_type; BOOST_MATH_STD_USING @@ -104,7 +102,7 @@ BOOST_MATH_GPU_ENABLED T log1p_imp(T const & x, const Policy& pol, const std::in if(a < tools::epsilon()) return x; detail::log1p_series s(x); - std::uintmax_t max_iter = policies::get_max_series_iterations(); + boost::math::uintmax_t max_iter = policies::get_max_series_iterations(); result_type result = tools::sum_series(s, policies::get_epsilon(), max_iter); @@ -113,7 +111,7 @@ BOOST_MATH_GPU_ENABLED T log1p_imp(T const & x, const Policy& pol, const std::in } template -BOOST_MATH_GPU_ENABLED T log1p_imp(T const& x, const Policy& pol, const std::integral_constant&) +BOOST_MATH_GPU_ENABLED T log1p_imp(T const& x, const Policy& pol, const boost::math::integral_constant&) { // The function returns the natural logarithm of 1 + x. BOOST_MATH_STD_USING @@ -166,7 +164,7 @@ BOOST_MATH_GPU_ENABLED T log1p_imp(T const& x, const Policy& pol, const std::int } template -BOOST_MATH_GPU_ENABLED T log1p_imp(T const& x, const Policy& pol, const std::integral_constant&) +BOOST_MATH_GPU_ENABLED T log1p_imp(T const& x, const Policy& pol, const boost::math::integral_constant&) { // The function returns the natural logarithm of 1 + x. BOOST_MATH_STD_USING @@ -221,7 +219,7 @@ BOOST_MATH_GPU_ENABLED T log1p_imp(T const& x, const Policy& pol, const std::int } template -BOOST_MATH_GPU_ENABLED T log1p_imp(T const& x, const Policy& pol, const std::integral_constant&) +BOOST_MATH_GPU_ENABLED T log1p_imp(T const& x, const Policy& pol, const boost::math::integral_constant&) { // The function returns the natural logarithm of 1 + x. BOOST_MATH_STD_USING @@ -276,8 +274,8 @@ struct log1p_initializer do_init(tag()); } template - BOOST_MATH_GPU_ENABLED static void do_init(const std::integral_constant&){} - BOOST_MATH_GPU_ENABLED static void do_init(const std::integral_constant&) + BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant&){} + BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant&) { boost::math::log1p(static_cast(0.25), Policy()); } @@ -286,7 +284,9 @@ struct log1p_initializer BOOST_MATH_STATIC const init initializer; BOOST_MATH_GPU_ENABLED static void force_instantiate() { + #ifndef BOOST_MATH_HAS_GPU_SUPPORT initializer.force_instantiate(); + #endif } }; @@ -309,7 +309,7 @@ BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type log1p(T x, c policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; - typedef std::integral_constant::type return -x * x / 2; boost::math::detail::log1p_series s(x); s(); - std::uintmax_t max_iter = policies::get_max_series_iterations(); + boost::math::uintmax_t max_iter = policies::get_max_series_iterations(); T result = boost::math::tools::sum_series(s, policies::get_epsilon(), max_iter); @@ -476,40 +476,6 @@ BOOST_MATH_GPU_ENABLED inline typename tools::promote_args::type log1pmx(T x) } // namespace math } // namespace boost -#else // Special handling for NVRTC platform - -namespace boost { -namespace math { - -template -BOOST_MATH_GPU_ENABLED auto log1p(T x) -{ - return ::log1p(x); -} - -template <> -BOOST_MATH_GPU_ENABLED auto log1p(float x) -{ - return ::log1pf(x); -} - -template -BOOST_MATH_GPU_ENABLED auto log1p(T x, const Policy&) -{ - return ::log1p(x); -} - -template -BOOST_MATH_GPU_ENABLED auto log1p(float x, const Policy&) -{ - return ::log1pf(x); -} - -} // namespace math -} // namespace boost - -#endif // BOOST_MATH_HAS_NVRTC - #ifdef _MSC_VER #pragma warning(pop) #endif diff --git a/include/boost/math/special_functions/math_fwd.hpp b/include/boost/math/special_functions/math_fwd.hpp index 21f51e5075..16ae3b61eb 100644 --- a/include/boost/math/special_functions/math_fwd.hpp +++ b/include/boost/math/special_functions/math_fwd.hpp @@ -512,28 +512,28 @@ namespace boost // gamma inverse. template - tools::promote_args_t gamma_p_inv(T1 a, T2 p); + BOOST_MATH_GPU_ENABLED tools::promote_args_t gamma_p_inv(T1 a, T2 p); template - tools::promote_args_t gamma_p_inva(T1 a, T2 p, const Policy&); + BOOST_MATH_GPU_ENABLED tools::promote_args_t gamma_p_inva(T1 a, T2 p, const Policy&); template - tools::promote_args_t gamma_p_inva(T1 a, T2 p); + BOOST_MATH_GPU_ENABLED tools::promote_args_t gamma_p_inva(T1 a, T2 p); template - tools::promote_args_t gamma_p_inv(T1 a, T2 p, const Policy&); + BOOST_MATH_GPU_ENABLED tools::promote_args_t gamma_p_inv(T1 a, T2 p, const Policy&); template - tools::promote_args_t gamma_q_inv(T1 a, T2 q); + BOOST_MATH_GPU_ENABLED tools::promote_args_t gamma_q_inv(T1 a, T2 q); template - tools::promote_args_t gamma_q_inv(T1 a, T2 q, const Policy&); + BOOST_MATH_GPU_ENABLED tools::promote_args_t gamma_q_inv(T1 a, T2 q, const Policy&); template - tools::promote_args_t gamma_q_inva(T1 a, T2 q); + BOOST_MATH_GPU_ENABLED tools::promote_args_t gamma_q_inva(T1 a, T2 q); template - tools::promote_args_t gamma_q_inva(T1 a, T2 q, const Policy&); + BOOST_MATH_GPU_ENABLED tools::promote_args_t gamma_q_inva(T1 a, T2 q, const Policy&); // digamma: template @@ -1447,16 +1447,16 @@ namespace boost BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t gamma_p_derivative(T1 a, T2 x){ return boost::math::gamma_p_derivative(a, x, Policy()); }\ \ template \ - inline boost::math::tools::promote_args_t gamma_p_inv(T1 a, T2 p){ return boost::math::gamma_p_inv(a, p, Policy()); }\ + BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t gamma_p_inv(T1 a, T2 p){ return boost::math::gamma_p_inv(a, p, Policy()); }\ \ template \ - inline boost::math::tools::promote_args_t gamma_p_inva(T1 a, T2 p){ return boost::math::gamma_p_inva(a, p, Policy()); }\ + BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t gamma_p_inva(T1 a, T2 p){ return boost::math::gamma_p_inva(a, p, Policy()); }\ \ template \ - inline boost::math::tools::promote_args_t gamma_q_inv(T1 a, T2 q){ return boost::math::gamma_q_inv(a, q, Policy()); }\ + BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t gamma_q_inv(T1 a, T2 q){ return boost::math::gamma_q_inv(a, q, Policy()); }\ \ template \ - inline boost::math::tools::promote_args_t gamma_q_inva(T1 a, T2 q){ return boost::math::gamma_q_inva(a, q, Policy()); }\ + BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t gamma_q_inva(T1 a, T2 q){ return boost::math::gamma_q_inva(a, q, Policy()); }\ \ template \ BOOST_MATH_GPU_ENABLED inline boost::math::tools::promote_args_t digamma(T x){ return boost::math::digamma(x, Policy()); }\ diff --git a/include/boost/math/tools/config.hpp b/include/boost/math/tools/config.hpp index 1f444c0048..a2af127630 100644 --- a/include/boost/math/tools/config.hpp +++ b/include/boost/math/tools/config.hpp @@ -677,6 +677,7 @@ namespace boost{ namespace math{ #include #include #include +#include # define BOOST_MATH_CUDA_ENABLED __host__ __device__ # define BOOST_MATH_HAS_GPU_SUPPORT @@ -802,6 +803,7 @@ BOOST_MATH_GPU_ENABLED constexpr T gpu_safe_max(const T& a, const T& b) { return #define BOOST_MATH_IF_CONSTEXPR if constexpr #define BOOST_MATH_IS_FLOAT(T) (boost::math::is_floating_point::value) #define BOOST_MATH_CONSTEXPR_TABLE_FUNCTION constexpr +#define BOOST_MATH_NO_EXCEPTIONS // This should be defined to nothing but since it is not specifically a math macro // we need to undef before proceeding diff --git a/include/boost/math/tools/fraction.hpp b/include/boost/math/tools/fraction.hpp index e5a31edf65..f36d024c40 100644 --- a/include/boost/math/tools/fraction.hpp +++ b/include/boost/math/tools/fraction.hpp @@ -18,9 +18,6 @@ #include #include #include -#include -#include -#include namespace boost{ namespace math{ namespace tools{ @@ -67,7 +64,7 @@ namespace detail template struct fraction_traits - : public std::conditional< + : public boost::math::conditional< is_pair::value, fraction_traits_pair, fraction_traits_simple>::type @@ -115,7 +112,7 @@ namespace detail { template BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type continued_fraction_b_impl(Gen& g, const U& factor, boost::math::uintmax_t& max_terms) noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits::result_type) - #ifndef BOOST_MATH_ENABLE_SYCL + #ifndef BOOST_MATH_HAS_GPU_SUPPORT // SYCL can not handle this condition so we only check float on that platform && noexcept(std::declval()()) #endif @@ -167,7 +164,7 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type template BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type continued_fraction_b(Gen& g, const U& factor, boost::math::uintmax_t& max_terms) noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits::result_type) - #ifndef BOOST_MATH_ENABLE_SYCL + #ifndef BOOST_MATH_HAS_GPU_SUPPORT && noexcept(std::declval()()) #endif ) @@ -178,19 +175,19 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type template BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type continued_fraction_b(Gen& g, const U& factor) noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits::result_type) - #ifndef BOOST_MATH_ENABLE_SYCL + #ifndef BOOST_MATH_HAS_GPU_SUPPORT && noexcept(std::declval()()) #endif ) { - boost::math::uintmax_t max_terms = (std::numeric_limits::max)(); + boost::math::uintmax_t max_terms = (boost::math::numeric_limits::max)(); return detail::continued_fraction_b_impl(g, factor, max_terms); } template BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type continued_fraction_b(Gen& g, int bits) noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits::result_type) - #ifndef BOOST_MATH_ENABLE_SYCL + #ifndef BOOST_MATH_HAS_GPU_SUPPORT && noexcept(std::declval()()) #endif ) @@ -201,14 +198,14 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type using result_type = typename traits::result_type; result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits); - boost::math::uintmax_t max_terms = (std::numeric_limits::max)(); + boost::math::uintmax_t max_terms = (boost::math::numeric_limits::max)(); return detail::continued_fraction_b_impl(g, factor, max_terms); } template BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type continued_fraction_b(Gen& g, int bits, boost::math::uintmax_t& max_terms) noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits::result_type) - #ifndef BOOST_MATH_ENABLE_SYCL + #ifndef BOOST_MATH_HAS_GPU_SUPPORT && noexcept(std::declval()()) #endif ) @@ -241,7 +238,7 @@ namespace detail { template BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type continued_fraction_a_impl(Gen& g, const U& factor, boost::math::uintmax_t& max_terms) noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits::result_type) - #ifndef BOOST_MATH_ENABLE_SYCL + #ifndef BOOST_MATH_HAS_GPU_SUPPORT && noexcept(std::declval()()) #endif ) @@ -294,7 +291,7 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type template BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type continued_fraction_a(Gen& g, const U& factor, boost::math::uintmax_t& max_terms) noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits::result_type) - #ifndef BOOST_MATH_ENABLE_SYCL + #ifndef BOOST_MATH_HAS_GPU_SUPPORT && noexcept(std::declval()()) #endif ) @@ -305,19 +302,19 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type template BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type continued_fraction_a(Gen& g, const U& factor) noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits::result_type) - #ifndef BOOST_MATH_ENABLE_SYCL + #ifndef BOOST_MATH_HAS_GPU_SUPPORT && noexcept(std::declval()()) #endif ) { - boost::math::uintmax_t max_iter = (std::numeric_limits::max)(); + boost::math::uintmax_t max_iter = (boost::math::numeric_limits::max)(); return detail::continued_fraction_a_impl(g, factor, max_iter); } template BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type continued_fraction_a(Gen& g, int bits) noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits::result_type) - #ifndef BOOST_MATH_ENABLE_SYCL + #ifndef BOOST_MATH_HAS_GPU_SUPPORT && noexcept(std::declval()()) #endif ) @@ -328,7 +325,7 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type typedef typename traits::result_type result_type; result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits); - boost::math::uintmax_t max_iter = (std::numeric_limits::max)(); + boost::math::uintmax_t max_iter = (boost::math::numeric_limits::max)(); return detail::continued_fraction_a_impl(g, factor, max_iter); } @@ -336,7 +333,7 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type template BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits::result_type continued_fraction_a(Gen& g, int bits, boost::math::uintmax_t& max_terms) noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits::result_type) - #ifndef BOOST_MATH_ENABLE_SYCL + #ifndef BOOST_MATH_HAS_GPU_SUPPORT && noexcept(std::declval()()) #endif ) diff --git a/include/boost/math/tools/precision.hpp b/include/boost/math/tools/precision.hpp index 1e4bf58ac0..662657732c 100644 --- a/include/boost/math/tools/precision.hpp +++ b/include/boost/math/tools/precision.hpp @@ -290,7 +290,7 @@ template BOOST_MATH_GPU_ENABLED inline T root_epsilon_imp(const T*, const Tag&) { BOOST_MATH_STD_USING - static const T r_eps = sqrt(tools::epsilon()); + BOOST_MATH_STATIC_LOCAL_VARIABLE const T r_eps = sqrt(tools::epsilon()); return r_eps; } diff --git a/include/boost/math/tools/series.hpp b/include/boost/math/tools/series.hpp index 50f2828bb4..4617ea3df7 100644 --- a/include/boost/math/tools/series.hpp +++ b/include/boost/math/tools/series.hpp @@ -14,6 +14,7 @@ #include #include #include +#include namespace boost{ namespace math{ namespace tools{ diff --git a/include/boost/math/tools/toms748_solve.hpp b/include/boost/math/tools/toms748_solve.hpp index ea93713224..dee2346853 100644 --- a/include/boost/math/tools/toms748_solve.hpp +++ b/include/boost/math/tools/toms748_solve.hpp @@ -1,4 +1,5 @@ // (C) Copyright John Maddock 2006. +// (C) Copyright Matt Borland 2024. // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) @@ -10,13 +11,13 @@ #pragma once #endif +#include #include +#include +#include +#include #include -#include #include -#include -#include -#include #ifdef BOOST_MATH_LOG_ROOT_ITERATIONS # define BOOST_MATH_LOGGER_INCLUDE @@ -32,29 +33,36 @@ template class eps_tolerance { public: - eps_tolerance() : eps(4 * tools::epsilon()) + BOOST_MATH_GPU_ENABLED eps_tolerance() : eps(4 * tools::epsilon()) { } - eps_tolerance(unsigned bits) + BOOST_MATH_GPU_ENABLED eps_tolerance(unsigned bits) { BOOST_MATH_STD_USING - eps = (std::max)(T(ldexp(1.0F, 1-bits)), T(4 * tools::epsilon())); + eps = BOOST_MATH_GPU_SAFE_MAX(T(ldexp(1.0F, 1-bits)), T(4 * tools::epsilon())); } - bool operator()(const T& a, const T& b) + BOOST_MATH_GPU_ENABLED bool operator()(const T& a, const T& b) { BOOST_MATH_STD_USING - return fabs(a - b) <= (eps * (std::min)(fabs(a), fabs(b))); + return fabs(a - b) <= (eps * BOOST_MATH_GPU_SAFE_MIN(fabs(a), fabs(b))); } private: T eps; }; +// CUDA warns about __host__ __device__ marker on defaulted constructor +// but the warning is benign +#ifdef BOOST_MATH_ENABLE_CUDA +# pragma nv_diag_suppress 20012 +#endif + struct equal_floor { - equal_floor()= default; + BOOST_MATH_GPU_ENABLED equal_floor() = default; + template - bool operator()(const T& a, const T& b) + BOOST_MATH_GPU_ENABLED bool operator()(const T& a, const T& b) { BOOST_MATH_STD_USING return (floor(a) == floor(b)) || (fabs((b-a)/b) < boost::math::tools::epsilon() * 2); @@ -63,9 +71,10 @@ struct equal_floor struct equal_ceil { - equal_ceil()= default; + BOOST_MATH_GPU_ENABLED equal_ceil() = default; + template - bool operator()(const T& a, const T& b) + BOOST_MATH_GPU_ENABLED bool operator()(const T& a, const T& b) { BOOST_MATH_STD_USING return (ceil(a) == ceil(b)) || (fabs((b - a) / b) < boost::math::tools::epsilon() * 2); @@ -74,19 +83,24 @@ struct equal_ceil struct equal_nearest_integer { - equal_nearest_integer()= default; + BOOST_MATH_GPU_ENABLED equal_nearest_integer() = default; + template - bool operator()(const T& a, const T& b) + BOOST_MATH_GPU_ENABLED bool operator()(const T& a, const T& b) { BOOST_MATH_STD_USING return (floor(a + 0.5f) == floor(b + 0.5f)) || (fabs((b - a) / b) < boost::math::tools::epsilon() * 2); } }; +#ifdef BOOST_MATH_ENABLE_CUDA +# pragma nv_diag_default 20012 +#endif + namespace detail{ template -void bracket(F f, T& a, T& b, T c, T& fa, T& fb, T& d, T& fd) +BOOST_MATH_GPU_ENABLED void bracket(F f, T& a, T& b, T c, T& fa, T& fb, T& d, T& fd) { // // Given a point c inside the existing enclosing interval @@ -150,7 +164,7 @@ void bracket(F f, T& a, T& b, T c, T& fa, T& fb, T& d, T& fd) } template -inline T safe_div(T num, T denom, T r) +BOOST_MATH_GPU_ENABLED inline T safe_div(T num, T denom, T r) { // // return num / denom without overflow, @@ -167,7 +181,7 @@ inline T safe_div(T num, T denom, T r) } template -inline T secant_interpolate(const T& a, const T& b, const T& fa, const T& fb) +BOOST_MATH_GPU_ENABLED inline T secant_interpolate(const T& a, const T& b, const T& fa, const T& fb) { // // Performs standard secant interpolation of [a,b] given @@ -188,9 +202,9 @@ inline T secant_interpolate(const T& a, const T& b, const T& fa, const T& fb) } template -T quadratic_interpolate(const T& a, const T& b, T const& d, - const T& fa, const T& fb, T const& fd, - unsigned count) +BOOST_MATH_GPU_ENABLED T quadratic_interpolate(const T& a, const T& b, T const& d, + const T& fa, const T& fb, T const& fd, + unsigned count) { // // Performs quadratic interpolation to determine the next point, @@ -244,9 +258,9 @@ T quadratic_interpolate(const T& a, const T& b, T const& d, } template -T cubic_interpolate(const T& a, const T& b, const T& d, - const T& e, const T& fa, const T& fb, - const T& fd, const T& fe) +BOOST_MATH_GPU_ENABLED T cubic_interpolate(const T& a, const T& b, const T& d, + const T& e, const T& fa, const T& fb, + const T& fd, const T& fe) { // // Uses inverse cubic interpolation of f(x) at points @@ -293,7 +307,7 @@ T cubic_interpolate(const T& a, const T& b, const T& d, } // namespace detail template -std::pair toms748_solve(F f, const T& ax, const T& bx, const T& fax, const T& fbx, Tol tol, std::uintmax_t& max_iter, const Policy& pol) +BOOST_MATH_GPU_ENABLED boost::math::pair toms748_solve(F f, const T& ax, const T& bx, const T& fax, const T& fbx, Tol tol, boost::math::uintmax_t& max_iter, const Policy& pol) { // // Main entry point and logic for Toms Algorithm 748 @@ -301,15 +315,15 @@ std::pair toms748_solve(F f, const T& ax, const T& bx, const T& fax, const // BOOST_MATH_STD_USING // For ADL of std math functions - static const char* function = "boost::math::tools::toms748_solve<%1%>"; + constexpr auto function = "boost::math::tools::toms748_solve<%1%>"; // // Sanity check - are we allowed to iterate at all? // if (max_iter == 0) - return std::make_pair(ax, bx); + return boost::math::make_pair(ax, bx); - std::uintmax_t count = max_iter; + boost::math::uintmax_t count = max_iter; T a, b, fa, fb, c, u, fu, a0, b0, d, fd, e, fe; static const T mu = 0.5f; @@ -330,7 +344,7 @@ std::pair toms748_solve(F f, const T& ax, const T& bx, const T& fax, const b = a; else if(fb == 0) a = b; - return std::make_pair(a, b); + return boost::math::make_pair(a, b); } if(boost::math::sign(fa) * boost::math::sign(fb) > 0) @@ -472,37 +486,37 @@ std::pair toms748_solve(F f, const T& ax, const T& bx, const T& fax, const a = b; } BOOST_MATH_LOG_COUNT(max_iter) - return std::make_pair(a, b); + return boost::math::make_pair(a, b); } template -inline std::pair toms748_solve(F f, const T& ax, const T& bx, const T& fax, const T& fbx, Tol tol, std::uintmax_t& max_iter) +BOOST_MATH_GPU_ENABLED inline boost::math::pair toms748_solve(F f, const T& ax, const T& bx, const T& fax, const T& fbx, Tol tol, boost::math::uintmax_t& max_iter) { return toms748_solve(f, ax, bx, fax, fbx, tol, max_iter, policies::policy<>()); } template -inline std::pair toms748_solve(F f, const T& ax, const T& bx, Tol tol, std::uintmax_t& max_iter, const Policy& pol) +BOOST_MATH_GPU_ENABLED inline boost::math::pair toms748_solve(F f, const T& ax, const T& bx, Tol tol, boost::math::uintmax_t& max_iter, const Policy& pol) { if (max_iter <= 2) - return std::make_pair(ax, bx); + return boost::math::make_pair(ax, bx); max_iter -= 2; - std::pair r = toms748_solve(f, ax, bx, f(ax), f(bx), tol, max_iter, pol); + boost::math::pair r = toms748_solve(f, ax, bx, f(ax), f(bx), tol, max_iter, pol); max_iter += 2; return r; } template -inline std::pair toms748_solve(F f, const T& ax, const T& bx, Tol tol, std::uintmax_t& max_iter) +BOOST_MATH_GPU_ENABLED inline boost::math::pair toms748_solve(F f, const T& ax, const T& bx, Tol tol, boost::math::uintmax_t& max_iter) { return toms748_solve(f, ax, bx, tol, max_iter, policies::policy<>()); } template -std::pair bracket_and_solve_root(F f, const T& guess, T factor, bool rising, Tol tol, std::uintmax_t& max_iter, const Policy& pol) +BOOST_MATH_GPU_ENABLED boost::math::pair bracket_and_solve_root(F f, const T& guess, T factor, bool rising, Tol tol, boost::math::uintmax_t& max_iter, const Policy& pol) { BOOST_MATH_STD_USING - static const char* function = "boost::math::tools::bracket_and_solve_root<%1%>"; + constexpr auto function = "boost::math::tools::bracket_and_solve_root<%1%>"; // // Set up initial brackets: // @@ -513,7 +527,7 @@ std::pair bracket_and_solve_root(F f, const T& guess, T factor, bool risin // // Set up invocation count: // - std::uintmax_t count = max_iter - 1; + boost::math::uintmax_t count = max_iter - 1; int step = 32; @@ -563,7 +577,7 @@ std::pair bracket_and_solve_root(F f, const T& guess, T factor, bool risin // Escape route just in case the answer is zero! max_iter -= count; max_iter += 1; - return a > 0 ? std::make_pair(T(0), T(a)) : std::make_pair(T(a), T(0)); + return a > 0 ? boost::math::make_pair(T(0), T(a)) : boost::math::make_pair(T(a), T(0)); } if(count == 0) return boost::math::detail::pair_from_single(policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", a, pol)); @@ -592,7 +606,7 @@ std::pair bracket_and_solve_root(F f, const T& guess, T factor, bool risin } max_iter -= count; max_iter += 1; - std::pair r = toms748_solve( + boost::math::pair r = toms748_solve( f, (a < 0 ? b : a), (a < 0 ? a : b), @@ -608,7 +622,7 @@ std::pair bracket_and_solve_root(F f, const T& guess, T factor, bool risin } template -inline std::pair bracket_and_solve_root(F f, const T& guess, const T& factor, bool rising, Tol tol, std::uintmax_t& max_iter) +BOOST_MATH_GPU_ENABLED inline boost::math::pair bracket_and_solve_root(F f, const T& guess, const T& factor, bool rising, Tol tol, boost::math::uintmax_t& max_iter) { return bracket_and_solve_root(f, guess, factor, rising, tol, max_iter, policies::policy<>()); } diff --git a/include/boost/math/tools/tuple.hpp b/include/boost/math/tools/tuple.hpp index 6a949d5e5a..82d23b8d77 100644 --- a/include/boost/math/tools/tuple.hpp +++ b/include/boost/math/tools/tuple.hpp @@ -13,6 +13,7 @@ #include #include +#include namespace boost { namespace math { diff --git a/test/cuda_jamfile b/test/cuda_jamfile index 64fac76fba..1fb55da197 100644 --- a/test/cuda_jamfile +++ b/test/cuda_jamfile @@ -34,6 +34,13 @@ run test_cauchy_quan_float.cu ; run test_cauchy_range_support_double.cu ; run test_cauchy_range_support_float.cu ; +run test_chi_squared_cdf_double.cu ; +run test_chi_squared_cdf_float.cu ; +run test_chi_squared_pdf_double.cu ; +run test_chi_squared_pdf_float.cu ; +run test_chi_squared_quan_double.cu ; +run test_chi_squared_quan_float.cu ; + run test_exponential_cdf_double.cu ; run test_exponential_cdf_float.cu ; run test_exponential_pdf_double.cu ; @@ -90,6 +97,13 @@ run test_saspoint5_pdf_float.cu ; run test_saspoint5_quan_double.cu ; run test_saspoint5_quan_float.cu ; +run test_weibull_cdf_double.cu ; +run test_weibull_cdf_float.cu ; +run test_weibull_pdf_double.cu ; +run test_weibull_pdf_float.cu ; +run test_weibull_quan_double.cu ; +run test_weibull_quan_float.cu ; + # Special Functions run test_beta_double.cu ; run test_beta_float.cu ; @@ -153,6 +167,12 @@ run test_lgamma_double.cu ; run test_lgamma_float.cu ; run test_tgamma_double.cu ; run test_tgamma_float.cu ; +run test_tgamma_ratio_double.cu ; +run test_tgamma_ratio_float.cu ; +run test_gamma_p_derivative_double.cu ; +run test_gamma_p_derivative_float.cu ; +run test_gamma_p_inv_double.cu ; +run test_gamma_p_inv_float.cu ; run test_log1p_double.cu ; run test_log1p_float.cu ; diff --git a/test/nvrtc_jamfile b/test/nvrtc_jamfile index de235822ef..4a37960a51 100644 --- a/test/nvrtc_jamfile +++ b/test/nvrtc_jamfile @@ -31,6 +31,13 @@ run test_cauchy_pdf_nvrtc_float.cpp ; run test_cauchy_quan_nvrtc_double.cpp ; run test_cauchy_quan_nvrtc_float.cpp ; +run test_chi_squared_cdf_nvrtc_double.cpp ; +run test_chi_squared_cdf_nvrtc_float.cpp ; +run test_chi_squared_pdf_nvrtc_double.cpp ; +run test_chi_squared_pdf_nvrtc_float.cpp ; +run test_chi_squared_quan_nvrtc_double.cpp ; +run test_chi_squared_quan_nvrtc_float.cpp ; + run test_exponential_cdf_nvrtc_double.cpp ; run test_exponential_cdf_nvrtc_float.cpp ; run test_exponential_pdf_nvrtc_double.cpp ; @@ -87,6 +94,13 @@ run test_saspoint5_pdf_nvrtc_float.cpp ; run test_saspoint5_quan_nvrtc_double.cpp ; run test_saspoint5_quan_nvrtc_float.cpp ; +run test_weibull_cdf_nvrtc_double.cpp ; +run test_weibull_cdf_nvrtc_float.cpp ; +run test_weibull_pdf_nvrtc_double.cpp ; +run test_weibull_pdf_nvrtc_float.cpp ; +run test_weibull_quan_nvrtc_double.cpp ; +run test_weibull_quan_nvrtc_float.cpp ; + # Special Functions run test_beta_nvrtc_double.cpp ; run test_beta_nvrtc_float.cpp ; @@ -148,6 +162,12 @@ run test_fpclassify_nvrtc_float.cpp ; run test_gamma_nvrtc_double.cpp ; run test_gamma_nvrtc_float.cpp ; +run test_gamma_p_derivative_nvrtc_double.cpp ; +run test_gamma_p_derivative_nvrtc_float.cpp ; +run test_gamma_p_inv_nvrtc_double.cpp ; +run test_gamma_p_inv_nvrtc_float.cpp ; +run test_tgamma_ratio_nvrtc_double.cpp ; +run test_tgamma_ratio_nvrtc_float.cpp ; run test_log1p_nvrtc_double.cpp ; run test_log1p_nvrtc_float.cpp ; diff --git a/test/sycl_jamfile b/test/sycl_jamfile index 97c48474ca..2fd5954ae1 100644 --- a/test/sycl_jamfile +++ b/test/sycl_jamfile @@ -13,6 +13,7 @@ project : requirements run test_arcsine.cpp ; run test_bernoulli.cpp ; run test_cauchy.cpp ; +run test_chi_squared.cpp ; run test_exponential_dist.cpp ; run test_extreme_value.cpp ; run test_holtsmark.cpp ; @@ -21,20 +22,35 @@ run test_laplace.cpp ; run test_logistic_dist.cpp ; run test_mapairy.cpp ; run test_saspoint5.cpp ; +run test_weibull.cpp ; # Special Functions run pow_test.cpp ; + run test_beta_simple.cpp ; + run test_bessel_i.cpp ; run test_bessel_j.cpp ; run test_bessel_k.cpp ; run test_bessel_y.cpp ; + run test_cbrt.cpp ; + run test_sign.cpp ; + run test_round.cpp ; + run test_expm1_simple.cpp ; + run test_log1p_simple.cpp ; + run test_digamma_simple.cpp ; + run test_trigamma.cpp ; + run test_erf.cpp ; + run test_gamma.cpp ; +run test_igamma.cpp ; +run test_igamma_inv.cpp ; +run test_igamma_inva.cpp ; diff --git a/test/test_bessel_i.cpp b/test/test_bessel_i.cpp index 70aac91e45..09487ddf1b 100644 --- a/test/test_bessel_i.cpp +++ b/test/test_bessel_i.cpp @@ -10,6 +10,14 @@ #include #endif +#ifdef __clang__ +# pragma clang diagnostic push +# pragma clang diagnostic ignored "-Wliteral-range" +#elif defined(__GNUC__) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wliteral-range" +#endif + #include "test_bessel_i.hpp" // diff --git a/test/test_bessel_j.cpp b/test/test_bessel_j.cpp index 516e34c29a..31a64bc579 100644 --- a/test/test_bessel_j.cpp +++ b/test/test_bessel_j.cpp @@ -10,6 +10,14 @@ #include #endif +#ifdef __clang__ +# pragma clang diagnostic push +# pragma clang diagnostic ignored "-Wliteral-range" +#elif defined(__GNUC__) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wliteral-range" +#endif + #include "test_bessel_j.hpp" // diff --git a/test/test_bessel_k.cpp b/test/test_bessel_k.cpp index d4ab7721f4..84ba0830f2 100644 --- a/test/test_bessel_k.cpp +++ b/test/test_bessel_k.cpp @@ -17,6 +17,14 @@ // Constants are too big for float case, but this doesn't matter for test. #endif +#ifdef __clang__ +# pragma clang diagnostic push +# pragma clang diagnostic ignored "-Wliteral-range" +#elif defined(__GNUC__) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wliteral-range" +#endif + #include "test_bessel_k.hpp" // diff --git a/test/test_bessel_y.cpp b/test/test_bessel_y.cpp index 0bbefba55f..232a903963 100644 --- a/test/test_bessel_y.cpp +++ b/test/test_bessel_y.cpp @@ -10,6 +10,14 @@ #include #endif +#ifdef __clang__ +# pragma clang diagnostic push +# pragma clang diagnostic ignored "-Wliteral-range" +#elif defined(__GNUC__) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wliteral-range" +#endif + #include "test_bessel_y.hpp" // diff --git a/test/test_chi_squared.cpp b/test/test_chi_squared.cpp index cc7747a6c0..bfd4b5f3a2 100644 --- a/test/test_chi_squared.cpp +++ b/test/test_chi_squared.cpp @@ -16,9 +16,13 @@ # pragma warning(disable: 4127) // conditional expression is constant #endif -#include // for real_concept +#include +#include "../include_private/boost/math/tools/test.hpp" + +#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS #include // for real_concept using ::boost::math::concepts::real_concept; +#endif #include // for chi_squared_distribution #include // for chi_squared_distribution diff --git a/test/test_chi_squared_cdf_double.cu b/test/test_chi_squared_cdf_double.cu new file mode 100644 index 0000000000..1b0c34ce6b --- /dev/null +++ b/test/test_chi_squared_cdf_double.cu @@ -0,0 +1,109 @@ +// Copyright John Maddock 2016. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef double float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = cdf(boost::math::chi_squared_distribution(1), in1[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + try{ + + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + boost::random::mt19937 gen; + boost::random::uniform_real_distribution dist; + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 512; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + cuda_test<<>>(input_vector1.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + if (err != cudaSuccess) + { + std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(cdf(boost::math::chi_squared_distribution(1), input_vector1[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 100.0) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + std::cerr << "Error rate was: " << boost::math::epsilon_difference(output_vector[i], results[i]) << "eps" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + } + return 0; +} \ No newline at end of file diff --git a/test/test_chi_squared_cdf_float.cu b/test/test_chi_squared_cdf_float.cu new file mode 100644 index 0000000000..8ca99ed2e0 --- /dev/null +++ b/test/test_chi_squared_cdf_float.cu @@ -0,0 +1,109 @@ +// Copyright John Maddock 2016. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef float float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = cdf(boost::math::chi_squared_distribution(1), in1[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + try{ + + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + boost::random::mt19937 gen; + boost::random::uniform_real_distribution dist; + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 512; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + cuda_test<<>>(input_vector1.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + if (err != cudaSuccess) + { + std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(cdf(boost::math::chi_squared_distribution(1), input_vector1[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 100.0) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + std::cerr << "Error rate was: " << boost::math::epsilon_difference(output_vector[i], results[i]) << "eps" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + } + return 0; +} \ No newline at end of file diff --git a/test/test_chi_squared_cdf_nvrtc_double.cpp b/test/test_chi_squared_cdf_nvrtc_double.cpp new file mode 100644 index 0000000000..0ad459fa67 --- /dev/null +++ b/test/test_chi_squared_cdf_nvrtc_double.cpp @@ -0,0 +1,191 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +typedef double float_type; + +const char* cuda_kernel = R"( +typedef double float_type; +#include +#include +extern "C" __global__ +void test_chi_squared_kernel(const float_type *in1, const float_type*, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = cdf(boost::math::chi_squared_distribution(1), in1[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_chi_squared_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_chi_squared_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_chi_squared_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + auto res = cdf(boost::math::chi_squared_distribution(1), h_in1[i]); + + if (boost::math::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_chi_squared_cdf_nvrtc_float.cpp b/test/test_chi_squared_cdf_nvrtc_float.cpp new file mode 100644 index 0000000000..1b26c5d6f2 --- /dev/null +++ b/test/test_chi_squared_cdf_nvrtc_float.cpp @@ -0,0 +1,191 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +typedef float float_type; + +const char* cuda_kernel = R"( +typedef float float_type; +#include +#include +extern "C" __global__ +void test_chi_squared_kernel(const float_type *in1, const float_type*, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = cdf(boost::math::chi_squared_distribution(1), in1[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_chi_squared_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_chi_squared_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_chi_squared_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + auto res = cdf(boost::math::chi_squared_distribution(1), h_in1[i]); + + if (boost::math::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_chi_squared_pdf_double.cu b/test/test_chi_squared_pdf_double.cu new file mode 100644 index 0000000000..ed45246d3c --- /dev/null +++ b/test/test_chi_squared_pdf_double.cu @@ -0,0 +1,109 @@ +// Copyright John Maddock 2016. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef double float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = pdf(boost::math::chi_squared_distribution(1), in1[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + try{ + + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + boost::random::mt19937 gen; + boost::random::uniform_real_distribution dist; + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 512; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + cuda_test<<>>(input_vector1.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + if (err != cudaSuccess) + { + std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(pdf(boost::math::chi_squared_distribution(1), input_vector1[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 100.0) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + std::cerr << "Error rate was: " << boost::math::epsilon_difference(output_vector[i], results[i]) << "eps" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + } + return 0; +} \ No newline at end of file diff --git a/test/test_chi_squared_pdf_float.cu b/test/test_chi_squared_pdf_float.cu new file mode 100644 index 0000000000..5a0f97db93 --- /dev/null +++ b/test/test_chi_squared_pdf_float.cu @@ -0,0 +1,109 @@ +// Copyright John Maddock 2016. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef float float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = pdf(boost::math::chi_squared_distribution(1), in1[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + try{ + + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + boost::random::mt19937 gen; + boost::random::uniform_real_distribution dist; + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 512; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + cuda_test<<>>(input_vector1.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + if (err != cudaSuccess) + { + std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(pdf(boost::math::chi_squared_distribution(1), input_vector1[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 100.0) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + std::cerr << "Error rate was: " << boost::math::epsilon_difference(output_vector[i], results[i]) << "eps" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + } + return 0; +} \ No newline at end of file diff --git a/test/test_chi_squared_pdf_nvrtc_double.cpp b/test/test_chi_squared_pdf_nvrtc_double.cpp new file mode 100644 index 0000000000..18d14a4b0e --- /dev/null +++ b/test/test_chi_squared_pdf_nvrtc_double.cpp @@ -0,0 +1,191 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +typedef double float_type; + +const char* cuda_kernel = R"( +typedef double float_type; +#include +#include +extern "C" __global__ +void test_chi_squared_kernel(const float_type *in1, const float_type*, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = pdf(boost::math::chi_squared_distribution(1), in1[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_chi_squared_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_chi_squared_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_chi_squared_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + auto res = pdf(boost::math::chi_squared_distribution(1), h_in1[i]); + + if (boost::math::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_chi_squared_pdf_nvrtc_float.cpp b/test/test_chi_squared_pdf_nvrtc_float.cpp new file mode 100644 index 0000000000..754cbf7fba --- /dev/null +++ b/test/test_chi_squared_pdf_nvrtc_float.cpp @@ -0,0 +1,191 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +typedef float float_type; + +const char* cuda_kernel = R"( +typedef float float_type; +#include +#include +extern "C" __global__ +void test_chi_squared_kernel(const float_type *in1, const float_type*, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = pdf(boost::math::chi_squared_distribution(1), in1[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_chi_squared_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_chi_squared_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_chi_squared_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + auto res = pdf(boost::math::chi_squared_distribution(1), h_in1[i]); + + if (boost::math::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_chi_squared_quan_double.cu b/test/test_chi_squared_quan_double.cu new file mode 100644 index 0000000000..3b7dad972b --- /dev/null +++ b/test/test_chi_squared_quan_double.cu @@ -0,0 +1,109 @@ +// Copyright John Maddock 2016. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef double float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = quantile(boost::math::chi_squared_distribution(1), in1[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + try{ + + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + boost::random::mt19937 gen; + boost::random::uniform_real_distribution dist; + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 512; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + cuda_test<<>>(input_vector1.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + if (err != cudaSuccess) + { + std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(quantile(boost::math::chi_squared_distribution(1), input_vector1[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 100.0) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + std::cerr << "Error rate was: " << boost::math::epsilon_difference(output_vector[i], results[i]) << "eps" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + } + return 0; +} \ No newline at end of file diff --git a/test/test_chi_squared_quan_float.cu b/test/test_chi_squared_quan_float.cu new file mode 100644 index 0000000000..3e779a090f --- /dev/null +++ b/test/test_chi_squared_quan_float.cu @@ -0,0 +1,109 @@ +// Copyright John Maddock 2016. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef float float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = quantile(boost::math::chi_squared_distribution(1), in1[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + try{ + + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + boost::random::mt19937 gen; + boost::random::uniform_real_distribution dist; + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 512; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + cuda_test<<>>(input_vector1.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + if (err != cudaSuccess) + { + std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(quantile(boost::math::chi_squared_distribution(1), input_vector1[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 100.0) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + std::cerr << "Error rate was: " << boost::math::epsilon_difference(output_vector[i], results[i]) << "eps" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + } + return 0; +} \ No newline at end of file diff --git a/test/test_chi_squared_quan_nvrtc_double.cpp b/test/test_chi_squared_quan_nvrtc_double.cpp new file mode 100644 index 0000000000..69b15b6cfd --- /dev/null +++ b/test/test_chi_squared_quan_nvrtc_double.cpp @@ -0,0 +1,191 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +typedef double float_type; + +const char* cuda_kernel = R"( +typedef double float_type; +#include +#include +extern "C" __global__ +void test_chi_squared_kernel(const float_type *in1, const float_type*, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = quantile(boost::math::chi_squared_distribution(1), in1[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_chi_squared_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_chi_squared_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_chi_squared_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + auto res = quantile(boost::math::chi_squared_distribution(1), h_in1[i]); + + if (boost::math::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_chi_squared_quan_nvrtc_float.cpp b/test/test_chi_squared_quan_nvrtc_float.cpp new file mode 100644 index 0000000000..d6e1b2a9b5 --- /dev/null +++ b/test/test_chi_squared_quan_nvrtc_float.cpp @@ -0,0 +1,191 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +typedef float float_type; + +const char* cuda_kernel = R"( +typedef float float_type; +#include +#include +extern "C" __global__ +void test_chi_squared_kernel(const float_type *in1, const float_type*, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = quantile(boost::math::chi_squared_distribution(1), in1[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_chi_squared_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_chi_squared_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_chi_squared_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + auto res = quantile(boost::math::chi_squared_distribution(1), h_in1[i]); + + if (boost::math::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_gamma_nvrtc_double.cpp b/test/test_gamma_nvrtc_double.cpp index 9302d1190e..9fe2933720 100644 --- a/test/test_gamma_nvrtc_double.cpp +++ b/test/test_gamma_nvrtc_double.cpp @@ -85,9 +85,9 @@ int main() nvrtcAddNameExpression(prog, "test_gamma_kernel"); #ifdef BOOST_MATH_NVRTC_CI_RUN - const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/"}; + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; #else - const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/"}; + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; #endif // Compile the program diff --git a/test/test_gamma_nvrtc_float.cpp b/test/test_gamma_nvrtc_float.cpp index 6e9ba72586..5d34b130ad 100644 --- a/test/test_gamma_nvrtc_float.cpp +++ b/test/test_gamma_nvrtc_float.cpp @@ -87,7 +87,7 @@ int main() #ifdef BOOST_MATH_NVRTC_CI_RUN const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/math/boost-root/libs/math/include/"}; #else - const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/"}; + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; #endif // Compile the program diff --git a/test/test_gamma_p_derivative_double.cu b/test/test_gamma_p_derivative_double.cu new file mode 100644 index 0000000000..566bc1657f --- /dev/null +++ b/test/test_gamma_p_derivative_double.cu @@ -0,0 +1,104 @@ + +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef double float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, const float_type *in2, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = boost::math::gamma_p_derivative(in1[i], in2[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed input vector B + cuda_managed_ptr input_vector2(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = rand()/(float_type)RAND_MAX; + input_vector2[i] = rand()/(float_type)RAND_MAX; + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 256; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + + cuda_test<<>>(input_vector1.get(), input_vector2.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + + std::cout << "CUDA kernal done in: " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + + if (err != cudaSuccess) + { + std::cerr << "Failed to launch CUDA kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(boost::math::gamma_p_derivative(input_vector1[i], input_vector2[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 10) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED, normal calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + + return 0; +} diff --git a/test/test_gamma_p_derivative_float.cu b/test/test_gamma_p_derivative_float.cu new file mode 100644 index 0000000000..f9fd52a50c --- /dev/null +++ b/test/test_gamma_p_derivative_float.cu @@ -0,0 +1,104 @@ + +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef float float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, const float_type *in2, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = boost::math::gamma_p_derivative(in1[i], in2[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed input vector B + cuda_managed_ptr input_vector2(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = rand()/(float_type)RAND_MAX; + input_vector2[i] = rand()/(float_type)RAND_MAX; + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 256; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + + cuda_test<<>>(input_vector1.get(), input_vector2.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + + std::cout << "CUDA kernal done in: " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + + if (err != cudaSuccess) + { + std::cerr << "Failed to launch CUDA kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(boost::math::gamma_p_derivative(input_vector1[i], input_vector2[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 10) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED, normal calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + + return 0; +} diff --git a/test/test_gamma_p_derivative_nvrtc_double.cpp b/test/test_gamma_p_derivative_nvrtc_double.cpp new file mode 100644 index 0000000000..53a752c2df --- /dev/null +++ b/test/test_gamma_p_derivative_nvrtc_double.cpp @@ -0,0 +1,190 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +typedef double float_type; + +const char* cuda_kernel = R"( +typedef double float_type; +#include +#include +extern "C" __global__ +void test_gamma_p_derivative_kernel(const float_type *in1, const float_type *in2, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = boost::math::gamma_p_derivative(in1[i], in2[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_gamma_p_derivative_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_gamma_p_derivative_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_gamma_p_derivative_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1000.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + const auto res = boost::math::gamma_p_derivative(h_in1[i], h_in2[i]); + + if (std::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_gamma_p_derivative_nvrtc_float.cpp b/test/test_gamma_p_derivative_nvrtc_float.cpp new file mode 100644 index 0000000000..da9c50855b --- /dev/null +++ b/test/test_gamma_p_derivative_nvrtc_float.cpp @@ -0,0 +1,190 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +typedef float float_type; + +const char* cuda_kernel = R"( +typedef float float_type; +#include +#include +extern "C" __global__ +void test_gamma_p_derivative_kernel(const float_type *in1, const float_type *in2, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = boost::math::gamma_p_derivative(in1[i], in2[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_gamma_p_derivative_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_gamma_p_derivative_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_gamma_p_derivative_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1000.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + const auto res = boost::math::gamma_p_derivative(h_in1[i], h_in2[i]); + + if (std::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_gamma_p_inv_double.cu b/test/test_gamma_p_inv_double.cu new file mode 100644 index 0000000000..4392f37d38 --- /dev/null +++ b/test/test_gamma_p_inv_double.cu @@ -0,0 +1,108 @@ + +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef double float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, const float_type *in2, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = boost::math::gamma_p_inv(in1[i], in2[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed input vector B + cuda_managed_ptr input_vector2(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + // Initialize the input vectors + std::mt19937_64 gen(42); + std::uniform_real_distribution dist(0, 1); + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + input_vector2[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 256; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + + cuda_test<<>>(input_vector1.get(), input_vector2.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + + std::cout << "CUDA kernal done in: " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + + if (err != cudaSuccess) + { + std::cerr << "Failed to launch CUDA kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(boost::math::gamma_p_inv(input_vector1[i], input_vector2[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 1000) + { + std::cerr << "Result verification failed at element " << i << "!\n" + << "Error found was: " << boost::math::epsilon_difference(output_vector[i], results[i]) << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED, normal calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + + return 0; +} diff --git a/test/test_gamma_p_inv_float.cu b/test/test_gamma_p_inv_float.cu new file mode 100644 index 0000000000..70033686c1 --- /dev/null +++ b/test/test_gamma_p_inv_float.cu @@ -0,0 +1,107 @@ + +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef float float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, const float_type *in2, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = boost::math::gamma_p_inv(in1[i], in2[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed input vector B + cuda_managed_ptr input_vector2(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + // Initialize the input vectors + std::mt19937_64 gen(42); + std::uniform_real_distribution dist(0, 1); + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + input_vector2[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 256; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + + cuda_test<<>>(input_vector1.get(), input_vector2.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + + std::cout << "CUDA kernal done in: " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + + if (err != cudaSuccess) + { + std::cerr << "Failed to launch CUDA kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(boost::math::gamma_p_inv(input_vector1[i], input_vector2[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 100) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED, normal calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + + return 0; +} diff --git a/test/test_gamma_p_inv_nvrtc_double.cpp b/test/test_gamma_p_inv_nvrtc_double.cpp new file mode 100644 index 0000000000..d270dbf901 --- /dev/null +++ b/test/test_gamma_p_inv_nvrtc_double.cpp @@ -0,0 +1,190 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +typedef double float_type; + +const char* cuda_kernel = R"( +typedef double float_type; +#include +#include +extern "C" __global__ +void test_gamma_p_inv_kernel(const float_type *in1, const float_type *in2, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = boost::math::gamma_p_inv(in1[i], in2[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_gamma_p_inv_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_gamma_p_inv_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_gamma_p_inv_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + const auto res = boost::math::gamma_p_inv(h_in1[i], h_in2[i]); + + if (std::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_gamma_p_inv_nvrtc_float.cpp b/test/test_gamma_p_inv_nvrtc_float.cpp new file mode 100644 index 0000000000..7c844eb682 --- /dev/null +++ b/test/test_gamma_p_inv_nvrtc_float.cpp @@ -0,0 +1,190 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +typedef float float_type; + +const char* cuda_kernel = R"( +typedef float float_type; +#include +#include +extern "C" __global__ +void test_gamma_p_inv_kernel(const float_type *in1, const float_type *in2, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = boost::math::gamma_p_inv(in1[i], in2[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_gamma_p_inv_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_gamma_p_inv_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_gamma_p_inv_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + const auto res = boost::math::gamma_p_inv(h_in1[i], h_in2[i]); + + if (std::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_holtsmark.cpp b/test/test_holtsmark.cpp index 6aa3577868..475f5400aa 100644 --- a/test/test_holtsmark.cpp +++ b/test/test_holtsmark.cpp @@ -16,6 +16,14 @@ # include #endif +#ifdef __clang__ +# pragma clang diagnostic push +# pragma clang diagnostic ignored "-Wliteral-range" +#elif defined(__GNUC__) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wliteral-range" +#endif + using boost::math::holtsmark_distribution; #ifndef BOOST_MATH_HAS_GPU_SUPPORT diff --git a/test/test_igamma.cpp b/test/test_igamma.cpp index 8e80c772c4..0ad7019963 100644 --- a/test/test_igamma.cpp +++ b/test/test_igamma.cpp @@ -3,7 +3,18 @@ // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +#ifndef SYCL_LANGUAGE_VERSION #include +#endif + +#ifdef __clang__ +# pragma clang diagnostic push +# pragma clang diagnostic ignored "-Wliteral-range" +#elif defined(__GNUC__) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wliteral-range" +#endif + #include "test_igamma.hpp" // diff --git a/test/test_igamma.hpp b/test/test_igamma.hpp index b434f727ee..bfe386d4de 100644 --- a/test/test_igamma.hpp +++ b/test/test_igamma.hpp @@ -8,11 +8,12 @@ #include #include +#include #define BOOST_TEST_MAIN #include #include +#include "../include_private/boost/math/tools/test.hpp" #include -#include #include #include #include diff --git a/test/test_igamma_inv.cpp b/test/test_igamma_inv.cpp index eafed0e1da..17e0bfb54f 100644 --- a/test/test_igamma_inv.cpp +++ b/test/test_igamma_inv.cpp @@ -3,7 +3,18 @@ // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +#ifndef SYCL_LANGUAGE_VERSION #include +#endif + +#ifdef __clang__ +# pragma clang diagnostic push +# pragma clang diagnostic ignored "-Wliteral-range" +#elif defined(__GNUC__) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wliteral-range" +#endif +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error #include "test_igamma_inv.hpp" #if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT) @@ -89,14 +100,22 @@ void expected_results() "linux.*", // platform largest_type, // test type(s) "[^|]*medium[^|]*", // test data group + #ifdef SYCL_LANGUAGE_VERSION + "[^|]*", 350, 50); + #else "[^|]*", 350, 5); // test function + #endif add_expected_result( "[^|]*", // compiler "[^|]*", // stdlib "linux.*", // platform largest_type, // test type(s) "[^|]*large[^|]*", // test data group + #ifdef SYCL_LANGUAGE_VERSION + "[^|]*", 150, 20); // test function + #else "[^|]*", 150, 5); // test function + #endif // diff --git a/test/test_igamma_inv.hpp b/test/test_igamma_inv.hpp index 7330e918a7..cf481537e7 100644 --- a/test/test_igamma_inv.hpp +++ b/test/test_igamma_inv.hpp @@ -6,13 +6,14 @@ #include #include +#include #define BOOST_TEST_MAIN #include #include #include #include #include -#include +#include "../include_private/boost/math/tools/test.hpp" #include #include #include diff --git a/test/test_igamma_inva.cpp b/test/test_igamma_inva.cpp index 047df11735..8d0e965962 100644 --- a/test/test_igamma_inva.cpp +++ b/test/test_igamma_inva.cpp @@ -3,7 +3,18 @@ // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +#ifndef SYCL_LANGUAGE_VERSION #include +#endif + +#ifdef __clang__ +# pragma clang diagnostic push +# pragma clang diagnostic ignored "-Wliteral-range" +#elif defined(__GNUC__) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wliteral-range" +#endif + #include "test_igamma_inva.hpp" #if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT) diff --git a/test/test_igamma_inva.hpp b/test/test_igamma_inva.hpp index 402ea2f8bc..d9d317da15 100644 --- a/test/test_igamma_inva.hpp +++ b/test/test_igamma_inva.hpp @@ -8,13 +8,14 @@ #include #include +#include #define BOOST_TEST_MAIN #include #include #include #include #include -#include +#include "../include_private/boost/math/tools/test.hpp" #include #include #include diff --git a/test/test_landau.cpp b/test/test_landau.cpp index 297589bc09..1625b21777 100644 --- a/test/test_landau.cpp +++ b/test/test_landau.cpp @@ -15,6 +15,14 @@ # include #endif +#ifdef __clang__ +# pragma clang diagnostic push +# pragma clang diagnostic ignored "-Wliteral-range" +#elif defined(__GNUC__) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wliteral-range" +#endif + using boost::math::landau_distribution; #ifndef BOOST_MATH_HAS_GPU_SUPPORT diff --git a/test/test_log1p_nvrtc_double.cpp b/test/test_log1p_nvrtc_double.cpp index 722194be6a..36b0771b1b 100644 --- a/test/test_log1p_nvrtc_double.cpp +++ b/test/test_log1p_nvrtc_double.cpp @@ -85,9 +85,9 @@ int main() nvrtcAddNameExpression(prog, "test_log1p_kernel"); #ifdef BOOST_MATH_NVRTC_CI_RUN - const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/"}; + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; #else - const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/"}; + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; #endif // Compile the program diff --git a/test/test_log1p_nvrtc_float.cpp b/test/test_log1p_nvrtc_float.cpp index 772c50cd10..7194ffb56a 100644 --- a/test/test_log1p_nvrtc_float.cpp +++ b/test/test_log1p_nvrtc_float.cpp @@ -7,6 +7,11 @@ #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error #define BOOST_MATH_PROMOTE_DOUBLE_POLICY false +// Must be included first +#include +#include +#include + #include #include #include @@ -14,9 +19,6 @@ #include #include #include -#include -#include -#include typedef float float_type; @@ -85,9 +87,9 @@ int main() nvrtcAddNameExpression(prog, "test_log1p_kernel"); #ifdef BOOST_MATH_NVRTC_CI_RUN - const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/"}; + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; #else - const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/"}; + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; #endif // Compile the program diff --git a/test/test_round.cpp b/test/test_round.cpp index e363efd565..e603aa510d 100644 --- a/test/test_round.cpp +++ b/test/test_round.cpp @@ -7,6 +7,11 @@ #include #endif +#ifdef __clang__ +# pragma clang diagnostic push +# pragma clang diagnostic ignored "-Wimplicit-const-int-float-conversion" +#endif + #include #define BOOST_TEST_MAIN #include diff --git a/test/test_tgamma_ratio_double.cu b/test/test_tgamma_ratio_double.cu new file mode 100644 index 0000000000..059e1c3c67 --- /dev/null +++ b/test/test_tgamma_ratio_double.cu @@ -0,0 +1,104 @@ + +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef double float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, const float_type *in2, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = boost::math::tgamma_ratio(in1[i], in2[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed input vector B + cuda_managed_ptr input_vector2(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = rand()/(float_type)RAND_MAX; + input_vector2[i] = rand()/(float_type)RAND_MAX; + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 256; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + + cuda_test<<>>(input_vector1.get(), input_vector2.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + + std::cout << "CUDA kernal done in: " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + + if (err != cudaSuccess) + { + std::cerr << "Failed to launch CUDA kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(boost::math::tgamma_ratio(input_vector1[i], input_vector2[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 10) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED, normal calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + + return 0; +} diff --git a/test/test_tgamma_ratio_float.cu b/test/test_tgamma_ratio_float.cu new file mode 100644 index 0000000000..dc669bd7fb --- /dev/null +++ b/test/test_tgamma_ratio_float.cu @@ -0,0 +1,104 @@ + +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef float float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, const float_type *in2, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = boost::math::tgamma_ratio(in1[i], in2[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed input vector B + cuda_managed_ptr input_vector2(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = rand()/(float_type)RAND_MAX; + input_vector2[i] = rand()/(float_type)RAND_MAX; + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 256; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + + cuda_test<<>>(input_vector1.get(), input_vector2.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + + std::cout << "CUDA kernal done in: " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + + if (err != cudaSuccess) + { + std::cerr << "Failed to launch CUDA kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(boost::math::tgamma_ratio(input_vector1[i], input_vector2[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 10) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED, normal calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + + return 0; +} diff --git a/test/test_tgamma_ratio_nvrtc_double.cpp b/test/test_tgamma_ratio_nvrtc_double.cpp new file mode 100644 index 0000000000..5b0c3b1e67 --- /dev/null +++ b/test/test_tgamma_ratio_nvrtc_double.cpp @@ -0,0 +1,190 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +typedef double float_type; + +const char* cuda_kernel = R"( +typedef double float_type; +#include +#include +extern "C" __global__ +void test_tgamma_ratio_kernel(const float_type *in1, const float_type *in2, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = boost::math::tgamma_ratio(in1[i], in2[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_tgamma_ratio_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_tgamma_ratio_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_tgamma_ratio_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1000.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + const auto res = boost::math::tgamma_ratio(h_in1[i], h_in2[i]); + + if (std::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_tgamma_ratio_nvrtc_float.cpp b/test/test_tgamma_ratio_nvrtc_float.cpp new file mode 100644 index 0000000000..ab1bf339b4 --- /dev/null +++ b/test/test_tgamma_ratio_nvrtc_float.cpp @@ -0,0 +1,190 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +typedef float float_type; + +const char* cuda_kernel = R"( +typedef float float_type; +#include +#include +extern "C" __global__ +void test_tgamma_ratio_kernel(const float_type *in1, const float_type *in2, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = boost::math::tgamma_ratio(in1[i], in2[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_tgamma_ratio_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_tgamma_ratio_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_tgamma_ratio_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1000.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + const auto res = boost::math::tgamma_ratio(h_in1[i], h_in2[i]); + + if (std::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_weibull.cpp b/test/test_weibull.cpp index 4b31a7f0b7..dc509b742b 100644 --- a/test/test_weibull.cpp +++ b/test/test_weibull.cpp @@ -12,15 +12,17 @@ # pragma warning (disable : 4127) // conditional expression is constant. #endif - +#include +#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS #include // for real_concept +#endif #define BOOST_TEST_MAIN #include // Boost.Test #include #include using boost::math::weibull_distribution; -#include +#include "../include_private/boost/math/tools/test.hpp" #include "test_out_of_range.hpp" #include diff --git a/test/test_weibull_cdf_double.cu b/test/test_weibull_cdf_double.cu new file mode 100644 index 0000000000..65efbe2528 --- /dev/null +++ b/test/test_weibull_cdf_double.cu @@ -0,0 +1,109 @@ +// Copyright John Maddock 2016. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef double float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = cdf(boost::math::weibull_distribution(1), in1[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + try{ + + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + boost::random::mt19937 gen; + boost::random::uniform_real_distribution dist; + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 512; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + cuda_test<<>>(input_vector1.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + if (err != cudaSuccess) + { + std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(cdf(boost::math::weibull_distribution(1), input_vector1[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 100.0) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + std::cerr << "Error rate was: " << boost::math::epsilon_difference(output_vector[i], results[i]) << "eps" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + } + return 0; +} \ No newline at end of file diff --git a/test/test_weibull_cdf_float.cu b/test/test_weibull_cdf_float.cu new file mode 100644 index 0000000000..65c3ce1ffc --- /dev/null +++ b/test/test_weibull_cdf_float.cu @@ -0,0 +1,109 @@ +// Copyright John Maddock 2016. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef float float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = cdf(boost::math::weibull_distribution(1), in1[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + try{ + + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + boost::random::mt19937 gen; + boost::random::uniform_real_distribution dist; + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 512; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + cuda_test<<>>(input_vector1.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + if (err != cudaSuccess) + { + std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(cdf(boost::math::weibull_distribution(1), input_vector1[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 100.0) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + std::cerr << "Error rate was: " << boost::math::epsilon_difference(output_vector[i], results[i]) << "eps" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + } + return 0; +} \ No newline at end of file diff --git a/test/test_weibull_cdf_nvrtc_double.cpp b/test/test_weibull_cdf_nvrtc_double.cpp new file mode 100644 index 0000000000..60d5ff5afb --- /dev/null +++ b/test/test_weibull_cdf_nvrtc_double.cpp @@ -0,0 +1,191 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +typedef double float_type; + +const char* cuda_kernel = R"( +typedef double float_type; +#include +#include +extern "C" __global__ +void test_weibull_kernel(const float_type *in1, const float_type*, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = cdf(boost::math::weibull_distribution(1), in1[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_weibull_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_weibull_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_weibull_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + auto res = cdf(boost::math::weibull_distribution(1), h_in1[i]); + + if (boost::math::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_weibull_cdf_nvrtc_float.cpp b/test/test_weibull_cdf_nvrtc_float.cpp new file mode 100644 index 0000000000..5085b2f7dd --- /dev/null +++ b/test/test_weibull_cdf_nvrtc_float.cpp @@ -0,0 +1,191 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +typedef float float_type; + +const char* cuda_kernel = R"( +typedef float float_type; +#include +#include +extern "C" __global__ +void test_weibull_kernel(const float_type *in1, const float_type*, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = cdf(boost::math::weibull_distribution(1), in1[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_weibull_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_weibull_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_weibull_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + auto res = cdf(boost::math::weibull_distribution(1), h_in1[i]); + + if (boost::math::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_weibull_pdf_double.cu b/test/test_weibull_pdf_double.cu new file mode 100644 index 0000000000..645df4c0ad --- /dev/null +++ b/test/test_weibull_pdf_double.cu @@ -0,0 +1,109 @@ +// Copyright John Maddock 2016. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef double float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = pdf(boost::math::weibull_distribution(1), in1[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + try{ + + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + boost::random::mt19937 gen; + boost::random::uniform_real_distribution dist; + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 512; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + cuda_test<<>>(input_vector1.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + if (err != cudaSuccess) + { + std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(pdf(boost::math::weibull_distribution(1), input_vector1[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 100.0) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + std::cerr << "Error rate was: " << boost::math::epsilon_difference(output_vector[i], results[i]) << "eps" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + } + return 0; +} \ No newline at end of file diff --git a/test/test_weibull_pdf_float.cu b/test/test_weibull_pdf_float.cu new file mode 100644 index 0000000000..f1e6917f03 --- /dev/null +++ b/test/test_weibull_pdf_float.cu @@ -0,0 +1,109 @@ +// Copyright John Maddock 2016. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef float float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = pdf(boost::math::weibull_distribution(1), in1[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + try{ + + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + boost::random::mt19937 gen; + boost::random::uniform_real_distribution dist; + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 512; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + cuda_test<<>>(input_vector1.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + if (err != cudaSuccess) + { + std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(pdf(boost::math::weibull_distribution(1), input_vector1[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 100.0) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + std::cerr << "Error rate was: " << boost::math::epsilon_difference(output_vector[i], results[i]) << "eps" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + } + return 0; +} \ No newline at end of file diff --git a/test/test_weibull_pdf_nvrtc_double.cpp b/test/test_weibull_pdf_nvrtc_double.cpp new file mode 100644 index 0000000000..2e5e237b20 --- /dev/null +++ b/test/test_weibull_pdf_nvrtc_double.cpp @@ -0,0 +1,191 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +typedef double float_type; + +const char* cuda_kernel = R"( +typedef double float_type; +#include +#include +extern "C" __global__ +void test_weibull_kernel(const float_type *in1, const float_type*, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = pdf(boost::math::weibull_distribution(1), in1[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_weibull_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_weibull_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_weibull_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + auto res = pdf(boost::math::weibull_distribution(1), h_in1[i]); + + if (boost::math::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_weibull_pdf_nvrtc_float.cpp b/test/test_weibull_pdf_nvrtc_float.cpp new file mode 100644 index 0000000000..6c3c5202c1 --- /dev/null +++ b/test/test_weibull_pdf_nvrtc_float.cpp @@ -0,0 +1,191 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +typedef float float_type; + +const char* cuda_kernel = R"( +typedef float float_type; +#include +#include +extern "C" __global__ +void test_weibull_kernel(const float_type *in1, const float_type*, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = pdf(boost::math::weibull_distribution(1), in1[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_weibull_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_weibull_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_weibull_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + auto res = pdf(boost::math::weibull_distribution(1), h_in1[i]); + + if (boost::math::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_weibull_quan_double.cu b/test/test_weibull_quan_double.cu new file mode 100644 index 0000000000..2f05006022 --- /dev/null +++ b/test/test_weibull_quan_double.cu @@ -0,0 +1,109 @@ +// Copyright John Maddock 2016. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef double float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = quantile(boost::math::weibull_distribution(1), in1[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + try{ + + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + boost::random::mt19937 gen; + boost::random::uniform_real_distribution dist; + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 512; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + cuda_test<<>>(input_vector1.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + if (err != cudaSuccess) + { + std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(quantile(boost::math::weibull_distribution(1), input_vector1[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 100.0) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + std::cerr << "Error rate was: " << boost::math::epsilon_difference(output_vector[i], results[i]) << "eps" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + } + return 0; +} \ No newline at end of file diff --git a/test/test_weibull_quan_float.cu b/test/test_weibull_quan_float.cu new file mode 100644 index 0000000000..3027e14dda --- /dev/null +++ b/test/test_weibull_quan_float.cu @@ -0,0 +1,109 @@ +// Copyright John Maddock 2016. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#include +#include +#include +#include +#include +#include "cuda_managed_ptr.hpp" +#include "stopwatch.hpp" + +// For the CUDA runtime routines (prefixed with "cuda_") +#include + +typedef float float_type; + +/** + * CUDA Kernel Device code + * + */ +__global__ void cuda_test(const float_type *in1, float_type *out, int numElements) +{ + using std::cos; + int i = blockDim.x * blockIdx.x + threadIdx.x; + + if (i < numElements) + { + out[i] = quantile(boost::math::weibull_distribution(1), in1[i]); + } +} + +/** + * Host main routine + */ +int main(void) +{ + try{ + + // Error code to check return values for CUDA calls + cudaError_t err = cudaSuccess; + + // Print the vector length to be used, and compute its size + int numElements = 50000; + std::cout << "[Vector operation on " << numElements << " elements]" << std::endl; + + // Allocate the managed input vector A + cuda_managed_ptr input_vector1(numElements); + + // Allocate the managed output vector C + cuda_managed_ptr output_vector(numElements); + + boost::random::mt19937 gen; + boost::random::uniform_real_distribution dist; + // Initialize the input vectors + for (int i = 0; i < numElements; ++i) + { + input_vector1[i] = dist(gen); + } + + // Launch the Vector Add CUDA Kernel + int threadsPerBlock = 512; + int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock; + std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl; + + watch w; + cuda_test<<>>(input_vector1.get(), output_vector.get(), numElements); + cudaDeviceSynchronize(); + std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl; + + err = cudaGetLastError(); + if (err != cudaSuccess) + { + std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl; + return EXIT_FAILURE; + } + + // Verify that the result vector is correct + std::vector results; + results.reserve(numElements); + w.reset(); + for(int i = 0; i < numElements; ++i) + results.push_back(quantile(boost::math::weibull_distribution(1), input_vector1[i])); + double t = w.elapsed(); + // check the results + for(int i = 0; i < numElements; ++i) + { + if (boost::math::epsilon_difference(output_vector[i], results[i]) > 100.0) + { + std::cerr << "Result verification failed at element " << i << "!" << std::endl; + std::cerr << "Error rate was: " << boost::math::epsilon_difference(output_vector[i], results[i]) << "eps" << std::endl; + return EXIT_FAILURE; + } + } + + std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl; + std::cout << "Done\n"; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + } + return 0; +} \ No newline at end of file diff --git a/test/test_weibull_quan_nvrtc_double.cpp b/test/test_weibull_quan_nvrtc_double.cpp new file mode 100644 index 0000000000..aed31865e8 --- /dev/null +++ b/test/test_weibull_quan_nvrtc_double.cpp @@ -0,0 +1,191 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +typedef double float_type; + +const char* cuda_kernel = R"( +typedef double float_type; +#include +#include +extern "C" __global__ +void test_weibull_kernel(const float_type *in1, const float_type*, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = quantile(boost::math::weibull_distribution(1), in1[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_weibull_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_weibull_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_weibull_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + auto res = quantile(boost::math::weibull_distribution(1), h_in1[i]); + + if (boost::math::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +} diff --git a/test/test_weibull_quan_nvrtc_float.cpp b/test/test_weibull_quan_nvrtc_float.cpp new file mode 100644 index 0000000000..98997b354b --- /dev/null +++ b/test/test_weibull_quan_nvrtc_float.cpp @@ -0,0 +1,191 @@ +// Copyright John Maddock 2016. +// Copyright Matt Borland 2024. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error +#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false + +// Must be included first +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +typedef float float_type; + +const char* cuda_kernel = R"( +typedef float float_type; +#include +#include +extern "C" __global__ +void test_weibull_kernel(const float_type *in1, const float_type*, float_type *out, int numElements) +{ + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < numElements) + { + out[i] = quantile(boost::math::weibull_distribution(1), in1[i]); + } +} +)"; + +void checkCUDAError(cudaError_t result, const char* msg) +{ + if (result != cudaSuccess) + { + std::cerr << msg << ": " << cudaGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkCUError(CUresult result, const char* msg) +{ + if (result != CUDA_SUCCESS) + { + const char* errorStr; + cuGetErrorString(result, &errorStr); + std::cerr << msg << ": " << errorStr << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkNVRTCError(nvrtcResult result, const char* msg) +{ + if (result != NVRTC_SUCCESS) + { + std::cerr << msg << ": " << nvrtcGetErrorString(result) << std::endl; + exit(EXIT_FAILURE); + } +} + +int main() +{ + try + { + // Initialize CUDA driver API + checkCUError(cuInit(0), "Failed to initialize CUDA"); + + // Create CUDA context + CUcontext context; + CUdevice device; + checkCUError(cuDeviceGet(&device, 0), "Failed to get CUDA device"); + checkCUError(cuCtxCreate(&context, 0, device), "Failed to create CUDA context"); + + nvrtcProgram prog; + nvrtcResult res; + + res = nvrtcCreateProgram(&prog, cuda_kernel, "test_weibull_kernel.cu", 0, nullptr, nullptr); + checkNVRTCError(res, "Failed to create NVRTC program"); + + nvrtcAddNameExpression(prog, "test_weibull_kernel"); + + #ifdef BOOST_MATH_NVRTC_CI_RUN + const char* opts[] = {"--std=c++14", "--gpu-architecture=compute_75", "--include-path=/home/runner/work/cuda-math/boost-root/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #else + const char* opts[] = {"--std=c++14", "--include-path=/home/mborland/Documents/boost/libs/cuda-math/include/", "-I/usr/local/cuda/include"}; + #endif + + // Compile the program + res = nvrtcCompileProgram(prog, sizeof(opts) / sizeof(const char*), opts); + if (res != NVRTC_SUCCESS) + { + size_t log_size; + nvrtcGetProgramLogSize(prog, &log_size); + char* log = new char[log_size]; + nvrtcGetProgramLog(prog, log); + std::cerr << "Compilation failed:\n" << log << std::endl; + delete[] log; + exit(EXIT_FAILURE); + } + + // Get PTX from the program + size_t ptx_size; + nvrtcGetPTXSize(prog, &ptx_size); + char* ptx = new char[ptx_size]; + nvrtcGetPTX(prog, ptx); + + // Load PTX into CUDA module + CUmodule module; + CUfunction kernel; + checkCUError(cuModuleLoadDataEx(&module, ptx, 0, 0, 0), "Failed to load module"); + checkCUError(cuModuleGetFunction(&kernel, module, "test_weibull_kernel"), "Failed to get kernel function"); + + int numElements = 5000; + float_type *h_in1, *h_in2, *h_out; + float_type *d_in1, *d_in2, *d_out; + + // Allocate memory on the host + h_in1 = new float_type[numElements]; + h_in2 = new float_type[numElements]; + h_out = new float_type[numElements]; + + // Initialize input arrays + std::mt19937_64 rng(42); + std::uniform_real_distribution dist(0.0f, 1.0f); + for (int i = 0; i < numElements; ++i) + { + h_in1[i] = static_cast(dist(rng)); + h_in2[i] = static_cast(dist(rng)); + } + + checkCUDAError(cudaMalloc(&d_in1, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in1"); + checkCUDAError(cudaMalloc(&d_in2, numElements * sizeof(float_type)), "Failed to allocate device memory for d_in2"); + checkCUDAError(cudaMalloc(&d_out, numElements * sizeof(float_type)), "Failed to allocate device memory for d_out"); + + checkCUDAError(cudaMemcpy(d_in1, h_in1, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in1"); + checkCUDAError(cudaMemcpy(d_in2, h_in2, numElements * sizeof(float_type), cudaMemcpyHostToDevice), "Failed to copy data to device for d_in2"); + + int blockSize = 256; + int numBlocks = (numElements + blockSize - 1) / blockSize; + void* args[] = { &d_in1, &d_in2, &d_out, &numElements }; + checkCUError(cuLaunchKernel(kernel, numBlocks, 1, 1, blockSize, 1, 1, 0, 0, args, 0), "Kernel launch failed"); + + checkCUDAError(cudaMemcpy(h_out, d_out, numElements * sizeof(float_type), cudaMemcpyDeviceToHost), "Failed to copy data back to host for h_out"); + + // Verify Result + for (int i = 0; i < numElements; ++i) + { + auto res = quantile(boost::math::weibull_distribution(1), h_in1[i]); + + if (boost::math::isfinite(res)) + { + if (boost::math::epsilon_difference(res, h_out[i]) > 300) + { + std::cout << "error at line: " << i + << "\nParallel: " << h_out[i] + << "\n Serial: " << res + << "\n Dist: " << boost::math::epsilon_difference(res, h_out[i]) << std::endl; + } + } + } + + cudaFree(d_in1); + cudaFree(d_in2); + cudaFree(d_out); + delete[] h_in1; + delete[] h_in2; + delete[] h_out; + + nvrtcDestroyProgram(&prog); + delete[] ptx; + + cuCtxDestroy(context); + + std::cout << "Kernel executed successfully." << std::endl; + return 0; + } + catch(const std::exception& e) + { + std::cerr << "Stopped with exception: " << e.what() << std::endl; + return EXIT_FAILURE; + } +}