Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_Random_Matrix: generate a random matrix
3 : //------------------------------------------------------------------------------
4 :
5 : // LAGraph, (c) 2019-2022 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 : // Constructs a sparse roughly uniformly distributed random matrix with roughly
19 : // density*nrows*ncols entries. If density == INFINITY then the matrix is
20 : // generated with all entries present.
21 :
22 : // If the type is GrB_FP32 or GrB_FP64, the values of A are returned in the
23 : // range [0,1]. If any duplicate entries are generated, the largest one is
24 : // take, so the distribution can be skewed towards 1 if the density is large.
25 : // This could be fixed by using the GxB_IGNORE_DUP operator, but this would
26 : // require SuiteSparse:GraphBLAS.
27 :
28 : #define LG_FREE_WORK \
29 : { \
30 : LAGraph_Free ((void **) &I, NULL) ; \
31 : LAGraph_Free ((void **) &J, NULL) ; \
32 : LAGraph_Free ((void **) &ignore, NULL) ; \
33 : LAGraph_Free (&X, NULL) ; \
34 : GrB_free (&Mod) ; \
35 : GrB_free (&Rows) ; \
36 : GrB_free (&Cols) ; \
37 : GrB_free (&Values) ; \
38 : GrB_free (&State) ; \
39 : GrB_free (&T) ; \
40 : }
41 :
42 : #define LG_FREE_ALL \
43 : { \
44 : LG_FREE_WORK ; \
45 : GrB_free (A) ; \
46 : }
47 :
48 : #include "LG_internal.h"
49 : #include "LAGraphX.h"
50 :
51 : // uncomment these to test vanilla case for just this file:
52 : // #undef LAGRAPH_SUITESPARSE
53 : // #define LAGRAPH_SUITESPARSE 0
54 :
55 : //------------------------------------------------------------------------------
56 : // mod function for uint64: z = x % y
57 : //------------------------------------------------------------------------------
58 :
59 220 : void mod_function (void *z, const void *x, const void *y)
60 : {
61 220 : uint64_t a = (*((uint64_t *) x)) ;
62 220 : uint64_t b = (*((uint64_t *) y)) ;
63 220 : (*((uint64_t *) z)) = a % b ;
64 220 : }
65 :
66 : #define MOD_FUNCTION_DEFN \
67 : "void mod_function (void *z, const void *x, const void *y) \n" \
68 : "{ \n" \
69 : " uint64_t a = (*((uint64_t *) x)) ; \n" \
70 : " uint64_t b = (*((uint64_t *) y)) ; \n" \
71 : " (*((uint64_t *) z)) = a % b ; \n" \
72 : "}"
73 :
74 : //------------------------------------------------------------------------------
75 : // LAGraph_Random_Matrix
76 : //------------------------------------------------------------------------------
77 :
78 80 : GrB_Info LAGraph_Random_Matrix // random matrix of any built-in type
79 : (
80 : // output
81 : GrB_Matrix *A, // A is constructed on output
82 : // input
83 : GrB_Type type, // type of matrix to construct
84 : GrB_Index nrows, // # of rows of A
85 : GrB_Index ncols, // # of columns of A
86 : double density, // density: build a sparse matrix with
87 : // density*nrows*cols values if not INFINITY;
88 : // build a dense matrix if INFINITY.
89 : uint64_t seed, // random number seed
90 : char *msg
91 : )
92 : {
93 :
94 : //--------------------------------------------------------------------------
95 : // check inputs
96 : //--------------------------------------------------------------------------
97 :
98 80 : LG_CLEAR_MSG ;
99 80 : GrB_BinaryOp Mod = NULL ;
100 80 : GrB_Vector Rows = NULL, Cols = NULL, Values = NULL, State = NULL ;
101 80 : GrB_Matrix T = NULL ;
102 80 : GrB_Index *I = NULL, *J = NULL, *ignore = NULL ;
103 80 : GrB_Index I_size = 0, J_size = 0, X_size = 0 ;
104 80 : void *X = NULL ;
105 80 : LG_ASSERT (A != NULL && type != NULL, GrB_NULL_POINTER) ;
106 80 : LG_ASSERT_MSG (density >= 0, GrB_INVALID_VALUE, "invalid density") ;
107 :
108 80 : LG_ASSERT_MSG (type == GrB_BOOL
109 : || type == GrB_INT8 || type == GrB_INT16 || type == GrB_INT32
110 : || type == GrB_INT64 || type == GrB_UINT8 || type == GrB_UINT16
111 : || type == GrB_UINT32 || type == GrB_UINT64
112 : || type == GrB_FP32 || type == GrB_FP64,
113 : GrB_NOT_IMPLEMENTED, "unsupported type") ;
114 :
115 79 : GRB_TRY (GrB_Matrix_new (A, type, nrows, ncols)) ;
116 79 : if (nrows == 0 || ncols == 0)
117 : {
118 : // nothing to do: return A as the requested empty matrix
119 1 : return (GrB_SUCCESS) ;
120 : }
121 :
122 : //--------------------------------------------------------------------------
123 : // create the Mod operator
124 : //--------------------------------------------------------------------------
125 :
126 : #if LAGRAPH_SUITESPARSE
127 78 : GRB_TRY (GxB_BinaryOp_new (&Mod, mod_function,
128 : GrB_UINT64, GrB_UINT64, GrB_UINT64,
129 : "mod_function", MOD_FUNCTION_DEFN)) ;
130 : #else
131 : GRB_TRY (GrB_BinaryOp_new (&Mod, mod_function,
132 : GrB_UINT64, GrB_UINT64, GrB_UINT64)) ;
133 : #endif
134 :
135 : //--------------------------------------------------------------------------
136 : // determine the number of entries to generate
137 : //--------------------------------------------------------------------------
138 :
139 78 : bool A_is_full = isinf (density) ;
140 : GrB_Index nvals ;
141 78 : if (A_is_full)
142 : {
143 : // determine number of tuples for building a random dense matrix
144 17 : double nx = (double) nrows * (double) ncols ;
145 17 : LG_ASSERT_MSG (nx < (double) GrB_INDEX_MAX, GrB_OUT_OF_MEMORY,
146 : "Problem too large") ;
147 17 : nvals = nrows * ncols ;
148 : }
149 : else
150 : {
151 : // determine number of tuples for building a random sparse matrix
152 61 : double nx = density * (double) nrows * (double) ncols ;
153 61 : nx = round (nx) ;
154 61 : nx = fmax (nx, (double) 0) ;
155 61 : nx = fmin (nx, (double) GrB_INDEX_MAX) ;
156 61 : nvals = (GrB_Index) nx ;
157 : }
158 :
159 : //--------------------------------------------------------------------------
160 : // allocate workspace
161 : //--------------------------------------------------------------------------
162 :
163 : #if !LAGRAPH_SUITESPARSE
164 : {
165 : LG_TRY (LAGraph_Malloc ((void **) &ignore, nvals, sizeof (GrB_Index),
166 : msg)) ;
167 : LG_TRY (LAGraph_Malloc ((void **) &I, nvals, sizeof (GrB_Index), msg)) ;
168 : LG_TRY (LAGraph_Malloc ((void **) &J, nvals, sizeof (GrB_Index), msg)) ;
169 : }
170 : #endif
171 :
172 : //--------------------------------------------------------------------------
173 : // construct the random dense State vector of size nvals
174 : //--------------------------------------------------------------------------
175 :
176 78 : GRB_TRY (GrB_Vector_new (&State, GrB_UINT64, nvals)) ;
177 78 : GRB_TRY (GrB_assign (State, NULL, NULL, 0, GrB_ALL, nvals, NULL)) ;
178 78 : LG_TRY (LAGraph_Random_Seed (State, seed, msg)) ;
179 :
180 : //--------------------------------------------------------------------------
181 : // construct the random indices if A is sparse, or all indices if full
182 : //--------------------------------------------------------------------------
183 :
184 78 : if (!A_is_full)
185 : {
186 :
187 : //----------------------------------------------------------------------
188 : // construct random indices for a sparse matrix
189 : //----------------------------------------------------------------------
190 :
191 : // Rows = mod (State, nrows) ;
192 61 : GRB_TRY (GrB_Vector_new (&Rows, GrB_UINT64, nvals)) ;
193 61 : GRB_TRY (GrB_apply (Rows, NULL, NULL, Mod, State, nrows, NULL)) ;
194 :
195 : // State = next (State)
196 61 : LG_TRY (LAGraph_Random_Next (State, msg)) ;
197 :
198 : // Cols = mod (State, ncols) ;
199 61 : GRB_TRY (GrB_Vector_new (&Cols, GrB_UINT64, nvals)) ;
200 61 : GRB_TRY (GrB_apply (Cols, NULL, NULL, Mod, State, ncols, NULL)) ;
201 :
202 : // State = next (State)
203 61 : LG_TRY (LAGraph_Random_Next (State, msg)) ;
204 :
205 : //----------------------------------------------------------------------
206 : // extract the indices
207 : //----------------------------------------------------------------------
208 :
209 : #if LAGRAPH_SUITESPARSE
210 : {
211 : // this takes O(1) time and space
212 61 : GRB_TRY (GxB_Vector_unpack_Full (Rows, (void **) &I, &I_size,
213 : NULL, NULL)) ;
214 61 : GRB_TRY (GxB_Vector_unpack_Full (Cols, (void **) &J, &J_size,
215 : NULL, NULL)) ;
216 : }
217 : #else
218 : {
219 : // this takes O(nvals) time and space
220 : GRB_TRY (GrB_Vector_extractTuples_UINT64 (ignore, I, &nvals,
221 : Rows)) ;
222 : GRB_TRY (GrB_Vector_extractTuples_UINT64 (ignore, J, &nvals,
223 : Cols)) ;
224 : }
225 : #endif
226 :
227 61 : GrB_free (&Rows) ;
228 61 : GrB_free (&Cols) ;
229 :
230 : }
231 : else
232 : {
233 :
234 : //----------------------------------------------------------------------
235 : // construct indices for a full matrix
236 : //----------------------------------------------------------------------
237 :
238 : #if !LAGRAPH_SUITESPARSE
239 : {
240 : // T = true (nrows, ncols) ;
241 : GRB_TRY (GrB_Matrix_new (&T, GrB_BOOL, nrows, ncols)) ;
242 : GRB_TRY (GrB_assign (T, NULL, NULL, true,
243 : GrB_ALL, nrows, GrB_ALL, ncols, NULL)) ;
244 : // extract the row and column indices from T
245 : GRB_TRY (GrB_Matrix_extractTuples (I, J, (bool *) ignore, &nvals,
246 : T)) ;
247 : GrB_free (&T) ;
248 : }
249 : #endif
250 : }
251 :
252 : //-------------------------------------------------------------------------
253 : // construct the random values
254 : //-------------------------------------------------------------------------
255 :
256 78 : if (type == GrB_BOOL)
257 : {
258 : // Values = (State < UINT64_MAX / 2)
259 5 : GRB_TRY (GrB_Vector_new (&Values, GrB_BOOL, nvals)) ;
260 5 : GRB_TRY (GrB_apply (Values, NULL, NULL,
261 : GrB_LT_UINT64, State, UINT64_MAX / 2, NULL)) ;
262 : }
263 73 : else if (type == GrB_UINT64)
264 : {
265 : // no need to allocate the Values vector; just use the State itself
266 5 : Values = State ;
267 5 : State = NULL ;
268 : }
269 : else
270 : {
271 : // Values = (type) State
272 68 : GRB_TRY (GrB_Vector_new (&Values, type, nvals)) ;
273 68 : GRB_TRY (GrB_assign (Values, NULL, NULL, State, GrB_ALL, nvals,
274 : NULL)) ;
275 : }
276 78 : GrB_free (&State) ;
277 :
278 : // scale the values to the range [0,1], if floating-point
279 78 : if (type == GrB_FP32)
280 : {
281 : // Values = Values / (float) UINT64_MAX
282 8 : GRB_TRY (GrB_apply (Values, NULL, NULL, GrB_DIV_FP32,
283 : Values, (float) UINT64_MAX, NULL)) ;
284 : }
285 70 : else if (type == GrB_FP64)
286 : {
287 : // Values = Values / (double) UINT64_MAX
288 25 : GRB_TRY (GrB_apply (Values, NULL, NULL, GrB_DIV_FP64,
289 : Values, (double) UINT64_MAX, NULL)) ;
290 : }
291 :
292 : //--------------------------------------------------------------------------
293 : // extract the values
294 : //--------------------------------------------------------------------------
295 :
296 : #if LAGRAPH_SUITESPARSE
297 : {
298 : // this takes O(1) time and space and works for any data type
299 78 : GRB_TRY (GxB_Vector_unpack_Full (Values, &X, &X_size, NULL, NULL)) ;
300 : }
301 : #else
302 : {
303 : // this takes O(nvals) time and space
304 : if (type == GrB_BOOL)
305 : {
306 : LG_TRY (LAGraph_Malloc ((void **) &X, nvals, sizeof (bool), msg)) ;
307 : GRB_TRY (GrB_Vector_extractTuples_BOOL (ignore, X, &nvals,
308 : Values)) ;
309 : }
310 : else if (type == GrB_INT8)
311 : {
312 : LG_TRY (LAGraph_Malloc ((void **) &X, nvals, sizeof (int8_t),
313 : msg)) ;
314 : GRB_TRY (GrB_Vector_extractTuples_INT8 (ignore, X, &nvals,
315 : Values)) ;
316 : }
317 : else if (type == GrB_INT16)
318 : {
319 : LG_TRY (LAGraph_Malloc ((void **) &X, nvals, sizeof (int16_t),
320 : msg)) ;
321 : GRB_TRY (GrB_Vector_extractTuples_INT16 (ignore, X, &nvals,
322 : Values)) ;
323 : }
324 : else if (type == GrB_INT32)
325 : {
326 : LG_TRY (LAGraph_Malloc ((void **) &X, nvals, sizeof (int32_t),
327 : msg)) ;
328 : GRB_TRY (GrB_Vector_extractTuples_INT32 (ignore, X, &nvals,
329 : Values)) ;
330 : }
331 : else if (type == GrB_INT64)
332 : {
333 : LG_TRY (LAGraph_Malloc ((void **) &X, nvals, sizeof (int64_t),
334 : msg)) ;
335 : GRB_TRY (GrB_Vector_extractTuples_INT64 (ignore, X, &nvals,
336 : Values)) ;
337 : }
338 : else if (type == GrB_UINT8)
339 : {
340 : LG_TRY (LAGraph_Malloc ((void **) &X, nvals, sizeof (uint8_t),
341 : msg)) ;
342 : GRB_TRY (GrB_Vector_extractTuples_UINT8 (ignore, X, &nvals,
343 : Values)) ;
344 : }
345 : else if (type == GrB_UINT16)
346 : {
347 : LG_TRY (LAGraph_Malloc ((void **) &X, nvals, sizeof (uint16_t),
348 : msg)) ;
349 : GRB_TRY (GrB_Vector_extractTuples_UINT16 (ignore, X, &nvals,
350 : Values)) ;
351 : }
352 : else if (type == GrB_UINT32)
353 : {
354 : LG_TRY (LAGraph_Malloc ((void **) &X, nvals, sizeof (uint32_t),
355 : msg)) ;
356 : GRB_TRY (GrB_Vector_extractTuples_UINT32 (ignore, X, &nvals,
357 : Values)) ;
358 : }
359 : else if (type == GrB_UINT64)
360 : {
361 : LG_TRY (LAGraph_Malloc ((void **) &X, nvals, sizeof (uint64_t),
362 : msg)) ;
363 : GRB_TRY (GrB_Vector_extractTuples_UINT64 (ignore, X, &nvals,
364 : Values)) ;
365 : }
366 : else if (type == GrB_FP32)
367 : {
368 : LG_TRY (LAGraph_Malloc ((void **) &X, nvals, sizeof (float),
369 : msg)) ;
370 : GRB_TRY (GrB_Vector_extractTuples_FP32 (ignore, X, &nvals,
371 : Values)) ;
372 : }
373 : else // if (type == GrB_FP64)
374 : {
375 : LG_TRY (LAGraph_Malloc ((void **) &X, nvals, sizeof (double),
376 : msg)) ;
377 : GRB_TRY (GrB_Vector_extractTuples_FP64 (ignore, X, &nvals,
378 : Values)) ;
379 : }
380 : LAGraph_Free ((void **) &ignore, NULL) ;
381 : }
382 : #endif
383 :
384 78 : GrB_free (&Values) ;
385 :
386 : //--------------------------------------------------------------------------
387 : // build the matrix
388 : //--------------------------------------------------------------------------
389 :
390 : // Using GxB_IGNORE_DUP for the dup operator would be faster for
391 : // SuiteSparse, but it would result in a different matrix as compared to
392 : // the pure GrB case.
393 :
394 : #if LAGRAPH_SUITESPARSE
395 78 : if (A_is_full)
396 : {
397 : // this takes O(1) time and space
398 17 : GRB_TRY (GxB_Matrix_pack_FullR (*A, &X, X_size, false, NULL)) ;
399 : }
400 : else
401 : #endif
402 61 : if (type == GrB_BOOL)
403 : {
404 4 : GRB_TRY (GrB_Matrix_build_BOOL (*A, I, J, X, nvals, GrB_LXOR)) ;
405 : }
406 57 : else if (type == GrB_INT8)
407 : {
408 4 : GRB_TRY (GrB_Matrix_build_INT8 (*A, I, J, X, nvals, GrB_PLUS_INT8)) ;
409 : }
410 53 : else if (type == GrB_INT16)
411 : {
412 4 : GRB_TRY (GrB_Matrix_build_INT16 (*A, I, J, X, nvals, GrB_PLUS_INT16)) ;
413 : }
414 49 : else if (type == GrB_INT32)
415 : {
416 4 : GRB_TRY (GrB_Matrix_build_INT32 (*A, I, J, X, nvals, GrB_PLUS_INT32)) ;
417 : }
418 45 : else if (type == GrB_INT64)
419 : {
420 4 : GRB_TRY (GrB_Matrix_build_INT64 (*A, I, J, X, nvals, GrB_PLUS_INT64)) ;
421 : }
422 41 : else if (type == GrB_UINT8)
423 : {
424 4 : GRB_TRY (GrB_Matrix_build_UINT8 (*A, I, J, X, nvals, GrB_PLUS_UINT8)) ;
425 : }
426 37 : else if (type == GrB_UINT16)
427 : {
428 4 : GRB_TRY (GrB_Matrix_build_UINT16 (*A, I, J, X, nvals, GrB_PLUS_UINT16));
429 : }
430 33 : else if (type == GrB_UINT32)
431 : {
432 4 : GRB_TRY (GrB_Matrix_build_UINT32 (*A, I, J, X, nvals, GrB_PLUS_UINT32));
433 : }
434 29 : else if (type == GrB_UINT64)
435 : {
436 4 : GRB_TRY (GrB_Matrix_build_UINT64 (*A, I, J, X, nvals, GrB_PLUS_UINT64));
437 : }
438 25 : else if (type == GrB_FP32)
439 : {
440 4 : GRB_TRY (GrB_Matrix_build_FP32 (*A, I, J, X, nvals, GrB_MAX_FP32)) ;
441 : }
442 : else // if (type == GrB_FP64)
443 : {
444 21 : GRB_TRY (GrB_Matrix_build_FP64 (*A, I, J, X, nvals, GrB_MAX_FP64)) ;
445 : }
446 :
447 : //--------------------------------------------------------------------------
448 : // free workspace and return result
449 : //--------------------------------------------------------------------------
450 :
451 78 : LG_FREE_WORK ;
452 78 : return (GrB_SUCCESS) ;
453 : }
|
