diff --git a/include/boost/math/distributions/arcsine.hpp b/include/boost/math/distributions/arcsine.hpp index 8aad5b2d0..ac5d6e565 100644 --- a/include/boost/math/distributions/arcsine.hpp +++ b/include/boost/math/distributions/arcsine.hpp @@ -54,7 +54,7 @@ namespace boost // Common error checking routines for arcsine distribution functions: // Duplicating for x_min and x_max provides specific error messages. template - inline bool check_x_min(const char* function, const RealType& x, RealType* result, const Policy& pol) + inline BOOST_GPU_ENABLED bool check_x_min(const char* function, const RealType& x, RealType* result, const Policy& pol) { if (!(boost::math::isfinite)(x)) { @@ -67,7 +67,7 @@ namespace boost } // bool check_x_min template - inline bool check_x_max(const char* function, const RealType& x, RealType* result, const Policy& pol) + inline BOOST_GPU_ENABLED bool check_x_max(const char* function, const RealType& x, RealType* result, const Policy& pol) { if (!(boost::math::isfinite)(x)) { @@ -81,11 +81,16 @@ namespace boost template - inline bool check_x_minmax(const char* function, const RealType& x_min, const RealType& x_max, RealType* result, const Policy& pol) + inline BOOST_GPU_ENABLED bool check_x_minmax(const char* function, const RealType& x_min, const RealType& x_max, RealType* result, const Policy& pol) { // Check x_min < x_max if (x_min >= x_max) { - std::string msg = "x_max argument is %1%, but must be > x_min = " + lexical_cast(x_min) + "!"; +#ifdef __CUDA_ARCH__ + *result = policies::raise_domain_error( + function, + "x_max value %1% is out of range", x_max, pol); +#else + std::string msg = "x_max argument is %1%, but must be > x_min = " + lexical_cast(x_min) + "!"; *result = policies::raise_domain_error( function, msg.c_str(), x_max, pol); @@ -94,13 +99,14 @@ namespace boost // But would require replication of all helpers functions in /policies/error_handling.hpp for two values, // as well as two value versions of raise_error, raise_domain_error and do_format ... // so use slightly hacky lexical_cast to string instead. +#endif return false; } return true; } // bool check_x_minmax template - inline bool check_prob(const char* function, const RealType& p, RealType* result, const Policy& pol) + inline BOOST_GPU_ENABLED bool check_prob(const char* function, const RealType& p, RealType* result, const Policy& pol) { if ((p < 0) || (p > 1) || !(boost::math::isfinite)(p)) { @@ -113,7 +119,7 @@ namespace boost } // bool check_prob template - inline bool check_x(const char* function, const RealType& x_min, const RealType& x_max, const RealType& x, RealType* result, const Policy& pol) + inline BOOST_GPU_ENABLED bool check_x(const char* function, const RealType& x_min, const RealType& x_max, const RealType& x, RealType* result, const Policy& pol) { // Check x finite and x_min < x < x_max. if (!(boost::math::isfinite)(x)) { @@ -137,7 +143,7 @@ namespace boost } // bool check_x template - inline bool check_dist(const char* function, const RealType& x_min, const RealType& x_max, RealType* result, const Policy& pol) + inline BOOST_GPU_ENABLED bool check_dist(const char* function, const RealType& x_min, const RealType& x_max, RealType* result, const Policy& pol) { // Check both x_min and x_max finite, and x_min < x_max. return check_x_min(function, x_min, result, pol) && check_x_max(function, x_max, result, pol) @@ -145,14 +151,14 @@ namespace boost } // bool check_dist template - inline bool check_dist_and_x(const char* function, const RealType& x_min, const RealType& x_max, RealType x, RealType* result, const Policy& pol) + inline BOOST_GPU_ENABLED bool check_dist_and_x(const char* function, const RealType& x_min, const RealType& x_max, RealType x, RealType* result, const Policy& pol) { return check_dist(function, x_min, x_max, result, pol) && arcsine_detail::check_x(function, x_min, x_max, x, result, pol); } // bool check_dist_and_x template - inline bool check_dist_and_prob(const char* function, const RealType& x_min, const RealType& x_max, RealType p, RealType* result, const Policy& pol) + inline BOOST_GPU_ENABLED bool check_dist_and_prob(const char* function, const RealType& x_min, const RealType& x_max, RealType p, RealType* result, const Policy& pol) { return check_dist(function, x_min, x_max, result, pol) && check_prob(function, p, result, pol); @@ -167,7 +173,7 @@ namespace boost typedef RealType value_type; typedef Policy policy_type; - arcsine_distribution(RealType x_min = 0, RealType x_max = 1) : m_x_min(x_min), m_x_max(x_max) + BOOST_GPU_ENABLED arcsine_distribution(RealType x_min = 0, RealType x_max = 1) : m_x_min(x_min), m_x_max(x_max) { // Default beta (alpha = beta = 0.5) is standard arcsine with x_min = 0, x_max = 1. // Generalized to allow x_min and x_max to be specified. RealType result; @@ -178,11 +184,11 @@ namespace boost &result, Policy()); } // arcsine_distribution constructor. // Accessor functions: - RealType x_min() const + BOOST_GPU_ENABLED RealType x_min() const { return m_x_min; } - RealType x_max() const + BOOST_GPU_ENABLED RealType x_max() const { return m_x_max; } @@ -197,21 +203,21 @@ namespace boost template - inline const std::pair range(const arcsine_distribution& dist) + inline BOOST_GPU_ENABLED const std::pair range(const arcsine_distribution& dist) { // Range of permissible values for random variable x. using boost::math::tools::max_value; return std::pair(static_cast(dist.x_min()), static_cast(dist.x_max())); } template - inline const std::pair support(const arcsine_distribution& dist) + inline BOOST_GPU_ENABLED const std::pair support(const arcsine_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(dist.x_min()), static_cast(dist.x_max())); } template - inline RealType mean(const arcsine_distribution& dist) + inline BOOST_GPU_ENABLED RealType mean(const arcsine_distribution& dist) { // Mean of arcsine distribution . RealType result; RealType x_min = dist.x_min(); @@ -230,7 +236,7 @@ namespace boost } // mean template - inline RealType variance(const arcsine_distribution& dist) + inline BOOST_GPU_ENABLED RealType variance(const arcsine_distribution& dist) { // Variance of standard arcsine distribution = (1-0)/8 = 0.125. RealType result; RealType x_min = dist.x_min(); @@ -248,7 +254,7 @@ namespace boost } // variance template - inline RealType mode(const arcsine_distribution& /* dist */) + inline BOOST_GPU_ENABLED RealType mode(const arcsine_distribution& /* dist */) { //There are always [*two] values for the mode, at ['x_min] and at ['x_max], default 0 and 1, // so instead we raise the exception domain_error. return policies::raise_domain_error( @@ -259,7 +265,7 @@ namespace boost } // mode template - inline RealType median(const arcsine_distribution& dist) + inline BOOST_GPU_ENABLED RealType median(const arcsine_distribution& dist) { // Median of arcsine distribution (a + b) / 2 == mean. RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); @@ -277,7 +283,7 @@ namespace boost } template - inline RealType skewness(const arcsine_distribution& dist) + inline BOOST_GPU_ENABLED RealType skewness(const arcsine_distribution& dist) { RealType result; RealType x_min = dist.x_min(); @@ -296,7 +302,7 @@ namespace boost } // skewness template - inline RealType kurtosis_excess(const arcsine_distribution& dist) + inline BOOST_GPU_ENABLED RealType kurtosis_excess(const arcsine_distribution& dist) { RealType result; RealType x_min = dist.x_min(); @@ -316,7 +322,7 @@ namespace boost } // kurtosis_excess template - inline RealType kurtosis(const arcsine_distribution& dist) + inline BOOST_GPU_ENABLED RealType kurtosis(const arcsine_distribution& dist) { RealType result; RealType x_min = dist.x_min(); @@ -336,12 +342,12 @@ namespace boost } // kurtosis template - inline RealType pdf(const arcsine_distribution& dist, const RealType& xx) + inline BOOST_GPU_ENABLED RealType pdf(const arcsine_distribution& dist, const RealType& xx) { // Probability Density/Mass Function arcsine. BOOST_FPU_EXCEPTION_GUARD BOOST_MATH_STD_USING // For ADL of std functions. - static const char* function = "boost::math::pdf(arcsine_distribution<%1%> const&, %1%)"; + BOOST_MATH_GPU_STATIC const char* function = "boost::math::pdf(arcsine_distribution<%1%> const&, %1%)"; RealType lo = dist.x_min(); RealType hi = dist.x_max(); @@ -362,11 +368,11 @@ namespace boost } // pdf template - inline RealType cdf(const arcsine_distribution& dist, const RealType& x) + inline BOOST_GPU_ENABLED RealType cdf(const arcsine_distribution& dist, const RealType& x) { // Cumulative Distribution Function arcsine. BOOST_MATH_STD_USING // For ADL of std functions. - static const char* function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; + BOOST_MATH_GPU_STATIC const char* function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; RealType x_min = dist.x_min(); RealType x_max = dist.x_max(); @@ -395,10 +401,10 @@ namespace boost } // arcsine cdf template - inline RealType cdf(const complemented2_type, RealType>& c) + inline BOOST_GPU_ENABLED RealType cdf(const complemented2_type, RealType>& c) { // Complemented Cumulative Distribution Function arcsine. BOOST_MATH_STD_USING // For ADL of std functions. - static const char* function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; + BOOST_MATH_GPU_STATIC const char* function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; RealType x = c.param; arcsine_distribution const& dist = c.dist; @@ -431,7 +437,7 @@ namespace boost } // arcine ccdf template - inline RealType quantile(const arcsine_distribution& dist, const RealType& p) + inline BOOST_GPU_ENABLED RealType quantile(const arcsine_distribution& dist, const RealType& p) { // Quantile or Percent Point arcsine function or // Inverse Cumulative probability distribution function CDF. @@ -445,7 +451,7 @@ namespace boost using boost::math::constants::half_pi; - static const char* function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; + BOOST_MATH_GPU_STATIC const char* function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; RealType result = 0; // of argument checks: RealType x_min = dist.x_min(); @@ -475,7 +481,7 @@ namespace boost } // quantile template - inline RealType quantile(const complemented2_type, RealType>& c) + inline BOOST_GPU_ENABLED RealType quantile(const complemented2_type, RealType>& c) { // Complement Quantile or Percent Point arcsine function. // Return the number of expected x for a given @@ -483,7 +489,7 @@ namespace boost BOOST_MATH_STD_USING // For ADL of std functions. using boost::math::constants::half_pi; - static const char* function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; + BOOST_MATH_GPU_STATIC const char* function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; // Error checks: RealType q = c.param; diff --git a/include/boost/math/distributions/cauchy.hpp b/include/boost/math/distributions/cauchy.hpp index 5a3a64f0f..74a052bf9 100644 --- a/include/boost/math/distributions/cauchy.hpp +++ b/include/boost/math/distributions/cauchy.hpp @@ -31,7 +31,7 @@ namespace detail { template -RealType cdf_imp(const cauchy_distribution& dist, const RealType& x, bool complement) +BOOST_GPU_ENABLED RealType cdf_imp(const cauchy_distribution& dist, const RealType& x, bool complement) { // // This calculates the cdf of the Cauchy distribution and/or its complement. @@ -55,7 +55,7 @@ RealType cdf_imp(const cauchy_distribution& dist, const RealTy // to get the result. // BOOST_MATH_STD_USING // for ADL of std functions - static const char* function = "boost::math::cdf(cauchy<%1%>&, %1%)"; + BOOST_MATH_GPU_STATIC const char* function = "boost::math::cdf(cauchy<%1%>&, %1%)"; RealType result = 0; RealType location = dist.location(); RealType scale = dist.scale(); @@ -67,6 +67,12 @@ RealType cdf_imp(const cauchy_distribution& dist, const RealTy { return result; } +#ifdef __CUDA_ARCH__ + if(x > tools::max_value()) + return static_cast((complement) ? 0 : 1); + if(x < -tools::max_value()) + return static_cast((complement) ? 1 : 0); +#else if(std::numeric_limits::has_infinity && x == std::numeric_limits::infinity()) { // cdf +infinity is unity. return static_cast((complement) ? 0 : 1); @@ -75,6 +81,7 @@ RealType cdf_imp(const cauchy_distribution& dist, const RealTy { // cdf -infinity is zero. return static_cast((complement) ? 1 : 0); } +#endif if(false == detail::check_x(function, x, &result, Policy())) { // Catches x == NaN return result; @@ -89,7 +96,7 @@ RealType cdf_imp(const cauchy_distribution& dist, const RealTy } // cdf template -RealType quantile_imp( +BOOST_GPU_ENABLED RealType quantile_imp( const cauchy_distribution& dist, const RealType& p, bool complement) @@ -102,7 +109,7 @@ RealType quantile_imp( // mid-point of the distribution. This is either added or subtracted // from the location parameter depending on whether `complement` is true. // - static const char* function = "boost::math::quantile(cauchy<%1%>&, %1%)"; + BOOST_MATH_GPU_STATIC const char* function = "boost::math::quantile(cauchy<%1%>&, %1%)"; BOOST_MATH_STD_USING // for ADL of std functions RealType result = 0; @@ -152,20 +159,20 @@ public: typedef RealType value_type; typedef Policy policy_type; - cauchy_distribution(RealType l_location = 0, RealType l_scale = 1) + BOOST_GPU_ENABLED cauchy_distribution(RealType l_location = 0, RealType l_scale = 1) : m_a(l_location), m_hg(l_scale) { - static const char* function = "boost::math::cauchy_distribution<%1%>::cauchy_distribution"; + BOOST_MATH_GPU_STATIC const char* function = "boost::math::cauchy_distribution<%1%>::cauchy_distribution"; RealType result; detail::check_location(function, l_location, &result, Policy()); detail::check_scale(function, l_scale, &result, Policy()); } // cauchy_distribution - RealType location()const + BOOST_GPU_ENABLED RealType location()const { return m_a; } - RealType scale()const + BOOST_GPU_ENABLED RealType scale()const { return m_hg; } @@ -178,7 +185,7 @@ private: typedef cauchy_distribution cauchy; template -inline const std::pair range(const cauchy_distribution&) +inline BOOST_GPU_ENABLED const std::pair range(const cauchy_distribution&) { // Range of permissible values for random variable x. if (std::numeric_limits::has_infinity) { @@ -192,7 +199,7 @@ inline const std::pair range(const cauchy_distribution -inline const std::pair support(const cauchy_distribution& ) +inline BOOST_GPU_ENABLED const std::pair support(const cauchy_distribution& ) { // 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. if (std::numeric_limits::has_infinity) @@ -207,11 +214,11 @@ inline const std::pair support(const cauchy_distribution -inline RealType pdf(const cauchy_distribution& dist, const RealType& x) +inline BOOST_GPU_ENABLED RealType pdf(const cauchy_distribution& dist, const RealType& x) { BOOST_MATH_STD_USING // for ADL of std functions - static const char* function = "boost::math::pdf(cauchy<%1%>&, %1%)"; + BOOST_MATH_GPU_STATIC const char* function = "boost::math::pdf(cauchy<%1%>&, %1%)"; RealType result = 0; RealType location = dist.location(); RealType scale = dist.scale(); @@ -244,31 +251,31 @@ inline RealType pdf(const cauchy_distribution& dist, const Rea } // pdf template -inline RealType cdf(const cauchy_distribution& dist, const RealType& x) +inline BOOST_GPU_ENABLED RealType cdf(const cauchy_distribution& dist, const RealType& x) { return detail::cdf_imp(dist, x, false); } // cdf template -inline RealType quantile(const cauchy_distribution& dist, const RealType& p) +inline BOOST_GPU_ENABLED RealType quantile(const cauchy_distribution& dist, const RealType& p) { return detail::quantile_imp(dist, p, false); } // quantile template -inline RealType cdf(const complemented2_type, RealType>& c) +inline BOOST_GPU_ENABLED RealType cdf(const complemented2_type, RealType>& c) { return detail::cdf_imp(c.dist, c.param, true); } // cdf complement template -inline RealType quantile(const complemented2_type, RealType>& c) +inline BOOST_GPU_ENABLED RealType quantile(const complemented2_type, RealType>& c) { return detail::quantile_imp(c.dist, c.param, true); } // quantile complement template -inline RealType mean(const cauchy_distribution&) +inline BOOST_GPU_ENABLED RealType mean(const cauchy_distribution&) { // There is no mean: typedef typename Policy::assert_undefined_type assert_type; BOOST_STATIC_ASSERT(assert_type::value == 0); @@ -281,7 +288,7 @@ inline RealType mean(const cauchy_distribution&) } template -inline RealType variance(const cauchy_distribution& /*dist*/) +inline BOOST_GPU_ENABLED RealType variance(const cauchy_distribution& /*dist*/) { // There is no variance: typedef typename Policy::assert_undefined_type assert_type; @@ -295,18 +302,18 @@ inline RealType variance(const cauchy_distribution& /*dist*/) } template -inline RealType mode(const cauchy_distribution& dist) +inline BOOST_GPU_ENABLED RealType mode(const cauchy_distribution& dist) { return dist.location(); } template -inline RealType median(const cauchy_distribution& dist) +inline BOOST_GPU_ENABLED RealType median(const cauchy_distribution& dist) { return dist.location(); } template -inline RealType skewness(const cauchy_distribution& /*dist*/) +inline BOOST_GPU_ENABLED RealType skewness(const cauchy_distribution& /*dist*/) { // There is no skewness: typedef typename Policy::assert_undefined_type assert_type; @@ -320,7 +327,7 @@ inline RealType skewness(const cauchy_distribution& /*dist*/) } template -inline RealType kurtosis(const cauchy_distribution& /*dist*/) +inline BOOST_GPU_ENABLED RealType kurtosis(const cauchy_distribution& /*dist*/) { // There is no kurtosis: typedef typename Policy::assert_undefined_type assert_type; @@ -334,7 +341,7 @@ inline RealType kurtosis(const cauchy_distribution& /*dist*/) } template -inline RealType kurtosis_excess(const cauchy_distribution& /*dist*/) +inline BOOST_GPU_ENABLED RealType kurtosis_excess(const cauchy_distribution& /*dist*/) { // There is no kurtosis excess: typedef typename Policy::assert_undefined_type assert_type; diff --git a/include/boost/math/distributions/detail/common_error_handling.hpp b/include/boost/math/distributions/detail/common_error_handling.hpp index 71a771e5a..5f3ba4305 100644 --- a/include/boost/math/distributions/detail/common_error_handling.hpp +++ b/include/boost/math/distributions/detail/common_error_handling.hpp @@ -18,7 +18,7 @@ namespace boost{ namespace math{ namespace detail { template -inline bool check_probability(const char* function, RealType const& prob, RealType* result, const Policy& pol) +inline BOOST_GPU_ENABLED bool check_probability(const char* function, RealType const& prob, RealType* result, const Policy& pol) { if((prob < 0) || (prob > 1) || !(boost::math::isfinite)(prob)) { @@ -31,7 +31,7 @@ inline bool check_probability(const char* function, RealType const& prob, RealTy } template -inline bool check_df(const char* function, RealType const& df, RealType* result, const Policy& pol) +inline BOOST_GPU_ENABLED bool check_df(const char* function, RealType const& df, RealType* result, const Policy& pol) { // df > 0 but NOT +infinity allowed. if((df <= 0) || !(boost::math::isfinite)(df)) { @@ -44,7 +44,7 @@ inline bool check_df(const char* function, RealType const& df, RealType* result, } template -inline bool check_df_gt0_to_inf(const char* function, RealType const& df, RealType* result, const Policy& pol) +inline BOOST_GPU_ENABLED bool check_df_gt0_to_inf(const char* function, RealType const& df, RealType* result, const Policy& pol) { // df > 0 or +infinity are allowed. if( (df <= 0) || (boost::math::isnan)(df) ) { // is bad df <= 0 or NaN or -infinity. @@ -58,7 +58,7 @@ inline bool check_df_gt0_to_inf(const char* function, RealType const& df, RealTy template -inline bool check_scale( +inline BOOST_GPU_ENABLED bool check_scale( const char* function, RealType scale, RealType* result, @@ -75,7 +75,7 @@ inline bool check_scale( } template -inline bool check_location( +inline BOOST_GPU_ENABLED bool check_location( const char* function, RealType location, RealType* result, @@ -92,7 +92,7 @@ inline bool check_location( } template -inline bool check_x( +inline BOOST_GPU_ENABLED bool check_x( const char* function, RealType x, RealType* result, @@ -113,7 +113,7 @@ inline bool check_x( } // bool check_x template -inline bool check_x_gt0( +inline BOOST_GPU_ENABLED bool check_x_gt0( const char* function, RealType x, RealType* result, @@ -134,7 +134,7 @@ inline bool check_x_gt0( } // bool check_x_gt0 template -inline bool check_positive_x( +inline BOOST_GPU_ENABLED bool check_positive_x( const char* function, RealType x, RealType* result, @@ -154,7 +154,7 @@ inline bool check_positive_x( } template -inline bool check_non_centrality( +inline BOOST_GPU_ENABLED bool check_non_centrality( const char* function, RealType ncp, RealType* result, @@ -171,7 +171,7 @@ inline bool check_non_centrality( } template -inline bool check_finite( +inline BOOST_GPU_ENABLED bool check_finite( const char* function, RealType x, RealType* result, diff --git a/include/boost/math/distributions/detail/derived_accessors.hpp b/include/boost/math/distributions/detail/derived_accessors.hpp index 00f5a9325..f9620677e 100644 --- a/include/boost/math/distributions/detail/derived_accessors.hpp +++ b/include/boost/math/distributions/detail/derived_accessors.hpp @@ -39,24 +39,24 @@ namespace boost{ namespace math{ template -typename Distribution::value_type variance(const Distribution& dist); +BOOST_GPU_ENABLED typename Distribution::value_type variance(const Distribution& dist); template -inline typename Distribution::value_type standard_deviation(const Distribution& dist) +inline BOOST_GPU_ENABLED typename Distribution::value_type standard_deviation(const Distribution& dist) { BOOST_MATH_STD_USING // ADL of sqrt. return sqrt(variance(dist)); } template -inline typename Distribution::value_type variance(const Distribution& dist) +inline BOOST_GPU_ENABLED typename Distribution::value_type variance(const Distribution& dist) { typename Distribution::value_type result = standard_deviation(dist); return result * result; } template -inline typename Distribution::value_type hazard(const Distribution& dist, const RealType& x) +inline BOOST_GPU_ENABLED typename Distribution::value_type hazard(const Distribution& dist, const RealType& x) { // hazard function // http://www.itl.nist.gov/div898/handbook/eda/section3/eda362.htm#HAZ typedef typename Distribution::value_type value_type; @@ -75,7 +75,7 @@ inline typename Distribution::value_type hazard(const Distribution& dist, const } template -inline typename Distribution::value_type chf(const Distribution& dist, const RealType& x) +inline BOOST_GPU_ENABLED typename Distribution::value_type chf(const Distribution& dist, const RealType& x) { // cumulative hazard function. // http://www.itl.nist.gov/div898/handbook/eda/section3/eda362.htm#HAZ BOOST_MATH_STD_USING @@ -83,7 +83,7 @@ inline typename Distribution::value_type chf(const Distribution& dist, const Rea } template -inline typename Distribution::value_type coefficient_of_variation(const Distribution& dist) +inline BOOST_GPU_ENABLED typename Distribution::value_type coefficient_of_variation(const Distribution& dist) { typedef typename Distribution::value_type value_type; typedef typename Distribution::policy_type policy_type; @@ -104,19 +104,19 @@ inline typename Distribution::value_type coefficient_of_variation(const Distribu // implementation with all arguments of the same type: // template -inline typename Distribution::value_type pdf(const Distribution& dist, const RealType& x) +inline BOOST_GPU_ENABLED typename Distribution::value_type pdf(const Distribution& dist, const RealType& x) { typedef typename Distribution::value_type value_type; return pdf(dist, static_cast(x)); } template -inline typename Distribution::value_type cdf(const Distribution& dist, const RealType& x) +inline BOOST_GPU_ENABLED typename Distribution::value_type cdf(const Distribution& dist, const RealType& x) { typedef typename Distribution::value_type value_type; return cdf(dist, static_cast(x)); } template -inline typename Distribution::value_type quantile(const Distribution& dist, const RealType& x) +inline BOOST_GPU_ENABLED typename Distribution::value_type quantile(const Distribution& dist, const RealType& x) { typedef typename Distribution::value_type value_type; return quantile(dist, static_cast(x)); @@ -130,21 +130,21 @@ inline typename Distribution::value_type chf(const Distribution& dist, const Rea } */ template -inline typename Distribution::value_type cdf(const complemented2_type& c) +inline BOOST_GPU_ENABLED typename Distribution::value_type cdf(const complemented2_type& c) { typedef typename Distribution::value_type value_type; return cdf(complement(c.dist, static_cast(c.param))); } template -inline typename Distribution::value_type quantile(const complemented2_type& c) +inline BOOST_GPU_ENABLED typename Distribution::value_type quantile(const complemented2_type& c) { typedef typename Distribution::value_type value_type; return quantile(complement(c.dist, static_cast(c.param))); } template -inline typename Dist::value_type median(const Dist& d) +inline BOOST_GPU_ENABLED typename Dist::value_type median(const Dist& d) { // median - default definition for those distributions for which a // simple closed form is not known, // and for which a domain_error and/or NaN generating function is NOT defined. diff --git a/test/cuda/arcsine_cdf_double.cu b/test/cuda/arcsine_cdf_double.cu new file mode 100644 index 000000000..21090130e --- /dev/null +++ b/test/cuda/arcsine_cdf_double.cu @@ -0,0 +1,108 @@ +// 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 "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::arcsine_distribution(), 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); + 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::arcsine_distribution(), 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; +} + diff --git a/test/cuda/arcsine_cdf_float.cu b/test/cuda/arcsine_cdf_float.cu new file mode 100644 index 000000000..a534a4770 --- /dev/null +++ b/test/cuda/arcsine_cdf_float.cu @@ -0,0 +1,108 @@ +// 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 "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::arcsine_distribution(), 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); + 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::arcsine_distribution(), 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; +} + diff --git a/test/cuda/arcsine_pdf_double.cu b/test/cuda/arcsine_pdf_double.cu new file mode 100644 index 000000000..c37ca3c99 --- /dev/null +++ b/test/cuda/arcsine_pdf_double.cu @@ -0,0 +1,108 @@ +// 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 "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::arcsine_distribution(), 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); + 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::arcsine_distribution(), 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; +} + diff --git a/test/cuda/arcsine_pdf_float.cu b/test/cuda/arcsine_pdf_float.cu new file mode 100644 index 000000000..e4042aa40 --- /dev/null +++ b/test/cuda/arcsine_pdf_float.cu @@ -0,0 +1,108 @@ +// 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 "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::arcsine_distribution(), 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); + 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::arcsine_distribution(), 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; +} + diff --git a/test/cuda/arcsine_quan_double.cu b/test/cuda/arcsine_quan_double.cu new file mode 100644 index 000000000..f95d9c1ea --- /dev/null +++ b/test/cuda/arcsine_quan_double.cu @@ -0,0 +1,108 @@ +// 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 "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::arcsine_distribution(), 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); + 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::arcsine_distribution(), 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; +} + diff --git a/test/cuda/arcsine_quan_float.cu b/test/cuda/arcsine_quan_float.cu new file mode 100644 index 000000000..a36ac0d65 --- /dev/null +++ b/test/cuda/arcsine_quan_float.cu @@ -0,0 +1,108 @@ +// 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 "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::arcsine_distribution(), 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); + 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::arcsine_distribution(), 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; +} + diff --git a/test/cuda/cauchy_cdf_double.cu b/test/cuda/cauchy_cdf_double.cu new file mode 100644 index 000000000..065fd3df0 --- /dev/null +++ b/test/cuda/cauchy_cdf_double.cu @@ -0,0 +1,108 @@ +// 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 "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::cauchy_distribution(), 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(-10000, 10000); + // 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); + 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::cauchy_distribution(), 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; +} + diff --git a/test/cuda/cauchy_cdf_float.cu b/test/cuda/cauchy_cdf_float.cu new file mode 100644 index 000000000..f93a6c5f8 --- /dev/null +++ b/test/cuda/cauchy_cdf_float.cu @@ -0,0 +1,108 @@ +// 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 "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::cauchy_distribution(), 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(-10000, 10000); + // 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); + 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::cauchy_distribution(), 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; +} + diff --git a/test/cuda/cauchy_pdf_double.cu b/test/cuda/cauchy_pdf_double.cu new file mode 100644 index 000000000..aa6e1ce4c --- /dev/null +++ b/test/cuda/cauchy_pdf_double.cu @@ -0,0 +1,108 @@ +// 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 "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::cauchy_distribution(), 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(-10000, 10000); + // 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); + 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::cauchy_distribution(), 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; +} + diff --git a/test/cuda/cauchy_pdf_float.cu b/test/cuda/cauchy_pdf_float.cu new file mode 100644 index 000000000..88ed9b84b --- /dev/null +++ b/test/cuda/cauchy_pdf_float.cu @@ -0,0 +1,108 @@ +// 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 "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::cauchy_distribution(), 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(-10000, 10000); + // 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); + 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::cauchy_distribution(), 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; +} + diff --git a/test/cuda/cauchy_quan_double.cu b/test/cuda/cauchy_quan_double.cu new file mode 100644 index 000000000..96ed2910a --- /dev/null +++ b/test/cuda/cauchy_quan_double.cu @@ -0,0 +1,108 @@ +// 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 "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::cauchy_distribution(), 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); + 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::cauchy_distribution(), 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; +} + diff --git a/test/cuda/cauchy_quan_float.cu b/test/cuda/cauchy_quan_float.cu new file mode 100644 index 000000000..51bfd06bc --- /dev/null +++ b/test/cuda/cauchy_quan_float.cu @@ -0,0 +1,108 @@ +// 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 "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::cauchy_distribution(), 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); + 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::cauchy_distribution(), 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; +} +