1 //------------------------------------------------------------------------------
2 // GraphBLAS/Demo/Source/prand: parallel random number generator
3 //------------------------------------------------------------------------------
4
5 // SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved.
6 // SPDX-License-Identifier: Apache-2.0
7
8 //------------------------------------------------------------------------------
9
10 // A simple thread-safe parallel pseudo-random nuumber generator.
11
12 #include "GraphBLAS.h"
13 #undef GB_PUBLIC
14 #define GB_LIBRARY
15 #include "graphblas_demos.h"
16
17 //------------------------------------------------------------------------------
18 // prand macros
19 //------------------------------------------------------------------------------
20
21 // Generate the next seed, and extract a random 15-bit value from a seed.
22
23 #define PRAND_RECURENCE(seed) ((seed) * 1103515245 + 12345)
24
25 #define PRAND_15_MAX 32767
26 #define PRAND_15(seed) (((seed)/65536) % (PRAND_15_MAX + 1))
27
28 //------------------------------------------------------------------------------
29 // global types and operators
30 //------------------------------------------------------------------------------
31
32 // These can be shared by all threads in a user application, and thus are
33 // safely declared as global objects.
34
35 GrB_Type prand_type = NULL ;
36 GrB_UnaryOp prand_next_op = NULL ;
37 GrB_UnaryOp prand_iget_op = NULL ;
38 GrB_UnaryOp prand_xget_op = NULL ;
39 GrB_BinaryOp prand_dup_op = NULL ;
40
41 //------------------------------------------------------------------------------
42 // prand_next_op: unary operator to construct the next seed
43 //------------------------------------------------------------------------------
44
45 // z = f(x), where x is the old seed and z is the new seed.
46
47 GB_PUBLIC
prand_next_f(prand_t * z,const prand_t * x)48 void prand_next_f (prand_t *z, const prand_t *x)
49 {
50 for (int k = 0 ; k < 5 ; k++)
51 {
52 z->seed [k] = PRAND_RECURENCE (x->seed [k]) ;
53 }
54 }
55
56 //------------------------------------------------------------------------------
57 // prand_iget: unary operator to construct get a random integer from the seed
58 //------------------------------------------------------------------------------
59
60 // z = f(x), where x is a random seed, and z is an unsigned 64-bit
61 // pseudo-random number constructed from the seed.
62
63 GB_PUBLIC
prand_iget_f(uint64_t * z,const prand_t * x)64 void prand_iget_f (uint64_t *z, const prand_t *x)
65 {
66 uint64_t i = 0 ;
67 for (int k = 0 ; k < 5 ; k++)
68 {
69 i = PRAND_15_MAX * i + PRAND_15 (x->seed [k]) ;
70 }
71 (*z) = i ;
72 }
73
74 //------------------------------------------------------------------------------
75 // prand_xget: unary operator to construct get a random double from the seed
76 //------------------------------------------------------------------------------
77
78 // z = f(x), where x is a random seed, and z is a double precision
79 // pseudo-random number constructed from the seed, in the range 0 to 1.
80
81 GB_PUBLIC
prand_xget_f(double * z,prand_t * x)82 void prand_xget_f (double *z, prand_t *x)
83 {
84 uint64_t i ;
85 prand_iget_f (&i, x) ;
86 (*z) = ((double) i) / ((double) UINT64_MAX) ;
87 }
88
89 //------------------------------------------------------------------------------
90 // prand_dup: binary operator to build a vector
91 //------------------------------------------------------------------------------
92
93 // This is required by GrB_Vector_build, but is never called since no
94 // duplicates are created. This is the SECOND operator for the prand_type.
95
96 #if defined ( __INTEL_COMPILER )
97 // disable icc warnings
98 // 869: unused parameters
99 #pragma warning (disable: 869 )
100 #elif defined ( __GNUC__ )
101 #pragma GCC diagnostic ignored "-Wunused-parameter"
102 #endif
103
104 GB_PUBLIC
prand_dup_f(prand_t * z,const prand_t * x,const prand_t * y)105 void prand_dup_f (prand_t *z, /* unused: */ const prand_t *x, const prand_t *y)
106 {
107 (*z) = (*y) ;
108 }
109
110 //------------------------------------------------------------------------------
111 // prand_init: create the random seed type and its operators
112 //------------------------------------------------------------------------------
113
114 #define PRAND_FREE_ALL \
115 { \
116 GrB_Type_free (&prand_type) ; \
117 GrB_UnaryOp_free (&prand_next_op) ; \
118 GrB_UnaryOp_free (&prand_iget_op) ; \
119 GrB_UnaryOp_free (&prand_xget_op) ; \
120 GrB_BinaryOp_free (&prand_dup_op) ; \
121 }
122
123 #undef OK
124 #define OK(method) \
125 { \
126 GrB_Info info = method ; \
127 if (info != GrB_SUCCESS) \
128 { \
129 PRAND_FREE_ALL ; \
130 printf ("GraphBLAS error: %d\n", info) ; \
131 return (info) ; \
132 } \
133 }
134
135 GB_PUBLIC
prand_init()136 GrB_Info prand_init ( )
137 {
138 prand_type = NULL ;
139 prand_next_op = NULL ;
140 prand_iget_op = NULL ;
141 prand_xget_op = NULL ;
142 prand_dup_op = NULL ;
143 OK (GrB_Type_new (&prand_type, sizeof (prand_t))) ;
144 OK (GrB_UnaryOp_new (&prand_next_op, (GxB_unary_function) prand_next_f,
145 prand_type, prand_type)) ;
146 OK (GrB_UnaryOp_new (&prand_iget_op, (GxB_unary_function) prand_iget_f,
147 GrB_UINT64, prand_type)) ;
148 OK (GrB_UnaryOp_new (&prand_xget_op, (GxB_unary_function) prand_xget_f,
149 GrB_FP64, prand_type)) ;
150 OK (GrB_BinaryOp_new (&prand_dup_op, (GxB_binary_function) prand_dup_f,
151 prand_type, prand_type, prand_type)) ;
152 return (GrB_SUCCESS) ;
153 }
154
155 //------------------------------------------------------------------------------
156 // prand_finalize: free the random seed type and its operators
157 //------------------------------------------------------------------------------
158
159 GB_PUBLIC
prand_finalize()160 GrB_Info prand_finalize ( )
161 {
162 PRAND_FREE_ALL ;
163 return (GrB_SUCCESS) ;
164 }
165
166 //------------------------------------------------------------------------------
167 // prand_next: get the next random numbers
168 //------------------------------------------------------------------------------
169
170 GB_PUBLIC
prand_next(GrB_Vector Seed)171 GrB_Info prand_next
172 (
173 GrB_Vector Seed
174 )
175 {
176 return (GrB_Vector_apply (Seed, NULL, NULL, prand_next_op, Seed, NULL)) ;
177 }
178
179 //------------------------------------------------------------------------------
180 // prand_seed: create a vector of random seeds
181 //------------------------------------------------------------------------------
182
183 // Returns a vector of random seed values.
184
185 #define PRAND_FREE_WORK \
186 { \
187 free (I) ; \
188 free (X) ; \
189 }
190
191 #undef PRAND_FREE_ALL
192 #define PRAND_FREE_ALL \
193 { \
194 PRAND_FREE_WORK ; \
195 GrB_Vector_free (Seed) ; \
196 }
197
198 GB_PUBLIC
prand_seed(GrB_Vector * Seed,int64_t seed,GrB_Index n,int nthreads)199 GrB_Info prand_seed
200 (
201 GrB_Vector *Seed, // vector of random number seeds
202 int64_t seed, // scalar input seed
203 GrB_Index n, // size of Seed to create
204 int nthreads // # of threads to use (OpenMP default if <= 0)
205 )
206 {
207
208 GrB_Index *I = NULL ;
209 prand_t *X = NULL ;
210
211 // allocate the Seed vector
212 OK (GrB_Vector_new (Seed, prand_type, n)) ;
213
214 // allocate the I and X arrays
215 I = (GrB_Index *) malloc ((n+1) * sizeof (GrB_Index)) ;
216 X = (prand_t *) malloc ((n+1) * sizeof (prand_t)) ;
217 if (I == NULL || X == NULL)
218 {
219 PRAND_FREE_ALL ;
220 return (GrB_OUT_OF_MEMORY) ;
221 }
222
223 // determine # of threads to use
224 int nthreads_max = 1 ;
225 #ifdef _OPENMP
226 nthreads_max = omp_get_max_threads ( ) ;
227 #endif
228 if (nthreads <= 0 || nthreads > nthreads_max)
229 {
230 nthreads = nthreads_max ;
231 }
232
233 // construct the tuples for the initial seeds
234 int64_t i, len = (int64_t) n ;
235 #pragma omp parallel for num_threads(nthreads) schedule(static)
236 for (i = 0 ; i < len ; i++)
237 {
238 I [i] = i ;
239 for (int k = 0 ; k < 5 ; k++)
240 {
241 X [i].seed [k] = (100000000*(seed) + 10*i + k + 1) ;
242 }
243 }
244
245 // build the Seed vector
246 OK (GrB_Vector_build_UDT (*Seed, I, X, n, prand_dup_op)) ;
247
248 // free workspace
249 PRAND_FREE_WORK ;
250
251 // advance to the first set of random numbers
252 OK (prand_next (*Seed)) ;
253
254 return (GrB_SUCCESS) ;
255 }
256
257 //------------------------------------------------------------------------------
258 // prand_print: print the Seed vector
259 //------------------------------------------------------------------------------
260
261 // This is meant for testing, not production use.
262
263 #undef PRAND_FREE_ALL
264 #define PRAND_FREE_ALL ;
265
266 GB_PUBLIC
prand_print(GrB_Vector Seed,int pr)267 GrB_Info prand_print
268 (
269 GrB_Vector Seed,
270 int pr // 0: print nothing, 1: print some, 2: print all
271 )
272 {
273 if (pr > 0)
274 {
275 GrB_Index n ;
276 OK (GrB_Vector_nvals (&n, Seed)) ;
277 printf ("\nSeed: length %g\n", (double) n) ;
278 prand_t x ;
279 for (int k = 0 ; k < 5 ; k++) x.seed [k] = -1 ;
280 for (int64_t i = 0 ; i < (int64_t) n ; i++)
281 {
282 if (GrB_Vector_extractElement_UDT (&x, Seed, i) == GrB_SUCCESS)
283 {
284 printf ("%g: ", (double) i) ;
285 for (int k = 0 ; k < 5 ; k++)
286 {
287 printf (" %.18g", (double) (x.seed [k])) ;
288 }
289 printf ("\n") ;
290 }
291 if (pr == 1 && i > 10) break ;
292 }
293 }
294 return (GrB_SUCCESS) ;
295 }
296
297 //------------------------------------------------------------------------------
298 // prand_iget: return a vector of random uint64 integers
299 //------------------------------------------------------------------------------
300
301 GB_PUBLIC
prand_iget(GrB_Vector X,GrB_Vector Seed)302 GrB_Info prand_iget
303 (
304 GrB_Vector X,
305 GrB_Vector Seed
306 )
307 {
308 OK (GrB_Vector_apply (X, NULL, NULL, prand_iget_op, Seed, NULL)) ;
309 return (prand_next (Seed)) ;
310 }
311
312 //------------------------------------------------------------------------------
313 // prand_xget: return a vector of random doubles, in range 0 to 1 inclusive
314 //------------------------------------------------------------------------------
315
316 GB_PUBLIC
prand_xget(GrB_Vector X,GrB_Vector Seed)317 GrB_Info prand_xget
318 (
319 GrB_Vector X,
320 GrB_Vector Seed
321 )
322 {
323 OK (GrB_Vector_apply (X, NULL, NULL, prand_xget_op, Seed, NULL)) ;
324 return (prand_next (Seed)) ;
325 }
326
327