Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_Random: generate a random vector (of any sparsity structure)
3 : //------------------------------------------------------------------------------
4 :
5 : // LAGraph, (c) 2019-2025 by The LAGraph Contributors, All Rights Reserved.
6 : // SPDX-License-Identifier: BSD-2-Clause
7 : //
8 : // For additional details (including references to third party source code and
9 : // other files) see the LICENSE file or contact permission@sei.cmu.edu. See
10 : // Contributors.txt for a full list of contributors. Created, in part, with
11 : // funding and support from the U.S. Government (see Acknowledgments.txt file).
12 : // DM22-0790
13 :
14 : // Contributed by Timothy A. Davis, Texas A&M University
15 :
16 : //------------------------------------------------------------------------------
17 :
18 : // A very simple thread-safe parallel pseudo-random number generator. These
19 : // have a period of 2^(64)-1. For longer periods, a better method is needed.
20 :
21 : // In the future, the State vector could be based on a 256-bit user-defined
22 : // type, created by LG_Random_Init below. The LAGraph_Random_Seed and
23 : // LAGraph_Random_Next API would not change. Instead, they would select the
24 : // operators based on the type of the State vector (GrB_UINT64 in the current
25 : // method, or the 256-bit type in a future method). Then we would need a
26 : // single new method, say LAGraph_Random_Value, which extracts the random
27 : // number from the State. It would compute R = f(State), where R is GrB_UINT64
28 : // and the State vector (with its 256-bit data type) is unchanged. With the
29 : // current method, R=State is implicit.
30 :
31 : #include "LG_internal.h"
32 :
33 : //------------------------------------------------------------------------------
34 : // global operator
35 : //------------------------------------------------------------------------------
36 :
37 : // These operators can be shared by all threads in a user application, and thus
38 : // are safely declared as global objects.
39 :
40 : GrB_UnaryOp LG_rand_next_op = NULL ;
41 : GrB_IndexUnaryOp LG_rand_init_op = NULL ;
42 :
43 : //------------------------------------------------------------------------------
44 : // unary and index-unary ops to construct the first and next states
45 : //------------------------------------------------------------------------------
46 :
47 : // z = f(x), where x is the old state and z is the new state.
48 :
49 : // using xorshift64, from https://en.wikipedia.org/wiki/Xorshift
50 : // with a state of uint64_t.
51 :
52 : // Reference: Marsaglia, George (July 2003). "Xorshift RNGs". Journal of
53 : // Statistical Software. 8 (14). https://doi.org/10.18637/jss.v008.i14 .
54 :
55 : // For this random number generator, the output random number is the same
56 : // as the state.
57 :
58 : #if 0
59 :
60 : The default initial state is given below, but is unused here:
61 : #define LG_RAND_MARSAGLIA_SEED 88172645463325252LL
62 :
63 : struct xorshift64_state {
64 : uint64_t a;
65 : };
66 :
67 : // the state must be initialized to nonzero
68 : uint64_t xorshift64(struct xorshift64_state *state)
69 : {
70 : uint64_t x = state->a;
71 : x ^= x << 13;
72 : x ^= x >> 7;
73 : x ^= x << 17;
74 : return state->a = x;
75 : }
76 :
77 : #endif
78 :
79 : // return a random uint64_t; for internal use in LAGraph
80 7571279 : uint64_t LG_Random64 (uint64_t *state)
81 : {
82 7571279 : (*state) ^= (*state) << 13 ;
83 7571279 : (*state) ^= (*state) >> 7 ;
84 7571279 : (*state) ^= (*state) << 17 ;
85 7571279 : return (*state) ;
86 : }
87 :
88 : // return a random uint64_t; as a unary operator
89 40354361 : void LG_rand_next_f2 (uint64_t *z, const uint64_t *x)
90 : {
91 40354361 : uint64_t state = (*x) ;
92 40354361 : state ^= state << 13 ;
93 40354361 : state ^= state >> 7 ;
94 40354361 : state ^= state << 17 ;
95 40354361 : (*z) = state ;
96 40354361 : }
97 :
98 : #define LG_RAND_NEXT_F2_DEFN \
99 : "void LG_rand_next_f2 (uint64_t *z, const uint64_t *x) \n" \
100 : "{ \n" \
101 : " uint64_t state = (*x) ; \n" \
102 : " state ^= state << 13 ; \n" \
103 : " state ^= state >> 7 ; \n" \
104 : " state ^= state << 17 ; \n" \
105 : " (*z) = state ; \n" \
106 : "}"
107 :
108 : // From these references, the recommendation is to create the initial state of
109 : // a random number generator with an entirely different random number
110 : // generator. splitmix64 is recommended, so we initialize the State(i) with
111 : // splitmix64 (i+seed). The method cannot return a value of zero, so it is
112 : // suitable as a seed for the xorshift64 generator, above.
113 :
114 : // References:
115 : //
116 : // David Blackman and Sebastiano Vigna. Scrambled linear pseudorandom number
117 : // generators. ACM Trans. Math. Softw., 47:1−32, 2021.
118 : //
119 : // Steele GL, Vigna S. Computationally easy, spectrally good multipliers for
120 : // congruential pseudorandom number generators. Software: Practice and
121 : // Experience 2022; 52(2): 443–458. https://doi.org/10.1002/spe.3030
122 : //
123 : // Guy L. Steele, Doug Lea, and Christine H. Flood. 2014. Fast splittable
124 : // pseudorandom number generators. SIGPLAN Not. 49, 10 (October 2014), 453–472.
125 : // https://doi.org/10.1145/2714064.2660195
126 : //
127 : // The splitmix64 below method is the mix64variant13 in the above paper.
128 :
129 : #define GOLDEN_GAMMA 0x9E3779B97F4A7C15LL
130 :
131 : #if 0
132 :
133 : struct splitmix64_state {
134 : uint64_t s;
135 : };
136 :
137 : uint64_t splitmix64(struct splitmix64_state *state)
138 : {
139 : uint64_t result = (state->s += 0x9E3779B97F4A7C15LL) ;
140 : result = (result ^ (result >> 30)) * 0xBF58476D1CE4E5B9LL ;
141 : result = (result ^ (result >> 27)) * 0x94D049BB133111EBLL ;
142 : return result ^ (result >> 31);
143 : }
144 :
145 : #endif
146 :
147 : // The init function computes z = splitmix64 (i + seed), but it does not
148 : // advance the seed value on return.
149 3586158 : void LG_rand_init_func (uint64_t *z, const void *x,
150 : GrB_Index i, GrB_Index j, const uint64_t *seed)
151 : {
152 3586158 : uint64_t state = i + (*seed) ;
153 3586158 : uint64_t result = (state += GOLDEN_GAMMA) ;
154 3586158 : result = (result ^ (result >> 30)) * 0xBF58476D1CE4E5B9LL ;
155 3586158 : result = (result ^ (result >> 27)) * 0x94D049BB133111EBLL ;
156 3586158 : result = (result ^ (result >> 31)) ;
157 3586158 : (*z) = result ;
158 3586158 : }
159 :
160 : #define LG_RAND_INIT_F_DEFN \
161 : "void LG_rand_init_func (uint64_t *z, const void *x, \n" \
162 : " GrB_Index i, GrB_Index j, const uint64_t *seed) \n" \
163 : "{ \n" \
164 : " uint64_t state = i + (*seed) ; \n" \
165 : " uint64_t result = (state += 0x9E3779B97F4A7C15LL) ; \n" \
166 : " result = (result ^ (result >> 30)) * 0xBF58476D1CE4E5B9LL ; \n" \
167 : " result = (result ^ (result >> 27)) * 0x94D049BB133111EBLL ; \n" \
168 : " result = (result ^ (result >> 31)) ; \n" \
169 : " (*z) = result ; \n" \
170 : "}"
171 :
172 : //------------------------------------------------------------------------------
173 : // LG_Random_Init: create the random state operator
174 : //------------------------------------------------------------------------------
175 :
176 : #undef LG_FREE_WORK
177 : #define LG_FREE_WORK \
178 : { \
179 : GrB_UnaryOp_free (&LG_rand_next_op) ; \
180 : GrB_IndexUnaryOp_free (&LG_rand_init_op) ; \
181 : }
182 :
183 : #if !LAGRAPH_SUITESPARSE
184 : typedef void (*grb_unary_function) (void *, const void *) ;
185 : typedef void (*grb_index_unary_function)
186 : (
187 : void *z, // output value z, of type ztype
188 : const void *x, // input value x of type xtype; value of v(i) or A(i,j)
189 : GrB_Index i, // row index of A(i,j)
190 : GrB_Index j, // column index of A(i,j), or zero for v(i)
191 : const void *y // input scalar y
192 : ) ;
193 : #endif
194 :
195 226 : int LG_Random_Init (char *msg)
196 : {
197 226 : LG_CLEAR_MSG ;
198 226 : LG_FREE_WORK ; // free the two ops in case LG_Random_Init is called twice
199 226 : LG_rand_next_op = NULL ;
200 226 : LG_rand_init_op = NULL ;
201 :
202 : #if LAGRAPH_SUITESPARSE
203 : {
204 : // give SuiteSparse:GraphBLAS the strings that define the functions
205 226 : GRB_TRY (GxB_UnaryOp_new (&LG_rand_next_op,
206 : (GxB_unary_function) LG_rand_next_f2,
207 : GrB_UINT64, GrB_UINT64,
208 : "LG_rand_next_f2", LG_RAND_NEXT_F2_DEFN)) ;
209 224 : GRB_TRY (GxB_IndexUnaryOp_new (&LG_rand_init_op,
210 : (GxB_index_unary_function) LG_rand_init_func,
211 : GrB_UINT64, GrB_UINT64, GrB_UINT64,
212 : "LG_rand_init_func", LG_RAND_INIT_F_DEFN)) ;
213 : }
214 : #else
215 : {
216 : // vanilla GraphBLAS, no strings to define the new operators
217 : GRB_TRY (GrB_UnaryOp_new (&LG_rand_next_op,
218 : (grb_unary_function) LG_rand_next_f2,
219 : GrB_UINT64, GrB_UINT64)) ;
220 : GRB_TRY (GrB_IndexUnaryOp_new (&LG_rand_init_op,
221 : (grb_index_unary_function) LG_rand_init_func,
222 : GrB_UINT64, GrB_UINT64, GrB_UINT64)) ;
223 : }
224 : #endif
225 :
226 222 : return (GrB_SUCCESS) ;
227 : }
228 :
229 : //------------------------------------------------------------------------------
230 : // LG_Random_Finalize: free the random state operator
231 : //------------------------------------------------------------------------------
232 :
233 282 : int LG_Random_Finalize (char *msg)
234 : {
235 282 : LG_CLEAR_MSG ;
236 282 : LG_FREE_WORK ;
237 282 : return (GrB_SUCCESS) ;
238 : }
239 :
240 : //------------------------------------------------------------------------------
241 : // LAGraph_Random_Seed: create a vector of random states
242 : //------------------------------------------------------------------------------
243 :
244 : #undef LG_FREE_WORK
245 : #define LG_FREE_WORK ;
246 :
247 : // Initializes a vector with random state values. The State vector must be
248 : // allocated on input, and should be of type GrB_UINT64. Its sparsity
249 : // structure is unchanged.
250 :
251 : #if defined ( COVERAGE )
252 : // for testing only
253 : bool random_hack = false ;
254 : #endif
255 :
256 5930 : int LAGraph_Random_Seed // construct a random State vector
257 : (
258 : // input/output:
259 : GrB_Vector State, // vector of random number States, normally GrB_UINT64
260 : // input:
261 : uint64_t seed, // scalar input seed
262 : char *msg
263 : )
264 : {
265 : // check inputs
266 5930 : LG_CLEAR_MSG ;
267 5930 : LG_ASSERT (State != NULL, GrB_NULL_POINTER) ;
268 :
269 : // State (i) = splitmix64 (i + seed) for all prior entries in State
270 5930 : GRB_TRY (GrB_apply (State, NULL, NULL, LG_rand_init_op, State, seed,
271 : NULL)) ;
272 :
273 : #if defined ( COVERAGE )
274 5924 : if (random_hack)
275 : {
276 : // Set all State values to 1, to break the random seed vector.
277 : // This is just for testing, to test algorithms that need to handle
278 : // extreme cases when the random number generator is non-random.
279 433 : GRB_TRY (GrB_apply (State, NULL, NULL, GrB_ONEB_UINT64, State, 0,
280 : NULL)) ;
281 : }
282 : #endif
283 :
284 5924 : return (GrB_SUCCESS) ;
285 : }
286 :
287 : //------------------------------------------------------------------------------
288 : // LAGraph_Random_Next: return next vector of random seeds
289 : //------------------------------------------------------------------------------
290 :
291 20583 : int LAGraph_Random_Next // advance to next random vector
292 : (
293 : // input/output:
294 : GrB_Vector State, // vector of random number States, normally GrB_UINT64
295 : char *msg
296 : )
297 : {
298 : // check inputs
299 20583 : LG_CLEAR_MSG ;
300 20583 : LG_ASSERT (State != NULL, GrB_NULL_POINTER) ;
301 : // State = xorshift64 (State)
302 20583 : GRB_TRY (GrB_apply (State, NULL, NULL, LG_rand_next_op, State, NULL)) ;
303 20583 : return (GrB_SUCCESS) ;
304 : }
305 :
|
