Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LG_CC_Boruvka.c: connected components using GrB* methods only
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 : // Modified by Timothy A. Davis, Texas A&M University
16 :
17 : //------------------------------------------------------------------------------
18 :
19 : // This is an Advanced algorithm (G->is_symmetric_structure must be known),
20 : // but it is not user-callable (see LAGr_ConnectedComponents instead).
21 :
22 : // Code is based on Boruvka's minimum spanning forest algorithm.
23 :
24 : // This method relies solely on GrB* methods in the V2.0 C API, but it much
25 : // slower in general than LG_CC_FastSV6, which uses GxB pack/unpack methods
26 : // for faster access to the contents of the matrices and vectors.
27 :
28 : #include <stdint.h>
29 :
30 : #include "LG_internal.h"
31 :
32 : //------------------------------------------------------------------------------
33 : // Reduce_assign
34 : //------------------------------------------------------------------------------
35 :
36 : // w[Px[i]] = min(w[Px[i]], s[i]) for i in [0..n-1].
37 :
38 40 : static GrB_Info Reduce_assign
39 : (
40 : GrB_Vector w, // input/output vector of size n
41 : GrB_Vector s, // input vector of size n
42 : GrB_Index *Px, // Px: array of size n
43 : GrB_Index *mem, // workspace of size 3*n
44 : GrB_Index n
45 : )
46 : {
47 40 : char *msg = NULL ;
48 40 : GrB_Index *ind = mem ;
49 40 : GrB_Index *sval = ind + n ;
50 40 : GrB_Index *wval = sval + n ;
51 40 : GRB_TRY (GrB_Vector_extractTuples (ind, wval, &n, w)) ;
52 40 : GRB_TRY (GrB_Vector_extractTuples (ind, sval, &n, s)) ;
53 34190 : for (GrB_Index j = 0 ; j < n ; j++)
54 : {
55 34150 : if (sval [j] < wval [Px [j]])
56 : {
57 12526 : wval [Px [j]] = sval [j] ;
58 : }
59 : }
60 40 : GRB_TRY (GrB_Vector_clear (w)) ;
61 40 : GRB_TRY (GrB_Vector_build (w, ind, wval, n, GrB_PLUS_UINT64)) ;
62 : LG_FREE_ALL ;
63 40 : return (GrB_SUCCESS) ;
64 : }
65 :
66 : //------------------------------------------------------------------------------
67 : // select_func: IndexUnaryOp for pruning entries from S
68 : //------------------------------------------------------------------------------
69 :
70 : // The pointer to the Px array is passed to the select function as a component
71 : // of a user-defined data type.
72 :
73 : typedef struct
74 : {
75 : GrB_Index *pointer ;
76 : }
77 : Parent_struct ;
78 :
79 507956 : void my_select_func (void *z, const void *x,
80 : const GrB_Index i, const GrB_Index j, const void *y)
81 : {
82 507956 : Parent_struct *Parent = (Parent_struct *) y ;
83 507956 : GrB_Index *Px = Parent->pointer ;
84 507956 : (*((bool *) z)) = (Px [i] != Px [j]) ;
85 507956 : }
86 :
87 : //------------------------------------------------------------------------------
88 : // LG_CC_Boruvka
89 : //------------------------------------------------------------------------------
90 :
91 : #undef LG_FREE_ALL
92 : #define LG_FREE_ALL \
93 : { \
94 : LG_FREE_WORK ; \
95 : GrB_free (&parent) ; \
96 : }
97 :
98 : #undef LG_FREE_WORK
99 : #define LG_FREE_WORK \
100 : { \
101 : LAGraph_Free ((void **) &I, NULL) ; \
102 : LAGraph_Free ((void **) &Px, NULL) ; \
103 : LAGraph_Free ((void **) &mem, NULL) ; \
104 : GrB_free (&S) ; \
105 : GrB_free (&Parent_Type) ; \
106 : GrB_free (&gp) ; \
107 : GrB_free (&mnp) ; \
108 : GrB_free (&ccmn) ; \
109 : GrB_free (&ramp) ; \
110 : GrB_free (&mask) ; \
111 : GrB_free (&select_op) ; \
112 : }
113 :
114 50 : int LG_CC_Boruvka
115 : (
116 : // output:
117 : GrB_Vector *component, // output: array of component identifiers
118 : // input:
119 : const LAGraph_Graph G, // input graph, not modified
120 : char *msg
121 : )
122 : {
123 :
124 : //--------------------------------------------------------------------------
125 : // check inputs
126 : //--------------------------------------------------------------------------
127 :
128 50 : GrB_Index n, *I = NULL, *Px = NULL, *mem = NULL ;
129 50 : GrB_Vector parent = NULL, gp = NULL, mnp = NULL, ccmn = NULL, ramp = NULL,
130 50 : mask = NULL ;
131 50 : GrB_IndexUnaryOp select_op = NULL ;
132 50 : GrB_Matrix S = NULL ;
133 50 : GrB_Type Parent_Type = NULL ;
134 : Parent_struct Parent ;
135 :
136 50 : LG_CLEAR_MSG ;
137 50 : LG_TRY (LAGraph_CheckGraph (G, msg)) ;
138 49 : LG_ASSERT (component != NULL, GrB_NULL_POINTER) ;
139 :
140 25 : LG_ASSERT_MSG ((G->kind == LAGraph_ADJACENCY_UNDIRECTED ||
141 : (G->kind == LAGraph_ADJACENCY_DIRECTED &&
142 : G->is_symmetric_structure == LAGraph_TRUE)),
143 : LAGRAPH_SYMMETRIC_STRUCTURE_REQUIRED,
144 : "G->A must be known to be symmetric") ;
145 :
146 : //--------------------------------------------------------------------------
147 : // initializations
148 : //--------------------------------------------------------------------------
149 :
150 : // S = structure of G->A
151 24 : LG_TRY (LAGraph_Matrix_Structure (&S, G->A, msg)) ;
152 :
153 24 : GRB_TRY (GrB_Matrix_nrows (&n, S)) ;
154 24 : GRB_TRY (GrB_Vector_new (&parent, GrB_UINT64, n)) ; // final result
155 24 : GRB_TRY (GrB_Vector_new (&gp, GrB_UINT64, n)) ; // grandparents
156 24 : GRB_TRY (GrB_Vector_new (&mnp, GrB_UINT64, n)) ; // min neighbor parent
157 24 : GRB_TRY (GrB_Vector_new (&ccmn, GrB_UINT64, n)) ; // cc's min neighbor
158 24 : GRB_TRY (GrB_Vector_new (&mask, GrB_BOOL, n)) ; // various uses
159 :
160 24 : LG_TRY (LAGraph_Malloc ((void **) &mem, 3*n, sizeof (GrB_Index), msg)) ;
161 24 : LG_TRY (LAGraph_Malloc ((void **) &Px, n, sizeof (GrB_Index), msg)) ;
162 24 : Parent.pointer = Px ;
163 :
164 24 : GRB_TRY (GrB_Type_new (&Parent_Type, sizeof (Parent_struct))) ;
165 :
166 : #if !LAGRAPH_SUITESPARSE
167 : // I is not needed for SuiteSparse and remains NULL
168 : LG_TRY (LAGraph_Malloc ((void **) &I, n, sizeof (GrB_Index), msg)) ;
169 : #endif
170 :
171 : // parent = 0:n-1, and copy to ramp
172 24 : GRB_TRY (GrB_assign (parent, NULL, NULL, 0, GrB_ALL, n, NULL)) ;
173 24 : GRB_TRY (GrB_apply (parent, NULL, NULL, GrB_ROWINDEX_INT64, parent, 0,
174 : NULL)) ;
175 24 : GRB_TRY (GrB_Vector_dup (&ramp, parent)) ;
176 :
177 : // Px is a non-opaque copy of the parent GrB_Vector
178 24 : GRB_TRY (GrB_Vector_extractTuples (I, Px, &n, parent)) ;
179 :
180 24 : GRB_TRY (GrB_IndexUnaryOp_new (&select_op, my_select_func, GrB_BOOL,
181 : /* aij: ignored */ GrB_BOOL, /* y: pointer to Px */ Parent_Type)) ;
182 :
183 : GrB_Index nvals ;
184 24 : GRB_TRY (GrB_Matrix_nvals (&nvals, S)) ;
185 :
186 : //--------------------------------------------------------------------------
187 : // find the connected components
188 : //--------------------------------------------------------------------------
189 :
190 64 : while (nvals > 0)
191 : {
192 :
193 : //----------------------------------------------------------------------
194 : // mnp[u] = u's minimum neighbor's parent for all nodes u
195 : //----------------------------------------------------------------------
196 :
197 : // every vertex points to a root vertex at the begining
198 40 : GRB_TRY (GrB_assign (mnp, NULL, NULL, n, GrB_ALL, n, NULL)) ;
199 40 : GRB_TRY (GrB_mxv (mnp, NULL, GrB_MIN_UINT64,
200 : GrB_MIN_SECOND_SEMIRING_UINT64, S, parent, NULL)) ;
201 :
202 : //----------------------------------------------------------------------
203 : // find the minimum neighbor
204 : //----------------------------------------------------------------------
205 :
206 : // ccmn[u] = connect component's minimum neighbor | if u is a root
207 : // = n | otherwise
208 40 : GRB_TRY (GrB_assign (ccmn, NULL, NULL, n, GrB_ALL, n, NULL)) ;
209 40 : GRB_TRY (Reduce_assign (ccmn, mnp, Px, mem, n)) ;
210 :
211 : //----------------------------------------------------------------------
212 : // parent[u] = ccmn[u] if ccmn[u] != n
213 : //----------------------------------------------------------------------
214 :
215 : // mask = (ccnm != n)
216 40 : GRB_TRY (GrB_apply (mask, NULL, NULL, GrB_NE_UINT64, ccmn, n, NULL)) ;
217 : // parent<mask> = ccmn
218 40 : GRB_TRY (GrB_assign (parent, mask, NULL, ccmn, GrB_ALL, n, NULL)) ;
219 :
220 : //----------------------------------------------------------------------
221 : // select new roots
222 : //----------------------------------------------------------------------
223 :
224 : // identify all pairs (u,v) where parent [u] == v and parent [v] == u
225 : // and then select the minimum of u, v as the new root;
226 : // if (parent [parent [i]] == i) parent [i] = min (parent [i], i)
227 :
228 : // compute grandparents: gp = parent (parent)
229 40 : GRB_TRY (GrB_Vector_extractTuples (I, Px, &n, parent)) ;
230 40 : GRB_TRY (GrB_extract (gp, NULL, NULL, parent, Px, n, NULL)) ;
231 :
232 : // mask = (gp == 0:n-1)
233 40 : GRB_TRY (GrB_eWiseMult (mask, NULL, NULL, GrB_EQ_UINT64, gp, ramp,
234 : NULL)) ;
235 : // parent<mask> = min (parent, ramp)
236 40 : GRB_TRY (GrB_assign (parent, mask, GrB_MIN_UINT64, ramp, GrB_ALL, n,
237 : NULL)) ;
238 :
239 : //----------------------------------------------------------------------
240 : // shortcutting: parent [i] = parent [parent [i]] until convergence
241 : //----------------------------------------------------------------------
242 :
243 40 : bool changing = true ;
244 : while (true)
245 : {
246 : // compute grandparents: gp = parent (parent)
247 122 : GRB_TRY (GrB_Vector_extractTuples (I, Px, &n, parent)) ;
248 122 : GRB_TRY (GrB_extract (gp, NULL, NULL, parent, Px, n, NULL)) ;
249 :
250 : // changing = or (parent != gp)
251 122 : GRB_TRY (GrB_eWiseMult (mask, NULL, NULL, GrB_NE_UINT64, parent, gp,
252 : NULL)) ;
253 122 : GRB_TRY (GrB_reduce (&changing, NULL, GrB_LOR_MONOID_BOOL, mask,
254 : NULL)) ;
255 122 : if (!changing) break ;
256 :
257 : // parent = gp
258 82 : GRB_TRY (GrB_assign (parent, NULL, NULL, gp, GrB_ALL, n, NULL)) ;
259 : }
260 :
261 : //----------------------------------------------------------------------
262 : // remove the edges inside each connected component
263 : //----------------------------------------------------------------------
264 :
265 : // This step is the costliest part of this algorithm when dealing with
266 : // large matrices.
267 40 : GRB_TRY (GrB_Matrix_select_UDT (S, NULL, NULL, select_op, S,
268 : (void *) (&Parent), NULL)) ;
269 40 : GRB_TRY (GrB_Matrix_nvals (&nvals, S)) ;
270 : }
271 :
272 : //--------------------------------------------------------------------------
273 : // free workspace and return result
274 : //--------------------------------------------------------------------------
275 :
276 24 : (*component) = parent ;
277 24 : LG_FREE_WORK ;
278 24 : return (GrB_SUCCESS) ;
279 : }
280 :
|
