From 39c6f7e8079fbca8961006bfe2732cd4afbbed9b Mon Sep 17 00:00:00 2001 From: GomezGab Date: Thu, 6 Aug 2026 17:45:31 -0500 Subject: [PATCH 01/13] Improve performance for calculating modularity metric --- .../algorithm/LAGr_ModularityMatrix.c | 45 +++++++++---------- 1 file changed, 21 insertions(+), 24 deletions(-) diff --git a/experimental/algorithm/LAGr_ModularityMatrix.c b/experimental/algorithm/LAGr_ModularityMatrix.c index 0482959717..78ccc2787b 100644 --- a/experimental/algorithm/LAGr_ModularityMatrix.c +++ b/experimental/algorithm/LAGr_ModularityMatrix.c @@ -27,6 +27,7 @@ // Current Test File: experimental/test/test_modularity.c +#include "GraphBLAS.h" #include "LG_internal.h" #include #include @@ -34,11 +35,7 @@ #define LG_FREE_MOD \ { \ GrB_free(&k); \ - GrB_free(&kk_); \ GrB_free(&B); \ - GrB_free(&BS); \ - GrB_free(&S_BS); \ - GrB_free(&Diag); \ } #undef LG_FREE_ALL #define LG_FREE_ALL \ @@ -78,29 +75,20 @@ int LAGr_AdjModularity( // GrB_set(GrB_GLOBAL, false, GxB_BURBLE); GrB_Index n; - GrB_Matrix k = NULL; - GrB_Matrix kk_ = NULL; + GrB_Vector k = NULL; GrB_Matrix B = NULL; - GrB_Matrix BS = NULL; - GrB_Matrix S_BS = NULL; - GrB_Matrix Diag = NULL; - GrB_Matrix mask = NULL; GRB_TRY(GrB_Matrix_nrows(&n, A)); GRB_TRY(GrB_Matrix_new(&B, GrB_FP64, n, n)); - GRB_TRY(GrB_Matrix_new(&k, GrB_FP64, n, 1)); - GRB_TRY(GrB_Matrix_new(&kk_, GrB_FP64, n, n)); - GRB_TRY(GrB_Matrix_new(&BS, GrB_FP64, n, n)); - GRB_TRY(GrB_Matrix_new(&S_BS, GrB_FP64, n, n)); - GRB_TRY(GrB_Matrix_new(&Diag, GrB_FP64, n, n)); + GRB_TRY(GrB_Vector_new(&k, GrB_FP64, n)); double m = 0.0; double Q_ = 0.0; - GRB_TRY(GrB_Matrix_reduce_Monoid((GrB_Vector)k, NULL, NULL, GrB_PLUS_MONOID_FP64, A, NULL)); + GRB_TRY (GrB_reduce (k, NULL, NULL, GrB_PLUS_MONOID_FP64, A, NULL)); + GRB_TRY (GrB_reduce (&m, NULL, GrB_PLUS_MONOID_FP64, k, NULL)); - GRB_TRY(GrB_Matrix_reduce_FP64(&m, GrB_PLUS_FP64, GrB_PLUS_MONOID_FP64, A, NULL)); m /= 2.0; if (m == 0.0) { @@ -108,16 +96,25 @@ int LAGr_AdjModularity( LG_FREE_ALL; return GrB_SUCCESS; } - GRB_TRY(GrB_mxm(kk_, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_FP64, k, k, GrB_DESC_T1)); double inv_m = -gamma / (2.0 * m); - GRB_TRY(GrB_Matrix_apply_BinaryOp2nd_FP64(kk_, NULL, NULL, GrB_TIMES_FP64, kk_, inv_m, GrB_DESC_R)); - GRB_TRY(GrB_eWiseAdd(B, NULL, NULL, GrB_PLUS_FP64, A, kk_, NULL)); - GRB_TRY(GrB_mxm(BS, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_FP64, B, S, NULL)); - GRB_TRY(GrB_mxm(S_BS, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_FP64, S, BS, GrB_DESC_T0)); - GRB_TRY(GrB_select(Diag, NULL, NULL, GrB_DIAG, S_BS, 0, NULL)); + // sum degrees per community + GRB_TRY (GrB_vxm (k, NULL, NULL, GxB_PLUS_FIRST_INT64, k, S, NULL)) ; + // square + GRB_TRY (GrB_assign (k, NULL, GrB_TIMES_INT64, k, GrB_ALL, n, NULL)) ; + + // sum all of the products + // this computed \sum_{i,j} (k_{i}k_{j}\delta(\sigma_i\sigma_j)) + GRB_TRY (GrB_reduce (&Q_, NULL, GrB_PLUS_MONOID_FP64, k, NULL)) ; + Q_ *= inv_m; + + // Count the number of edges per node in the cluster originating at a node + // in the cluster + GRB_TRY (GrB_mxm (B, S, NULL, GxB_PLUS_FIRST_FP64, A, S, GrB_DESC_S)); + + // sum the intra-cluster edges per node for all nodes + GRB_TRY (GrB_reduce (&Q_, GrB_PLUS_FP64, GrB_PLUS_MONOID_FP64, B, NULL)) ; - GRB_TRY(GrB_Matrix_reduce_FP64(&Q_, NULL, GrB_PLUS_MONOID_FP64, Diag, NULL)); Q_ *= -inv_m; *Q = Q_; From 7522000064ad194f1b612deab09fa32c4ae17b5a Mon Sep 17 00:00:00 2001 From: GomezGab Date: Fri, 7 Aug 2026 23:06:28 -0500 Subject: [PATCH 02/13] [WIP] Leiden implementation with GraphBLAS --- experimental/algorithm/LAGraph_Leiden.c | 290 ++++++++++++++++-------- experimental/test/test_leiden.c | 104 +++++++-- 2 files changed, 277 insertions(+), 117 deletions(-) diff --git a/experimental/algorithm/LAGraph_Leiden.c b/experimental/algorithm/LAGraph_Leiden.c index 56a644703a..683ac67138 100644 --- a/experimental/algorithm/LAGraph_Leiden.c +++ b/experimental/algorithm/LAGraph_Leiden.c @@ -11,11 +11,14 @@ // funding and support from the U.S. Government (see Acknowledgments.txt file). // DM22-0790 +#include "GraphBLAS.h" #include "LG_internal.h" #include #include +#include #include #include +#include #define LEIDEN_MAX_ITER 100 @@ -79,7 +82,7 @@ GrB_free (&S_mat) ; \ GrB_free (&A_temp) ; \ GrB_free (&one_scalar) ; \ - LAGR_LEIDEN_FREE_CONTAINER ; \ + GrB_free (&cont) ; \ LAGraph_Free ((void **) &k_arr, NULL) ; \ LAGraph_Free ((void **) &c_arr, NULL) ; \ LAGraph_Free ((void **) &k_comm, NULL) ; \ @@ -109,12 +112,185 @@ if (c_handle != NULL) GrB_free (c_handle) ; \ } -#if LG_SUITESPARSE_GRAPHBLAS_V10 -#define LAGR_LEIDEN_FREE_CONTAINER GxB_Container_free (&cont) +# if 1 +int LAGraph_Leiden +( + // output: + GrB_Vector *c_handle, // c[i] = community label (0..K-1) for node i + // input: + LAGraph_Graph G, // input graph (must be symmetric, no self-loops) + uint64_t seed, // random seed (reserved; not yet used) + char *msg +) +{ +#ifdef LG_SUITESPARSE_GRAPHBLAS_V10 + // queue: circular buffer with next nodes to look up + bool *enqueued = NULL; // true if node is in the queue + uint64_t *community = NULL ; // array: i in community c[i] + GrB_Matrix C = NULL ; // contains readonly pointer to community array + GrB_Matrix C_t = NULL ; // transpose of C + GrB_Matrix mask = NULL ; // mask each community + + uint64_t *sub_com = NULL ; // array: i in sub-community sub_com[i] + GrB_Matrix S = NULL ; // contains readonly pointer to sub-community array + + GrB_Vector deg = NULL ; // degree of each node + GrB_Vector c_deg = NULL ; // degree of each community + GrB_Vector c_list = NULL ; // active communities + + GrB_Matrix A = NULL ; // G->A + GrB_Vector x = NULL ; // used for BFS on the current node + GrB_Matrix X = NULL ; // used for BFS on the current node + GrB_Vector gain = NULL ; // gain from moving into a community + GrB_Vector k = NULL; + + GrB_Descriptor desc; + + uint64_t n; + + GxB_Iterator neighbor_it = NULL; + + A = G->A; + GrB_Matrix_nvals(&n, A); + // TODO: Asserts + // no self edges + // degree is cached + // undirected (although could be weighted) + + // TODO: initialize everything + + + double m = 0; + GrB_reduce (&m, NULL, GrB_PLUS_MONOID_FP64, A, NULL) ; + double m_inv2 = -1 / (2 * m); + // TODO: all of this should be wrapped in one more while loop + // TODO: make queue in random order + // TODO: load community array into C's j vector as readonly + + //-------------------------------------------------------------------------- + // Phase 1: move nodes + //-------------------------------------------------------------------------- + while (true) { // TODO: while queue not empty + uint64_t node_id, com_id; + + // TODO: pop node_id from the queue + enqueued [node_id] = false; + double com_deg, n_deg; + com_id = community[node_id] ; + GrB_Vector_setElement_FP64 (x, 0.0, node_id); + GrB_vxm (x, NULL, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, x, A, NULL) ; + // give gain the number of edges connecting x to community i + GrB_vxm (gain, NULL, NULL, GxB_PLUS_SECOND_FP64, x, C, NULL) ; + + // momentarily remove node degree from the community total + GrB_Vector_extractElement_FP64 (&n_deg, deg, node_id); + + // c_deg [com_id] -= ndeg + GrB_assign (c_deg, NULL, GrB_MINUS_FP64, n_deg, &com_id, 1, NULL); + + // Calculate gain in each neighboring community + GrB_apply (gain, gain, GrB_PLUS_FP64, GrB_TIMES_FP64, c_deg, m_inv2 * n_deg, GrB_DESC_S); + uint64_t max_gain_c; + // TODO: find which community has the max gain + + community [node_id] = max_gain_c; + // c_deg [node_id] -= ndeg + GrB_assign (c_deg, NULL, GrB_PLUS_FP64, n_deg, &max_gain_c, 1, NULL); + + if (max_gain_c == com_id) continue; // no change in community + + // push every neighbor in a different community that is not already in + // the queue to the back of the queue + GxB_Vector_Iterator_attach (neighbor_it, x, NULL); + GrB_Info info = GxB_Vector_Iterator_seek (neighbor_it, 0); + while (info == GrB_SUCCESS) { + uint64_t neighbor_id = GxB_Vector_Iterator_getIndex (neighbor_it); + if (max_gain_c != community [neighbor_id] + && !enqueued [neighbor_id]) { + // TODO: push to queue + } + + info = GxB_Vector_Iterator_next (neighbor_it); + } + } + + //-------------------------------------------------------------------------- + // Phase 2: refine + //-------------------------------------------------------------------------- + // TODO: resets + // TODO: check if transpose is even needed + + GrB_transpose (C_t, NULL, NULL, C, NULL) ; + // find number of in-cluster connections for each node in the cluster + // TODO: make sure C_t is 0s + // TODO: double check this + GrB_mxm (C_t, C_t, GrB_PLUS_FP64, GxB_PLUS_FIRST_FP64, A, C_t, GrB_DESC_S) ; + // TODO: discriminate the nodes with too few in-cluster connections + // (do not add them to the queue, and select them out of C_t) + GrB_mxv (c_list, NULL, NULL, GxB_ANY_PAIR_BOOL, C_t, x, NULL) ; + // TODO: interpret row list index descriptor + uint64_t n_communities; + GrB_Vector_nvals (&n_communities, c_list) ; + + GrB_Vector_new (&k, GrB_INT64, n_communities) ; + GrB_Matrix_new (&mask, GrB_INT64, n_communities, n) ; + GrB_Matrix_new (&X, GrB_INT64, n_communities, n) ; + + GxB_Matrix_extract_Vector (mask, NULL, NULL, C_t, c_list, NULL, desc); + + // NOTE: I am wiring this so that it will be relatively easy to parallelize + // this phase in the FUTURE + while (true) { // TODO: while queue not empty + uint64_t node_id, com_id; + + // TODO: pop node_id from the queue + enqueued [node_id] = false; + double com_deg, n_deg; + com_id = community[node_id] ; + uint64_t sub_com_id = sub_com [node_id] ; + + GrB_Matrix_setElement_FP64 (X, 0.0, com_id, node_id); + GrB_mxm (X, mask, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, X, A, GrB_DESC_S) ; + // give X the number of edges connecting x to community i + GrB_mxm (X, X, NULL, GxB_PLUS_SECOND_FP64, X, S, NULL) ; + + // momentarily remove node degree from the community total + GrB_Vector_extractElement_FP64 (&n_deg, deg, node_id); + + // c_deg [com_id] -= ndeg + GrB_assign (c_deg, NULL, GrB_MINUS_FP64, n_deg, &com_id, 1, NULL) ; + + GrB_Vector_setElement_FP64(k, m_inv2 * n_deg, com_id) ; + // Calculate gain in each neighboring community + GrB_mxm (X, X, GrB_PLUS_FP64, GrB_PLUS_TIMES_SEMIRING_FP64, (GrB_Matrix) c_deg, (GrB_Matrix) k, GrB_DESC_ST1) ; + + uint64_t new_c; + // TODO: find which community to go into via random selection + + sub_com [node_id] = new_c; + // c_deg [node_id] -= ndeg + GrB_assign (s_deg, NULL, GrB_PLUS_FP64, n_deg, &new_c, 1, NULL); + } + + //-------------------------------------------------------------------------- + // Phase 3: aggregate + //-------------------------------------------------------------------------- + + GrB_mxv (c_list, NULL, NULL, GxB_ANY_PAIR_BOOL, S_t, x, NULL) ; + // FIXME: this modifies the input matrix + GrB_mxm (A, NULL, NULL, GxB_PLUS_SECOND_FP64, S_t, A, NULL) ; + GrB_mxm (A, NULL, NULL, GxB_PLUS_FIRST_FP64, A, S, NULL) ; + GrB_select (A, NULL, NULL, GrB_OFFDIAG, A, 0, NULL) ; + GxB_Matrix_extract_Vector (A, NULL, NULL, A, c_list, c_list, desc) ; + // TODO: carry over the community vector and matrix and start from the top + + + return GrB_NOT_IMPLEMENTED; #else -#define LAGR_LEIDEN_FREE_CONTAINER ((void) 0) + return GrB_NOT_IMPLEMENTED; #endif - +} +#else int LAGraph_Leiden ( // output: @@ -137,9 +313,7 @@ int LAGraph_Leiden GrB_Matrix S_mat = NULL ; // temporary membership matrix (Phase 3) GrB_Matrix A_temp = NULL ; // temporary for mxm (Phase 3) GrB_Scalar one_scalar = NULL ; // FP64 scalar with value 1.0 for build_Scalar -#if LG_SUITESPARSE_GRAPHBLAS_V10 GxB_Container cont = NULL ; // for unloading A_cur into raw CSR arrays -#endif double *k_arr = NULL ; // k_arr[i] = degree of node i (current level) int64_t *c_arr = NULL ; // c_arr[i] = Phase-1 community label double *k_comm = NULL ; // k_comm[l] = total degree of community l @@ -200,6 +374,10 @@ int LAGraph_Leiden // allocate workspace (all arrays sized n; used for indices 0..n_cur-1) //-------------------------------------------------------------------------- + // void *workspace = NULL, p_workspace = NULL; + // LG_TRY (LAGraph_Malloc (&workspace, n, sizeof (double) * 4 + + // sizeof (int64_t) * 4 + sizeof (int8_t) + sizeof (GrB_Index) *4, msg)) ; + LG_TRY (LAGraph_Malloc ((void **) &k_arr, n, sizeof (double), msg)) ; LG_TRY (LAGraph_Malloc ((void **) &c_arr, n, sizeof (int64_t), msg)) ; LG_TRY (LAGraph_Malloc ((void **) &k_comm, n, sizeof (double), msg)) ; @@ -213,7 +391,6 @@ int LAGraph_Leiden LG_TRY (LAGraph_Malloc ((void **) &o_comm, n, sizeof (int64_t), msg)) ; LG_TRY (LAGraph_Malloc ((void **) &init_comm, n, sizeof (GrB_Index), msg)) ; LG_TRY (LAGraph_Malloc ((void **) &iota, n, sizeof (GrB_Index), msg)) ; - for (GrB_Index i = 0 ; i < n ; i++) iota[i] = i ; // Reusable FP64 scalar with value 1.0 for GxB_Matrix_build_Scalar. GRB_TRY (GrB_Scalar_new (&one_scalar, GrB_FP64)) ; @@ -221,9 +398,7 @@ int LAGraph_Leiden // Reusable container for unloading A_cur into raw CSR arrays per level // (only available with the SuiteSparse v10+ Container API). -#if LG_SUITESPARSE_GRAPHBLAS_V10 GRB_TRY (GxB_Container_new (&cont)) ; -#endif //-------------------------------------------------------------------------- // compute m = total edge weight / 2 from G->A (invariant under aggregation) @@ -233,15 +408,16 @@ int LAGraph_Leiden GRB_TRY (GrB_Matrix_reduce_FP64 (&m, NULL, GrB_PLUS_MONOID_FP64, A, NULL)) ; m /= 2.0 ; - double two_m = 2.0 * m ; // denominator of the modularity penalty // Empty graph: return a singleton partition. if (m == 0.0) { // c[i] = i for all i, built in one call instead of n setElement calls. GRB_TRY (GrB_Vector_new (c_handle, GrB_INT64, n)) ; - GRB_TRY (GrB_Vector_build_INT64 (*c_handle, iota, iota, n, - GrB_FIRST_INT64)) ; + GRB_TRY (GrB_Vector_assign_INT64 ( + *c_handle, NULL, NULL, 0, GrB_ALL, 0, NULL)) ; + GRB_TRY (GrB_Vector_apply_IndexOp_INT64 ( + *c_handle, NULL, NULL, GrB_ROWINDEX_INT64, *c_handle, 0, NULL)) ; LG_FREE_WORK ; return (GrB_SUCCESS) ; } @@ -257,6 +433,7 @@ int LAGraph_Leiden for (GrB_Index i = 0 ; i < n ; i++) { + iota[i] = i ; o_comm[i] = (int64_t) i ; init_comm[i] = i ; } @@ -301,7 +478,6 @@ int LAGraph_Leiden GRB_TRY (GrB_Matrix_reduce_Monoid (k_vec, NULL, GrB_PLUS_FP64, GrB_PLUS_MONOID_FP64, A_cur, NULL)) ; -#if LG_SUITESPARSE_GRAPHBLAS_V10 // Unload k_vec into k_arr (dense FP64 array of length n_cur). // The previous k_arr workspace allocation is replaced by a pointer // owned by GraphBLAS until LAGraph_Free reclaims it. @@ -366,75 +542,6 @@ int LAGraph_Leiden Ap = (GrB_Index *) pv ; Aj = (GrB_Index *) iv ; Ax = (double *) xv ; -#else - // Fallback for SuiteSparse:GraphBLAS < v10.0.0 (no Container API): - // copy degrees out of k_vec, then materialize CSR via extractTuples - // + counting-sort scatter. k_vec contains entries only for non-zero - // rows, so zero-fill k_arr first then scatter. - for (GrB_Index i = 0 ; i < n_cur ; i++) k_arr[i] = 0.0 ; - { - GrB_Index nvk = n_cur ; - GrB_Index *Ik = NULL ; - double *Xk = NULL ; - LG_TRY (LAGraph_Malloc ((void **) &Ik, n_cur, sizeof (GrB_Index), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &Xk, n_cur, sizeof (double), msg)) ; - GRB_TRY (GrB_Vector_extractTuples_FP64 (Ik, Xk, &nvk, k_vec)) ; - for (GrB_Index t = 0 ; t < nvk ; t++) k_arr[Ik[t]] = Xk[t] ; - LAGraph_Free ((void **) &Ik, NULL) ; - LAGraph_Free ((void **) &Xk, NULL) ; - } - - GrB_Index Anz ; - GRB_TRY (GrB_Matrix_nvals (&Anz, A_cur)) ; - if (Ap_cap < n_cur + 1) - { - LAGraph_Free ((void **) &Ap, NULL) ; - LAGraph_Free ((void **) &cursor, NULL) ; - LG_TRY (LAGraph_Malloc ((void **) &Ap, n_cur + 1, - sizeof (GrB_Index), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &cursor, n_cur, - sizeof (GrB_Index), msg)) ; - Ap_cap = n_cur + 1 ; - } - if (Anz_cap < Anz) - { - GrB_Index newcap = (Anz < 16) ? 16 : Anz ; - LAGraph_Free ((void **) &Aj, NULL) ; - LAGraph_Free ((void **) &Ax, NULL) ; - LAGraph_Free ((void **) &I_tup, NULL) ; - LAGraph_Free ((void **) &J_tup, NULL) ; - LAGraph_Free ((void **) &X_tup, NULL) ; - LG_TRY (LAGraph_Malloc ((void **) &Aj, newcap, - sizeof (GrB_Index), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &Ax, newcap, - sizeof (double), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &I_tup, newcap, - sizeof (GrB_Index), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &J_tup, newcap, - sizeof (GrB_Index), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &X_tup, newcap, - sizeof (double), msg)) ; - Anz_cap = newcap ; - } - - memset (Ap, 0, (n_cur + 1) * sizeof (GrB_Index)) ; - if (Anz > 0) - { - GrB_Index nout = Anz ; - GRB_TRY (GrB_Matrix_extractTuples_FP64 (I_tup, J_tup, X_tup, - &nout, A_cur)) ; - for (GrB_Index t = 0 ; t < Anz ; t++) Ap[I_tup[t] + 1]++ ; - for (GrB_Index r = 0 ; r < n_cur ; r++) Ap[r + 1] += Ap[r] ; - memcpy (cursor, Ap, n_cur * sizeof (GrB_Index)) ; - for (GrB_Index t = 0 ; t < Anz ; t++) - { - GrB_Index r = I_tup[t] ; - GrB_Index dst = cursor[r]++ ; - Aj[dst] = J_tup[t] ; - Ax[dst] = X_tup[t] ; - } - } -#endif //---------------------------------------------------------------------- // PHASE 1: Local Move Phase @@ -487,7 +594,7 @@ int LAGraph_Leiden } double T_ci = dirty[ci] ? T_local[ci] : 0.0 ; - double score_ci = T_ci - ki * k_comm[ci] / two_m ; + double score_ci = T_ci - ki * k_comm[ci] / (m * 2) ; double best_score = score_ci ; int64_t best_c = ci ; @@ -495,7 +602,7 @@ int LAGraph_Leiden { int64_t c_cand = (int64_t) dirty_list[d] ; if (c_cand == ci) continue ; - double score = T_local[c_cand] - ki * k_comm[c_cand] / two_m ; + double score = T_local[c_cand] - ki * k_comm[c_cand] / (m * 2) ; if (score > best_score) { best_score = score ; @@ -572,7 +679,7 @@ int LAGraph_Leiden } double T_ci_ref = dirty[ci_ref] ? T_local[ci_ref] : 0.0 ; - double score_ci_ref = T_ci_ref - ki * k_ref_comm[ci_ref] / two_m ; + double score_ci_ref = T_ci_ref - ki * k_ref_comm[ci_ref] / (m * 2) ; double best_score = score_ci_ref ; int64_t best_c_ref = ci_ref ; @@ -580,7 +687,7 @@ int LAGraph_Leiden { int64_t c_cand = (int64_t) dirty_list[d] ; if (c_cand == ci_ref) continue ; - double score = T_local[c_cand] - ki * k_ref_comm[c_cand] / two_m ; + double score = T_local[c_cand] - ki * k_ref_comm[c_cand] / (m * 2) ; if (score > best_score) { best_score = score ; @@ -656,7 +763,6 @@ int LAGraph_Leiden // No-op on the v9 fallback (A_cur was never unloaded). //---------------------------------------------------------------------- -#if LG_SUITESPARSE_GRAPHBLAS_V10 GRB_TRY (GxB_Vector_load (cont->p, (void **) &Ap, pty, pn, psz, ph, NULL)) ; Ap = NULL ; @@ -667,7 +773,6 @@ int LAGraph_Leiden xn, xsz, xh, NULL)) ; Ax = NULL ; GRB_TRY (GxB_load_Matrix_from_Container (A_cur, cont, NULL)) ; -#endif //---------------------------------------------------------------------- // PHASE 3: Aggregation — build coarsened graph if communities merged @@ -719,16 +824,13 @@ int LAGraph_Leiden //-------------------------------------------------------------------------- GRB_TRY (GrB_Vector_new (c_handle, GrB_INT64, n)) ; -#if LG_SUITESPARSE_GRAPHBLAS_V10 GRB_TRY (GrB_set (*c_handle, GxB_FULL, GxB_SPARSITY_CONTROL)) ; GRB_TRY (GxB_Vector_load (*c_handle, (void **) &o_comm, GrB_INT64, n, n * sizeof (int64_t), GrB_DEFAULT, NULL)) ; o_comm = NULL ; // ownership transferred to *c_handle -#else - GRB_TRY (GrB_Vector_build_INT64 (*c_handle, iota, o_comm, n, - GrB_FIRST_INT64)) ; -#endif LG_FREE_WORK ; return (GrB_SUCCESS) ; } +#endif + diff --git a/experimental/test/test_leiden.c b/experimental/test/test_leiden.c index 3c9d1a65e9..b2dd3c3728 100644 --- a/experimental/test/test_leiden.c +++ b/experimental/test/test_leiden.c @@ -8,6 +8,7 @@ #include #include +#include "GraphBLAS.h" #include "LG_Xtest.h" #include #include @@ -22,12 +23,17 @@ typedef struct { const char *matrix_file ; bool expect_communities ; // true iff non-trivial Q is expected + double min_modularity ; // minimum acceptable modularity + double min_coverage ; // minimum intra-cluster edge ratio + double min_performance ; // minimum partition performance } matrix_info ; const matrix_info files[] = { - { "karate.mtx", true }, // Zachary karate club (Q > 0 expected) - { "", false } + // matrix_file, expect_communities, min_Q, min_coverage, min_performance + { "karate.mtx", true, 0.35, 0.15, 0.30 }, + { "comm0.mtx", true, 0.25, 0.10, 0.20 }, + { "", false, -1.0, -1.0, -1.0 } } ; //------------------------------------------------------------------------------ @@ -38,12 +44,6 @@ void test_Leiden (void) { LAGraph_Init (msg) ; -#if LAGRAPH_SUITESPARSE - // Disable JIT before any GraphBLAS operations; the JIT compiler may not - // be available in all build environments. - OK (LG_SET_JIT (LG_JIT_OFF)) ; -#endif - for (int k = 0 ;; k++) { const char *aname = files[k].matrix_file ; @@ -65,7 +65,7 @@ void test_Leiden (void) // Ensure symmetry cache is populated (required by some checks). OK (LAGraph_Cached_IsSymmetricStructure (G, msg)) ; - uint64_t seed = 42 ; + uint64_t seed = 0 ; //unused GrB_Vector c = NULL ; OK (LAGraph_Leiden (&c, G, seed, msg)) ; @@ -80,28 +80,86 @@ void test_Leiden (void) (unsigned long long) n, (unsigned long long) nvals) ; // Community labels must be in [0, n-1]. - for (GrB_Index i = 0 ; i < n ; i++) - { - int64_t label = 0 ; - OK (GrB_Vector_extractElement_INT64 (&label, c, i)) ; - TEST_CHECK (label >= 0 && (GrB_Index) label < n) ; - } + int64_t min_label = 0, max_label = 0 ; + + OK (GrB_Vector_reduce_INT64 ( + &min_label, NULL, GxB_MIN_INT64_MONOID, c, NULL)) ; + OK (GrB_Vector_reduce_INT64 ( + &max_label, NULL, GxB_MAX_INT64_MONOID, c, NULL)) ; + + TEST_CHECK (min_label >= 0) ; + TEST_CHECK (max_label < n) ; + + GrB_Index n_communities = max_label + 1 ; -#if LAGRAPH_SUITESPARSE // Compute modularity Q (requires SuiteSparse:GraphBLAS). double Q = 0.0 ; OK (LAGr_Modularity (&Q, 1.0, c, G, msg)) ; - printf (" Modularity Q = %f (K_ref <= %llu)\n", - Q, (unsigned long long) n) ; + printf (" Modularity Q = %f\n", Q) ; if (files[k].expect_communities) { - // Multi-level Leiden on karate.mtx achieves Q ≈ 0.42, matching - // published Louvain/Leiden benchmarks on this graph. - TEST_CHECK (Q > 0.37) ; - TEST_MSG ("Expected Q > 0.37 for %s, got Q = %f", aname, Q) ; + // Validate modularity meets threshold for this graph + TEST_CHECK (Q > files[k].min_modularity) ; + TEST_MSG ("Expected Q > %f for %s, got Q = %f", + files[k].min_modularity, aname, Q) ; } -#endif + + // Compute per-community statistics for robustness validation + printf (" Community statistics:\n") ; + GrB_Index total_community_edges = 0 ; + GrB_Index total_edges = 0 ; + OK (GrB_Matrix_nvals (&total_edges, G->A)) ; + + for (GrB_Index com = 0 ; com < n_communities ; com++) + { + GrB_Index com_size = 0 ; + GrB_Index com_edges = 0 ; + + // Count nodes and edges in this community + for (GrB_Index i = 0 ; i < n ; i++) + { + int64_t label ; + GrB_Info info = GrB_Vector_extractElement_INT64 (&label, c, i) ; + if (info == GrB_SUCCESS && label == (int64_t)com) + { + com_size++ ; + // Count internal edges: edges from i to other nodes in same community + for (GrB_Index j = 0 ; j < n ; j++) + { + int64_t label_j ; + GrB_Info info_j = GrB_Vector_extractElement_INT64 (&label_j, c, j) ; + if (info_j == GrB_SUCCESS && label_j == (int64_t)com) + { + bool has_edge = false ; + GrB_Info info_e = GrB_Matrix_extractElement_BOOL (&has_edge, G->A, i, j) ; + if (info_e == GrB_SUCCESS && has_edge) + { + com_edges++ ; + } + } + } + } + } + + if (com_size > 0) + { + double density = (com_size > 1) ? + (double)com_edges / (com_size * (com_size - 1)) : 0.0 ; + printf (" Community %llu: size=%llu, edges=%llu, density=%f\n", + (unsigned long long) com, (unsigned long long) com_size, + (unsigned long long) com_edges, density) ; + total_community_edges += com_edges ; + } + } + + // Compute edge-cut ratio (edges crossing communities / total edges) + GrB_Index cut_edges = total_edges - total_community_edges ; + double cut_ratio = (total_edges > 0) ? + (double)cut_edges / total_edges : 0.0 ; + printf (" Edge-cut ratio = %f (cut=%llu, total=%llu)\n", + cut_ratio, (unsigned long long) cut_edges, + (unsigned long long) total_edges) ; GrB_free (&c) ; OK (LAGraph_Delete (&G, msg)) ; From 124c836919babc49a49ab7b1402a15843de1745d Mon Sep 17 00:00:00 2001 From: GomezGab Date: Mon, 10 Aug 2026 18:25:21 -0500 Subject: [PATCH 03/13] [WIP] Break up Leiden phases into helper functions --- experimental/algorithm/LAGraph_Leiden.c | 1044 +++++++++-------------- 1 file changed, 397 insertions(+), 647 deletions(-) diff --git a/experimental/algorithm/LAGraph_Leiden.c b/experimental/algorithm/LAGraph_Leiden.c index 683ac67138..afae8284d6 100644 --- a/experimental/algorithm/LAGraph_Leiden.c +++ b/experimental/algorithm/LAGraph_Leiden.c @@ -19,6 +19,7 @@ #include #include #include +#include #define LEIDEN_MAX_ITER 100 @@ -76,221 +77,317 @@ #undef LG_FREE_WORK #define LG_FREE_WORK \ { \ - GrB_free (&k_vec) ; \ - GrB_free (&A_agg) ; \ - GrB_free (&A_new) ; \ - GrB_free (&S_mat) ; \ - GrB_free (&A_temp) ; \ - GrB_free (&one_scalar) ; \ - GrB_free (&cont) ; \ - LAGraph_Free ((void **) &k_arr, NULL) ; \ - LAGraph_Free ((void **) &c_arr, NULL) ; \ - LAGraph_Free ((void **) &k_comm, NULL) ; \ - LAGraph_Free ((void **) &c_p1, NULL) ; \ - LAGraph_Free ((void **) &c_ref, NULL) ; \ - LAGraph_Free ((void **) &k_ref_comm, NULL) ; \ - LAGraph_Free ((void **) &T_local, NULL) ; \ - LAGraph_Free ((void **) &dirty, NULL) ; \ - LAGraph_Free ((void **) &dirty_list, NULL) ; \ - LAGraph_Free ((void **) &remap, NULL) ; \ - LAGraph_Free ((void **) &o_comm, NULL) ; \ - LAGraph_Free ((void **) &init_comm, NULL) ; \ - LAGraph_Free ((void **) &Ap, NULL) ; \ - LAGraph_Free ((void **) &Aj, NULL) ; \ - LAGraph_Free ((void **) &Ax, NULL) ; \ - LAGraph_Free ((void **) &I_tup, NULL) ; \ - LAGraph_Free ((void **) &J_tup, NULL) ; \ - LAGraph_Free ((void **) &X_tup, NULL) ; \ - LAGraph_Free ((void **) &cursor, NULL) ; \ - LAGraph_Free ((void **) &iota, NULL) ; \ + GrB_free (&it); \ + GrB_free (&c_deg); \ + GrB_free (&x); \ + GrB_free (&gain); \ } #undef LG_FREE_ALL #define LG_FREE_ALL \ { \ LG_FREE_WORK ; \ - if (c_handle != NULL) GrB_free (c_handle) ; \ } -# if 1 -int LAGraph_Leiden +// helper function: phase 1 of leiden. Output already allocated +int LG_Leiden_move_nodes ( // output: - GrB_Vector *c_handle, // c[i] = community label (0..K-1) for node i + GrB_Matrix C, // community matrix + uint64_t *community, // community array (inside of commmunity matrix) // input: - LAGraph_Graph G, // input graph (must be symmetric, no self-loops) - uint64_t seed, // random seed (reserved; not yet used) - char *msg -) -{ -#ifdef LG_SUITESPARSE_GRAPHBLAS_V10 - // queue: circular buffer with next nodes to look up - bool *enqueued = NULL; // true if node is in the queue - uint64_t *community = NULL ; // array: i in community c[i] - GrB_Matrix C = NULL ; // contains readonly pointer to community array - GrB_Matrix C_t = NULL ; // transpose of C - GrB_Matrix mask = NULL ; // mask each community - - uint64_t *sub_com = NULL ; // array: i in sub-community sub_com[i] - GrB_Matrix S = NULL ; // contains readonly pointer to sub-community array - - GrB_Vector deg = NULL ; // degree of each node - GrB_Vector c_deg = NULL ; // degree of each community - GrB_Vector c_list = NULL ; // active communities - - GrB_Matrix A = NULL ; // G->A - GrB_Vector x = NULL ; // used for BFS on the current node - GrB_Matrix X = NULL ; // used for BFS on the current node - GrB_Vector gain = NULL ; // gain from moving into a community - GrB_Vector k = NULL; - - GrB_Descriptor desc; - - uint64_t n; - - GxB_Iterator neighbor_it = NULL; - - A = G->A; - GrB_Matrix_nvals(&n, A); - // TODO: Asserts - // no self edges - // degree is cached - // undirected (although could be weighted) - - // TODO: initialize everything - - - double m = 0; - GrB_reduce (&m, NULL, GrB_PLUS_MONOID_FP64, A, NULL) ; - double m_inv2 = -1 / (2 * m); - // TODO: all of this should be wrapped in one more while loop - // TODO: make queue in random order - // TODO: load community array into C's j vector as readonly - - //-------------------------------------------------------------------------- - // Phase 1: move nodes - //-------------------------------------------------------------------------- - while (true) { // TODO: while queue not empty - uint64_t node_id, com_id; - - // TODO: pop node_id from the queue + const GrB_Matrix A, // adjacency matrix + const GrB_Vector deg, // degree vector + uint64_t *queue, // queue to use (contains all nodes) + bool *enqueued, // nodes in queue (all at the start) + uint64_t m_inv2, // 1 / (2 * m) + char* msg +) { + uint64_t node_id, com_id, n_deg, queue_head = 0, queue_tail, n; + GrB_Vector x = NULL, gain = NULL; + GrB_Vector c_deg = NULL ; + GxB_Iterator it = NULL; + GrB_Matrix_nrows (&n, A) ; + queue_tail = n; + + + + GRB_TRY (GxB_Iterator_new (&it)) ; + GRB_TRY (GrB_Vector_new (&c_deg, GrB_FP64, n)) ; + GRB_TRY (GrB_Vector_new (&x, GrB_FP64, n)) ; + GRB_TRY (GrB_Vector_new (&gain, GrB_FP64, n)) ; + + GRB_TRY (GrB_mxv(c_deg, NULL, NULL, GxB_PLUS_SECOND_FP64, C, deg, NULL)); + // TODO: decide if queue should be made inside this function. + while (queue_head != queue_tail) { // while queue not empty + node_id = queue[queue_head]; + queue_head = (queue_head + 1) % n; enqueued [node_id] = false; double com_deg, n_deg; com_id = community[node_id] ; - GrB_Vector_setElement_FP64 (x, 0.0, node_id); - GrB_vxm (x, NULL, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, x, A, NULL) ; + GRB_TRY (GrB_Vector_setElement_FP64 (x, 0.0, node_id)); + GRB_TRY (GrB_vxm (x, NULL, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, x, A, NULL)) ; // give gain the number of edges connecting x to community i - GrB_vxm (gain, NULL, NULL, GxB_PLUS_SECOND_FP64, x, C, NULL) ; + GRB_TRY (GrB_vxm (gain, NULL, NULL, GxB_PLUS_SECOND_FP64, x, C, NULL)) ; // momentarily remove node degree from the community total - GrB_Vector_extractElement_FP64 (&n_deg, deg, node_id); + GRB_TRY (GrB_Vector_extractElement_FP64 (&n_deg, deg, node_id)); // c_deg [com_id] -= ndeg - GrB_assign (c_deg, NULL, GrB_MINUS_FP64, n_deg, &com_id, 1, NULL); + GRB_TRY (GrB_assign (c_deg, NULL, GrB_MINUS_FP64, n_deg, &com_id, 1, NULL)); // Calculate gain in each neighboring community - GrB_apply (gain, gain, GrB_PLUS_FP64, GrB_TIMES_FP64, c_deg, m_inv2 * n_deg, GrB_DESC_S); + GRB_TRY (GrB_apply (gain, gain, GrB_PLUS_FP64, GrB_TIMES_FP64, c_deg, m_inv2 * n_deg, GrB_DESC_S)); uint64_t max_gain_c; - // TODO: find which community has the max gain + double max_gain_val = -1e100; + for (uint64_t c = 0; c < n; c++) { + double g; + GRB_TRY (GrB_Vector_extractElement_FP64(&g, gain, c)); + if (g > max_gain_val) { + max_gain_val = g; + max_gain_c = c; + } + } community [node_id] = max_gain_c; // c_deg [node_id] -= ndeg - GrB_assign (c_deg, NULL, GrB_PLUS_FP64, n_deg, &max_gain_c, 1, NULL); + GRB_TRY (GrB_assign (c_deg, NULL, GrB_PLUS_FP64, n_deg, &max_gain_c, 1, NULL)); if (max_gain_c == com_id) continue; // no change in community // push every neighbor in a different community that is not already in // the queue to the back of the queue - GxB_Vector_Iterator_attach (neighbor_it, x, NULL); - GrB_Info info = GxB_Vector_Iterator_seek (neighbor_it, 0); + GRB_TRY (GxB_Vector_Iterator_attach (it, x, NULL)); + GrB_Info info = GxB_Vector_Iterator_seek (it, 0); while (info == GrB_SUCCESS) { - uint64_t neighbor_id = GxB_Vector_Iterator_getIndex (neighbor_it); + uint64_t neighbor_id = GxB_Vector_Iterator_getIndex (it); if (max_gain_c != community [neighbor_id] && !enqueued [neighbor_id]) { - // TODO: push to queue + queue[queue_tail] = neighbor_id; + queue_tail = (queue_tail + 1) % n; + enqueued[neighbor_id] = true; } - info = GxB_Vector_Iterator_next (neighbor_it); + info = GxB_Vector_Iterator_next (it); } } + LG_FREE_WORK ; + return GrB_SUCCESS ; +} + +#undef LG_FREE_WORK +#define LG_FREE_WORK \ +{ \ + GrB_free (&s_deg); \ + GrB_free (&c_list); \ + GrB_free (&x_count); \ + GrB_free (&k); \ + GrB_free (&X); \ + GrB_free (&C_t); \ + GrB_free (&mask); \ +} + +#undef LG_FREE_ALL +#define LG_FREE_ALL \ +{ \ + LG_FREE_WORK ; \ +} + +// helper function: phase 2 of leiden. Output already allocated +int LG_Leiden_refinement +( + // output: + GrB_Matrix S, // sub-community matrix + uint64_t *sub_com, // sub-community array (inside of S matrix) + // input: + const GrB_Matrix C, // community matrix + const uint64_t *community, // community array (inside of C matrix) + const GrB_Vector c_deg, // degree vector + const GrB_Matrix A, // adjacency matrix + const GrB_Vector deg, // degree vector + uint64_t *queue, // queue to use (contains all nodes) + uint64_t m_inv2, // 1 / (2 * m) + char* msg +) { + uint64_t node_id, com_id, n_deg, queue_head = 0, queue_tail, n; + GrB_Matrix X = NULL, mask = NULL; + GrB_Vector s_deg = NULL; + GrB_Vector c_list = NULL; + GrB_Vector k = NULL; + GrB_Matrix C_t = NULL ; // transpose of C + GrB_Vector x_count = NULL; + + GRB_TRY (GrB_Matrix_nrows (&n, A)) ; + queue_tail = n; + + // TODO: initialize x. decide if queue should be made inside this function. + GRB_TRY (GrB_Matrix_new (&C_t, GrB_FP64, n, n)) ; + GRB_TRY (GrB_transpose (C_t, NULL, NULL, C, NULL)) ; + + // Find number of in-cluster connections for each node in the cluster + GRB_TRY (GrB_mxm (C_t, C_t, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, C_t, A, GrB_DESC_ST1)) ; + + // Identify nodes with in-cluster connections (more than 0) + // TODO: use a smarter selection using c_deg + GRB_TRY (GrB_Vector_new (&x_count, GrB_FP64, n)); + GRB_TRY (GrB_assign (x_count, NULL, NULL, 0.0, GrB_ALL, n, NULL)) ; + GRB_TRY (GrB_reduce (x_count, NULL, GrB_PLUS_FP64, GrB_PLUS_MONOID_FP64, C_t, GrB_DESC_T0)) ; + GRB_TRY (GrB_mxv (c_list, NULL, NULL, GxB_ANY_PAIR_BOOL, C_t, deg, NULL)) ; - //-------------------------------------------------------------------------- - // Phase 2: refine - //-------------------------------------------------------------------------- - // TODO: resets - // TODO: check if transpose is even needed - - GrB_transpose (C_t, NULL, NULL, C, NULL) ; - // find number of in-cluster connections for each node in the cluster - // TODO: make sure C_t is 0s - // TODO: double check this - GrB_mxm (C_t, C_t, GrB_PLUS_FP64, GxB_PLUS_FIRST_FP64, A, C_t, GrB_DESC_S) ; - // TODO: discriminate the nodes with too few in-cluster connections - // (do not add them to the queue, and select them out of C_t) - GrB_mxv (c_list, NULL, NULL, GxB_ANY_PAIR_BOOL, C_t, x, NULL) ; - // TODO: interpret row list index descriptor uint64_t n_communities; - GrB_Vector_nvals (&n_communities, c_list) ; - GrB_Vector_new (&k, GrB_INT64, n_communities) ; - GrB_Matrix_new (&mask, GrB_INT64, n_communities, n) ; - GrB_Matrix_new (&X, GrB_INT64, n_communities, n) ; + GRB_TRY (GrB_Vector_nvals (&n_communities, c_list)) ; + GRB_TRY (GrB_Vector_new (&k, GrB_FP64, n_communities)) ; + GRB_TRY (GrB_Matrix_new (&mask, GrB_FP64, n_communities, n)) ; + GRB_TRY (GrB_Matrix_new (&X, GrB_FP64, n_communities, n)) ; - GxB_Matrix_extract_Vector (mask, NULL, NULL, C_t, c_list, NULL, desc); + GRB_TRY (GxB_Matrix_extract_Vector (mask, NULL, NULL, C_t, c_list, NULL, NULL)) ; + GRB_TRY (GrB_select (mask, NULL, NULL, GrB_VALUEGT_FP64, mask, 0.0, NULL)) ; // NOTE: I am wiring this so that it will be relatively easy to parallelize // this phase in the FUTURE - while (true) { // TODO: while queue not empty - uint64_t node_id, com_id; + while (queue_head != queue_tail) { // while queue not empty + double count; + node_id = queue[queue_head]; + GRB_TRY (GrB_Vector_extractElement_FP64 (&count, x_count, node_id)); + if (count <= 0.0) continue; // skip singletons - // TODO: pop node_id from the queue - enqueued [node_id] = false; + queue_head = (queue_head + 1) % n; double com_deg, n_deg; - com_id = community[node_id] ; + com_id = community [node_id] ; uint64_t sub_com_id = sub_com [node_id] ; - GrB_Matrix_setElement_FP64 (X, 0.0, com_id, node_id); - GrB_mxm (X, mask, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, X, A, GrB_DESC_S) ; + GRB_TRY (GrB_Matrix_setElement_FP64 (X, 0.0, com_id, node_id)); + GRB_TRY (GrB_mxm (X, mask, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, X, A, GrB_DESC_S)) ; // give X the number of edges connecting x to community i - GrB_mxm (X, X, NULL, GxB_PLUS_SECOND_FP64, X, S, NULL) ; + GRB_TRY (GrB_mxm (X, NULL, NULL, GxB_PLUS_SECOND_FP64, X, S, NULL)) ; // momentarily remove node degree from the community total - GrB_Vector_extractElement_FP64 (&n_deg, deg, node_id); + GRB_TRY (GrB_Vector_extractElement_FP64 (&n_deg, deg, node_id)); - // c_deg [com_id] -= ndeg - GrB_assign (c_deg, NULL, GrB_MINUS_FP64, n_deg, &com_id, 1, NULL) ; + // s_deg [sub_com_id] -= ndeg + GRB_TRY (GrB_assign (s_deg, NULL, GrB_MINUS_FP64, n_deg, &sub_com_id, 1, NULL)) ; - GrB_Vector_setElement_FP64(k, m_inv2 * n_deg, com_id) ; + GRB_TRY (GrB_Vector_setElement_FP64(k, m_inv2 * n_deg, sub_com_id)) ; // Calculate gain in each neighboring community - GrB_mxm (X, X, GrB_PLUS_FP64, GrB_PLUS_TIMES_SEMIRING_FP64, (GrB_Matrix) c_deg, (GrB_Matrix) k, GrB_DESC_ST1) ; + GRB_TRY (GrB_mxm (X, X, GrB_PLUS_FP64, GrB_PLUS_TIMES_SEMIRING_FP64, + (GrB_Matrix) s_deg, (GrB_Matrix) k, GrB_DESC_ST1)) ; uint64_t new_c; - // TODO: find which community to go into via random selection + // TODO: select from X (those with gain >= current and choose randomly) sub_com [node_id] = new_c; - // c_deg [node_id] -= ndeg - GrB_assign (s_deg, NULL, GrB_PLUS_FP64, n_deg, &new_c, 1, NULL); + // s_deg [new_c] += ndeg + GRB_TRY (GrB_assign (s_deg, NULL, GrB_PLUS_FP64, n_deg, &new_c, 1, NULL)); } + LG_FREE_WORK ; + return GrB_SUCCESS ; +} - //-------------------------------------------------------------------------- - // Phase 3: aggregate - //-------------------------------------------------------------------------- +#undef LG_FREE_WORK +#define LG_FREE_WORK \ +{ \ + GrB_free (&desc) ; \ + GrB_free (&s_list) ; \ + GrB_free (&S_squished) ; \ + GrB_free (&C_squished) ; \ + GrB_free (&A_squished) ; \ +} - GrB_mxv (c_list, NULL, NULL, GxB_ANY_PAIR_BOOL, S_t, x, NULL) ; - // FIXME: this modifies the input matrix - GrB_mxm (A, NULL, NULL, GxB_PLUS_SECOND_FP64, S_t, A, NULL) ; - GrB_mxm (A, NULL, NULL, GxB_PLUS_FIRST_FP64, A, S, NULL) ; - GrB_select (A, NULL, NULL, GrB_OFFDIAG, A, 0, NULL) ; - GxB_Matrix_extract_Vector (A, NULL, NULL, A, c_list, c_list, desc) ; - // TODO: carry over the community vector and matrix and start from the top +#undef LG_FREE_ALL +#define LG_FREE_ALL \ +{ \ + GrB_free (S) ; \ + GrB_free (C) ; \ + if (free_A) GrB_free (A) ; \ + LG_FREE_WORK ; \ +} +// helper function: phase 3 of leiden. input freed and output allocated +int LG_Leiden_aggregate +( + // input / output: + GrB_Matrix *S, // sub-communities to agregate + GrB_Matrix *C, // communities to preserve + GrB_Matrix *A, // adjacency matrix + bool free_A, // can A be freed by this function? + char* msg +) { + uint64_t node_id, com_id, n_deg, queue_head = 0, queue_tail = 0; + GrB_Vector s_list = NULL; + GrB_Matrix S_squished = NULL, S_new = NULL ; + GrB_Matrix C_squished = NULL, C_new = NULL ; + GrB_Matrix A_squished = NULL, A_new = NULL ; + + GrB_Descriptor desc = NULL ; + uint64_t n, n_new ; + GRB_TRY (GrB_Matrix_nrows (&n, *A)) ; + + GRB_TRY (GrB_Descriptor_new (&desc)) ; + GRB_TRY (GrB_set (desc, GxB_USE_INDICES, GxB_ROWINDEX_LIST)) ; + GRB_TRY (GrB_set (desc, GxB_USE_INDICES, GxB_COLINDEX_LIST)) ; + + GRB_TRY (GxB_Vector_diag (s_list, *S, 0, NULL)) ; + GRB_TRY (GrB_Vector_nvals (&n_new, s_list)) ; + + GRB_TRY (GrB_Matrix_new (&S_new, GrB_BOOL, n_new, n)) ; + GRB_TRY (GrB_Matrix_new (&C_new, GrB_BOOL, n_new, n_new)) ; + GRB_TRY (GrB_Matrix_new (&A_new, GrB_BOOL, n_new, n_new)) ; + + GRB_TRY (GrB_Vector_new (&s_list, GrB_BOOL, n)) ; + GRB_TRY (GrB_Matrix_new (&S_squished, GrB_BOOL, n, n_new)) ; + GRB_TRY (GrB_Matrix_new (&A_squished, GrB_BOOL, n, n_new)) ; + GRB_TRY (GrB_Matrix_new (&C_squished, GrB_BOOL, n, n_new)) ; + + GRB_TRY (GxB_Matrix_extract_Vector (S_squished, NULL, NULL, *S, NULL, s_list, desc)) ; + GRB_TRY (GrB_free (S)) ; + + GRB_TRY (GrB_mxm (C_squished, NULL, NULL, GxB_PLUS_FIRST_FP64, *C, S_squished, NULL)) ; + GRB_TRY (GrB_free (C)) ; + + GRB_TRY (GrB_mxm (A_squished, NULL, NULL, GxB_PLUS_FIRST_FP64, *A, S_squished, NULL)) ; + if (free_A) + GRB_TRY (GrB_free (A)) ; + + GRB_TRY (GrB_transpose (S_new, NULL, NULL, S_squished, NULL)) ; + GRB_TRY (GrB_free (&S_squished)) ; + + GRB_TRY (GrB_mxm (C_new, NULL, NULL, GxB_PLUS_FIRST_FP64, S_new, C_squished, NULL)) ; + GRB_TRY (GrB_free (&C_squished)) ; + + GRB_TRY (GrB_mxm (A_new, NULL, NULL, GxB_PLUS_FIRST_FP64, S_new, A_squished, NULL)) ; + GRB_TRY (GrB_free (&A_squished)) ; + + *S = S_new ; + *C = C_new ; + *A = A_new ; + + LG_FREE_WORK; + return GrB_SUCCESS ; +} +#undef LG_FREE_WORK +#define LG_FREE_WORK \ +{ \ + GrB_free (&C); \ + GrB_free (&S); \ + if (free_A) GrB_free (&A) ; \ + GrB_free (&mask); \ + GrB_free (°); \ + GrB_free (&c_deg); \ + GrB_free (&gain); \ + GrB_free (&x); \ + LAGraph_Free ((void **) queue, NULL); \ + LAGraph_Free ((void **) enqueued, NULL); \ + LAGraph_Free ((void **) community, NULL); \ + LAGraph_Free ((void **) sub_com, NULL); \ +} - return GrB_NOT_IMPLEMENTED; -#else - return GrB_NOT_IMPLEMENTED; -#endif +#undef LG_FREE_ALL +#define LG_FREE_ALL \ +{ \ + LG_FREE_WORK ; \ } -#else + int LAGraph_Leiden ( // output: @@ -301,536 +398,189 @@ int LAGraph_Leiden char *msg ) { +#ifdef LG_SUITESPARSE_GRAPHBLAS_V10 + uint64_t *queue = NULL; // circular buffer with next nodes to look up + bool *enqueued = NULL; // true if node is in the queue + uint64_t *community = NULL ; // array: i in community c[i] + GrB_Matrix C = NULL ; // contains readonly pointer to community array + GrB_Matrix mask = NULL ; // mask each community + GxB_Container cont = NULL; - //-------------------------------------------------------------------------- - // declare all workspace (must precede any LG_TRY/GRB_TRY calls so that - // LG_FREE_ALL can safely free them even on early exit) - //-------------------------------------------------------------------------- - - GrB_Vector k_vec = NULL ; - GrB_Matrix A_agg = NULL ; // owned coarsened graph (Phase 3) - GrB_Matrix A_new = NULL ; // next-level aggregate before ownership transfer - GrB_Matrix S_mat = NULL ; // temporary membership matrix (Phase 3) - GrB_Matrix A_temp = NULL ; // temporary for mxm (Phase 3) - GrB_Scalar one_scalar = NULL ; // FP64 scalar with value 1.0 for build_Scalar - GxB_Container cont = NULL ; // for unloading A_cur into raw CSR arrays - double *k_arr = NULL ; // k_arr[i] = degree of node i (current level) - int64_t *c_arr = NULL ; // c_arr[i] = Phase-1 community label - double *k_comm = NULL ; // k_comm[l] = total degree of community l - int64_t *c_p1 = NULL ; // c_p1[i] = Phase-1 parent community - int64_t *c_ref = NULL ; // c_ref[i] = refined sub-community label - double *k_ref_comm = NULL ; // k_ref_comm[l] = total degree of sub-community l - double *T_local = NULL ; // scratch: edge sums from node i to each community - int8_t *dirty = NULL ; // dirty[l] = 1 if T_local[l] was written - GrB_Index *dirty_list = NULL ; // list of community labels touched this node - GrB_Index *remap = NULL ; // remap[old_label] -> new contiguous label - int64_t *o_comm = NULL ; // o_comm[i] = community of original node i - GrB_Index *init_comm = NULL ; // init_comm[r] = initial c_arr for aggregate node r - - // Raw CSR pointers for inner-loop walks. On v10+ they are obtained by - // unloading A_cur into the SuiteSparse Container (zero-copy); ownership - // returns to GraphBLAS on reload (we then null them so LG_FREE_WORK - // doesn't double-free). On older versions they are allocated by us - // and (re)filled per level via GrB_Matrix_extractTuples + counting-sort. - GrB_Index *Ap = NULL ; // row pointers, size n_cur+1 - GrB_Index *Aj = NULL ; // column indices, size Anz - double *Ax = NULL ; // values, size Anz (only if !iso on v10) - GrB_Index *I_tup = NULL ; // raw row indices from extractTuples (v9 fallback) - GrB_Index *J_tup = NULL ; // raw col indices from extractTuples (v9 fallback) - double *X_tup = NULL ; // raw values from extractTuples (v9 fallback) - GrB_Index *cursor = NULL ; // scatter cursor for CSR build (v9 fallback) - GrB_Index *iota = NULL ; // [0,1,...,n-1] for vector/matrix build - GrB_Index Ap_cap = 0 ; // current allocated capacity of Ap (v9 fallback) - GrB_Index Anz_cap = 0 ; // current allocated capacity of Aj/Ax/tuples (v9 fallback) - - //-------------------------------------------------------------------------- - // check inputs - //-------------------------------------------------------------------------- - - LG_CLEAR_MSG ; - (void) seed ; // reserved for future randomized refinement - LG_ASSERT (c_handle != NULL, GrB_NULL_POINTER) ; - (*c_handle) = NULL ; - LG_TRY (LAGraph_CheckGraph (G, msg)) ; - LG_ASSERT_MSG ( - G->kind == LAGraph_ADJACENCY_UNDIRECTED || - (G->kind == LAGraph_ADJACENCY_DIRECTED && - G->is_symmetric_structure == LAGraph_TRUE), - LAGRAPH_NOT_CACHED, - "G must be undirected or have symmetric structure") ; - - GrB_Matrix A = G->A ; - GrB_Index n ; - GRB_TRY (GrB_Matrix_nrows (&n, A)) ; + uint64_t *sub_com = NULL ; // array: i in sub-community sub_com[i] + GrB_Matrix S = NULL ; // contains readonly pointer to sub-community array - // Degenerate: return immediately for empty (0-node) graph. - if (n == 0) - { - GRB_TRY (GrB_Vector_new (c_handle, GrB_INT64, 0)) ; - return (GrB_SUCCESS) ; - } + GrB_Vector deg = NULL ; // degree of each node + GrB_Vector c_deg = NULL ; // degree of each community - //-------------------------------------------------------------------------- - // allocate workspace (all arrays sized n; used for indices 0..n_cur-1) - //-------------------------------------------------------------------------- - - // void *workspace = NULL, p_workspace = NULL; - // LG_TRY (LAGraph_Malloc (&workspace, n, sizeof (double) * 4 + - // sizeof (int64_t) * 4 + sizeof (int8_t) + sizeof (GrB_Index) *4, msg)) ; - - LG_TRY (LAGraph_Malloc ((void **) &k_arr, n, sizeof (double), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &c_arr, n, sizeof (int64_t), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &k_comm, n, sizeof (double), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &c_p1, n, sizeof (int64_t), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &c_ref, n, sizeof (int64_t), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &k_ref_comm, n, sizeof (double), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &T_local, n, sizeof (double), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &dirty, n, sizeof (int8_t), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &dirty_list, n, sizeof (GrB_Index), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &remap, n, sizeof (GrB_Index), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &o_comm, n, sizeof (int64_t), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &init_comm, n, sizeof (GrB_Index), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &iota, n, sizeof (GrB_Index), msg)) ; - - // Reusable FP64 scalar with value 1.0 for GxB_Matrix_build_Scalar. - GRB_TRY (GrB_Scalar_new (&one_scalar, GrB_FP64)) ; - GRB_TRY (GrB_Scalar_setElement_FP64 (one_scalar, 1.0)) ; - - // Reusable container for unloading A_cur into raw CSR arrays per level - // (only available with the SuiteSparse v10+ Container API). - GRB_TRY (GxB_Container_new (&cont)) ; + GrB_Matrix A = NULL ; // G->A + GrB_Vector gain = NULL ; // gain from moving into a community + GrB_Vector x = NULL; // full boolean vector + GrB_Matrix node_map = NULL ; // maps aggregate node to nodes it contains + GrB_Matrix new_node_map = NULL ; // maps aggregate node to nodes it contains - //-------------------------------------------------------------------------- - // compute m = total edge weight / 2 from G->A (invariant under aggregation) - //-------------------------------------------------------------------------- - - double m = 0.0 ; - GRB_TRY (GrB_Matrix_reduce_FP64 (&m, NULL, GrB_PLUS_MONOID_FP64, - A, NULL)) ; - m /= 2.0 ; - - // Empty graph: return a singleton partition. - if (m == 0.0) - { - // c[i] = i for all i, built in one call instead of n setElement calls. - GRB_TRY (GrB_Vector_new (c_handle, GrB_INT64, n)) ; - GRB_TRY (GrB_Vector_assign_INT64 ( - *c_handle, NULL, NULL, 0, GrB_ALL, 0, NULL)) ; - GRB_TRY (GrB_Vector_apply_IndexOp_INT64 ( - *c_handle, NULL, NULL, GrB_ROWINDEX_INT64, *c_handle, 0, NULL)) ; - LG_FREE_WORK ; - return (GrB_SUCCESS) ; - } + uint64_t n, original_n; - //-------------------------------------------------------------------------- - // initialise multi-level state - // - // o_comm[i] = community of original node i (tracks the composition of - // all levels' refined partitions). Starts as identity. - // init_comm[r] = initial Phase-1 community for aggregate node r at the - // next level. First level: singleton start (c_arr[i] = i). - //-------------------------------------------------------------------------- - - for (GrB_Index i = 0 ; i < n ; i++) - { - iota[i] = i ; - o_comm[i] = (int64_t) i ; - init_comm[i] = i ; - } + GxB_Iterator neighbor_it = NULL; - // Duplicate G->A as FP64 so we own A_cur and may unload it via the - // container API expecting FP64 values. Done after the m == 0 early - // return to avoid an unused copy for empty-edge graphs. G->A may be - // any numeric type (BOOL on pattern-only matrices, INT*, FP32, ...); - // we typecast once here so the inner-loop CSR walks always read double. - GRB_TRY (GrB_Matrix_new (&A_agg, GrB_FP64, n, n)) ; - GRB_TRY (GrB_Matrix_assign (A_agg, NULL, NULL, A, - GrB_ALL, n, GrB_ALL, n, NULL)) ; - GrB_Matrix A_cur = A_agg ; - GrB_Index n_cur = n ; - - //========================================================================== - // OUTER AGGREGATION LOOP - //========================================================================== - - bool outer_changed = true ; - while (outer_changed) - { - outer_changed = false ; + bool free_A = false ; + A = G->A; + GRB_TRY (GrB_Matrix_nvals(&n, A)); + original_n = n; // Store original number of nodes + + // Input checking + LG_ASSERT_MSG (G->out_degree != NULL, GrB_EMPTY_OBJECT, + "G->out_degree must be defined") ; + LG_ASSERT_MSG( + G->kind == LAGraph_ADJACENCY_UNDIRECTED, GrB_INVALID_VALUE, + "G->A must be symmetric") ; + + // initialize queue + LG_TRY (LAGraph_Malloc ((void **) &enqueued, n, sizeof(bool), msg)) ; + LG_TRY (LAGraph_Malloc ((void **) &queue, n, sizeof (uint64_t), msg)) ; + + // Initialize GraphBLAS data structures + GRB_TRY (GxB_Iterator_new(&neighbor_it)); + GRB_TRY (GxB_Container_new (&cont)) ; + GRB_TRY (GrB_Vector_new(&x, GrB_FP64, n)); + GRB_TRY (GrB_Vector_new(°, GrB_FP64, n)); + GRB_TRY (GrB_Vector_new(&gain, GrB_FP64, n)); - //---------------------------------------------------------------------- - // Compute degrees k_arr[i] = sum of row i (includes self-loops in - // A_agg) and unload A_cur into raw CSR arrays Ap/Aj/Ax via the - // SuiteSparse Container API. This replaces O(n_cur) GrB_Col_extract - // + extractTuples calls per inner iteration with direct pointer walks. - // - // Force A_cur into a deterministic state before unloading: dense - // reduce target for k, sparse + row-major + non-iso + 64-bit indices - // for A. Container fields are then asserted to match expectations. - //---------------------------------------------------------------------- + GRB_TRY (GrB_Vector_new(&x, GrB_BOOL, n)); + GRB_TRY (GrB_assign (x, NULL, NULL, (bool) 0, GrB_ALL, n, NULL)); + GRB_TRY (GrB_Matrix_diag (&C, x, 0)) ; + GRB_TRY (GrB_Matrix_dup (&S, C)) ; - // 1) Compute degrees with GraphBLAS reduce *before* unloading the - // matrix. Zero-fill k_vec first so isolated rows produce 0.0 - // (otherwise reduce leaves them as missing entries). - GRB_TRY (GrB_Vector_new (&k_vec, GrB_FP64, n_cur)) ; - GRB_TRY (GrB_assign (k_vec, NULL, NULL, (double) 0.0, - GrB_ALL, n_cur, NULL)) ; - GRB_TRY (GrB_Matrix_reduce_Monoid (k_vec, NULL, GrB_PLUS_FP64, - GrB_PLUS_MONOID_FP64, A_cur, NULL)) ; - - // Unload k_vec into k_arr (dense FP64 array of length n_cur). - // The previous k_arr workspace allocation is replaced by a pointer - // owned by GraphBLAS until LAGraph_Free reclaims it. - LAGraph_Free ((void **) &k_arr, NULL) ; - do { - GrB_Type k_type = NULL ; - uint64_t k_n = 0, k_size = 0 ; - int k_handling = GrB_DEFAULT ; - void *k_void = NULL ; - GRB_TRY (GxB_Vector_unload (k_vec, &k_void, &k_type, &k_n, - &k_size, &k_handling, NULL)) ; - LG_ASSERT_MSG (k_type == GrB_FP64 && k_n == n_cur, - GrB_INVALID_VALUE, - "k_vec unload: unexpected type or length") ; - k_arr = (double *) k_void ; - } while (0); - - // 2) Force A_cur into the format we want, then unload into container. - // Hints: sparse, row-major, non-iso, 64-bit row pointers/indices. - // GxB_unload_Matrix_into_Container materializes pending work. - GRB_TRY (GrB_Matrix_set_INT32 ( - A_cur, GxB_SPARSE, GxB_SPARSITY_CONTROL)) ; - GRB_TRY (GrB_Matrix_set_INT32 ( - A_cur, GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT)) ; - GRB_TRY (GrB_Matrix_set_INT32 (A_cur, false, GxB_ISO)) ; - GRB_TRY (GrB_Matrix_set_INT32 (A_cur, 64, GxB_OFFSET_INTEGER_HINT)) ; - GRB_TRY (GrB_Matrix_set_INT32 (A_cur, 64, GxB_ROWINDEX_INTEGER_HINT)) ; - GRB_TRY (GrB_Matrix_set_INT32 (A_cur, 64, GxB_COLINDEX_INTEGER_HINT)) ; - GRB_TRY (GrB_wait (A_cur, GrB_MATERIALIZE)) ; - - GRB_TRY (GxB_unload_Matrix_into_Container (A_cur, cont, NULL)) ; - LG_ASSERT_MSG (cont->format == GxB_SPARSE, - GrB_INVALID_VALUE, "A_cur container is not sparse CSR") ; - LG_ASSERT_MSG (cont->orientation == GrB_ROWMAJOR, - GrB_INVALID_VALUE, "A_cur container is not row-major") ; - LG_ASSERT_MSG (!cont->iso, - GrB_INVALID_VALUE, "A_cur container unexpectedly iso") ; - - // Unload row pointers, column indices, and values from the container's - // internal vectors into raw arrays for the inner-loop CSR walks. - // We hold these arrays as our own until reload below. - GrB_Type pty = NULL, ity = NULL, xty = NULL ; - uint64_t pn = 0, in_ = 0, xn = 0 ; - uint64_t psz = 0, isz = 0, xsz = 0 ; - int ph = GrB_DEFAULT, ih = GrB_DEFAULT, xh = GrB_DEFAULT ; - void *pv = NULL, *iv = NULL, *xv = NULL ; - - GRB_TRY (GxB_Vector_unload (cont->p, &pv, &pty, &pn, &psz, &ph, NULL)); - GRB_TRY (GxB_Vector_unload (cont->i, &iv, &ity, &in_, &isz, &ih, NULL)); - GRB_TRY (GxB_Vector_unload (cont->x, &xv, &xty, &xn, &xsz, &xh, NULL)); - // Offsets must be 64-bit unsigned (we forced via INTEGER_HINT). - // Column indices may come back as either UINT64 or INT64 depending - // on SuiteSparse's internal choice; both have identical bit width - // and represent non-negative indices, so reinterpret cast is safe. - // Values must be FP64 since A_agg was constructed as FP64. - LG_ASSERT_MSG (pty == GrB_UINT64, - GrB_INVALID_VALUE, "container offsets are not 64-bit unsigned") ; - LG_ASSERT_MSG (ity == GrB_UINT64 || ity == GrB_INT64, - GrB_INVALID_VALUE, "container indices are not 64-bit") ; - LG_ASSERT_MSG (xty == GrB_FP64, - GrB_INVALID_VALUE, "container values are not FP64") ; - Ap = (GrB_Index *) pv ; - Aj = (GrB_Index *) iv ; - Ax = (double *) xv ; + // Initialize full degree vector from graph + GRB_TRY (GrB_assign (deg, NULL, NULL, 0, GrB_ALL, n, NULL)); + GRB_TRY (GrB_assign (deg, NULL, GrB_PLUS_FP64, G->out_degree, GrB_ALL, n, NULL)); - //---------------------------------------------------------------------- - // PHASE 1: Local Move Phase - // - // Initialise partition from init_comm (singletons on first level; - // induced Phase-1 partition on subsequent levels). - // Score: score(i->c) = T[c] - k[i]*k_comm[c]/m (self-loops skipped). - //---------------------------------------------------------------------- - memset (dirty, 0, n * sizeof (int8_t)) ; - memset (T_local, 0, n * sizeof (double)) ; - memset (k_comm, 0, n * sizeof (double)) ; - for (GrB_Index i = 0 ; i < n_cur ; i++) - { - c_arr[i] = (int64_t) init_comm[i] ; - k_comm[init_comm[i]] += k_arr[i] ; - } + double m = 0; + GRB_TRY (GrB_reduce (&m, NULL, GrB_PLUS_MONOID_FP64, deg, NULL)) ; + double m_inv2 = -1 / (2 * m); - bool changed = true ; - for (int p1_iter = 0 ; changed && p1_iter < LEIDEN_MAX_ITER ; p1_iter++) - { - changed = false ; - for (GrB_Index i = 0 ; i < n_cur ; i++) - { - double ki = k_arr[i] ; - if (ki == 0.0) continue ; - - int64_t ci = c_arr[i] ; - - GrB_Index row_begin = Ap[i] ; - GrB_Index row_end = Ap[i + 1] ; - if (row_begin == row_end) continue ; - - // Temporarily remove i from community ci. - k_comm[ci] -= ki ; - - GrB_Index ndirty = 0 ; - for (GrB_Index t = row_begin ; t < row_end ; t++) - { - GrB_Index j = Aj[t] ; - if (j == i) continue ; // skip self-loop (in A_agg) - int64_t cj = c_arr[j] ; - if (!dirty[cj]) - { - dirty[cj] = 1 ; - dirty_list[ndirty++] = (GrB_Index) cj ; - T_local[cj] = 0.0 ; - } - T_local[cj] += Ax[t] ; - } - - double T_ci = dirty[ci] ? T_local[ci] : 0.0 ; - double score_ci = T_ci - ki * k_comm[ci] / (m * 2) ; - double best_score = score_ci ; - int64_t best_c = ci ; - - for (GrB_Index d = 0 ; d < ndirty ; d++) - { - int64_t c_cand = (int64_t) dirty_list[d] ; - if (c_cand == ci) continue ; - double score = T_local[c_cand] - ki * k_comm[c_cand] / (m * 2) ; - if (score > best_score) - { - best_score = score ; - best_c = c_cand ; - } - } - - c_arr[i] = best_c ; - if (best_c == ci) - { - k_comm[ci] += ki ; - } - else - { - k_comm[best_c] += ki ; - changed = true ; - } - - for (GrB_Index d = 0 ; d < ndirty ; d++) - { - dirty[dirty_list[d]] = 0 ; - } - } - } + srand(seed); - //---------------------------------------------------------------------- - // PHASE 2: Refinement Phase (key Leiden addition) - // - // Save Phase-1 result. Restart each node in a singleton sub-community. - // Only allow moves within the same Phase-1 parent community. - //---------------------------------------------------------------------- + // Outer loop: repeat phases until convergence + double prev_modularity = 0; + double modularity_epsilon = 1e-6; - memcpy (c_p1, c_arr, n_cur * sizeof (int64_t)) ; - memcpy (k_ref_comm, k_arr, n_cur * sizeof (double)) ; + while (true) { + // Initialize queue in random order (seed-based) + GRB_TRY (GrB_Vector_extractTuples_BOOL (queue, NULL, &n, deg)) ; - for (GrB_Index i = 0 ; i < n_cur ; i++) - { - c_ref[i] = (int64_t) i ; + // TODO: shuffle with LAGraph Random instead. + for (uint64_t i = n - 1; i > 0; i--) { + uint64_t j = rand() % (i + 1); + uint64_t temp = queue[i]; + queue[i] = queue[j]; + queue[j] = temp; } + memset (enqueued, true, n) ; - changed = true ; - for (int p2_iter = 0 ; changed && p2_iter < LEIDEN_MAX_ITER ; p2_iter++) - { - changed = false ; - for (GrB_Index i = 0 ; i < n_cur ; i++) - { - double ki = k_arr[i] ; - if (ki == 0.0) continue ; - - int64_t pi = c_p1[i] ; - int64_t ci_ref = c_ref[i] ; - - GrB_Index row_begin = Ap[i] ; - GrB_Index row_end = Ap[i + 1] ; - if (row_begin == row_end) continue ; - - k_ref_comm[ci_ref] -= ki ; - - GrB_Index ndirty = 0 ; - for (GrB_Index t = row_begin ; t < row_end ; t++) - { - GrB_Index j = Aj[t] ; - if (j == i) continue ; // skip self-loop - if (c_p1[j] != pi) continue ; // cross-parent: skip - - int64_t cj_ref = c_ref[j] ; - if (!dirty[cj_ref]) - { - dirty[cj_ref] = 1 ; - dirty_list[ndirty++] = (GrB_Index) cj_ref ; - T_local[cj_ref] = 0.0 ; - } - T_local[cj_ref] += Ax[t] ; - } - - double T_ci_ref = dirty[ci_ref] ? T_local[ci_ref] : 0.0 ; - double score_ci_ref = T_ci_ref - ki * k_ref_comm[ci_ref] / (m * 2) ; - double best_score = score_ci_ref ; - int64_t best_c_ref = ci_ref ; - - for (GrB_Index d = 0 ; d < ndirty ; d++) - { - int64_t c_cand = (int64_t) dirty_list[d] ; - if (c_cand == ci_ref) continue ; - double score = T_local[c_cand] - ki * k_ref_comm[c_cand] / (m * 2) ; - if (score > best_score) - { - best_score = score ; - best_c_ref = c_cand ; - } - } - - c_ref[i] = best_c_ref ; - if (best_c_ref == ci_ref) - { - k_ref_comm[ci_ref] += ki ; - } - else - { - k_ref_comm[best_c_ref] += ki ; - changed = true ; - } - - for (GrB_Index d = 0 ; d < ndirty ; d++) - { - dirty[dirty_list[d]] = 0 ; - } - } - } + // get i vector from C and S to use as community and sub_com arrays + // keep inside C and S as read-only + do { + GrB_Type type; + uint64_t n_community, X_memsize; + int handling; + + GRB_TRY (GrB_set (C, GxB_SPARSE, GxB_SPARSITY_CONTROL)) ; + GRB_TRY (GrB_set (C, GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT)) ; + GRB_TRY (GxB_unload_Matrix_into_Container(C, cont, NULL)) ; + GRB_TRY (GxB_Vector_unload (cont->i, (void **) &community, &type, &n_community, + &X_memsize, &handling, NULL)) ; + GRB_TRY (GxB_Vector_load (cont->i, (void **) &community, type, n_community, + X_memsize, GxB_IS_READONLY, NULL)) ; + GRB_TRY (GxB_load_Matrix_from_Container(C, cont, NULL)) ; + + GRB_TRY (GrB_set (S, GxB_SPARSE, GxB_SPARSITY_CONTROL)) ; + GRB_TRY (GrB_set (S, GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT)) ; + GRB_TRY (GxB_unload_Matrix_into_Container (S, cont, NULL)) ; + GRB_TRY (GxB_Vector_unload (cont->i, (void **) &sub_com, &type, &n_community, + &X_memsize, &handling, NULL)) ; + GRB_TRY (GxB_Vector_load (cont->i, (void **) &sub_com, type, n_community, + X_memsize, GxB_IS_READONLY, NULL)) ; + GRB_TRY (GxB_load_Matrix_from_Container(S, cont, NULL)) ; + } while (0) ; //---------------------------------------------------------------------- - // Relabel c_ref to contiguous integers 0..K_ref-1 + // Phase 1: move nodes //---------------------------------------------------------------------- - - // Use n as sentinel ("not yet assigned"); safe because c_ref values - // are in 0..n_cur-1 < n. - for (GrB_Index i = 0 ; i < n_cur ; i++) remap[i] = n ; - - GrB_Index K_ref = 0 ; - for (GrB_Index i = 0 ; i < n_cur ; i++) - { - GrB_Index old_label = (GrB_Index) c_ref[i] ; - if (remap[old_label] == n) remap[old_label] = K_ref++ ; - c_ref[i] = (int64_t) remap[old_label] ; - } + LG_TRY (LG_Leiden_move_nodes (C, community, A, deg, queue, enqueued, m_inv2, msg)) ; //---------------------------------------------------------------------- - // Compute init_comm for next level. - // - // Aggregate node r (0..K_ref-1) is the refined community c_ref[i] for - // any i with that label. All such nodes have the same Phase-1 parent - // c_arr[i] (Leiden invariant), so we record that as the initial - // community for aggregate node r in the next outer iteration. + // Prep phase 2 //---------------------------------------------------------------------- - - for (GrB_Index i = 0 ; i < n_cur ; i++) - { - // c_ref[i] is in 0..K_ref-1 and init_comm is size n >= K_ref. - init_comm[(GrB_Index) c_ref[i]] = (GrB_Index) c_arr[i] ; + // Initialize queue + GRB_TRY (GrB_Vector_extractTuples_BOOL (queue, NULL, &n, deg)) ; + + // TODO: shuffle with LAGraph Random instead. + for (uint64_t i = n - 1; i > 0; i--) { + uint64_t j = rand() % (i + 1); + uint64_t temp = queue[i]; + queue[i] = queue[j]; + queue[j] = temp; } //---------------------------------------------------------------------- - // Compose o_comm: original node i now maps to aggregate community - // c_ref[o_comm[i]]. o_comm[i] is always a valid index in c_ref - // because it was set to some value in 0..n_cur-1 in the previous - // iteration (or to i on the first iteration). + // Phase 2: refine //---------------------------------------------------------------------- - - for (GrB_Index i = 0 ; i < n ; i++) - { - o_comm[i] = c_ref[o_comm[i]] ; - } + // Save Phase 1 communities as parent communities + LG_TRY (LG_Leiden_refinement (S, sub_com, C, community, c_deg, A, deg, queue, m_inv2, msg)) ; //---------------------------------------------------------------------- - // Reload A_cur from the container before any further GraphBLAS use - // (Phase 3 mxm or next-level unload). Ownership of Ap/Aj/Ax returns - // to GraphBLAS; we null our pointers so LG_FREE_WORK won't double-free. - // No-op on the v9 fallback (A_cur was never unloaded). + // Phase 3: aggregate //---------------------------------------------------------------------- - GRB_TRY (GxB_Vector_load (cont->p, (void **) &Ap, pty, - pn, psz, ph, NULL)) ; - Ap = NULL ; - GRB_TRY (GxB_Vector_load (cont->i, (void **) &Aj, ity, - in_, isz, ih, NULL)) ; - Aj = NULL ; - GRB_TRY (GxB_Vector_load (cont->x, (void **) &Ax, xty, - xn, xsz, xh, NULL)) ; - Ax = NULL ; - GRB_TRY (GxB_load_Matrix_from_Container (A_cur, cont, NULL)) ; + // Compute modularity to check convergence + double current_modularity; + LG_TRY (LAGr_AdjModularity(¤t_modularity, 1.0, A, S, msg)); - //---------------------------------------------------------------------- - // PHASE 3: Aggregation — build coarsened graph if communities merged - //---------------------------------------------------------------------- - - if (K_ref < n_cur) - { - outer_changed = true ; - - // S_mat: n_cur × K_ref indicator matrix; S[i, c_ref[i]] = 1 for - // every i. All values are 1.0, so use GxB_Matrix_build_Scalar - // and a shared scalar instead of materializing a 1.0-array. - // rows: iota (precomputed [0..n-1], reused) - // cols: c_ref reinterpret-cast to GrB_Index*. c_ref values - // are non-negative community labels in [0, K_ref); on - // all targeted platforms int64_t and uint64_t share the - // same width and representation for non-negative values. - GRB_TRY (GrB_Matrix_new (&S_mat, GrB_FP64, n_cur, K_ref)) ; - GRB_TRY (GxB_Matrix_build_Scalar (S_mat, iota, - (GrB_Index *) c_ref, one_scalar, n_cur)) ; - - // A_temp = A_cur * S (n_cur × K_ref) - GRB_TRY (GrB_Matrix_new (&A_temp, GrB_FP64, n_cur, K_ref)) ; - GRB_TRY (GrB_mxm (A_temp, NULL, NULL, - GrB_PLUS_TIMES_SEMIRING_FP64, A_cur, S_mat, NULL)) ; - - // A_new = S^T * A_temp (K_ref × K_ref) - GRB_TRY (GrB_Matrix_new (&A_new, GrB_FP64, K_ref, K_ref)) ; - GRB_TRY (GrB_mxm (A_new, NULL, NULL, - GrB_PLUS_TIMES_SEMIRING_FP64, S_mat, A_temp, GrB_DESC_T0)) ; - - GrB_free (&S_mat) ; S_mat = NULL ; - GrB_free (&A_temp) ; A_temp = NULL ; - GrB_free (&A_agg) ; // free previous level's aggregate graph - A_agg = A_new ; - A_new = NULL ; // ownership transferred to A_agg - A_cur = A_agg ; - n_cur = K_ref ; + if (fabs(current_modularity - prev_modularity) < modularity_epsilon) { + break; // converged + } + prev_modularity = current_modularity; + + LG_TRY (LG_Leiden_aggregate (&S, &C, &A, free_A, msg)) ; + // free internal pointers that were not freed + LG_TRY (LAGraph_Free ((void **) &community, msg)) ; + LG_TRY (LAGraph_Free ((void **) &sub_com, msg)) ; + GRB_TRY (GrB_Matrix_nvals (&n, A)) ; + + free_A = true; + if (node_map == NULL) { + node_map = S; + S = NULL; + } else { + GRB_TRY (GrB_Matrix_new (&new_node_map, GrB_BOOL, n, original_n)) ; + GRB_TRY (GrB_mxm (new_node_map, NULL, NULL, GxB_ANY_PAIR_BOOL, S, node_map, NULL)) ; + + GRB_TRY (GrB_free (&node_map)) ; + node_map = new_node_map; new_node_map = NULL; + GRB_TRY (GrB_free (&S)) ; } - // K_ref == n_cur: no communities merged this level → converged. - } - //-------------------------------------------------------------------------- - // Build output GrB_Vector from o_comm with move semantics: hand the - // o_comm buffer directly to GraphBLAS (no copy) and null our pointer so - // LG_FREE_WORK doesn't double-free. o_comm values are already relabeled - // 0..K_final-1 from the last iteration; the loaded vector is "full" - // (every index has a value), so set sparsity hint accordingly. - //-------------------------------------------------------------------------- + // Re-initialize degree from coarsened graph + GRB_TRY (GrB_free (°)) ; + GRB_TRY (GrB_Vector_new (°, GrB_FP64, n)) ; + GRB_TRY (GrB_assign(deg, NULL, NULL, 0.0, GrB_ALL, n, NULL)); + GRB_TRY (GrB_reduce(deg, NULL, NULL, GrB_PLUS_MONOID_FP64, A, NULL)); + GRB_TRY (GrB_Vector_resize (x, n)) ; + GRB_TRY (GrB_Matrix_diag (&S, x, 0)) ; + } // end outer while loop - GRB_TRY (GrB_Vector_new (c_handle, GrB_INT64, n)) ; - GRB_TRY (GrB_set (*c_handle, GxB_FULL, GxB_SPARSITY_CONTROL)) ; - GRB_TRY (GxB_Vector_load (*c_handle, (void **) &o_comm, GrB_INT64, - n, n * sizeof (int64_t), GrB_DEFAULT, NULL)) ; - o_comm = NULL ; // ownership transferred to *c_handle - LG_FREE_WORK ; - return (GrB_SUCCESS) ; -} + // Map refined sub-communities back to original nodes + GRB_TRY (GrB_Vector_new(c_handle, GrB_INT64, original_n)); + for (uint64_t i = 0; i < original_n; i++) { + GRB_TRY (GrB_Vector_setElement_INT64(*c_handle, community[i], i)); + } + + LG_FREE_WORK; + return GrB_SUCCESS; +#else + return GrB_NOT_IMPLEMENTED; #endif +} From d052ed783c7a6ff592155b09b20437f8b824b03e Mon Sep 17 00:00:00 2001 From: GomezGab Date: Mon, 10 Aug 2026 22:16:28 -0500 Subject: [PATCH 04/13] Leiden initial tests passing --- experimental/algorithm/LAGraph_Leiden.c | 116 +++++++++++++++--------- experimental/test/test_leiden.c | 1 + 2 files changed, 76 insertions(+), 41 deletions(-) diff --git a/experimental/algorithm/LAGraph_Leiden.c b/experimental/algorithm/LAGraph_Leiden.c index afae8284d6..ae690d2bd9 100644 --- a/experimental/algorithm/LAGraph_Leiden.c +++ b/experimental/algorithm/LAGraph_Leiden.c @@ -100,17 +100,15 @@ int LG_Leiden_move_nodes const GrB_Vector deg, // degree vector uint64_t *queue, // queue to use (contains all nodes) bool *enqueued, // nodes in queue (all at the start) - uint64_t m_inv2, // 1 / (2 * m) + double m_inv2, // -1 / (2 * m) char* msg ) { - uint64_t node_id, com_id, n_deg, queue_head = 0, queue_tail, n; + uint64_t node_id, com_id, n_deg, n; GrB_Vector x = NULL, gain = NULL; GrB_Vector c_deg = NULL ; GxB_Iterator it = NULL; GrB_Matrix_nrows (&n, A) ; - queue_tail = n; - - + uint64_t queue_head = 0, queue_tail = n, queue_size = n + 1; GRB_TRY (GxB_Iterator_new (&it)) ; GRB_TRY (GrB_Vector_new (&c_deg, GrB_FP64, n)) ; @@ -121,14 +119,18 @@ int LG_Leiden_move_nodes // TODO: decide if queue should be made inside this function. while (queue_head != queue_tail) { // while queue not empty node_id = queue[queue_head]; - queue_head = (queue_head + 1) % n; + queue_head = (queue_head + 1) % queue_size; enqueued [node_id] = false; double com_deg, n_deg; - com_id = community[node_id] ; + uint64_t n_neighbors; + com_id = community [node_id] ; + GRB_TRY (GrB_Vector_clear (x)); GRB_TRY (GrB_Vector_setElement_FP64 (x, 0.0, node_id)); GRB_TRY (GrB_vxm (x, NULL, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, x, A, NULL)) ; + GRB_TRY (GrB_Vector_nvals (&n_neighbors, x)) ; + if (n_neighbors == 1) continue; // skip singletons // give gain the number of edges connecting x to community i - GRB_TRY (GrB_vxm (gain, NULL, NULL, GxB_PLUS_SECOND_FP64, x, C, NULL)) ; + GRB_TRY (GrB_vxm (gain, NULL, NULL, GxB_PLUS_FIRST_FP64, x, C, NULL)) ; // momentarily remove node degree from the community total GRB_TRY (GrB_Vector_extractElement_FP64 (&n_deg, deg, node_id)); @@ -140,13 +142,16 @@ int LG_Leiden_move_nodes GRB_TRY (GrB_apply (gain, gain, GrB_PLUS_FP64, GrB_TIMES_FP64, c_deg, m_inv2 * n_deg, GrB_DESC_S)); uint64_t max_gain_c; double max_gain_val = -1e100; - for (uint64_t c = 0; c < n; c++) { - double g; - GRB_TRY (GrB_Vector_extractElement_FP64(&g, gain, c)); + GRB_TRY (GxB_Vector_Iterator_attach (it, gain, NULL)); + GrB_Info info = GxB_Vector_Iterator_seek (it, 0); + while (info == GrB_SUCCESS) { + uint64_t c = GxB_Vector_Iterator_getIndex (it); + double g = GxB_Iterator_get_FP64 (it); if (g > max_gain_val) { max_gain_val = g; max_gain_c = c; } + info = GxB_Vector_Iterator_next (it); } community [node_id] = max_gain_c; @@ -158,13 +163,13 @@ int LG_Leiden_move_nodes // push every neighbor in a different community that is not already in // the queue to the back of the queue GRB_TRY (GxB_Vector_Iterator_attach (it, x, NULL)); - GrB_Info info = GxB_Vector_Iterator_seek (it, 0); + info = GxB_Vector_Iterator_seek (it, 0); while (info == GrB_SUCCESS) { uint64_t neighbor_id = GxB_Vector_Iterator_getIndex (it); if (max_gain_c != community [neighbor_id] && !enqueued [neighbor_id]) { queue[queue_tail] = neighbor_id; - queue_tail = (queue_tail + 1) % n; + queue_tail = (queue_tail + 1) % queue_size; enqueued[neighbor_id] = true; } @@ -179,12 +184,11 @@ int LG_Leiden_move_nodes #define LG_FREE_WORK \ { \ GrB_free (&s_deg); \ - GrB_free (&c_list); \ GrB_free (&x_count); \ GrB_free (&k); \ GrB_free (&X); \ GrB_free (&C_t); \ - GrB_free (&mask); \ + GrB_free (&it); \ } #undef LG_FREE_ALL @@ -206,21 +210,21 @@ int LG_Leiden_refinement const GrB_Matrix A, // adjacency matrix const GrB_Vector deg, // degree vector uint64_t *queue, // queue to use (contains all nodes) - uint64_t m_inv2, // 1 / (2 * m) + double m_inv2, // -1 / (2 * m) char* msg ) { - uint64_t node_id, com_id, n_deg, queue_head = 0, queue_tail, n; - GrB_Matrix X = NULL, mask = NULL; + uint64_t node_id, com_id, n_deg, n; + GrB_Matrix X = NULL; GrB_Vector s_deg = NULL; - GrB_Vector c_list = NULL; GrB_Vector k = NULL; GrB_Matrix C_t = NULL ; // transpose of C GrB_Vector x_count = NULL; + GxB_Iterator it = NULL ; GRB_TRY (GrB_Matrix_nrows (&n, A)) ; - queue_tail = n; + uint64_t queue_head = 0, queue_tail = n, queue_size = n + 1; - // TODO: initialize x. decide if queue should be made inside this function. + // TODO: decide if queue should be made inside this function. GRB_TRY (GrB_Matrix_new (&C_t, GrB_FP64, n, n)) ; GRB_TRY (GrB_transpose (C_t, NULL, NULL, C, NULL)) ; @@ -232,33 +236,35 @@ int LG_Leiden_refinement GRB_TRY (GrB_Vector_new (&x_count, GrB_FP64, n)); GRB_TRY (GrB_assign (x_count, NULL, NULL, 0.0, GrB_ALL, n, NULL)) ; GRB_TRY (GrB_reduce (x_count, NULL, GrB_PLUS_FP64, GrB_PLUS_MONOID_FP64, C_t, GrB_DESC_T0)) ; - GRB_TRY (GrB_mxv (c_list, NULL, NULL, GxB_ANY_PAIR_BOOL, C_t, deg, NULL)) ; uint64_t n_communities; - GRB_TRY (GrB_Vector_nvals (&n_communities, c_list)) ; - GRB_TRY (GrB_Vector_new (&k, GrB_FP64, n_communities)) ; - GRB_TRY (GrB_Matrix_new (&mask, GrB_FP64, n_communities, n)) ; - GRB_TRY (GrB_Matrix_new (&X, GrB_FP64, n_communities, n)) ; + GRB_TRY (GrB_Vector_new (&k, GrB_FP64, n)) ; + GRB_TRY (GrB_Vector_new (&s_deg, GrB_FP64, n)) ; + GRB_TRY (GrB_Matrix_new (&X, GrB_FP64, n, n)) ; + GRB_TRY (GxB_Iterator_new(&it)) ; - GRB_TRY (GxB_Matrix_extract_Vector (mask, NULL, NULL, C_t, c_list, NULL, NULL)) ; - GRB_TRY (GrB_select (mask, NULL, NULL, GrB_VALUEGT_FP64, mask, 0.0, NULL)) ; + GRB_TRY (GrB_select (C_t, NULL, NULL, GrB_VALUEGT_FP64, C_t, 0.0, NULL)) ; // NOTE: I am wiring this so that it will be relatively easy to parallelize // this phase in the FUTURE + while (queue_head != queue_tail) { // while queue not empty double count; node_id = queue[queue_head]; GRB_TRY (GrB_Vector_extractElement_FP64 (&count, x_count, node_id)); if (count <= 0.0) continue; // skip singletons - queue_head = (queue_head + 1) % n; + queue_head = (queue_head + 1) % queue_size; double com_deg, n_deg; com_id = community [node_id] ; uint64_t sub_com_id = sub_com [node_id] ; + // TODO: this is probably much easier with extract and if we select + // nodes via their communities. + GRB_TRY (GrB_Matrix_clear (X)); GRB_TRY (GrB_Matrix_setElement_FP64 (X, 0.0, com_id, node_id)); - GRB_TRY (GrB_mxm (X, mask, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, X, A, GrB_DESC_S)) ; + GRB_TRY (GrB_mxm (X, C_t, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, X, A, GrB_DESC_S)) ; // give X the number of edges connecting x to community i GRB_TRY (GrB_mxm (X, NULL, NULL, GxB_PLUS_SECOND_FP64, X, S, NULL)) ; @@ -273,8 +279,33 @@ int LG_Leiden_refinement GRB_TRY (GrB_mxm (X, X, GrB_PLUS_FP64, GrB_PLUS_TIMES_SEMIRING_FP64, (GrB_Matrix) s_deg, (GrB_Matrix) k, GrB_DESC_ST1)) ; - uint64_t new_c; - // TODO: select from X (those with gain >= current and choose randomly) + uint64_t new_c = 0; + + // FUTURE: This part will have to change if we want parrallel clusters + // to work + + // Select from X (those with gain >= current and choose randomly) + double current_gain, total_gain; + GRB_TRY (GrB_Matrix_extractElement_FP64 (¤t_gain, X, com_id, node_id)) ; + GRB_TRY (GrB_apply (X, NULL, NULL, GrB_MINUS_FP64, X, current_gain, NULL)) ; + GRB_TRY (GrB_select (X, NULL, NULL, GrB_VALUELT_FP64, X, 0.0, NULL)) ; + GRB_TRY (GrB_apply (X, NULL, NULL, GxB_POW_FP64, M_E, X, NULL)) ; + GRB_TRY (GrB_reduce (&total_gain, NULL, GrB_PLUS_MONOID_FP64, X, NULL)) ; + double rand_f = drand48 () * total_gain ; + + GRB_TRY (GxB_Matrix_Iterator_attach(it, X, NULL)) ; + GrB_Info info = GxB_Matrix_Iterator_seek (it, 0); + total_gain = 0; + while (info == GrB_SUCCESS) { + total_gain += GxB_Iterator_get_FC64 (it) ; + + if (total_gain > rand_f) { + GxB_Matrix_Iterator_getIndex(it, NULL, &new_c) ; + break; + } + + info = GxB_Vector_Iterator_next (it); + } sub_com [node_id] = new_c; // s_deg [new_c] += ndeg @@ -376,10 +407,10 @@ int LG_Leiden_aggregate GrB_free (&c_deg); \ GrB_free (&gain); \ GrB_free (&x); \ - LAGraph_Free ((void **) queue, NULL); \ - LAGraph_Free ((void **) enqueued, NULL); \ - LAGraph_Free ((void **) community, NULL); \ - LAGraph_Free ((void **) sub_com, NULL); \ + LAGraph_Free ((void **) &queue, NULL); \ + LAGraph_Free ((void **) &enqueued, NULL); \ + LAGraph_Free ((void **) &community, NULL); \ + LAGraph_Free ((void **) &sub_com, NULL); \ } #undef LG_FREE_ALL @@ -424,19 +455,20 @@ int LAGraph_Leiden bool free_A = false ; A = G->A; - GRB_TRY (GrB_Matrix_nvals(&n, A)); + GRB_TRY (GrB_Matrix_nrows (&n, A)); original_n = n; // Store original number of nodes // Input checking - LG_ASSERT_MSG (G->out_degree != NULL, GrB_EMPTY_OBJECT, + LG_ASSERT_MSG (G->out_degree != NULL, LAGRAPH_NOT_CACHED, "G->out_degree must be defined") ; LG_ASSERT_MSG( - G->kind == LAGraph_ADJACENCY_UNDIRECTED, GrB_INVALID_VALUE, + G->is_symmetric_structure == LAGraph_TRUE, + LAGRAPH_SYMMETRIC_STRUCTURE_REQUIRED, "G->A must be symmetric") ; // initialize queue LG_TRY (LAGraph_Malloc ((void **) &enqueued, n, sizeof(bool), msg)) ; - LG_TRY (LAGraph_Malloc ((void **) &queue, n, sizeof (uint64_t), msg)) ; + LG_TRY (LAGraph_Malloc ((void **) &queue, n + 1, sizeof (uint64_t), msg)) ; // Initialize GraphBLAS data structures GRB_TRY (GxB_Iterator_new(&neighbor_it)); @@ -451,7 +483,7 @@ int LAGraph_Leiden GRB_TRY (GrB_Matrix_dup (&S, C)) ; // Initialize full degree vector from graph - GRB_TRY (GrB_assign (deg, NULL, NULL, 0, GrB_ALL, n, NULL)); + GRB_TRY (GrB_assign (deg, NULL, NULL, 0.0, GrB_ALL, n, NULL)); GRB_TRY (GrB_assign (deg, NULL, GrB_PLUS_FP64, G->out_degree, GrB_ALL, n, NULL)); @@ -487,6 +519,7 @@ int LAGraph_Leiden GRB_TRY (GrB_set (C, GxB_SPARSE, GxB_SPARSITY_CONTROL)) ; GRB_TRY (GrB_set (C, GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT)) ; + GRB_TRY (GrB_set (C, 64, GxB_COLINDEX_INTEGER_HINT)) ; GRB_TRY (GxB_unload_Matrix_into_Container(C, cont, NULL)) ; GRB_TRY (GxB_Vector_unload (cont->i, (void **) &community, &type, &n_community, &X_memsize, &handling, NULL)) ; @@ -496,6 +529,7 @@ int LAGraph_Leiden GRB_TRY (GrB_set (S, GxB_SPARSE, GxB_SPARSITY_CONTROL)) ; GRB_TRY (GrB_set (S, GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT)) ; + GRB_TRY (GrB_set (S, 64, GxB_COLINDEX_INTEGER_HINT)) ; GRB_TRY (GxB_unload_Matrix_into_Container (S, cont, NULL)) ; GRB_TRY (GxB_Vector_unload (cont->i, (void **) &sub_com, &type, &n_community, &X_memsize, &handling, NULL)) ; diff --git a/experimental/test/test_leiden.c b/experimental/test/test_leiden.c index b2dd3c3728..2a429ecb81 100644 --- a/experimental/test/test_leiden.c +++ b/experimental/test/test_leiden.c @@ -64,6 +64,7 @@ void test_Leiden (void) // Ensure symmetry cache is populated (required by some checks). OK (LAGraph_Cached_IsSymmetricStructure (G, msg)) ; + OK (LAGraph_Cached_OutDegree (G, msg)) ; uint64_t seed = 0 ; //unused GrB_Vector c = NULL ; From 7314b9996b2244a2eb2ede93fdf316245b23ee76 Mon Sep 17 00:00:00 2001 From: GomezGab Date: Tue, 11 Aug 2026 13:24:55 -0500 Subject: [PATCH 05/13] bug fixes for Leiden --- experimental/algorithm/LAGraph_Leiden.c | 95 ++++++++++++++++--------- 1 file changed, 62 insertions(+), 33 deletions(-) diff --git a/experimental/algorithm/LAGraph_Leiden.c b/experimental/algorithm/LAGraph_Leiden.c index ae690d2bd9..393775197e 100644 --- a/experimental/algorithm/LAGraph_Leiden.c +++ b/experimental/algorithm/LAGraph_Leiden.c @@ -21,7 +21,7 @@ #include #include -#define LEIDEN_MAX_ITER 100 +#define LEIDEN_MAX_ITER 20 //------------------------------------------------------------------------------ // The Leiden algorithm is a modularity-based community detection method that @@ -252,10 +252,10 @@ int LG_Leiden_refinement while (queue_head != queue_tail) { // while queue not empty double count; node_id = queue[queue_head]; + queue_head = (queue_head + 1) % queue_size; GRB_TRY (GrB_Vector_extractElement_FP64 (&count, x_count, node_id)); if (count <= 0.0) continue; // skip singletons - queue_head = (queue_head + 1) % queue_size; double com_deg, n_deg; com_id = community [node_id] ; uint64_t sub_com_id = sub_com [node_id] ; @@ -279,7 +279,7 @@ int LG_Leiden_refinement GRB_TRY (GrB_mxm (X, X, GrB_PLUS_FP64, GrB_PLUS_TIMES_SEMIRING_FP64, (GrB_Matrix) s_deg, (GrB_Matrix) k, GrB_DESC_ST1)) ; - uint64_t new_c = 0; + uint64_t new_c = com_id, row = 0; // FUTURE: This part will have to change if we want parrallel clusters // to work @@ -288,7 +288,7 @@ int LG_Leiden_refinement double current_gain, total_gain; GRB_TRY (GrB_Matrix_extractElement_FP64 (¤t_gain, X, com_id, node_id)) ; GRB_TRY (GrB_apply (X, NULL, NULL, GrB_MINUS_FP64, X, current_gain, NULL)) ; - GRB_TRY (GrB_select (X, NULL, NULL, GrB_VALUELT_FP64, X, 0.0, NULL)) ; + GRB_TRY (GrB_select (X, NULL, NULL, GrB_VALUEGE_FP64, X, 0.0, NULL)) ; GRB_TRY (GrB_apply (X, NULL, NULL, GxB_POW_FP64, M_E, X, NULL)) ; GRB_TRY (GrB_reduce (&total_gain, NULL, GrB_PLUS_MONOID_FP64, X, NULL)) ; double rand_f = drand48 () * total_gain ; @@ -297,14 +297,14 @@ int LG_Leiden_refinement GrB_Info info = GxB_Matrix_Iterator_seek (it, 0); total_gain = 0; while (info == GrB_SUCCESS) { - total_gain += GxB_Iterator_get_FC64 (it) ; + total_gain += GxB_Iterator_get_FP64 (it) ; if (total_gain > rand_f) { - GxB_Matrix_Iterator_getIndex(it, NULL, &new_c) ; + GxB_Matrix_Iterator_getIndex(it, &row, &new_c) ; break; } - info = GxB_Vector_Iterator_next (it); + info = GxB_Matrix_Iterator_next (it); } sub_com [node_id] = new_c; @@ -319,6 +319,7 @@ int LG_Leiden_refinement #define LG_FREE_WORK \ { \ GrB_free (&desc) ; \ + GrB_free (&x) ; \ GrB_free (&s_list) ; \ GrB_free (&S_squished) ; \ GrB_free (&C_squished) ; \ @@ -345,7 +346,7 @@ int LG_Leiden_aggregate char* msg ) { uint64_t node_id, com_id, n_deg, queue_head = 0, queue_tail = 0; - GrB_Vector s_list = NULL; + GrB_Vector s_list = NULL, x = NULL; GrB_Matrix S_squished = NULL, S_new = NULL ; GrB_Matrix C_squished = NULL, C_new = NULL ; GrB_Matrix A_squished = NULL, A_new = NULL ; @@ -358,22 +359,25 @@ int LG_Leiden_aggregate GRB_TRY (GrB_set (desc, GxB_USE_INDICES, GxB_ROWINDEX_LIST)) ; GRB_TRY (GrB_set (desc, GxB_USE_INDICES, GxB_COLINDEX_LIST)) ; - GRB_TRY (GxB_Vector_diag (s_list, *S, 0, NULL)) ; + GRB_TRY (GrB_Vector_new (&s_list, GrB_BOOL, n)) ; + GRB_TRY (GrB_Vector_new (&x, GrB_BOOL, n)) ; + GRB_TRY (GrB_assign (x, NULL, NULL, (bool) 0, GrB_ALL, n, NULL)) ; + GRB_TRY (GrB_vxm (s_list, NULL, NULL, GxB_ANY_PAIR_BOOL, x, *S, NULL)) ; + GRB_TRY (GrB_Vector_nvals (&n_new, s_list)) ; GRB_TRY (GrB_Matrix_new (&S_new, GrB_BOOL, n_new, n)) ; GRB_TRY (GrB_Matrix_new (&C_new, GrB_BOOL, n_new, n_new)) ; - GRB_TRY (GrB_Matrix_new (&A_new, GrB_BOOL, n_new, n_new)) ; + GRB_TRY (GrB_Matrix_new (&A_new, GrB_FP64, n_new, n_new)) ; - GRB_TRY (GrB_Vector_new (&s_list, GrB_BOOL, n)) ; GRB_TRY (GrB_Matrix_new (&S_squished, GrB_BOOL, n, n_new)) ; - GRB_TRY (GrB_Matrix_new (&A_squished, GrB_BOOL, n, n_new)) ; GRB_TRY (GrB_Matrix_new (&C_squished, GrB_BOOL, n, n_new)) ; + GRB_TRY (GrB_Matrix_new (&A_squished, GrB_FP64, n, n_new)) ; GRB_TRY (GxB_Matrix_extract_Vector (S_squished, NULL, NULL, *S, NULL, s_list, desc)) ; GRB_TRY (GrB_free (S)) ; - GRB_TRY (GrB_mxm (C_squished, NULL, NULL, GxB_PLUS_FIRST_FP64, *C, S_squished, NULL)) ; + GRB_TRY (GrB_mxm (C_squished, NULL, NULL, GxB_ANY_PAIR_BOOL, *C, S_squished, NULL)) ; GRB_TRY (GrB_free (C)) ; GRB_TRY (GrB_mxm (A_squished, NULL, NULL, GxB_PLUS_FIRST_FP64, *A, S_squished, NULL)) ; @@ -383,10 +387,10 @@ int LG_Leiden_aggregate GRB_TRY (GrB_transpose (S_new, NULL, NULL, S_squished, NULL)) ; GRB_TRY (GrB_free (&S_squished)) ; - GRB_TRY (GrB_mxm (C_new, NULL, NULL, GxB_PLUS_FIRST_FP64, S_new, C_squished, NULL)) ; + GRB_TRY (GrB_mxm (C_new, NULL, NULL, GxB_ANY_PAIR_BOOL, S_new, C_squished, NULL)) ; GRB_TRY (GrB_free (&C_squished)) ; - GRB_TRY (GrB_mxm (A_new, NULL, NULL, GxB_PLUS_FIRST_FP64, S_new, A_squished, NULL)) ; + GRB_TRY (GrB_mxm (A_new, NULL, NULL, GxB_PLUS_SECOND_FP64, S_new, A_squished, NULL)) ; GRB_TRY (GrB_free (&A_squished)) ; *S = S_new ; @@ -497,16 +501,20 @@ int LAGraph_Leiden double prev_modularity = 0; double modularity_epsilon = 1e-6; - while (true) { + for (int count = 0; count < LEIDEN_MAX_ITER; count++) { + printf ("Iteration %d\n", count) ; // Initialize queue in random order (seed-based) - GRB_TRY (GrB_Vector_extractTuples_BOOL (queue, NULL, &n, deg)) ; + uint64_t queue_len = n; + GRB_TRY (GrB_Vector_extractTuples_FP64 (queue, NULL, &queue_len, deg)) ; // TODO: shuffle with LAGraph Random instead. - for (uint64_t i = n - 1; i > 0; i--) { - uint64_t j = rand() % (i + 1); - uint64_t temp = queue[i]; - queue[i] = queue[j]; - queue[j] = temp; + if (queue_len > 1) { + for (uint64_t i = queue_len - 1; i > 0; i--) { + uint64_t j = rand() % (i + 1); + uint64_t temp = queue[i]; + queue[i] = queue[j]; + queue[j] = temp; + } } memset (enqueued, true, n) ; @@ -520,6 +528,7 @@ int LAGraph_Leiden GRB_TRY (GrB_set (C, GxB_SPARSE, GxB_SPARSITY_CONTROL)) ; GRB_TRY (GrB_set (C, GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT)) ; GRB_TRY (GrB_set (C, 64, GxB_COLINDEX_INTEGER_HINT)) ; + GrB_wait (C, GrB_MATERIALIZE) ; GRB_TRY (GxB_unload_Matrix_into_Container(C, cont, NULL)) ; GRB_TRY (GxB_Vector_unload (cont->i, (void **) &community, &type, &n_community, &X_memsize, &handling, NULL)) ; @@ -547,14 +556,17 @@ int LAGraph_Leiden // Prep phase 2 //---------------------------------------------------------------------- // Initialize queue - GRB_TRY (GrB_Vector_extractTuples_BOOL (queue, NULL, &n, deg)) ; + queue_len = n; + GRB_TRY (GrB_Vector_extractTuples_FP64 (queue, NULL, &queue_len, deg)) ; // TODO: shuffle with LAGraph Random instead. - for (uint64_t i = n - 1; i > 0; i--) { - uint64_t j = rand() % (i + 1); - uint64_t temp = queue[i]; - queue[i] = queue[j]; - queue[j] = temp; + if (queue_len > 1) { + for (uint64_t i = queue_len - 1; i > 0; i--) { + uint64_t j = rand() % (i + 1); + uint64_t temp = queue[i]; + queue[i] = queue[j]; + queue[j] = temp; + } } //---------------------------------------------------------------------- @@ -569,8 +581,9 @@ int LAGraph_Leiden // Compute modularity to check convergence double current_modularity; - LG_TRY (LAGr_AdjModularity(¤t_modularity, 1.0, A, S, msg)); + LG_TRY (LAGr_AdjModularity (¤t_modularity, 1.0, A, S, msg)); + ASSERT (current_modularity >= prev_modularity) ; if (fabs(current_modularity - prev_modularity) < modularity_epsilon) { break; // converged } @@ -580,7 +593,7 @@ int LAGraph_Leiden // free internal pointers that were not freed LG_TRY (LAGraph_Free ((void **) &community, msg)) ; LG_TRY (LAGraph_Free ((void **) &sub_com, msg)) ; - GRB_TRY (GrB_Matrix_nvals (&n, A)) ; + GRB_TRY (GrB_Matrix_nrows (&n, A)) ; free_A = true; if (node_map == NULL) { @@ -602,13 +615,29 @@ int LAGraph_Leiden GRB_TRY (GrB_reduce(deg, NULL, NULL, GrB_PLUS_MONOID_FP64, A, NULL)); GRB_TRY (GrB_Vector_resize (x, n)) ; GRB_TRY (GrB_Matrix_diag (&S, x, 0)) ; + #ifndef NDEBUG + LG_TRY (LAGr_AdjModularity (¤t_modularity, 1.0, A, S, msg)); + ASSERT (current_modularity == prev_modularity) ; + #endif } // end outer while loop - // Map refined sub-communities back to original nodes + // Map refined communities back to original nodes GRB_TRY (GrB_Vector_new(c_handle, GrB_INT64, original_n)); - for (uint64_t i = 0; i < original_n; i++) { - GRB_TRY (GrB_Vector_setElement_INT64(*c_handle, community[i], i)); + GRB_TRY (GrB_free (&C)) ; + if (node_map == NULL) { + ASSERT (n == original_n) ; + GRB_TRY (GxB_Vector_load (*c_handle, (void **) &community, GrB_INT64, n, + sizeof (int64_t) * n, GrB_DEFAULT, NULL)); + } else { + GRB_TRY (GxB_Matrix_Iterator_attach (neighbor_it, node_map, NULL)) ; + GrB_Info info = GxB_Matrix_Iterator_seek (neighbor_it, 0) ; + while (info == GrB_SUCCESS) { + uint64_t coarse_id, orig_id; + GxB_Matrix_Iterator_getIndex (neighbor_it, &coarse_id, &orig_id) ; + GRB_TRY (GrB_Vector_setElement_INT64 (*c_handle, community[coarse_id], orig_id)) ; + info = GxB_Matrix_Iterator_next (neighbor_it) ; + } } LG_FREE_WORK; From 621f5ba451074138370418033acc1daec560f577 Mon Sep 17 00:00:00 2001 From: GomezGab Date: Tue, 11 Aug 2026 19:24:20 -0500 Subject: [PATCH 06/13] more bug fixes --- .../algorithm/LAGr_ModularityMatrix.c | 4 +- experimental/algorithm/LAGraph_Leiden.c | 24 ++- experimental/benchmark/leiden_demo.c | 90 ++++++++ experimental/test/test_leiden.c | 196 ++++++++++++------ 4 files changed, 241 insertions(+), 73 deletions(-) create mode 100644 experimental/benchmark/leiden_demo.c diff --git a/experimental/algorithm/LAGr_ModularityMatrix.c b/experimental/algorithm/LAGr_ModularityMatrix.c index 78ccc2787b..d14f473865 100644 --- a/experimental/algorithm/LAGr_ModularityMatrix.c +++ b/experimental/algorithm/LAGr_ModularityMatrix.c @@ -99,9 +99,9 @@ int LAGr_AdjModularity( double inv_m = -gamma / (2.0 * m); // sum degrees per community - GRB_TRY (GrB_vxm (k, NULL, NULL, GxB_PLUS_FIRST_INT64, k, S, NULL)) ; + GRB_TRY (GrB_vxm (k, NULL, NULL, GxB_PLUS_FIRST_FP64, k, S, NULL)) ; // square - GRB_TRY (GrB_assign (k, NULL, GrB_TIMES_INT64, k, GrB_ALL, n, NULL)) ; + GRB_TRY (GrB_assign (k, NULL, GrB_TIMES_FP64, k, GrB_ALL, n, NULL)) ; // sum all of the products // this computed \sum_{i,j} (k_{i}k_{j}\delta(\sigma_i\sigma_j)) diff --git a/experimental/algorithm/LAGraph_Leiden.c b/experimental/algorithm/LAGraph_Leiden.c index 393775197e..db2b628257 100644 --- a/experimental/algorithm/LAGraph_Leiden.c +++ b/experimental/algorithm/LAGraph_Leiden.c @@ -117,6 +117,7 @@ int LG_Leiden_move_nodes GRB_TRY (GrB_mxv(c_deg, NULL, NULL, GxB_PLUS_SECOND_FP64, C, deg, NULL)); // TODO: decide if queue should be made inside this function. + int64_t count = 0; while (queue_head != queue_tail) { // while queue not empty node_id = queue[queue_head]; queue_head = (queue_head + 1) % queue_size; @@ -155,7 +156,7 @@ int LG_Leiden_move_nodes } community [node_id] = max_gain_c; - // c_deg [node_id] -= ndeg + // c_deg [node_id] += ndeg GRB_TRY (GrB_assign (c_deg, NULL, GrB_PLUS_FP64, n_deg, &max_gain_c, 1, NULL)); if (max_gain_c == com_id) continue; // no change in community @@ -469,6 +470,14 @@ int LAGraph_Leiden G->is_symmetric_structure == LAGraph_TRUE, LAGRAPH_SYMMETRIC_STRUCTURE_REQUIRED, "G->A must be symmetric") ; + LG_ASSERT_MSG (G->emin_state != LAGraph_STATE_UNKNOWN, LAGRAPH_NOT_CACHED, + "G->emin must be defined") ; + // cast to FP64 if needed + double min_val = 0.0; + GRB_TRY (GrB_Scalar_extractElement_FP64 (&min_val, G->emin)) ; + LG_ASSERT_MSG (min_val >= 0.0, GrB_INVALID_VALUE, + "G->emin must be non-negative") ; + // TODO: check numeric, and should we allow 0 weight edges? // initialize queue LG_TRY (LAGraph_Malloc ((void **) &enqueued, n, sizeof(bool), msg)) ; @@ -488,11 +497,16 @@ int LAGraph_Leiden // Initialize full degree vector from graph GRB_TRY (GrB_assign (deg, NULL, NULL, 0.0, GrB_ALL, n, NULL)); - GRB_TRY (GrB_assign (deg, NULL, GrB_PLUS_FP64, G->out_degree, GrB_ALL, n, NULL)); + GRB_TRY (GrB_reduce(deg, NULL, NULL, GrB_PLUS_MONOID_FP64, A, NULL)); + // don't use out degree unless A is bool adj matrix + // GRB_TRY (GrB_assign (deg, NULL, GrB_PLUS_FP64, G->out_degree, GrB_ALL, n, NULL)); double m = 0; GRB_TRY (GrB_reduce (&m, NULL, GrB_PLUS_MONOID_FP64, deg, NULL)) ; + LG_ASSERT_MSG (isfinite(m), GrB_INVALID_VALUE, + "Matrix must reduce to a finite value") ; + double m_inv2 = -1 / (2 * m); srand(seed); @@ -502,7 +516,7 @@ int LAGraph_Leiden double modularity_epsilon = 1e-6; for (int count = 0; count < LEIDEN_MAX_ITER; count++) { - printf ("Iteration %d\n", count) ; + printf ("Iteration %d, node count %lu \n", count, n) ; // Initialize queue in random order (seed-based) uint64_t queue_len = n; GRB_TRY (GrB_Vector_extractTuples_FP64 (queue, NULL, &queue_len, deg)) ; @@ -583,7 +597,7 @@ int LAGraph_Leiden double current_modularity; LG_TRY (LAGr_AdjModularity (¤t_modularity, 1.0, A, S, msg)); - ASSERT (current_modularity >= prev_modularity) ; + ASSERT (current_modularity >= prev_modularity - modularity_epsilon) ; if (fabs(current_modularity - prev_modularity) < modularity_epsilon) { break; // converged } @@ -617,7 +631,7 @@ int LAGraph_Leiden GRB_TRY (GrB_Matrix_diag (&S, x, 0)) ; #ifndef NDEBUG LG_TRY (LAGr_AdjModularity (¤t_modularity, 1.0, A, S, msg)); - ASSERT (current_modularity == prev_modularity) ; + ASSERT (fabs (current_modularity - prev_modularity) < modularity_epsilon) ; #endif } // end outer while loop diff --git a/experimental/benchmark/leiden_demo.c b/experimental/benchmark/leiden_demo.c new file mode 100644 index 0000000000..f61c55f406 --- /dev/null +++ b/experimental/benchmark/leiden_demo.c @@ -0,0 +1,90 @@ +//------------------------------------------------------------------------------ +// LAGraph/experimental/benchmark/leiden_demo.c: +// run Leiden on a sanitized input graph and report modularity +//------------------------------------------------------------------------------ + +// LAGraph, (c) 2026 by The LAGraph Contributors, All Rights Reserved. +// SPDX-License-Identifier: BSD-2-Clause + +// Usage: +// leiden_demo < matrixmarketfile.mtx +// leiden_demo matrixmarketfile.mtx +// leiden_demo matrixmarketfile.grb + +#include "../../src/benchmark/LAGraph_demo.h" +#include "LAGraphX.h" + +#define LG_FREE_ALL \ +{ \ + LAGraph_Delete (&G, NULL) ; \ + GrB_free (&c) ; \ +} + +int main (int argc, char **argv) +{ +#if LAGRAPH_SUITESPARSE + + char msg [LAGRAPH_MSG_LEN] ; + msg [0] = '\0' ; + LAGraph_Graph G = NULL ; + GrB_Vector c = NULL ; + + bool burble = false ; + LAGRAPH_TRY (demo_init (burble)) ; + + // read and sanitize: + // - make symmetric (A = A + A') + // - remove self edges + // - keep numeric values (not structural) + // - cast to FP64 + // - ensure all values are nonnegative + LAGRAPH_TRY (readproblem (&G, NULL, + true, true, false, GrB_FP64, true, argc, argv)) ; + + LAGRAPH_TRY (LAGraph_Cached_IsSymmetricStructure (G, msg)) ; + LAGRAPH_TRY (LAGraph_Cached_OutDegree (G, msg)) ; + LAGRAPH_TRY (LAGraph_Cached_EMin (G, msg)) ; + + double emin = 0 ; + GRB_TRY (GrB_Scalar_extractElement_FP64 (&emin, G->emin)) ; + if (emin < 0) + { + printf ("error: sanitized graph has negative edge weights (emin=%g)\n", + emin) ; + LG_FREE_ALL ; + LAGRAPH_TRY (LAGraph_Finalize (msg)) ; + return (GrB_INVALID_VALUE) ; + } + + uint64_t seed = 0 ; + double t = LAGraph_WallClockTime ( ) ; + LAGRAPH_TRY (LAGraph_Leiden (&c, G, seed, msg)) ; + t = LAGraph_WallClockTime ( ) - t ; + + double Q = 0 ; + LAGRAPH_TRY (LAGr_Modularity (&Q, 1.0, c, G, msg)) ; + + GrB_Index n = 0, nvals = 0 ; + GRB_TRY (GrB_Matrix_nrows (&n, G->A)) ; + GRB_TRY (GrB_Vector_nvals (&nvals, c)) ; + + int64_t max_label = -1 ; + GRB_TRY (GrB_Vector_reduce_INT64 ( + &max_label, NULL, GxB_MAX_INT64_MONOID, c, NULL)) ; + GrB_Index n_communities = (max_label < 0) ? 0 : (GrB_Index) (max_label + 1) ; + + printf ("Leiden complete\n") ; + printf (" n : %llu\n", (unsigned long long) n) ; + printf (" labeled nodes : %llu\n", (unsigned long long) nvals) ; + printf (" communities : %llu\n", (unsigned long long) n_communities) ; + printf (" modularity (gamma=1): %.12g\n", Q) ; + printf (" run time (sec) : %.6g\n", t) ; + + LG_FREE_ALL ; + LAGRAPH_TRY (LAGraph_Finalize (msg)) ; + return (GrB_SUCCESS) ; +#else + return (GrB_NOT_IMPLEMENTED) ; +#endif +} + diff --git a/experimental/test/test_leiden.c b/experimental/test/test_leiden.c index 2a429ecb81..7d4ca66a2a 100644 --- a/experimental/test/test_leiden.c +++ b/experimental/test/test_leiden.c @@ -6,6 +6,7 @@ // SPDX-License-Identifier: BSD-2-Clause #include +#include #include #include "GraphBLAS.h" @@ -22,20 +23,58 @@ char filename[LEN + 1] ; typedef struct { const char *matrix_file ; - bool expect_communities ; // true iff non-trivial Q is expected - double min_modularity ; // minimum acceptable modularity - double min_coverage ; // minimum intra-cluster edge ratio - double min_performance ; // minimum partition performance + LAGraph_Kind kind ; + bool force_symmetric ; // make A symmetric via A + A' + bool force_positive ; // make all edge weights positive + bool require_fp64 ; // assert resulting matrix is FP64 + double min_modularity ; // set < -1 to skip threshold } matrix_info ; const matrix_info files[] = { - // matrix_file, expect_communities, min_Q, min_coverage, min_performance - { "karate.mtx", true, 0.35, 0.15, 0.30 }, - { "comm0.mtx", true, 0.25, 0.10, 0.20 }, - { "", false, -1.0, -1.0, -1.0 } + // matrix_file, kind, force_sym, force_positive, require_fp64, min_Q + { "karate.mtx", LAGraph_ADJACENCY_UNDIRECTED, false, true, false, 0.35 }, + { "comm0.mtx", LAGraph_ADJACENCY_UNDIRECTED, false, true, false, 0.25 }, + { "west0067.mtx", LAGraph_ADJACENCY_DIRECTED, true, true, true, -2.0 }, + { "jagmesh7.mtx", LAGraph_ADJACENCY_UNDIRECTED, false, true, false, -2.0 }, + { "bcsstk13.mtx", LAGraph_ADJACENCY_UNDIRECTED, false, true, true, -2.0 }, + { "cryg2500.mtx", LAGraph_ADJACENCY_DIRECTED, true, true, true, -2.0 }, + { "", LAGraph_ADJACENCY_UNDIRECTED, false, false, false, -2.0 } } ; +const char *nonfinite_files[] = +{ + "matrix_fp64.mtx", + "skew_fp64.mtx", + "" +} ; + +static int check_all_finite (const GrB_Matrix A, bool *all_finite) +{ + GxB_Iterator it = NULL ; + GrB_Info info ; + *all_finite = true ; + info = GxB_Iterator_new (&it) ; + if (info != GrB_SUCCESS) return info ; + info = GxB_Matrix_Iterator_attach (it, A, NULL) ; + if (info != GrB_SUCCESS) goto done ; + info = GxB_Matrix_Iterator_seek (it, 0) ; + while (info == GrB_SUCCESS) + { + double aij = GxB_Iterator_get_FP64 (it) ; + if (!isfinite (aij)) + { + *all_finite = false ; + break ; + } + info = GxB_Matrix_Iterator_next (it) ; + } +done: + GrB_free (&it) ; + if (info == GxB_EXHAUSTED || info == GrB_SUCCESS) return GrB_SUCCESS ; + return info ; +} + //------------------------------------------------------------------------------ // test_Leiden //------------------------------------------------------------------------------ @@ -59,17 +98,55 @@ void test_Leiden (void) OK (LAGraph_MMRead (&A, f, msg)) ; fclose (f) ; - OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_UNDIRECTED, msg)) ; + OK (LAGraph_New (&G, &A, files[k].kind, msg)) ; TEST_CHECK (A == NULL) ; // LAGraph_New takes ownership - // Ensure symmetry cache is populated (required by some checks). - OK (LAGraph_Cached_IsSymmetricStructure (G, msg)) ; + // Ensure symmetry as needed for Leiden. + if (files[k].force_symmetric) + { + OK (LAGraph_Cached_AT (G, msg)) ; + OK (LAGraph_Cached_IsSymmetricStructure (G, msg)) ; + if (G->is_symmetric_structure == LAGraph_FALSE) + { + OK (GrB_eWiseAdd (G->A, NULL, NULL, GrB_ONEB_FP64, G->A, G->AT, NULL)) ; + G->is_symmetric_structure = LAGraph_TRUE ; + } + } + else + { + OK (LAGraph_Cached_IsSymmetricStructure (G, msg)) ; + } + TEST_CHECK (G->is_symmetric_structure == LAGraph_TRUE) ; + + if (files[k].force_positive) + { + OK (GrB_apply (G->A, NULL, NULL, GrB_ABS_FP64, G->A, NULL)) ; + } + + if (files[k].require_fp64) + { + GrB_Type atype = NULL ; + OK (GxB_Matrix_type (&atype, G->A)) ; + TEST_CHECK (atype == GrB_FP64) ; + } + + OK (LAGraph_Cached_EMin (G, msg)) ; + double min_val = 0.0 ; + OK (GrB_Scalar_extractElement_FP64 (&min_val, G->emin)) ; + TEST_CHECK (min_val >= 0.0) ; OK (LAGraph_Cached_OutDegree (G, msg)) ; uint64_t seed = 0 ; //unused GrB_Vector c = NULL ; - OK (LAGraph_Leiden (&c, G, seed, msg)) ; + GrB_Info info = LAGraph_Leiden (&c, G, seed, msg) ; + TEST_CHECK (info == GrB_SUCCESS) ; + if (info != GrB_SUCCESS) + { + GrB_free (&c) ; + OK (LAGraph_Delete (&G, msg)) ; + continue ; + } TEST_CHECK (c != NULL) ; // Every node must have a community label. @@ -91,78 +168,64 @@ void test_Leiden (void) TEST_CHECK (min_label >= 0) ; TEST_CHECK (max_label < n) ; - GrB_Index n_communities = max_label + 1 ; - // Compute modularity Q (requires SuiteSparse:GraphBLAS). double Q = 0.0 ; OK (LAGr_Modularity (&Q, 1.0, c, G, msg)) ; printf (" Modularity Q = %f\n", Q) ; - if (files[k].expect_communities) + if (files[k].min_modularity > -1.0) { // Validate modularity meets threshold for this graph TEST_CHECK (Q > files[k].min_modularity) ; TEST_MSG ("Expected Q > %f for %s, got Q = %f", files[k].min_modularity, aname, Q) ; } + TEST_CHECK (isfinite (Q)) ; + TEST_CHECK (Q >= -1.0 && Q <= 1.0) ; - // Compute per-community statistics for robustness validation - printf (" Community statistics:\n") ; - GrB_Index total_community_edges = 0 ; - GrB_Index total_edges = 0 ; - OK (GrB_Matrix_nvals (&total_edges, G->A)) ; + GrB_free (&c) ; + OK (LAGraph_Delete (&G, msg)) ; + } - for (GrB_Index com = 0 ; com < n_communities ; com++) - { - GrB_Index com_size = 0 ; - GrB_Index com_edges = 0 ; + LAGraph_Finalize (msg) ; +} - // Count nodes and edges in this community - for (GrB_Index i = 0 ; i < n ; i++) - { - int64_t label ; - GrB_Info info = GrB_Vector_extractElement_INT64 (&label, c, i) ; - if (info == GrB_SUCCESS && label == (int64_t)com) - { - com_size++ ; - // Count internal edges: edges from i to other nodes in same community - for (GrB_Index j = 0 ; j < n ; j++) - { - int64_t label_j ; - GrB_Info info_j = GrB_Vector_extractElement_INT64 (&label_j, c, j) ; - if (info_j == GrB_SUCCESS && label_j == (int64_t)com) - { - bool has_edge = false ; - GrB_Info info_e = GrB_Matrix_extractElement_BOOL (&has_edge, G->A, i, j) ; - if (info_e == GrB_SUCCESS && has_edge) - { - com_edges++ ; - } - } - } - } - } +void test_Leiden_NonfiniteInputs (void) +{ + LAGraph_Init (msg) ; - if (com_size > 0) - { - double density = (com_size > 1) ? - (double)com_edges / (com_size * (com_size - 1)) : 0.0 ; - printf (" Community %llu: size=%llu, edges=%llu, density=%f\n", - (unsigned long long) com, (unsigned long long) com_size, - (unsigned long long) com_edges, density) ; - total_community_edges += com_edges ; - } + for (int k = 0 ;; k++) + { + const char *aname = nonfinite_files [k] ; + if (strlen (aname) == 0) break ; + + printf ("\n====== nonfinite %s ======\n", aname) ; + snprintf (filename, LEN, LG_DATA_DIR "%s", aname) ; + + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + TEST_MSG ("Cannot open %s", filename) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + fclose (f) ; + + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; + OK (LAGraph_Cached_AT (G, msg)) ; + OK (LAGraph_Cached_IsSymmetricStructure (G, msg)) ; + if (G->is_symmetric_structure == LAGraph_FALSE) + { + OK (GrB_eWiseAdd (G->A, NULL, NULL, GrB_PLUS_FP64, G->A, G->AT, NULL)) ; + G->is_symmetric_structure = LAGraph_TRUE ; } - // Compute edge-cut ratio (edges crossing communities / total edges) - GrB_Index cut_edges = total_edges - total_community_edges ; - double cut_ratio = (total_edges > 0) ? - (double)cut_edges / total_edges : 0.0 ; - printf (" Edge-cut ratio = %f (cut=%llu, total=%llu)\n", - cut_ratio, (unsigned long long) cut_edges, - (unsigned long long) total_edges) ; + GrB_Type atype = NULL ; + OK (GxB_Matrix_type (&atype, G->A)) ; + TEST_CHECK (atype == GrB_FP64) ; + + bool all_finite = true ; + OK (check_all_finite (G->A, &all_finite)) ; + TEST_CHECK (!all_finite) ; - GrB_free (&c) ; OK (LAGraph_Delete (&G, msg)) ; } @@ -176,5 +239,6 @@ void test_Leiden (void) TEST_LIST = { { "Leiden", test_Leiden }, + { "Leiden nonfinite inputs", test_Leiden_NonfiniteInputs }, { NULL, NULL } } ; From 20b48f761dfa5bcf078207f2827e2a17a611143d Mon Sep 17 00:00:00 2001 From: GomezGab Date: Sat, 15 Aug 2026 22:56:34 -0500 Subject: [PATCH 07/13] improve preformance by bypassing costly ops --- .../algorithm/LAGr_ModularityMatrix.c | 27 +- experimental/algorithm/LAGraph_Leiden.c | 265 ++++++++++++++---- experimental/test/test_leiden.c | 13 +- experimental/test/test_modularity.c | 102 +++++++ 4 files changed, 333 insertions(+), 74 deletions(-) diff --git a/experimental/algorithm/LAGr_ModularityMatrix.c b/experimental/algorithm/LAGr_ModularityMatrix.c index d14f473865..af0580c847 100644 --- a/experimental/algorithm/LAGr_ModularityMatrix.c +++ b/experimental/algorithm/LAGr_ModularityMatrix.c @@ -36,6 +36,7 @@ { \ GrB_free(&k); \ GrB_free(&B); \ + GrB_free(&S_t); \ } #undef LG_FREE_ALL #define LG_FREE_ALL \ @@ -62,25 +63,28 @@ #endif int LAGr_AdjModularity( - double *Q, + double *Q, // TODO: scalar double gamma, - GrB_Matrix A, - GrB_Matrix S, + const GrB_Matrix A, + const GrB_Matrix S, char *msg) { #if LG_SUITESPARSE_GRAPHBLAS_V10_2 LG_CLEAR_MSG; char MATRIX_TYPE[LAGRAPH_MSG_LEN]; - // GrB_set(GrB_GLOBAL, false, GxB_BURBLE); + //TODO: add checks + GrB_set(GrB_GLOBAL, false, GxB_BURBLE); GrB_Index n; GrB_Vector k = NULL; GrB_Matrix B = NULL; + GrB_Matrix S_t = NULL; GRB_TRY(GrB_Matrix_nrows(&n, A)); GRB_TRY(GrB_Matrix_new(&B, GrB_FP64, n, n)); + GRB_TRY(GrB_Matrix_new(&S_t, GrB_BOOL, n, n)); GRB_TRY(GrB_Vector_new(&k, GrB_FP64, n)); double m = 0.0; @@ -93,32 +97,35 @@ int LAGr_AdjModularity( if (m == 0.0) { *Q = 0.0; - LG_FREE_ALL; + LG_FREE_WORK; return GrB_SUCCESS; } double inv_m = -gamma / (2.0 * m); - // sum degrees per community - GRB_TRY (GrB_vxm (k, NULL, NULL, GxB_PLUS_FIRST_FP64, k, S, NULL)) ; + GRB_TRY (GrB_transpose (S_t, NULL, NULL, S, NULL)) ; + // GRB_TRY (GrB_vxm (k, NULL, NULL, GxB_PLUS_FIRST_FP64, k, S, NULL)) ; + GRB_TRY (GrB_mxv (k, NULL, NULL, GxB_PLUS_SECOND_FP64, S_t, k, NULL)) ; // square - GRB_TRY (GrB_assign (k, NULL, GrB_TIMES_FP64, k, GrB_ALL, n, NULL)) ; + GRB_TRY (GrB_apply (k, NULL, NULL, GxB_POW_FP64, k, 2.0, NULL)) ; // sum all of the products // this computed \sum_{i,j} (k_{i}k_{j}\delta(\sigma_i\sigma_j)) GRB_TRY (GrB_reduce (&Q_, NULL, GrB_PLUS_MONOID_FP64, k, NULL)) ; Q_ *= inv_m; + GRB_TRY (GrB_free (&k)) ; // Count the number of edges per node in the cluster originating at a node // in the cluster - GRB_TRY (GrB_mxm (B, S, NULL, GxB_PLUS_FIRST_FP64, A, S, GrB_DESC_S)); + GRB_TRY (GrB_mxm (B, S, NULL, GxB_PLUS_FIRST_FP64, A, S_t, GrB_DESC_ST1)); // sum the intra-cluster edges per node for all nodes GRB_TRY (GrB_reduce (&Q_, GrB_PLUS_FP64, GrB_PLUS_MONOID_FP64, B, NULL)) ; + GrB_set(GrB_GLOBAL, false, GxB_BURBLE); Q_ *= -inv_m; *Q = Q_; - LG_FREE_ALL; + LG_FREE_WORK; return (GrB_SUCCESS); #else return (GrB_NOT_IMPLEMENTED); diff --git a/experimental/algorithm/LAGraph_Leiden.c b/experimental/algorithm/LAGraph_Leiden.c index db2b628257..1364be036a 100644 --- a/experimental/algorithm/LAGraph_Leiden.c +++ b/experimental/algorithm/LAGraph_Leiden.c @@ -23,6 +23,24 @@ #define LEIDEN_MAX_ITER 20 +#ifndef LG_LEIDEN_TIMING +#define LG_LEIDEN_TIMING 0 +#endif + +#if LG_LEIDEN_TIMING +#define LG_LEIDEN_TIC(t) double t = LAGraph_WallClockTime ( ) +#define LG_LEIDEN_ELAPSED(t) (LAGraph_WallClockTime ( ) - (t)) +#define LG_LEIDEN_PRINTF(...) printf (__VA_ARGS__) +#define LG_LEIDEN_BURBLE_ON GRB_TRY (LG_SET_BURBLE (true)) +#define LG_LEIDEN_BURBLE_OFF GRB_TRY (LG_SET_BURBLE (false)) +#else +#define LG_LEIDEN_TIC(t) +#define LG_LEIDEN_ELAPSED(t) (0.0) +#define LG_LEIDEN_PRINTF(...) +#define LG_LEIDEN_BURBLE_ON +#define LG_LEIDEN_BURBLE_OFF +#endif + //------------------------------------------------------------------------------ // The Leiden algorithm is a modularity-based community detection method that // guarantees well-connected communities by introducing a Refinement phase @@ -80,7 +98,6 @@ GrB_free (&it); \ GrB_free (&c_deg); \ GrB_free (&x); \ - GrB_free (&gain); \ } #undef LG_FREE_ALL @@ -95,6 +112,9 @@ int LG_Leiden_move_nodes // output: GrB_Matrix C, // community matrix uint64_t *community, // community array (inside of commmunity matrix) + uint64_t *nodes_popped_handle, // queue pops in this phase + uint64_t *nodes_evaluated_handle, // non-singleton nodes evaluated + uint64_t *nodes_moved_handle, // nodes moved to a new community // input: const GrB_Matrix A, // adjacency matrix const GrB_Vector deg, // degree vector @@ -104,21 +124,26 @@ int LG_Leiden_move_nodes char* msg ) { uint64_t node_id, com_id, n_deg, n; - GrB_Vector x = NULL, gain = NULL; + GrB_Vector x = NULL; GrB_Vector c_deg = NULL ; - GxB_Iterator it = NULL; + GxB_Iterator it = NULL, neighbor_it = NULL; GrB_Matrix_nrows (&n, A) ; uint64_t queue_head = 0, queue_tail = n, queue_size = n + 1; GRB_TRY (GxB_Iterator_new (&it)) ; GRB_TRY (GrB_Vector_new (&c_deg, GrB_FP64, n)) ; GRB_TRY (GrB_Vector_new (&x, GrB_FP64, n)) ; - GRB_TRY (GrB_Vector_new (&gain, GrB_FP64, n)) ; + GRB_TRY (GxB_Iterator_new (&neighbor_it)) ; + // FIXME: could modify A + GrB_set (A, GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT) ; + GRB_TRY (GxB_rowIterator_attach (neighbor_it, A, NULL)); - GRB_TRY (GrB_mxv(c_deg, NULL, NULL, GxB_PLUS_SECOND_FP64, C, deg, NULL)); + GRB_TRY (GrB_assign (c_deg, NULL, NULL, 0.0, GrB_ALL, n, NULL)) ; + GRB_TRY (GrB_vxm(c_deg, NULL, NULL, GxB_PLUS_FIRST_FP64, deg, C, NULL)); // TODO: decide if queue should be made inside this function. - int64_t count = 0; + uint64_t nodes_popped = 0, nodes_evaluated = 0, nodes_moved = 0; while (queue_head != queue_tail) { // while queue not empty + nodes_popped++; node_id = queue[queue_head]; queue_head = (queue_head + 1) % queue_size; enqueued [node_id] = false; @@ -126,12 +151,17 @@ int LG_Leiden_move_nodes uint64_t n_neighbors; com_id = community [node_id] ; GRB_TRY (GrB_Vector_clear (x)); - GRB_TRY (GrB_Vector_setElement_FP64 (x, 0.0, node_id)); - GRB_TRY (GrB_vxm (x, NULL, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, x, A, NULL)) ; + // GRB_TRY (GrB_Vector_setElement_FP64 (x, 0.0, node_id)) ; + // GRB_TRY (GrB_vxm (x, NULL, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, x, A, NULL)) ; + // TODO: deg and c_deg as C_arrays? + GRB_TRY (GrB_Vector_setElement_FP64 (x, 0.0, node_id)) ; + GRB_TRY (GrB_Col_extract (x, NULL, GrB_PLUS_FP64, A, GrB_ALL, n, node_id, GrB_DESC_T0)) ; + GRB_TRY (GrB_Vector_nvals (&n_neighbors, x)) ; if (n_neighbors == 1) continue; // skip singletons + nodes_evaluated++; // give gain the number of edges connecting x to community i - GRB_TRY (GrB_vxm (gain, NULL, NULL, GxB_PLUS_FIRST_FP64, x, C, NULL)) ; + GRB_TRY (GrB_vxm (x, NULL, NULL, GxB_PLUS_FIRST_FP64, x, C, NULL)) ; // momentarily remove node degree from the community total GRB_TRY (GrB_Vector_extractElement_FP64 (&n_deg, deg, node_id)); @@ -140,16 +170,19 @@ int LG_Leiden_move_nodes GRB_TRY (GrB_assign (c_deg, NULL, GrB_MINUS_FP64, n_deg, &com_id, 1, NULL)); // Calculate gain in each neighboring community - GRB_TRY (GrB_apply (gain, gain, GrB_PLUS_FP64, GrB_TIMES_FP64, c_deg, m_inv2 * n_deg, GrB_DESC_S)); uint64_t max_gain_c; double max_gain_val = -1e100; - GRB_TRY (GxB_Vector_Iterator_attach (it, gain, NULL)); + double y = m_inv2 * n_deg; + GRB_TRY (GxB_Vector_Iterator_attach (it, x, NULL)); GrB_Info info = GxB_Vector_Iterator_seek (it, 0); while (info == GrB_SUCCESS) { uint64_t c = GxB_Vector_Iterator_getIndex (it); - double g = GxB_Iterator_get_FP64 (it); - if (g > max_gain_val) { - max_gain_val = g; + double gain = GxB_Iterator_get_FP64 (it); + double com_deg; + GRB_TRY (GrB_Vector_extractElement_FP64 (&com_deg, c_deg, c)) ; + gain += y * com_deg ; + if (gain > max_gain_val) { + max_gain_val = gain; max_gain_c = c; } info = GxB_Vector_Iterator_next (it); @@ -160,23 +193,27 @@ int LG_Leiden_move_nodes GRB_TRY (GrB_assign (c_deg, NULL, GrB_PLUS_FP64, n_deg, &max_gain_c, 1, NULL)); if (max_gain_c == com_id) continue; // no change in community + nodes_moved++; // push every neighbor in a different community that is not already in // the queue to the back of the queue - GRB_TRY (GxB_Vector_Iterator_attach (it, x, NULL)); - info = GxB_Vector_Iterator_seek (it, 0); + info = GxB_rowIterator_seekRow (neighbor_it, node_id); while (info == GrB_SUCCESS) { - uint64_t neighbor_id = GxB_Vector_Iterator_getIndex (it); + uint64_t neighbor_id = GxB_rowIterator_getColIndex (neighbor_it); if (max_gain_c != community [neighbor_id] && !enqueued [neighbor_id]) { + // TODO: make a macro for enque deque queue[queue_tail] = neighbor_id; queue_tail = (queue_tail + 1) % queue_size; enqueued[neighbor_id] = true; } - info = GxB_Vector_Iterator_next (it); + info = GxB_rowIterator_nextCol (neighbor_it); } } + if (nodes_popped_handle != NULL) (*nodes_popped_handle) = nodes_popped; + if (nodes_evaluated_handle != NULL) (*nodes_evaluated_handle) = nodes_evaluated; + if (nodes_moved_handle != NULL) (*nodes_moved_handle) = nodes_moved; LG_FREE_WORK ; return GrB_SUCCESS ; } @@ -184,12 +221,12 @@ int LG_Leiden_move_nodes #undef LG_FREE_WORK #define LG_FREE_WORK \ { \ - GrB_free (&s_deg); \ GrB_free (&x_count); \ GrB_free (&k); \ GrB_free (&X); \ GrB_free (&C_t); \ GrB_free (&it); \ + LAGraph_Free ((void **) &s_deg, NULL) ; \ } #undef LG_FREE_ALL @@ -198,12 +235,28 @@ int LG_Leiden_move_nodes LG_FREE_WORK ; \ } +typedef struct { + double *s_deg; + double km_inv; +} leiden_ctx ; + +void LG_Leiden_gain +( + double *z, + const double *x, + GrB_Index i, + GrB_Index j, + const leiden_ctx *thunk +) { + *z = *x + thunk->km_inv * thunk->s_deg [j]; +} // helper function: phase 2 of leiden. Output already allocated int LG_Leiden_refinement ( // output: GrB_Matrix S, // sub-community matrix uint64_t *sub_com, // sub-community array (inside of S matrix) + uint64_t *nodes_evaluated_handle, // non-singleton nodes evaluated // input: const GrB_Matrix C, // community matrix const uint64_t *community, // community array (inside of C matrix) @@ -216,15 +269,21 @@ int LG_Leiden_refinement ) { uint64_t node_id, com_id, n_deg, n; GrB_Matrix X = NULL; - GrB_Vector s_deg = NULL; + double *s_deg = NULL; GrB_Vector k = NULL; GrB_Matrix C_t = NULL ; // transpose of C GrB_Vector x_count = NULL; GxB_Iterator it = NULL ; + GrB_IndexUnaryOp gain_op = NULL; + GrB_Type leiden_ctx_t = NULL; + GrB_Scalar ctx_s = NULL; GRB_TRY (GrB_Matrix_nrows (&n, A)) ; uint64_t queue_head = 0, queue_tail = n, queue_size = n + 1; + // TODO: make C_t compact, use an different queue per each community. + // This makes it easier to do the per-comunity masking. + // TODO: decide if queue should be made inside this function. GRB_TRY (GrB_Matrix_new (&C_t, GrB_FP64, n, n)) ; GRB_TRY (GrB_transpose (C_t, NULL, NULL, C, NULL)) ; @@ -233,29 +292,35 @@ int LG_Leiden_refinement GRB_TRY (GrB_mxm (C_t, C_t, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, C_t, A, GrB_DESC_ST1)) ; // Identify nodes with in-cluster connections (more than 0) - // TODO: use a smarter selection using c_deg GRB_TRY (GrB_Vector_new (&x_count, GrB_FP64, n)); GRB_TRY (GrB_assign (x_count, NULL, NULL, 0.0, GrB_ALL, n, NULL)) ; GRB_TRY (GrB_reduce (x_count, NULL, GrB_PLUS_FP64, GrB_PLUS_MONOID_FP64, C_t, GrB_DESC_T0)) ; uint64_t n_communities; + LG_TRY (LAGraph_Malloc ((void **) &s_deg, n, sizeof (uint64_t), msg)) ; + GRB_TRY (GrB_Vector_extractTuples_FP64 (NULL, s_deg, &n, deg)) ; + GRB_TRY (GrB_Type_new (&leiden_ctx_t, sizeof(leiden_ctx))) ; GRB_TRY (GrB_Vector_new (&k, GrB_FP64, n)) ; - GRB_TRY (GrB_Vector_new (&s_deg, GrB_FP64, n)) ; GRB_TRY (GrB_Matrix_new (&X, GrB_FP64, n, n)) ; + GRB_TRY (GrB_Scalar_new (&ctx_s, leiden_ctx_t)) ; + GRB_TRY (GrB_IndexUnaryOp_new ( + &gain_op, (GxB_index_unary_function) LG_Leiden_gain, + GrB_FP64, GrB_FP64, leiden_ctx_t)) ; GRB_TRY (GxB_Iterator_new(&it)) ; GRB_TRY (GrB_select (C_t, NULL, NULL, GrB_VALUEGT_FP64, C_t, 0.0, NULL)) ; + leiden_ctx ctx = {.s_deg = s_deg, .km_inv = 0.0} ; - // NOTE: I am wiring this so that it will be relatively easy to parallelize - // this phase in the FUTURE + uint64_t nodes_evaluated = 0; while (queue_head != queue_tail) { // while queue not empty double count; node_id = queue[queue_head]; queue_head = (queue_head + 1) % queue_size; GRB_TRY (GrB_Vector_extractElement_FP64 (&count, x_count, node_id)); if (count <= 0.0) continue; // skip singletons + nodes_evaluated++; double com_deg, n_deg; com_id = community [node_id] ; @@ -272,13 +337,15 @@ int LG_Leiden_refinement // momentarily remove node degree from the community total GRB_TRY (GrB_Vector_extractElement_FP64 (&n_deg, deg, node_id)); - // s_deg [sub_com_id] -= ndeg - GRB_TRY (GrB_assign (s_deg, NULL, GrB_MINUS_FP64, n_deg, &sub_com_id, 1, NULL)) ; + s_deg [sub_com_id] -= n_deg; - GRB_TRY (GrB_Vector_setElement_FP64(k, m_inv2 * n_deg, sub_com_id)) ; + ctx.km_inv = m_inv2 * n_deg; + GRB_TRY (GrB_Scalar_setElement_UDT (ctx_s, &leiden_ctx_t)) ; // Calculate gain in each neighboring community - GRB_TRY (GrB_mxm (X, X, GrB_PLUS_FP64, GrB_PLUS_TIMES_SEMIRING_FP64, - (GrB_Matrix) s_deg, (GrB_Matrix) k, GrB_DESC_ST1)) ; + // FIXME: killing preformance, this isn't using the mask. + GRB_TRY (GrB_apply (X, NULL, NULL, gain_op, X, ctx_s, NULL)) ; + // GRB_TRY (GrB_mxm (X, X, GrB_PLUS_FP64, GrB_PLUS_TIMES_SEMIRING_FP64, + // (GrB_Matrix) k, (GrB_Matrix) s_deg, GrB_DESC_ST1)) ; uint64_t new_c = com_id, row = 0; @@ -287,6 +354,7 @@ int LG_Leiden_refinement // Select from X (those with gain >= current and choose randomly) double current_gain, total_gain; + // FIXME: just do this inside the loop below. GRB_TRY (GrB_Matrix_extractElement_FP64 (¤t_gain, X, com_id, node_id)) ; GRB_TRY (GrB_apply (X, NULL, NULL, GrB_MINUS_FP64, X, current_gain, NULL)) ; GRB_TRY (GrB_select (X, NULL, NULL, GrB_VALUEGE_FP64, X, 0.0, NULL)) ; @@ -309,9 +377,9 @@ int LG_Leiden_refinement } sub_com [node_id] = new_c; - // s_deg [new_c] += ndeg - GRB_TRY (GrB_assign (s_deg, NULL, GrB_PLUS_FP64, n_deg, &new_c, 1, NULL)); + s_deg [new_c] += n_deg; } + if (nodes_evaluated_handle != NULL) (*nodes_evaluated_handle) = nodes_evaluated; LG_FREE_WORK ; return GrB_SUCCESS ; } @@ -348,6 +416,7 @@ int LG_Leiden_aggregate ) { uint64_t node_id, com_id, n_deg, queue_head = 0, queue_tail = 0; GrB_Vector s_list = NULL, x = NULL; + // TODO:rename S_new, maybe make it pure output instead of i/o GrB_Matrix S_squished = NULL, S_new = NULL ; GrB_Matrix C_squished = NULL, C_new = NULL ; GrB_Matrix A_squished = NULL, A_new = NULL ; @@ -407,10 +476,8 @@ int LG_Leiden_aggregate GrB_free (&C); \ GrB_free (&S); \ if (free_A) GrB_free (&A) ; \ - GrB_free (&mask); \ GrB_free (°); \ GrB_free (&c_deg); \ - GrB_free (&gain); \ GrB_free (&x); \ LAGraph_Free ((void **) &queue, NULL); \ LAGraph_Free ((void **) &enqueued, NULL); \ @@ -429,17 +496,17 @@ int LAGraph_Leiden // output: GrB_Vector *c_handle, // c[i] = community label (0..K-1) for node i // input: - LAGraph_Graph G, // input graph (must be symmetric, no self-loops) - uint64_t seed, // random seed (reserved; not yet used) + const LAGraph_Graph G, // input graph (must be symmetric, no self-loops, + // have positive numerical weights) + uint64_t seed, // random seed char *msg ) { -#ifdef LG_SUITESPARSE_GRAPHBLAS_V10 +#if LG_SUITESPARSE_GRAPHBLAS_V10 uint64_t *queue = NULL; // circular buffer with next nodes to look up bool *enqueued = NULL; // true if node is in the queue uint64_t *community = NULL ; // array: i in community c[i] GrB_Matrix C = NULL ; // contains readonly pointer to community array - GrB_Matrix mask = NULL ; // mask each community GxB_Container cont = NULL; uint64_t *sub_com = NULL ; // array: i in sub-community sub_com[i] @@ -449,7 +516,6 @@ int LAGraph_Leiden GrB_Vector c_deg = NULL ; // degree of each community GrB_Matrix A = NULL ; // G->A - GrB_Vector gain = NULL ; // gain from moving into a community GrB_Vector x = NULL; // full boolean vector GrB_Matrix node_map = NULL ; // maps aggregate node to nodes it contains GrB_Matrix new_node_map = NULL ; // maps aggregate node to nodes it contains @@ -464,9 +530,7 @@ int LAGraph_Leiden original_n = n; // Store original number of nodes // Input checking - LG_ASSERT_MSG (G->out_degree != NULL, LAGRAPH_NOT_CACHED, - "G->out_degree must be defined") ; - LG_ASSERT_MSG( + LG_ASSERT_MSG ( G->is_symmetric_structure == LAGraph_TRUE, LAGRAPH_SYMMETRIC_STRUCTURE_REQUIRED, "G->A must be symmetric") ; @@ -486,9 +550,7 @@ int LAGraph_Leiden // Initialize GraphBLAS data structures GRB_TRY (GxB_Iterator_new(&neighbor_it)); GRB_TRY (GxB_Container_new (&cont)) ; - GRB_TRY (GrB_Vector_new(&x, GrB_FP64, n)); GRB_TRY (GrB_Vector_new(°, GrB_FP64, n)); - GRB_TRY (GrB_Vector_new(&gain, GrB_FP64, n)); GRB_TRY (GrB_Vector_new(&x, GrB_BOOL, n)); GRB_TRY (GrB_assign (x, NULL, NULL, (bool) 0, GrB_ALL, n, NULL)); @@ -497,7 +559,7 @@ int LAGraph_Leiden // Initialize full degree vector from graph GRB_TRY (GrB_assign (deg, NULL, NULL, 0.0, GrB_ALL, n, NULL)); - GRB_TRY (GrB_reduce(deg, NULL, NULL, GrB_PLUS_MONOID_FP64, A, NULL)); + GRB_TRY (GrB_reduce (deg, NULL, GrB_PLUS_FP64, GrB_PLUS_MONOID_FP64, A, NULL)); // don't use out degree unless A is bool adj matrix // GRB_TRY (GrB_assign (deg, NULL, GrB_PLUS_FP64, G->out_degree, GrB_ALL, n, NULL)); @@ -506,19 +568,27 @@ int LAGraph_Leiden GRB_TRY (GrB_reduce (&m, NULL, GrB_PLUS_MONOID_FP64, deg, NULL)) ; LG_ASSERT_MSG (isfinite(m), GrB_INVALID_VALUE, "Matrix must reduce to a finite value") ; - - double m_inv2 = -1 / (2 * m); + double m_inv2 = -1.0 / (m); srand(seed); // Outer loop: repeat phases until convergence - double prev_modularity = 0; + double prev_modularity = -1; double modularity_epsilon = 1e-6; + double total_move_time = 0, total_refine_time = 0, total_agg_time = 0; + double total_modularity_time = 0, total_unpack_time = 0, total_queue_time = 0; + uint64_t total_move_evaluated = 0, total_move_popped = 0, total_move_moved = 0; + uint64_t total_refine_evaluated = 0; + for (int count = 0; count < LEIDEN_MAX_ITER; count++) { - printf ("Iteration %d, node count %lu \n", count, n) ; + double current_modularity; + LG_LEIDEN_PRINTF ("[Leiden timing] iteration %d, coarse nodes=%llu\n", + count, (unsigned long long) n) ; // Initialize queue in random order (seed-based) + LG_LEIDEN_TIC (t_queue_phase1); uint64_t queue_len = n; + // TODO: store values and memcopy in a different array GRB_TRY (GrB_Vector_extractTuples_FP64 (queue, NULL, &queue_len, deg)) ; // TODO: shuffle with LAGraph Random instead. @@ -531,10 +601,14 @@ int LAGraph_Leiden } } memset (enqueued, true, n) ; + double queue_phase1_time = LG_LEIDEN_ELAPSED (t_queue_phase1); + total_queue_time += queue_phase1_time; // get i vector from C and S to use as community and sub_com arrays // keep inside C and S as read-only - do { + LG_LEIDEN_TIC (t_unpack); + + { GrB_Type type; uint64_t n_community, X_memsize; int handling; @@ -553,23 +627,50 @@ int LAGraph_Leiden GRB_TRY (GrB_set (S, GxB_SPARSE, GxB_SPARSITY_CONTROL)) ; GRB_TRY (GrB_set (S, GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT)) ; GRB_TRY (GrB_set (S, 64, GxB_COLINDEX_INTEGER_HINT)) ; + GrB_wait (S, GrB_MATERIALIZE) ; GRB_TRY (GxB_unload_Matrix_into_Container (S, cont, NULL)) ; GRB_TRY (GxB_Vector_unload (cont->i, (void **) &sub_com, &type, &n_community, &X_memsize, &handling, NULL)) ; GRB_TRY (GxB_Vector_load (cont->i, (void **) &sub_com, type, n_community, X_memsize, GxB_IS_READONLY, NULL)) ; GRB_TRY (GxB_load_Matrix_from_Container(S, cont, NULL)) ; - } while (0) ; + } + + double unpack_time = LG_LEIDEN_ELAPSED (t_unpack); + total_unpack_time += unpack_time; //---------------------------------------------------------------------- // Phase 1: move nodes //---------------------------------------------------------------------- - LG_TRY (LG_Leiden_move_nodes (C, community, A, deg, queue, enqueued, m_inv2, msg)) ; + uint64_t move_popped = 0, move_evaluated = 0, move_moved = 0; + LG_LEIDEN_TIC (t_move); + // LG_LEIDEN_BURBLE_ON; + LG_TRY (LG_Leiden_move_nodes (C, community, &move_popped, &move_evaluated, + &move_moved, A, deg, queue, enqueued, m_inv2, msg)) ; + LG_LEIDEN_BURBLE_OFF; + double move_time = LG_LEIDEN_ELAPSED (t_move); + total_move_time += move_time; + total_move_popped += move_popped; + total_move_evaluated += move_evaluated; + total_move_moved += move_moved; + LG_LEIDEN_PRINTF ( + " move_nodes: t=%.6fs, queue_pops=%llu, evaluated=%llu, moved=%llu, amortized=%.3fus/evaluated\n", + move_time, + (unsigned long long) move_popped, + (unsigned long long) move_evaluated, + (unsigned long long) move_moved, + (move_evaluated == 0) ? 0.0 : (1e6 * move_time / (double) move_evaluated)) ; + + #ifndef NDEBUG + LG_TRY (LAGr_AdjModularity (¤t_modularity, 1.0, A, C, msg)); + ASSERT (current_modularity >= prev_modularity - modularity_epsilon) ; + #endif //---------------------------------------------------------------------- // Prep phase 2 //---------------------------------------------------------------------- // Initialize queue + LG_LEIDEN_TIC (t_queue_phase2); queue_len = n; GRB_TRY (GrB_Vector_extractTuples_FP64 (queue, NULL, &queue_len, deg)) ; @@ -582,20 +683,41 @@ int LAGraph_Leiden queue[j] = temp; } } + double queue_phase2_time = LG_LEIDEN_ELAPSED (t_queue_phase2); + total_queue_time += queue_phase2_time; //---------------------------------------------------------------------- // Phase 2: refine //---------------------------------------------------------------------- // Save Phase 1 communities as parent communities - LG_TRY (LG_Leiden_refinement (S, sub_com, C, community, c_deg, A, deg, queue, m_inv2, msg)) ; + uint64_t refine_evaluated = 0; + LG_LEIDEN_TIC (t_refine); + LG_LEIDEN_BURBLE_ON; + LG_TRY (LG_Leiden_refinement (S, sub_com, &refine_evaluated, C, community, + c_deg, A, deg, queue, m_inv2, msg)) ; + LG_LEIDEN_BURBLE_OFF; + double refine_time = LG_LEIDEN_ELAPSED (t_refine); + total_refine_time += refine_time; + total_refine_evaluated += refine_evaluated; + LG_LEIDEN_PRINTF ( + " refinement: t=%.6fs, evaluated=%llu, amortized=%.3fus/evaluated\n", + refine_time, + (unsigned long long) refine_evaluated, + (refine_evaluated == 0) ? 0.0 : (1e6 * refine_time / (double) refine_evaluated)) ; //---------------------------------------------------------------------- // Phase 3: aggregate //---------------------------------------------------------------------- // Compute modularity to check convergence - double current_modularity; + LG_LEIDEN_TIC (t_modularity); + // LG_LEIDEN_BURBLE_ON; LG_TRY (LAGr_AdjModularity (¤t_modularity, 1.0, A, S, msg)); + LG_LEIDEN_BURBLE_OFF; + double modularity_time = LG_LEIDEN_ELAPSED (t_modularity); + total_modularity_time += modularity_time; + LG_LEIDEN_PRINTF (" modularity: Q=%.12g, t=%.6fs\n", + current_modularity, modularity_time) ; ASSERT (current_modularity >= prev_modularity - modularity_epsilon) ; if (fabs(current_modularity - prev_modularity) < modularity_epsilon) { @@ -603,7 +725,13 @@ int LAGraph_Leiden } prev_modularity = current_modularity; + LG_LEIDEN_TIC (t_agg); + // LG_LEIDEN_BURBLE_ON; LG_TRY (LG_Leiden_aggregate (&S, &C, &A, free_A, msg)) ; + LG_LEIDEN_BURBLE_OFF; + double agg_time = LG_LEIDEN_ELAPSED (t_agg); + total_agg_time += agg_time; + LG_LEIDEN_PRINTF (" aggregate: t=%.6fs\n", agg_time) ; // free internal pointers that were not freed LG_TRY (LAGraph_Free ((void **) &community, msg)) ; LG_TRY (LAGraph_Free ((void **) &sub_com, msg)) ; @@ -635,6 +763,23 @@ int LAGraph_Leiden #endif } // end outer while loop + LG_LEIDEN_PRINTF ( + "[Leiden timing] totals: queue=%.6fs, unpack=%.6fs, move=%.6fs, refine=%.6fs, modularity=%.6fs, aggregate=%.6fs\n", + total_queue_time, total_unpack_time, total_move_time, total_refine_time, + total_modularity_time, total_agg_time) ; + LG_LEIDEN_PRINTF ( + "[Leiden timing] move totals: queue_pops=%llu, evaluated=%llu, moved=%llu, amortized=%.3fus/evaluated\n", + (unsigned long long) total_move_popped, + (unsigned long long) total_move_evaluated, + (unsigned long long) total_move_moved, + (total_move_evaluated == 0) ? 0.0 : + (1e6 * total_move_time / (double) total_move_evaluated)) ; + LG_LEIDEN_PRINTF ( + "[Leiden timing] refinement totals: evaluated=%llu, amortized=%.3fus/evaluated\n", + (unsigned long long) total_refine_evaluated, + (total_refine_evaluated == 0) ? 0.0 : + (1e6 * total_refine_time / (double) total_refine_evaluated)) ; + // Map refined communities back to original nodes GRB_TRY (GrB_Vector_new(c_handle, GrB_INT64, original_n)); @@ -644,14 +789,11 @@ int LAGraph_Leiden GRB_TRY (GxB_Vector_load (*c_handle, (void **) &community, GrB_INT64, n, sizeof (int64_t) * n, GrB_DEFAULT, NULL)); } else { - GRB_TRY (GxB_Matrix_Iterator_attach (neighbor_it, node_map, NULL)) ; - GrB_Info info = GxB_Matrix_Iterator_seek (neighbor_it, 0) ; - while (info == GrB_SUCCESS) { - uint64_t coarse_id, orig_id; - GxB_Matrix_Iterator_getIndex (neighbor_it, &coarse_id, &orig_id) ; - GRB_TRY (GrB_Vector_setElement_INT64 (*c_handle, community[coarse_id], orig_id)) ; - info = GxB_Matrix_Iterator_next (neighbor_it) ; - } + GRB_TRY (GrB_Vector_resize (x, n)) ; + GRB_TRY (GrB_assign ( + *c_handle, NULL, NULL, (int64_t) 0, GrB_ALL, original_n, NULL)) ; + GRB_TRY (GrB_vxm (*c_handle, NULL, GrB_PLUS_INT64, + GxB_PLUS_FIRSTJ_INT64, x, node_map, NULL)) ; } LG_FREE_WORK; @@ -660,4 +802,3 @@ int LAGraph_Leiden return GrB_NOT_IMPLEMENTED; #endif } - diff --git a/experimental/test/test_leiden.c b/experimental/test/test_leiden.c index 7d4ca66a2a..c845b4e9b3 100644 --- a/experimental/test/test_leiden.c +++ b/experimental/test/test_leiden.c @@ -17,6 +17,8 @@ char msg[LAGRAPH_MSG_LEN] ; LAGraph_Graph G = NULL ; GrB_Matrix A = NULL ; +GrB_Matrix C = NULL ; +GrB_Scalar zero_bool = NULL ; #define LEN 512 char filename[LEN + 1] ; @@ -168,9 +170,16 @@ void test_Leiden (void) TEST_CHECK (min_label >= 0) ; TEST_CHECK (max_label < n) ; - // Compute modularity Q (requires SuiteSparse:GraphBLAS). + // Compute modularity Q double Q = 0.0 ; - OK (LAGr_Modularity (&Q, 1.0, c, G, msg)) ; + GrB_Descriptor desc = NULL; + OK (GrB_Matrix_new (&C, GrB_BOOL, n, n)) ; + OK (GrB_Scalar_new (&zero_bool, GrB_BOOL)) ; + OK (GrB_Scalar_setElement_BOOL (zero_bool, false)) ; + OK (GrB_Descriptor_new (&desc)) ; + OK (GrB_set (desc, GxB_USE_INDICES, GxB_ROWINDEX_LIST)) ; + OK (GxB_Matrix_build_Scalar_Vector(C, c, c, zero_bool, desc)) ; + OK (LAGr_AdjModularity (&Q, 1.0, G->A, C, msg)) ; printf (" Modularity Q = %f\n", Q) ; if (files[k].min_modularity > -1.0) diff --git a/experimental/test/test_modularity.c b/experimental/test/test_modularity.c index 7a22a5e73f..01a379852b 100644 --- a/experimental/test/test_modularity.c +++ b/experimental/test/test_modularity.c @@ -26,6 +26,20 @@ const matrix_info files[] = { {"comm0.mtx", "comm0_Sa.mtx", 1, 0.35714}, // {"com-Amazon.mtx", "comm0_Sa.mtx",1, -1}, {"", "", -1, -1}}; + +typedef struct +{ + const char *matrix_file; + const char *cluster_file; + const double gamma; +} large_matrix_info; + +const large_matrix_info large_files[] = { + {"jagmesh7.mtx", "jagmesh7_cluster.mtx", 1}, + {"bcsstk13.mtx", "bcsstk13_cluster.mtx", 1}, + {"west0067.mtx", "west0067_cluster.mtx", 1}, + {"", "", -1}}; + void test_modularity(void) { #if LG_SUITESPARSE_GRAPHBLAS_V10_2 @@ -60,6 +74,94 @@ void test_modularity(void) OK(GrB_free(&S)); } + for (int k = 0;; k++) + { + if (strlen(large_files[k].matrix_file) == 0) + break; + + snprintf(filename, LEN, LG_DATA_DIR "%s", large_files[k].matrix_file); + FILE *f = fopen(filename, "r"); + TEST_CHECK(f != NULL); + OK(LAGraph_MMRead(&A, f, msg)); + OK(fclose(f)); + + GrB_Vector c = NULL; + snprintf(filename2, LEN, LG_DATA_DIR "%s", large_files[k].cluster_file); + FILE *t = fopen(filename2, "r"); + TEST_CHECK(t != NULL); + OK(LAGraph_MMRead((GrB_Matrix *)&c, t, msg)); + OK(fclose(t)); + + GrB_Index n = 0, csize = 0; + OK(GrB_Matrix_nrows(&n, A)); + OK(GrB_Vector_size(&csize, c)); + TEST_CHECK(n == csize); + + GrB_Index cnvals = 0; + OK(GrB_Vector_nvals(&cnvals, c)); + + GrB_Index *I = NULL; + int64_t *X = NULL; + int64_t *labels = NULL; + int64_t *unique_labels = NULL; + GrB_Index *comm_ids = NULL; + + OK(LAGraph_Malloc((void **)&I, cnvals, sizeof(GrB_Index), msg)); + OK(LAGraph_Malloc((void **)&X, cnvals, sizeof(int64_t), msg)); + OK(LAGraph_Calloc((void **)&labels, n, sizeof(int64_t), msg)); + OK(LAGraph_Malloc((void **)&unique_labels, n, sizeof(int64_t), msg)); + OK(LAGraph_Malloc((void **)&comm_ids, n, sizeof(GrB_Index), msg)); + + OK(GrB_Vector_extractTuples_INT64(I, X, &cnvals, c)); + for (GrB_Index p = 0; p < cnvals; p++) + { + labels[I[p]] = X[p]; + } + + GrB_Index ncomms = 0; + for (GrB_Index i = 0; i < n; i++) + { + int64_t label = labels[i]; + GrB_Index cidx = 0; + bool found = false; + for (GrB_Index j = 0; j < ncomms; j++) + { + if (unique_labels[j] == label) + { + cidx = j; + found = true; + break; + } + } + if (!found) + { + unique_labels[ncomms] = label; + cidx = ncomms; + ncomms++; + } + comm_ids[i] = cidx; + } + + OK(GrB_Matrix_new(&S, GrB_BOOL, n, n)); + for (GrB_Index i = 0; i < n; i++) + { + OK(GrB_Matrix_setElement_BOOL(S, true, i, comm_ids[i])); + } + + double Q; + OK(LAGr_AdjModularity(&Q, large_files[k].gamma, A, S, msg)); + TEST_CHECK(isfinite(Q)); + + OK(LAGraph_Free((void **)&I, msg)); + OK(LAGraph_Free((void **)&X, msg)); + OK(LAGraph_Free((void **)&labels, msg)); + OK(LAGraph_Free((void **)&unique_labels, msg)); + OK(LAGraph_Free((void **)&comm_ids, msg)); + OK(GrB_free(&c)); + OK(GrB_free(&A)); + OK(GrB_free(&S)); + } + OK(GrB_free(&A)); OK(GrB_free(&S)); From 733623ba85ad74b43a53ab1c6c1d18d39f5236e4 Mon Sep 17 00:00:00 2001 From: GomezGab Date: Wed, 26 Aug 2026 20:45:37 -0500 Subject: [PATCH 08/13] Change ModularityMatrix, update tests, and update leiden --- .../algorithm/LAGr_ModularityMatrix.c | 73 +++--- experimental/algorithm/LAGraph_Leiden.c | 231 +++++++++++++----- experimental/test/test_leiden.c | 54 +++- experimental/test/test_louvain.c | 4 +- experimental/test/test_modularity.c | 3 +- 5 files changed, 248 insertions(+), 117 deletions(-) diff --git a/experimental/algorithm/LAGr_ModularityMatrix.c b/experimental/algorithm/LAGr_ModularityMatrix.c index af0580c847..8387217d60 100644 --- a/experimental/algorithm/LAGr_ModularityMatrix.c +++ b/experimental/algorithm/LAGr_ModularityMatrix.c @@ -30,40 +30,22 @@ #include "GraphBLAS.h" #include "LG_internal.h" #include -#include -#define LG_FREE_MOD \ - { \ - GrB_free(&k); \ - GrB_free(&B); \ - GrB_free(&S_t); \ +#undef LG_FREE_WORK +#define LG_FREE_WORK \ + { \ + GrB_free (&k); \ + GrB_free (&B); \ + GrB_free (&S_t); \ } #undef LG_FREE_ALL #define LG_FREE_ALL \ { \ - LG_FREE_MOD; \ + LG_FREE_WORK; \ } -#define DEBUG 0 -#if DEBUG -#define check() printf("here") -#define dbg(x) \ - if (DEBUG) \ - GxB_print(x, 5) -#define err(x, info) \ - if (!(info == GrB_SUCCESS || info == GrB_NO_VALUE)) \ - { \ - char **err; \ - GrB_error(err, x); \ - printf("\ninfo: %d error: %s\n", info, err); \ - } -#else -#define check() -#define dbg(x) -#define err(x, info) -#endif -int LAGr_AdjModularity( - double *Q, // TODO: scalar +int LAGr_AdjModularity ( + double *Q, double gamma, const GrB_Matrix A, const GrB_Matrix S, @@ -71,21 +53,28 @@ int LAGr_AdjModularity( { #if LG_SUITESPARSE_GRAPHBLAS_V10_2 LG_CLEAR_MSG; - char MATRIX_TYPE[LAGRAPH_MSG_LEN]; - //TODO: add checks - GrB_set(GrB_GLOBAL, false, GxB_BURBLE); - - GrB_Index n; - GrB_Vector k = NULL; + GrB_Index nrows, ncols, nrows_s, ncols_s; + GrB_Vector k = NULL, com_deg = NULL; GrB_Matrix B = NULL; GrB_Matrix S_t = NULL; - GRB_TRY(GrB_Matrix_nrows(&n, A)); + LG_ASSERT (A != NULL, GrB_NULL_POINTER) ; + LG_ASSERT (S != NULL, GrB_NULL_POINTER) ; + + GxB_print(A, 2); + GRB_TRY (GrB_Matrix_nrows (&nrows, A)); + GRB_TRY (GrB_Matrix_ncols (&ncols, A)); + GRB_TRY (GrB_Matrix_nrows (&nrows_s, S)); + GRB_TRY (GrB_Matrix_ncols (&ncols_s, S)); + + LG_ASSERT (nrows == ncols, GrB_DIMENSION_MISMATCH) ; + LG_ASSERT (ncols == nrows_s, GrB_DIMENSION_MISMATCH) ; - GRB_TRY(GrB_Matrix_new(&B, GrB_FP64, n, n)); - GRB_TRY(GrB_Matrix_new(&S_t, GrB_BOOL, n, n)); - GRB_TRY(GrB_Vector_new(&k, GrB_FP64, n)); + GRB_TRY (GrB_Matrix_new (&B, GrB_FP64, nrows, ncols)); + GRB_TRY (GrB_Matrix_new (&S_t, GrB_BOOL, ncols_s, nrows_s)); + GRB_TRY (GrB_Vector_new (&k, GrB_FP64, ncols)); + GRB_TRY (GrB_Vector_new (&com_deg, GrB_FP64, ncols_s)); double m = 0.0; double Q_ = 0.0; @@ -103,16 +92,15 @@ int LAGr_AdjModularity( double inv_m = -gamma / (2.0 * m); // sum degrees per community GRB_TRY (GrB_transpose (S_t, NULL, NULL, S, NULL)) ; - // GRB_TRY (GrB_vxm (k, NULL, NULL, GxB_PLUS_FIRST_FP64, k, S, NULL)) ; - GRB_TRY (GrB_mxv (k, NULL, NULL, GxB_PLUS_SECOND_FP64, S_t, k, NULL)) ; + GRB_TRY (GrB_mxv (com_deg, NULL, NULL, GxB_PLUS_SECOND_FP64, S_t, k, NULL)) ; // square - GRB_TRY (GrB_apply (k, NULL, NULL, GxB_POW_FP64, k, 2.0, NULL)) ; + GRB_TRY (GrB_apply (com_deg, NULL, NULL, GxB_POW_FP64, com_deg, 2.0, NULL)) ; // sum all of the products // this computed \sum_{i,j} (k_{i}k_{j}\delta(\sigma_i\sigma_j)) - GRB_TRY (GrB_reduce (&Q_, NULL, GrB_PLUS_MONOID_FP64, k, NULL)) ; + GRB_TRY (GrB_reduce (&Q_, NULL, GrB_PLUS_MONOID_FP64, com_deg, NULL)) ; Q_ *= inv_m; - GRB_TRY (GrB_free (&k)) ; + GRB_TRY (GrB_free (&com_deg)) ; // Count the number of edges per node in the cluster originating at a node // in the cluster @@ -120,7 +108,6 @@ int LAGr_AdjModularity( // sum the intra-cluster edges per node for all nodes GRB_TRY (GrB_reduce (&Q_, GrB_PLUS_FP64, GrB_PLUS_MONOID_FP64, B, NULL)) ; - GrB_set(GrB_GLOBAL, false, GxB_BURBLE); Q_ *= -inv_m; *Q = Q_; diff --git a/experimental/algorithm/LAGraph_Leiden.c b/experimental/algorithm/LAGraph_Leiden.c index 1364be036a..4f81f9263e 100644 --- a/experimental/algorithm/LAGraph_Leiden.c +++ b/experimental/algorithm/LAGraph_Leiden.c @@ -24,14 +24,14 @@ #define LEIDEN_MAX_ITER 20 #ifndef LG_LEIDEN_TIMING -#define LG_LEIDEN_TIMING 0 +#define LG_LEIDEN_TIMING 1 #endif #if LG_LEIDEN_TIMING #define LG_LEIDEN_TIC(t) double t = LAGraph_WallClockTime ( ) #define LG_LEIDEN_ELAPSED(t) (LAGraph_WallClockTime ( ) - (t)) #define LG_LEIDEN_PRINTF(...) printf (__VA_ARGS__) -#define LG_LEIDEN_BURBLE_ON GRB_TRY (LG_SET_BURBLE (true)) +#define LG_LEIDEN_BURBLE_ON GRB_TRY (LG_SET_BURBLE (false)) #define LG_LEIDEN_BURBLE_OFF GRB_TRY (LG_SET_BURBLE (false)) #else #define LG_LEIDEN_TIC(t) @@ -41,6 +41,18 @@ #define LG_LEIDEN_BURBLE_OFF #endif +#define LG_QUEUE_ENQUEUE(queue, queue_tail, queue_size, value) \ +{ \ + (queue) [queue_tail] = (value); \ + (queue_tail) = ((queue_tail) + 1) % (queue_size); \ +} + +#define LG_QUEUE_DEQUEUE(queue, queue_head, queue_size, value) \ +{ \ + (value) = (queue) [queue_head]; \ + (queue_head) = ((queue_head) + 1) % (queue_size); \ +} + //------------------------------------------------------------------------------ // The Leiden algorithm is a modularity-based community detection method that // guarantees well-connected communities by introducing a Refinement phase @@ -95,6 +107,7 @@ #undef LG_FREE_WORK #define LG_FREE_WORK \ { \ + GrB_free (&neighbor_it); \ GrB_free (&it); \ GrB_free (&c_deg); \ GrB_free (&x); \ @@ -127,15 +140,13 @@ int LG_Leiden_move_nodes GrB_Vector x = NULL; GrB_Vector c_deg = NULL ; GxB_Iterator it = NULL, neighbor_it = NULL; - GrB_Matrix_nrows (&n, A) ; + GRB_TRY (GrB_Matrix_nrows (&n, A)) ; uint64_t queue_head = 0, queue_tail = n, queue_size = n + 1; GRB_TRY (GxB_Iterator_new (&it)) ; GRB_TRY (GrB_Vector_new (&c_deg, GrB_FP64, n)) ; GRB_TRY (GrB_Vector_new (&x, GrB_FP64, n)) ; GRB_TRY (GxB_Iterator_new (&neighbor_it)) ; - // FIXME: could modify A - GrB_set (A, GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT) ; GRB_TRY (GxB_rowIterator_attach (neighbor_it, A, NULL)); GRB_TRY (GrB_assign (c_deg, NULL, NULL, 0.0, GrB_ALL, n, NULL)) ; @@ -144,8 +155,7 @@ int LG_Leiden_move_nodes uint64_t nodes_popped = 0, nodes_evaluated = 0, nodes_moved = 0; while (queue_head != queue_tail) { // while queue not empty nodes_popped++; - node_id = queue[queue_head]; - queue_head = (queue_head + 1) % queue_size; + LG_QUEUE_DEQUEUE (queue, queue_head, queue_size, node_id) ; enqueued [node_id] = false; double com_deg, n_deg; uint64_t n_neighbors; @@ -202,9 +212,7 @@ int LG_Leiden_move_nodes uint64_t neighbor_id = GxB_rowIterator_getColIndex (neighbor_it); if (max_gain_c != community [neighbor_id] && !enqueued [neighbor_id]) { - // TODO: make a macro for enque deque - queue[queue_tail] = neighbor_id; - queue_tail = (queue_tail + 1) % queue_size; + LG_QUEUE_ENQUEUE (queue, queue_tail, queue_size, neighbor_id) ; enqueued[neighbor_id] = true; } @@ -221,12 +229,14 @@ int LG_Leiden_move_nodes #undef LG_FREE_WORK #define LG_FREE_WORK \ { \ - GrB_free (&x_count); \ - GrB_free (&k); \ + GrB_free (&gain_op); \ + GrB_free (&ctx_s); \ + GrB_free (&leiden_ctx_t); \ GrB_free (&X); \ GrB_free (&C_t); \ GrB_free (&it); \ LAGraph_Free ((void **) &s_deg, NULL) ; \ + LAGraph_Free ((void **) &is_rep, NULL) ; \ } #undef LG_FREE_ALL @@ -270,9 +280,8 @@ int LG_Leiden_refinement uint64_t node_id, com_id, n_deg, n; GrB_Matrix X = NULL; double *s_deg = NULL; - GrB_Vector k = NULL; + bool *is_rep = NULL; GrB_Matrix C_t = NULL ; // transpose of C - GrB_Vector x_count = NULL; GxB_Iterator it = NULL ; GrB_IndexUnaryOp gain_op = NULL; GrB_Type leiden_ctx_t = NULL; @@ -285,41 +294,32 @@ int LG_Leiden_refinement // This makes it easier to do the per-comunity masking. // TODO: decide if queue should be made inside this function. - GRB_TRY (GrB_Matrix_new (&C_t, GrB_FP64, n, n)) ; + GRB_TRY (GrB_Matrix_new (&C_t, GrB_BOOL, n, n)) ; GRB_TRY (GrB_transpose (C_t, NULL, NULL, C, NULL)) ; - // Find number of in-cluster connections for each node in the cluster - GRB_TRY (GrB_mxm (C_t, C_t, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, C_t, A, GrB_DESC_ST1)) ; - - // Identify nodes with in-cluster connections (more than 0) - GRB_TRY (GrB_Vector_new (&x_count, GrB_FP64, n)); - GRB_TRY (GrB_assign (x_count, NULL, NULL, 0.0, GrB_ALL, n, NULL)) ; - GRB_TRY (GrB_reduce (x_count, NULL, GrB_PLUS_FP64, GrB_PLUS_MONOID_FP64, C_t, GrB_DESC_T0)) ; - uint64_t n_communities; + LG_TRY (LAGraph_Malloc ((void **) &is_rep, n, sizeof (bool), msg)) ; + memset (is_rep, 0, n) ; LG_TRY (LAGraph_Malloc ((void **) &s_deg, n, sizeof (uint64_t), msg)) ; GRB_TRY (GrB_Vector_extractTuples_FP64 (NULL, s_deg, &n, deg)) ; GRB_TRY (GrB_Type_new (&leiden_ctx_t, sizeof(leiden_ctx))) ; - GRB_TRY (GrB_Vector_new (&k, GrB_FP64, n)) ; GRB_TRY (GrB_Matrix_new (&X, GrB_FP64, n, n)) ; GRB_TRY (GrB_Scalar_new (&ctx_s, leiden_ctx_t)) ; GRB_TRY (GrB_IndexUnaryOp_new ( &gain_op, (GxB_index_unary_function) LG_Leiden_gain, GrB_FP64, GrB_FP64, leiden_ctx_t)) ; GRB_TRY (GxB_Iterator_new(&it)) ; - - GRB_TRY (GrB_select (C_t, NULL, NULL, GrB_VALUEGT_FP64, C_t, 0.0, NULL)) ; leiden_ctx ctx = {.s_deg = s_deg, .km_inv = 0.0} ; - uint64_t nodes_evaluated = 0; while (queue_head != queue_tail) { // while queue not empty double count; - node_id = queue[queue_head]; - queue_head = (queue_head + 1) % queue_size; - GRB_TRY (GrB_Vector_extractElement_FP64 (&count, x_count, node_id)); - if (count <= 0.0) continue; // skip singletons + LG_QUEUE_DEQUEUE (queue, queue_head, queue_size, node_id) ; + + if (is_rep [node_id]) { + continue; + } nodes_evaluated++; double com_deg, n_deg; @@ -330,9 +330,9 @@ int LG_Leiden_refinement // nodes via their communities. GRB_TRY (GrB_Matrix_clear (X)); GRB_TRY (GrB_Matrix_setElement_FP64 (X, 0.0, com_id, node_id)); - GRB_TRY (GrB_mxm (X, C_t, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, X, A, GrB_DESC_S)) ; + GRB_TRY (GrB_mxm (X, C_t, NULL, GxB_PLUS_SECOND_FP64, X, A, GrB_DESC_S)) ; // give X the number of edges connecting x to community i - GRB_TRY (GrB_mxm (X, NULL, NULL, GxB_PLUS_SECOND_FP64, X, S, NULL)) ; + GRB_TRY (GrB_mxm (X, NULL, NULL, GxB_PLUS_FIRST_FP64, X, S, NULL)) ; // momentarily remove node degree from the community total GRB_TRY (GrB_Vector_extractElement_FP64 (&n_deg, deg, node_id)); @@ -340,42 +340,28 @@ int LG_Leiden_refinement s_deg [sub_com_id] -= n_deg; ctx.km_inv = m_inv2 * n_deg; - GRB_TRY (GrB_Scalar_setElement_UDT (ctx_s, &leiden_ctx_t)) ; - // Calculate gain in each neighboring community - // FIXME: killing preformance, this isn't using the mask. + GRB_TRY (GrB_Scalar_setElement_UDT (ctx_s, &ctx)) ; + // Calculate gain in each neighboring community. GRB_TRY (GrB_apply (X, NULL, NULL, gain_op, X, ctx_s, NULL)) ; - // GRB_TRY (GrB_mxm (X, X, GrB_PLUS_FP64, GrB_PLUS_TIMES_SEMIRING_FP64, - // (GrB_Matrix) k, (GrB_Matrix) s_deg, GrB_DESC_ST1)) ; - - uint64_t new_c = com_id, row = 0; // FUTURE: This part will have to change if we want parrallel clusters // to work - - // Select from X (those with gain >= current and choose randomly) - double current_gain, total_gain; - // FIXME: just do this inside the loop below. - GRB_TRY (GrB_Matrix_extractElement_FP64 (¤t_gain, X, com_id, node_id)) ; - GRB_TRY (GrB_apply (X, NULL, NULL, GrB_MINUS_FP64, X, current_gain, NULL)) ; - GRB_TRY (GrB_select (X, NULL, NULL, GrB_VALUEGE_FP64, X, 0.0, NULL)) ; - GRB_TRY (GrB_apply (X, NULL, NULL, GxB_POW_FP64, M_E, X, NULL)) ; - GRB_TRY (GrB_reduce (&total_gain, NULL, GrB_PLUS_MONOID_FP64, X, NULL)) ; - double rand_f = drand48 () * total_gain ; - + uint64_t new_c = sub_com_id, row = 0; + // Select the best gain candidate (temporary greedy refinement). GRB_TRY (GxB_Matrix_Iterator_attach(it, X, NULL)) ; - GrB_Info info = GxB_Matrix_Iterator_seek (it, 0); - total_gain = 0; + GrB_Info info = GxB_Matrix_Iterator_seek (it, 0) ; + double best_gain = 0.0 ; while (info == GrB_SUCCESS) { - total_gain += GxB_Iterator_get_FP64 (it) ; - - if (total_gain > rand_f) { + double gain = GxB_Iterator_get_FP64 (it) ; + if (gain > best_gain) { + best_gain = gain ; GxB_Matrix_Iterator_getIndex(it, &row, &new_c) ; - break; } info = GxB_Matrix_Iterator_next (it); } + is_rep [new_c] = true; sub_com [node_id] = new_c; s_deg [new_c] += n_deg; } @@ -384,6 +370,90 @@ int LG_Leiden_refinement return GrB_SUCCESS ; } +#undef LG_FREE_WORK +#define LG_FREE_WORK \ +{ \ + GrB_free (&row_parent_count); \ + GrB_free (&row_split_count); \ + GrB_free (&all_ones); \ + GrB_free (&C_tS_in_C_t); \ + GrB_free (&C_tS); \ + GrB_free (&C_t); \ +} + +#undef LG_FREE_ALL +#define LG_FREE_ALL \ +{ \ + LG_FREE_WORK ; \ +} + +// helper function: phase 2 metrics (split counts + parent containment check) +int LG_Leiden_refinement_stats +( + // output: + uint64_t *split_total_handle, // total extra splits beyond 1 per parent + double *split_avg_handle, // average extra splits per parent + bool *within_parent_handle, // true iff C_tS is a submatrix of C_t + // input: + const GrB_Matrix C, // parent-community matrix + const GrB_Matrix S, // refined sub-community matrix + char *msg +) +{ + (void) msg; + GrB_Index n; + GrB_Matrix C_t = NULL, C_tS = NULL, C_tS_in_C_t = NULL; + GrB_Vector all_ones = NULL, row_split_count = NULL, row_parent_count = NULL; + + GRB_TRY (GrB_Matrix_nrows (&n, C)) ; + GRB_TRY (GrB_Matrix_new (&C_t, GrB_BOOL, n, n)) ; + GRB_TRY (GrB_transpose (C_t, NULL, NULL, C, NULL)) ; + GRB_TRY (GrB_assign ( + C_t, C_t, NULL, (bool) true, GrB_ALL, n, GrB_ALL, n, GrB_DESC_S)) ; + + GRB_TRY (GrB_Matrix_new (&C_tS, GrB_BOOL, n, n)) ; + GRB_TRY (GrB_mxm (C_tS, NULL, NULL, LAGraph_any_one_bool, C_t, S, NULL)) ; + GRB_TRY (GrB_assign ( + C_tS, C_tS, NULL, (bool) true, GrB_ALL, n, GrB_ALL, n, GrB_DESC_S)) ; + + GRB_TRY (GrB_Vector_new (&all_ones, GrB_BOOL, n)) ; + GRB_TRY (GrB_assign (all_ones, NULL, NULL, (bool) true, GrB_ALL, n, NULL)) ; + GRB_TRY (GrB_Vector_new (&row_split_count, GrB_INT32, n)) ; + GRB_TRY (GrB_Vector_new (&row_parent_count, GrB_INT32, n)) ; + GRB_TRY (GrB_mxv ( + row_split_count, NULL, NULL, LAGraph_plus_one_int32, C_tS, all_ones, NULL)) ; + GRB_TRY (GrB_mxv ( + row_parent_count, NULL, NULL, LAGraph_plus_one_int32, C_t, all_ones, NULL)) ; + + int64_t split_total_raw = 0; + GRB_TRY (GrB_reduce ( + &split_total_raw, NULL, GrB_PLUS_MONOID_INT64, row_split_count, NULL)) ; + GrB_Index n_parent = 0; + GRB_TRY (GrB_Vector_nvals (&n_parent, row_parent_count)) ; + int64_t split_total = split_total_raw - (int64_t) n_parent; + if (split_total < 0) split_total = 0; + + if (split_total_handle != NULL) (*split_total_handle) = (uint64_t) split_total; + if (split_avg_handle != NULL) { + (*split_avg_handle) = + (n_parent == 0) ? 0.0 : ((double) split_total / (double) n_parent); + } + + // Verify that every (parent, sub-community) incidence remains in C_t. + GRB_TRY (GrB_Matrix_new (&C_tS_in_C_t, GrB_BOOL, n, n)) ; + GRB_TRY (GrB_eWiseMult ( + C_tS_in_C_t, NULL, NULL, GrB_LAND, C_tS, C_t, NULL)) ; + GrB_Index nvals_C_tS = 0, nvals_intersection = 0; + GRB_TRY (GrB_Matrix_nvals (&nvals_C_tS, C_tS)) ; + GRB_TRY (GrB_Matrix_nvals (&nvals_intersection, C_tS_in_C_t)) ; + if (within_parent_handle != NULL) { + (*within_parent_handle) = (nvals_C_tS == nvals_intersection); + } + + LG_FREE_WORK ; + return GrB_SUCCESS ; +} + #undef LG_FREE_WORK #define LG_FREE_WORK \ { \ @@ -439,6 +509,7 @@ int LG_Leiden_aggregate GRB_TRY (GrB_Matrix_new (&S_new, GrB_BOOL, n_new, n)) ; GRB_TRY (GrB_Matrix_new (&C_new, GrB_BOOL, n_new, n_new)) ; GRB_TRY (GrB_Matrix_new (&A_new, GrB_FP64, n_new, n_new)) ; + GrB_set (A_new, GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT) ; GRB_TRY (GrB_Matrix_new (&S_squished, GrB_BOOL, n, n_new)) ; GRB_TRY (GrB_Matrix_new (&C_squished, GrB_BOOL, n, n_new)) ; @@ -473,12 +544,17 @@ int LG_Leiden_aggregate #undef LG_FREE_WORK #define LG_FREE_WORK \ { \ + GrB_free (&neighbor_it); \ + GrB_free (&cont); \ + GrB_free (&C_to_orig); \ GrB_free (&C); \ GrB_free (&S); \ if (free_A) GrB_free (&A) ; \ GrB_free (°); \ GrB_free (&c_deg); \ GrB_free (&x); \ + GrB_free (&node_map); \ + GrB_free (&new_node_map); \ LAGraph_Free ((void **) &queue, NULL); \ LAGraph_Free ((void **) &enqueued, NULL); \ LAGraph_Free ((void **) &community, NULL); \ @@ -519,6 +595,7 @@ int LAGraph_Leiden GrB_Vector x = NULL; // full boolean vector GrB_Matrix node_map = NULL ; // maps aggregate node to nodes it contains GrB_Matrix new_node_map = NULL ; // maps aggregate node to nodes it contains + GrB_Matrix C_to_orig = NULL ; // maps final C communities to original nodes uint64_t n, original_n; @@ -571,6 +648,7 @@ int LAGraph_Leiden double m_inv2 = -1.0 / (m); srand(seed); + srand48 ((long int) seed) ; // Outer loop: repeat phases until convergence double prev_modularity = -1; @@ -580,6 +658,8 @@ int LAGraph_Leiden double total_modularity_time = 0, total_unpack_time = 0, total_queue_time = 0; uint64_t total_move_evaluated = 0, total_move_popped = 0, total_move_moved = 0; uint64_t total_refine_evaluated = 0; + uint64_t total_refine_splits = 0, total_refine_calls = 0; + double total_refine_split_avg = 0.0; for (int count = 0; count < LEIDEN_MAX_ITER; count++) { double current_modularity; @@ -664,6 +744,7 @@ int LAGraph_Leiden #ifndef NDEBUG LG_TRY (LAGr_AdjModularity (¤t_modularity, 1.0, A, C, msg)); ASSERT (current_modularity >= prev_modularity - modularity_epsilon) ; + LG_LEIDEN_PRINTF ( " modularity %f\n", current_modularity) ; #endif //---------------------------------------------------------------------- @@ -691,18 +772,31 @@ int LAGraph_Leiden //---------------------------------------------------------------------- // Save Phase 1 communities as parent communities uint64_t refine_evaluated = 0; + uint64_t refine_split_total = 0; + double refine_split_avg = 0.0; + bool refine_within_parent = false; LG_LEIDEN_TIC (t_refine); LG_LEIDEN_BURBLE_ON; + // memcpy (sub_com, community, n * sizeof (uint64_t)) ; LG_TRY (LG_Leiden_refinement (S, sub_com, &refine_evaluated, C, community, - c_deg, A, deg, queue, m_inv2, msg)) ; + c_deg, A, deg, queue, m_inv2, msg)) ; + LG_TRY (LG_Leiden_refinement_stats ( + &refine_split_total, &refine_split_avg, &refine_within_parent, + C, S, msg)) ; LG_LEIDEN_BURBLE_OFF; double refine_time = LG_LEIDEN_ELAPSED (t_refine); total_refine_time += refine_time; total_refine_evaluated += refine_evaluated; + total_refine_splits += refine_split_total; + total_refine_split_avg += refine_split_avg; + total_refine_calls++; LG_LEIDEN_PRINTF ( - " refinement: t=%.6fs, evaluated=%llu, amortized=%.3fus/evaluated\n", + " refinement: t=%.6fs, evaluated=%llu, splits_total=%llu, splits_avg=%.6g, within_parent=%s, amortized=%.3fus/evaluated\n", refine_time, (unsigned long long) refine_evaluated, + (unsigned long long) refine_split_total, + refine_split_avg, + refine_within_parent ? "true" : "false", (refine_evaluated == 0) ? 0.0 : (1e6 * refine_time / (double) refine_evaluated)) ; //---------------------------------------------------------------------- @@ -712,7 +806,7 @@ int LAGraph_Leiden // Compute modularity to check convergence LG_LEIDEN_TIC (t_modularity); // LG_LEIDEN_BURBLE_ON; - LG_TRY (LAGr_AdjModularity (¤t_modularity, 1.0, A, S, msg)); + LG_TRY (LAGr_AdjModularity (¤t_modularity, 1.0, A, C, msg)); LG_LEIDEN_BURBLE_OFF; double modularity_time = LG_LEIDEN_ELAPSED (t_modularity); total_modularity_time += modularity_time; @@ -758,8 +852,8 @@ int LAGraph_Leiden GRB_TRY (GrB_Vector_resize (x, n)) ; GRB_TRY (GrB_Matrix_diag (&S, x, 0)) ; #ifndef NDEBUG - LG_TRY (LAGr_AdjModularity (¤t_modularity, 1.0, A, S, msg)); - ASSERT (fabs (current_modularity - prev_modularity) < modularity_epsilon) ; + LG_TRY (LAGr_AdjModularity (¤t_modularity, 1.0, A, C, msg)); + ASSERT (isfinite (current_modularity)) ; #endif } // end outer while loop @@ -775,25 +869,30 @@ int LAGraph_Leiden (total_move_evaluated == 0) ? 0.0 : (1e6 * total_move_time / (double) total_move_evaluated)) ; LG_LEIDEN_PRINTF ( - "[Leiden timing] refinement totals: evaluated=%llu, amortized=%.3fus/evaluated\n", + "[Leiden timing] refinement totals: evaluated=%llu, splits_total=%llu, splits_avg=%.6g, amortized=%.3fus/evaluated\n", (unsigned long long) total_refine_evaluated, + (unsigned long long) total_refine_splits, + (total_refine_calls == 0) ? 0.0 : + (total_refine_split_avg / (double) total_refine_calls), (total_refine_evaluated == 0) ? 0.0 : (1e6 * total_refine_time / (double) total_refine_evaluated)) ; - // Map refined communities back to original nodes + // Map communities in C back to original nodes GRB_TRY (GrB_Vector_new(c_handle, GrB_INT64, original_n)); - GRB_TRY (GrB_free (&C)) ; if (node_map == NULL) { ASSERT (n == original_n) ; GRB_TRY (GxB_Vector_load (*c_handle, (void **) &community, GrB_INT64, n, sizeof (int64_t) * n, GrB_DEFAULT, NULL)); } else { + GRB_TRY (GrB_Matrix_new (&C_to_orig, GrB_BOOL, n, original_n)) ; + GRB_TRY (GrB_mxm ( + C_to_orig, NULL, NULL, GxB_ANY_PAIR_BOOL, C, node_map, GrB_DESC_T0)) ; GRB_TRY (GrB_Vector_resize (x, n)) ; GRB_TRY (GrB_assign ( *c_handle, NULL, NULL, (int64_t) 0, GrB_ALL, original_n, NULL)) ; GRB_TRY (GrB_vxm (*c_handle, NULL, GrB_PLUS_INT64, - GxB_PLUS_FIRSTJ_INT64, x, node_map, NULL)) ; + GxB_PLUS_FIRSTJ_INT64, x, C_to_orig, NULL)) ; } LG_FREE_WORK; diff --git a/experimental/test/test_leiden.c b/experimental/test/test_leiden.c index c845b4e9b3..72a2c319b0 100644 --- a/experimental/test/test_leiden.c +++ b/experimental/test/test_leiden.c @@ -17,8 +17,6 @@ char msg[LAGRAPH_MSG_LEN] ; LAGraph_Graph G = NULL ; GrB_Matrix A = NULL ; -GrB_Matrix C = NULL ; -GrB_Scalar zero_bool = NULL ; #define LEN 512 char filename[LEN + 1] ; @@ -172,13 +170,15 @@ void test_Leiden (void) // Compute modularity Q double Q = 0.0 ; - GrB_Descriptor desc = NULL; + GrB_Matrix C = NULL ; + GrB_Scalar zero_bool = NULL ; + GrB_Descriptor desc = NULL ; OK (GrB_Matrix_new (&C, GrB_BOOL, n, n)) ; OK (GrB_Scalar_new (&zero_bool, GrB_BOOL)) ; OK (GrB_Scalar_setElement_BOOL (zero_bool, false)) ; OK (GrB_Descriptor_new (&desc)) ; OK (GrB_set (desc, GxB_USE_INDICES, GxB_ROWINDEX_LIST)) ; - OK (GxB_Matrix_build_Scalar_Vector(C, c, c, zero_bool, desc)) ; + OK (GxB_Matrix_build_Scalar_Vector (C, c, c, zero_bool, desc)) ; OK (LAGr_AdjModularity (&Q, 1.0, G->A, C, msg)) ; printf (" Modularity Q = %f\n", Q) ; @@ -192,6 +192,9 @@ void test_Leiden (void) TEST_CHECK (isfinite (Q)) ; TEST_CHECK (Q >= -1.0 && Q <= 1.0) ; + GrB_free (&desc) ; + GrB_free (&zero_bool) ; + GrB_free (&C) ; GrB_free (&c) ; OK (LAGraph_Delete (&G, msg)) ; } @@ -241,6 +244,46 @@ void test_Leiden_NonfiniteInputs (void) LAGraph_Finalize (msg) ; } +#if LG_BRUTAL_TESTS +void test_Leiden_brutal (void) +{ + OK (LG_brutal_setup (msg)) ; + OK (GxB_Global_Option_set (GxB_JIT_C_CONTROL, GxB_JIT_OFF)) ; + + snprintf (filename, LEN, LG_DATA_DIR "%s", "comm0.mtx") ; + uint64_t seed = 0 ; + GrB_Vector c = NULL ; + LAGraph_Graph H = NULL ; + GrB_Matrix B = NULL ; + + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + TEST_MSG ("Cannot open %s", filename) ; + OK (LAGraph_MMRead (&B, f, msg)) ; + fclose (f) ; + OK (LAGraph_New (&H, &B, LAGraph_ADJACENCY_UNDIRECTED, msg)) ; + OK (LAGraph_Cached_IsSymmetricStructure (H, msg)) ; + TEST_CHECK (H->is_symmetric_structure == LAGraph_TRUE) ; + OK (GrB_apply (H->A, NULL, NULL, GrB_ABS_FP64, H->A, NULL)) ; + OK (LAGraph_Cached_EMin (H, msg)) ; + + LG_brutal = INT64_MAX ; + OK (LAGraph_Leiden (&c, H, seed, msg)) ; + LG_brutal = -1 ; + + GrB_Index n, nvals ; + OK (GrB_Matrix_nrows (&n, H->A)) ; + OK (GrB_Vector_nvals (&nvals, c)) ; + TEST_CHECK (nvals == n) ; + + GrB_free (&c) ; + OK (LAGraph_Delete (&H, msg)) ; + GrB_free (&B) ; + + OK (LG_brutal_teardown (msg)) ; +} +#endif + //------------------------------------------------------------------------------ // test list //------------------------------------------------------------------------------ @@ -249,5 +292,8 @@ TEST_LIST = { { "Leiden", test_Leiden }, { "Leiden nonfinite inputs", test_Leiden_NonfiniteInputs }, + #if LG_BRUTAL_TESTS + { "Leiden brutal", test_Leiden_brutal }, + #endif { NULL, NULL } } ; diff --git a/experimental/test/test_louvain.c b/experimental/test/test_louvain.c index 6fc4a1638b..21ffef646b 100644 --- a/experimental/test/test_louvain.c +++ b/experimental/test/test_louvain.c @@ -25,10 +25,10 @@ typedef struct } matrix_info; const matrix_info files[] = { - {"empty.mtx",0}, {"comm0.mtx", 0.357142857142857}, {"karate.mtx", .42}, - {"", -1}}; + {"", -1} +}; void test_LouvainSeq(void) diff --git a/experimental/test/test_modularity.c b/experimental/test/test_modularity.c index 01a379852b..d4c67856a4 100644 --- a/experimental/test/test_modularity.c +++ b/experimental/test/test_modularity.c @@ -20,8 +20,7 @@ typedef struct const double mod; // expected modularity from Q = 1/2m (S^TBS) } matrix_info; -const matrix_info files[] = { - {"empty.mtx","empty.mtx",1,0}, +const matrix_info files[] = { {"comm0.mtx", "comm0_S.mtx", 1, -0.17347}, {"comm0.mtx", "comm0_Sa.mtx", 1, 0.35714}, // {"com-Amazon.mtx", "comm0_Sa.mtx",1, -1}, From 2ec46f3cd63b32c1b23bb56bef02b056832f5af3 Mon Sep 17 00:00:00 2001 From: GomezGab Date: Sat, 29 Aug 2026 16:29:25 -0500 Subject: [PATCH 09/13] Allow splitting in Leiden Phase 1 --- experimental/algorithm/LAGraph_Leiden.c | 133 ++++++++++++++++-------- experimental/test/test_leiden.c | 60 +++++++++++ 2 files changed, 150 insertions(+), 43 deletions(-) diff --git a/experimental/algorithm/LAGraph_Leiden.c b/experimental/algorithm/LAGraph_Leiden.c index 4f81f9263e..050049d012 100644 --- a/experimental/algorithm/LAGraph_Leiden.c +++ b/experimental/algorithm/LAGraph_Leiden.c @@ -107,10 +107,15 @@ #undef LG_FREE_WORK #define LG_FREE_WORK \ { \ + GrB_free (&c_deg_vec); \ + GrB_free (&c_size_vec); \ GrB_free (&neighbor_it); \ GrB_free (&it); \ - GrB_free (&c_deg); \ GrB_free (&x); \ + LAGraph_Free ((void **) &c_deg, NULL) ; \ + LAGraph_Free ((void **) &node_deg, NULL) ; \ + LAGraph_Free ((void **) &community_size, NULL) ; \ + LAGraph_Free ((void **) &empty_stack, NULL) ; \ } #undef LG_FREE_ALL @@ -136,30 +141,62 @@ int LG_Leiden_move_nodes double m_inv2, // -1 / (2 * m) char* msg ) { - uint64_t node_id, com_id, n_deg, n; - GrB_Vector x = NULL; - GrB_Vector c_deg = NULL ; - GxB_Iterator it = NULL, neighbor_it = NULL; + uint64_t node_id, com_id, n ; + GrB_Vector x = NULL ; + GrB_Vector c_deg_vec = NULL ; + GrB_Vector c_size_vec = NULL ; + double *c_deg = NULL ; + double *node_deg = NULL ; + uint64_t *community_size = NULL ; + uint64_t *empty_stack = NULL ; + uint64_t empty_top = 0 ; + GxB_Iterator it = NULL, neighbor_it = NULL ; GRB_TRY (GrB_Matrix_nrows (&n, A)) ; uint64_t queue_head = 0, queue_tail = n, queue_size = n + 1; GRB_TRY (GxB_Iterator_new (&it)) ; - GRB_TRY (GrB_Vector_new (&c_deg, GrB_FP64, n)) ; GRB_TRY (GrB_Vector_new (&x, GrB_FP64, n)) ; + GRB_TRY (GrB_Vector_new (&c_deg_vec, GrB_FP64, n)) ; + GRB_TRY (GrB_Vector_new (&c_size_vec, GrB_UINT64, n)) ; GRB_TRY (GxB_Iterator_new (&neighbor_it)) ; - GRB_TRY (GxB_rowIterator_attach (neighbor_it, A, NULL)); + GRB_TRY (GxB_rowIterator_attach (neighbor_it, A, NULL)) ; + LG_TRY (LAGraph_Malloc ((void **) &empty_stack, n, sizeof (uint64_t), msg)) ; + + // initialize degree and community degree vectors, then unpack to C arrays + GRB_TRY (GrB_assign (c_deg_vec, NULL, NULL, 0.0, GrB_ALL, n, NULL)) ; + GRB_TRY (GrB_assign (c_size_vec, NULL, NULL, (uint64_t) 0, GrB_ALL, n, NULL)) ; + GRB_TRY (GrB_vxm (c_deg_vec, NULL, GrB_PLUS_FP64, GxB_PLUS_FIRST_FP64, deg, C, NULL)) ; + GRB_TRY (GrB_vxm (c_size_vec, NULL, GrB_PLUS_UINT64, GxB_PLUS_PAIR_UINT64, deg, C, NULL)) ; + + GrB_Type node_type = NULL, cdeg_type = NULL, csize_type = NULL ; + uint64_t node_nheld = 0, cdeg_nheld = 0, csize_nheld = 0 ; + uint64_t node_size = 0, cdeg_size = 0, csize_size = 0 ; + int node_handling = 0, cdeg_handling = 0, csize_handling = 0 ; + GRB_TRY (GxB_Vector_unload (deg, (void **) &node_deg, &node_type, + &node_nheld, &node_size, &node_handling, NULL)) ; + GRB_TRY (GxB_Vector_unload (c_deg_vec, (void **) &c_deg, &cdeg_type, + &cdeg_nheld, &cdeg_size, &cdeg_handling, NULL)) ; + GRB_TRY (GxB_Vector_unload (c_size_vec, (void **) &community_size, &csize_type, + &csize_nheld, &csize_size, &csize_handling, NULL)) ; + + for (uint64_t c = 0 ; c < n ; c++) + { + if (community_size [c] == 0) + { + empty_stack [empty_top++] = c ; + } + } - GRB_TRY (GrB_assign (c_deg, NULL, NULL, 0.0, GrB_ALL, n, NULL)) ; - GRB_TRY (GrB_vxm(c_deg, NULL, NULL, GxB_PLUS_FIRST_FP64, deg, C, NULL)); // TODO: decide if queue should be made inside this function. uint64_t nodes_popped = 0, nodes_evaluated = 0, nodes_moved = 0; while (queue_head != queue_tail) { // while queue not empty nodes_popped++; LG_QUEUE_DEQUEUE (queue, queue_head, queue_size, node_id) ; - enqueued [node_id] = false; - double com_deg, n_deg; + enqueued [node_id] = false ; uint64_t n_neighbors; com_id = community [node_id] ; + const uint64_t old_size = community_size [com_id] ; + const double degree_i = node_deg [node_id] ; GRB_TRY (GrB_Vector_clear (x)); // GRB_TRY (GrB_Vector_setElement_FP64 (x, 0.0, node_id)) ; // GRB_TRY (GrB_vxm (x, NULL, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, x, A, NULL)) ; @@ -174,51 +211,65 @@ int LG_Leiden_move_nodes GRB_TRY (GrB_vxm (x, NULL, NULL, GxB_PLUS_FIRST_FP64, x, C, NULL)) ; // momentarily remove node degree from the community total - GRB_TRY (GrB_Vector_extractElement_FP64 (&n_deg, deg, node_id)); - - // c_deg [com_id] -= ndeg - GRB_TRY (GrB_assign (c_deg, NULL, GrB_MINUS_FP64, n_deg, &com_id, 1, NULL)); + c_deg [com_id] -= degree_i ; + community_size [com_id]-- ; // Calculate gain in each neighboring community - uint64_t max_gain_c; - double max_gain_val = -1e100; - double y = m_inv2 * n_deg; - GRB_TRY (GxB_Vector_Iterator_attach (it, x, NULL)); - GrB_Info info = GxB_Vector_Iterator_seek (it, 0); + uint64_t max_gain_c = com_id ; + double max_gain_val = -1.0 ; + double y = m_inv2 * degree_i ; + GRB_TRY (GxB_Vector_Iterator_attach (it, x, NULL)) ; + GrB_Info info = GxB_Vector_Iterator_seek (it, 0) ; while (info == GrB_SUCCESS) { - uint64_t c = GxB_Vector_Iterator_getIndex (it); - double gain = GxB_Iterator_get_FP64 (it); - double com_deg; - GRB_TRY (GrB_Vector_extractElement_FP64 (&com_deg, c_deg, c)) ; - gain += y * com_deg ; + uint64_t c = GxB_Vector_Iterator_getIndex (it) ; + double gain = GxB_Iterator_get_FP64 (it) ; + gain += y * c_deg [c] ; if (gain > max_gain_val) { - max_gain_val = gain; - max_gain_c = c; + max_gain_val = gain ; + max_gain_c = c ; } - info = GxB_Vector_Iterator_next (it); + info = GxB_Vector_Iterator_next (it) ; + } + + uint64_t new_com_id = com_id ; + if (max_gain_val >= 0.0) + { + new_com_id = max_gain_c ; + } + else if (old_size > 1 && empty_top > 0) + { + new_com_id = empty_stack [--empty_top] ; } - community [node_id] = max_gain_c; - // c_deg [node_id] += ndeg - GRB_TRY (GrB_assign (c_deg, NULL, GrB_PLUS_FP64, n_deg, &max_gain_c, 1, NULL)); + community [node_id] = new_com_id ; + c_deg [new_com_id] += degree_i ; + community_size [new_com_id]++ ; - if (max_gain_c == com_id) continue; // no change in community - nodes_moved++; + if (new_com_id == com_id) continue; // no change in community + nodes_moved++ ; + if (community_size [com_id] == 0) + { + empty_stack [empty_top++] = com_id ; + } // push every neighbor in a different community that is not already in // the queue to the back of the queue - info = GxB_rowIterator_seekRow (neighbor_it, node_id); + info = GxB_rowIterator_seekRow (neighbor_it, node_id) ; while (info == GrB_SUCCESS) { - uint64_t neighbor_id = GxB_rowIterator_getColIndex (neighbor_it); - if (max_gain_c != community [neighbor_id] + uint64_t neighbor_id = GxB_rowIterator_getColIndex (neighbor_it) ; + if (new_com_id != community [neighbor_id] && !enqueued [neighbor_id]) { LG_QUEUE_ENQUEUE (queue, queue_tail, queue_size, neighbor_id) ; - enqueued[neighbor_id] = true; + enqueued[neighbor_id] = true ; } - info = GxB_rowIterator_nextCol (neighbor_it); + info = GxB_rowIterator_nextCol (neighbor_it) ; } } + + GRB_TRY (GxB_Vector_load (deg, (void **) &node_deg, node_type, + node_nheld, node_size, node_handling, NULL)) ; + if (nodes_popped_handle != NULL) (*nodes_popped_handle) = nodes_popped; if (nodes_evaluated_handle != NULL) (*nodes_evaluated_handle) = nodes_evaluated; if (nodes_moved_handle != NULL) (*nodes_moved_handle) = nodes_moved; @@ -270,7 +321,6 @@ int LG_Leiden_refinement // input: const GrB_Matrix C, // community matrix const uint64_t *community, // community array (inside of C matrix) - const GrB_Vector c_deg, // degree vector const GrB_Matrix A, // adjacency matrix const GrB_Vector deg, // degree vector uint64_t *queue, // queue to use (contains all nodes) @@ -551,7 +601,6 @@ int LG_Leiden_aggregate GrB_free (&S); \ if (free_A) GrB_free (&A) ; \ GrB_free (°); \ - GrB_free (&c_deg); \ GrB_free (&x); \ GrB_free (&node_map); \ GrB_free (&new_node_map); \ @@ -589,8 +638,6 @@ int LAGraph_Leiden GrB_Matrix S = NULL ; // contains readonly pointer to sub-community array GrB_Vector deg = NULL ; // degree of each node - GrB_Vector c_deg = NULL ; // degree of each community - GrB_Matrix A = NULL ; // G->A GrB_Vector x = NULL; // full boolean vector GrB_Matrix node_map = NULL ; // maps aggregate node to nodes it contains @@ -779,7 +826,7 @@ int LAGraph_Leiden LG_LEIDEN_BURBLE_ON; // memcpy (sub_com, community, n * sizeof (uint64_t)) ; LG_TRY (LG_Leiden_refinement (S, sub_com, &refine_evaluated, C, community, - c_deg, A, deg, queue, m_inv2, msg)) ; + A, deg, queue, m_inv2, msg)) ; LG_TRY (LG_Leiden_refinement_stats ( &refine_split_total, &refine_split_avg, &refine_within_parent, C, S, msg)) ; diff --git a/experimental/test/test_leiden.c b/experimental/test/test_leiden.c index 72a2c319b0..e3fe699959 100644 --- a/experimental/test/test_leiden.c +++ b/experimental/test/test_leiden.c @@ -20,6 +20,21 @@ GrB_Matrix A = NULL ; #define LEN 512 char filename[LEN + 1] ; +int LG_Leiden_move_nodes +( + GrB_Matrix C, + uint64_t *community, + uint64_t *nodes_popped_handle, + uint64_t *nodes_evaluated_handle, + uint64_t *nodes_moved_handle, + const GrB_Matrix A, + const GrB_Vector deg, + uint64_t *queue, + bool *enqueued, + double m_inv2, + char *msg +) ; + typedef struct { const char *matrix_file ; @@ -141,6 +156,7 @@ void test_Leiden (void) GrB_Info info = LAGraph_Leiden (&c, G, seed, msg) ; TEST_CHECK (info == GrB_SUCCESS) ; + if (info != GrB_SUCCESS) { GrB_free (&c) ; @@ -244,6 +260,49 @@ void test_Leiden_NonfiniteInputs (void) LAGraph_Finalize (msg) ; } +void test_Leiden_Phase1SingletonSplit (void) +{ + LAGraph_Init (msg) ; + + const GrB_Index n = 4 ; + GrB_Matrix B = NULL, C = NULL ; + GrB_Vector deg = NULL ; + OK (GrB_Matrix_new (&B, GrB_FP64, n, n)) ; + OK (GrB_Matrix_new (&C, GrB_BOOL, n, n)) ; + OK (GrB_Vector_new (°, GrB_FP64, n)) ; + + GrB_Index Ai [4] = { 0, 2, 0, 3 } ; + GrB_Index Aj [4] = { 2, 0, 3, 0 } ; + double Ax [4] = { 1.0, 1.0, 1.0, 1.0 } ; + OK (GrB_Matrix_build_FP64 (B, Ai, Aj, Ax, 4, GrB_PLUS_FP64)) ; + + uint64_t community [4] = { 0, 0, 2, 3 } ; + for (GrB_Index i = 0 ; i < n ; i++) + { + OK (GrB_Matrix_setElement_BOOL (C, true, i, community [i])) ; + } + + // Deliberately chosen to force negative neighbor gains for node 0. + OK (GrB_Vector_setElement_FP64 (deg, 100.0, 0)) ; + OK (GrB_Vector_setElement_FP64 (deg, 100.0, 1)) ; + OK (GrB_Vector_setElement_FP64 (deg, 100.0, 2)) ; + OK (GrB_Vector_setElement_FP64 (deg, 100.0, 3)) ; + + uint64_t queue [5] = { 0, 1, 2, 3, 0 } ; + bool enqueued [4] = { true, true, true, true } ; + uint64_t popped = 0, evaluated = 0, moved = 0 ; + OK (LG_Leiden_move_nodes (C, community, &popped, &evaluated, &moved, + B, deg, queue, enqueued, -1.0 / 400.0, msg)) ; + + // Community 1 starts empty; node 0 should split into it. + TEST_CHECK (community [0] == 1) ; + + OK (GrB_free (°)) ; + OK (GrB_free (&C)) ; + OK (GrB_free (&B)) ; + LAGraph_Finalize (msg) ; +} + #if LG_BRUTAL_TESTS void test_Leiden_brutal (void) { @@ -292,6 +351,7 @@ TEST_LIST = { { "Leiden", test_Leiden }, { "Leiden nonfinite inputs", test_Leiden_NonfiniteInputs }, + { "Leiden phase1 singleton split", test_Leiden_Phase1SingletonSplit }, #if LG_BRUTAL_TESTS { "Leiden brutal", test_Leiden_brutal }, #endif From f81b4025a923b308f0f36ebee349dfecb67eb4b9 Mon Sep 17 00:00:00 2001 From: GomezGab Date: Sat, 29 Aug 2026 22:10:17 -0500 Subject: [PATCH 10/13] work around 10.3.0 crash --- experimental/algorithm/LAGraph_Leiden.c | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/experimental/algorithm/LAGraph_Leiden.c b/experimental/algorithm/LAGraph_Leiden.c index 050049d012..f37b583169 100644 --- a/experimental/algorithm/LAGraph_Leiden.c +++ b/experimental/algorithm/LAGraph_Leiden.c @@ -355,6 +355,14 @@ int LG_Leiden_refinement GRB_TRY (GrB_Vector_extractTuples_FP64 (NULL, s_deg, &n, deg)) ; GRB_TRY (GrB_Type_new (&leiden_ctx_t, sizeof(leiden_ctx))) ; GRB_TRY (GrB_Matrix_new (&X, GrB_FP64, n, n)) ; + // Work around SuiteSparse:GraphBLAS bug (67) present in v10.3.0 and + // fixed in v10.3.1: incorrect JIT kernel for R=masker(C,M,Z) when R is + // hypersparse. Apply only to affected versions by disabling hypersparsity + // for X on this path. + #if LAGRAPH_SUITESPARSE && (GxB_IMPLEMENTATION < GxB_VERSION (10,3,1)) + GRB_TRY (GxB_Matrix_Option_set_FP64 (X, GxB_HYPER_SWITCH, 0.0)) ; + #endif + GRB_TRY (GrB_set (X, GxB_SPARSE, GxB_SPARSITY_CONTROL)) ; GRB_TRY (GrB_Scalar_new (&ctx_s, leiden_ctx_t)) ; GRB_TRY (GrB_IndexUnaryOp_new ( &gain_op, (GxB_index_unary_function) LG_Leiden_gain, From 80efae4e7e30707f9e09adb439a7843ea0405f30 Mon Sep 17 00:00:00 2001 From: Roi Lipman Date: Wed, 26 Aug 2026 10:27:27 +0300 Subject: [PATCH 11/13] [wip] louvain --- experimental/algorithm/LAGraph_louvain.c | 360 +++++++++++++++++++++++ experimental/test/test_louvain.c | 280 +++++++++--------- experimental/test/test_louvain_seqis.c | 152 ++++++++++ 3 files changed, 661 insertions(+), 131 deletions(-) create mode 100644 experimental/algorithm/LAGraph_louvain.c create mode 100644 experimental/test/test_louvain_seqis.c diff --git a/experimental/algorithm/LAGraph_louvain.c b/experimental/algorithm/LAGraph_louvain.c new file mode 100644 index 0000000000..cdc51bf6be --- /dev/null +++ b/experimental/algorithm/LAGraph_louvain.c @@ -0,0 +1,360 @@ +//------------------------------------------------------------------------------ +// LAGraph_louvain.c: Louvain method +//------------------------------------------------------------------------------ + +// LAGraph, (c) 2019-2026 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 Roi Lipman and Gabriel Gomez, FalkorDB + +//------------------------------------------------------------------------------ + +#include +#include +#include + +#undef LG_FREE_ALL +#define LG_FREE_ALL \ +{ \ + GrB_free(&V); \ + GrB_free(&cont); \ +} + +// construct C, C[i:] = [i] +// and N, N[i] = i's community id +static GrB_Info _initCommunities +( + GrB_Matrix *C, // C [i] = i + uint64_t **I, // C's indices array + GrB_Index node_count // node count +) +{ + //-------------------------------------------------------------------------- + // initialize C + //-------------------------------------------------------------------------- + + // create full vector V[i] = 1 + char *msg = NULL ; + GrB_Vector V = NULL ; + GxB_Container cont = NULL ; + + GRB_TRY (GrB_Vector_new (&V, GrB_BOOL, node_count)) ; + GRB_TRY (GrB_assign (V, NULL, NULL, true, GrB_ALL, node_count, NULL)) ; + + // create diagonal matrix C from V + GRB_TRY (GrB_Matrix_diag (C, V, 0)) ; + GRB_TRY (GrB_free (&V)) ; + + // C should be CSR with 64 bit indicies + GRB_TRY (GrB_set (*C, 64, GxB_COLINDEX_INTEGER_HINT )) ; + GRB_TRY (GrB_set (*C, GxB_SPARSE, GxB_SPARSITY_CONTROL )) ; + GRB_TRY (GrB_set (*C, GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT)) ; + + //-------------------------------------------------------------------------- + // get a handle to C's indicies array + //-------------------------------------------------------------------------- + + // unload C into a container + GRB_TRY (GxB_Container_new (&cont)) ; + + // GrB_Matrix -> GxB_Container + GRB_TRY (GxB_unload_Matrix_into_Container (*C, cont, NULL)) ; + + int handling ; + GrB_Type type ; + uint64_t n, X_memsize ; + + // unload I + GRB_TRY (GxB_Vector_unload (cont->i, (void**)I, &type, &n, &X_memsize, + &handling, NULL)) ; + + // load I, mark it as READONLY + GRB_TRY (GxB_Vector_load (cont->i, (void **)I, type, n, X_memsize, + GxB_IS_READONLY, NULL)) ; + + // GxB_Container -> GrB_Matrix + GRB_TRY (GxB_load_Matrix_from_Container (*C, cont, NULL)) ; + GRB_TRY (GrB_free (&cont)) ; + + LG_FREE_ALL ; + + return GrB_SUCCESS ; +} + +#undef LG_FREE_WORK +#define LG_FREE_WORK \ +{ \ + GrB_free(&C); \ + GrB_free(&D); \ + GrB_free(&it); \ + GrB_free(&Ni); \ + GrB_free(&desc); \ + LAGraph_Free((void**)°ree, NULL); \ + LAGraph_Free((void**)&community_degree, NULL); \ +} + +#undef LG_FREE_ALL +#define LG_FREE_ALL \ +{ \ + LG_FREE_WORK ; \ +} + +// computes the modularity contribution of attaching a node of degree +// i_degree, with kin edges into a community, to that community -- where +// sigma_tot is the total degree of the community's members, *excluding* +// the node itself if it happens to already be a member (the caller is +// responsible for that adjustment; see base_gain) +static inline double gain +( + uint64_t i_degree, // node i's degree + uint64_t kjn, // number of edges connecting node i to community j + uint64_t sigma_tot, // total degree of community j members + double M2 // number of edges * 2 +) +{ + return (2.0 * (double) kjn) / M2 - + ((double) sigma_tot * (double) i_degree) / (M2 * M2) ; +} + +// compute clustering by running the Louvain algorithm against the graph's +// adjacency matrix A +GrB_Info LAGraph_louvain +( + GrB_Vector *com, // output communities + LAGraph_Graph G, // graph adjacency matrix + int itermax, // max number of modularity improvements sweeps per level + int levelmax, // max number of modularity improve and cluster condense + float e, // min change in modularity considered an improvement + char *msg // error message +) +{ + //-------------------------------------------------------------------------- + // check inputs + //-------------------------------------------------------------------------- + + if (com == NULL || G == NULL || msg == NULL) { + return (GrB_NULL_POINTER) ; + } + + GrB_Matrix C = NULL ; // map between node to community id + GrB_Vector D = NULL ; // map between nodes to community ids + GxB_Iterator it = NULL ; + GrB_Vector Ni = NULL ; // node i's neighbors + GrB_Descriptor desc = NULL ; + + uint64_t *degree = NULL ; // node degree + uint64_t *community_degree = NULL ; // sigma tot + + // find out if graph is symmetric, compute G->out_degree, and G->nself_edges + LG_TRY (LAGraph_Cached_IsSymmetricStructure (G, msg)) ; + LG_TRY (LAGraph_Cached_OutDegree (G, msg)) ; + + LG_TRY (LAGraph_Cached_NSelfEdges (G, msg)) ; + LG_ASSERT_MSG (G->nself_edges == 0, GrB_INVALID_VALUE, + "G->nself_edges must be zero") ; + + GrB_Matrix A = G->A ; + + // TODO: make sure `A` is row-wise + + GrB_Index nrows ; + GrB_Index ncols ; + GRB_TRY (GrB_Matrix_nrows (&nrows, A)) ; + GRB_TRY (GrB_Matrix_ncols (&ncols, A)) ; + + // expecting a square matrix + LG_TRY (LAGraph_CheckGraph (G, msg)) ; + LG_ASSERT_MSG (nrows == ncols, LAGRAPH_INVALID_GRAPH, + "adjacency matrix must be square") ; + + GrB_Type t ; + GRB_TRY (GxB_Matrix_type (&t, A)) ; + LG_ASSERT_MSG (t == GrB_BOOL, LAGRAPH_INVALID_GRAPH, "A must be boolean") ; + + GrB_Index nvals ; + GRB_TRY (GrB_Matrix_nvals (&nvals, A)) ; + double M2 = (double) nvals ; + + //-------------------------------------------------------------------------- + // compute nodes degree + //-------------------------------------------------------------------------- + + GRB_TRY (GrB_Vector_new (&D, GrB_UINT64, nrows)) ; + GRB_TRY (GrB_assign (D, NULL, NULL, 0, GrB_ALL, nrows, NULL)) ; + GRB_TRY (GrB_assign (D, NULL, GrB_PLUS_UINT64, G->out_degree, GrB_ALL, + nrows, NULL)) ; + + // TODO: enable once we support weighted graphs + //GRB_TRY (GrB_mxv (D, NULL, GrB_PLUS_UINT64, GxB_PLUS_PAIR_UINT64, A, D, + //NULL)) ; + + // unpack D to an array for direct access + uint64_t degree_size ; + GRB_TRY (GxB_Vector_unpack_Full (D, (void**)°ree, °ree_size, NULL, + NULL)) ; + GRB_TRY (GrB_free (&D)) ; + + // community_degree [i] == degree [i] + // as each node is in its own community + LG_TRY (LAGraph_Malloc ((void **) &community_degree, 1, degree_size, msg)) ; + memcpy (community_degree, degree, degree_size) ; + + // initialize C + uint64_t *I = NULL ; + GRB_TRY (_initCommunities (&C, &I, nrows)) ; + LG_ASSERT (I != NULL, GrB_NULL_POINTER) ; + + GRB_TRY (GrB_Vector_new (&Ni, GrB_UINT64, ncols)) ; + GRB_TRY (GxB_Iterator_new (&it)) ; + + GRB_TRY (GrB_Descriptor_new (&desc)) ; + GRB_TRY (GrB_set (desc, GrB_STRUCTURE, GrB_MASK_FIELD)) ; + GRB_TRY (GrB_set (desc, GxB_USE_INDICES, GxB_ROWINDEX_LIST)) ; + + bool improved ; + double modularity_gain ; + + // start sweep + do + { + modularity_gain = 0 ; + + // for each node + for (GrB_Index i = 0 ; i < nrows ; i++) { + // determine i's current community + uint64_t ic = I [i] ; // i's community + + //------------------------------------------------------------------ + // get a set of communities i can migrate to + //------------------------------------------------------------------ + + // get i's neighbors; A[i:] + GRB_TRY (GrB_Col_extract (Ni, NULL, NULL, A, GrB_ALL, ncols, + i, GrB_DESC_T0)) ; + + // get neighbors communities + // Ni * C = X[i]=j i community ID, j #neighbors in the ith community + GRB_TRY (GrB_vxm (Ni, NULL, NULL, GxB_PLUS_PAIR_UINT64, Ni, C, NULL)) ; + + GRB_TRY (GxB_Vector_Iterator_attach (it, Ni, NULL)) ; + GrB_Info info = GxB_Vector_Iterator_seek (it, 0) ; + + //------------------------------------------------------------------ + // compute base score removing i from its current community + //------------------------------------------------------------------ + + // i's degree + uint64_t i_degree = degree [i] ; + + // number of edges connecting i to its current community + uint64_t kin = 0 ; + GRB_TRY (GrB_Vector_extractElement (&kin, Ni, ic)) ; + + // adjusted sigma tot + uint64_t sigma_tot = community_degree [ic] ; + //LG_ASSERT (sigma_tot >= i_degree) ; + sigma_tot -= i_degree ; + + double base_gain = gain (i_degree, kin, sigma_tot, M2) ; + double max_modularity = base_gain ; + uint64_t best_community = ic ; + + while (info != GxB_EXHAUSTED) + { + // candidate community + GrB_Index jc = GxB_Vector_Iterator_getIndex (it) ; + + // number of edges from node i to candidate community + uint64_t kjn = GxB_Iterator_get_UINT64 (it) ; + + // move to the next entry in Ni + info = GxB_Vector_Iterator_next (it) ; + + // skip i's community + if (jc == ic) + { + continue ; + } + + // total degree of nodes in community + sigma_tot = community_degree [jc] ; + + // modularity gain of i joining the candidate community + double g = gain (i_degree, kjn, sigma_tot, M2) ; + + if (g > max_modularity) + { + best_community = jc ; + max_modularity = g ; + } + } + + if (best_community != ic) { + // migrate i to its new community j + I [i] = best_community ; + + //-------------------------------------------------------------- + // update community degree + // community_degree [ic] -= degree [i] + // community_degree [best_community] += degree [i] + //-------------------------------------------------------------- + + community_degree [ic] -= i_degree ; + community_degree [best_community] += i_degree ; + + // accumulate modularity + // TODO: at the end of the sweep compute modularity from A and C + // compare that againt the initial modularity before the sweep + modularity_gain += max_modularity - base_gain ; + } + } + + improved = (modularity_gain > e) ; + } while (improved && --itermax > 0) ; + + + for (uint l = 0 ; l < levelmax ; l++) + { + //---------------------------------------------------------------------- + // pick a representative for each community + //---------------------------------------------------------------------- + + GrB_Vector R ; // representatives R [i] = n; i community ID, n node ID + GRB_TRY (GrB_Vector_new (&R, GrB_UINT64, nrows)) ; + GRB_TRY (GrB_Matrix_reduce_Monoid (R, NULL, NULL, GrB_MAX_MONOID_UINT64, + C, GrB_DESC_T0)) ; + + //---------------------------------------------------------------------- + // map representative to its members + //---------------------------------------------------------------------- + + // MAP [:i] nodes represented by i + // MAP = MAP * (C * Rdiag) + + //---------------------------------------------------------------------- + // compute A + //---------------------------------------------------------------------- + + // A = CT * C; A[i,j] = x community i is connected to community j with x + // different connections + } + + //-------------------------------------------------------------------------- + // set output + //-------------------------------------------------------------------------- + + GRB_TRY (GrB_Vector_new (com, GrB_UINT64, nrows)) ; + GRB_TRY (GxB_Vector_load (*com, (void **)(&I), GrB_UINT64, nrows, + sizeof (uint64_t) * nrows, GrB_DEFAULT, NULL)) ; + + LG_FREE_WORK ; + + return GrB_SUCCESS ; +} + diff --git a/experimental/test/test_louvain.c b/experimental/test/test_louvain.c index 21ffef646b..b295d888ed 100644 --- a/experimental/test/test_louvain.c +++ b/experimental/test/test_louvain.c @@ -1,152 +1,170 @@ +//---------------------------------------------------------------------------- +// LAGraph/experimental/test/test_louvain.c: test cases for LAGraph_louvain +//---------------------------------------------------------------------------- + +// LAGraph, (c) 2019-2026 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 Roi Lipman and Gabriel Gomez, FalkorDB + +//----------------------------------------------------------------------------- + +// Smoke test for the Louvain community-detection method: run the algorithm on +// a few small graphs and confirm the result is a well-formed community vector. + #include #include #include #include -#include "LG_Xtest.h" - -#define dbg(x) GxB_print(x, 5) -#define err(x, info) \ - if (!(info == GrB_SUCCESS || info == GrB_NO_VALUE)) \ - { \ - char **err; \ - GrB_error(err, x); \ - printf("\ninfo: %d error: %s\n", info, err); \ - } -char msg[LAGRAPH_MSG_LEN]; -LAGraph_Graph G = NULL; -GrB_Matrix A = NULL; + +char msg [LAGRAPH_MSG_LEN] ; +GrB_Matrix A = NULL ; +LAGraph_Graph G = NULL ; +GrB_Vector com = NULL ; + #define LEN 512 -char filename[LEN + 1]; -typedef struct -{ - const char *matrix_file; // Adjeancy matrix or graph - const double mod; -} matrix_info; +char filename [LEN+1] ; -const matrix_info files[] = { - {"comm0.mtx", 0.357142857142857}, - {"karate.mtx", .42}, - {"", -1} -}; +// Louvain parameters used throughout the smoke test. +#define ITERMAX 6 // max modularity-improvement sweeps per level +#define LEVELMAX 2 // max improve-and-condense levels +#define EPSILON 1e-5f // min modularity change considered an improvement +typedef struct +{ + const char *name ; +} +matrix_info ; -void test_LouvainSeq(void) +// pattern / boolean symmetric matrices: LAGraph_louvain requires G->A to be a +// square boolean adjacency matrix, which is exactly what these files load as. +const matrix_info files [ ] = { - #if LG_SUITESPARSE_GRAPHBLAS_V10_2 - LAGraph_Init(msg); - // LG_SET_BURBLE (false) ; + { "A.mtx" }, // tiny 7x7 graph + { "karate.mtx" }, // Zachary's karate club, the classic clustering example + { "" }, +} ; - for (int k = 0;; k++) - { - uint64_t seed = 1224; - - const char *aname = files[k].matrix_file; - if (strlen(aname) == 0) - break; - printf("\n================================== %s:\n", aname); - snprintf(filename, LEN, LG_DATA_DIR "%s", files[k].matrix_file); - FILE *f = fopen(filename, "r"); - TEST_CHECK(f != NULL); - OK(LAGraph_MMRead(&A, f, msg)); - fclose(f); - - OK(LAGraph_New(&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)); - TEST_CHECK(A == NULL); - - - OK(LAGraph_Cached_AT(G, msg)); - // check if the pattern is symmetric - if it isn't make it. - OK(LAGraph_Cached_IsSymmetricStructure(G, msg)); - GrB_Matrix S = NULL; - printf("Sequntial Louvain"); - for (int jit = 0 ; jit <= 1 ; jit++) - { - OK (LG_SET_JIT (jit ? GxB_JIT_ON : GxB_JIT_OFF)) ; - - double tsimple = LAGraph_WallClockTime(); - OK(LAGraph_LouvainSeq(&S, G,seed, msg)); - - // OK(LAGraph_Louvain_res(&S,G,.3,msg)); - tsimple = LAGraph_WallClockTime() - tsimple; - double Q = 0.0; - OK(LAGr_AdjModularity(&Q, 1.0, G->A, S, msg)); - printf("Q:%f\n", Q); - // printf("Number of Communities: %d",comms); - printf(" time: %f\n", tsimple); - } - - OK(LAGraph_Delete(&G, msg)); - } - LAGraph_Finalize(msg); - #endif -} -void test_LouvainIS(void) +//**************************************************************************** +void test_louvain (void) { - #if LG_SUITESPARSE_GRAPHBLAS_V10_2 - LAGraph_Init(msg); - // LG_SET_BURBLE (true) ; + OK (LAGraph_Init (msg)) ; - for (int k = 0;; k++) + for (int k = 0 ; ; k++) { - uint64_t seed = 1249141465; - const char *aname = files[k].matrix_file; - if (strlen(aname) == 0) - break; - printf("\n================================== %s:\n", aname); - snprintf(filename, LEN, LG_DATA_DIR "%s", files[k].matrix_file); - FILE *f = fopen(filename, "r"); - TEST_CHECK(f != NULL); - OK(LAGraph_MMRead(&A, f, msg)); - fclose(f); - - OK(LAGraph_New(&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)); - TEST_CHECK(A == NULL); - - - OK(LAGraph_Cached_AT(G, msg)); - // check if the pattern is symmetric - if it isn't make it. - OK(LAGraph_Cached_IsSymmetricStructure(G, msg)); - GrB_Matrix S = NULL; - printf("Isolate Sets Louvain\n"); - for (int jit = 0 ; jit <= 1 ; jit++) - { - OK (LG_SET_JIT (jit ? GxB_JIT_ON : GxB_JIT_OFF)) ; - - double tsimple = LAGraph_WallClockTime(); - OK(LAGraph_LouvainIS(&S,seed, G, msg)); - // GrB_Info info = (LAGraph_LouvainIS(&S,seed, G, msg)); - // printf ("info %d\nmsg: %s\n", info, msg) ; - // OK (info) ; - - // OK(LAGraph_Louvain_res(&S,G,.3,msg)); - tsimple = LAGraph_WallClockTime() - tsimple; - double Q = 0.0; - double tsimple2 = LAGraph_WallClockTime(); - OK(LAGr_AdjModularity(&Q, 1.0, G->A, S, msg)); - tsimple2 = LAGraph_WallClockTime() - tsimple2; - - printf("Q:%f time to calc Q: %f\n", Q,tsimple2); - // printf("Number of Communities: %d",comms); - printf(" time: %f\n", tsimple); - GrB_free(&S); - } - - OK(LAGraph_Delete(&G, msg)); - + // load the matrix as A + const char *aname = files [k].name ; + if (strlen (aname) == 0) + { + break ; + } + + printf ("\n================================== %s:\n", aname) ; + TEST_CASE (aname) ; + snprintf (filename, LEN, LG_DATA_DIR "%s", aname) ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + + // construct an undirected graph G with adjacency matrix A + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_UNDIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; // A has been moved into G->A + + // Louvain operates on a simple graph + OK (LAGraph_DeleteSelfEdges (G, msg)) ; + + GrB_Index n ; + OK (GrB_Matrix_nrows (&n, G->A)) ; + + // compute the communities with LAGraph_louvain + int result = LAGraph_louvain (&com, G, ITERMAX, LEVELMAX, EPSILON, msg) ; + TEST_CHECK (result == GrB_SUCCESS) ; + + //---------------------------------------------------------------------- + // check the result is a well-formed community vector + //---------------------------------------------------------------------- + + // one community id per node + TEST_CHECK (com != NULL) ; + GrB_Index com_size ; + OK (GrB_Vector_size (&com_size, com)) ; + TEST_CHECK (com_size == n) ; + TEST_MSG ("community vector size %g, expected %g", + (double) com_size, (double) n) ; + + // the result must be a full vector: one community id per node + GrB_Index nvals ; + OK (GrB_Vector_nvals (&nvals, com)) ; + TEST_CHECK (nvals == n) ; + TEST_MSG ("community vector has %g entries, expected %g", + (double) nvals, (double) n) ; + + // print the result + LAGraph_PrintLevel pr = (n <= 100) ? LAGraph_COMPLETE : LAGraph_SHORT ; + printf ("\nlouvain (computed communities):\n") ; + OK (LAGraph_Vector_Print (com, pr, stdout, msg)) ; + + // community ids are node indices, so the largest id must be < n; + // ids are unsigned, so this bounds every id to the valid range [0, n) + uint64_t max_id = 0 ; + OK (GrB_reduce (&max_id, NULL, GrB_MAX_MONOID_UINT64, com, NULL)) ; + TEST_CHECK (max_id < n) ; + TEST_MSG ("largest community id %g, must be < %g", + (double) max_id, (double) n) ; + + OK (GrB_free (&com)) ; + OK (LAGraph_Delete (&G, msg)) ; } - LAGraph_Finalize(msg); - #endif + OK (LAGraph_Finalize (msg)) ; } -TEST_LIST = { - - {"LouvainSeq", test_LouvainSeq}, - {"LouvainIS", test_LouvainIS}, - - {NULL, NULL}}; +//------------------------------------------------------------------------------ +// test_errors +//------------------------------------------------------------------------------ +void test_louvain_errors (void) +{ + OK (LAGraph_Init (msg)) ; + + snprintf (filename, LEN, LG_DATA_DIR "%s", "karate.mtx") ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + TEST_MSG ("Loading of adjacency matrix failed") ; + OK (fclose (f)) ; + + // construct an undirected graph G with adjacency matrix A + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_UNDIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; + OK (LAGraph_DeleteSelfEdges (G, msg)) ; + + // com is NULL + int result = LAGraph_louvain (NULL, G, ITERMAX, LEVELMAX, EPSILON, msg) ; + printf ("\nresult: %d\n", result) ; + TEST_CHECK (result == GrB_NULL_POINTER) ; + + // G is NULL + result = LAGraph_louvain (&com, NULL, ITERMAX, LEVELMAX, EPSILON, msg) ; + printf ("result: %d\n", result) ; + TEST_CHECK (result == GrB_NULL_POINTER) ; + + OK (LAGraph_Delete (&G, msg)) ; + OK (LAGraph_Finalize (msg)) ; +} +//**************************************************************************** - +TEST_LIST = { + {"louvain", test_louvain}, + {"louvain_errors", test_louvain_errors}, + {NULL, NULL} +} ; diff --git a/experimental/test/test_louvain_seqis.c b/experimental/test/test_louvain_seqis.c new file mode 100644 index 0000000000..21ffef646b --- /dev/null +++ b/experimental/test/test_louvain_seqis.c @@ -0,0 +1,152 @@ +#include +#include + +#include +#include +#include "LG_Xtest.h" + +#define dbg(x) GxB_print(x, 5) +#define err(x, info) \ + if (!(info == GrB_SUCCESS || info == GrB_NO_VALUE)) \ + { \ + char **err; \ + GrB_error(err, x); \ + printf("\ninfo: %d error: %s\n", info, err); \ + } +char msg[LAGRAPH_MSG_LEN]; +LAGraph_Graph G = NULL; +GrB_Matrix A = NULL; +#define LEN 512 +char filename[LEN + 1]; +typedef struct +{ + const char *matrix_file; // Adjeancy matrix or graph + const double mod; +} matrix_info; + +const matrix_info files[] = { + {"comm0.mtx", 0.357142857142857}, + {"karate.mtx", .42}, + {"", -1} +}; + + +void test_LouvainSeq(void) +{ + #if LG_SUITESPARSE_GRAPHBLAS_V10_2 + LAGraph_Init(msg); + // LG_SET_BURBLE (false) ; + + for (int k = 0;; k++) + { + uint64_t seed = 1224; + + const char *aname = files[k].matrix_file; + if (strlen(aname) == 0) + break; + printf("\n================================== %s:\n", aname); + snprintf(filename, LEN, LG_DATA_DIR "%s", files[k].matrix_file); + FILE *f = fopen(filename, "r"); + TEST_CHECK(f != NULL); + OK(LAGraph_MMRead(&A, f, msg)); + fclose(f); + + OK(LAGraph_New(&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)); + TEST_CHECK(A == NULL); + + + OK(LAGraph_Cached_AT(G, msg)); + // check if the pattern is symmetric - if it isn't make it. + OK(LAGraph_Cached_IsSymmetricStructure(G, msg)); + GrB_Matrix S = NULL; + printf("Sequntial Louvain"); + for (int jit = 0 ; jit <= 1 ; jit++) + { + OK (LG_SET_JIT (jit ? GxB_JIT_ON : GxB_JIT_OFF)) ; + + double tsimple = LAGraph_WallClockTime(); + OK(LAGraph_LouvainSeq(&S, G,seed, msg)); + + // OK(LAGraph_Louvain_res(&S,G,.3,msg)); + tsimple = LAGraph_WallClockTime() - tsimple; + double Q = 0.0; + OK(LAGr_AdjModularity(&Q, 1.0, G->A, S, msg)); + printf("Q:%f\n", Q); + // printf("Number of Communities: %d",comms); + printf(" time: %f\n", tsimple); + } + + OK(LAGraph_Delete(&G, msg)); + } + LAGraph_Finalize(msg); + #endif +} +void test_LouvainIS(void) +{ + #if LG_SUITESPARSE_GRAPHBLAS_V10_2 + LAGraph_Init(msg); + // LG_SET_BURBLE (true) ; + + for (int k = 0;; k++) + { + uint64_t seed = 1249141465; + const char *aname = files[k].matrix_file; + if (strlen(aname) == 0) + break; + printf("\n================================== %s:\n", aname); + snprintf(filename, LEN, LG_DATA_DIR "%s", files[k].matrix_file); + FILE *f = fopen(filename, "r"); + TEST_CHECK(f != NULL); + OK(LAGraph_MMRead(&A, f, msg)); + fclose(f); + + OK(LAGraph_New(&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)); + TEST_CHECK(A == NULL); + + + OK(LAGraph_Cached_AT(G, msg)); + // check if the pattern is symmetric - if it isn't make it. + OK(LAGraph_Cached_IsSymmetricStructure(G, msg)); + GrB_Matrix S = NULL; + printf("Isolate Sets Louvain\n"); + for (int jit = 0 ; jit <= 1 ; jit++) + { + OK (LG_SET_JIT (jit ? GxB_JIT_ON : GxB_JIT_OFF)) ; + + double tsimple = LAGraph_WallClockTime(); + OK(LAGraph_LouvainIS(&S,seed, G, msg)); + // GrB_Info info = (LAGraph_LouvainIS(&S,seed, G, msg)); + // printf ("info %d\nmsg: %s\n", info, msg) ; + // OK (info) ; + + // OK(LAGraph_Louvain_res(&S,G,.3,msg)); + tsimple = LAGraph_WallClockTime() - tsimple; + double Q = 0.0; + double tsimple2 = LAGraph_WallClockTime(); + OK(LAGr_AdjModularity(&Q, 1.0, G->A, S, msg)); + tsimple2 = LAGraph_WallClockTime() - tsimple2; + + printf("Q:%f time to calc Q: %f\n", Q,tsimple2); + // printf("Number of Communities: %d",comms); + printf(" time: %f\n", tsimple); + GrB_free(&S); + } + + OK(LAGraph_Delete(&G, msg)); + + } + + LAGraph_Finalize(msg); + #endif +} + +TEST_LIST = { + + {"LouvainSeq", test_LouvainSeq}, + {"LouvainIS", test_LouvainIS}, + + {NULL, NULL}}; + + + + From 00b6a9e49b828666278fa5b83ca4162b2740876d Mon Sep 17 00:00:00 2001 From: GomezGab Date: Fri, 4 Sep 2026 19:34:59 -0500 Subject: [PATCH 12/13] [WIP] Louvain --- experimental/algorithm/LAGraph_Leiden.c | 8 +- experimental/algorithm/LAGraph_louvain.c | 783 ++++++++++++++--------- include/LAGraphX.h | 11 + 3 files changed, 504 insertions(+), 298 deletions(-) diff --git a/experimental/algorithm/LAGraph_Leiden.c b/experimental/algorithm/LAGraph_Leiden.c index f37b583169..6b10fd2f4c 100644 --- a/experimental/algorithm/LAGraph_Leiden.c +++ b/experimental/algorithm/LAGraph_Leiden.c @@ -135,13 +135,14 @@ int LG_Leiden_move_nodes uint64_t *nodes_moved_handle, // nodes moved to a new community // input: const GrB_Matrix A, // adjacency matrix - const GrB_Vector deg, // degree vector + const GrB_Vector deg, // weighteed degree vector uint64_t *queue, // queue to use (contains all nodes) bool *enqueued, // nodes in queue (all at the start) - double m_inv2, // -1 / (2 * m) + double m_inv2, // -1 / (2 * m) char* msg ) { uint64_t node_id, com_id, n ; + // TODO: comment GrB_Vector x = NULL ; GrB_Vector c_deg_vec = NULL ; GrB_Vector c_size_vec = NULL ; @@ -192,6 +193,7 @@ int LG_Leiden_move_nodes while (queue_head != queue_tail) { // while queue not empty nodes_popped++; LG_QUEUE_DEQUEUE (queue, queue_head, queue_size, node_id) ; + // TODO: put in macro enqueued [node_id] = false ; uint64_t n_neighbors; com_id = community [node_id] ; @@ -200,7 +202,7 @@ int LG_Leiden_move_nodes GRB_TRY (GrB_Vector_clear (x)); // GRB_TRY (GrB_Vector_setElement_FP64 (x, 0.0, node_id)) ; // GRB_TRY (GrB_vxm (x, NULL, GrB_PLUS_FP64, GxB_PLUS_SECOND_FP64, x, A, NULL)) ; - // TODO: deg and c_deg as C_arrays? + // TODO: is this set needed? GRB_TRY (GrB_Vector_setElement_FP64 (x, 0.0, node_id)) ; GRB_TRY (GrB_Col_extract (x, NULL, GrB_PLUS_FP64, A, GrB_ALL, n, node_id, GrB_DESC_T0)) ; diff --git a/experimental/algorithm/LAGraph_louvain.c b/experimental/algorithm/LAGraph_louvain.c index cdc51bf6be..e5ac908686 100644 --- a/experimental/algorithm/LAGraph_louvain.c +++ b/experimental/algorithm/LAGraph_louvain.c @@ -22,339 +22,532 @@ #undef LG_FREE_ALL #define LG_FREE_ALL \ { \ - GrB_free(&V); \ - GrB_free(&cont); \ + GrB_free (&V) ; \ + GrB_free (&cont) ; \ } -// construct C, C[i:] = [i] -// and N, N[i] = i's community id +// construct C, C[i,:] = [i] +// and I, where I[i] is i's community id (readonly handle to C's col-indices) static GrB_Info _initCommunities ( - GrB_Matrix *C, // C [i] = i - uint64_t **I, // C's indices array - GrB_Index node_count // node count + GrB_Matrix *C, // C [i] = i + uint64_t **I, // C's indices array + GrB_Index node_count // node count ) { - //-------------------------------------------------------------------------- - // initialize C - //-------------------------------------------------------------------------- + //-------------------------------------------------------------------------- + // initialize C + //-------------------------------------------------------------------------- - // create full vector V[i] = 1 - char *msg = NULL ; - GrB_Vector V = NULL ; - GxB_Container cont = NULL ; + char *msg = NULL ; + GrB_Vector V = NULL ; + GxB_Container cont = NULL ; - GRB_TRY (GrB_Vector_new (&V, GrB_BOOL, node_count)) ; - GRB_TRY (GrB_assign (V, NULL, NULL, true, GrB_ALL, node_count, NULL)) ; + GRB_TRY (GrB_Vector_new (&V, GrB_BOOL, node_count)) ; + GRB_TRY (GrB_assign (V, NULL, NULL, true, GrB_ALL, node_count, NULL)) ; + GRB_TRY (GrB_Matrix_diag (C, V, 0)) ; + GRB_TRY (GrB_free (&V)) ; - // create diagonal matrix C from V - GRB_TRY (GrB_Matrix_diag (C, V, 0)) ; - GRB_TRY (GrB_free (&V)) ; + // C should be CSR with 64-bit indices + GRB_TRY (GrB_set (*C, 64, GxB_COLINDEX_INTEGER_HINT )) ; + GRB_TRY (GrB_set (*C, GxB_SPARSE, GxB_SPARSITY_CONTROL )) ; + GRB_TRY (GrB_set (*C, GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT)) ; - // C should be CSR with 64 bit indicies - GRB_TRY (GrB_set (*C, 64, GxB_COLINDEX_INTEGER_HINT )) ; - GRB_TRY (GrB_set (*C, GxB_SPARSE, GxB_SPARSITY_CONTROL )) ; - GRB_TRY (GrB_set (*C, GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT)) ; + //-------------------------------------------------------------------------- + // get a handle to C's indices array + //-------------------------------------------------------------------------- - //-------------------------------------------------------------------------- - // get a handle to C's indicies array - //-------------------------------------------------------------------------- + GRB_TRY (GxB_Container_new (&cont)) ; + GRB_TRY (GxB_unload_Matrix_into_Container (*C, cont, NULL)) ; - // unload C into a container - GRB_TRY (GxB_Container_new (&cont)) ; + int handling ; + GrB_Type type ; + uint64_t n, X_memsize ; - // GrB_Matrix -> GxB_Container - GRB_TRY (GxB_unload_Matrix_into_Container (*C, cont, NULL)) ; + GRB_TRY (GxB_Vector_unload (cont->i, (void**) I, &type, &n, &X_memsize, + &handling, NULL)) ; - int handling ; - GrB_Type type ; - uint64_t n, X_memsize ; + // load I and mark it read-only + GRB_TRY (GxB_Vector_load (cont->i, (void**) I, type, n, X_memsize, + GxB_IS_READONLY, NULL)) ; - // unload I - GRB_TRY (GxB_Vector_unload (cont->i, (void**)I, &type, &n, &X_memsize, - &handling, NULL)) ; + GRB_TRY (GxB_load_Matrix_from_Container (*C, cont, NULL)) ; + GRB_TRY (GrB_free (&cont)) ; - // load I, mark it as READONLY - GRB_TRY (GxB_Vector_load (cont->i, (void **)I, type, n, X_memsize, - GxB_IS_READONLY, NULL)) ; - - // GxB_Container -> GrB_Matrix - GRB_TRY (GxB_load_Matrix_from_Container (*C, cont, NULL)) ; - GRB_TRY (GrB_free (&cont)) ; - - LG_FREE_ALL ; + LG_FREE_ALL ; + return GrB_SUCCESS ; +} - return GrB_SUCCESS ; +// computes the modularity contribution of attaching a node of degree i_degree, +// with kin edges into a community, to that community. +static inline double gain +( + uint64_t i_degree, + uint64_t kjn, + uint64_t sigma_tot, + double M2 +) +{ + return (2.0 * (double) kjn) / M2 - + ((double) sigma_tot * (double) i_degree) / (M2 * M2) ; } #undef LG_FREE_WORK -#define LG_FREE_WORK \ -{ \ - GrB_free(&C); \ - GrB_free(&D); \ - GrB_free(&it); \ - GrB_free(&Ni); \ - GrB_free(&desc); \ - LAGraph_Free((void**)°ree, NULL); \ - LAGraph_Free((void**)&community_degree, NULL); \ +#define LG_FREE_WORK \ +{ \ + GrB_free (&desc) ; \ + GrB_free (&c_list) ; \ + GrB_free (&x) ; \ + GrB_free (&S_squished) ; \ + GrB_free (&S_new) ; \ + GrB_free (&A_squished) ; \ + GrB_free (&A_new) ; \ } #undef LG_FREE_ALL -#define LG_FREE_ALL \ -{ \ - LG_FREE_WORK ; \ +#define LG_FREE_ALL \ +{ \ + if (free_A) GrB_free (&A) ; \ + LG_FREE_WORK ; \ } -// computes the modularity contribution of attaching a node of degree -// i_degree, with kin edges into a community, to that community -- where -// sigma_tot is the total degree of the community's members, *excluding* -// the node itself if it happens to already be a member (the caller is -// responsible for that adjustment; see base_gain) -static inline double gain +// helper function: phase 2 of Louvain. +// Input C is node->community (possibly with empty community columns). +// Output S maps coarse nodes to previous level nodes, and A_new is coarsened. +static GrB_Info LG_Louvain_aggregate ( - uint64_t i_degree, // node i's degree - uint64_t kjn, // number of edges connecting node i to community j - uint64_t sigma_tot, // total degree of community j members - double M2 // number of edges * 2 + GrB_Matrix *S_handle, // output: coarse->prev level map (n_new x n_old) + GrB_Matrix *A_handle, // output: coarsened adjacency (n_new x n_new) + const GrB_Matrix C, // input node->community map (n_old x n_old) + GrB_Matrix A, // input adjacency (n_old x n_old) + bool free_A, // if true, helper may free input A + char *msg ) { - return (2.0 * (double) kjn) / M2 - - ((double) sigma_tot * (double) i_degree) / (M2 * M2) ; + GrB_Vector c_list = NULL, x = NULL ; + GrB_Matrix S_squished = NULL, S_new = NULL ; + GrB_Matrix A_squished = NULL, A_new = NULL ; + GrB_Descriptor desc = NULL ; + + uint64_t n, n_new ; + GRB_TRY (GrB_Matrix_nrows (&n, A)) ; + + GrB_Type type_A ; + GrB_Semiring plus_first, plus_second ; + int code = 0 ; + GRB_TRY (GrB_get (A, &code, GrB_EL_TYPE_CODE)) ; + switch (code) + { + case GrB_BOOL_CODE: + case GrB_INT8_CODE: + case GrB_INT16_CODE: + case GrB_INT32_CODE: + case GrB_INT64_CODE: + case GrB_UINT8_CODE: + case GrB_UINT16_CODE: + case GrB_UINT32_CODE: + case GrB_UINT64_CODE: + type_A = GrB_INT64 ; + plus_first = GxB_PLUS_FIRST_INT64 ; + plus_second = GxB_PLUS_SECOND_INT64 ; + break ; + + case GrB_FP32_CODE: + case GrB_FP64_CODE: + type_A = GrB_FP64 ; + plus_first = GxB_PLUS_FIRST_FP64 ; + plus_second = GxB_PLUS_SECOND_FP64 ; + break ; + + default: + LG_ERROR_MSG ("LAGraph failed (file %s, line %d): unsupported " + "adjacency matrix type for Louvain aggregation", + __FILE__, __LINE__) ; + LG_FREE_ALL ; + return (GrB_NOT_IMPLEMENTED) ; + } + + GRB_TRY (GrB_Descriptor_new (&desc)) ; + GRB_TRY (GrB_set (desc, GxB_USE_INDICES, GxB_ROWINDEX_LIST)) ; + GRB_TRY (GrB_set (desc, GxB_USE_INDICES, GxB_COLINDEX_LIST)) ; + + // Find active communities as column indices in C. + GRB_TRY (GrB_Vector_new (&c_list, GrB_BOOL, n)) ; + GRB_TRY (GrB_Vector_new (&x, GrB_BOOL, n)) ; + GRB_TRY (GrB_assign (x, NULL, NULL, (bool) 0, GrB_ALL, n, NULL)) ; + GRB_TRY (GrB_vxm (c_list, NULL, NULL, GxB_ANY_PAIR_BOOL, x, C, NULL)) ; + GRB_TRY (GrB_Vector_nvals (&n_new, c_list)) ; + LG_ASSERT_MSG (n_new > 0, GrB_INVALID_VALUE, "No active communities") ; + + // S_squished: previous node -> active-community-id (contiguous columns) + GRB_TRY (GrB_Matrix_new (&S_squished, GrB_BOOL, n, n_new)) ; + GRB_TRY (GxB_Matrix_extract_Vector ( + S_squished, NULL, NULL, C, NULL, c_list, desc)) ; + + // A_squished: aggregate destination by community + GRB_TRY (GrB_Matrix_new (&A_squished, type_A, n, n_new)) ; + GRB_TRY (GrB_mxm (A_squished, NULL, NULL, plus_first, A, S_squished, NULL)) ; + + // S_new: active-community-id -> previous nodes + GRB_TRY (GrB_Matrix_new (&S_new, GrB_BOOL, n_new, n)) ; + GRB_TRY (GrB_transpose (S_new, NULL, NULL, S_squished, NULL)) ; + + // A_new: aggregate source by community + GRB_TRY (GrB_Matrix_new (&A_new, type_A, n_new, n_new)) ; + GRB_TRY (GrB_set (A_new, GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT)) ; + GRB_TRY (GrB_mxm (A_new, NULL, NULL, plus_second, S_new, A_squished, NULL)) ; + + if (free_A) + { + GRB_TRY (GrB_free (&A)) ; + } + + *S_handle = S_new ; S_new = NULL ; + *A_handle = A_new ; A_new = NULL ; + + LG_FREE_WORK ; + return GrB_SUCCESS ; +} + +#undef LG_FREE_WORK +#define LG_FREE_WORK \ +{ \ + GrB_free (&C) ; \ + GrB_free (&D) ; \ + GrB_free (&it) ; \ + GrB_free (&Ni) ; \ + GrB_free (&desc) ; \ + GrB_free (&node_map) ; \ + GrB_free (&new_node_map) ; \ + GrB_free (&S_level) ; \ + GrB_free (&A_next) ; \ + GrB_free (&C_to_orig) ; \ + GrB_free (&active_comms) ; \ + GrB_free (&x_active) ; \ + LAGraph_Free ((void**) °ree, NULL) ; \ + LAGraph_Free ((void**) &community_degree, NULL); \ + LAGraph_Free ((void**) &tuple_i, NULL) ; \ + LAGraph_Free ((void**) &tuple_j, NULL) ; \ + LAGraph_Free ((void**) &tuple_x, NULL) ; \ + LAGraph_Free ((void**) &com_i, NULL) ; \ + LAGraph_Free ((void**) &com_x, NULL) ; \ + LAGraph_Free ((void**) &d_i, NULL) ; \ + LAGraph_Free ((void**) &d_x, NULL) ; \ +} + +#undef LG_FREE_ALL +#define LG_FREE_ALL \ +{ \ + if (free_A) GrB_free (&A) ; \ + LG_FREE_WORK ; \ } // compute clustering by running the Louvain algorithm against the graph's // adjacency matrix A GrB_Info LAGraph_louvain ( - GrB_Vector *com, // output communities - LAGraph_Graph G, // graph adjacency matrix - int itermax, // max number of modularity improvements sweeps per level - int levelmax, // max number of modularity improve and cluster condense - float e, // min change in modularity considered an improvement - char *msg // error message + GrB_Vector *com, // output communities + LAGraph_Graph G, // graph adjacency matrix + int itermax, // max number of modularity improvements sweeps per level + int levelmax, // max number of modularity improve and cluster condense + float e, // min change in modularity considered an improvement + char *msg // error message ) { - //-------------------------------------------------------------------------- - // check inputs - //-------------------------------------------------------------------------- - - if (com == NULL || G == NULL || msg == NULL) { - return (GrB_NULL_POINTER) ; - } - - GrB_Matrix C = NULL ; // map between node to community id - GrB_Vector D = NULL ; // map between nodes to community ids - GxB_Iterator it = NULL ; - GrB_Vector Ni = NULL ; // node i's neighbors - GrB_Descriptor desc = NULL ; - - uint64_t *degree = NULL ; // node degree - uint64_t *community_degree = NULL ; // sigma tot - - // find out if graph is symmetric, compute G->out_degree, and G->nself_edges + if (com == NULL || G == NULL || msg == NULL) + { + return (GrB_NULL_POINTER) ; + } + + GrB_Matrix C = NULL ; // current level node->community map + GrB_Vector D = NULL ; // current level degree vector + GxB_Iterator it = NULL ; + GrB_Vector Ni = NULL ; // node i's neighbors aggregated by comm + GrB_Descriptor desc = NULL ; + + GrB_Matrix node_map = NULL ; // current-level node -> original node map + GrB_Matrix new_node_map = NULL ; + GrB_Matrix S_level = NULL ; // coarse->current node map from aggregate + GrB_Matrix A_next = NULL ; // coarsened adjacency + GrB_Matrix C_to_orig = NULL ; + GrB_Vector active_comms = NULL ; + GrB_Vector x_active = NULL ; + + uint64_t *degree = NULL ; // node degree for current level + uint64_t *community_degree = NULL ; // sigma_tot per community id + bool *tuple_x = NULL ; // values for C_to_orig extraction + GrB_Index *tuple_i = NULL, *tuple_j = NULL ; + GrB_Index *com_i = NULL ; + uint64_t *com_x = NULL ; + GrB_Index *d_i = NULL ; + uint64_t *d_x = NULL ; + + uint64_t *I = NULL ; // readonly handle to C col-indices + GrB_Index original_n = 0 ; + GrB_Index nrows = 0, ncols = 0 ; + GrB_Index final_nrows = 0 ; + + GrB_Matrix A = G->A ; // current adjacency over levels + bool free_A = false ; // true once A becomes internally allocated + + // find out if graph is symmetric, compute cached values, and check loops LG_TRY (LAGraph_Cached_IsSymmetricStructure (G, msg)) ; LG_TRY (LAGraph_Cached_OutDegree (G, msg)) ; - LG_TRY (LAGraph_Cached_NSelfEdges (G, msg)) ; LG_ASSERT_MSG (G->nself_edges == 0, GrB_INVALID_VALUE, - "G->nself_edges must be zero") ; - - GrB_Matrix A = G->A ; - - // TODO: make sure `A` is row-wise - - GrB_Index nrows ; - GrB_Index ncols ; - GRB_TRY (GrB_Matrix_nrows (&nrows, A)) ; - GRB_TRY (GrB_Matrix_ncols (&ncols, A)) ; - - // expecting a square matrix - LG_TRY (LAGraph_CheckGraph (G, msg)) ; - LG_ASSERT_MSG (nrows == ncols, LAGRAPH_INVALID_GRAPH, - "adjacency matrix must be square") ; - - GrB_Type t ; - GRB_TRY (GxB_Matrix_type (&t, A)) ; - LG_ASSERT_MSG (t == GrB_BOOL, LAGRAPH_INVALID_GRAPH, "A must be boolean") ; - - GrB_Index nvals ; - GRB_TRY (GrB_Matrix_nvals (&nvals, A)) ; - double M2 = (double) nvals ; - - //-------------------------------------------------------------------------- - // compute nodes degree - //-------------------------------------------------------------------------- - - GRB_TRY (GrB_Vector_new (&D, GrB_UINT64, nrows)) ; - GRB_TRY (GrB_assign (D, NULL, NULL, 0, GrB_ALL, nrows, NULL)) ; - GRB_TRY (GrB_assign (D, NULL, GrB_PLUS_UINT64, G->out_degree, GrB_ALL, - nrows, NULL)) ; - - // TODO: enable once we support weighted graphs - //GRB_TRY (GrB_mxv (D, NULL, GrB_PLUS_UINT64, GxB_PLUS_PAIR_UINT64, A, D, - //NULL)) ; - - // unpack D to an array for direct access - uint64_t degree_size ; - GRB_TRY (GxB_Vector_unpack_Full (D, (void**)°ree, °ree_size, NULL, - NULL)) ; - GRB_TRY (GrB_free (&D)) ; - - // community_degree [i] == degree [i] - // as each node is in its own community - LG_TRY (LAGraph_Malloc ((void **) &community_degree, 1, degree_size, msg)) ; - memcpy (community_degree, degree, degree_size) ; - - // initialize C - uint64_t *I = NULL ; - GRB_TRY (_initCommunities (&C, &I, nrows)) ; - LG_ASSERT (I != NULL, GrB_NULL_POINTER) ; - - GRB_TRY (GrB_Vector_new (&Ni, GrB_UINT64, ncols)) ; - GRB_TRY (GxB_Iterator_new (&it)) ; - - GRB_TRY (GrB_Descriptor_new (&desc)) ; - GRB_TRY (GrB_set (desc, GrB_STRUCTURE, GrB_MASK_FIELD)) ; - GRB_TRY (GrB_set (desc, GxB_USE_INDICES, GxB_ROWINDEX_LIST)) ; - - bool improved ; - double modularity_gain ; - - // start sweep - do - { - modularity_gain = 0 ; - - // for each node - for (GrB_Index i = 0 ; i < nrows ; i++) { - // determine i's current community - uint64_t ic = I [i] ; // i's community - - //------------------------------------------------------------------ - // get a set of communities i can migrate to - //------------------------------------------------------------------ - - // get i's neighbors; A[i:] - GRB_TRY (GrB_Col_extract (Ni, NULL, NULL, A, GrB_ALL, ncols, - i, GrB_DESC_T0)) ; - - // get neighbors communities - // Ni * C = X[i]=j i community ID, j #neighbors in the ith community - GRB_TRY (GrB_vxm (Ni, NULL, NULL, GxB_PLUS_PAIR_UINT64, Ni, C, NULL)) ; - - GRB_TRY (GxB_Vector_Iterator_attach (it, Ni, NULL)) ; - GrB_Info info = GxB_Vector_Iterator_seek (it, 0) ; - - //------------------------------------------------------------------ - // compute base score removing i from its current community - //------------------------------------------------------------------ - - // i's degree - uint64_t i_degree = degree [i] ; - - // number of edges connecting i to its current community - uint64_t kin = 0 ; - GRB_TRY (GrB_Vector_extractElement (&kin, Ni, ic)) ; - - // adjusted sigma tot - uint64_t sigma_tot = community_degree [ic] ; - //LG_ASSERT (sigma_tot >= i_degree) ; - sigma_tot -= i_degree ; - - double base_gain = gain (i_degree, kin, sigma_tot, M2) ; - double max_modularity = base_gain ; - uint64_t best_community = ic ; - - while (info != GxB_EXHAUSTED) - { - // candidate community - GrB_Index jc = GxB_Vector_Iterator_getIndex (it) ; - - // number of edges from node i to candidate community - uint64_t kjn = GxB_Iterator_get_UINT64 (it) ; - - // move to the next entry in Ni - info = GxB_Vector_Iterator_next (it) ; - - // skip i's community - if (jc == ic) - { - continue ; - } - - // total degree of nodes in community - sigma_tot = community_degree [jc] ; - - // modularity gain of i joining the candidate community - double g = gain (i_degree, kjn, sigma_tot, M2) ; - - if (g > max_modularity) - { - best_community = jc ; - max_modularity = g ; - } - } - - if (best_community != ic) { - // migrate i to its new community j - I [i] = best_community ; - - //-------------------------------------------------------------- - // update community degree - // community_degree [ic] -= degree [i] - // community_degree [best_community] += degree [i] - //-------------------------------------------------------------- - - community_degree [ic] -= i_degree ; - community_degree [best_community] += i_degree ; - - // accumulate modularity - // TODO: at the end of the sweep compute modularity from A and C - // compare that againt the initial modularity before the sweep - modularity_gain += max_modularity - base_gain ; - } - } - - improved = (modularity_gain > e) ; - } while (improved && --itermax > 0) ; - - - for (uint l = 0 ; l < levelmax ; l++) - { - //---------------------------------------------------------------------- - // pick a representative for each community - //---------------------------------------------------------------------- - - GrB_Vector R ; // representatives R [i] = n; i community ID, n node ID - GRB_TRY (GrB_Vector_new (&R, GrB_UINT64, nrows)) ; - GRB_TRY (GrB_Matrix_reduce_Monoid (R, NULL, NULL, GrB_MAX_MONOID_UINT64, - C, GrB_DESC_T0)) ; - - //---------------------------------------------------------------------- - // map representative to its members - //---------------------------------------------------------------------- - - // MAP [:i] nodes represented by i - // MAP = MAP * (C * Rdiag) - - //---------------------------------------------------------------------- - // compute A - //---------------------------------------------------------------------- - - // A = CT * C; A[i,j] = x community i is connected to community j with x - // different connections - } - - //-------------------------------------------------------------------------- - // set output - //-------------------------------------------------------------------------- - - GRB_TRY (GrB_Vector_new (com, GrB_UINT64, nrows)) ; - GRB_TRY (GxB_Vector_load (*com, (void **)(&I), GrB_UINT64, nrows, - sizeof (uint64_t) * nrows, GrB_DEFAULT, NULL)) ; - - LG_FREE_WORK ; - - return GrB_SUCCESS ; + "G->nself_edges must be zero") ; + + GRB_TRY (GrB_Matrix_nrows (&nrows, A)) ; + GRB_TRY (GrB_Matrix_ncols (&ncols, A)) ; + original_n = nrows ; + + LG_TRY (LAGraph_CheckGraph (G, msg)) ; + LG_ASSERT_MSG (nrows == ncols, LAGRAPH_INVALID_GRAPH, + "adjacency matrix must be square") ; + + GrB_Type t ; + GRB_TRY (GxB_Matrix_type (&t, A)) ; + LG_ASSERT_MSG (t == GrB_BOOL, LAGRAPH_INVALID_GRAPH, "A must be boolean") ; + + bool improved = true ; + + for (int level = 0 ; level < levelmax ; level++) + { + GRB_TRY (GrB_Matrix_nrows (&nrows, A)) ; + GRB_TRY (GrB_Matrix_ncols (&ncols, A)) ; + LG_ASSERT_MSG (nrows == ncols, LAGRAPH_INVALID_GRAPH, + "coarsened adjacency matrix must be square") ; + final_nrows = nrows ; + + //---------------------------------------------------------------------- + // initialize degree for this level + //---------------------------------------------------------------------- + + GRB_TRY (GrB_Vector_new (&D, GrB_UINT64, nrows)) ; + GRB_TRY (GrB_assign (D, NULL, NULL, (uint64_t) 0, GrB_ALL, nrows, NULL)) ; + GRB_TRY (GrB_reduce (D, NULL, GrB_PLUS_UINT64, GrB_PLUS_MONOID_UINT64, + A, NULL)) ; + + LG_TRY (LAGraph_Malloc ((void**) °ree, nrows, sizeof (uint64_t), msg)) ; + memset (degree, 0, nrows * sizeof (uint64_t)) ; + + GrB_Index d_nvals = 0 ; + GRB_TRY (GrB_Vector_nvals (&d_nvals, D)) ; + LG_TRY (LAGraph_Malloc ((void**) &d_i, d_nvals, sizeof (GrB_Index), msg)) ; + LG_TRY (LAGraph_Malloc ((void**) &d_x, d_nvals, sizeof (uint64_t), msg)) ; + GRB_TRY (GrB_Vector_extractTuples_UINT64 (d_i, d_x, &d_nvals, D)) ; + for (GrB_Index k = 0 ; k < d_nvals ; k++) + { + degree [d_i [k]] = d_x [k] ; + } + LG_TRY (LAGraph_Free ((void**) &d_i, msg)) ; + LG_TRY (LAGraph_Free ((void**) &d_x, msg)) ; + GRB_TRY (GrB_free (&D)) ; + + LG_TRY (LAGraph_Malloc ((void**) &community_degree, nrows, + sizeof (uint64_t), msg)) ; + memcpy (community_degree, degree, nrows * sizeof (uint64_t)) ; + + //---------------------------------------------------------------------- + // initialize communities C and readonly community id handle I + //---------------------------------------------------------------------- + + GRB_TRY (_initCommunities (&C, &I, nrows)) ; + LG_ASSERT (I != NULL, GrB_NULL_POINTER) ; + + GRB_TRY (GrB_Vector_new (&Ni, GrB_UINT64, ncols)) ; + GRB_TRY (GxB_Iterator_new (&it)) ; + GRB_TRY (GrB_Descriptor_new (&desc)) ; + GRB_TRY (GrB_set (desc, GrB_STRUCTURE, GrB_MASK_FIELD)) ; + GRB_TRY (GrB_set (desc, GxB_USE_INDICES, GxB_ROWINDEX_LIST)) ; + + //---------------------------------------------------------------------- + // phase 1: local moving + //---------------------------------------------------------------------- + + double modularity_gain = 0 ; + improved = true ; + + double M2 = 0.0 ; + for (GrB_Index i = 0 ; i < nrows ; i++) + { + M2 += (double) degree [i] ; + } + LG_ASSERT_MSG (M2 > 0.0, GrB_INVALID_VALUE, + "sum of degrees must be positive") ; + + for (int iter = 0 ; iter < itermax && improved ; iter++) + { + modularity_gain = 0 ; + + for (GrB_Index i = 0 ; i < nrows ; i++) + { + uint64_t ic = I [i] ; + + GRB_TRY (GrB_Col_extract ( + Ni, NULL, NULL, A, GrB_ALL, ncols, i, GrB_DESC_T0)) ; + + // Aggregate node i's edge-weight to each candidate community. + GRB_TRY (GrB_vxm ( + Ni, NULL, NULL, GxB_PLUS_FIRST_UINT64, Ni, C, NULL)) ; + + uint64_t i_degree = degree [i] ; + + uint64_t kin = 0 ; + GrB_Info kin_info = GrB_Vector_extractElement (&kin, Ni, ic) ; + if (kin_info != GrB_SUCCESS && kin_info != GrB_NO_VALUE) + { + GRB_TRY (kin_info) ; + } + if (kin_info == GrB_NO_VALUE) kin = 0 ; + + uint64_t sigma_tot = community_degree [ic] ; + sigma_tot -= i_degree ; + + double base_gain = gain (i_degree, kin, sigma_tot, M2) ; + double max_modularity = base_gain ; + uint64_t best_community = ic ; + + GRB_TRY (GxB_Vector_Iterator_attach (it, Ni, NULL)) ; + GrB_Info info = GxB_Vector_Iterator_seek (it, 0) ; + while (info != GxB_EXHAUSTED) + { + GrB_Index jc = GxB_Vector_Iterator_getIndex (it) ; + uint64_t kjn = GxB_Iterator_get_UINT64 (it) ; + info = GxB_Vector_Iterator_next (it) ; + + if (jc == ic) + { + continue ; + } + + sigma_tot = community_degree [jc] ; + double g = gain (i_degree, kjn, sigma_tot, M2) ; + if (g > max_modularity) + { + best_community = jc ; + max_modularity = g ; + } + } + + if (best_community != ic) + { + I [i] = best_community ; + community_degree [ic] -= i_degree ; + community_degree [best_community] += i_degree ; + modularity_gain += max_modularity - base_gain ; + } + } + + improved = (modularity_gain > (double) e) ; + } + + //---------------------------------------------------------------------- + // decide whether to stop at this level before aggregation + //---------------------------------------------------------------------- + + if (!improved || (level + 1) >= levelmax || nrows <= 1) + { + break ; + } + + // Count active communities at this level. + GRB_TRY (GrB_Vector_new (&active_comms, GrB_BOOL, nrows)) ; + GRB_TRY (GrB_Vector_new (&x_active, GrB_BOOL, nrows)) ; + GRB_TRY (GrB_assign (x_active, NULL, NULL, (bool) 0, GrB_ALL, nrows, NULL)) ; + GRB_TRY (GrB_vxm (active_comms, NULL, NULL, GxB_ANY_PAIR_BOOL, + x_active, C, NULL)) ; + GrB_Index n_communities = 0 ; + GRB_TRY (GrB_Vector_nvals (&n_communities, active_comms)) ; + GRB_TRY (GrB_free (&active_comms)) ; + GRB_TRY (GrB_free (&x_active)) ; + + // No further coarsening is possible. + if (n_communities == nrows || n_communities <= 1) + { + break ; + } + + //---------------------------------------------------------------------- + // phase 2: aggregate + //---------------------------------------------------------------------- + + GRB_TRY (LG_Louvain_aggregate (&S_level, &A_next, C, A, free_A, msg)) ; + + // Compose level maps so node_map always maps current nodes -> original. + if (node_map == NULL) + { + node_map = S_level ; + S_level = NULL ; + } + else + { + GRB_TRY (GrB_Matrix_new (&new_node_map, GrB_BOOL, + n_communities, original_n)) ; + GRB_TRY (GrB_mxm (new_node_map, NULL, NULL, GxB_ANY_PAIR_BOOL, + S_level, node_map, NULL)) ; + GRB_TRY (GrB_free (&node_map)) ; + GRB_TRY (GrB_free (&S_level)) ; + node_map = new_node_map ; + new_node_map = NULL ; + } + + // Prepare next level. + GRB_TRY (GrB_free (&C)) ; + I = NULL ; + A = A_next ; + A_next = NULL ; + free_A = true ; + + GRB_TRY (GrB_free (&it)) ; + GRB_TRY (GrB_free (&Ni)) ; + GRB_TRY (GrB_free (&desc)) ; + LG_TRY (LAGraph_Free ((void**) °ree, msg)) ; + LG_TRY (LAGraph_Free ((void**) &community_degree, msg)) ; + } + + //---------------------------------------------------------------------- + // set output + //---------------------------------------------------------------------- + + LG_ASSERT (C != NULL, GrB_NULL_POINTER) ; + GRB_TRY (GrB_Vector_new (com, GrB_UINT64, original_n)) ; + + if (node_map == NULL) + { + LG_ASSERT_MSG (final_nrows == original_n, GrB_INVALID_VALUE, + "unexpected shape mismatch without node aggregation") ; + GRB_TRY (GxB_Vector_load (*com, (void**) (&I), GrB_UINT64, original_n, + sizeof (uint64_t) * original_n, GrB_DEFAULT, NULL)) ; + } + else + { + // C_to_orig = C' * node_map: community rows, original-node columns + GRB_TRY (GrB_Matrix_new (&C_to_orig, GrB_BOOL, final_nrows, original_n)) ; + GRB_TRY (GrB_mxm (C_to_orig, NULL, NULL, GxB_ANY_PAIR_BOOL, + C, node_map, GrB_DESC_T0)) ; + + GrB_Index cnvals = 0 ; + GRB_TRY (GrB_Matrix_nvals (&cnvals, C_to_orig)) ; + LG_TRY (LAGraph_Malloc ((void**) &tuple_i, cnvals, sizeof (GrB_Index), msg)) ; + LG_TRY (LAGraph_Malloc ((void**) &tuple_j, cnvals, sizeof (GrB_Index), msg)) ; + LG_TRY (LAGraph_Malloc ((void**) &tuple_x, cnvals, sizeof (bool), msg)) ; + LG_TRY (LAGraph_Malloc ((void**) &com_i, cnvals, sizeof (GrB_Index), msg)) ; + LG_TRY (LAGraph_Malloc ((void**) &com_x, cnvals, sizeof (uint64_t), msg)) ; + + GRB_TRY (GrB_Matrix_extractTuples_BOOL ( + tuple_i, tuple_j, tuple_x, &cnvals, C_to_orig)) ; + + for (GrB_Index k = 0 ; k < cnvals ; k++) + { + com_i [k] = tuple_j [k] ; + com_x [k] = (uint64_t) tuple_i [k] ; + } + + GRB_TRY (GrB_Vector_build_UINT64 ( + *com, com_i, com_x, cnvals, GrB_SECOND_UINT64)) ; + } + + LG_FREE_ALL ; + return GrB_SUCCESS ; } - diff --git a/include/LAGraphX.h b/include/LAGraphX.h index 689df55ad9..e3c954f1ae 100644 --- a/include/LAGraphX.h +++ b/include/LAGraphX.h @@ -1493,6 +1493,17 @@ int LAGr_PartitionQuality( char *msg ); +LAGRAPHX_PUBLIC +GrB_Info LAGraph_louvain +( + GrB_Vector *com, // output communities + LAGraph_Graph G, // graph adjacency matrix + int itermax, // max number of modularity improvements sweeps per level + int levelmax, // max number of modularity improve and cluster condense + float e, // min change in modularity considered an improvement + char *msg // error message +) ; + LAGRAPHX_PUBLIC int LAGr_Modularity( // Outputs From a630bcb70e37365ac7cc2d502c163409224dafd2 Mon Sep 17 00:00:00 2001 From: Gabriel Gomez Date: Mon, 7 Sep 2026 13:10:53 -0500 Subject: [PATCH 13/13] use more graphblas functions for louvain --- experimental/algorithm/LAGraph_louvain.c | 72 +++++++++--------------- 1 file changed, 26 insertions(+), 46 deletions(-) diff --git a/experimental/algorithm/LAGraph_louvain.c b/experimental/algorithm/LAGraph_louvain.c index e5ac908686..4399444073 100644 --- a/experimental/algorithm/LAGraph_louvain.c +++ b/experimental/algorithm/LAGraph_louvain.c @@ -224,13 +224,9 @@ static GrB_Info LG_Louvain_aggregate GrB_free (&C_to_orig) ; \ GrB_free (&active_comms) ; \ GrB_free (&x_active) ; \ + GrB_free (&x_com) ; \ LAGraph_Free ((void**) °ree, NULL) ; \ LAGraph_Free ((void**) &community_degree, NULL); \ - LAGraph_Free ((void**) &tuple_i, NULL) ; \ - LAGraph_Free ((void**) &tuple_j, NULL) ; \ - LAGraph_Free ((void**) &tuple_x, NULL) ; \ - LAGraph_Free ((void**) &com_i, NULL) ; \ - LAGraph_Free ((void**) &com_x, NULL) ; \ LAGraph_Free ((void**) &d_i, NULL) ; \ LAGraph_Free ((void**) &d_x, NULL) ; \ } @@ -272,13 +268,10 @@ GrB_Info LAGraph_louvain GrB_Matrix C_to_orig = NULL ; GrB_Vector active_comms = NULL ; GrB_Vector x_active = NULL ; + GrB_Vector x_com = NULL ; // dummy dense vector for C_to_orig -> com uint64_t *degree = NULL ; // node degree for current level uint64_t *community_degree = NULL ; // sigma_tot per community id - bool *tuple_x = NULL ; // values for C_to_orig extraction - GrB_Index *tuple_i = NULL, *tuple_j = NULL ; - GrB_Index *com_i = NULL ; - uint64_t *com_x = NULL ; GrB_Index *d_i = NULL ; uint64_t *d_x = NULL ; @@ -290,6 +283,12 @@ GrB_Info LAGraph_louvain GrB_Matrix A = G->A ; // current adjacency over levels bool free_A = false ; // true once A becomes internally allocated + // variables for unloading vector D + GrB_Type unload_type = NULL ; + uint64_t unload_n = 0 ; + uint64_t unload_size = 0 ; + int unload_handling = 0 ; + // find out if graph is symmetric, compute cached values, and check loops LG_TRY (LAGraph_Cached_IsSymmetricStructure (G, msg)) ; LG_TRY (LAGraph_Cached_OutDegree (G, msg)) ; @@ -328,20 +327,13 @@ GrB_Info LAGraph_louvain GRB_TRY (GrB_reduce (D, NULL, GrB_PLUS_UINT64, GrB_PLUS_MONOID_UINT64, A, NULL)) ; - LG_TRY (LAGraph_Malloc ((void**) °ree, nrows, sizeof (uint64_t), msg)) ; - memset (degree, 0, nrows * sizeof (uint64_t)) ; + // Calculate M2 using GrB_reduce before unloading D + uint64_t M2_uint64 = 0 ; + GRB_TRY (GrB_reduce (&M2_uint64, NULL, GrB_PLUS_MONOID_UINT64, D, NULL)) ; - GrB_Index d_nvals = 0 ; - GRB_TRY (GrB_Vector_nvals (&d_nvals, D)) ; - LG_TRY (LAGraph_Malloc ((void**) &d_i, d_nvals, sizeof (GrB_Index), msg)) ; - LG_TRY (LAGraph_Malloc ((void**) &d_x, d_nvals, sizeof (uint64_t), msg)) ; - GRB_TRY (GrB_Vector_extractTuples_UINT64 (d_i, d_x, &d_nvals, D)) ; - for (GrB_Index k = 0 ; k < d_nvals ; k++) - { - degree [d_i [k]] = d_x [k] ; - } - LG_TRY (LAGraph_Free ((void**) &d_i, msg)) ; - LG_TRY (LAGraph_Free ((void**) &d_x, msg)) ; + // Unload vector D into the degree array + GRB_TRY (GxB_Vector_unload (D, (void **) °ree, &unload_type, + &unload_n, &unload_size, &unload_handling, NULL)) ; GRB_TRY (GrB_free (&D)) ; LG_TRY (LAGraph_Malloc ((void**) &community_degree, nrows, @@ -368,11 +360,8 @@ GrB_Info LAGraph_louvain double modularity_gain = 0 ; improved = true ; - double M2 = 0.0 ; - for (GrB_Index i = 0 ; i < nrows ; i++) - { - M2 += (double) degree [i] ; - } + // M2 is the sum of all degrees (computed via GrB_reduce above) + double M2 = (double) M2_uint64 ; LG_ASSERT_MSG (M2 > 0.0, GrB_INVALID_VALUE, "sum of degrees must be positive") ; @@ -527,25 +516,16 @@ GrB_Info LAGraph_louvain GRB_TRY (GrB_mxm (C_to_orig, NULL, NULL, GxB_ANY_PAIR_BOOL, C, node_map, GrB_DESC_T0)) ; - GrB_Index cnvals = 0 ; - GRB_TRY (GrB_Matrix_nvals (&cnvals, C_to_orig)) ; - LG_TRY (LAGraph_Malloc ((void**) &tuple_i, cnvals, sizeof (GrB_Index), msg)) ; - LG_TRY (LAGraph_Malloc ((void**) &tuple_j, cnvals, sizeof (GrB_Index), msg)) ; - LG_TRY (LAGraph_Malloc ((void**) &tuple_x, cnvals, sizeof (bool), msg)) ; - LG_TRY (LAGraph_Malloc ((void**) &com_i, cnvals, sizeof (GrB_Index), msg)) ; - LG_TRY (LAGraph_Malloc ((void**) &com_x, cnvals, sizeof (uint64_t), msg)) ; - - GRB_TRY (GrB_Matrix_extractTuples_BOOL ( - tuple_i, tuple_j, tuple_x, &cnvals, C_to_orig)) ; - - for (GrB_Index k = 0 ; k < cnvals ; k++) - { - com_i [k] = tuple_j [k] ; - com_x [k] = (uint64_t) tuple_i [k] ; - } - - GRB_TRY (GrB_Vector_build_UINT64 ( - *com, com_i, com_x, cnvals, GrB_SECOND_UINT64)) ; + // Each column of C_to_orig has exactly one nonzero, at the row index + // of the final community for that original node. Use a positional + // semiring to pull out that row index without ever leaving GrB. + GRB_TRY (GrB_Vector_new (&x_com, GrB_BOOL, final_nrows)) ; + GRB_TRY (GrB_assign (x_com, NULL, NULL, (bool) 0, GrB_ALL, final_nrows, + NULL)) ; + GRB_TRY (GrB_assign (*com, NULL, NULL, (uint64_t) 0, GrB_ALL, + original_n, NULL)) ; + GRB_TRY (GrB_vxm (*com, NULL, NULL, GxB_PLUS_FIRSTJ_INT64, + x_com, C_to_orig, NULL)) ; } LG_FREE_ALL ;