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