Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGr_TriangleCount: Triangle counting using various methods
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 : // Count the number of triangles in a graph,
19 :
20 : // This is an Advanced algorithm (G->nself_edges, G->out_degree,
21 : // G->is_symmetric_structure are required).
22 :
23 : // Given a symmetric graph A with no-self edges, LAGr_TriangleCount counts the
24 : // number of triangles in the graph. A triangle is a clique of size three,
25 : // that is, 3 nodes that are all pairwise connected.
26 :
27 : // One of 6 methods are used, defined below where L and U are the strictly
28 : // lower and strictly upper triangular parts of the symmetrix matrix A,
29 : // respectively. Each method computes the same result, ntri:
30 :
31 : // 0: default: use the default method (currently method Sandia_LUT)
32 : // 1: Burkhardt: ntri = sum (sum ((A^2) .* A)) / 6
33 : // 2: Cohen: ntri = sum (sum ((L * U) .* A)) / 2
34 : // 3: Sandia_LL: ntri = sum (sum ((L * L) .* L))
35 : // 4: Sandia_UU: ntri = sum (sum ((U * U) .* U))
36 : // 5: Sandia_LUT: ntri = sum (sum ((L * U') .* L)). Note that L=U'.
37 : // 6: Sandia_ULT: ntri = sum (sum ((U * L') .* U)). Note that U=L'.
38 :
39 : // A is a square symmetric matrix, of any type. Its values are ignored.
40 : // Results are undefined for methods 1 and 2 if self-edges exist in A. Results
41 : // are undefined for all methods if A is unsymmetric.
42 :
43 : // The Sandia_* methods all tend to be faster than the Burkhardt or Cohen
44 : // methods. For the largest graphs, Sandia_LUT tends to be fastest, except for
45 : // the GAP-urand matrix, where the saxpy-based Sandia_LL method (L*L.*L) is
46 : // fastest. For many small graphs, the saxpy-based Sandia_LL and Sandia_UU
47 : // methods are often faster that the dot-product-based methods.
48 :
49 : // Reference for the Burkhardt method: Burkhardt, Paul. "Graphing Trillions of
50 : // Triangles." Information Visualization 16, no. 3 (July 2017): 157–66.
51 : // https://doi.org/10.1177/1473871616666393.
52 :
53 : // Reference for the Cohen method: J. Cohen, "Graph twiddling in a mapreduce
54 : // world," Computing in Science & Engineering, vol. 11, no. 4, pp. 29–41, 2009.
55 : // https://doi.org/10.1109/MCSE.2009.120
56 :
57 : // Reference for the "Sandia_*" methods: Wolf, Deveci, Berry, Hammond,
58 : // Rajamanickam, "Fast linear algebra- based triangle counting with
59 : // KokkosKernels", IEEE HPEC'17, https://dx.doi.org/10.1109/HPEC.2017.8091043
60 :
61 : #define LG_FREE_ALL \
62 : { \
63 : GrB_free (L) ; \
64 : GrB_free (U) ; \
65 : }
66 :
67 : #include "LG_internal.h"
68 :
69 : //------------------------------------------------------------------------------
70 : // tricount_prep: construct L and U for LAGr_TriangleCount
71 : //------------------------------------------------------------------------------
72 :
73 4235 : static int tricount_prep
74 : (
75 : GrB_Matrix *L, // if present, compute L = tril (A,-1)
76 : GrB_Matrix *U, // if present, compute U = triu (A, 1)
77 : GrB_Matrix A, // input matrix
78 : char *msg
79 : )
80 : {
81 : GrB_Index n ;
82 4235 : GRB_TRY (GrB_Matrix_nrows (&n, A)) ;
83 :
84 4235 : if (L != NULL)
85 : {
86 : // L = tril (A,-1)
87 3746 : GRB_TRY (GrB_Matrix_new (L, GrB_BOOL, n, n)) ;
88 3263 : GRB_TRY (GrB_select (*L, NULL, NULL, GrB_TRIL, A, (int64_t) (-1),
89 : NULL)) ;
90 2378 : GRB_TRY (GrB_Matrix_wait (*L, GrB_MATERIALIZE)) ;
91 : }
92 :
93 2867 : if (U != NULL)
94 : {
95 : // U = triu (A,1)
96 2626 : GRB_TRY (GrB_Matrix_new (U, GrB_BOOL, n, n)) ;
97 2116 : GRB_TRY (GrB_select (*U, NULL, NULL, GrB_TRIU, A, (int64_t) 1, NULL)) ;
98 1217 : GRB_TRY (GrB_Matrix_wait (*U, GrB_MATERIALIZE)) ;
99 : }
100 1458 : return (GrB_SUCCESS) ;
101 : }
102 :
103 : //------------------------------------------------------------------------------
104 : // LAGraph_tricount: count the number of triangles in a graph
105 : //------------------------------------------------------------------------------
106 :
107 : #undef LG_FREE_ALL
108 : #define LG_FREE_ALL \
109 : { \
110 : GrB_free (&C) ; \
111 : GrB_free (&L) ; \
112 : GrB_free (&T) ; \
113 : GrB_free (&U) ; \
114 : LAGraph_Free ((void **) &P, NULL) ; \
115 : }
116 :
117 5856 : int LAGr_TriangleCount
118 : (
119 : // output:
120 : uint64_t *ntriangles,
121 : // input:
122 : const LAGraph_Graph G,
123 : LAGr_TriangleCount_Method *p_method,
124 : LAGr_TriangleCount_Presort *p_presort,
125 : char *msg
126 : )
127 : {
128 :
129 : //--------------------------------------------------------------------------
130 : // check inputs
131 : //--------------------------------------------------------------------------
132 :
133 5856 : LG_CLEAR_MSG ;
134 5856 : GrB_Matrix C = NULL, L = NULL, U = NULL, T = NULL ;
135 5856 : int64_t *P = NULL ;
136 :
137 : // get the method
138 : LAGr_TriangleCount_Method method ;
139 5856 : method = (p_method == NULL) ? LAGr_TriangleCount_AutoMethod : (*p_method) ;
140 5856 : LG_ASSERT_MSG (
141 : method == LAGr_TriangleCount_AutoMethod || // 0: use auto method
142 : method == LAGr_TriangleCount_Burkhardt || // 1: sum (sum ((A^2) .* A))/6
143 : method == LAGr_TriangleCount_Cohen || // 2: sum (sum ((L * U) .*A))/2
144 : method == LAGr_TriangleCount_Sandia_LL || // 3: sum (sum ((L * L) .* L))
145 : method == LAGr_TriangleCount_Sandia_UU || // 4: sum (sum ((U * U) .* U))
146 : method == LAGr_TriangleCount_Sandia_LUT || // 5: sum (sum ((L * U') .* L))
147 : method == LAGr_TriangleCount_Sandia_ULT, // 6: sum (sum ((U * L') .* U))
148 : GrB_INVALID_VALUE, "method is invalid") ;
149 :
150 : // get the presort
151 : LAGr_TriangleCount_Presort presort ;
152 5846 : presort = (p_presort == NULL) ? LAGr_TriangleCount_AutoSort : (*p_presort) ;
153 5846 : LG_ASSERT_MSG (
154 : presort == LAGr_TriangleCount_NoSort ||
155 : presort == LAGr_TriangleCount_Ascending ||
156 : presort == LAGr_TriangleCount_Descending ||
157 : presort == LAGr_TriangleCount_AutoSort,
158 : GrB_INVALID_VALUE, "presort is invalid") ;
159 :
160 5836 : LG_TRY (LAGraph_CheckGraph (G, msg)) ;
161 5836 : LG_ASSERT (ntriangles != NULL, GrB_NULL_POINTER) ;
162 5836 : LG_ASSERT (G->nself_edges == 0, LAGRAPH_NO_SELF_EDGES_ALLOWED) ;
163 :
164 5836 : LG_ASSERT_MSG ((G->kind == LAGraph_ADJACENCY_UNDIRECTED ||
165 : (G->kind == LAGraph_ADJACENCY_DIRECTED &&
166 : G->is_symmetric_structure == LAGraph_TRUE)),
167 : LAGRAPH_SYMMETRIC_STRUCTURE_REQUIRED,
168 : "G->A must be known to be symmetric") ;
169 :
170 5836 : if (method == LAGr_TriangleCount_AutoMethod)
171 : {
172 : // AutoMethod: use default, Sandia_LUT: sum (sum ((L * U') .* L))
173 1433 : method = LAGr_TriangleCount_Sandia_LUT ;
174 : }
175 :
176 : // only the Sandia_* methods can benefit from the presort
177 5836 : bool method_can_use_presort =
178 5118 : method == LAGr_TriangleCount_Sandia_LL || // sum (sum ((L * L) .* L))
179 4402 : method == LAGr_TriangleCount_Sandia_UU || // sum (sum ((U * U) .* U))
180 10954 : method == LAGr_TriangleCount_Sandia_LUT || // sum (sum ((L * U') .* L))
181 : method == LAGr_TriangleCount_Sandia_ULT ; // sum (sum ((U * L') .* U))
182 :
183 5836 : GrB_Matrix A = G->A ;
184 5836 : GrB_Vector Degree = G->out_degree ;
185 :
186 5836 : bool auto_sort = (presort == LAGr_TriangleCount_AutoSort) ;
187 5836 : if (auto_sort && method_can_use_presort)
188 : {
189 1796 : LG_ASSERT_MSG (Degree != NULL,
190 : LAGRAPH_NOT_CACHED, "G->out_degree is required") ;
191 : }
192 :
193 : //--------------------------------------------------------------------------
194 : // initializations
195 : //--------------------------------------------------------------------------
196 :
197 : GrB_Index n ;
198 5832 : GRB_TRY (GrB_Matrix_nrows (&n, A)) ;
199 5832 : GRB_TRY (GrB_Matrix_new (&C, GrB_INT64, n, n)) ;
200 5162 : GrB_Semiring semiring = LAGraph_plus_one_int64 ;
201 5162 : GrB_Monoid monoid = GrB_PLUS_MONOID_INT64 ;
202 :
203 : //--------------------------------------------------------------------------
204 : // heuristic sort rule
205 : //--------------------------------------------------------------------------
206 :
207 5162 : if (!method_can_use_presort)
208 : {
209 : // no sorting for the Burkhardt and Cohen methods: presort parameter
210 : // is ignored.
211 929 : presort = LAGr_TriangleCount_NoSort ;
212 : }
213 4233 : else if (auto_sort)
214 : {
215 : // auto selection of sorting method for Sandia_* methods
216 1602 : presort = LAGr_TriangleCount_NoSort ; // default is not to sort
217 :
218 1602 : if (method_can_use_presort)
219 : {
220 : // This rule is very similar to Scott Beamer's rule in the GAP TC
221 : // benchmark, except that it is extended to handle the ascending
222 : // sort needed by methods 3 and 5. It also uses a stricter rule,
223 : // since the performance of triangle counting in SuiteSparse:
224 : // GraphBLAS is less sensitive to the sorting as compared to the
225 : // GAP algorithm. This is because the dot products in SuiteSparse:
226 : // GraphBLAS use binary search if one vector is very sparse
227 : // compared to the other. As a result, SuiteSparse:GraphBLAS needs
228 : // the sort for fewer matrices, as compared to the GAP algorithm.
229 :
230 : // With this rule, the GAP-kron and GAP-twitter matrices are
231 : // sorted, and the others remain unsorted. With the rule in the
232 : // GAP tc.cc benchmark, GAP-kron and GAP-twitter are sorted, and so
233 : // is GAP-web, but GAP-web is not sorted here.
234 :
235 : #define NSAMPLES 2000
236 : GrB_Index nvals ;
237 1614 : GRB_TRY (GrB_Matrix_nvals (&nvals, A)) ;
238 1602 : if (n > NSAMPLES && ((double) nvals / ((double) n)) >= 10)
239 : {
240 : // estimate the mean and median degrees
241 : double mean, median ;
242 479 : LG_TRY (LAGr_SampleDegree (&mean, &median,
243 : G, true, NSAMPLES, n, msg)) ;
244 : // sort if the average degree is very high vs the median
245 467 : if (mean > 3 * median)
246 : {
247 292 : switch (method)
248 : {
249 40 : case LAGr_TriangleCount_Sandia_LL:
250 : // 3:sum (sum ((L * L) .* L))
251 40 : presort = LAGr_TriangleCount_Ascending ;
252 40 : break ;
253 38 : case LAGr_TriangleCount_Sandia_UU:
254 : // 4: sum (sum ((U * U) .* U))
255 38 : presort = LAGr_TriangleCount_Descending ;
256 38 : break ;
257 172 : default:
258 : case LAGr_TriangleCount_Sandia_LUT:
259 : // 5: sum (sum ((L * U') .* L))
260 172 : presort = LAGr_TriangleCount_Ascending ;
261 172 : break ;
262 42 : case LAGr_TriangleCount_Sandia_ULT:
263 : // 6: sum (sum ((U * L') .* U))
264 42 : presort = LAGr_TriangleCount_Descending ;
265 42 : break ;
266 : }
267 : }
268 : }
269 : }
270 : }
271 :
272 : //--------------------------------------------------------------------------
273 : // sort the input matrix, if requested
274 : //--------------------------------------------------------------------------
275 :
276 5150 : if (presort != LAGr_TriangleCount_NoSort)
277 : {
278 : // P = permutation that sorts the rows by their degree
279 1888 : LG_TRY (LAGr_SortByDegree (&P, G, true,
280 : presort == LAGr_TriangleCount_Ascending, msg)) ;
281 :
282 : // T = A (P,P) and typecast to boolean
283 1709 : GRB_TRY (GrB_Matrix_new (&T, GrB_BOOL, n, n)) ;
284 1538 : GRB_TRY (GrB_extract (T, NULL, NULL, A, (GrB_Index *) P, n,
285 : (GrB_Index *) P, n, NULL)) ;
286 1184 : A = T ;
287 :
288 : // free workspace
289 1184 : LG_TRY (LAGraph_Free ((void **) &P, NULL)) ;
290 : }
291 :
292 : //--------------------------------------------------------------------------
293 : // count triangles
294 : //--------------------------------------------------------------------------
295 :
296 : int64_t ntri ;
297 :
298 4446 : switch (method)
299 : {
300 :
301 211 : case LAGr_TriangleCount_Burkhardt: // 1: sum (sum ((A^2) .* A)) / 6
302 :
303 211 : GRB_TRY (GrB_mxm (C, A, NULL, semiring, A, A, GrB_DESC_S)) ;
304 61 : GRB_TRY (GrB_reduce (&ntri, NULL, monoid, C, NULL)) ;
305 61 : ntri /= 6 ;
306 61 : break ;
307 :
308 718 : case LAGr_TriangleCount_Cohen: // 2: sum (sum ((L * U) .* A)) / 2
309 :
310 718 : LG_TRY (tricount_prep (&L, &U, A, msg)) ;
311 214 : GRB_TRY (GrB_mxm (C, A, NULL, semiring, L, U, GrB_DESC_S)) ;
312 61 : GRB_TRY (GrB_reduce (&ntri, NULL, monoid, C, NULL)) ;
313 61 : ntri /= 2 ;
314 61 : break ;
315 :
316 489 : case LAGr_TriangleCount_Sandia_LL: // 3: sum (sum ((L * L) .* L))
317 :
318 : // using the masked saxpy3 method
319 489 : LG_TRY (tricount_prep (&L, NULL, A, msg)) ;
320 241 : GRB_TRY (GrB_mxm (C, L, NULL, semiring, L, L, GrB_DESC_S)) ;
321 62 : GRB_TRY (GrB_reduce (&ntri, NULL, monoid, C, NULL)) ;
322 62 : break ;
323 :
324 489 : case LAGr_TriangleCount_Sandia_UU: // 4: sum (sum ((U * U) .* U))
325 :
326 : // using the masked saxpy3 method
327 489 : LG_TRY (tricount_prep (NULL, &U, A, msg)) ;
328 241 : GRB_TRY (GrB_mxm (C, U, NULL, semiring, U, U, GrB_DESC_S)) ;
329 62 : GRB_TRY (GrB_reduce (&ntri, NULL, monoid, C, NULL)) ;
330 62 : break ;
331 :
332 1837 : default:
333 : case LAGr_TriangleCount_Sandia_LUT: // 5: sum (sum ((L * U') .* L))
334 :
335 : // This tends to be the fastest method for most large matrices, but
336 : // the Sandia_ULT method is also very fast.
337 :
338 : // using the masked dot product
339 1837 : LG_TRY (tricount_prep (&L, &U, A, msg)) ;
340 554 : GRB_TRY (GrB_mxm (C, L, NULL, semiring, L, U, GrB_DESC_ST1)) ;
341 164 : GRB_TRY (GrB_reduce (&ntri, NULL, monoid, C, NULL)) ;
342 164 : break ;
343 :
344 702 : case LAGr_TriangleCount_Sandia_ULT: // 6: sum (sum ((U * L') .* U))
345 :
346 : // using the masked dot product
347 702 : LG_TRY (tricount_prep (&L, &U, A, msg)) ;
348 208 : GRB_TRY (GrB_mxm (C, U, NULL, semiring, U, L, GrB_DESC_ST1)) ;
349 62 : GRB_TRY (GrB_reduce (&ntri, NULL, monoid, C, NULL)) ;
350 62 : break ;
351 : }
352 :
353 : //--------------------------------------------------------------------------
354 : // return result
355 : //--------------------------------------------------------------------------
356 :
357 472 : LG_FREE_ALL ;
358 472 : if (p_method != NULL) (*p_method) = method ;
359 472 : if (p_presort != NULL) (*p_presort) = presort ;
360 472 : (*ntriangles) = (uint64_t) ntri ;
361 472 : return (GrB_SUCCESS) ;
362 : }
|
