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