diff --git a/experimental/algorithm/LAGr_MaximumMatching.c b/experimental/algorithm/LAGr_MaximumMatching.c new file mode 100644 index 0000000000..f95738f55a --- /dev/null +++ b/experimental/algorithm/LAGr_MaximumMatching.c @@ -0,0 +1,959 @@ +//------------------------------------------------------------------------------ +// LAGr_MaximumMatching: maximum matching between nodes of disjoint sets +// in bipartite graphs +//------------------------------------------------------------------------------ + +// LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved. +// SPDX-License-Identifier: BSD-2-Clause +// +// For additional details (including references to third party source code and +// other files) see the LICENSE file or contact permission@sei.cmu.edu. See +// Contributors.txt for a full list of contributors. Created, in part, with +// funding and support from the U.S. Government (see Acknowledgments.txt file). +// DM22-0790 + +// Contributed by Christina Koutsou, Aristotle University of Thessaloniki + +//------------------------------------------------------------------------------ + +// This implmentation is based on the algorithm described in the following +// paper: "Distributed-Memory Algorithms for Maximum Cardinality Matching in +// Bipartite Graphs" by A. Azad and A. Buluç. +// +// In brief, this algorithm explores alternate paths by reversing an already +// traversed path to see if one of the already matched columns encountered in +// the path has at least one free children to be matched with instead. If so, it +// concedes its previous match to another previously matched parent or to the +// unmatched root of the path. + +// More detailed explanation of the algorithm and its implementation can be +// found in the LAGraph "papers" folder, in file +// "MaximumMatching_Report.pdf". + +// This is an Advanced method, G->A is required as input instead of G due to +// lack of a Bipartite graph kind in GraphBLAS. + +//------------------------------------------------------------------------------ + +#include "LAGraphX.h" +#include "LG_internal.h" + +//------------------------------------------------------------------------------ +// the Vertex tuple: (parentC, rootC) +//------------------------------------------------------------------------------ + +typedef struct +{ + uint64_t parentC; + uint64_t rootC; +} vertex; + +// repeat the typedef as a string, to give to GraphBLAS +#define VERTEX_DEFN \ + "typedef struct " \ + "{ " \ + "uint64_t parentC; " \ + "uint64_t rootC; " \ + "} " \ + "vertex; " + +void initFrontier(vertex *z, void *x, uint64_t i, uint64_t j, const void *y) +{ + z->parentC = i; + z->rootC = i; +} + +#define INIT_FRONTIER_DEFN \ + "void initFrontier(vertex *z, void *x, uint64_t i, uint64_t j, const " \ + "void " \ + "*y) " \ + "{ " \ + "z->parentC = i; " \ + "z->rootC = i; " \ + "} " + +void minparent(vertex *z, vertex *x, vertex *y) +{ + *z = x->parentC < y->parentC ? *x : *y; +} + +#define MIN_PARENT_DEFN \ + "void minparent(vertex *z, vertex *x, vertex *y) " \ + "{ " \ + "*z = x->parentC < y->parentC ? *x : *y; " \ + "} " + +// FIXME: revise GraphBLAS so we can tell it that the select2nd operator +// does not use the 'x' input. +void select2nd(vertex *z, bool *x, vertex *y) +{ + z->parentC = y->parentC; + z->rootC = y->rootC; +} + +#define SELECT_2ND_DEFN \ + "void select2nd(vertex *z, bool *x, vertex *y) " \ + "{ " \ + "z->parentC = y->parentC; " \ + "z->rootC = y->rootC;" \ + "} " + +void select1st(vertex *z, vertex *x, bool *y) +{ + z->parentC = x->parentC; + z->rootC = x->rootC; +} + +#define SELECT_1ST_DEFN \ + "void select1st(vertex *z, vertex *x, bool *y) " \ + "{ " \ + "z->parentC = x->parentC; " \ + "z->rootC = x->rootC;" \ + "} " + +void keepParents(uint64_t *z, vertex *x) { *z = x->parentC; } + +#define KEEP_PARENTS_DEFN \ + "void keepParents(uint64_t *z, vertex *x) { *z = x->parentC; } " + +void keepRoots(uint64_t *z, vertex *x) { *z = x->rootC; } + +#define KEEP_ROOTS_DEFN \ + "void keepRoots(uint64_t *z, vertex *x) { *z = x->rootC; } " + +void buildfCTuples(vertex *z, uint64_t *x, uint64_t i, uint64_t j, + const void *y) +{ + z->parentC = i; + z->rootC = *x; +} + +#define BUILT_FC_TUPLES_DEFN \ + "void buildfCTuples(vertex *z, uint64_t *x, uint64_t i, uint64_t j, " \ + "const void *y) " \ + "{ " \ + "z->parentC = i; " \ + "z->rootC = *x; " \ + "} " + +void vertexTypecast(vertex *z, uint64_t *x) +{ + z->parentC = *x; + z->rootC = *x; +} + +#define VERTEX_TYPECAST_DEFN \ + "void vertexTypecast(vertex *z, uint64_t *x) " \ + "{ " \ + "z->parentC = *x; " \ + "z->rootC = *x; " \ + "} " + +void setParentsMates(vertex *z, vertex *x, vertex *y) +{ + z->parentC = y->parentC; + z->rootC = x->rootC; +} + +#define SET_PARENTS_MATES_DEFN \ + "void setParentsMates(vertex *z, vertex *x, vertex *y) " \ + "{ " \ + "z->parentC = y->parentC; " \ + "z->rootC = x->rootC; " \ + "} " + +//------------------------------------------------------------------------------ +// invert +//------------------------------------------------------------------------------ + +// This function "inverts" an input vector by swapping its row indices +// and its values, returning the result in an output vector. +// +// For example, for the indices/values of an input vector (in) with 5 entries +// and length 100: +// +// indices: 0 3 5 42 99 +// values: 4 98 1 3 12 +// +// on output, the out vector will contain: +// +// indices: 4 98 1 3 12 +// values: 0 3 5 42 99 +// +// The output vector will normally be jumbled since the values will not appear +// in any particular order. The method assumes that the input values are in +// range 0 to n-1 where n = length(out). The values in the input vector +// may be duplicated and this argument of the function must be set accordingly. +// Both the in vector and out vector must have the same type (GrB_UINT64). The +// lengths of the two vectors need not be the same, so long as the indices +// remain in range. Results are undefined if these conditions do not hold. +// +// The in and out vectors may be aliased. If not aliased, the input vector is +// cleared of all entries on output. If in and out are aliased, then the +// inversion is performed in-place. +// +// In SuiteSparse:GraphBLAS, this method takes O(1) time if the in vector is in +// CSC (sparse, by column) format. Otherwise it can take O(e) time if e = +// nvals(in), because the unpack below will convert the in vector to CSC and +// then unpack it. + +#undef LG_FREE_ALL +#define LG_FREE_ALL \ + { \ + LAGraph_Free((void *)&I, NULL); \ + LAGraph_Free((void *)&X1, NULL); \ + LAGraph_Free((void *)&X2, NULL); \ + } + +static inline GrB_Info invert_nondestructive( + GrB_Vector out, // input/output. On input, only the size and type are + // kept; any entries in the 'out' vector cleared. It is + // then replaced with the inversion of the input vector. + GrB_Vector in, // input vector, not modified. There must be no duplicate + // values in the input vector. + char *msg) +{ + bool jumbled = 1; + GrB_Index *I = NULL; + GrB_Index *X1 = NULL; + GrB_Index *X2 = NULL; // not used + GrB_Index IBytes = 0, XBytes = 0; + uint64_t nvals = 0; + + // the input and output vectors cannot be the same vector + ASSERT(in != out); + + // All input/output vectors must be of type GrB_UINT64. +#if LAGRAPH_SUITESPARSE + GRB_TRY( + GxB_Vector_unpack_CSC(in, (GrB_Index **)&I, (void **)&X1, &IBytes, + &XBytes, NULL, &nvals, &jumbled, + NULL)); // the output and input should have no + // duplicates, so the order doesn't matter +#else + GRB_TRY(GrB_Vector_nvals(&nvals, in)); + LG_TRY(LAGraph_Malloc((void **)&I, nvals, sizeof(GrB_Index), msg)); + LG_TRY(LAGraph_Malloc((void **)&X1, nvals, sizeof(GrB_Index), msg)); + GRB_TRY(GrB_Vector_wait(in, GrB_MATERIALIZE)); + GRB_TRY(GrB_Vector_extractTuples_UINT64( + I, X1, &nvals, in)); // the output and input should have no + // duplicates, so the order doesn't matter +#endif + + GRB_TRY(GrB_Vector_clear(out)); // clear the output first as a prerequisite + // of the build method + GRB_TRY(GrB_Vector_build_UINT64( + out, X1, I, nvals, + NULL)); // build does not take ownership of the lists I and X, + // but only copies them, these lists will be given + // again to the input + // the input should have no duplicates in the + // values list, so dups are not handled +#if LAGRAPH_SUITESPARSE + GRB_TRY(GxB_Vector_pack_CSC(in, (GrB_Index **)&I, (void **)&X1, IBytes, + XBytes, NULL, nvals, jumbled, NULL)); +#endif +} + +static inline GrB_Info +invert(GrB_Vector out, // input/output. Same as invert_nondescructive above. + GrB_Vector in, // input vector, empty on output (unless in == out) + const bool dups, // flag that indicates if duplicates exist in the input + // vector's values + char *msg) +{ + // The input and output vectors can be the same vector + // that is, in == out is OK. + // All input/output vectors must be of type GrB_UINT64. + + // the output vector will normally be returned in a jumbled state + bool jumbled = dups ? 0 : 1; // if there are duplicates, we want the indices + // to be ordered so we can keep the min child + GrB_Index *I = NULL; // unpack allocates space for these lists + GrB_Index *X1 = NULL; + GrB_Index *X2 = NULL; // not used + GrB_Index IBytes = 0, XBytes = 0; + uint64_t nvals = 0; + +#if LAGRAPH_SUITESPARSE + GRB_TRY(GxB_Vector_unpack_CSC(in, (GrB_Index **)&I, (void **)&X1, &IBytes, + &XBytes, NULL, &nvals, &jumbled, NULL)); +#else + // vanilla case using extractTuples and build: + // allocate I and X for GrB_extractTuples + GRB_TRY(GrB_Vector_nvals(&nvals, in)); + LG_TRY(LAGraph_Malloc((void **)&I, nvals, sizeof(GrB_Index), msg)); + LG_TRY(LAGraph_Malloc((void **)&X1, nvals, sizeof(GrB_Index), msg)); + GRB_TRY(GrB_Vector_wait(in, GrB_MATERIALIZE)); + GRB_TRY(GrB_Vector_extractTuples_UINT64( + I, X1, &nvals, in)); // the output and input should have no + // duplicates, so the order doesn't matter + GRB_TRY(GrB_Vector_clear(in)); +#endif + if (!dups) + { +#if LAGRAPH_SUITESPARSE + GRB_TRY(GxB_Vector_pack_CSC(out, (GrB_Index **)&X1, (void **)&I, XBytes, + IBytes, NULL, nvals, true, NULL)); +#else + GRB_TRY(GrB_Vector_clear(out)); + // GrB_MIN_UINT64 is used instead of first because first is an extension + GRB_TRY(GrB_Vector_build_UINT64(out, X1, I, nvals, GrB_MIN_UINT64)); + // build copies the lists so they need to be freed in LG_FREE_ALL + LG_FREE_ALL; +#endif + } + else + { + GRB_TRY(GrB_Vector_clear(out)); +#if LAGRAPH_SUITESPARSE + GRB_TRY(GrB_Vector_build_UINT64(out, X1, I, nvals, GrB_FIRST_UINT64)); +#else + GRB_TRY(GrB_Vector_build_UINT64(out, X1, I, nvals, GrB_MIN_UINT64)); +#endif + // build copies the lists so they need to be freed in LG_FREE_ALL + LG_FREE_ALL; + } +} + +static inline GrB_Info +invert_2(GrB_Vector out, // input/output + GrB_Vector in1, // input vector, empty on output (unless in1 == out) + GrB_Vector in2, // input vector, empty on output (unless in2 == out) + const bool dups, // flag that indicates if duplicates exist in the + // input vector's values + const bool udt, // type of the first input, in1, and out vectors + char *msg) +{ + // The input vectors cannot be aliased. However in1==out or in2==out is + // OK. The two input vectors must have the same # of entries. + // All input/output vectors must be of type GrB_UINT64. + ASSERT(in1 != in2); + + GrB_Index *I = NULL; + GrB_Index *X1 = NULL; + GrB_Index *X2 = NULL; + GrB_Index IBytes = 0, X1Bytes = 0, X2Bytes = 0; + uint64_t nvals1 = 0, nvals2 = 0; + +#if LAGRAPH_SUITESPARSE + GRB_TRY(GxB_Vector_unpack_CSC(in1, (GrB_Index **)&I, (void **)&X1, &IBytes, + &X1Bytes, NULL, &nvals1, NULL, NULL)); + LAGraph_Free((void *)&I, NULL); + GRB_TRY(GxB_Vector_unpack_CSC(in2, (GrB_Index **)&I, (void **)&X2, &IBytes, + &X2Bytes, NULL, &nvals2, NULL, NULL)); + ASSERT(nvals1 == nvals2); + if (!dups) + { + LAGraph_Free((void *)&I, NULL); + GRB_TRY(GxB_Vector_pack_CSC(out, (GrB_Index **)&X2, (void **)&X1, + X2Bytes, X1Bytes, NULL, nvals2, true, + NULL)); // the values are not ordered, + // so the indices are jumbled + } + else + { + GRB_TRY(GrB_Vector_clear(out)); + GRB_TRY(GrB_Vector_build_UINT64(out, X2, X1, nvals2, GrB_FIRST_UINT64)); + LG_FREE_ALL; + } +#else + // vanilla case using extractTuples and build: + if (dups) + { + // invert in2 to eliminate dups and get the final indices of out (we + // cannot invert into out if out is of type udt, because we would have a + // domain mismatch with in2) + GrB_Vector in_rev = NULL; + GrB_Index out_size = 0; + GRB_TRY(GrB_Vector_size(&out_size, out)); + GRB_TRY(GrB_Vector_new(&in_rev, GrB_UINT64, out_size)); + LG_TRY(invert(in_rev, in2, true, msg)); + GRB_TRY(GrB_Vector_nvals( + &nvals2, + in_rev)); // out may have less entries than in2 and, thus, in1 + LG_TRY(LAGraph_Malloc((void **)&I, nvals2, sizeof(GrB_Index), msg)); + LG_TRY(LAGraph_Malloc((void **)&X2, nvals2, sizeof(GrB_Index), msg)); + // get indices of out + GRB_TRY(GrB_Vector_extractTuples_UINT64( + X2, I, &nvals2, + in_rev)); // X2 are the final values of in2 and will be the indices + // of out. I are the indices of in2 (which are a subset of + // the indices of in1) + GRB_TRY(GrB_Vector_free(&in_rev)); + } + else + { + GRB_TRY(GrB_Vector_nvals(&nvals2, in2)); + LG_TRY(LAGraph_Malloc((void **)&I, nvals2, sizeof(GrB_Index), msg)); + LG_TRY(LAGraph_Malloc((void **)&X2, nvals2, sizeof(GrB_Index), msg)); + GRB_TRY(GrB_Vector_wait(in2, GrB_MATERIALIZE)); + GRB_TRY(GrB_Vector_extractTuples_UINT64( + I, X2, &nvals2, in2)); // X2 will be the indices of out and I are + // the indices of in2 (which are the indices + // of in1 as well) + } + + GRB_TRY(GrB_Vector_clear(out)); + GRB_TRY(GrB_Vector_wait(in1, GrB_MATERIALIZE)); + if (udt) + { + vertex *X1_V = NULL; + LG_TRY(LAGraph_Malloc((void **)&X1_V, nvals2, sizeof(vertex), msg)); + for (uint64_t i = 0; i < nvals2; i++) + { + GRB_TRY(GrB_Vector_extractElement_UDT( + (void *)&X1_V[i], in1, + (GrB_Index)I[i])); // value I[i] corresponds to pos X2[i] in + // out due to the tuple extraction and X[i] + // will replace it + } + GRB_TRY(GrB_Vector_clear(in1)); + GRB_TRY(GrB_Vector_build_UDT(out, X2, (void *)X1_V, nvals2, + NULL)); // no dups are left + LAGRAPH_TRY(LAGraph_Free((void **)&X1_V, msg)); + } + else + { + LG_TRY(LAGraph_Malloc((void **)&X1, nvals2, sizeof(GrB_UINT64), msg)); + for (GrB_Index i = 0; i < nvals2; i++) + { + GRB_TRY(GrB_Vector_extractElement_UINT64(&X1[i], in1, I[i])); + } + GRB_TRY(GrB_Vector_clear(in1)); + GRB_TRY(GrB_Vector_build_UINT64(out, X2, X1, nvals2, NULL)); + } + LG_FREE_ALL; +#endif +} + +//------------------------------------------------------------------------------ +// LAGr_MaximumMatching +//------------------------------------------------------------------------------ + +#undef LG_FREE_WORK +#define LG_FREE_WORK \ + { \ + GrB_free(&pathC); \ + GrB_free(&parentsR); \ + GrB_free(&Vertex); \ + GrB_free(&frontierC); \ + GrB_free(&frontierR); \ + GrB_free(&initFrontierOp); \ + GrB_free(&I); \ + GrB_free(&MinParent); \ + GrB_free(&MinParent_Monoid); \ + GrB_free(&Select2ndOp); \ + GrB_free(&Select1stOp); \ + GrB_free(&MinParent_2nd_Semiring); \ + GrB_free(&MinParent_1st_Semiring); \ + GrB_free(&getParentsOp); \ + GrB_free(&getRootsOp); \ + GrB_free(&parentsUpdate); \ + GrB_free(&ufrontierR); \ + GrB_free(&mateR); \ + GrB_free(&rootsufR); \ + GrB_free(&pathUpdate); \ + GrB_free(&rootufRIndexes); \ + GrB_free(&rootsfR); \ + GrB_free(&rootfRIndexes); \ + GrB_free(&buildfCTuplesOp); \ + GrB_free(&vertexTypecastOp); \ + GrB_free(&setParentsMatesOp); \ + GrB_free(&vr); \ + GrB_free(&pathCopy); \ + GrB_free(¤tMatesR); \ + } + +#undef LG_FREE_ALL +#define LG_FREE_ALL \ + { \ + LG_FREE_WORK; \ + GrB_free(&mateC); \ + } + +int LAGr_MaximumMatching( + // outputs: + GrB_Vector + *mateC_handle, // mateC(j) = i : Column j of the C subset is matched + // to row i of the R subset (ignored on input) + GrB_Vector *mateR_handle, // mateR(i) = j : Row i of the R subset is matched + // to column j of the C subset (ignored on input) + // inputs: + GrB_Matrix A, // input adjacency matrix, TODO: this should be a LAGraph + // of a BIPARTITE kind + GrB_Matrix AT, // transpose of the input adjacency matrix, necessary to + // perform push-pull optimization + // if NULL, the push-pull optimization is not performed + GrB_Vector mate_init, // input only, not modified, ignored if NULL + bool col_init, // flag to indicate if the initial matching is provided from + // the columns' or from the rows' perspective, ignored if + // mate_init is NULL + char *msg) +{ + + //-------------------------------------------------------------------------- + // check inputs + //-------------------------------------------------------------------------- + + GrB_Vector pathC = NULL; // make this bitmap/sparse, + // if dense I would have to give + // all the entries and make the matrix 1-based + GrB_Vector parentsR = NULL; // parents of row nodes that are reachable + // from paths of the initial column frontier + GrB_Type Vertex = NULL; + GrB_Vector frontierC = NULL; + GrB_Vector frontierR = NULL; + GrB_IndexUnaryOp initFrontierOp = NULL; + GrB_Vector I = NULL; // dense vector of 1's + GrB_BinaryOp MinParent = NULL; + GrB_Monoid MinParent_Monoid = NULL; + GrB_BinaryOp Select2ndOp = NULL; + GrB_BinaryOp Select1stOp = NULL; + GrB_Semiring MinParent_2nd_Semiring = NULL; + GrB_Semiring MinParent_1st_Semiring = NULL; + GrB_UnaryOp getParentsOp = NULL; + GrB_UnaryOp getRootsOp = NULL; + GrB_Vector parentsUpdate = NULL; // unmatched rows of R frontier + GrB_Vector ufrontierR = NULL; // unmatched rows of R frontier + GrB_Vector mateR = NULL; // mateR(i) = j : Row i of the R subset is + // matched to column j of the C subset + GrB_Vector rootsufR = NULL; + GrB_Vector pathUpdate = NULL; + GrB_Vector rootufRIndexes = NULL; + GrB_Vector rootsfR = NULL; + GrB_Vector rootfRIndexes = NULL; + GrB_IndexUnaryOp buildfCTuplesOp = NULL; + GrB_UnaryOp vertexTypecastOp = NULL; + GrB_BinaryOp setParentsMatesOp = NULL; + GrB_Vector vr = NULL; + GrB_Vector pathCopy = NULL; + GrB_Vector currentMatesR = NULL; + + GrB_Vector mateC = NULL; + + LG_CLEAR_MSG; + + LG_ASSERT_MSG(mateC_handle != NULL || mateR_handle != NULL, + GrB_NULL_POINTER, "At least one output must be provided\n"); + LG_ASSERT_MSG(A != NULL || AT != NULL, GrB_NULL_POINTER, + "A matrix is NULL\n"); + + if (mateC_handle != NULL) + { + (*mateC_handle) = NULL; + } + if (mateR_handle != NULL) + { + (*mateR_handle) = NULL; + } + + bool do_pushpull = (AT != NULL) && (A != NULL); + + uint64_t ncols = 0; + uint64_t nrows = 0; + + if (A != NULL) + { + GRB_TRY(GrB_Matrix_ncols(&ncols, A)); + GRB_TRY(GrB_Matrix_nrows(&nrows, A)); + } + else + { + GRB_TRY(GrB_Matrix_nrows(&ncols, AT)); + GRB_TRY(GrB_Matrix_ncols(&nrows, AT)); + } + + GRB_TRY(GrB_Vector_new(&pathC, GrB_UINT64, ncols)); + + GRB_TRY(GrB_Vector_new(&parentsR, GrB_UINT64, nrows)); + +#if LAGRAPH_SUITESPARSE + GRB_TRY(GxB_Type_new(&Vertex, sizeof(vertex), "vertex", VERTEX_DEFN)); +#else + GRB_TRY(GrB_Type_new(&Vertex, sizeof(vertex))); +#endif + + GRB_TRY(GrB_Vector_new(&frontierC, Vertex, ncols)); + + GRB_TRY(GrB_Vector_new(&frontierR, Vertex, nrows)); + +#if LAGRAPH_SUITESPARSE + GRB_TRY(GxB_IndexUnaryOp_new(&initFrontierOp, (void *)initFrontier, Vertex, + GrB_BOOL, GrB_BOOL, "initFrontier", + INIT_FRONTIER_DEFN)); +#else + GRB_TRY(GrB_IndexUnaryOp_new(&initFrontierOp, (void *)initFrontier, Vertex, + GrB_BOOL, GrB_BOOL)); +#endif + + GRB_TRY(GrB_Vector_new(&I, GrB_BOOL, ncols)); + GRB_TRY(GrB_Vector_assign_BOOL(I, NULL, NULL, 1, GrB_ALL, ncols, NULL)); + +#if LAGRAPH_SUITESPARSE + GRB_TRY(GxB_BinaryOp_new(&MinParent, (void *)minparent, Vertex, Vertex, + Vertex, "minparent", MIN_PARENT_DEFN)); +#else + GRB_TRY(GrB_BinaryOp_new(&MinParent, (void *)minparent, Vertex, Vertex, + Vertex)); +#endif + vertex infinityParent = {GrB_INDEX_MAX + 1, 0}; + GRB_TRY(GrB_Monoid_new_UDT(&MinParent_Monoid, MinParent, &infinityParent)); + +#if LAGRAPH_SUITESPARSE + GRB_TRY(GxB_BinaryOp_new(&Select2ndOp, (void *)select2nd, Vertex, GrB_BOOL, + Vertex, "select2nd", SELECT_2ND_DEFN)); +#else + GRB_TRY(GrB_BinaryOp_new(&Select2ndOp, (void *)select2nd, Vertex, GrB_BOOL, + Vertex)); +#endif + + GRB_TRY(GrB_Semiring_new(&MinParent_2nd_Semiring, MinParent_Monoid, + Select2ndOp)); + +#if LAGRAPH_SUITESPARSE + GRB_TRY(GxB_BinaryOp_new(&Select1stOp, (void *)select1st, Vertex, Vertex, + GrB_BOOL, "select1st", SELECT_1ST_DEFN)); +#else + GRB_TRY(GrB_BinaryOp_new(&Select1stOp, (void *)select1st, Vertex, Vertex, + GrB_BOOL)); +#endif + + GRB_TRY(GrB_Semiring_new(&MinParent_1st_Semiring, MinParent_Monoid, + Select1stOp)); + +#if LAGRAPH_SUITESPARSE + GRB_TRY(GxB_UnaryOp_new(&getParentsOp, (void *)keepParents, GrB_UINT64, + Vertex, "keepParents", KEEP_PARENTS_DEFN)); +#else + GRB_TRY(GrB_UnaryOp_new(&getParentsOp, (void *)keepParents, GrB_UINT64, + Vertex)); +#endif + +#if LAGRAPH_SUITESPARSE + GRB_TRY(GxB_UnaryOp_new(&getRootsOp, (void *)keepRoots, GrB_UINT64, Vertex, + "keepRoots", KEEP_ROOTS_DEFN)); +#else + GRB_TRY( + GrB_UnaryOp_new(&getRootsOp, (void *)keepRoots, GrB_UINT64, Vertex)); +#endif + + GRB_TRY(GrB_Vector_new(&parentsUpdate, GrB_UINT64, nrows)); + + GRB_TRY(GrB_Vector_new(&ufrontierR, Vertex, nrows)); + + GRB_TRY(GrB_Vector_new(&rootsufR, GrB_UINT64, nrows)); + + GRB_TRY(GrB_Vector_new(&pathUpdate, GrB_UINT64, ncols)); + + GRB_TRY(GrB_Vector_new(&rootufRIndexes, GrB_UINT64, ncols)); + + GRB_TRY(GrB_Vector_new(&rootsfR, GrB_UINT64, nrows)); + + GRB_TRY(GrB_Vector_new(&rootfRIndexes, GrB_UINT64, ncols)); + +#if LAGRAPH_SUITESPARSE + GRB_TRY(GxB_IndexUnaryOp_new(&buildfCTuplesOp, (void *)buildfCTuples, + Vertex, GrB_UINT64, GrB_BOOL, "buildfCTuples", + BUILT_FC_TUPLES_DEFN)); +#else + GRB_TRY(GrB_IndexUnaryOp_new(&buildfCTuplesOp, (void *)buildfCTuples, + Vertex, GrB_UINT64, GrB_BOOL)); +#endif + +#if LAGRAPH_SUITESPARSE + GRB_TRY(GxB_UnaryOp_new(&vertexTypecastOp, (void *)vertexTypecast, Vertex, + GrB_UINT64, "vertexTypecast", + VERTEX_TYPECAST_DEFN)); +#else + GRB_TRY(GrB_UnaryOp_new(&vertexTypecastOp, (void *)vertexTypecast, Vertex, + GrB_UINT64)); +#endif + +#if LAGRAPH_SUITESPARSE + GRB_TRY(GxB_BinaryOp_new(&setParentsMatesOp, (void *)setParentsMates, + Vertex, Vertex, Vertex, "setParentsMates", + SET_PARENTS_MATES_DEFN)); +#else + GRB_TRY(GrB_BinaryOp_new(&setParentsMatesOp, (void *)setParentsMates, + Vertex, Vertex, Vertex)); +#endif + + GRB_TRY(GrB_Vector_new(&vr, GrB_UINT64, nrows)); + + GRB_TRY(GrB_Vector_new(&pathCopy, GrB_UINT64, ncols)); + + GRB_TRY(GrB_Vector_new(¤tMatesR, GrB_UINT64, nrows)); + + uint64_t npath = 0; + bool y = 0; + + GRB_TRY(GrB_Vector_new(&mateC, GrB_UINT64, ncols)); + GRB_TRY(GrB_Vector_new(&mateR, GrB_UINT64, nrows)); + + if (mate_init != NULL) + { + uint64_t nmatched = 0; + GRB_TRY(GrB_Vector_nvals(&nmatched, mate_init)); + if (nmatched) + { + if (col_init) + { + // mateC = (uint64_t) mate_init + GRB_TRY(GrB_assign(mateC, NULL, NULL, mate_init, GrB_ALL, ncols, + NULL)); + // mateR = invert (mateC), but do not clear the input + LAGRAPH_TRY(invert_nondestructive(mateR, mateC, msg)); + } + else + { + // mateR = (uint64_t) mate_init + GRB_TRY(GrB_assign(mateR, NULL, NULL, mate_init, GrB_ALL, nrows, + NULL)); + // mateC = invert (mateR), but do not clear the input + LAGRAPH_TRY(invert_nondestructive(mateC, mateR, msg)); + } + } + } + + do + { + GRB_TRY(GrB_Vector_clear(parentsR)); + // for every col j not matched, assign f(j) = VERTEX(j,j) + GRB_TRY(GrB_Vector_apply_IndexOp_UDT( + frontierC, mateC, NULL, initFrontierOp, I, &y, GrB_DESC_RSC)); + + uint64_t nfC = 0; + + do + { + //---------------------------------------------------------------------- + // STEPS 1,2: Explore neighbors of column frontier (one step of + // BFS), keeping only unvisited rows in the frontierR + //---------------------------------------------------------------------- + if (do_pushpull) + { + int32_t kind; + LAGRAPH_TRY(LG_GET_FORMAT_HINT(frontierC, &kind)); + if (kind == LG_BITMAP || kind == LG_FULL) + { + // the frontierC vector is bitmap or full + // pull (vector's values are pulled into A) + GRB_TRY(GrB_mxv(frontierR, parentsR, NULL, + MinParent_2nd_Semiring, A, frontierC, + GrB_DESC_RSC)); + } + else + { + // the frontierC vector is sparse or hypersparse + // push (vector's values are pushed to A) + GRB_TRY(GrB_vxm(frontierR, NULL, NULL, + MinParent_1st_Semiring, frontierC, AT, + GrB_DESC_R)); + GRB_TRY(GrB_Vector_assign(frontierR, parentsR, NULL, + frontierR, GrB_ALL, nrows, + GrB_DESC_RSC)); + } + } + else + { + if (A != NULL) + { + // Only the pull method can be used if AT is not given + GRB_TRY(GrB_mxv(frontierR, parentsR, NULL, + MinParent_2nd_Semiring, A, frontierC, + GrB_DESC_RSC)); + } + else + { + // Only the push method can be used if A is not given + GRB_TRY(GrB_vxm(frontierR, NULL, NULL, + MinParent_1st_Semiring, frontierC, AT, + GrB_DESC_R)); + GRB_TRY(GrB_Vector_assign(frontierR, parentsR, NULL, + frontierR, GrB_ALL, nrows, + GrB_DESC_RSC)); + } + } + + //---------------------------------------------------------------------- + // STEPS 3,4: Select univisited, matched and unmatched row vertices + //---------------------------------------------------------------------- + // set parents of row frontier + GRB_TRY(GrB_Vector_apply( + parentsR, frontierR, NULL, getParentsOp, frontierR, + GrB_DESC_S)); // use input as mask to only update or insert + // parents without deleting the ones not + // updated + + // select unmatched rows of the R frontier + GRB_TRY(GrB_Vector_assign(ufrontierR, mateR, NULL, frontierR, + GrB_ALL, nrows, GrB_DESC_RSC)); + // select matched rows of the R frontier + GRB_TRY(GrB_Vector_assign(frontierR, mateR, NULL, frontierR, + GrB_ALL, nrows, GrB_DESC_RS)); + + // keep only mates of rows in frontierR + GRB_TRY(GrB_Vector_assign(currentMatesR, frontierR, NULL, mateR, + GrB_ALL, nrows, GrB_DESC_RS)); + + uint64_t nUfR = 0, nfR = 0; + GRB_TRY(GrB_Vector_nvals(&nUfR, ufrontierR)); + GRB_TRY(GrB_Vector_nvals(&nfR, frontierR)); + + if (nUfR) + { + //---------------------------------------------------------------------- + // STEP 5: Store endpoints of newly discovered augmenting paths + //---------------------------------------------------------------------- + // get roots of unmatched row nodes in the R frontier + GRB_TRY(GrB_Vector_apply(rootsufR, NULL, NULL, getRootsOp, + ufrontierR, NULL)); + + // pathUpdate = invert (rootsufR), but need to handle + // duplicates + LAGRAPH_TRY(invert(pathUpdate, rootsufR, /*dups=*/true, msg)); + + GRB_TRY(GrB_Vector_assign( + pathC, pathUpdate, NULL, pathUpdate, GrB_ALL, ncols, + GrB_DESC_S)); // update path without deleting the values + // not updated when GrB_ALL is used, ni is + // the number of rows of the vector + + //---------------------------------------------------------------------- + // STEP 6: Prune vertices in trees yielding augmenting paths + //---------------------------------------------------------------------- + GRB_TRY(GrB_Vector_clear(rootfRIndexes)); + + if (nfR) + { + // get roots of row nodes in the current R frontier + GRB_TRY(GrB_Vector_apply(rootsfR, NULL, NULL, getRootsOp, + frontierR, NULL)); + + // keep mates and roots of the R frontier (ordered indices) + LAGRAPH_TRY( + invert_2(rootfRIndexes, currentMatesR, rootsfR, + /*dups=*/true, /*udt=*/false, + msg)); // rootsfRIndexes(j) = i, where i + // is the col mate of the first + // row included in the current R + // frontier with a col root of j + + // keep only col roots that are not included in ufR + GRB_TRY(GrB_Vector_assign(rootfRIndexes, pathUpdate, NULL, + rootfRIndexes, GrB_ALL, ncols, + GrB_DESC_RSC)); + + //---------------------------------------------------------------------- + // STEP 7a (ufrontierR not empty): Move values in the + // correct positions for the C frontier + //---------------------------------------------------------------------- + // rootfRIndexes = invert (rootfRIndexes), so that: + // if LAGRAPH_SUITESPARSE: rootfRIndexes(i) = j, + // where (i,j) = (parentC, rootC) of the new frontier C + // else: rootfRIndexes(i) = j, where i is the matched child + // of the current frontier R that stems from root path of j + LAGRAPH_TRY(invert(rootfRIndexes, rootfRIndexes, + /*dups=*/false, msg)); + } + + //---------------------------------------------------------------------- + // STEP 7b (ufrontierR not empty): Build tuple of (parentC, + // rootC) + //---------------------------------------------------------------------- + GRB_TRY(GrB_Vector_clear(frontierC)); + GRB_TRY(GrB_Vector_apply_IndexOp_UDT(frontierC, NULL, NULL, + buildfCTuplesOp, + rootfRIndexes, &y, NULL)); + } + else + { + //---------------------------------------------------------------------- + // STEP 7a (ufrontierR is empty): Set parents of the R frontier + // to their mates + //---------------------------------------------------------------------- + // typecast mateR to ensure domain match with frontier R and + // apply op on frontier to set parents to mates + GRB_TRY( + GrB_Vector_apply(frontierR, NULL, setParentsMatesOp, + vertexTypecastOp, currentMatesR, + NULL)); // fR(i) = (column mate of i, + // rootC) add the structural mask + + //---------------------------------------------------------------------- + // STEP 7b (ufrontierR is empty): Move values in the correct + // positions for the C frontier + //---------------------------------------------------------------------- + // invert fr and assign to fC + // (currentMatesR already contains only the rows of fR) + LAGRAPH_TRY(invert_2(frontierC, frontierR, currentMatesR, + /*dups=*/false, /*udt=*/true, msg)); + } + + GRB_TRY(GrB_Vector_nvals(&nfC, frontierC)); + + } while (nfC); + + //---------------------------------------------------------------------- + // STEP 8: Augment matching by all augmenting paths discovered in + // this phase + //---------------------------------------------------------------------- + GRB_TRY(GrB_Vector_nvals(&npath, pathC)); + uint64_t npathCopy = npath; + while (npath) + { + // vr = invert (pathC), leaving pathC empty + // pathC doesn't have dup values as it stems from an invertion + LAGRAPH_TRY(invert(vr, pathC, /*dups=*/false, msg)); + + // assign parents of rows to rows + GRB_TRY(GrB_Vector_assign( + vr, vr, NULL, parentsR, GrB_ALL, nrows, + GrB_DESC_S)); // update the values of vr (descriptor needed + // to use mask's structure and not values) + + // update mateR: mateR = vr + GRB_TRY(GrB_Vector_assign(mateR, vr, NULL, vr, GrB_ALL, nrows, + GrB_DESC_S)); + + // pathC = invert (vr), leaving vr empty + LAGRAPH_TRY(invert(pathC, vr, /*dups=*/false, msg)); + + // keep a copy of the previous row matches of the matched cols + // that will alter mates + GRB_TRY(GrB_Vector_assign(pathCopy, pathC, NULL, mateC, GrB_ALL, + ncols, GrB_DESC_RS)); + + // update mateC + GRB_TRY(GrB_Vector_assign(mateC, pathC, NULL, pathC, GrB_ALL, ncols, + GrB_DESC_S)); + + // At this point, mateR and mateC are in sync. That is, they + // are inversions of each other (mateR == invert(mateC) and + // mateC == invert (mateR) both hold). + + // swap path and pathCopy + GrB_Vector tmp = pathC; + pathC = pathCopy; + pathCopy = tmp; + + GRB_TRY(GrB_Vector_nvals(&npath, pathC)); + } + + npath = npathCopy; + } while (npath); // only in the first and last iteration should this + // condition be false + + if (mateC_handle != NULL) + { + (*mateC_handle) = mateC; + } + if (mateR_handle != NULL) + { + (*mateR_handle) = mateR; + } + LG_FREE_WORK; + + return (GrB_SUCCESS); +} diff --git a/experimental/benchmark/mcm_demo.c b/experimental/benchmark/mcm_demo.c new file mode 100644 index 0000000000..6e4c72ec36 --- /dev/null +++ b/experimental/benchmark/mcm_demo.c @@ -0,0 +1,273 @@ +//------------------------------------------------------------------------------ +// LAGraph/experimental/benchmark/matching_demo.c: benchmarks for +// LAGr_MaximumMatching +//------------------------------------------------------------------------------ + +// LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved. +// SPDX-License-Identifier: BSD-2-Clause +// +// For additional details (including references to third party source code and +// other files) see the LICENSE file or contact permission@sei.cmu.edu. See +// Contributors.txt for a full list of contributors. Created, in part, with +// funding and support from the U.S. Government (see Acknowledgments.txt file). +// DM22-0790 + +// Contributed by Christina Koutsou, Aristotle University of Thessaloniki + +/* +Usage: + +*/ + +//------------------------------------------------------------------------------ + +#include "../../src/benchmark/LAGraph_demo.h" +#include "LAGraphX.h" +#include "LG_internal.h" +#include + +// #define VERBOSE + +#define NTHREAD_LIST 1 +#define THREAD_LIST 0 + +#define LG_FREE_ALL \ + { \ + GrB_free(&A); \ + GrB_free(&M); \ + GrB_free(&mateC); \ + LAGraph_Free((void *)&I, NULL); \ + LAGraph_Free((void *)&X, NULL); \ + } + +GrB_Info check_matching(GrB_Matrix A, GrB_Vector mateC, char *msg) +{ + GrB_Index nmatched = 0; + GrB_Vector mateR = NULL; + GrB_Matrix M = NULL; + GrB_Index *I = NULL, *X = NULL; + + uint64_t ncols = 0, nrows = 0; + GRB_TRY(GrB_Matrix_ncols(&ncols, A)); + GRB_TRY(GrB_Matrix_nrows(&nrows, A)); + + // invert to check for dups + GrB_Index IBytes = 0, XBytes = 0; + bool jumbled = 1; + GRB_TRY(GrB_Vector_new(&mateR, GrB_UINT64, nrows)); + GRB_TRY(GxB_Vector_unpack_CSC(mateC, (GrB_Index **)&I, (void **)&X, &IBytes, + &XBytes, NULL, &nmatched, &jumbled, NULL)); + GRB_TRY(GrB_Vector_build_UINT64(mateR, X, I, nmatched, GrB_FIRST_UINT64)); + GrB_Index nmateR = 0; + GRB_TRY(GrB_Vector_nvals(&nmateR, mateR)); + // if nvals of mateC and mateR don't match, then there's at least + // one row that is used in at least one matching + if (nmatched != nmateR) + { + printf("Duplicates in mateC"); + fflush(stdout); + abort(); + } + + // pack matched values in a matrix + bool *val; + LAGRAPH_TRY(LAGraph_Malloc((void **)&val, nmatched, sizeof(bool), msg)); + for (uint64_t i = 0; i < nmatched; i++) + val[i] = 1; + GRB_TRY(GrB_Matrix_new(&M, GrB_BOOL, nrows, ncols)); + GRB_TRY(GrB_Matrix_build_BOOL(M, X, I, val, nmatched, NULL)); + LAGRAPH_TRY(LAGraph_Free((void **)&val, msg)); + // mask with matrix A to check if all edges are present in A + GRB_TRY(GrB_Matrix_assign(M, M, NULL, A, GrB_ALL, nrows, GrB_ALL, ncols, + GrB_DESC_S)); + GrB_Index nvalsM = 0; + GRB_TRY(GrB_Matrix_nvals(&nvalsM, M)); + // if values have been eliminated then edges do not exist in A + if (nvalsM != nmatched) + { + printf("mateC invalid!\n"); + fflush(stdout); + abort(); + } + + GRB_TRY(GxB_Vector_pack_CSC(mateC, (GrB_Index **)&I, (void **)&X, IBytes, + XBytes, NULL, nmatched, jumbled, NULL)); + + GrB_Vector_free(&mateR); + GrB_Matrix_free(&M); +} + +#undef LG_FREE_ALL +#define LG_FREE_ALL \ + { \ + LAGraph_Delete(&G, NULL); \ + GrB_free(&A); \ + GrB_free(&AT); \ + GrB_free(&mateC_init); \ + GrB_free(&mateC); \ + } + +int main(int argc, char **argv) +{ + //-------------------------------------------------------------------------- + // declare inputs and outputs + //-------------------------------------------------------------------------- + + char msg[LAGRAPH_MSG_LEN]; + + LAGraph_Graph G = NULL; + GrB_Matrix A = NULL; + GrB_Matrix AT = NULL; + GrB_Vector mateC_init = NULL; + GrB_Vector mateC = NULL; + + //-------------------------------------------------------------------------- + // startup LAGraph and GraphBLAS + //-------------------------------------------------------------------------- + + bool burble = false; + demo_init(burble); + + //-------------------------------------------------------------------------- + // read in the graph + //-------------------------------------------------------------------------- + + if (argc < 2) + { + printf("Invalid usage, please read comments\n"); + fflush(stdout); + return 0; + } + char *matrix_name = (argc > 1) ? argv[1] : "stdin"; + + LAGRAPH_TRY(LAGraph_Random_Init(msg)); + bool make_symmetric = false, remove_self_edges = false, structural = true, + ensure_positive = false; + LAGRAPH_TRY(readproblem(&G, NULL, make_symmetric, remove_self_edges, + structural, NULL, ensure_positive, argc, argv)); + + A = G->A; + // compute AT to be able to use push-pull optimization + if (G->is_symmetric_structure) + AT = A; + else + { + LAGRAPH_TRY(LAGraph_Cached_AT(G, msg)); + AT = G->AT; + } + + //-------------------------------------------------------------------------- + // determine the number of threads to run the algorithm with + //-------------------------------------------------------------------------- + + int nt = NTHREAD_LIST; + int Nthreads[20] = {0, THREAD_LIST}; + int nthreads_max, nthreads_outer, nthreads_inner; + LAGRAPH_TRY(LAGraph_GetNumThreads(&nthreads_outer, &nthreads_inner, msg)); +#ifdef VERBOSE + printf("nthreads_outer: %d, nthreads_inner: %d\n", nthreads_outer, + nthreads_inner); + fflush(stdout); +#endif + nthreads_max = nthreads_outer * nthreads_inner; + if (Nthreads[1] == 0) // THREAD_LIST == 0 + { + // create thread list automatically + Nthreads[1] = nthreads_max; + for (int t = 2; t <= nt; t++) + { + Nthreads[t] = Nthreads[t - 1] / 2; + if (Nthreads[t] == 0) + nt = t - 1; + } + } +#ifdef VERBOSE + printf("threads to test: "); + for (int t = 1; t <= nt; t++) + { + int nthreads = Nthreads[t]; + if (nthreads > nthreads_max) + continue; + printf(" %d", nthreads); + } + printf("\n"); + fflush(stdout); +#endif + + //-------------------------------------------------------------------------- + // warmup before benchmarking + //-------------------------------------------------------------------------- + + double t = LAGraph_WallClockTime(); + LAGRAPH_TRY( + LAGr_MaximumMatching(&mateC, NULL, A, AT, mateC_init, true, msg)); + t = LAGraph_WallClockTime() - t; + LAGRAPH_TRY(check_matching(A, mateC, msg)); + uint64_t sprank = 0; + GRB_TRY(GrB_Vector_nvals(&sprank, mateC)); + printf("number of matches: %ld\n", sprank); + fflush(stdout); + GRB_TRY(GrB_free(&mateC)); +#ifdef VERBOSE + printf("warmup time %g sec\n", t); + fflush(stdout); +#endif + + //-------------------------------------------------------------------------- + // benchmark + //-------------------------------------------------------------------------- + + // the GAP benchmark requires 16 trials + int ntrials = 16; +#ifdef VERBOSE + printf("# of trials: %d\n", ntrials); + fflush(stdout); +#endif + + for (int kk = 1; kk <= nt; kk++) + { + int nthreads = Nthreads[kk]; + if (nthreads > nthreads_max) + continue; + LAGRAPH_TRY(LAGraph_SetNumThreads(1, nthreads, msg)); + +#ifdef VERBOSE + printf("\n--------------------------- nthreads: %2d\n", nthreads); + fflush(stdout); +#endif + + double total_time = 0; + + for (int trial = 0; trial < ntrials; trial++) + { + t = LAGraph_WallClockTime(); + LAGRAPH_TRY(LAGr_MaximumMatching(&mateC, NULL, A, AT, mateC_init, + true, msg)); + t = LAGraph_WallClockTime() - t; + GRB_TRY(GrB_free(&mateC)); +#ifdef VERBOSE + printf("trial: %2d time: %10.7f sec\n", trial, t); + fflush(stdout); +#endif + total_time += t; + } + + double total_time_per_trial = total_time / ntrials; + +#ifndef VERBOSE + printf("%.7f\n", total_time_per_trial); +#else + printf("maximum matching: %3d: avg time: %10.7f (sec) matrix: %s\n", + nthreads, total_time_per_trial, matrix_name); +#endif + fflush(stdout); + } + + //-------------------------------------------------------------------------- + // free all workspace and finish + //-------------------------------------------------------------------------- + LG_FREE_ALL; + + LAGRAPH_TRY(LAGraph_Finalize(msg)); + return (GrB_SUCCESS); +} diff --git a/experimental/test/test_MaximumMatching.c b/experimental/test/test_MaximumMatching.c new file mode 100644 index 0000000000..0cdbf7bc07 --- /dev/null +++ b/experimental/test/test_MaximumMatching.c @@ -0,0 +1,142 @@ +#include +#include + +#include "LG_internal.h" +#include +#include +#include + +char msg[LAGRAPH_MSG_LEN]; +LAGraph_Graph G = NULL; + +#define LEN 512 +char filename[LEN + 1]; + +#define NTESTS 5 + +const char *filenames[NTESTS] = {"random_weighted_bipartite2.mtx", + "test_FW_2500.mtx", "LFAT5_hypersparse.mtx", + "lp_afiro_structure.mtx", "sources_7.mtx"}; +const uint64_t spranks[NTESTS] = {298, 2009, 14, 27, 1}; + +void test_MCM(void) +{ + LAGraph_Init(msg); + + OK(LG_SET_BURBLE(1)); + + for (uint8_t jit = 0; jit < 2; jit++) + { + uint8_t JIT_flag = jit * 4; // JIT_OFF = 0 and JIT_ON = 4 + OK(GxB_Global_Option_set(GxB_JIT_C_CONTROL, JIT_flag)); + for (uint64_t test = 0; test < NTESTS; test++) + { + + GrB_Matrix A = NULL; + GrB_Matrix AT = NULL; + snprintf(filename, LEN, LG_DATA_DIR "%s", filenames[test]); + FILE *f = fopen(filename, "r"); + TEST_CHECK(f != NULL); + OK(LAGraph_MMRead(&A, f, msg)); + OK(fclose(f)); + GrB_Index nrows = 0, ncols = 0, nvals = 0; + OK(GrB_Matrix_nrows(&nrows, A)); + OK(GrB_Matrix_ncols(&ncols, A)); + OK(GrB_Matrix_nvals(&nvals, A)); + + // make A a bool matrix and iso-valued + GrB_Index *I, *J, *X; + double *dummy; + bool *iso_value; + + OK(LAGraph_Malloc((void **)&I, nvals, sizeof(GrB_Index), msg)); + OK(LAGraph_Malloc((void **)&J, nvals, sizeof(GrB_Index), msg)); + OK(LAGraph_Malloc((void **)&dummy, nvals, sizeof(double), msg)); + OK(LAGraph_Malloc((void **)&iso_value, nvals, sizeof(bool), msg)); + + for (uint64_t i = 0; i < nvals; i++) + iso_value[i] = 1; + OK(GrB_Matrix_extractTuples_FP64(I, J, dummy, &nvals, A)); + TEST_CHECK(I != NULL); + OK(GrB_Matrix_new(&A, GrB_BOOL, nrows, ncols)); + OK(GrB_Matrix_build_BOOL(A, I, J, iso_value, nvals, + GrB_FIRST_BOOL)); + + OK(LAGraph_Free((void **)&I, msg)); + OK(LAGraph_Free((void **)&J, msg)); + OK(LAGraph_Free((void **)&dummy, msg)); + OK(LAGraph_Free((void **)&iso_value, msg)); + + GrB_Vector mateC = NULL; + OK(GrB_Vector_new(&mateC, GrB_UINT64, ncols)); + + GrB_Vector mateC_init = NULL; + + if (filenames[test] == "lp_afiro_structure.mtx") + { + OK(GrB_Vector_new(&mateC_init, GrB_UINT64, ncols)); + OK(GrB_Vector_setElement_UINT64( + mateC_init, 0, 19)); // col 20 matched with row 1 (1-based) + OK(GrB_Matrix_new(&AT, GrB_BOOL, ncols, + nrows)); // transpose matrix has the reverse + // dimensions from the original + OK(GrB_transpose(AT, NULL, NULL, A, NULL)); + } + + OK(LAGr_MaximumMatching(&mateC, NULL, A, AT, mateC_init, true, + msg)); + printf("\nmsg: %s\n", msg); + + GrB_Index nmatched = 0; + + GrB_Vector mateR = NULL; + OK(GrB_Vector_new(&mateR, GrB_UINT64, nrows)); + + // invert to check for dups + GrB_Index Ibytes = 0, Jbytes = 0, Xbytes = 0; + bool jumbled = 1; + OK(GxB_Vector_unpack_CSC(mateC, (GrB_Index **)&J, (void **)&X, + &Jbytes, &Xbytes, NULL, &nmatched, + &jumbled, NULL)); + OK(GrB_Vector_build_UINT64(mateR, X, J, nmatched, + GrB_FIRST_UINT64)); + GrB_Index nmateR = 0; + OK(GrB_Vector_nvals(&nmateR, mateR)); + // if nvals of mateC and mateR don't match, then there's at least + // one row that is used in at least one matching + TEST_CHECK(nmatched == nmateR); + + // pack matched values in a matrix + GrB_Matrix M = NULL; + bool *val; + OK(LAGraph_Malloc((void **)&val, nmatched, sizeof(bool), msg)); + for (uint64_t i = 0; i < nmatched; i++) + val[i] = 1; + OK(GrB_Matrix_new(&M, GrB_BOOL, nrows, ncols)); + OK(GrB_Matrix_build_BOOL(M, X, J, val, nmatched, NULL)); + OK(LAGraph_Free((void **)&val, msg)); + OK(LAGraph_Free((void **)&J, msg)); + OK(LAGraph_Free((void **)&X, msg)); + // mask with matrix A to check if all edges are present in A + OK(GrB_Matrix_assign(M, M, NULL, A, GrB_ALL, nrows, GrB_ALL, ncols, + GrB_DESC_S)); + GrB_Index nvalsM = 0; + OK(GrB_Matrix_nvals(&nvalsM, M)); + // if values have been eliminated then edges do not exist in A + TEST_CHECK(nvalsM == nmatched); + + // sprank must be equal to nvals of mateC (nmatched) + TEST_CHECK(nmatched == spranks[test]); + + OK(GrB_Vector_free(&mateC)); + OK(GrB_Vector_free(&mateR)); + OK(GrB_Matrix_free(&M)); + OK(GrB_Matrix_free(&A)); + OK(GrB_Matrix_free(&AT)); + } + } + LAGraph_Finalize(msg); +} + +TEST_LIST = {{"MaximumMatching", test_MCM}, // just one test in this example + {NULL, NULL}}; \ No newline at end of file diff --git a/include/LAGraphX.h b/include/LAGraphX.h index 9ad988c166..2d0c40901f 100644 --- a/include/LAGraphX.h +++ b/include/LAGraphX.h @@ -1229,8 +1229,26 @@ int LAGraph_argminmax int dim, // dim=1: cols of A, dim=2: rows of A bool is_min, char *msg -); +); +LAGRAPHX_PUBLIC +int LAGr_MaximumMatching( + // outputs + GrB_Vector + *mateC_handle, // mateC(j) = i : Column j of the C subset is matched to + // row i of the R subset (ignored on input) + GrB_Vector *mateR_handle, // mateR(i) = j : Row i of the R subset is matched + // to column j of the C subset (ignored on input) + // inputs + GrB_Matrix A, // input adjacency matrix, TODO: this should be a LAGraph of a + // BIPARTITE kind + GrB_Matrix + AT, // trasnpose of the input adjacency matrix, NULL if not provided + GrB_Vector mate_init, // input only, not modified, ignored if NULL + bool col_init, // flag to indicate if the initial matching is provided from + // the columns' or from the rows' perspective, ignored if + // mate_init is NULL + char *msg); #if defined ( __cplusplus ) } diff --git a/papers/Maximum-matching-doc.pdf b/papers/Maximum-matching-doc.pdf new file mode 100644 index 0000000000..47ba343349 Binary files /dev/null and b/papers/Maximum-matching-doc.pdf differ diff --git a/papers/Maximum-matching-poster.pdf b/papers/Maximum-matching-poster.pdf new file mode 100644 index 0000000000..f2c4c7cf10 Binary files /dev/null and b/papers/Maximum-matching-poster.pdf differ