Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_scc.c
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 Yongzhe Zhang (zyz915@gmail.com)
15 :
16 : //------------------------------------------------------------------------------
17 :
18 : // TODO: not ready for src; uses global variables
19 :
20 : /**
21 : * Code is based on the Min-Label algorithm described in the following paper:
22 : * D. Yan, J. Cheng, K. Xin, Y. Lu, W. Ng, Y. Bu, "Pregel Algorithms for Graph
23 : * Connectivity Problems with Performance Guarantees"
24 : * Proc. VLDB Endow. 7, 14 (October 2014), 1821–1832.
25 : * DOI: https://doi.org/10.14778/2733085.2733089
26 : **/
27 :
28 : #define LG_FREE_ALL ;
29 :
30 : #include "LG_internal.h"
31 : #include <LAGraph.h>
32 : #include <LAGraphX.h>
33 :
34 : #if LAGRAPH_SUITESPARSE
35 :
36 : //****************************************************************************
37 : //arrays used in SelectOp
38 : typedef struct
39 : {
40 : uint64_t *F, *B;
41 : bool *M;
42 : } LG_SCC_Context;
43 : #define SCCCONTEXT \
44 : "typedef struct \n" \
45 : "{\n" \
46 : " uint64_t *F, *B;\n" \
47 : " bool *M;\n" \
48 : "} LG_SCC_Context;\n"
49 :
50 :
51 : // LG_SCC_edge_removal:
52 : // - remove the edges connected to newly identified SCCs (vertices u with M[u]==1)
53 : // - remove the edges (u, v) where u and v can never be in the same SCC.
54 : //
55 : // Here's a brief explanation of the second case. After the forward and backward
56 : // propagation, each vertex u has two labels
57 : // - F[u]: the smallest vertex that can reach u
58 : // - B[u]: the smallest vertex that is reachable from u
59 : // If two vertices u and v are in the same SCC, then F[u]==F[v] and B[u]==B[v] must
60 : // hold. The converse is not true unless F[u]==B[u]. However, we can safely remove
61 : // an edge (u, v) if either F[u]!=F[v] or B[u]!=B[v] holds, which can accelerate
62 : // the SCC computation in the future rounds.
63 :
64 : void LG_SCC_edge_removal (bool *z, const void *x, GrB_Index i, GrB_Index j, const LG_SCC_Context *thunk) ;
65 1474596 : void LG_SCC_edge_removal (bool *z, const void *x, GrB_Index i, GrB_Index j, const LG_SCC_Context *thunk)
66 : {
67 1194 : (*z) = (!thunk->M[i] && !thunk->M[j]
68 974 : && thunk->F[i] == thunk->F[j]
69 1475790 : && thunk->B[i] == thunk->B[j]) ;
70 1474596 : }
71 : #define EDGE_REMOVAL \
72 : "void LG_SCC_edge_removal \n" \
73 : "(bool *z, const void *x, GrB_Index i, GrB_Index j, const LG_SCC_Context *thunk)\n" \
74 : "{\n" \
75 : " (*z) = (!thunk->M[i] && !thunk->M[j] \n" \
76 : " && thunk->F[i] == thunk->F[j] \n" \
77 : " && thunk->B[i] == thunk->B[j]) ;\n" \
78 : "}\n"
79 :
80 : //****************************************************************************
81 : // LG_SCC_trim_one: remove the edges connected to trivial SCCs
82 : // - A vertex is a trivial SCC if it has no incoming or outgoing edges.
83 : // - M[i] = i | if vertex i is a trivial SCC
84 : // M[i] = n | otherwise
85 :
86 : void LG_SCC_trim_one (bool *z, const void *x, GrB_Index i, GrB_Index j, const LG_SCC_Context *thunk) ;
87 127480 : void LG_SCC_trim_one (bool *z, const void *x, GrB_Index i, GrB_Index j, const LG_SCC_Context *thunk)
88 : {
89 127480 : (*z) = (thunk->F[i] == thunk->F[j]) ;
90 127480 : }
91 : #define TRIM_ONE \
92 : "void LG_SCC_trim_one\n" \
93 : "(bool *z, const void *x, GrB_Index i, GrB_Index j, const LG_SCC_Context *thunk)\n" \
94 : "{\n" \
95 : " (*z) = (thunk->F[i] == thunk->F[j]) ;\n" \
96 : "}\n"
97 :
98 : //****************************************************************************
99 : // label propagation
100 : // - label : (input/output) labels
101 : // - mask : (input) mask
102 : // - A : (input) original matrix
103 : // - AT : (input) transposed matrix
104 : // - n : (input) number of vertices
105 :
106 : #undef LG_FREE_ALL
107 : #define LG_FREE_ALL \
108 : { \
109 : GrB_free (&s) ; \
110 : GrB_free (&t) ; \
111 : }
112 :
113 296 : static GrB_Info propagate (GrB_Vector label, GrB_Vector mask,
114 : const GrB_Matrix A, const GrB_Matrix AT, GrB_Index n, char *msg)
115 : {
116 : GrB_Info info;
117 : // semirings
118 296 : GrB_Vector s = NULL, t = NULL;
119 296 : GRB_TRY (GrB_Vector_new (&s, GrB_UINT64, n));
120 296 : GRB_TRY (GrB_Vector_new (&t, GrB_UINT64, n));
121 296 : GRB_TRY (GrB_assign (s, mask, 0, label, GrB_ALL, 0, 0));
122 : // GxB_fprint(s, GxB_SHORT, stdout);
123 296 : GRB_TRY (GrB_assign (t, 0, 0, label, GrB_ALL, 0, 0));
124 296 : GRB_TRY (GrB_wait(A, GrB_MATERIALIZE));
125 :
126 : bool active;
127 : while (true)
128 : {
129 : // GRB_TRY (GrB_mxv
130 : // (t, 0, GrB_MIN_UINT64, GrB_MIN_SECOND_SEMIRING_UINT64, AT, s, 0));
131 14184 : GRB_TRY (GrB_vxm (t, 0, GrB_MIN_UINT64,
132 : GrB_MIN_FIRST_SEMIRING_UINT64, s, A, 0));
133 14184 : GRB_TRY (GrB_eWiseMult (mask, 0, 0, GrB_NE_UINT64, t, label, 0));
134 14184 : GRB_TRY (GrB_assign (label, NULL, NULL, t, GrB_ALL, n, NULL));
135 14184 : GRB_TRY (GrB_reduce (&active, 0, GrB_LOR_MONOID_BOOL, mask, 0));
136 14184 : if (!active) break;
137 13888 : GRB_TRY (GrB_Vector_clear(s));
138 13888 : GRB_TRY (GrB_assign (s, mask, 0, label, GrB_ALL, 0, 0));
139 13888 : GRB_TRY (GrB_wait(s, GrB_MATERIALIZE));
140 : }
141 :
142 296 : LG_FREE_ALL ;
143 296 : return GrB_SUCCESS;
144 : }
145 : //****************************************************************************
146 :
147 : #undef LG_FREE_ALL
148 : #define LG_FREE_ALL \
149 : LAGraph_Free ((void **) &contx.F, msg); \
150 : LAGraph_Free ((void **) &contx.B, msg); \
151 : LAGraph_Free ((void **) &contx.M, msg); \
152 : GrB_free (&ind); \
153 : GrB_free (&inf); \
154 : GrB_free (&f); \
155 : GrB_free (&b); \
156 : GrB_free (&D); \
157 : GrB_free (&x); \
158 : GrB_free (&mask); \
159 : GrB_free (&m2); \
160 : GrB_free (&FW); \
161 : GrB_free (&BW); \
162 : GrB_free (&sel1); \
163 : GrB_free (&sel2); \
164 : GrB_free (&contx_type); \
165 : GrB_free (&scc);
166 :
167 : #endif
168 :
169 : //****************************************************************************
170 104 : int LAGraph_scc
171 : (
172 : GrB_Vector *result, // output: array of component identifiers
173 : GrB_Matrix A, // input matrix
174 : char *msg
175 : )
176 : {
177 : #if LAGRAPH_SUITESPARSE
178 :
179 104 : LG_CLEAR_MSG ;
180 104 : LG_SCC_Context contx = {NULL, NULL, NULL};
181 104 : GrB_Info info = GrB_SUCCESS;
182 104 : GrB_Type contx_type = NULL;
183 104 : GrB_Type type_F = NULL, type_B = NULL, type_M = NULL;
184 104 : int hand_F = GrB_DEFAULT, hand_B = GrB_DEFAULT, hand_M = GrB_DEFAULT;
185 104 : uint64_t n_F = 0, n_B = 0, n_M = 0, size_F = 0, size_B = 0, size_M = 0;
186 104 : GrB_Vector scc = NULL ;
187 104 : GrB_Vector ind = NULL ;
188 104 : GrB_Vector inf = NULL ;
189 104 : GrB_Vector x = NULL ;
190 104 : GrB_Vector f = NULL, b = NULL, mask = NULL, m2 = NULL;
191 104 : GrB_IndexUnaryOp sel1 = NULL, sel2 = NULL ;
192 104 : GrB_Monoid Add = NULL ;
193 104 : GrB_Matrix FW = NULL, BW = NULL, D = NULL;
194 104 : LG_ASSERT(result != NULL, GrB_NULL_POINTER);
195 103 : LG_ASSERT(A != NULL, GrB_NULL_POINTER);
196 :
197 : GrB_Index n, ncols, nvals;
198 103 : GRB_TRY (GrB_Matrix_nrows (&n, A));
199 103 : GRB_TRY (GrB_Matrix_ncols (&ncols, A));
200 103 : LG_ASSERT(n == ncols, GrB_DIMENSION_MISMATCH);
201 :
202 : #if !LG_SUITESPARSE_GRAPHBLAS_V10
203 : LG_TRY (LAGraph_Malloc ((void **) &contx.F, n, sizeof (uint64_t), msg)) ;
204 : LG_TRY (LAGraph_Malloc ((void **) &contx.B, n, sizeof (uint64_t), msg)) ;
205 : LG_TRY (LAGraph_Malloc ((void **) &contx.M, n, sizeof (bool), msg)) ;
206 : #endif
207 : // scc: the SCC identifier for each vertex
208 : // scc[u] == n: not assigned yet
209 102 : GRB_TRY (GrB_Vector_new (&scc, GrB_UINT64, n));
210 : // vector of indices: ind[i] == i
211 102 : GRB_TRY (GrB_Vector_new (&ind, GrB_UINT64, n));
212 102 : GRB_TRY (GrB_Vector_assign_UINT64 (
213 : ind, NULL, NULL, (uint64_t) 0, GrB_ALL, n, NULL)) ;
214 102 : GRB_TRY (GrB_Vector_apply_IndexOp_UINT64 (
215 : ind, NULL, NULL, GrB_ROWINDEX_INT64, ind, 0, NULL)) ;
216 : // vector of infinite value: inf[i] == n
217 102 : GRB_TRY (GrB_Vector_new (&inf, GrB_UINT64, n));
218 102 : GRB_TRY (GrB_assign (inf, NULL, NULL, n, GrB_ALL, 0, NULL));
219 : // other vectors
220 102 : GRB_TRY (GrB_Vector_new (&f, GrB_UINT64, n));
221 102 : GRB_TRY (GrB_Vector_new (&b, GrB_UINT64, n));
222 102 : GRB_TRY (GrB_Vector_new (&mask, GrB_BOOL, n));
223 102 : GRB_TRY (GrB_Vector_new (&m2, GrB_BOOL, n));
224 102 : GRB_TRY (GrB_Vector_new (&x, GrB_BOOL, n));
225 102 : GRB_TRY (GxB_Type_new (
226 : &contx_type, sizeof(LG_SCC_Context), "LG_SCC_Context", SCCCONTEXT)) ;
227 :
228 102 : GRB_TRY (GxB_IndexUnaryOp_new (
229 : &sel1, (GxB_index_unary_function) LG_SCC_trim_one,
230 : GrB_BOOL, GrB_UINT64, contx_type,
231 : // NULL, NULL
232 : "LG_SCC_trim_one", TRIM_ONE
233 : ));
234 102 : GRB_TRY (GxB_IndexUnaryOp_new (
235 : &sel2, (GxB_index_unary_function) LG_SCC_edge_removal,
236 : GrB_BOOL, GrB_UINT64, contx_type,
237 : // NULL, NULL
238 : "LG_SCC_edge_removal", EDGE_REMOVAL
239 : ));
240 :
241 : // store the graph in both directions (forward / backward)
242 102 : GRB_TRY (GrB_Matrix_new (&FW, GrB_BOOL, n, n));
243 102 : GRB_TRY (GrB_Matrix_new (&BW, GrB_BOOL, n, n));
244 102 : GRB_TRY (GrB_Matrix_assign_BOOL(
245 : FW, A, NULL, true, GrB_ALL, n, GrB_ALL, n, GrB_DESC_S)) ;
246 102 : GRB_TRY (GrB_transpose (BW, NULL, NULL, FW, NULL)); // BW = FW'
247 :
248 : // check format
249 : int32_t A_format, AT_format;
250 102 : GRB_TRY (GrB_get (FW, &A_format , GrB_STORAGE_ORIENTATION_HINT));
251 102 : GRB_TRY (GrB_get (BW, &AT_format, GrB_STORAGE_ORIENTATION_HINT));
252 :
253 102 : bool is_csr = (A_format == GrB_ROWMAJOR && AT_format == GrB_ROWMAJOR);
254 102 : LG_ASSERT (is_csr, GrB_INVALID_VALUE) ;
255 :
256 : // remove trivial SCCs
257 102 : GRB_TRY (GrB_Vector_assign_BOOL (x, NULL, NULL, 0, GrB_ALL, n, NULL)) ;
258 102 : GRB_TRY (GrB_mxv (m2, NULL, NULL, GxB_ANY_PAIR_BOOL, FW, x, NULL)) ;
259 102 : GRB_TRY (GrB_mxv (mask, m2, NULL, GxB_ANY_PAIR_BOOL, BW, x, GrB_DESC_S)) ;
260 102 : GRB_TRY (GrB_Vector_nvals (&nvals, mask));
261 :
262 102 : GRB_TRY (GrB_assign (scc, NULL, NULL, ind, GrB_ALL, 0, NULL));
263 102 : GRB_TRY (GrB_assign (scc, mask, NULL, n, GrB_ALL, 0, NULL));
264 :
265 102 : if (nvals < n)
266 : {
267 : // No reason for context.
268 : #if LG_SUITESPARSE_GRAPHBLAS_V10
269 18 : GRB_TRY(GxB_Vector_unload(
270 : scc, (void **) &contx.F, &type_F, &n_F, &size_F, &hand_F, NULL)) ;
271 : #else
272 : GRB_TRY (GrB_Vector_extractTuples_UINT64 (NULL, contx.F, &n, scc));
273 : #endif
274 18 : GRB_TRY (GrB_Matrix_select_UDT (
275 : FW, NULL, NULL, sel1, FW, &contx, NULL)) ;
276 18 : GRB_TRY (GrB_Matrix_select_UDT (
277 : BW, NULL, NULL, sel1, BW, &contx, NULL)) ;
278 : #if LG_SUITESPARSE_GRAPHBLAS_V10
279 18 : GRB_TRY(GxB_Vector_load(
280 : scc, (void **) &contx.F, type_F, n_F, size_F, hand_F, NULL)) ;
281 : #endif
282 : }
283 :
284 102 : GRB_TRY (GrB_Matrix_nvals (&nvals, FW));
285 250 : while (nvals > 0)
286 : {
287 148 : GRB_TRY (GrB_Vector_apply_BinaryOp2nd_UINT64 (
288 : mask, NULL, NULL, GrB_EQ_UINT64, scc, n, NULL));
289 148 : GRB_TRY (GrB_assign (f, NULL, NULL, ind, GrB_ALL, 0, NULL));
290 148 : LG_TRY (propagate (f, mask, FW, BW, n, msg));
291 :
292 148 : GRB_TRY (GrB_eWiseMult (mask, NULL, NULL, GrB_EQ_UINT64, f, ind, NULL));
293 148 : GRB_TRY (GrB_Vector_assign_UINT64 (
294 : b, NULL, NULL, n, GrB_ALL, 0, NULL)) ;
295 148 : GRB_TRY (GrB_assign (b, mask, NULL, ind, GrB_ALL, 0, NULL));
296 148 : LG_TRY (propagate (b, mask, BW, FW, n, msg));
297 :
298 148 : GRB_TRY (GrB_eWiseMult (mask, NULL, NULL, GrB_EQ_UINT64, f, b, NULL));
299 148 : GRB_TRY (GrB_assign (scc, mask, GrB_MIN_UINT64, f, GrB_ALL, 0, NULL));
300 :
301 : #if LG_SUITESPARSE_GRAPHBLAS_V10
302 148 : GRB_TRY(GxB_Vector_unload(
303 : f, (void **) &contx.F, &type_F, &n_F, &size_F, &hand_F, NULL)) ;
304 148 : GRB_TRY(GxB_Vector_unload(
305 : b, (void **) &contx.B, &type_B, &n_B, &size_B, &hand_B, NULL)) ;
306 148 : GRB_TRY(GxB_Vector_unload(
307 : mask, (void **) &contx.M, &type_M, &n_M, &size_M, &hand_M, NULL
308 : )) ;
309 : #else
310 : GRB_TRY (GrB_Vector_extractTuples_UINT64 (NULL, contx.F, &n, f));
311 : GRB_TRY (GrB_Vector_extractTuples_UINT64 (NULL, contx.B, &n, b));
312 : GRB_TRY (GrB_Vector_extractTuples_BOOL (NULL, contx.M, &n, mask));
313 : #endif
314 :
315 148 : GRB_TRY (GrB_Matrix_select_UDT (
316 : FW, NULL, NULL, sel2, FW, &contx, NULL)) ;
317 148 : GRB_TRY (GrB_Matrix_select_UDT (
318 : BW, NULL, NULL, sel2, BW, &contx, NULL)) ;
319 : #if LG_SUITESPARSE_GRAPHBLAS_V10
320 148 : GRB_TRY(GxB_Vector_load(
321 : f, (void **) &contx.F, type_F, n_F, size_F, hand_F, NULL)) ;
322 148 : GRB_TRY(GxB_Vector_load(
323 : b, (void **) &contx.B, type_B, n_B, size_B, hand_B, NULL)) ;
324 148 : GRB_TRY(GxB_Vector_load(
325 : mask, (void **) &contx.M, type_M, n_M, size_M, hand_M, NULL
326 : )) ;
327 : #endif
328 148 : GRB_TRY (GrB_Matrix_nvals (&nvals, FW));
329 : }
330 102 : GRB_TRY (GrB_Vector_apply_BinaryOp2nd_UINT64 (
331 : mask, NULL, NULL, GrB_EQ_UINT64, scc, n, NULL));
332 102 : GRB_TRY (GrB_assign (scc, mask, NULL, ind, GrB_ALL, 0, NULL));
333 :
334 102 : *result = scc;
335 102 : scc = NULL;
336 :
337 102 : LG_FREE_ALL ;
338 102 : return (GrB_SUCCESS) ;
339 : #else
340 : return (GrB_NOT_IMPLEMENTED) ;
341 : #endif
342 : }
|
