Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/RNG/MKLRNG.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Interface to Intel MKL VSL random number generators
9 : : \details Interface to Intel MKL VSL random number generators.
10 : : */
11 : : // *****************************************************************************
12 : : #ifndef MKLRNG_h
13 : : #define MKLRNG_h
14 : :
15 : : #include <memory>
16 : : #include <mkl_vsl.h>
17 : :
18 : : #include "Exception.hpp"
19 : : #include "Keywords.hpp"
20 : :
21 : : namespace tk {
22 : :
23 : : //! MKL-based random number generator used polymorphically with tk::RNG
24 : : class MKLRNG {
25 : :
26 : : using ncomp_t = kw::ncomp::info::expect::type;
27 : :
28 : : public:
29 : : //! Constructor
30 : : //! \param[in] n Initialize RNG using this many independent streams
31 : : //! \param[in] brng Index of the basic generator to initialize the stream
32 : : //! \param[in] seed RNG seed
33 : : //! \param[in] uniform_method MKL ID of the method to use for uniform RNGs
34 : : //! \param[in] gaussian_method MKL ID of the method to use for Gaussian RNGs
35 : : //! \param[in] gaussianmv_method MKL ID of the method to use for
36 : : //! multi-variate Gaussian RNGs
37 : : //! \param[in] beta_method MKL ID of the method to use for beta RNGs
38 : : //! \param[in] gamma_method MKL ID of the method to use for gamma RNGs
39 : 174 : explicit MKLRNG(
40 : : int n = 1,
41 : : int brng = VSL_BRNG_MCG59,
42 : : unsigned int seed = 0,
43 : : int uniform_method = VSL_RNG_METHOD_UNIFORM_STD,
44 : : int gaussian_method = VSL_RNG_METHOD_GAUSSIAN_BOXMULLER,
45 : : int gaussianmv_method = VSL_RNG_METHOD_GAUSSIANMV_BOXMULLER2,
46 : : int beta_method = VSL_RNG_METHOD_BETA_CJA,
47 : 174 : int gamma_method = VSL_RNG_METHOD_GAMMA_GNORM ) :
48 : : m_brng( brng ),
49 : : m_seed( seed ),
50 : : m_uniform_method( uniform_method ),
51 : : m_gaussian_method( gaussian_method ),
52 : : m_gaussianmv_method( gaussianmv_method ),
53 : : m_beta_method( beta_method ),
54 : : m_gamma_method( gamma_method ),
55 : : m_nthreads( n ),
56 [ + - ]: 174 : m_stream()
57 : : {
58 : : Assert( n > 0, "Need at least one thread" );
59 : : Assert( brng > 0, "Basic RNG MKL parameter must be positive" );
60 : : // Allocate array of stream-pointers for threads
61 : 174 : m_stream = std::make_unique< VSLStreamStatePtr[] >(
62 [ + - ]: 174 : static_cast<std::size_t>(n) );
63 : : // Initialize thread-streams for block-splitting. These MKL VSL functions
64 : : // dynamically allocate memory, so these calls being in a constructor are
65 : : // a potential memory leak hazard in the presence of exceptions. However,
66 : : // thankfully, the MKL functions below only emit warnings if they
67 : : // encounter errors and always continue. As a result, the constructor
68 : : // finishes, the MKLRNG object gets created, so the destructor will also
69 : : // get called when leaving scope.
70 [ + + ]: 174 : if (n == 1)
71 [ + - ][ + - ]: 7 : errchk( vslNewStream( &m_stream[0], brng, seed ) );
72 : : else
73 [ + + ]: 815 : for (int i=0; i<n; ++i) {
74 : 649 : auto I = static_cast< std::size_t >( i );
75 [ + - ][ + - ]: 649 : errchk( vslNewStream( &m_stream[I], brng, seed ) );
76 [ + - ][ + + ]: 649 : errchk( vslLeapfrogStream( m_stream[I], i, n ) );
77 : : }
78 : 173 : }
79 : :
80 : : //! Destructor
81 : 580 : ~MKLRNG() noexcept { deleteStreams(); }
82 : :
83 : : //! Uniform RNG: Generate uniform random numbers
84 : : //! \param[in] tid Thread (or more precisely stream) ID
85 : : //! \param[in] num Number of RNGs to generate
86 : : //! \param[in,out] r Pointer to memory to write the random numbers to
87 : : void uniform( int tid, ncomp_t num, double* r ) const {
88 : 907326242 : vdRngUniform( m_uniform_method,
89 [ + - ]: 907326242 : m_stream[ static_cast<std::size_t>(tid) ],
90 : : static_cast< long long >( num ),
91 : : r,
92 : : 0.0, 1.0 );
93 : 8 : }
94 : :
95 : : //! Gaussian RNG: Generate Gaussian random numbers
96 : : //! \param[in] tid Thread (or more precisely stream) ID
97 : : //! \param[in] num Number of RNGs to generate
98 : : //! \param[in,out] r Pointer to memory to write the random numbers to
99 : : void gaussian( int tid, ncomp_t num, double* r ) const {
100 : 442000089 : vdRngGaussian( m_gaussian_method,
101 [ + - ]: 442000089 : m_stream[ static_cast<std::size_t>(tid) ],
102 : : static_cast< long long >( num ),
103 : : r,
104 : : 0.0, 1.0 );
105 : 61 : }
106 : :
107 : : //! \brief Multi-variate Gaussian RNG: Generate multi-variate Gaussian
108 : : //! random numbers
109 : : //! \param[in] tid Thread (or more precisely stream) ID
110 : : //! \param[in] num Number of RNGs to generate
111 : : //! \param[in] d Dimension d ( d ≥ 1) of output random vectors
112 : : //! \param[in] mean Mean vector of dimension d
113 : : //! \param[in] cov Lower triangle of covariance matrix, stored as a vector
114 : : //! of length d(d+1)/2
115 : : //! \param[in,out] r Pointer to memory to write the random numbers to
116 : : void gaussianmv( int tid, ncomp_t num, ncomp_t d, const double* const mean,
117 : : const double* const cov, double* r ) const
118 : : {
119 : : Assert( d > 0,
120 : : "Dimension of multi-variate Gaussian RNGs must be positive" );
121 : 16 : vdRngGaussianMV( m_gaussianmv_method,
122 [ + - ][ + - ]: 16 : m_stream[ static_cast<std::size_t>(tid) ],
123 : : static_cast< long long >( num ),
124 : : r,
125 : : static_cast< int >( d ),
126 : : VSL_MATRIX_STORAGE_PACKED,
127 : : mean,
128 : : cov );
129 : 16 : }
130 : :
131 : : //! Beta RNG: Generate beta random numbers
132 : : //! \param[in] tid Thread (or more precisely stream) ID
133 : : //! \param[in] num Number of RNGs to generate
134 : : //! \param[in] p First beta shape parameter
135 : : //! \param[in] q Second beta shape parameter
136 : : //! \param[in] a Beta displacement parameter
137 : : //! \param[in] b Beta scale factor
138 : : //! \param[in,out] r Pointer to memory to write the random numbers to
139 : : void beta( int tid, ncomp_t num, double p, double q, double a, double b,
140 : : double* r ) const
141 : : {
142 : 32 : vdRngBeta( m_beta_method,
143 [ + - ][ + - ]: 32 : m_stream[ static_cast<std::size_t>(tid) ],
144 : : static_cast< long long >( num ),
145 : : r,
146 : : p, q, a, b );
147 : 24 : }
148 : :
149 : : //! Gamma RNG: Generate gamma random numbers
150 : : //! \param[in] tid Thread (or more precisely stream) ID
151 : : //! \param[in] num Number of RNGs to generate
152 : : //! \param[in] a Gamma shape parameter
153 : : //! \param[in] b Gamma scale factor
154 : : //! \param[in,out] r Pointer to memory to write the random numbers to
155 : : void gamma( int tid, ncomp_t num, double a, double b, double* r ) const
156 : : {
157 : 0 : vdRngGamma( m_beta_method,
158 : 0 : m_stream[ static_cast<std::size_t>(tid) ],
159 : : static_cast< long long >( num ),
160 : : r,
161 : : a, 0.0, b ); // displacement = 0.0
162 : : }
163 : :
164 : : //! Copy assignment
165 : 19 : MKLRNG& operator=( const MKLRNG& x ) {
166 : 19 : m_brng = x.m_brng;
167 : 19 : m_seed = x.m_seed;
168 : 19 : m_uniform_method = x.m_uniform_method;
169 : 19 : m_gaussian_method = x.m_gaussian_method;
170 : 19 : m_gaussianmv_method = x.m_gaussianmv_method;
171 : 19 : m_beta_method = x.m_beta_method;
172 : 19 : m_gamma_method = x.m_gamma_method;
173 : 19 : m_nthreads = x.m_nthreads;
174 : 19 : m_stream = std::make_unique< VSLStreamStatePtr[] >(
175 : 19 : static_cast<std::size_t>(x.m_nthreads) );
176 [ + + ]: 19 : if (m_nthreads == 1)
177 : 4 : errchk( vslNewStream( &m_stream[0], x.m_brng, x.m_seed ) );
178 : : else
179 [ + + ]: 75 : for (int i=0; i<x.m_nthreads; ++i) {
180 : 60 : auto I = static_cast< std::size_t >( i );
181 : 60 : errchk( vslNewStream( &m_stream[I], x.m_brng, x.m_seed ) );
182 : 60 : errchk( vslLeapfrogStream( m_stream[I], i, x.m_nthreads ) );
183 : : }
184 : 19 : return *this;
185 : : }
186 : :
187 : : //! Copy constructor: in terms of copy assignment
188 [ + - ]: 16 : MKLRNG( const MKLRNG& x ) { operator=(x); }
189 : :
190 : : //! Move assignment
191 : 391 : MKLRNG& operator=( MKLRNG&& x ) {
192 : 391 : deleteStreams();
193 : 391 : m_brng = x.m_brng;
194 : 391 : m_seed = x.m_seed;
195 : 391 : m_uniform_method = x.m_uniform_method;
196 : 391 : m_gaussian_method = x.m_gaussian_method;
197 : 391 : m_gaussianmv_method = x.m_gaussianmv_method;
198 : 391 : m_beta_method = x.m_beta_method;
199 : 391 : m_gamma_method = x.m_gamma_method;
200 : 391 : m_nthreads = x.m_nthreads;
201 : 391 : m_stream = std::make_unique< VSLStreamStatePtr[] >(
202 : 391 : static_cast<std::size_t>(x.m_nthreads) );
203 [ + + ]: 1896 : for (int i=0; i<x.m_nthreads; ++i) {
204 : 1505 : auto I = static_cast< std::size_t >( i );
205 : 1505 : m_stream[I] = x.m_stream[I];
206 : 1505 : x.m_stream[I] = nullptr;
207 : : }
208 : 391 : x.m_brng = 0;
209 : 391 : x.m_seed = 0;
210 : 391 : x.m_uniform_method = 0;
211 : 391 : x.m_gaussian_method = 0;
212 : 391 : x.m_gaussianmv_method = 0;
213 : 391 : x.m_beta_method = 0;
214 : 391 : x.m_gamma_method = 0;
215 [ + - ]: 391 : x.m_nthreads = 0;
216 : : x.m_stream.reset( nullptr );
217 : 391 : return *this;
218 : : }
219 : :
220 : : //! Move constructor: in terms of move assignment
221 : 391 : MKLRNG( MKLRNG&& x ) :
222 : : m_brng( 0 ),
223 : : m_seed( 0 ),
224 : : m_uniform_method( 0 ),
225 : : m_gaussian_method( 0 ),
226 : : m_gaussianmv_method( 0 ),
227 : : m_beta_method( 0 ),
228 : : m_gamma_method( 0 ),
229 : : m_nthreads( 0 ),
230 [ + - ]: 391 : m_stream( nullptr )
231 [ + - ]: 391 : { *this = std::move(x); }
232 : :
233 : : //! Accessor to the number of threads we operate on
234 : : std::size_t nthreads() const noexcept
235 : 39 : { return static_cast< std::size_t >( m_nthreads); }
236 : :
237 : : private:
238 : : //! Delete all thread streams
239 : 971 : void deleteStreams() {
240 [ + + ]: 1681 : for (int i=0; i<m_nthreads; ++i) {
241 [ + - ]: 710 : auto I = static_cast< std::size_t >( i );
242 [ + - ]: 710 : if (m_stream[I]) {
243 : 710 : vslDeleteStream( &m_stream[I] );
244 : 710 : m_stream[I] = nullptr;
245 : : }
246 : : }
247 : 971 : }
248 : :
249 : : //! MKL VSL error check
250 : : //! \param[in] err MKL VSL error code as returned from MKL VSL functions
251 : : //! \details This calls ErrChk(), i.e., it is not compiled away in Release
252 : : //! mode as an error here can result due to user input incompatible with
253 : : //! the MKL library.
254 : 1429 : void errchk( int err ) {
255 [ + + ][ + - ]: 1430 : ErrChk( err == VSL_STATUS_OK, "MKL VSL Error Code: " +
[ + - ][ + - ]
[ + - ][ + - ]
[ - + ][ - + ]
[ - + ][ - + ]
[ - - ][ - - ]
[ - - ]
256 : : std::to_string(err) + ", see mkl_vsl_defines.h for more info" );
257 : 1428 : }
258 : :
259 : : int m_brng; //!< MKL RNG id
260 : : unsigned int m_seed; //!< Seed
261 : : int m_uniform_method; //!< Uniform method to use
262 : : int m_gaussian_method; //!< Gaussian method to use
263 : : int m_gaussianmv_method; //!< Multi-variate Gaussian method to use
264 : : int m_beta_method; //!< Beta method to use
265 : : int m_gamma_method; //!< Gamma method to use
266 : : int m_nthreads; //!< Number of threads
267 : : std::unique_ptr< VSLStreamStatePtr[] > m_stream; //!< Random number streams
268 : : };
269 : :
270 : : } // tk::
271 : :
272 : : #endif // MKLRNG_h
|
