Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_FastGraphletTransform: fast graphlet transform
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 Tanner Hoke, Texas A&M University
15 :
16 : //------------------------------------------------------------------------------
17 :
18 : // LAGraph_FastGraphletTransform: computes the Fast Graphlet Transform of
19 : // an undirected graph. No self edges are allowed on the input graph.
20 :
21 : // TODO: rename this, and add to src (but it uses GxB)
22 :
23 : // https://arxiv.org/pdf/2007.11111.pdf
24 :
25 : //------------------------------------------------------------------------------
26 : #define F_UNARY(f) ((void (*)(void *, const void *)) f)
27 :
28 : #define LG_FREE_WORK \
29 : { \
30 : GrB_free (&C_3) ; \
31 : GrB_free (&d_0) ; \
32 : GrB_free (&d_1) ; \
33 : GrB_free (&d_2) ; \
34 : GrB_free (&d_3) ; \
35 : GrB_free (&d_4) ; \
36 : GrB_free (&d_5) ; \
37 : GrB_free (&d_6) ; \
38 : GrB_free (&d_7) ; \
39 : GrB_free (&d_8) ; \
40 : GrB_free (&d_9) ; \
41 : GrB_free (&d_10) ; \
42 : GrB_free (&d_11) ; \
43 : GrB_free (&d_12) ; \
44 : GrB_free (&d_13) ; \
45 : GrB_free (&d_14) ; \
46 : GrB_free (&d_15) ; \
47 : GrB_free (&d_2) ; \
48 : GrB_free (&v) ; \
49 : GrB_free (&p_1_minus_one) ; \
50 : GrB_free (&p_1_minus_two) ; \
51 : GrB_free (&two_c_3) ; \
52 : GrB_free (&p_1_p_1_had) ; \
53 : GrB_free (&C_42) ; \
54 : GrB_free (&P_2) ; \
55 : GrB_free (&D_1) ; \
56 : GrB_free (&D_4c) ; \
57 : GrB_free (&D_43) ; \
58 : GrB_free (&U_inv) ; \
59 : GrB_free (&F_raw) ; \
60 : GrB_free (&C_4) ; \
61 : GrB_free (&Sub_one_mult) ; \
62 : GrB_free (&T) ; \
63 : if (A_Tiles != NULL) \
64 : { \
65 : for (int i = 0; i < tile_cnt; ++i) GrB_free (&A_Tiles [i]) ; \
66 : } \
67 : if (D_Tiles != NULL) \
68 : { \
69 : for (int i = 0; i < tile_cnt; ++i) GrB_free (&D_Tiles [i]) ; \
70 : } \
71 : if (C_Tiles != NULL) \
72 : { \
73 : for (int i = 0; i < tile_cnt; ++i) GrB_free (&C_Tiles [i]) ; \
74 : } \
75 : LAGraph_Free ((void **) &neighbors, msg) ; \
76 : LAGraph_Free ((void **) &k4cmn, msg) ; \
77 : LAGraph_Free ((void **) &f15, msg) ; \
78 : LAGraph_Free ((void **) &I, msg) ; \
79 : LAGraph_Free ((void **) &isNeighbor, msg) ; \
80 : LAGraph_Free ((void **) &J, msg) ; \
81 : LAGraph_Free ((void **) &vals, msg) ; \
82 : LAGraph_Free ((void **) &A_Tiles, msg) ; \
83 : LAGraph_Free ((void **) &D_Tiles, msg) ; \
84 : LAGraph_Free ((void **) &C_Tiles, msg) ; \
85 : LAGraph_Free ((void **) &Tile_nrows, msg) ; \
86 : }
87 :
88 : #define LG_FREE_ALL \
89 : { \
90 : LG_FREE_WORK ; \
91 : }
92 :
93 : #define F_UNARY(f) ((void (*)(void *, const void *)) f)
94 :
95 : #include "LG_internal.h"
96 : #include "LAGraphX.h"
97 :
98 993 : void sub_one_mult (int64_t *z, const int64_t *x) { (*z) = (*x) * ((*x)-1) ; }
99 :
100 2 : int LAGraph_FastGraphletTransform
101 : (
102 : // outputs:
103 : GrB_Matrix *F_net, // 16-by-n matrix of graphlet counts
104 : // inputs:
105 : LAGraph_Graph G,
106 : bool compute_d_15, // probably this makes most sense
107 : char *msg
108 : )
109 : {
110 : #if LAGRAPH_SUITESPARSE
111 :
112 2 : LG_CLEAR_MSG ;
113 2 : GrB_Index const U_inv_I[] = {0, 1, 2, 2, 3, 3, 4, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 12, 12, 12, 12, 13, 13, 14, 14, 15} ;
114 2 : GrB_Index const U_inv_J[] = {0, 1, 2, 4, 3, 4, 4, 5, 9, 10, 12, 13, 14, 15, 6, 10, 11, 12, 13, 14, 15, 7, 9, 10, 13, 14, 15, 8, 11, 14, 15, 9, 13, 15, 10, 13, 14, 15, 11, 14, 15, 12, 13, 14, 15, 13, 15, 14, 15, 15} ;
115 2 : int64_t const U_inv_X[] = {1, 1, 1, -2, 1, -1, 1, 1, -2, -1, -2, 4, 2, -6, 1, -1, -2, -2, 2, 4, -6, 1, -1, -1, 2, 1, -3, 1, -1, 1, -1, 1, -2, 3, 1, -2, -2, 6, 1, -2, 3, 1, -1, -1, 3, 1, -3, 1, -3, 1} ;
116 2 : GrB_Index const U_inv_nvals = 50;
117 2 : GrB_UnaryOp Sub_one_mult = NULL ;
118 2 : int tile_cnt = 0 ;
119 2 : GrB_Matrix *A_Tiles = NULL ;
120 2 : GrB_Matrix *D_Tiles = NULL ;
121 2 : GrB_Matrix *C_Tiles = NULL ;
122 2 : GrB_Index *Tile_nrows = NULL ;
123 2 : GrB_Matrix T = NULL ;
124 :
125 2 : GrB_Index *neighbors = NULL ;
126 2 : GrB_Index *k4cmn = NULL ;
127 2 : int64_t *f15 = NULL ;
128 2 : GrB_Index *I = NULL ;
129 2 : int *isNeighbor = NULL ;
130 2 : GrB_Index *J = NULL ;
131 2 : int64_t *vals = NULL ;
132 :
133 2 : GrB_Matrix C_3 = NULL,
134 2 : A = NULL,
135 2 : C_42 = NULL,
136 2 : P_2 = NULL,
137 2 : D_1 = NULL,
138 2 : D_4c = NULL,
139 2 : D_43 = NULL,
140 2 : U_inv = NULL,
141 2 : F_raw = NULL,
142 2 : C_4 = NULL ;
143 :
144 2 : GrB_Vector d_0 = NULL,
145 2 : d_1 = NULL,
146 2 : d_2 = NULL,
147 2 : d_3 = NULL,
148 2 : d_4 = NULL,
149 2 : d_5 = NULL,
150 2 : d_6 = NULL,
151 2 : d_7 = NULL,
152 2 : d_8 = NULL,
153 2 : d_9 = NULL,
154 2 : d_10 = NULL,
155 2 : d_11 = NULL,
156 2 : d_12 = NULL,
157 2 : d_13 = NULL,
158 2 : d_14 = NULL,
159 2 : d_15 = NULL;
160 :
161 2 : GrB_Vector v = NULL,
162 2 : two_c_3 = NULL,
163 2 : p_1_minus_one = NULL,
164 2 : p_1_minus_two = NULL,
165 2 : p_1_p_1_had = NULL ;
166 :
167 : GrB_Index nvals ;
168 : int64_t ntri ;
169 :
170 2 : A = G->A ;
171 :
172 : GrB_Index n ;
173 2 : GRB_TRY (GrB_Matrix_nrows (&n, A)) ;
174 :
175 : //--------------------------------------------------------------------------
176 : // compute d_0 = e
177 : //--------------------------------------------------------------------------
178 :
179 : // d_0 = e
180 2 : GRB_TRY (GrB_Vector_new (&d_0, GrB_INT64, n)) ;
181 2 : GRB_TRY (GrB_assign (d_0, NULL, NULL, 1, GrB_ALL, n, NULL)) ;
182 :
183 : //--------------------------------------------------------------------------
184 : // compute d_1 = Ae (in_degree)
185 : //--------------------------------------------------------------------------
186 :
187 : // GRB_TRY (GrB_Vector_new (&d_1, GrB_INT64, n)) ;
188 :
189 : // d_1 = Ae (in_degree)
190 2 : LG_TRY (LAGraph_Cached_OutDegree (G, msg)) ;
191 :
192 2 : GRB_TRY (GrB_Vector_dup (&d_1, G->out_degree)) ;
193 :
194 : //--------------------------------------------------------------------------
195 : // compute d_2 = p_2
196 : //--------------------------------------------------------------------------
197 :
198 2 : GRB_TRY (GrB_Vector_new (&d_2, GrB_INT64, n)) ;
199 :
200 : // TODO: use LAGraph_plus_second_int64
201 : // d_2 = p_2 = A*p_1 - c_2 = A*d_1 - d_1
202 2 : GRB_TRY (GrB_mxv (d_2, NULL, NULL, GxB_PLUS_SECOND_INT64, A, d_1, NULL)) ;
203 2 : GRB_TRY (GrB_eWiseMult (d_2, NULL, NULL, GrB_MINUS_INT64, d_2, d_1, NULL)) ;
204 :
205 : //--------------------------------------------------------------------------
206 : // compute d_3 = hadamard(p_1, p_1 - 1) / 2
207 : //--------------------------------------------------------------------------
208 :
209 2 : GRB_TRY (GrB_Vector_new (&d_3, GrB_INT64, n)) ;
210 :
211 2 : GRB_TRY (GrB_UnaryOp_new (&Sub_one_mult, F_UNARY (sub_one_mult), GrB_INT64, GrB_INT64)) ;
212 :
213 2 : GRB_TRY (GrB_apply (d_3, NULL, NULL, Sub_one_mult, d_1, NULL)) ;
214 2 : GRB_TRY (GrB_apply (d_3, NULL, NULL, GrB_DIV_INT64, d_3, (int64_t) 2, NULL)) ;
215 :
216 : //--------------------------------------------------------------------------
217 : // compute d_4 = C_3e/2
218 : //--------------------------------------------------------------------------
219 :
220 2 : GRB_TRY (GrB_Matrix_new (&C_3, GrB_INT64, n, n)) ;
221 2 : GRB_TRY (GrB_Vector_new (&d_4, GrB_INT64, n)) ;
222 :
223 : // TODO: use LAGraph_plus_first_int64
224 : // C_3 = hadamard(A, A^2)
225 2 : GRB_TRY (GrB_mxm (C_3, A, NULL, GxB_PLUS_FIRST_INT64, A, A, GrB_DESC_ST1)) ;
226 :
227 : // d_4 = c_3 = C_3e/2
228 2 : GRB_TRY (GrB_reduce (d_4, NULL, NULL, GrB_PLUS_MONOID_INT64, C_3, NULL)) ;
229 2 : GRB_TRY (GrB_apply (d_4, NULL, NULL, GrB_DIV_INT64, d_4, (int64_t) 2, NULL)) ;
230 :
231 : //--------------------------------------------------------------------------
232 : // compute d_5 = p_3 = A*d_2 - hadamard(p_1, p_1 - 1) - 2c_3
233 : //--------------------------------------------------------------------------
234 :
235 2 : GRB_TRY (GrB_Vector_new (&v, GrB_INT64, n)) ;
236 2 : GRB_TRY (GrB_Vector_new (&two_c_3, GrB_INT64, n)) ;
237 2 : GRB_TRY (GrB_Vector_new (&d_5, GrB_INT64, n)) ;
238 :
239 : // v = hadamard(p_1, p_1 - 1)
240 2 : GRB_TRY (GrB_apply (v, NULL, NULL, Sub_one_mult, d_1, NULL)) ;
241 :
242 : // two_c_3 = 2 * c_3 = 2 * d_4
243 2 : GRB_TRY (GrB_apply (two_c_3, NULL, NULL, GrB_TIMES_INT64, 2, d_4, NULL)) ;
244 :
245 : // d_5 = A * d_2
246 : // TODO: use LAGraph_plus_second_int64
247 2 : GRB_TRY (GrB_mxv (d_5, NULL, NULL, GxB_PLUS_SECOND_INT64, A, d_2, NULL)) ;
248 :
249 : // d_5 -= hadamard(p_1, p_1 - 1)
250 2 : GRB_TRY (GrB_eWiseAdd (d_5, NULL, NULL, GrB_MINUS_INT64, d_5, v, NULL)) ;
251 :
252 : // d_5 -= two_c_3
253 2 : GRB_TRY (GrB_eWiseAdd (d_5, NULL, NULL, GrB_MINUS_INT64, d_5, two_c_3, NULL)) ;
254 :
255 : //--------------------------------------------------------------------------
256 : // compute d_6 = hadamard(d_2, p_1-1) - 2c_3
257 : //--------------------------------------------------------------------------
258 :
259 2 : GRB_TRY (GrB_Vector_new (&p_1_minus_one, GrB_INT64, n)) ;
260 2 : GRB_TRY (GrB_Vector_new (&d_6, GrB_INT64, n)) ;
261 :
262 : // p_1_minus_one = p_1 - 1
263 2 : GRB_TRY (GrB_apply (p_1_minus_one, NULL, NULL, GrB_MINUS_INT64, d_1, (int64_t) 1, NULL)) ;
264 :
265 : // d_6 = hadamard(d_2, p_1-1)
266 2 : GRB_TRY (GrB_eWiseMult (d_6, NULL, NULL, GrB_TIMES_INT64, d_2, p_1_minus_one, NULL)) ;
267 :
268 : // d_6 -= 2c_3
269 2 : GRB_TRY (GrB_eWiseAdd (d_6, NULL, NULL, GrB_MINUS_INT64, d_6, two_c_3, NULL)) ;
270 :
271 : //--------------------------------------------------------------------------
272 : // compute d_7 = A*hadamard(p_1-1, p_1-2) / 2
273 : //--------------------------------------------------------------------------
274 :
275 2 : GRB_TRY (GrB_Vector_new (&p_1_minus_two, GrB_INT64, n)) ;
276 2 : GRB_TRY (GrB_Vector_new (&p_1_p_1_had, GrB_INT64, n)) ;
277 2 : GRB_TRY (GrB_Vector_new (&d_7, GrB_INT64, n)) ;
278 :
279 2 : GRB_TRY (GrB_apply (p_1_minus_two, NULL, NULL, GrB_MINUS_INT64, d_1, (int64_t) 2, NULL)) ;
280 2 : GRB_TRY (GrB_eWiseMult (p_1_p_1_had, NULL, NULL, GrB_TIMES_INT64, p_1_minus_one, p_1_minus_two, NULL)) ;
281 :
282 : // TODO: use LAGraph_plus_second_int64
283 2 : GRB_TRY (GrB_mxv (d_7, NULL, NULL, GxB_PLUS_SECOND_INT64, A, p_1_p_1_had, NULL)) ;
284 2 : GRB_TRY (GrB_apply (d_7, NULL, NULL, GrB_DIV_INT64, d_7, (int64_t) 2, NULL)) ;
285 :
286 : //--------------------------------------------------------------------------
287 : // compute d_8 = hadamard(p_1, p_1_p_1_had) / 6
288 : //--------------------------------------------------------------------------
289 :
290 2 : GRB_TRY (GrB_Vector_new (&d_8, GrB_INT64, n)) ;
291 :
292 2 : GRB_TRY (GrB_eWiseMult (d_8, NULL, NULL, GrB_TIMES_INT64, d_1, p_1_p_1_had, NULL)) ;
293 2 : GRB_TRY (GrB_apply (d_8, NULL, NULL, GrB_DIV_INT64, d_8, (int64_t) 6, NULL)) ;
294 :
295 : //--------------------------------------------------------------------------
296 : // compute d_9 = A*c_3 - 2*c_3
297 : //--------------------------------------------------------------------------
298 :
299 2 : GRB_TRY (GrB_Vector_new (&d_9, GrB_INT64, n)) ;
300 :
301 : // TODO: use LAGraph_plus_second_int64
302 2 : GRB_TRY (GrB_mxv (d_9, NULL, NULL, GxB_PLUS_SECOND_INT64, A, d_4, NULL)) ;
303 2 : GRB_TRY (GrB_eWiseAdd (d_9, NULL, NULL, GrB_MINUS_INT64, d_9, two_c_3, NULL)) ;
304 :
305 : //--------------------------------------------------------------------------
306 : // compute d_10 = C_3 * (p_1 - 2)
307 : //--------------------------------------------------------------------------
308 :
309 2 : GRB_TRY (GrB_Vector_new (&d_10, GrB_INT64, n)) ;
310 :
311 : // TODO: use GrB_PLUS_TIMES_INT64
312 2 : GRB_TRY (GrB_mxv (d_10, NULL, NULL, GxB_PLUS_TIMES_INT64, C_3, p_1_minus_two, NULL)) ;
313 :
314 : //--------------------------------------------------------------------------
315 : // compute d_11 = hadamard(p_1 - 2, c_3)
316 : //--------------------------------------------------------------------------
317 :
318 2 : GRB_TRY (GrB_Vector_new (&d_11, GrB_INT64, n)) ;
319 :
320 2 : GRB_TRY (GrB_eWiseMult (d_11, NULL, NULL, GrB_TIMES_INT64, p_1_minus_two, d_4, NULL)) ;
321 :
322 : //--------------------------------------------------------------------------
323 : // compute d_12 = c_4 = C_{4,2}e/2
324 : //--------------------------------------------------------------------------
325 :
326 2 : GRB_TRY (GrB_Matrix_new (&C_4, GrB_INT64, n, 1)) ;
327 2 : GRB_TRY (GrB_Matrix_new (&D_1, GrB_INT64, n, n)) ;
328 :
329 2 : GRB_TRY (GrB_Vector_new (&d_12, GrB_INT64, n)) ;
330 :
331 : // D_1 = diag(d_1)
332 : // TODO: use GrB_Matrix_diag
333 2 : GRB_TRY (GxB_Matrix_diag (D_1, d_1, (int64_t) 0, NULL)) ;
334 :
335 : // TODO: true GxB starts here. For vanilla method, put this
336 : // in a single call to GrB_mxm to compute C_4
337 :
338 2 : GRB_TRY (GrB_Matrix_nvals (&nvals, A));
339 :
340 2 : const GrB_Index entries_per_tile = 1000;
341 2 : GrB_Index ntiles = (nvals + entries_per_tile - 1) / entries_per_tile ;
342 :
343 2 : LG_TRY (LAGraph_Calloc ((void **) &A_Tiles, ntiles, sizeof (GrB_Matrix),
344 : msg)) ;
345 2 : LG_TRY (LAGraph_Calloc ((void **) &D_Tiles, ntiles, sizeof (GrB_Matrix),
346 : msg)) ;
347 2 : LG_TRY (LAGraph_Calloc ((void **) &C_Tiles, ntiles, sizeof (GrB_Matrix),
348 : msg)) ;
349 2 : LG_TRY (LAGraph_Calloc ((void **) &Tile_nrows, ntiles, sizeof (GrB_Index),
350 : msg)) ;
351 :
352 2 : GrB_Index Tile_ncols [1] = {n} ;
353 :
354 2 : int64_t tot_deg = 0 ;
355 2 : GrB_Index last_row = -1 ;
356 43 : for (GrB_Index i = 0; i < n; ++i)
357 : {
358 41 : int64_t deg = 0 ;
359 41 : GRB_TRY (GrB_Vector_extractElement (°, d_1, i)) ;
360 41 : if (i == n - 1 || (tot_deg / entries_per_tile != (tot_deg + deg) / entries_per_tile))
361 : {
362 2 : Tile_nrows [tile_cnt++] = i - last_row ;
363 2 : last_row = i ;
364 : }
365 41 : tot_deg += deg ;
366 : }
367 :
368 2 : GRB_TRY (GxB_Matrix_split (A_Tiles, tile_cnt, 1, Tile_nrows, Tile_ncols, A, NULL)) ;
369 2 : GRB_TRY (GxB_Matrix_split (D_Tiles, tile_cnt, 1, Tile_nrows, Tile_ncols, D_1, NULL)) ;
370 :
371 : #define TRY(method) \
372 : { \
373 : GrB_Info info2 = method ; \
374 : if (info2 != GrB_SUCCESS) \
375 : { \
376 : GrB_free (&A_i) ; \
377 : GrB_free (&C_Tiles [i_tile]) ; \
378 : GrB_free (&e) ; \
379 : info1 = info2 ; \
380 : continue ; \
381 : } \
382 : }
383 :
384 : int save_nthreads_outer, save_nthreads_inner ;
385 2 : LG_TRY (LAGraph_GetNumThreads (&save_nthreads_outer, &save_nthreads_inner, msg)) ;
386 2 : LG_TRY (LAGraph_SetNumThreads (1, 1, msg)) ;
387 :
388 : int i_tile;
389 2 : GrB_Info info1 = GrB_SUCCESS ;
390 : #pragma omp parallel for num_threads(omp_get_max_threads()) schedule(dynamic,1)
391 4 : for (i_tile = 0; i_tile < tile_cnt; ++i_tile)
392 : {
393 2 : GrB_Matrix A_i = NULL, e = NULL ;
394 :
395 2 : TRY (GrB_Matrix_new (&e, GrB_INT64, n, 1)) ;
396 2 : TRY (GrB_assign (e, NULL, NULL, (int64_t) 1, GrB_ALL, n, GrB_ALL, 1, NULL)) ;
397 :
398 2 : TRY (GrB_Matrix_new (&A_i, GrB_INT64, Tile_nrows [i_tile], n)) ;
399 2 : TRY (GrB_Matrix_new (&C_Tiles [i_tile], GrB_INT64, Tile_nrows [i_tile], 1)) ;
400 :
401 2 : TRY (GrB_mxm (A_i, NULL, NULL, GxB_PLUS_PAIR_INT64, A_Tiles [i_tile], A, NULL)) ;
402 2 : TRY (GrB_eWiseAdd (A_i, NULL, NULL, GrB_MINUS_INT64, A_i, D_Tiles [i_tile], NULL)) ;
403 2 : TRY (GrB_apply (A_i, NULL, NULL, Sub_one_mult, A_i, NULL)) ;
404 :
405 : // multiply A_i by it on the right
406 2 : TRY (GrB_mxm (C_Tiles [i_tile], NULL, NULL, GxB_PLUS_FIRST_INT64, A_i, e, NULL)) ;
407 :
408 2 : GrB_free (&A_i) ;
409 2 : GrB_free (&e) ;
410 : }
411 :
412 2 : GRB_TRY (info1) ;
413 :
414 2 : LG_TRY (LAGraph_SetNumThreads (save_nthreads_outer, save_nthreads_inner, msg)) ;
415 :
416 2 : GRB_TRY (GxB_Matrix_concat (C_4, C_Tiles, tile_cnt, 1, NULL)) ;
417 :
418 2 : if (A_Tiles != NULL)
419 : {
420 4 : for (int i = 0; i < tile_cnt; ++i) GrB_free (&A_Tiles [i]) ;
421 : }
422 2 : if (D_Tiles != NULL)
423 : {
424 4 : for (int i = 0; i < tile_cnt; ++i) GrB_free (&D_Tiles [i]) ;
425 : }
426 2 : if (C_Tiles != NULL)
427 : {
428 4 : for (int i = 0; i < tile_cnt; ++i) GrB_free (&C_Tiles [i]) ;
429 : }
430 2 : LAGraph_Free ((void **) &Tile_nrows, msg) ;
431 2 : LAGraph_Free ((void **) &A_Tiles, msg) ;
432 2 : LAGraph_Free ((void **) &D_Tiles, msg) ;
433 2 : LAGraph_Free ((void **) &C_Tiles, msg) ;
434 :
435 : // d_12 = sum (C_4)
436 2 : GRB_TRY (GrB_reduce (d_12, NULL, NULL, GrB_PLUS_MONOID_INT64, C_4, NULL)) ;
437 2 : GRB_TRY (GrB_apply (d_12, NULL, NULL, GrB_DIV_INT64, d_12, 2, NULL)) ;
438 :
439 : //--------------------------------------------------------------------------
440 : // compute d_13 = D_{4,c}e/2
441 : //--------------------------------------------------------------------------
442 :
443 2 : GRB_TRY (GrB_Matrix_new (&D_4c, GrB_INT64, n, n)) ;
444 2 : GRB_TRY (GrB_Vector_new (&d_13, GrB_INT64, n)) ;
445 :
446 2 : GRB_TRY (GrB_eWiseMult (D_4c, NULL, NULL, GrB_MINUS_INT64, C_3, A, NULL)) ; // can be mult because we mask with A next
447 2 : GRB_TRY (GrB_mxm (D_4c, A, NULL, GxB_PLUS_SECOND_INT64, A, D_4c, GrB_DESC_S)) ;
448 :
449 : // d_13 = D_{4,c}*e/2
450 2 : GRB_TRY (GrB_reduce (d_13, NULL, NULL, GrB_PLUS_INT64, D_4c, NULL)) ;
451 2 : GRB_TRY (GrB_apply (d_13, NULL, NULL, GrB_DIV_INT64, d_13, (int64_t) 2, NULL)) ;
452 :
453 : //--------------------------------------------------------------------------
454 : // compute d_14 = D_{4,3}e/2 = hadamard(A, C_42)e/2
455 : //--------------------------------------------------------------------------
456 :
457 2 : GRB_TRY (GrB_Matrix_new (&D_43, GrB_INT64, n, n)) ;
458 2 : GRB_TRY (GrB_Vector_new (&d_14, GrB_INT64, n)) ;
459 2 : GRB_TRY (GrB_Matrix_new (&C_42, GrB_INT64, n, n)) ;
460 2 : GRB_TRY (GrB_Matrix_new (&P_2, GrB_INT64, n, n)) ;
461 :
462 : // P_2 = A*A - diag(d_1)
463 2 : GRB_TRY (GrB_eWiseAdd (P_2, A, NULL, GrB_MINUS_INT64, C_3, D_1, NULL)) ;
464 :
465 : // C_42 = hadamard(P_2, P_2 - 1)
466 2 : GRB_TRY (GrB_apply (C_42, A, NULL, Sub_one_mult, P_2, NULL)) ;
467 :
468 2 : GRB_TRY (GrB_eWiseMult (D_43, NULL, NULL, GrB_TIMES_INT64, A, C_42, NULL)) ;
469 :
470 : // d_14 = D_{4,3}*e/2
471 2 : GRB_TRY (GrB_reduce (d_14, NULL, NULL, GrB_PLUS_INT64, D_43, NULL)) ;
472 2 : GRB_TRY (GrB_apply (d_14, NULL, NULL, GrB_DIV_INT64, d_14, (int64_t) 2, NULL)) ;
473 :
474 : //--------------------------------------------------------------------------
475 : // compute d_15 = Te/6
476 : //--------------------------------------------------------------------------
477 :
478 : // TODO: computing d_15 requires GxB iterators. For vanilla GrB,
479 : // export the CSR matrix and use it directly.
480 :
481 2 : if (compute_d_15)
482 : {
483 2 : LG_TRY (LAGraph_KTruss (&T, G, 4, msg)) ;
484 2 : GRB_TRY (GrB_Vector_new (&d_15, GrB_INT64, n)) ;
485 :
486 2 : int nthreads = 1 ;
487 : // todo: parallelize this...
488 : //#pragma omp parallel for num_threads(nthreads)
489 : //for (int tid = 0 ; tid < nthreads ; tid++)
490 : {
491 :
492 : // allocate workspcae
493 2 : LG_TRY (LAGraph_Malloc ((void **) &neighbors, n,
494 : sizeof (GrB_Index), msg)) ;
495 2 : LG_TRY (LAGraph_Malloc ((void **) &k4cmn, n,
496 : sizeof (GrB_Index), msg)) ;
497 2 : LG_TRY (LAGraph_Malloc ((void **) &f15, n,
498 : sizeof (int64_t), msg)) ;
499 2 : LG_TRY (LAGraph_Malloc ((void **) &I, n,
500 : sizeof (GrB_Index), msg)) ;
501 2 : LG_TRY (LAGraph_Malloc ((void **) &isNeighbor, n,
502 : sizeof (int), msg)) ;
503 :
504 43 : for (int i = 0; i < n; ++i) {
505 41 : neighbors [i] = k4cmn [i] = f15 [i] = isNeighbor [i] = 0 ;
506 41 : I [i] = i ;
507 : }
508 :
509 : // thread tid operates on T(row1:row2-1,:)
510 2 : GrB_Index row1 = 0;//tid * (n / nthreads) ;
511 2 : GrB_Index row2 = n;//(tid == nthreads - 1) ? n : ((tid+1) * (n / nthreads)) ;
512 :
513 : GxB_Iterator riterator ;
514 2 : GxB_Iterator_new (&riterator) ;
515 2 : GRB_TRY (GxB_rowIterator_attach (riterator, T, NULL)) ;
516 :
517 : GxB_Iterator iterator ;
518 2 : GxB_Iterator_new (&iterator) ;
519 2 : GRB_TRY (GxB_rowIterator_attach (iterator, T, NULL)) ;
520 :
521 : // seek to T(row1,:)
522 2 : GrB_Info info = GxB_rowIterator_seekRow (iterator, row1) ;
523 43 : while (info != GxB_EXHAUSTED)
524 : {
525 : // iterate over entries in T(i,:)
526 41 : GrB_Index idx2 = GxB_rowIterator_getRowIndex (iterator) ;
527 41 : if (idx2 >= row2) break ;
528 41 : int neighbor_cnt = 0 ;
529 109 : while (info == GrB_SUCCESS)
530 : {
531 : // working with edge (idx2, j)
532 68 : GrB_Index j = GxB_rowIterator_getColIndex (iterator) ;
533 :
534 68 : if (j > idx2) {
535 34 : neighbors [neighbor_cnt++] = j ;
536 34 : isNeighbor [j] = 1 ;
537 : }
538 :
539 68 : info = GxB_rowIterator_nextCol (iterator) ;
540 : }
541 :
542 75 : for (int neighbor_id = 0 ; neighbor_id < neighbor_cnt ; ++neighbor_id) {
543 34 : GrB_Index j = neighbors [neighbor_id] ;
544 34 : int cmn_cnt = 0 ;
545 34 : info = GxB_rowIterator_seekRow(riterator, j) ;
546 :
547 180 : while (info == GrB_SUCCESS) { // iterate over neighbors of j
548 146 : GrB_Index k = GxB_rowIterator_getColIndex (riterator) ;
549 146 : if (k > j && isNeighbor [k]) {
550 31 : k4cmn [cmn_cnt++] = k ;
551 31 : isNeighbor [k] = -1 ;
552 : }
553 146 : info = GxB_rowIterator_nextCol (riterator) ;
554 : }
555 : // check every combination
556 65 : for (int k_1 = 0 ; k_1 < cmn_cnt ; k_1++) {
557 31 : GrB_Index k = k4cmn [k_1] ;
558 31 : info = GxB_rowIterator_seekRow(riterator, k) ;
559 :
560 164 : while (info == GrB_SUCCESS) { // iterate over neighbors of k
561 133 : GrB_Index l = GxB_rowIterator_getColIndex (riterator) ;
562 133 : if (l > k && isNeighbor [l] == -1) {
563 13 : f15[idx2]++ ;
564 13 : f15[j]++ ;
565 13 : f15[k]++ ;
566 13 : f15[l]++ ;
567 : }
568 133 : info = GxB_rowIterator_nextCol (riterator) ;
569 : }
570 : }
571 65 : for (int k_1 = 0 ; k_1 < cmn_cnt ; k_1++) {
572 31 : isNeighbor[k4cmn[k_1]] = 1 ;
573 : }
574 : }
575 :
576 75 : for (int neighbor_id = 0 ; neighbor_id < neighbor_cnt ; ++neighbor_id) {
577 34 : GrB_Index j = neighbors [neighbor_id] ;
578 34 : isNeighbor [j] = 0 ;
579 : }
580 :
581 : // move to the next row, T(i+1,:)
582 41 : info = GxB_rowIterator_nextRow (iterator) ;
583 : }
584 2 : GrB_free (&iterator) ;
585 2 : GrB_free (&riterator) ;
586 2 : GrB_free (&T) ;
587 2 : GRB_TRY (GrB_Vector_build (d_15, I, f15, n, NULL)) ;
588 :
589 : // free workspace
590 2 : LAGraph_Free ((void **) &neighbors, msg) ;
591 2 : LAGraph_Free ((void **) &k4cmn, msg) ;
592 2 : LAGraph_Free ((void **) &f15, msg) ;
593 2 : LAGraph_Free ((void **) &I, msg) ;
594 2 : LAGraph_Free ((void **) &isNeighbor, msg) ;
595 : }
596 : }
597 :
598 : //--------------------------------------------------------------------------
599 : // construct raw frequencies matrix F_raw
600 : //--------------------------------------------------------------------------
601 :
602 2 : GRB_TRY (GrB_Matrix_new (&F_raw, GrB_INT64, 16, n)) ;
603 :
604 2 : GrB_Vector d[16] = {d_0, d_1, d_2, d_3, d_4, d_5, d_6, d_7, d_8, d_9, d_10, d_11, d_12, d_13, d_14, d_15} ;
605 :
606 34 : for (int i = 0; i < 15 + (compute_d_15 ? 1 : 0); ++i)
607 : {
608 32 : GRB_TRY (GrB_Vector_nvals (&nvals, d[i]));
609 :
610 : // allocate workspace
611 32 : LG_TRY (LAGraph_Malloc ((void **) &J, nvals, sizeof (GrB_Index), msg)) ;
612 32 : LG_TRY (LAGraph_Malloc ((void **) &vals, nvals, sizeof (int64_t),
613 : msg)) ;
614 :
615 32 : GRB_TRY (GrB_Vector_extractTuples (J, vals, &nvals, d[i])) ;
616 678 : for (int j = 0; j < nvals; ++j) {
617 646 : GRB_TRY (GrB_Matrix_setElement (F_raw, vals[j], i, J[j])) ;
618 : }
619 :
620 : // free workspace
621 32 : LAGraph_Free ((void **) &J, msg) ;
622 32 : LAGraph_Free ((void **) &vals, msg) ;
623 : }
624 :
625 : //--------------------------------------------------------------------------
626 : // construct U_inv
627 : //--------------------------------------------------------------------------
628 :
629 2 : GRB_TRY (GrB_Matrix_new (&U_inv, GrB_INT64, 16, 16)) ;
630 2 : GRB_TRY (GrB_Matrix_build (U_inv, U_inv_I, U_inv_J, U_inv_X, U_inv_nvals, GrB_PLUS_INT64)) ;
631 :
632 : //--------------------------------------------------------------------------
633 : // construct net frequencies matrix F_net
634 : //--------------------------------------------------------------------------
635 :
636 2 : GRB_TRY (GrB_Matrix_new (F_net, GrB_INT64, 16, n)) ;
637 2 : GRB_TRY (GrB_mxm (*F_net, NULL, NULL, GxB_PLUS_TIMES_INT64, U_inv, F_raw, NULL)) ;
638 2 : GrB_Vector f_net = NULL ;
639 2 : GRB_TRY (GrB_Vector_new (&f_net, GrB_INT64, 16)) ;
640 2 : GRB_TRY (GrB_reduce (f_net, NULL, NULL, GrB_PLUS_INT64, *F_net, NULL)) ;
641 2 : GRB_TRY (GrB_free (&f_net)) ;
642 :
643 : //--------------------------------------------------------------------------
644 : // free work
645 : //--------------------------------------------------------------------------
646 :
647 2 : LG_FREE_WORK ;
648 2 : return (GrB_SUCCESS) ;
649 : #else
650 : return (GrB_NOT_IMPLEMENTED) ;
651 : #endif
652 : }
|
