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