1 // =============================================================================
2 // PROJECT CHRONO - http://projectchrono.org
3 //
4 // Copyright (c) 2016 projectchrono.org
5 // All rights reserved.
6 //
7 // Use of this source code is governed by a BSD-style license that can be found
8 // in the LICENSE file at the top level of the distribution and at
9 // http://projectchrono.org/license-chrono.txt.
10 //
11 // =============================================================================
12 // Authors: Hammad Mazhar, Radu Serban
13 // =============================================================================
14 //
15 // Description: definition of some convenience functions for math operations
16 // =============================================================================
17
18 #pragma once
19
20 #include "chrono/multicore_math/real.h"
21 #include "chrono/multicore_math/real2.h"
22 #include "chrono/multicore_math/real3.h"
23 #include "chrono/multicore_math/real4.h"
24 #include "chrono/multicore_math/other_types.h"
25 #include "chrono/multicore_math/utility.h"
26
27 #undef _GLIBCXX_ATOMIC_BUILTINS
28 #undef _GLIBCXX_USE_INT128
29
30 #include <iostream>
31
32 // -----------------------------------------------------------------------------
33 // Thrust related defines
34 // -----------------------------------------------------------------------------
35
36 // Always include ChConfig.h *before* any Thrust headers!
37 #include "chrono/ChConfig.h"
38 #include <thrust/reduce.h>
39 #include <thrust/gather.h>
40 #include <thrust/scan.h>
41 #include <thrust/fill.h>
42 #include <thrust/copy.h>
43 #include <thrust/iterator/counting_iterator.h>
44
45 #if defined(CHRONO_OPENMP_ENABLED)
46 #include <thrust/system/omp/execution_policy.h>
47 #elif defined(CHRONO_TBB_ENABLED)
48 #include <thrust/system/tbb/execution_policy.h>
49 #endif
50
51 #ifndef _MSC_VER
52 #include <cfenv>
53 #endif
54
55 namespace chrono {
56
57 /// @addtogroup chrono_mc_math
58 /// @{
59
60 #if defined(CHRONO_OPENMP_ENABLED)
61 #define THRUST_PAR thrust::omp::par,
62 #elif defined(CHRONO_TBB_ENABLED)
63 #define THRUST_PAR thrust::tbb::par,
64 #else
65 #define THRUST_PAR
66 #endif
67
68 #define Thrust_Inclusive_Scan_Sum(x, y) \
69 thrust::inclusive_scan(THRUST_PAR x.begin(), x.end(), x.begin()); \
70 y = x.back();
71 #define Thrust_Sort_By_Key(x, y) thrust::sort_by_key(THRUST_PAR x.begin(), x.end(), y.begin())
72
73 #define Run_Length_Encode(y, z, w) \
74 (thrust::reduce_by_key(THRUST_PAR y.begin(), y.end(), thrust::constant_iterator<uint>(1), z.begin(), w.begin()) \
75 .second) - \
76 w.begin()
77
78 #define Thrust_Inclusive_Scan(x) thrust::inclusive_scan(THRUST_PAR x.begin(), x.end(), x.begin())
79 #define Thrust_Exclusive_Scan(x) thrust::exclusive_scan(THRUST_PAR x.begin(), x.end(), x.begin())
80 #define Thrust_Fill(x, y) thrust::fill(x.begin(), x.end(), y)
81 #define Thrust_Sort(x) thrust::sort(THRUST_PAR x.begin(), x.end())
82 #define Thrust_Count(x, y) thrust::count(THRUST_PAR x.begin(), x.end(), y)
83 #define Thrust_Sequence(x) thrust::sequence(x.begin(), x.end())
84 #define Thrust_Equal(x, y) thrust::equal(THRUST_PAR x.begin(), x.end(), y.begin())
85 #define Thrust_Max(x) x[thrust::max_element(THRUST_PAR x.begin(), x.end()) - x.begin()]
86 #define Thrust_Min(x) x[thrust::min_element(THRUST_PAR x.begin(), x.end()) - x.begin()]
87 #define Thrust_Total(x) thrust::reduce(THRUST_PAR x.begin(), x.end())
88 #define Thrust_Unique(x) thrust::unique(THRUST_PAR x.begin(), x.end()) - x.begin();
89 #define DBG(x) printf(x);
90
91 /// Explicit conversion of scoped enumeration to int (e.g. for streaming).
92 template <typename Enumeration>
93 auto as_integer(Enumeration const value) -> typename std::underlying_type<Enumeration>::type {
94 return static_cast<typename std::underlying_type<Enumeration>::type>(value);
95 }
96
97 /// Utility to expand an input sequence by
98 /// replicating each element a variable number of times. For example,
99 /// - expand([2,2,2],[A,B,C]) -> [A,A,B,B,C,C]
100 /// - expand([3,0,1],[A,B,C]) -> [A,A,A,C]
101 /// - expand([1,3,2],[A,B,C]) -> [A,B,B,B,C,C]
102 ///
103 /// The element counts are assumed to be non-negative integers
104 /// (from Thrust's example expand.cu)
105 template <typename InputIterator1, typename InputIterator2, typename OutputIterator>
Thrust_Expand(InputIterator1 first1,InputIterator1 last1,InputIterator2 first2,OutputIterator output)106 OutputIterator Thrust_Expand(InputIterator1 first1,
107 InputIterator1 last1,
108 InputIterator2 first2,
109 OutputIterator output) {
110 typedef typename thrust::iterator_difference<InputIterator1>::type difference_type;
111
112 difference_type input_size = thrust::distance(first1, last1);
113 difference_type output_size = thrust::reduce(THRUST_PAR first1, last1);
114
115 // scan the counts to obtain output offsets for each input element
116 std::vector<difference_type> output_offsets(input_size, 0);
117 thrust::exclusive_scan(THRUST_PAR first1, last1, output_offsets.begin());
118
119 // scatter the nonzero counts into their corresponding output positions
120 std::vector<difference_type> output_indices(output_size, 0);
121 thrust::scatter_if(THRUST_PAR thrust::counting_iterator<difference_type>(0),
122 thrust::counting_iterator<difference_type>(input_size), output_offsets.begin(), first1,
123 output_indices.begin());
124
125 // compute max-scan over the output indices, filling in the holes
126 thrust::inclusive_scan(THRUST_PAR output_indices.begin(), output_indices.end(), output_indices.begin(),
127 thrust::maximum<difference_type>());
128
129 // gather input values according to index array (output = first2[output_indices])
130 OutputIterator output_end = output;
131 thrust::advance(output_end, output_size);
132 thrust::gather(THRUST_PAR output_indices.begin(), output_indices.end(), first2, output);
133
134 // return output + output_size
135 thrust::advance(output, output_size);
136 return output;
137 }
138
139 /// @} chrono_mc_math
140
141 } // end namespace chrono
142