Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGr_MaximumMatching: maximum matching between nodes of disjoint sets
3 : // in bipartite graphs
4 : //------------------------------------------------------------------------------
5 :
6 : // LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved.
7 : // SPDX-License-Identifier: BSD-2-Clause
8 : //
9 : // For additional details (including references to third party source code and
10 : // other files) see the LICENSE file or contact permission@sei.cmu.edu. See
11 : // Contributors.txt for a full list of contributors. Created, in part, with
12 : // funding and support from the U.S. Government (see Acknowledgments.txt file).
13 : // DM22-0790
14 :
15 : // Contributed by Christina Koutsou, Aristotle University of Thessaloniki
16 :
17 : //------------------------------------------------------------------------------
18 :
19 : // This implmentation is based on the algorithm described in the following
20 : // paper: "Distributed-Memory Algorithms for Maximum Cardinality Matching in
21 : // Bipartite Graphs" by A. Azad and A. Buluç.
22 : //
23 : // In brief, this algorithm explores alternate paths by reversing an already
24 : // traversed path to see if one of the already matched columns encountered in
25 : // the path has at least one free children to be matched with instead. If so, it
26 : // concedes its previous match to another previously matched parent or to the
27 : // unmatched root of the path.
28 :
29 : // More detailed explanation of the algorithm and its implementation can be
30 : // found in the LAGraph "papers" folder, in file
31 : // "MaximumMatching_Report.pdf".
32 :
33 : // This is an Advanced method, G->A is required as input instead of G due to
34 : // lack of a Bipartite graph kind in GraphBLAS.
35 :
36 : //------------------------------------------------------------------------------
37 :
38 : #include "LAGraphX.h"
39 : #include "LG_internal.h"
40 :
41 : //------------------------------------------------------------------------------
42 : // the Vertex tuple: (parentC, rootC)
43 : //------------------------------------------------------------------------------
44 :
45 : typedef struct
46 : {
47 : uint64_t parentC;
48 : uint64_t rootC;
49 : } vertex;
50 :
51 : // repeat the typedef as a string, to give to GraphBLAS
52 : #define VERTEX_DEFN \
53 : "typedef struct \n" \
54 : "{ \n" \
55 : " uint64_t parentC; \n" \
56 : " uint64_t rootC; \n" \
57 : "} vertex;"
58 :
59 73362 : void initFrontier(vertex *z, const bool *x, uint64_t i, uint64_t j, const bool *y)
60 : {
61 73362 : z->parentC = i;
62 73362 : z->rootC = i;
63 73362 : }
64 :
65 : #define INIT_FRONTIER_DEFN \
66 : "void initFrontier(vertex *z, const bool *x, uint64_t i, uint64_t j, const bool *y) \n" \
67 : "{ \n" \
68 : " z->parentC = i; \n" \
69 : " z->rootC = i; \n" \
70 : "}"
71 :
72 32261 : void minparent(vertex *z, const vertex *x, const vertex *y)
73 : {
74 32261 : *z = x->parentC < y->parentC ? *x : *y;
75 32261 : }
76 :
77 : #define MIN_PARENT_DEFN \
78 : "void minparent(vertex *z, const vertex *x, const vertex *y) \n" \
79 : "{ \n" \
80 : " *z = x->parentC < y->parentC ? *x : *y; \n" \
81 : "}"
82 :
83 : // TODO: revise GraphBLAS so we can tell it that the select2nd operator
84 : // does not use the 'x' input.
85 33420 : void select2nd(vertex *z, const bool *x, const vertex *y)
86 : {
87 33420 : z->parentC = y->parentC;
88 33420 : z->rootC = y->rootC;
89 33420 : }
90 :
91 : #define SELECT_2ND_DEFN \
92 : "void select2nd(vertex *z, const bool *x, const vertex *y) \n" \
93 : "{ \n" \
94 : " z->parentC = y->parentC; \n" \
95 : " z->rootC = y->rootC; \n" \
96 : "}"
97 :
98 23430 : void select1st(vertex *z, const vertex *x, const bool *y)
99 : {
100 23430 : z->parentC = x->parentC;
101 23430 : z->rootC = x->rootC;
102 23430 : }
103 :
104 : #define SELECT_1ST_DEFN \
105 : "void select1st(vertex *z, const vertex *x, const bool *y) \n" \
106 : "{ \n" \
107 : " z->parentC = x->parentC; \n" \
108 : " z->rootC = x->rootC; \n" \
109 : "}"
110 :
111 20280 : void keepParents(uint64_t *z, const vertex *x) { *z = x->parentC; }
112 :
113 : #define KEEP_PARENTS_DEFN \
114 : "void keepParents(uint64_t *z, const vertex *x) { *z = x->parentC; }\n"
115 :
116 14199 : void keepRoots(uint64_t *z, const vertex *x) { *z = x->rootC; }
117 :
118 : #define KEEP_ROOTS_DEFN \
119 : "void keepRoots(uint64_t *z, const vertex *x) { *z = x->rootC; }\n"
120 :
121 2277 : void buildfCTuples(vertex *z, const uint64_t *x, uint64_t i, uint64_t j,
122 : const bool *y)
123 : {
124 2277 : z->parentC = i;
125 2277 : z->rootC = *x;
126 2277 : }
127 :
128 : #define BUILT_FC_TUPLES_DEFN \
129 : "void buildfCTuples(vertex *z, const uint64_t *x, uint64_t i, uint64_t j, \n" \
130 : " const bool *y) \n" \
131 : "{ \n" \
132 : " z->parentC = i; \n" \
133 : " z->rootC = *x; \n" \
134 : "}"
135 :
136 6081 : void vertexTypecast(vertex *z, const uint64_t *x)
137 : {
138 6081 : z->parentC = *x;
139 6081 : z->rootC = *x;
140 6081 : }
141 :
142 : #define VERTEX_TYPECAST_DEFN \
143 : "void vertexTypecast(vertex *z, const uint64_t *x) \n" \
144 : "{ \n" \
145 : " z->parentC = *x; \n" \
146 : " z->rootC = *x; \n" \
147 : "}"
148 :
149 6081 : void setParentsMates(vertex *z, const vertex *x, const vertex *y)
150 : {
151 6081 : z->parentC = y->parentC;
152 6081 : z->rootC = x->rootC;
153 6081 : }
154 :
155 : #define SET_PARENTS_MATES_DEFN \
156 : "void setParentsMates(vertex *z, const vertex *x, const vertex *y) \n" \
157 : "{ \n" \
158 : " z->parentC = y->parentC; \n" \
159 : " z->rootC = x->rootC; \n" \
160 : "}"
161 :
162 : //------------------------------------------------------------------------------
163 : // invert
164 : //------------------------------------------------------------------------------
165 :
166 : // This function "inverts" an input vector by swapping its row indices
167 : // and its values, returning the result in an output vector.
168 : //
169 : // For example, for the indices/values of an input vector (in) with 5 entries
170 : // and length 100:
171 : //
172 : // indices: 0 3 5 42 99
173 : // values: 4 98 1 3 12
174 : //
175 : // on output, the out vector will contain:
176 : //
177 : // indices: 4 98 1 3 12
178 : // values: 0 3 5 42 99
179 : //
180 : // The output vector will normally be jumbled since the values will not appear
181 : // in any particular order. The method assumes that the input values are in
182 : // range 0 to n-1 where n = length(out). The values in the input vector
183 : // may be duplicated and this argument of the function must be set accordingly.
184 : // Both the in vector and out vector must have the same type (GrB_UINT64). The
185 : // lengths of the two vectors need not be the same, so long as the indices
186 : // remain in range. Results are undefined if these conditions do not hold.
187 : //
188 : // The in and out vectors may be aliased. If not aliased, the input vector is
189 : // cleared of all entries on output. If in and out are aliased, then the
190 : // inversion is performed in-place.
191 : //
192 : // In SuiteSparse:GraphBLAS, this method takes O(1) time if the in vector is in
193 : // CSC (sparse, by column) format. Otherwise it can take O(e) time if e =
194 : // nvals(in), because the unpack below will convert the in vector to CSC and
195 : // then unpack it.
196 :
197 : #undef LG_FREE_ALL
198 : #define LG_FREE_ALL \
199 : { \
200 : LAGraph_Free((void *)&I, NULL); \
201 : LAGraph_Free((void *)&X1, NULL); \
202 : LAGraph_Free((void *)&X2, NULL); \
203 : }
204 :
205 6 : static inline GrB_Info invert_nondestructive(
206 : GrB_Vector out, // input/output. On input, only the size and type are
207 : // kept; any entries in the 'out' vector cleared. It is
208 : // then replaced with the inversion of the input vector.
209 : GrB_Vector in, // input vector, not modified. There must be no duplicate
210 : // values in the input vector.
211 : char *msg)
212 : {
213 6 : bool jumbled = 1;
214 6 : GrB_Index *I = NULL;
215 6 : GrB_Index *X1 = NULL;
216 6 : GrB_Index *X2 = NULL; // not used
217 6 : GrB_Index IBytes = 0, XBytes = 0;
218 6 : uint64_t nvals = 0;
219 :
220 : // the input and output vectors cannot be the same vector
221 : ASSERT(in != out);
222 :
223 : // All input/output vectors must be of type GrB_UINT64.
224 : #if LAGRAPH_SUITESPARSE
225 6 : GRB_TRY(
226 : GxB_Vector_unpack_CSC(in, (GrB_Index **)&I, (void **)&X1, &IBytes,
227 : &XBytes, NULL, &nvals, &jumbled,
228 : NULL)); // the output and input should have no
229 : // duplicates, so the order doesn't matter
230 : #else
231 : GRB_TRY(GrB_Vector_nvals(&nvals, in));
232 : LG_TRY(LAGraph_Malloc((void **)&I, nvals, sizeof(GrB_Index), msg));
233 : LG_TRY(LAGraph_Malloc((void **)&X1, nvals, sizeof(GrB_Index), msg));
234 : GRB_TRY(GrB_Vector_wait(in, GrB_MATERIALIZE));
235 : GRB_TRY(GrB_Vector_extractTuples_UINT64(
236 : I, X1, &nvals, in)); // the output and input should have no
237 : // duplicates, so the order doesn't matter
238 : #endif
239 :
240 6 : GRB_TRY(GrB_Vector_clear(out)); // clear the output first as a prerequisite
241 : // of the build method
242 6 : GRB_TRY(GrB_Vector_build_UINT64(
243 : out, X1, I, nvals,
244 : NULL)); // build does not take ownership of the lists I and X,
245 : // but only copies them, these lists will be given
246 : // again to the input
247 : // the input should have no duplicates in the
248 : // values list, so dups are not handled
249 : #if LAGRAPH_SUITESPARSE
250 6 : GRB_TRY(GxB_Vector_pack_CSC(in, (GrB_Index **)&I, (void **)&X1, IBytes,
251 : XBytes, false, nvals, jumbled, NULL));
252 : #endif
253 6 : return (GrB_SUCCESS);
254 : }
255 :
256 : static inline GrB_Info
257 1650 : invert(GrB_Vector out, // input/output. Same as invert_nondescructive above.
258 : GrB_Vector in, // input vector, empty on output (unless in == out)
259 : const bool dups, // flag that indicates if duplicates exist in the input
260 : // vector's values
261 : char *msg)
262 : {
263 : // The input and output vectors can be the same vector
264 : // that is, in == out is OK.
265 : // All input/output vectors must be of type GrB_UINT64.
266 :
267 : // the output vector will normally be returned in a jumbled state
268 1650 : bool jumbled = dups ? 0 : 1; // if there are duplicates, we want the indices
269 : // to be ordered so we can keep the min child
270 1650 : GrB_Index *I = NULL; // unpack allocates space for these lists
271 1650 : GrB_Index *X1 = NULL;
272 1650 : GrB_Index *X2 = NULL; // not used
273 1650 : GrB_Index IBytes = 0, XBytes = 0;
274 1650 : uint64_t nvals = 0;
275 :
276 : #if LAGRAPH_SUITESPARSE
277 1650 : GRB_TRY(GxB_Vector_unpack_CSC(in, (GrB_Index **)&I, (void **)&X1, &IBytes,
278 : &XBytes, NULL, &nvals, &jumbled, NULL));
279 : #else
280 : // vanilla case using extractTuples and build:
281 : // allocate I and X for GrB_extractTuples
282 : GRB_TRY(GrB_Vector_nvals(&nvals, in));
283 : LG_TRY(LAGraph_Malloc((void **)&I, nvals, sizeof(GrB_Index), msg));
284 : LG_TRY(LAGraph_Malloc((void **)&X1, nvals, sizeof(GrB_Index), msg));
285 : GRB_TRY(GrB_Vector_wait(in, GrB_MATERIALIZE));
286 : GRB_TRY(GrB_Vector_extractTuples_UINT64(
287 : I, X1, &nvals, in)); // the output and input should have no
288 : // duplicates, so the order doesn't matter
289 : GRB_TRY(GrB_Vector_clear(in));
290 : #endif
291 1650 : if (!dups)
292 : {
293 : #if LAGRAPH_SUITESPARSE
294 1386 : GRB_TRY(GxB_Vector_pack_CSC(out, (GrB_Index **)&X1, (void **)&I, XBytes,
295 : IBytes, false, nvals, true, NULL));
296 : #else
297 : GRB_TRY(GrB_Vector_clear(out));
298 : // GrB_MIN_UINT64 is used instead of first because first is an extension
299 : GRB_TRY(GrB_Vector_build_UINT64(out, X1, I, nvals, GrB_MIN_UINT64));
300 : // build copies the lists so they need to be freed in LG_FREE_ALL
301 : LG_FREE_ALL;
302 : #endif
303 : }
304 : else
305 : {
306 264 : GRB_TRY(GrB_Vector_clear(out));
307 : #if LAGRAPH_SUITESPARSE
308 264 : GRB_TRY(GrB_Vector_build_UINT64(out, X1, I, nvals, GrB_FIRST_UINT64));
309 : #else
310 : GRB_TRY(GrB_Vector_build_UINT64(out, X1, I, nvals, GrB_MIN_UINT64));
311 : #endif
312 : // build copies the lists so they need to be freed in LG_FREE_ALL
313 264 : LG_FREE_ALL;
314 : }
315 1650 : return (GrB_SUCCESS);
316 : }
317 :
318 : static inline GrB_Info
319 810 : invert_2(GrB_Vector out, // input/output
320 : GrB_Vector in1, // input vector, empty on output (unless in1 == out)
321 : GrB_Vector in2, // input vector, empty on output (unless in2 == out)
322 : const bool dups, // flag that indicates if duplicates exist in the
323 : // input vector's values
324 : const bool udt, // type of the first input, in1, and out vectors
325 : char *msg)
326 : {
327 : // The input vectors cannot be aliased. However in1==out or in2==out is
328 : // OK. The two input vectors must have the same # of entries.
329 : // All input/output vectors must be of type GrB_UINT64.
330 : ASSERT(in1 != in2);
331 :
332 810 : GrB_Index *I = NULL;
333 810 : GrB_Index *X1 = NULL;
334 810 : GrB_Index *X2 = NULL;
335 810 : GrB_Index IBytes = 0, X1Bytes = 0, X2Bytes = 0;
336 810 : uint64_t nvals1 = 0, nvals2 = 0;
337 :
338 : #if LAGRAPH_SUITESPARSE
339 810 : GRB_TRY(GxB_Vector_unpack_CSC(in1, (GrB_Index **)&I, (void **)&X1, &IBytes,
340 : &X1Bytes, NULL, &nvals1, NULL, NULL));
341 810 : LAGraph_Free((void *)&I, NULL);
342 810 : GRB_TRY(GxB_Vector_unpack_CSC(in2, (GrB_Index **)&I, (void **)&X2, &IBytes,
343 : &X2Bytes, NULL, &nvals2, NULL, NULL));
344 : ASSERT(nvals1 == nvals2);
345 810 : if (!dups)
346 : {
347 576 : LAGraph_Free((void *)&I, NULL);
348 576 : GRB_TRY(GxB_Vector_pack_CSC(out, (GrB_Index **)&X2, (void **)&X1,
349 : X2Bytes, X1Bytes, false, nvals2, true,
350 : NULL)); // the values are not ordered,
351 : // so the indices are jumbled
352 : }
353 : else
354 : {
355 234 : GRB_TRY(GrB_Vector_clear(out));
356 234 : GRB_TRY(GrB_Vector_build_UINT64(out, X2, X1, nvals2, GrB_FIRST_UINT64));
357 234 : LG_FREE_ALL;
358 : }
359 : #else
360 : // vanilla case using extractTuples and build:
361 : if (dups)
362 : {
363 : // invert in2 to eliminate dups and get the final indices of out (we
364 : // cannot invert into out if out is of type udt, because we would have a
365 : // domain mismatch with in2)
366 : GrB_Vector in_rev = NULL;
367 : GrB_Index out_size = 0;
368 : GRB_TRY(GrB_Vector_size(&out_size, out));
369 : GRB_TRY(GrB_Vector_new(&in_rev, GrB_UINT64, out_size));
370 : LG_TRY(invert(in_rev, in2, true, msg));
371 : GRB_TRY(GrB_Vector_nvals(
372 : &nvals2,
373 : in_rev)); // out may have less entries than in2 and, thus, in1
374 : LG_TRY(LAGraph_Malloc((void **)&I, nvals2, sizeof(GrB_Index), msg));
375 : LG_TRY(LAGraph_Malloc((void **)&X2, nvals2, sizeof(GrB_Index), msg));
376 : // get indices of out
377 : GRB_TRY(GrB_Vector_extractTuples_UINT64(
378 : X2, I, &nvals2,
379 : in_rev)); // X2 are the final values of in2 and will be the indices
380 : // of out. I are the indices of in2 (which are a subset of
381 : // the indices of in1)
382 : GRB_TRY(GrB_Vector_free(&in_rev));
383 : }
384 : else
385 : {
386 : GRB_TRY(GrB_Vector_nvals(&nvals2, in2));
387 : LG_TRY(LAGraph_Malloc((void **)&I, nvals2, sizeof(GrB_Index), msg));
388 : LG_TRY(LAGraph_Malloc((void **)&X2, nvals2, sizeof(GrB_Index), msg));
389 : GRB_TRY(GrB_Vector_wait(in2, GrB_MATERIALIZE));
390 : GRB_TRY(GrB_Vector_extractTuples_UINT64(
391 : I, X2, &nvals2, in2)); // X2 will be the indices of out and I are
392 : // the indices of in2 (which are the indices
393 : // of in1 as well)
394 : }
395 :
396 : GRB_TRY(GrB_Vector_clear(out));
397 : GRB_TRY(GrB_Vector_wait(in1, GrB_MATERIALIZE));
398 : if (udt)
399 : {
400 : vertex *X1_V = NULL;
401 : LG_TRY(LAGraph_Malloc((void **)&X1_V, nvals2, sizeof(vertex), msg));
402 : for (uint64_t i = 0; i < nvals2; i++)
403 : {
404 : GRB_TRY(GrB_Vector_extractElement_UDT(
405 : (void *)&X1_V[i], in1,
406 : (GrB_Index)I[i])); // value I[i] corresponds to pos X2[i] in
407 : // out due to the tuple extraction and X[i]
408 : // will replace it
409 : }
410 : GRB_TRY(GrB_Vector_clear(in1));
411 : GRB_TRY(GrB_Vector_build_UDT(out, X2, (void *)X1_V, nvals2,
412 : NULL)); // no dups are left
413 : LAGRAPH_TRY(LAGraph_Free((void **)&X1_V, msg));
414 : }
415 : else
416 : {
417 : LG_TRY(LAGraph_Malloc((void **)&X1, nvals2, sizeof(GrB_UINT64), msg));
418 : for (GrB_Index i = 0; i < nvals2; i++)
419 : {
420 : GRB_TRY(GrB_Vector_extractElement_UINT64(&X1[i], in1, I[i]));
421 : }
422 : GRB_TRY(GrB_Vector_clear(in1));
423 : GRB_TRY(GrB_Vector_build_UINT64(out, X2, X1, nvals2, NULL));
424 : }
425 : LG_FREE_ALL;
426 : #endif
427 810 : return (GrB_SUCCESS);
428 : }
429 :
430 : //------------------------------------------------------------------------------
431 : // LAGr_MaximumMatching
432 : //------------------------------------------------------------------------------
433 :
434 : #undef LG_FREE_WORK
435 : #define LG_FREE_WORK \
436 : { \
437 : GrB_free(&pathC); \
438 : GrB_free(&parentsR); \
439 : GrB_free(&Vertex); \
440 : GrB_free(&frontierC); \
441 : GrB_free(&frontierR); \
442 : GrB_free(&initFrontierOp); \
443 : GrB_free(&I); \
444 : GrB_free(&MinParent); \
445 : GrB_free(&MinParent_Monoid); \
446 : GrB_free(&Select2ndOp); \
447 : GrB_free(&Select1stOp); \
448 : GrB_free(&MinParent_2nd_Semiring); \
449 : GrB_free(&MinParent_1st_Semiring); \
450 : GrB_free(&getParentsOp); \
451 : GrB_free(&getRootsOp); \
452 : GrB_free(&parentsUpdate); \
453 : GrB_free(&ufrontierR); \
454 : GrB_free(&mateC); \
455 : GrB_free(&mateR); \
456 : GrB_free(&rootsufR); \
457 : GrB_free(&pathUpdate); \
458 : GrB_free(&rootufRIndexes); \
459 : GrB_free(&rootsfR); \
460 : GrB_free(&rootfRIndexes); \
461 : GrB_free(&buildfCTuplesOp); \
462 : GrB_free(&vertexTypecastOp); \
463 : GrB_free(&setParentsMatesOp); \
464 : GrB_free(&vr); \
465 : GrB_free(&pathCopy); \
466 : GrB_free(¤tMatesR); \
467 : }
468 :
469 : #undef LG_FREE_ALL
470 : #define LG_FREE_ALL \
471 : { \
472 : LG_FREE_WORK; \
473 : GrB_free(&mateC); \
474 : }
475 :
476 30 : int LAGr_MaximumMatching(
477 : // outputs:
478 : GrB_Vector
479 : *mateC_handle, // mateC(j) = i : Column j of the C subset is matched
480 : // to row i of the R subset (ignored on input)
481 : GrB_Vector *mateR_handle, // mateR(i) = j : Row i of the R subset is matched
482 : // to column j of the C subset (ignored on input)
483 : // inputs:
484 : GrB_Matrix A, // input adjacency matrix, TODO: this should be a LAGraph
485 : // of a BIPARTITE kind
486 : GrB_Matrix AT, // transpose of the input adjacency matrix, necessary to
487 : // perform push-pull optimization
488 : // if NULL, the push-pull optimization is not performed
489 : GrB_Vector mate_init, // input only, not modified, ignored if NULL
490 : bool col_init, // flag to indicate if the initial matching is provided from
491 : // the columns' or from the rows' perspective, ignored if
492 : // mate_init is NULL
493 : char *msg)
494 : {
495 :
496 : //--------------------------------------------------------------------------
497 : // check inputs
498 : //--------------------------------------------------------------------------
499 :
500 30 : GrB_Vector pathC = NULL; // make this bitmap/sparse,
501 : // if dense I would have to give
502 : // all the entries and make the matrix 1-based
503 30 : GrB_Vector parentsR = NULL; // parents of row nodes that are reachable
504 : // from paths of the initial column frontier
505 30 : GrB_Type Vertex = NULL;
506 30 : GrB_Vector frontierC = NULL;
507 30 : GrB_Vector frontierR = NULL;
508 30 : GrB_IndexUnaryOp initFrontierOp = NULL;
509 30 : GrB_Vector I = NULL; // dense vector of 1's
510 30 : GrB_BinaryOp MinParent = NULL;
511 30 : GrB_Monoid MinParent_Monoid = NULL;
512 30 : GrB_BinaryOp Select2ndOp = NULL;
513 30 : GrB_BinaryOp Select1stOp = NULL;
514 30 : GrB_Semiring MinParent_2nd_Semiring = NULL;
515 30 : GrB_Semiring MinParent_1st_Semiring = NULL;
516 30 : GrB_UnaryOp getParentsOp = NULL;
517 30 : GrB_UnaryOp getRootsOp = NULL;
518 30 : GrB_Vector parentsUpdate = NULL; // unmatched rows of R frontier
519 30 : GrB_Vector ufrontierR = NULL; // unmatched rows of R frontier
520 30 : GrB_Vector mateR = NULL; // mateR(i) = j : Row i of the R subset is
521 : // matched to column j of the C subset
522 30 : GrB_Vector rootsufR = NULL;
523 30 : GrB_Vector pathUpdate = NULL;
524 30 : GrB_Vector rootufRIndexes = NULL;
525 30 : GrB_Vector rootsfR = NULL;
526 30 : GrB_Vector rootfRIndexes = NULL;
527 30 : GrB_IndexUnaryOp buildfCTuplesOp = NULL;
528 30 : GrB_UnaryOp vertexTypecastOp = NULL;
529 30 : GrB_BinaryOp setParentsMatesOp = NULL;
530 30 : GrB_Vector vr = NULL;
531 30 : GrB_Vector pathCopy = NULL;
532 30 : GrB_Vector currentMatesR = NULL;
533 :
534 30 : GrB_Vector mateC = NULL;
535 :
536 30 : LG_CLEAR_MSG;
537 :
538 30 : LG_ASSERT_MSG(mateC_handle != NULL || mateR_handle != NULL,
539 : GrB_NULL_POINTER, "At least one output must be provided\n");
540 30 : LG_ASSERT_MSG(A != NULL || AT != NULL, GrB_NULL_POINTER,
541 : "A matrix is NULL\n");
542 :
543 30 : if (mateC_handle != NULL)
544 : {
545 30 : (*mateC_handle) = NULL;
546 : }
547 30 : if (mateR_handle != NULL)
548 : {
549 20 : (*mateR_handle) = NULL;
550 : }
551 :
552 30 : bool do_pushpull = (AT != NULL) && (A != NULL);
553 :
554 30 : uint64_t ncols = 0;
555 30 : uint64_t nrows = 0;
556 :
557 30 : if (A != NULL)
558 : {
559 20 : GRB_TRY(GrB_Matrix_ncols(&ncols, A));
560 20 : GRB_TRY(GrB_Matrix_nrows(&nrows, A));
561 : }
562 : else
563 : {
564 10 : GRB_TRY(GrB_Matrix_nrows(&ncols, AT));
565 10 : GRB_TRY(GrB_Matrix_ncols(&nrows, AT));
566 : }
567 :
568 30 : GRB_TRY(GrB_Vector_new(&pathC, GrB_UINT64, ncols));
569 :
570 30 : GRB_TRY(GrB_Vector_new(&parentsR, GrB_UINT64, nrows));
571 :
572 : #if LAGRAPH_SUITESPARSE
573 30 : GRB_TRY(GxB_Type_new(&Vertex, sizeof(vertex), "vertex", VERTEX_DEFN));
574 : #else
575 : GRB_TRY(GrB_Type_new(&Vertex, sizeof(vertex)));
576 : #endif
577 :
578 30 : GRB_TRY(GrB_Vector_new(&frontierC, Vertex, ncols));
579 :
580 30 : GRB_TRY(GrB_Vector_new(&frontierR, Vertex, nrows));
581 :
582 : #if LAGRAPH_SUITESPARSE
583 30 : GRB_TRY(GxB_IndexUnaryOp_new(&initFrontierOp, (void *)initFrontier, Vertex,
584 : GrB_BOOL, GrB_BOOL, "initFrontier",
585 : INIT_FRONTIER_DEFN));
586 : #else
587 : GRB_TRY(GrB_IndexUnaryOp_new(&initFrontierOp, (void *)initFrontier, Vertex,
588 : GrB_BOOL, GrB_BOOL));
589 : #endif
590 :
591 30 : GRB_TRY(GrB_Vector_new(&I, GrB_BOOL, ncols));
592 30 : GRB_TRY(GrB_Vector_assign_BOOL(I, NULL, NULL, 1, GrB_ALL, ncols, NULL));
593 :
594 : #if LAGRAPH_SUITESPARSE
595 30 : GRB_TRY(GxB_BinaryOp_new(&MinParent, (void *)minparent, Vertex, Vertex,
596 : Vertex, "minparent", MIN_PARENT_DEFN));
597 : #else
598 : GRB_TRY(GrB_BinaryOp_new(&MinParent, (void *)minparent, Vertex, Vertex,
599 : Vertex));
600 : #endif
601 30 : vertex infinityParent = {GrB_INDEX_MAX + 1, 0};
602 30 : GRB_TRY(GrB_Monoid_new_UDT(&MinParent_Monoid, MinParent, &infinityParent));
603 :
604 : #if LAGRAPH_SUITESPARSE
605 30 : GRB_TRY(GxB_BinaryOp_new(&Select2ndOp, (void *)select2nd, Vertex, GrB_BOOL,
606 : Vertex, "select2nd", SELECT_2ND_DEFN));
607 : #else
608 : GRB_TRY(GrB_BinaryOp_new(&Select2ndOp, (void *)select2nd, Vertex, GrB_BOOL,
609 : Vertex));
610 : #endif
611 :
612 30 : GRB_TRY(GrB_Semiring_new(&MinParent_2nd_Semiring, MinParent_Monoid,
613 : Select2ndOp));
614 :
615 : #if LAGRAPH_SUITESPARSE
616 30 : GRB_TRY(GxB_BinaryOp_new(&Select1stOp, (void *)select1st, Vertex, Vertex,
617 : GrB_BOOL, "select1st", SELECT_1ST_DEFN));
618 : #else
619 : GRB_TRY(GrB_BinaryOp_new(&Select1stOp, (void *)select1st, Vertex, Vertex,
620 : GrB_BOOL));
621 : #endif
622 :
623 30 : GRB_TRY(GrB_Semiring_new(&MinParent_1st_Semiring, MinParent_Monoid,
624 : Select1stOp));
625 :
626 : #if LAGRAPH_SUITESPARSE
627 30 : GRB_TRY(GxB_UnaryOp_new(&getParentsOp, (void *)keepParents, GrB_UINT64,
628 : Vertex, "keepParents", KEEP_PARENTS_DEFN));
629 : #else
630 : GRB_TRY(GrB_UnaryOp_new(&getParentsOp, (void *)keepParents, GrB_UINT64,
631 : Vertex));
632 : #endif
633 :
634 : #if LAGRAPH_SUITESPARSE
635 30 : GRB_TRY(GxB_UnaryOp_new(&getRootsOp, (void *)keepRoots, GrB_UINT64, Vertex,
636 : "keepRoots", KEEP_ROOTS_DEFN));
637 : #else
638 : GRB_TRY(
639 : GrB_UnaryOp_new(&getRootsOp, (void *)keepRoots, GrB_UINT64, Vertex));
640 : #endif
641 :
642 30 : GRB_TRY(GrB_Vector_new(&parentsUpdate, GrB_UINT64, nrows));
643 :
644 30 : GRB_TRY(GrB_Vector_new(&ufrontierR, Vertex, nrows));
645 :
646 30 : GRB_TRY(GrB_Vector_new(&rootsufR, GrB_UINT64, nrows));
647 :
648 30 : GRB_TRY(GrB_Vector_new(&pathUpdate, GrB_UINT64, ncols));
649 :
650 30 : GRB_TRY(GrB_Vector_new(&rootufRIndexes, GrB_UINT64, ncols));
651 :
652 30 : GRB_TRY(GrB_Vector_new(&rootsfR, GrB_UINT64, nrows));
653 :
654 30 : GRB_TRY(GrB_Vector_new(&rootfRIndexes, GrB_UINT64, ncols));
655 :
656 : #if LAGRAPH_SUITESPARSE
657 30 : GRB_TRY(GxB_IndexUnaryOp_new(&buildfCTuplesOp, (void *)buildfCTuples,
658 : Vertex, GrB_UINT64, GrB_BOOL, "buildfCTuples",
659 : BUILT_FC_TUPLES_DEFN));
660 : #else
661 : GRB_TRY(GrB_IndexUnaryOp_new(&buildfCTuplesOp, (void *)buildfCTuples,
662 : Vertex, GrB_UINT64, GrB_BOOL));
663 : #endif
664 :
665 : #if LAGRAPH_SUITESPARSE
666 30 : GRB_TRY(GxB_UnaryOp_new(&vertexTypecastOp, (void *)vertexTypecast, Vertex,
667 : GrB_UINT64, "vertexTypecast",
668 : VERTEX_TYPECAST_DEFN));
669 : #else
670 : GRB_TRY(GrB_UnaryOp_new(&vertexTypecastOp, (void *)vertexTypecast, Vertex,
671 : GrB_UINT64));
672 : #endif
673 :
674 : #if LAGRAPH_SUITESPARSE
675 30 : GRB_TRY(GxB_BinaryOp_new(&setParentsMatesOp, (void *)setParentsMates,
676 : Vertex, Vertex, Vertex, "setParentsMates",
677 : SET_PARENTS_MATES_DEFN));
678 : #else
679 : GRB_TRY(GrB_BinaryOp_new(&setParentsMatesOp, (void *)setParentsMates,
680 : Vertex, Vertex, Vertex));
681 : #endif
682 :
683 30 : GRB_TRY(GrB_Vector_new(&vr, GrB_UINT64, nrows));
684 :
685 30 : GRB_TRY(GrB_Vector_new(&pathCopy, GrB_UINT64, ncols));
686 :
687 30 : GRB_TRY(GrB_Vector_new(¤tMatesR, GrB_UINT64, nrows));
688 :
689 30 : uint64_t npath = 0;
690 30 : bool y = 0;
691 :
692 30 : GRB_TRY(GrB_Vector_new(&mateC, GrB_UINT64, ncols));
693 30 : GRB_TRY(GrB_Vector_new(&mateR, GrB_UINT64, nrows));
694 :
695 30 : if (mate_init != NULL)
696 : {
697 6 : uint64_t nmatched = 0;
698 6 : GRB_TRY(GrB_Vector_nvals(&nmatched, mate_init));
699 6 : if (nmatched)
700 : {
701 6 : if (col_init)
702 : {
703 : // mateC = (uint64_t) mate_init
704 4 : GRB_TRY(GrB_assign(mateC, NULL, NULL, mate_init, GrB_ALL, ncols,
705 : NULL));
706 : // mateR = invert (mateC), but do not clear the input
707 4 : LAGRAPH_TRY(invert_nondestructive(mateR, mateC, msg));
708 : }
709 : else
710 : {
711 : // mateR = (uint64_t) mate_init
712 2 : GRB_TRY(GrB_assign(mateR, NULL, NULL, mate_init, GrB_ALL, nrows,
713 : NULL));
714 : // mateC = invert (mateR), but do not clear the input
715 2 : LAGRAPH_TRY(invert_nondestructive(mateC, mateR, msg));
716 : }
717 : }
718 : }
719 :
720 : do
721 : {
722 126 : GRB_TRY(GrB_Vector_clear(parentsR));
723 : // for every col j not matched, assign f(j) = VERTEX(j,j)
724 126 : GRB_TRY(GrB_Vector_apply_IndexOp_BOOL(
725 : frontierC, mateC, NULL, initFrontierOp, I, true, GrB_DESC_RSC));
726 :
727 126 : uint64_t nfC = 0;
728 :
729 : do
730 : {
731 : //----------------------------------------------------------------------
732 : // STEPS 1,2: Explore neighbors of column frontier (one step of
733 : // BFS), keeping only unvisited rows in the frontierR
734 : //----------------------------------------------------------------------
735 840 : if (do_pushpull)
736 : {
737 : int32_t kind;
738 36 : LAGRAPH_TRY(LG_GET_FORMAT_HINT(frontierC, &kind));
739 36 : if (kind == LG_BITMAP || kind == LG_FULL)
740 : {
741 : // the frontierC vector is bitmap or full
742 : // pull (vector's values are pulled into A)
743 16 : GRB_TRY(GrB_mxv(frontierR, parentsR, NULL,
744 : MinParent_2nd_Semiring, A, frontierC,
745 : GrB_DESC_RSC));
746 : }
747 : else
748 : {
749 : // the frontierC vector is sparse or hypersparse
750 : // push (vector's values are pushed to A)
751 20 : GRB_TRY(GrB_vxm(frontierR, NULL, NULL,
752 : MinParent_1st_Semiring, frontierC, AT,
753 : GrB_DESC_R));
754 20 : GRB_TRY(GrB_Vector_assign(frontierR, parentsR, NULL,
755 : frontierR, GrB_ALL, nrows,
756 : GrB_DESC_RSC));
757 : }
758 : }
759 : else
760 : {
761 804 : if (A != NULL)
762 : {
763 : // Only the pull method can be used if AT is not given
764 524 : GRB_TRY(GrB_mxv(frontierR, parentsR, NULL,
765 : MinParent_2nd_Semiring, A, frontierC,
766 : GrB_DESC_RSC));
767 : }
768 : else
769 : {
770 : // Only the push method can be used if A is not given
771 280 : GRB_TRY(GrB_vxm(frontierR, NULL, NULL,
772 : MinParent_1st_Semiring, frontierC, AT,
773 : GrB_DESC_R));
774 280 : GRB_TRY(GrB_Vector_assign(frontierR, parentsR, NULL,
775 : frontierR, GrB_ALL, nrows,
776 : GrB_DESC_RSC));
777 : }
778 : }
779 :
780 : //----------------------------------------------------------------------
781 : // STEPS 3,4: Select univisited, matched and unmatched row vertices
782 : //----------------------------------------------------------------------
783 : // set parents of row frontier
784 840 : GRB_TRY(GrB_Vector_apply(
785 : parentsR, frontierR, NULL, getParentsOp, frontierR,
786 : GrB_DESC_S)); // use input as mask to only update or insert
787 : // parents without deleting the ones not
788 : // updated
789 :
790 : // select unmatched rows of the R frontier
791 840 : GRB_TRY(GrB_Vector_assign(ufrontierR, mateR, NULL, frontierR,
792 : GrB_ALL, nrows, GrB_DESC_RSC));
793 : // select matched rows of the R frontier
794 840 : GRB_TRY(GrB_Vector_assign(frontierR, mateR, NULL, frontierR,
795 : GrB_ALL, nrows, GrB_DESC_RS));
796 :
797 : // keep only mates of rows in frontierR
798 840 : GRB_TRY(GrB_Vector_assign(currentMatesR, frontierR, NULL, mateR,
799 : GrB_ALL, nrows, GrB_DESC_RS));
800 :
801 840 : uint64_t nUfR = 0, nfR = 0;
802 840 : GRB_TRY(GrB_Vector_nvals(&nUfR, ufrontierR));
803 840 : GRB_TRY(GrB_Vector_nvals(&nfR, frontierR));
804 :
805 840 : if (nUfR)
806 : {
807 : //----------------------------------------------------------------------
808 : // STEP 5: Store endpoints of newly discovered augmenting paths
809 : //----------------------------------------------------------------------
810 : // get roots of unmatched row nodes in the R frontier
811 264 : GRB_TRY(GrB_Vector_apply(rootsufR, NULL, NULL, getRootsOp,
812 : ufrontierR, NULL));
813 :
814 : // pathUpdate = invert (rootsufR), but need to handle
815 : // duplicates
816 264 : LAGRAPH_TRY(invert(pathUpdate, rootsufR, /*dups=*/true, msg));
817 :
818 264 : GRB_TRY(GrB_Vector_assign(
819 : pathC, pathUpdate, NULL, pathUpdate, GrB_ALL, ncols,
820 : GrB_DESC_S)); // update path without deleting the values
821 : // not updated when GrB_ALL is used, ni is
822 : // the number of rows of the vector
823 :
824 : //----------------------------------------------------------------------
825 : // STEP 6: Prune vertices in trees yielding augmenting paths
826 : //----------------------------------------------------------------------
827 264 : GRB_TRY(GrB_Vector_clear(rootfRIndexes));
828 :
829 264 : if (nfR)
830 : {
831 : // get roots of row nodes in the current R frontier
832 234 : GRB_TRY(GrB_Vector_apply(rootsfR, NULL, NULL, getRootsOp,
833 : frontierR, NULL));
834 :
835 : // keep mates and roots of the R frontier (ordered indices)
836 234 : LAGRAPH_TRY(
837 : invert_2(rootfRIndexes, currentMatesR, rootsfR,
838 : /*dups=*/true, /*udt=*/false,
839 : msg)); // rootsfRIndexes(j) = i, where i
840 : // is the col mate of the first
841 : // row included in the current R
842 : // frontier with a col root of j
843 :
844 : // keep only col roots that are not included in ufR
845 234 : GRB_TRY(GrB_Vector_assign(rootfRIndexes, pathUpdate, NULL,
846 : rootfRIndexes, GrB_ALL, ncols,
847 : GrB_DESC_RSC));
848 :
849 : //----------------------------------------------------------------------
850 : // STEP 7a (ufrontierR not empty): Move values in the
851 : // correct positions for the C frontier
852 : //----------------------------------------------------------------------
853 : // rootfRIndexes = invert (rootfRIndexes), so that:
854 : // if LAGRAPH_SUITESPARSE: rootfRIndexes(i) = j,
855 : // where (i,j) = (parentC, rootC) of the new frontier C
856 : // else: rootfRIndexes(i) = j, where i is the matched child
857 : // of the current frontier R that stems from root path of j
858 234 : LAGRAPH_TRY(invert(rootfRIndexes, rootfRIndexes,
859 : /*dups=*/false, msg));
860 : }
861 :
862 : //----------------------------------------------------------------------
863 : // STEP 7b (ufrontierR not empty): Build tuple of (parentC,
864 : // rootC)
865 : //----------------------------------------------------------------------
866 264 : GRB_TRY(GrB_Vector_clear(frontierC));
867 264 : GRB_TRY(GrB_Vector_apply_IndexOp_BOOL(frontierC, NULL, NULL,
868 : buildfCTuplesOp,
869 : rootfRIndexes, true, NULL));
870 : }
871 : else
872 : {
873 : //----------------------------------------------------------------------
874 : // STEP 7a (ufrontierR is empty): Set parents of the R frontier
875 : // to their mates
876 : //----------------------------------------------------------------------
877 : // typecast mateR to ensure domain match with frontier R and
878 : // apply op on frontier to set parents to mates
879 576 : GRB_TRY(
880 : GrB_Vector_apply(frontierR, NULL, setParentsMatesOp,
881 : vertexTypecastOp, currentMatesR,
882 : NULL)); // fR(i) = (column mate of i,
883 : // rootC) add the structural mask
884 :
885 : //----------------------------------------------------------------------
886 : // STEP 7b (ufrontierR is empty): Move values in the correct
887 : // positions for the C frontier
888 : //----------------------------------------------------------------------
889 : // invert fr and assign to fC
890 : // (currentMatesR already contains only the rows of fR)
891 576 : LAGRAPH_TRY(invert_2(frontierC, frontierR, currentMatesR,
892 : /*dups=*/false, /*udt=*/true, msg));
893 : }
894 :
895 840 : GRB_TRY(GrB_Vector_nvals(&nfC, frontierC));
896 :
897 840 : } while (nfC);
898 :
899 : //----------------------------------------------------------------------
900 : // STEP 8: Augment matching by all augmenting paths discovered in
901 : // this phase
902 : //----------------------------------------------------------------------
903 126 : GRB_TRY(GrB_Vector_nvals(&npath, pathC));
904 126 : uint64_t npathCopy = npath;
905 702 : while (npath)
906 : {
907 : // vr = invert (pathC), leaving pathC empty
908 : // pathC doesn't have dup values as it stems from an invertion
909 576 : LAGRAPH_TRY(invert(vr, pathC, /*dups=*/false, msg));
910 :
911 : // assign parents of rows to rows
912 576 : GRB_TRY(GrB_Vector_assign(
913 : vr, vr, NULL, parentsR, GrB_ALL, nrows,
914 : GrB_DESC_S)); // update the values of vr (descriptor needed
915 : // to use mask's structure and not values)
916 :
917 : // update mateR: mateR<vr> = vr
918 576 : GRB_TRY(GrB_Vector_assign(mateR, vr, NULL, vr, GrB_ALL, nrows,
919 : GrB_DESC_S));
920 :
921 : // pathC = invert (vr), leaving vr empty
922 576 : LAGRAPH_TRY(invert(pathC, vr, /*dups=*/false, msg));
923 :
924 : // keep a copy of the previous row matches of the matched cols
925 : // that will alter mates
926 576 : GRB_TRY(GrB_Vector_assign(pathCopy, pathC, NULL, mateC, GrB_ALL,
927 : ncols, GrB_DESC_RS));
928 :
929 : // update mateC
930 576 : GRB_TRY(GrB_Vector_assign(mateC, pathC, NULL, pathC, GrB_ALL, ncols,
931 : GrB_DESC_S));
932 :
933 : // At this point, mateR and mateC are in sync. That is, they
934 : // are inversions of each other (mateR == invert(mateC) and
935 : // mateC == invert (mateR) both hold).
936 :
937 : // swap path and pathCopy
938 576 : GrB_Vector tmp = pathC;
939 576 : pathC = pathCopy;
940 576 : pathCopy = tmp;
941 :
942 576 : GRB_TRY(GrB_Vector_nvals(&npath, pathC));
943 : }
944 :
945 126 : npath = npathCopy;
946 126 : } while (npath); // only in the first and last iteration should this
947 : // condition be false
948 :
949 : // return result
950 30 : if (mateC_handle != NULL)
951 : {
952 30 : (*mateC_handle) = mateC;
953 30 : mateC = NULL ;
954 : }
955 30 : if (mateR_handle != NULL)
956 : {
957 20 : (*mateR_handle) = mateR;
958 20 : mateR = NULL ;
959 : }
960 30 : LG_FREE_WORK;
961 :
962 30 : return (GrB_SUCCESS);
963 : }
|
