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