Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_msf.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 Boruvka's minimum spanning forest algorithm
20 : */
21 :
22 : // TODO: is this ready for src? It uses global values, so not yet ready.
23 : // TODO: Reduce_assign is slow. See src/algorithm/LG_CC_FastSV6/7.
24 :
25 : #include "LG_internal.h"
26 : #include <LAGraph.h>
27 : #include <LAGraphX.h>
28 :
29 : //****************************************************************************
30 : // encode each edge into a single uint64_t
31 103852 : static void combine (void *z, const void *x, const void *y)
32 : {
33 103852 : *(uint64_t*)z = ((*(uint64_t*)x) << 32) + (*(uint64_t*)y);
34 103852 : }
35 :
36 9771 : static void get_fst (void *y, const void *x)
37 : {
38 9771 : *(uint64_t*)y = (*(uint64_t*)x) >> 32;
39 9771 : }
40 :
41 16258 : static void get_snd (void *y, const void *x)
42 : {
43 16258 : *(uint64_t*)y = (*(uint64_t*)x) & INT_MAX;
44 16258 : }
45 :
46 : //****************************************************************************
47 : // TODO: Reduce_assign is slow. See src/algorithm/LG_CC_FastSV6.
48 :
49 : #undef LG_FREE_ALL
50 : #define LG_FREE_ALL LAGraph_Free ((void **) &mem, msg) ;
51 :
52 : // w[index[i]] = min(w[index[i]], s[i]) for i in [0..n-1]
53 30 : static GrB_Info Reduce_assign (GrB_Vector w,
54 : GrB_Vector s, GrB_Index *index, GrB_Index n, char *msg)
55 : {
56 30 : GrB_Index *mem = NULL ;
57 30 : LG_TRY (LAGraph_Malloc ((void **) &mem, n*3, sizeof (GrB_Index), msg)) ;
58 30 : GrB_Index *ind = mem, *sval = mem + n, *wval = sval + n;
59 30 : LG_TRY (GrB_Vector_extractTuples(ind, wval, &n, w));
60 30 : LG_TRY (GrB_Vector_extractTuples(ind, sval, &n, s));
61 13004 : for (GrB_Index i = 0; i < n; i++)
62 12974 : if (sval[i] < wval[index[i]])
63 6603 : wval[index[i]] = sval[i];
64 30 : LG_TRY (GrB_Vector_clear(w));
65 30 : LG_TRY (GrB_Vector_build(w, ind, wval, n, GrB_PLUS_UINT64));
66 30 : LG_FREE_ALL ;
67 30 : return GrB_SUCCESS;
68 : }
69 :
70 : //****************************************************************************
71 : // global C arrays (for implementing various GrB_IndexUnaryOps)
72 : static GrB_Index *weight = NULL, *parent = NULL, *partner = NULL;
73 :
74 : // generate solution:
75 : // for each element A(i, j), it is selected if
76 : // 1. weight[i] == A(i, j) -- where weight[i] stores i's minimum edge weight
77 : // 2. parent[j] == partner[i] -- j belongs to the specified connected component
78 :
79 192530 : void f1 (bool *z, const void *x, GrB_Index i, GrB_Index j, const void *thunk)
80 : {
81 192530 : uint64_t *aij = (uint64_t*) x;
82 192530 : (*z) = (weight[i] == *aij) && (parent[j] == partner[i]);
83 192530 : }
84 :
85 : // edge removal:
86 : // A(i, j) is removed when parent[i] == parent[j]
87 :
88 192530 : void f2 (bool *z, const void *x, GrB_Index i, GrB_Index j, const void *thunk)
89 : {
90 192530 : (*z) = (parent[i] != parent[j]);
91 192530 : }
92 :
93 : //****************************************************************************
94 :
95 : #undef LG_FREE_ALL
96 : #define LG_FREE_ALL \
97 : { \
98 : GrB_free (&S); \
99 : GrB_free (&T); \
100 : LAGraph_Free ((void **) &I, msg); \
101 : LAGraph_Free ((void **) &V, msg); \
102 : LAGraph_Free ((void **) &SI, msg); \
103 : LAGraph_Free ((void **) &SJ, msg); \
104 : LAGraph_Free ((void **) &SX, msg); \
105 : LAGraph_Free ((void **) &parent, msg); \
106 : LAGraph_Free ((void **) &partner, msg); \
107 : LAGraph_Free ((void **) &weight, msg); \
108 : GrB_free (&f); \
109 : GrB_free (&i); \
110 : GrB_free (&t); \
111 : GrB_free (&edge); \
112 : GrB_free (&cedge); \
113 : GrB_free (&mask); \
114 : GrB_free (&index); \
115 : GrB_free (&comb); \
116 : GrB_free (&combMin); \
117 : GrB_free (&fst); \
118 : GrB_free (&snd); \
119 : GrB_free (&s1); \
120 : GrB_free (&s2); \
121 : }
122 :
123 : //****************************************************************************
124 12 : int LAGraph_msf
125 : (
126 : GrB_Matrix *result, // output: an unsymmetrical matrix, the spanning forest
127 : GrB_Matrix A, // input matrix
128 : bool sanitize, // if true, ensure A is symmetric
129 : char *msg
130 : )
131 : {
132 : #if LAGRAPH_SUITESPARSE
133 :
134 12 : LG_CLEAR_MSG ;
135 :
136 : GrB_Info info;
137 : GrB_Index n;
138 12 : GrB_Matrix S = NULL, T = NULL;
139 12 : GrB_Vector f = NULL, i = NULL, t = NULL,
140 12 : edge = NULL, cedge = NULL, mask = NULL, index = NULL;
141 12 : GrB_Index *I = NULL, *V = NULL, *SI = NULL, *SJ = NULL, *SX = NULL;
142 :
143 12 : GrB_BinaryOp comb = NULL;
144 12 : GrB_Semiring combMin = NULL;
145 12 : GrB_UnaryOp fst = NULL, snd = NULL;
146 :
147 12 : GrB_IndexUnaryOp s1 = NULL, s2 = NULL;
148 12 : if (result == NULL || A == NULL) return (GrB_NULL_POINTER) ;
149 :
150 : GrB_Index ncols ;
151 11 : GRB_TRY (GrB_Matrix_nrows (&n, A));
152 11 : GRB_TRY (GrB_Matrix_ncols (&ncols, A));
153 11 : if (n != ncols) return (GrB_DIMENSION_MISMATCH) ;
154 :
155 10 : if (sanitize)
156 : {
157 : // S = A+A'
158 1 : GRB_TRY (GrB_Matrix_new (&S, GrB_UINT64, n, n));
159 1 : GRB_TRY (GrB_eWiseAdd (S, 0, 0, GrB_PLUS_UINT64, A, A, GrB_DESC_T1));
160 : }
161 : else
162 : {
163 : // Use the input as-is, and assume it is GrB_UINT64 and symmetric
164 9 : GRB_TRY (GrB_Matrix_dup (&S, A));
165 : }
166 :
167 10 : GRB_TRY (GrB_Matrix_new (&T, GrB_UINT64, n, n));
168 :
169 10 : GRB_TRY (GrB_Vector_new (&t, GrB_UINT64, n));
170 10 : GRB_TRY (GrB_Vector_new (&f, GrB_UINT64, n));
171 10 : GRB_TRY (GrB_Vector_new (&i, GrB_UINT64, n));
172 10 : GRB_TRY (GrB_Vector_new (&edge, GrB_UINT64, n));
173 10 : GRB_TRY (GrB_Vector_new (&cedge, GrB_UINT64, n));
174 10 : GRB_TRY (GrB_Vector_new (&mask, GrB_BOOL, n));
175 10 : GRB_TRY (GrB_Vector_new (&index, GrB_UINT64, n));
176 :
177 : // temporary arrays
178 10 : LG_TRY (LAGraph_Malloc ((void **) &I, n, sizeof (GrB_Index), msg)) ;
179 10 : LG_TRY (LAGraph_Malloc ((void **) &V, n, sizeof (GrB_Index), msg)) ;
180 10 : LG_TRY (LAGraph_Malloc ((void **) &SI, 2*n, sizeof (GrB_Index), msg)) ;
181 10 : LG_TRY (LAGraph_Malloc ((void **) &SJ, 2*n, sizeof (GrB_Index), msg)) ;
182 10 : LG_TRY (LAGraph_Malloc ((void **) &SX, 2*n, sizeof (GrB_Index), msg)) ;
183 :
184 : // global arrays
185 10 : LG_TRY (LAGraph_Malloc ((void **) &parent, n, sizeof (GrB_Index), msg)) ;
186 10 : LG_TRY (LAGraph_Malloc ((void **) &weight, n, sizeof (GrB_Index), msg)) ;
187 10 : LG_TRY (LAGraph_Malloc ((void **) &partner, n, sizeof (GrB_Index), msg)) ;
188 :
189 : // prepare vectors
190 3304 : for (GrB_Index i = 0; i < n; i++)
191 3294 : I[i] = parent[i] = i;
192 10 : GRB_TRY (GrB_Vector_build (f, I, parent, n, GrB_PLUS_UINT64));
193 10 : GRB_TRY (GrB_assign (i, 0, 0, f, GrB_ALL, 0, 0));
194 :
195 : // semiring & monoid
196 10 : GrB_Index inf = ((uint64_t) INT_MAX << 32) ^ INT_MAX;
197 10 : GRB_TRY (GrB_BinaryOp_new (&comb, combine, GrB_UINT64, GrB_UINT64, GrB_UINT64));
198 10 : GRB_TRY (GrB_Semiring_new (&combMin, GrB_MIN_MONOID_UINT64, comb));
199 10 : GRB_TRY (GrB_UnaryOp_new (&fst, get_fst, GrB_UINT64, GrB_UINT64));
200 10 : GRB_TRY (GrB_UnaryOp_new (&snd, get_snd, GrB_UINT64, GrB_UINT64));
201 :
202 : // ops for GrB_select
203 10 : GrB_IndexUnaryOp_new (&s1, (void *) f1, GrB_BOOL, GrB_UINT64, GrB_UINT64);
204 10 : GrB_IndexUnaryOp_new (&s2, (void *) f2, GrB_BOOL, GrB_UINT64, GrB_UINT64);
205 :
206 : // the main computation
207 10 : GrB_Index nvals, diff, ntuples = 0, num;
208 10 : GRB_TRY (GrB_Matrix_nvals (&nvals, S));
209 15 : for (int iters = 1; nvals > 0; iters++)
210 : {
211 : // every vertex points to a root vertex at the beginning
212 : // edge[u] = u's minimum edge (weight and index are encoded together)
213 15 : GRB_TRY (GrB_assign (edge, 0, 0, inf, GrB_ALL, 0, 0));
214 15 : GRB_TRY (GrB_mxv (edge, 0, GrB_MIN_UINT64, combMin, S, f, 0));
215 : // cedge[u] = children's minimum edge | if u is a root
216 : // = (INT_MAX, u) | otherwise
217 15 : GRB_TRY (GrB_assign (t, 0, 0, (uint64_t) INT_MAX, GrB_ALL, 0, 0));
218 15 : GRB_TRY (GrB_eWiseMult (cedge, 0, 0, comb, t, i, 0));
219 15 : LG_TRY (Reduce_assign (cedge, edge, parent, n, msg));
220 : // if (f[u] == u) f[u] := snd(cedge[u]) -- the index part of the edge
221 15 : GRB_TRY (GrB_eWiseMult (mask, 0, 0, GrB_EQ_UINT64, f, i, 0));
222 15 : GRB_TRY (GrB_apply (f, mask, GrB_SECOND_UINT64, snd, cedge, 0));
223 : // identify all the vertex pairs (u, v) where f[u] == v and f[v] == u
224 : // and then select the minimum of u, v as the new root;
225 : // if (f[f[i]] == i) f[i] = min(f[i], i)
226 15 : GRB_TRY (GrB_Vector_extractTuples (I, V, &n, f));
227 15 : GRB_TRY (GrB_extract (t, 0, 0, f, V, n, 0));
228 15 : GRB_TRY (GrB_eWiseMult (mask, 0, 0, GrB_EQ_UINT64, i, t, 0));
229 15 : GRB_TRY (GrB_assign (f, mask, GrB_MIN_UINT64, i, GrB_ALL, 0, 0));
230 :
231 : // five steps to generate the solution
232 : // 1. new roots (f[i] == i) revise their entries in cedge
233 15 : GRB_TRY (GrB_eWiseMult (mask, 0, 0, GrB_EQ_UINT64, i, f, 0));
234 15 : GRB_TRY (GrB_assign (cedge, mask, 0, inf, GrB_ALL, 0, 0));
235 :
236 : // 2. every vertex tries to know whether one of its edges is selected
237 15 : GRB_TRY (GrB_extract (t, 0, 0, cedge, parent, n, 0));
238 15 : GRB_TRY (GrB_eWiseMult (mask ,0, 0, GrB_EQ_UINT64, edge, t, 0));
239 :
240 : // 3. each root picks a vertex from its children to generate the solution
241 15 : GRB_TRY (GrB_assign (index, 0, 0, n, GrB_ALL, 0, 0));
242 15 : GRB_TRY (GrB_assign (index, mask, 0, i, GrB_ALL, 0, 0));
243 15 : GRB_TRY (GrB_assign (t, 0, 0, n, GrB_ALL, 0, 0));
244 15 : LG_TRY (Reduce_assign (t, index, parent, n, msg));
245 15 : GRB_TRY (GrB_extract (index, 0, 0, t, parent, n, 0));
246 15 : GRB_TRY (GrB_eWiseMult (mask ,0, 0, GrB_EQ_UINT64, i, index, 0));
247 :
248 : // 4. generate the select function (set the global pointers)
249 15 : GRB_TRY (GrB_assign (t, 0, 0, inf, GrB_ALL, 0, 0));
250 15 : GRB_TRY (GrB_apply (t, mask, 0, fst, edge, 0));
251 15 : GRB_TRY (GrB_Vector_extractTuples (I, weight, &n, t));
252 15 : GRB_TRY (GrB_assign (t, 0, 0, inf, GrB_ALL, 0, 0));
253 15 : GRB_TRY (GrB_apply (t, mask, 0, snd, edge, 0));
254 15 : GRB_TRY (GrB_Vector_extractTuples (I, partner, &n, t));
255 15 : GRB_TRY (GrB_select (T, 0, 0, s1, S, 0, 0));
256 15 : GRB_TRY (GrB_Vector_clear (t));
257 :
258 : // 5. the generated matrix may still have redundant edges
259 : // remove the duplicates by GrB_mxv() and store them as tuples
260 15 : GRB_TRY (GrB_Vector_clear (edge));
261 15 : GRB_TRY (GrB_mxv (edge, mask, GrB_MIN_UINT64, combMin, T, i, 0));
262 15 : GRB_TRY (GrB_Vector_nvals (&num, edge));
263 15 : GRB_TRY (GrB_apply (t, 0, 0, snd, edge, 0));
264 15 : GRB_TRY (GrB_Vector_extractTuples (SI + ntuples, SJ + ntuples, &num, t));
265 15 : GRB_TRY (GrB_apply (t, 0, 0, fst, edge, 0));
266 15 : GRB_TRY (GrB_Vector_extractTuples (SI + ntuples, SX + ntuples, &num, t));
267 15 : GRB_TRY (GrB_Vector_clear (t));
268 15 : ntuples += num;
269 :
270 : // path halving until every vertex points on a root
271 : do {
272 45 : GRB_TRY (GrB_Vector_extractTuples (I, V, &n, f));
273 45 : GRB_TRY (GrB_extract (t, 0, 0, f, V, n, 0));
274 45 : GRB_TRY (GrB_eWiseMult (mask, 0, 0, GrB_NE_UINT64, f, t, 0));
275 45 : GRB_TRY (GrB_assign (f, 0, 0, t, GrB_ALL, 0, 0));
276 45 : GRB_TRY (GrB_reduce (&diff, 0, GrB_PLUS_MONOID_UINT64, mask, 0));
277 45 : } while (diff != 0);
278 :
279 : // remove the edges in the same connected component
280 15 : GRB_TRY (GrB_Vector_extractTuples (I, parent, &n, f));
281 15 : GRB_TRY (GrB_select (S, 0, 0, s2, S, 0, 0));
282 15 : GrB_Matrix_nvals (&nvals, S);
283 15 : if (nvals == 0) break;
284 : }
285 10 : GRB_TRY (GrB_Matrix_clear (T));
286 10 : GRB_TRY (GrB_Matrix_build (T, SI, SJ, SX, ntuples, GrB_SECOND_UINT64));
287 10 : *result = T;
288 10 : T = NULL ;
289 :
290 10 : LG_FREE_ALL;
291 10 : return (GrB_SUCCESS) ;
292 : #else
293 : return (GrB_NOT_IMPLEMENTED) ;
294 : #endif
295 : }
|
