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 : // global C arrays used in SelectOp
38 : GrB_Index *I = NULL, *V = NULL, *F = NULL, *B = NULL, *M = NULL;
39 :
40 : // edge_removal:
41 : // - remove the edges connected to newly identified SCCs (vertices u with M[u]==1)
42 : // - remove the edges (u, v) where u and v can never be in the same SCC.
43 : //
44 : // Here's a brief explanation of the second case. After the forward and backward
45 : // propagation, each vertex u has two labels
46 : // - F[u]: the smallest vertex that can reach u
47 : // - B[u]: the smallest vertex that is reachable from u
48 : // If two vertices u and v are in the same SCC, then F[u]==F[v] and B[u]==B[v] must
49 : // hold. The converse is not true unless F[u]==B[u]. However, we can safely remove
50 : // an edge (u, v) if either F[u]!=F[v] or B[u]!=B[v] holds, which can accelerate
51 : // the SCC computation in the future rounds.
52 :
53 : void edge_removal (bool *z, const void *x, GrB_Index i, GrB_Index j, const void *thunk) ;
54 1515400 : void edge_removal (bool *z, const void *x, GrB_Index i, GrB_Index j, const void *thunk)
55 : {
56 1515400 : (*z) = (!M[i] && !M[j] && F[i] == F[j] && B[i] == B[j]) ;
57 1515400 : }
58 :
59 : //****************************************************************************
60 : // trim_one: remove the edges connected to trivial SCCs
61 : // - A vertex is a trivial SCC if it has no incoming or outgoing edges.
62 : // - M[i] = i | if vertex i is a trivial SCC
63 : // M[i] = n | otherwise
64 :
65 : void trim_one (bool *z, const void *x, GrB_Index i, GrB_Index j, const void *thunk) ;
66 127480 : void trim_one (bool *z, const void *x, GrB_Index i, GrB_Index j, const void *thunk)
67 : {
68 127480 : (*z) = (M[i] == M[j]) ;
69 127480 : }
70 :
71 : //****************************************************************************
72 : // label propagation
73 : // - label : (input/output) labels
74 : // - mask : (input) mask
75 : // - A : (input) original matrix
76 : // - AT : (input) transposed matrix
77 : // - n : (input) number of vertices
78 :
79 : #undef LG_FREE_ALL
80 : #define LG_FREE_ALL \
81 : GrB_free (&s) ; \
82 : GrB_free (&t) ;
83 :
84 160 : static GrB_Info propagate (GrB_Vector label, GrB_Vector mask,
85 : GrB_Matrix A, GrB_Matrix AT, GrB_Index n, char *msg)
86 : {
87 : GrB_Info info;
88 : // semirings
89 :
90 160 : GrB_Vector s = NULL, t = NULL;
91 160 : GRB_TRY (GrB_Vector_new (&s, GrB_UINT64, n));
92 160 : GRB_TRY (GrB_Vector_new (&t, GrB_UINT64, n));
93 160 : GRB_TRY (GrB_assign (s, mask, 0, label, GrB_ALL, 0, 0));
94 160 : GRB_TRY (GrB_assign (t, 0, 0, label, GrB_ALL, 0, 0));
95 :
96 : GrB_Index active;
97 : while (true)
98 : {
99 7186 : GRB_TRY (GrB_vxm (t, 0, GrB_MIN_UINT64,
100 : GrB_MIN_FIRST_SEMIRING_UINT64, s, A, 0));
101 7186 : GRB_TRY (GrB_eWiseMult (mask, 0, 0, GxB_ISNE_UINT64, t, label, 0));
102 7186 : GRB_TRY (GrB_assign (label, mask, 0, t, GrB_ALL, 0, 0));
103 7186 : GRB_TRY (GrB_reduce (&active, 0, GrB_PLUS_MONOID_UINT64, mask, 0));
104 7186 : if (active == 0) break;
105 7026 : GRB_TRY (GrB_Vector_clear (s));
106 7026 : GRB_TRY (GrB_assign (s, mask, 0, label, GrB_ALL, 0, 0));
107 : }
108 :
109 160 : LG_FREE_ALL ;
110 160 : return GrB_SUCCESS;
111 : }
112 :
113 : //****************************************************************************
114 :
115 : #undef LG_FREE_ALL
116 : #define LG_FREE_ALL \
117 : LAGraph_Free ((void **) &I, msg); \
118 : LAGraph_Free ((void **) &V, msg); \
119 : LAGraph_Free ((void **) &F, msg); \
120 : LAGraph_Free ((void **) &B, msg); \
121 : LAGraph_Free ((void **) &M, msg); \
122 : GrB_free (&ind); \
123 : GrB_free (&inf); \
124 : GrB_free (&f); \
125 : GrB_free (&b); \
126 : GrB_free (&mask); \
127 : GrB_free (&FW); \
128 : GrB_free (&BW); \
129 : GrB_free (&sel1); \
130 : GrB_free (&sel2); \
131 : GrB_free (&scc);
132 :
133 : #endif
134 :
135 : //****************************************************************************
136 53 : int LAGraph_scc
137 : (
138 : GrB_Vector *result, // output: array of component identifiers
139 : GrB_Matrix A, // input matrix
140 : char *msg
141 : )
142 : {
143 : #if LAGRAPH_SUITESPARSE
144 :
145 53 : LG_CLEAR_MSG ;
146 :
147 : GrB_Info info;
148 53 : GrB_Vector scc = NULL ;
149 53 : GrB_Vector ind = NULL ;
150 53 : GrB_Vector inf = NULL ;
151 53 : GrB_Vector f = NULL, b = NULL, mask = NULL ;
152 53 : GrB_IndexUnaryOp sel1 = NULL, sel2 = NULL ;
153 53 : GrB_Monoid Add = NULL ;
154 53 : GrB_Matrix FW = NULL, BW = NULL;
155 :
156 53 : if (result == NULL || A == NULL) return (GrB_NULL_POINTER) ;
157 :
158 : GrB_Index n, ncols, nvals;
159 52 : GRB_TRY (GrB_Matrix_nrows (&n, A));
160 52 : GRB_TRY (GrB_Matrix_ncols (&ncols, A));
161 52 : if (n != ncols) return (GrB_DIMENSION_MISMATCH) ;
162 :
163 : // store the graph in both directions (forward / backward)
164 51 : GRB_TRY (GrB_Matrix_new (&FW, GrB_BOOL, n, n));
165 51 : GRB_TRY (GrB_Matrix_new (&BW, GrB_BOOL, n, n));
166 51 : GRB_TRY (GrB_transpose (FW, 0, 0, A, GrB_DESC_T0)); // FW = A
167 51 : GRB_TRY (GrB_transpose (BW, 0, 0, A, 0)); // BW = A'
168 :
169 : // check format
170 : int32_t A_format, AT_format;
171 51 : GRB_TRY (GrB_get (FW, &A_format , GrB_STORAGE_ORIENTATION_HINT));
172 51 : GRB_TRY (GrB_get (BW, &AT_format, GrB_STORAGE_ORIENTATION_HINT));
173 :
174 51 : bool is_csr = (A_format == GrB_ROWMAJOR && AT_format == GrB_ROWMAJOR);
175 51 : if (!is_csr) return (GrB_INVALID_VALUE) ;
176 :
177 51 : LG_TRY (LAGraph_Malloc ((void **) &I, n, sizeof (GrB_Index), msg)) ;
178 51 : LG_TRY (LAGraph_Malloc ((void **) &V, n, sizeof (GrB_Index), msg)) ;
179 51 : LG_TRY (LAGraph_Malloc ((void **) &F, n, sizeof (GrB_Index), msg)) ;
180 51 : LG_TRY (LAGraph_Malloc ((void **) &B, n, sizeof (GrB_Index), msg)) ;
181 51 : LG_TRY (LAGraph_Malloc ((void **) &M, n, sizeof (GrB_Index), msg)) ;
182 :
183 19605 : for (GrB_Index i = 0; i < n; i++)
184 19554 : I[i] = V[i] = i;
185 :
186 : // scc: the SCC identifier for each vertex
187 : // scc[u] == n: not assigned yet
188 51 : GRB_TRY (GrB_Vector_new (&scc, GrB_UINT64, n));
189 : // vector of indices: ind[i] == i
190 51 : GRB_TRY (GrB_Vector_new (&ind, GrB_UINT64, n));
191 51 : GRB_TRY (GrB_Vector_build (ind, I, V, n, GrB_PLUS_UINT64));
192 : // vector of infinite value: inf[i] == n
193 51 : GRB_TRY (GrB_Vector_new (&inf, GrB_UINT64, n));
194 51 : GRB_TRY (GrB_assign (inf, 0, 0, n, GrB_ALL, 0, 0));
195 : // other vectors
196 51 : GRB_TRY (GrB_Vector_new (&f, GrB_UINT64, n));
197 51 : GRB_TRY (GrB_Vector_new (&b, GrB_UINT64, n));
198 51 : GRB_TRY (GrB_Vector_new (&mask, GrB_UINT64, n));
199 51 : GRB_TRY (GrB_IndexUnaryOp_new (&sel1, (void *) trim_one, GrB_BOOL, GrB_UINT64, GrB_UINT64));
200 51 : GRB_TRY (GrB_IndexUnaryOp_new (&sel2, (void *) edge_removal, GrB_BOOL, GrB_UINT64, GrB_UINT64));
201 :
202 : // remove trivial SCCs
203 51 : GRB_TRY (GrB_reduce (f, 0, GrB_PLUS_UINT64, GrB_PLUS_UINT64, FW, 0));
204 51 : GRB_TRY (GrB_reduce (b, 0, GrB_PLUS_UINT64, GrB_PLUS_UINT64, BW, 0));
205 51 : GRB_TRY (GrB_eWiseMult (mask, 0, GxB_LAND_UINT64, GxB_LAND_UINT64, f, b, 0));
206 51 : GRB_TRY (GrB_Vector_nvals (&nvals, mask));
207 :
208 51 : GRB_TRY (GrB_assign (scc, 0, 0, ind, GrB_ALL, 0, 0));
209 51 : GRB_TRY (GrB_assign (scc, mask, 0, n, GrB_ALL, 0, 0));
210 51 : GRB_TRY (GrB_Vector_clear (mask));
211 :
212 51 : if (nvals < n)
213 : {
214 9 : GRB_TRY (GrB_Vector_extractTuples (I, M, &n, scc));
215 9 : GRB_TRY (GrB_select (FW, 0, 0, sel1, FW, 0, 0));
216 9 : GRB_TRY (GrB_select (BW, 0, 0, sel1, BW, 0, 0));
217 : }
218 :
219 51 : GRB_TRY (GrB_Matrix_nvals (&nvals, FW));
220 131 : while (nvals > 0)
221 : {
222 80 : GRB_TRY (GrB_eWiseMult (mask, 0, 0, GxB_ISEQ_UINT64, scc, inf, 0));
223 80 : GRB_TRY (GrB_assign (f, 0, 0, ind, GrB_ALL, 0, 0));
224 80 : LG_TRY (propagate (f, mask, FW, BW, n, msg));
225 :
226 80 : GRB_TRY (GrB_eWiseMult (mask, 0, 0, GxB_ISEQ_UINT64, f, ind, 0));
227 80 : GRB_TRY (GrB_assign (b, 0, 0, inf, GrB_ALL, 0, 0));
228 80 : GRB_TRY (GrB_assign (b, mask, 0, ind, GrB_ALL, 0, 0));
229 80 : LG_TRY (propagate (b, mask, BW, FW, n, msg));
230 :
231 80 : GRB_TRY (GrB_eWiseMult (mask, 0, 0, GxB_ISEQ_UINT64, f, b, 0));
232 80 : GRB_TRY (GrB_assign (scc, mask, GrB_MIN_UINT64, f, GrB_ALL, 0, 0));
233 :
234 80 : GRB_TRY (GrB_Vector_extractTuples (I, F, &n, f));
235 80 : GRB_TRY (GrB_Vector_extractTuples (I, B, &n, b));
236 80 : GRB_TRY (GrB_Vector_extractTuples (I, M, &n, mask));
237 :
238 80 : GRB_TRY (GrB_select (FW, 0, 0, sel2, FW, 0, 0));
239 80 : GRB_TRY (GrB_select (BW, 0, 0, sel2, BW, 0, 0));
240 :
241 80 : GRB_TRY (GrB_Matrix_nvals (&nvals, FW));
242 : }
243 51 : GRB_TRY (GrB_eWiseMult (mask, 0, 0, GxB_ISEQ_UINT64, scc, inf, 0));
244 51 : GRB_TRY (GrB_assign (scc, mask, 0, ind, GrB_ALL, 0, 0));
245 :
246 51 : GRB_TRY (GrB_eWiseMult (mask, 0, 0, GxB_ISEQ_UINT64, scc, ind, 0));
247 51 : GRB_TRY (GrB_reduce (&nvals, 0, GrB_PLUS_MONOID_UINT64, mask, 0));
248 :
249 51 : *result = scc;
250 51 : scc = NULL;
251 :
252 51 : LG_FREE_ALL ;
253 51 : return (GrB_SUCCESS) ;
254 : #else
255 : return (GrB_NOT_IMPLEMENTED) ;
256 : #endif
257 : }
|
