From 08e313e4f160f8e038616be74425bad10fbf32a6 Mon Sep 17 00:00:00 2001 From: Gabriel Gomez Date: Tue, 7 Apr 2026 14:36:43 -0500 Subject: [PATCH 1/3] RCC biop --- .../algorithm/LAGraph_RichClubCoefficient.c | 130 ++++++++++++++---- experimental/test/test_RichClubCoefficient.c | 5 +- 2 files changed, 103 insertions(+), 32 deletions(-) diff --git a/experimental/algorithm/LAGraph_RichClubCoefficient.c b/experimental/algorithm/LAGraph_RichClubCoefficient.c index d20ed9dac1..55216f6052 100755 --- a/experimental/algorithm/LAGraph_RichClubCoefficient.c +++ b/experimental/algorithm/LAGraph_RichClubCoefficient.c @@ -49,8 +49,14 @@ GrB_free(&iseq_2lt) ; \ GrB_free(&plus_2le) ; \ GrB_free(&rcCalculation) ; \ + GrB_free(&node_compare) ; \ GrB_free(&ramp_v) ; \ + GrB_free(&ctx_s) ; \ + GrB_free(&ctx_type) ; \ + GrB_free(&node_compare_idxop) ; \ + GrB_free(&node_compare) ; \ LAGraph_Free(&a_space, NULL) ; \ + LAGraph_Free((void **) &ctx.arr, NULL) ; \ } @@ -77,6 +83,29 @@ void LG_RCC_rich_club_formula(double *z, const int64_t *x, const int64_t *y) (*z) = ((double)(*x)) / (((double)(*y)) * (((double)(*y)) - 1.0)); }, LG_RCC_RICH_CLUB_FORMULA) +LG_JIT_STRING( +typedef struct { + int64_t *arr; +} LG_RCC_context ; +, LG_RCC_CONTEXT) + +LG_JIT_STRING( +void LG_RCC_compare_row_col ( + int64_t *z, + const bool *x, + GrB_Index ix, + GrB_Index jx, + const bool *y, + GrB_Index iy, + GrB_Index jy, + const LG_RCC_context *contx +) { + int64_t row_deg = contx->arr[ix] ; + int64_t col_deg = contx->arr[iy] ; + (*z) = (int64_t)((row_deg < col_deg) + (row_deg <= col_deg)) ; +} +, LG_RCC_COMPARE_ROW_COL) + int LAGraph_RichClubCoefficient ( // output: @@ -134,6 +163,15 @@ int LAGraph_RichClubCoefficient // 2E_K / (N_k (N_k -1)) GrB_BinaryOp rcCalculation = NULL; + // TODO + GrB_BinaryOp node_compare = NULL; + GxB_IndexBinaryOp node_compare_idxop = NULL; + + // context type + LG_RCC_context ctx = {0}; + GrB_Scalar ctx_s = NULL; + GrB_Type ctx_type = NULL; + GrB_Matrix A = NULL; // G->A, the adjacency matrix // Matrix used for row reduction @@ -147,9 +185,8 @@ int LAGraph_RichClubCoefficient void *a_space = NULL; - int64_t *node_edges_arr = NULL, *deg_arr = NULL, - *epd_arr = NULL, *ones = NULL, - *vpd_arr = NULL; + int64_t *node_edges_arr = NULL, *deg_arr = NULL, *epd_arr = NULL, + *ones = NULL, *vpd_arr = NULL; GrB_Type epd_type = NULL, vpd_type = NULL; uint64_t epd_n = 0, vpd_n = 0, epd_size = 0, vpd_size = 0; int epd_h = 0, vpd_h = 0; @@ -178,31 +215,33 @@ int LAGraph_RichClubCoefficient GRB_TRY (GrB_Vector_new(°rees, GrB_INT64, n)) ; GRB_TRY (GrB_Vector_new(&node_edges, GrB_INT64, n)) ; -#if LAGRAPH_SUITESPARSE - GRB_TRY (GxB_BinaryOp_new( - &iseq_2lt, (GxB_binary_function) (&LG_RCC_iseq_2islt), - GrB_INT64, GrB_INT64, GrB_INT64, "LG_RCC_iseq_2islt", -LG_RCC_ISEQ_2ISLT)) ; + +#if LG_SUITESPARSE_GRAPHBLAS_V10 + GRB_TRY (GxB_Type_new ( + &ctx_type, sizeof(LG_RCC_context), "LG_RCC_context", LG_RCC_CONTEXT)) ; + GRB_TRY (GxB_IndexBinaryOp_new ( + &node_compare_idxop, (GxB_index_binary_function) (LG_RCC_compare_row_col), + GrB_INT64, GrB_BOOL, GrB_BOOL, ctx_type, "LG_RCC_compare_row_col", + LG_RCC_COMPARE_ROW_COL)) ; GRB_TRY (GxB_BinaryOp_new( - &rcCalculation, (GxB_binary_function) (&LG_RCC_rich_club_formula), + &rcCalculation, (GxB_binary_function) (LG_RCC_rich_club_formula), GrB_FP64, GrB_INT64, GrB_INT64, "LG_RCC_rich_club_formula", LG_RCC_RICH_CLUB_FORMULA)) ; #else GRB_TRY (GrB_BinaryOp_new( - &iseq_2lt, (GxB_binary_function) (&LG_RCC_iseq_2islt), + &iseq_2lt, (GxB_binary_function) (LG_RCC_iseq_2islt), GrB_INT64, GrB_INT64, GrB_INT64)) ; GRB_TRY (GrB_BinaryOp_new( - &rcCalculation, (GxB_binary_function) (&LG_RCC_rich_club_formula), + &rcCalculation, (GxB_binary_function) (LG_RCC_rich_club_formula), GrB_FP64, GrB_INT64, GrB_INT64 )) ; #endif - GRB_TRY (GrB_Semiring_new(&plus_2le, GrB_PLUS_MONOID_INT64, iseq_2lt)) ; GRB_TRY (GrB_Vector_reduce_INT64( &max_deg, NULL, GrB_MAX_MONOID_INT64, G->out_degree, NULL)) ; - GRB_TRY (GrB_Vector_new(&edges_per_deg, GrB_INT64, max_deg)) ; - GRB_TRY (GrB_Vector_new(&verts_per_deg, GrB_INT64, max_deg)) ; - GRB_TRY (GrB_Vector_new(rccs, GrB_FP64, max_deg)) ; + GRB_TRY (GrB_Vector_new (&edges_per_deg, GrB_INT64, max_deg)) ; + GRB_TRY (GrB_Vector_new (&verts_per_deg, GrB_INT64, max_deg)) ; + GRB_TRY (GrB_Vector_new (rccs, GrB_FP64, max_deg)) ; //-------------------------------------------------------------------------- // Calculations @@ -211,28 +250,59 @@ LG_RCC_ISEQ_2ISLT)) ; // degrees = G->out_degree - 1 // Fill out degree vector, to target col_scale mxm on graphs // with singletons, scalar value irrelevant. - GRB_TRY (GrB_Vector_assign_INT64( + GRB_TRY (GrB_Vector_assign_INT64 ( degrees, NULL, NULL, (int64_t) -1, GrB_ALL, 0, NULL)) ; - GRB_TRY (GrB_Vector_assign( + GRB_TRY (GrB_Vector_assign ( degrees, NULL, GrB_PLUS_INT64, G->out_degree, GrB_ALL, 0, NULL)) ; - GRB_TRY (GrB_Matrix_diag(&D, degrees, 0)) ; - // Each edge in the graph gets the value of the degree of its row node + #if LAGRAPH_SUITESPARSE - GRB_TRY (GrB_mxm( - A_deg, NULL, NULL, GxB_ANY_FIRST_INT64, D, A, NULL)) ; + GrB_Type type; + uint64_t deg_n; + uint64_t deg_size; + int handling; + GRB_TRY (GxB_Vector_unload ( + degrees, (void **) &ctx.arr, &type, °_n, °_size, &handling, NULL + )) ; + GRB_TRY (GrB_Scalar_new (&ctx_s, ctx_type)); + GRB_TRY (GrB_Scalar_setElement_UDT (ctx_s, &ctx)); + GRB_TRY (GxB_BinaryOp_new_IndexOp ( + &node_compare, node_compare_idxop, ctx_s)); + GRB_TRY (GrB_Semiring_new(&plus_2le, GrB_PLUS_MONOID_INT64, node_compare)) ; + GRB_TRY (GrB_free (°rees)); + GRB_TRY (GrB_Vector_new (°rees, GrB_BOOL, deg_n)); + GRB_TRY (GrB_Vector_assign_INT64 ( + degrees, NULL, NULL, (bool) 0, GrB_ALL, 0, NULL)) ; + + // Sum the number of edges each node is "responsible" for. + GRB_TRY (GrB_mxv ( + node_edges, NULL, GrB_PLUS_INT64, plus_2le, A, degrees, NULL)) ; + + // The rest of this is indexing the number of edges and number of nodes at + // each degree and then doing a cummulative sum to know the amount of edges + // and nodes at degree geq k. + GRB_TRY (GrB_Vector_nvals (&edge_vec_nvals, node_edges)) ; + + GRB_TRY (GxB_Vector_load ( + degrees, (void **) &ctx.arr, type, deg_n, deg_size, handling, NULL + )) ; #else - GRB_TRY (GrB_mxm( + // Each edge in the graph gets the value of the degree of its row node + GRB_TRY (GrB_Matrix_diag (&D, degrees, 0)) ; + GRB_TRY (GrB_mxm ( A_deg, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_INT64, D, A, NULL)) ; -#endif + // Sum the number of edges each node is "responsible" for. - GRB_TRY (GrB_mxv( - node_edges, NULL, GrB_PLUS_INT64, plus_2le, A_deg, degrees, NULL)) ; + GRB_TRY (GrB_mxv ( + node_edges, NULL, GrB_PLUS_INT64, plus_2le, A, degrees, NULL)) ; // The rest of this is indexing the number of edges and number of nodes at // each degree and then doing a cummulative sum to know the amount of edges // and nodes at degree geq k. GRB_TRY (GrB_Vector_nvals (&edge_vec_nvals, node_edges)) ; +#endif + + #if LG_SUITESPARSE_GRAPHBLAS_V10 if(n == edge_vec_nvals) { @@ -242,12 +312,13 @@ LG_RCC_ISEQ_2ISLT)) ; node_edges = NULL; } else -{ + { GRB_TRY (GrB_Vector_assign( degrees, G->out_degree, NULL, degrees, GrB_ALL, 0, GrB_DESC_RS )) ; - GRB_TRY (GrB_Vector_new(°_x, GrB_BOOL, 0)) ; - GRB_TRY (GrB_Vector_new(&node_edges_x, GrB_BOOL, 0)) ; + + GRB_TRY (GrB_Vector_new(°_x, GrB_BOOL, 0)) ; + GRB_TRY (GrB_Vector_new(&node_edges_x, GrB_BOOL, 0)) ; GRB_TRY (GxB_Vector_extractTuples_Vector( NULL, deg_x, degrees, NULL )) ; @@ -256,7 +327,7 @@ LG_RCC_ISEQ_2ISLT)) ; )) ; } GRB_TRY (GrB_Vector_nvals(&edge_vec_nvals, node_edges_x)) - GRB_TRY (GrB_Vector_new(&ones_v, GrB_INT64, edge_vec_nvals)) ; + GRB_TRY (GrB_Vector_new(&ones_v, GrB_INT64, edge_vec_nvals)) ; GRB_TRY (GrB_Vector_assign_INT64( @@ -273,7 +344,6 @@ LG_RCC_ISEQ_2ISLT)) ; GRB_TRY (GrB_apply ( ramp_v, NULL, NULL, GrB_ROWINDEX_INT64, ramp_v, 0, NULL)) ; #endif - LG_TRY (LAGraph_FastAssign_Semiring ( edges_per_deg, NULL, GrB_PLUS_INT64, deg_x, node_edges_x, ramp_v, LAGraph_plus_second_int64, NULL, msg diff --git a/experimental/test/test_RichClubCoefficient.c b/experimental/test/test_RichClubCoefficient.c index 63866fb4d3..1a082e3853 100644 --- a/experimental/test/test_RichClubCoefficient.c +++ b/experimental/test/test_RichClubCoefficient.c @@ -111,7 +111,6 @@ void test_RichClubCoefficient (void) // start LAGraph //-------------------------------------------------------------------------- OK (LAGraph_Init (msg)) ; - for (int k = 0 ; ; k++) { @@ -170,8 +169,9 @@ void test_RichClubCoefficient (void) printf ("RCC computation begins:\n") ; GrB_set (GrB_GLOBAL, (int32_t) (true), GxB_BURBLE) ; - OK(LAGraph_RichClubCoefficient ( &rcc, G, msg)); + GrB_Info res = LAGraph_RichClubCoefficient ( &rcc, G, msg); printf("%s\n", msg); + OK(res); GrB_set (GrB_GLOBAL, (int32_t) (false), GxB_BURBLE) ; printf ("RCC computation ends:\n") ; @@ -200,6 +200,7 @@ void iseq(bool *z, const double *x, const double *y) { (*z) = (isnan(*x) && isnan(*y)) ||*x == *y ; } + //------------------------------------------------------------------------------ // test RichClubCoefficient vs C code //------------------------------------------------------------------------------ From ebde9ffc6fecd1d86fa3036e74b44b8c8dfb0442 Mon Sep 17 00:00:00 2001 From: Gabriel Gomez Date: Tue, 7 Apr 2026 15:30:23 -0500 Subject: [PATCH 2/3] solve issues for out_degree - 1 + // Fill out degree vector, to target col_scale mxm on graphs + // with singletons, scalar value irrelevant. + GRB_TRY (GrB_Vector_assign_INT64 ( + degrees, NULL, NULL, (int64_t) -1, GrB_ALL, 0, NULL)) ; + GRB_TRY (GrB_Vector_assign ( + degrees, NULL, GrB_PLUS_INT64, G->out_degree, GrB_ALL, 0, NULL)) ; + #if LG_SUITESPARSE_GRAPHBLAS_V10 GRB_TRY (GxB_Type_new ( &ctx_type, sizeof(LG_RCC_context), "LG_RCC_context", LG_RCC_CONTEXT)) ; @@ -227,48 +239,40 @@ int LAGraph_RichClubCoefficient &rcCalculation, (GxB_binary_function) (LG_RCC_rich_club_formula), GrB_FP64, GrB_INT64, GrB_INT64, "LG_RCC_rich_club_formula", LG_RCC_RICH_CLUB_FORMULA)) ; + + GrB_Type type; + uint64_t deg_n; + uint64_t deg_size; + int handling; + GRB_TRY (GxB_Vector_unload ( + degrees, (void **) &ctx.arr, &type, °_n, °_size, &handling, NULL + )) ; + GRB_TRY (GrB_Scalar_new (&ctx_s, ctx_type)); + GRB_TRY (GrB_Scalar_setElement_UDT (ctx_s, &ctx)); + GRB_TRY (GxB_BinaryOp_new_IndexOp ( + &node_compare, node_compare_idxop, ctx_s)); #else GRB_TRY (GrB_BinaryOp_new( - &iseq_2lt, (GxB_binary_function) (LG_RCC_iseq_2islt), + &node_compare, (GxB_binary_function) (LG_RCC_iseq_2islt), GrB_INT64, GrB_INT64, GrB_INT64)) ; GRB_TRY (GrB_BinaryOp_new( &rcCalculation, (GxB_binary_function) (LG_RCC_rich_club_formula), GrB_FP64, GrB_INT64, GrB_INT64 )) ; #endif + GRB_TRY (GrB_Semiring_new(&plus_2le, GrB_PLUS_MONOID_INT64, node_compare)) ; + //-------------------------------------------------------------------------- + // Calculations + //-------------------------------------------------------------------------- GRB_TRY (GrB_Vector_reduce_INT64( &max_deg, NULL, GrB_MAX_MONOID_INT64, G->out_degree, NULL)) ; GRB_TRY (GrB_Vector_new (&edges_per_deg, GrB_INT64, max_deg)) ; GRB_TRY (GrB_Vector_new (&verts_per_deg, GrB_INT64, max_deg)) ; GRB_TRY (GrB_Vector_new (rccs, GrB_FP64, max_deg)) ; - //-------------------------------------------------------------------------- - // Calculations - //-------------------------------------------------------------------------- - // degrees = G->out_degree - 1 - // Fill out degree vector, to target col_scale mxm on graphs - // with singletons, scalar value irrelevant. - GRB_TRY (GrB_Vector_assign_INT64 ( - degrees, NULL, NULL, (int64_t) -1, GrB_ALL, 0, NULL)) ; - GRB_TRY (GrB_Vector_assign ( - degrees, NULL, GrB_PLUS_INT64, G->out_degree, GrB_ALL, 0, NULL)) ; - - -#if LAGRAPH_SUITESPARSE - GrB_Type type; - uint64_t deg_n; - uint64_t deg_size; - int handling; - GRB_TRY (GxB_Vector_unload ( - degrees, (void **) &ctx.arr, &type, °_n, °_size, &handling, NULL - )) ; - GRB_TRY (GrB_Scalar_new (&ctx_s, ctx_type)); - GRB_TRY (GrB_Scalar_setElement_UDT (ctx_s, &ctx)); - GRB_TRY (GxB_BinaryOp_new_IndexOp ( - &node_compare, node_compare_idxop, ctx_s)); - GRB_TRY (GrB_Semiring_new(&plus_2le, GrB_PLUS_MONOID_INT64, node_compare)) ; +#if LG_SUITESPARSE_GRAPHBLAS_V10 GRB_TRY (GrB_free (°rees)); GRB_TRY (GrB_Vector_new (°rees, GrB_BOOL, deg_n)); GRB_TRY (GrB_Vector_assign_INT64 ( @@ -286,24 +290,7 @@ int LAGraph_RichClubCoefficient GRB_TRY (GxB_Vector_load ( degrees, (void **) &ctx.arr, type, deg_n, deg_size, handling, NULL )) ; -#else - // Each edge in the graph gets the value of the degree of its row node - GRB_TRY (GrB_Matrix_diag (&D, degrees, 0)) ; - GRB_TRY (GrB_mxm ( - A_deg, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_INT64, D, A, NULL)) ; - // Sum the number of edges each node is "responsible" for. - GRB_TRY (GrB_mxv ( - node_edges, NULL, GrB_PLUS_INT64, plus_2le, A, degrees, NULL)) ; - - // The rest of this is indexing the number of edges and number of nodes at - // each degree and then doing a cummulative sum to know the amount of edges - // and nodes at degree geq k. - GRB_TRY (GrB_Vector_nvals (&edge_vec_nvals, node_edges)) ; -#endif - - -#if LG_SUITESPARSE_GRAPHBLAS_V10 if(n == edge_vec_nvals) { deg_x = degrees; @@ -374,6 +361,20 @@ int LAGraph_RichClubCoefficient verts_per_deg, (void **) &vpd_arr, vpd_type, vpd_n, vpd_size, vpd_h, NULL)) ; #else + // Each edge in the graph gets the value of the degree of its row node + GRB_TRY (GrB_Matrix_diag (&D, degrees, 0)) ; + GRB_TRY (GrB_mxm ( + A_deg, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_INT64, D, A, NULL)) ; + + // Sum the number of edges each node is "responsible" for. + GRB_TRY (GrB_mxv ( + node_edges, NULL, GrB_PLUS_INT64, plus_2le, A_deg, degrees, NULL)) ; + + // The rest of this is indexing the number of edges and number of nodes at + // each degree and then doing a cummulative sum to know the amount of edges + // and nodes at degree geq k. + GRB_TRY (GrB_Vector_nvals (&edge_vec_nvals, node_edges)) ; + LG_TRY (LAGraph_Malloc( &a_space, edge_vec_nvals * 3 + max_deg * 4, sizeof(int64_t), NULL )) ; @@ -458,8 +459,7 @@ int LAGraph_RichClubCoefficient //Computes the RCC of a matrix GRB_TRY(GrB_eWiseMult( - *rccs, NULL, NULL, rcCalculation, - edges_per_deg, verts_per_deg, NULL + *rccs, NULL, NULL, rcCalculation, edges_per_deg, verts_per_deg, NULL )) ; LG_FREE_WORK ; From c64135101da83089fb88c24a5e27623a45394f9a Mon Sep 17 00:00:00 2001 From: Gabriel Gomez Date: Tue, 7 Apr 2026 16:13:17 -0500 Subject: [PATCH 3/3] split off vanilla function --- .../algorithm/LAGraph_RichClubCoefficient.c | 468 +++++++++++------- 1 file changed, 282 insertions(+), 186 deletions(-) diff --git a/experimental/algorithm/LAGraph_RichClubCoefficient.c b/experimental/algorithm/LAGraph_RichClubCoefficient.c index 97b59d0d8b..3dc92c6d89 100755 --- a/experimental/algorithm/LAGraph_RichClubCoefficient.c +++ b/experimental/algorithm/LAGraph_RichClubCoefficient.c @@ -29,45 +29,10 @@ // References: -// Julian J. McAuley, Luciano da Fontoura Costa, and Tibério S. Caetano, “The -// rich-club phenomenon across complex network hierarchies”, Applied Physics +// Julian J. McAuley, Luciano da Fontoura Costa, and Tibério S. Caetano, "The +// rich-club phenomenon across complex network hierarchies", Applied Physics // Letters Vol 91 Issue 8, August 2007. https://arxiv.org/abs/physics/0701290 -#define LG_FREE_WORK \ -{ \ - /* free any workspace used here */ \ - GrB_free(&D) ; \ - GrB_free(&P) ; \ - GrB_free(&A_deg) ; \ - GrB_free(°rees) ; \ - GrB_free(°_x) ; \ - GrB_free(&node_edges) ; \ - GrB_free(&node_edges_x) ; \ - GrB_free(&ones_v) ; \ - GrB_free(&edges_per_deg) ; \ - GrB_free(&verts_per_deg) ; \ - GrB_free(&iseq_2lt) ; \ - GrB_free(&plus_2le) ; \ - GrB_free(&rcCalculation) ; \ - GrB_free(&node_compare) ; \ - GrB_free(&ramp_v) ; \ - GrB_free(&ctx_s) ; \ - GrB_free(&ctx_type) ; \ - GrB_free(&node_compare_idxop) ; \ - GrB_free(&node_compare) ; \ - LAGraph_Free(&a_space, NULL) ; \ - LAGraph_Free((void **) &ctx.arr, NULL) ; \ -} - - -#define LG_FREE_ALL \ -{ \ - /* free any workspace used here */ \ - LG_FREE_WORK ; \ - /* free all the output variable(s) */ \ - GrB_free (rccs) ; \ -} - #include "LG_internal.h" #include "LAGraphX.h" @@ -110,125 +75,114 @@ void LG_RCC_compare_row_col ( typedef GrB_BinaryOp GxB_IndexBinaryOp; #endif -int LAGraph_RichClubCoefficient +//------------------------------------------------------------------------------ +// LG_RCC_compute_v10: SuiteSparse GraphBLAS v10+ implementation +//------------------------------------------------------------------------------ + +// Uses GxB_IndexBinaryOp (LG_RCC_compare_row_col) which reads the row and +// column node indices and looks up their degrees from a context array passed +// via a GrB_Scalar. This avoids building the intermediate D and A_deg +// matrices. The scatter of per-node edge counts into per-degree buckets uses +// LAGraph_FastAssign_Semiring. + +#undef LG_FREE_WORK +#define LG_FREE_WORK \ +{ \ + GrB_free(°rees) ; \ + GrB_free(°_x) ; \ + GrB_free(&node_edges) ; \ + GrB_free(&node_edges_x) ; \ + GrB_free(&ones_v) ; \ + GrB_free(&edges_per_deg) ; \ + GrB_free(&verts_per_deg) ; \ + GrB_free(&plus_2le) ; \ + GrB_free(&rcCalculation) ; \ + GrB_free(&node_compare) ; \ + GrB_free(&ramp_v) ; \ + GrB_free(&ctx_s) ; \ + GrB_free(&ctx_type) ; \ + GrB_free(&node_compare_idxop) ; \ + LAGraph_Free((void **) &ctx.arr, NULL) ; \ +} + +#undef LG_FREE_ALL +#define LG_FREE_ALL \ +{ \ + LG_FREE_WORK ; \ + GrB_free (rccs) ; \ +} + +#if LG_SUITESPARSE_GRAPHBLAS_V10 +static int LG_RCC_compute_v10 ( - // output: - //rccs(i): rich club coefficent of i GrB_Vector *rccs, - - // input: - LAGraph_Graph G, //input graph + LAGraph_Graph G, char *msg ) { //-------------------------------------------------------------------------- - // Declorations + // Declarations //-------------------------------------------------------------------------- - LG_CLEAR_MSG ; + GrB_Matrix A = NULL ; // G->A, the adjacency matrix - // n x n Adjacency Matrix - // With values cooresponding to the degree of its column - GrB_Matrix A_deg = NULL; + // n degrees vector (degrees = out_degree - 1) + GrB_Vector degrees = NULL, deg_x = NULL ; - // n x n Diagonal Matrix - // entries corresponding to degrees. - GrB_Matrix D = NULL; + // n x 1: number of edges each node is "responsible" for + // (smallest-degree end * 2 + same-degree edges, to handle double-counting) + GrB_Vector node_edges = NULL, node_edges_x = NULL ; - // n degrees vector - GrB_Vector degrees = NULL, deg_x = NULL; + // max_degree x 1: edge count / vertex count bucketed by degree + GrB_Vector edges_per_deg = NULL ; + GrB_Vector verts_per_deg = NULL ; - // n x 1 - // contains the number of edges for which the ith node is - // the smallest degree node * 2 + # edges w/ same degree as the other node - // to account for double counting of edges w/ same degree as the other node. - GrB_Vector node_edges = NULL, node_edges_x = NULL; + // edge_vec_nvals x 1: vector of ones + GrB_Vector ones_v = NULL ; - // max_degree x 1 - // the ith entry contains the number of edges whose lowest degree is i. - GrB_Vector edges_per_deg = NULL; + // ramp vector for LAGraph_FastAssign_Semiring + GrB_Vector ramp_v = NULL ; - // max_degree x 1 - // the ith entry contains the number of verticies whose degree is i. - GrB_Vector verts_per_deg = NULL; + // [+].[LG_RCC_compare_row_col]: compares row vs. col degree via index op + GrB_Semiring plus_2le = NULL ; - // edge_vec_nvals x 1 - // Vector of ones - GrB_Vector ones_v = NULL; + // 2E_K / (N_k (N_k - 1)) + GrB_BinaryOp rcCalculation = NULL ; - // Ramp vector - GrB_Vector ramp_v = NULL; + // index-binary-op wrapper carrying the degree context + GrB_BinaryOp node_compare = NULL ; + GxB_IndexBinaryOp node_compare_idxop = NULL ; - // 2 * (x < y) + (x == y) - GrB_BinaryOp iseq_2lt = NULL; - - // [+].[iseq_2lt] - GrB_Semiring plus_2le = NULL; - - // 2E_K / (N_k (N_k -1)) - GrB_BinaryOp rcCalculation = NULL; - - // TODO - GrB_BinaryOp node_compare = NULL; - GxB_IndexBinaryOp node_compare_idxop = NULL; - - // context type - LG_RCC_context ctx = {0}; - GrB_Scalar ctx_s = NULL; - GrB_Type ctx_type = NULL; - - GrB_Matrix A = NULL; // G->A, the adjacency matrix - - // Matrix used for row reduction - GrB_Matrix P = NULL; + // context: holds pointer to the raw degree array + LG_RCC_context ctx = {0} ; + GrB_Scalar ctx_s = NULL ; + GrB_Type ctx_type = NULL ; GrB_Index n ; + GrB_Index edge_vec_nvals ; + int64_t max_deg ; - GrB_Index edge_vec_nvals; - int64_t max_deg; - bool iso = false; - - void *a_space = NULL; - - int64_t *node_edges_arr = NULL, *deg_arr = NULL, *epd_arr = NULL, - *ones = NULL, *vpd_arr = NULL; - GrB_Type epd_type = NULL, vpd_type = NULL; - uint64_t epd_n = 0, vpd_n = 0, epd_size = 0, vpd_size = 0; - int epd_h = 0, vpd_h = 0; - GrB_Index *epd_index = NULL, *vpd_index = NULL; - - //-------------------------------------------------------------------------- - // Check inputs - //-------------------------------------------------------------------------- - LG_TRY (LAGraph_CheckGraph (G, msg)) ; - LG_ASSERT (rccs != NULL, GrB_NULL_POINTER); - - LG_ASSERT_MSG( - G->kind == LAGraph_ADJACENCY_UNDIRECTED, GrB_INVALID_VALUE, - "G->A must be symmetric") ; - LG_ASSERT_MSG (G->out_degree != NULL, GrB_EMPTY_OBJECT, - "G->out_degree must be defined") ; - LG_ASSERT_MSG (G->nself_edges == 0, GrB_INVALID_VALUE, - "G->nself_edges must be zero") ; + GrB_Type epd_type = NULL, vpd_type = NULL ; + uint64_t epd_n = 0, vpd_n = 0, epd_size = 0, vpd_size = 0 ; + int epd_h = 0, vpd_h = 0 ; + int64_t *epd_arr = NULL, *vpd_arr = NULL ; //-------------------------------------------------------------------------- // Initializations //-------------------------------------------------------------------------- A = G->A ; GRB_TRY (GrB_Matrix_nrows (&n, A)) ; - GRB_TRY (GrB_Matrix_new(&A_deg, GrB_INT64,n,n)) ; - GRB_TRY (GrB_Vector_new(°rees, GrB_INT64, n)) ; - GRB_TRY (GrB_Vector_new(&node_edges, GrB_INT64, n)) ; + GRB_TRY (GrB_Vector_new (°rees, GrB_INT64, n)) ; + GRB_TRY (GrB_Vector_new (&node_edges, GrB_INT64, n)) ; // degrees = G->out_degree - 1 - // Fill out degree vector, to target col_scale mxm on graphs + // Fill out degree vector, to target col_scale mxm on graphs // with singletons, scalar value irrelevant. GRB_TRY (GrB_Vector_assign_INT64 ( degrees, NULL, NULL, (int64_t) -1, GrB_ALL, 0, NULL)) ; GRB_TRY (GrB_Vector_assign ( degrees, NULL, GrB_PLUS_INT64, G->out_degree, GrB_ALL, 0, NULL)) ; -#if LG_SUITESPARSE_GRAPHBLAS_V10 GRB_TRY (GxB_Type_new ( &ctx_type, sizeof(LG_RCC_context), "LG_RCC_context", LG_RCC_CONTEXT)) ; GRB_TRY (GxB_IndexBinaryOp_new ( @@ -240,27 +194,19 @@ int LAGraph_RichClubCoefficient GrB_FP64, GrB_INT64, GrB_INT64, "LG_RCC_rich_club_formula", LG_RCC_RICH_CLUB_FORMULA)) ; - GrB_Type type; - uint64_t deg_n; - uint64_t deg_size; - int handling; + GrB_Type type ; + uint64_t deg_n ; + uint64_t deg_size ; + int handling ; GRB_TRY (GxB_Vector_unload ( degrees, (void **) &ctx.arr, &type, °_n, °_size, &handling, NULL )) ; - GRB_TRY (GrB_Scalar_new (&ctx_s, ctx_type)); - GRB_TRY (GrB_Scalar_setElement_UDT (ctx_s, &ctx)); + GRB_TRY (GrB_Scalar_new (&ctx_s, ctx_type)) ; + GRB_TRY (GrB_Scalar_setElement_UDT (ctx_s, &ctx)) ; GRB_TRY (GxB_BinaryOp_new_IndexOp ( - &node_compare, node_compare_idxop, ctx_s)); -#else - GRB_TRY (GrB_BinaryOp_new( - &node_compare, (GxB_binary_function) (LG_RCC_iseq_2islt), - GrB_INT64, GrB_INT64, GrB_INT64)) ; - GRB_TRY (GrB_BinaryOp_new( - &rcCalculation, (GxB_binary_function) (LG_RCC_rich_club_formula), - GrB_FP64, GrB_INT64, GrB_INT64 )) ; -#endif + &node_compare, node_compare_idxop, ctx_s)) ; - GRB_TRY (GrB_Semiring_new(&plus_2le, GrB_PLUS_MONOID_INT64, node_compare)) ; + GRB_TRY (GrB_Semiring_new (&plus_2le, GrB_PLUS_MONOID_INT64, node_compare)) ; //-------------------------------------------------------------------------- // Calculations @@ -271,10 +217,8 @@ int LAGraph_RichClubCoefficient GRB_TRY (GrB_Vector_new (&verts_per_deg, GrB_INT64, max_deg)) ; GRB_TRY (GrB_Vector_new (rccs, GrB_FP64, max_deg)) ; - -#if LG_SUITESPARSE_GRAPHBLAS_V10 - GRB_TRY (GrB_free (°rees)); - GRB_TRY (GrB_Vector_new (°rees, GrB_BOOL, deg_n)); + GRB_TRY (GrB_free (°rees)) ; + GRB_TRY (GrB_Vector_new (°rees, GrB_BOOL, deg_n)) ; GRB_TRY (GrB_Vector_assign_INT64 ( degrees, NULL, NULL, (bool) 0, GrB_ALL, 0, NULL)) ; @@ -293,10 +237,10 @@ int LAGraph_RichClubCoefficient if(n == edge_vec_nvals) { - deg_x = degrees; - degrees = NULL; - node_edges_x = node_edges; - node_edges = NULL; + deg_x = degrees ; + degrees = NULL ; + node_edges_x = node_edges ; + node_edges = NULL ; } else { @@ -304,8 +248,8 @@ int LAGraph_RichClubCoefficient degrees, G->out_degree, NULL, degrees, GrB_ALL, 0, GrB_DESC_RS )) ; - GRB_TRY (GrB_Vector_new(°_x, GrB_BOOL, 0)) ; - GRB_TRY (GrB_Vector_new(&node_edges_x, GrB_BOOL, 0)) ; + GRB_TRY (GrB_Vector_new (°_x, GrB_BOOL, 0)) ; + GRB_TRY (GrB_Vector_new (&node_edges_x, GrB_BOOL, 0)) ; GRB_TRY (GxB_Vector_extractTuples_Vector( NULL, deg_x, degrees, NULL )) ; @@ -313,9 +257,8 @@ int LAGraph_RichClubCoefficient NULL, node_edges_x, node_edges, NULL )) ; } - GRB_TRY (GrB_Vector_nvals(&edge_vec_nvals, node_edges_x)) - GRB_TRY (GrB_Vector_new(&ones_v, GrB_INT64, edge_vec_nvals)) ; - + GRB_TRY (GrB_Vector_nvals (&edge_vec_nvals, node_edges_x)) + GRB_TRY (GrB_Vector_new (&ones_v, GrB_INT64, edge_vec_nvals)) ; GRB_TRY (GrB_Vector_assign_INT64( edges_per_deg, NULL, NULL, (int64_t) 0, GrB_ALL, 0, NULL)) ; @@ -325,12 +268,13 @@ int LAGraph_RichClubCoefficient ones_v, NULL, NULL, (int64_t) 0, GrB_ALL, 0, NULL)) ; #ifndef COVERAGE - GRB_TRY (GrB_Vector_new(&ramp_v, GrB_INT64, edge_vec_nvals + 1)) ; + GRB_TRY (GrB_Vector_new (&ramp_v, GrB_INT64, edge_vec_nvals + 1)) ; GRB_TRY (GrB_Vector_assign_INT64( ramp_v, NULL, NULL, (int64_t) 0, GrB_ALL, 0, NULL)) ; GRB_TRY (GrB_apply ( ramp_v, NULL, NULL, GrB_ROWINDEX_INT64, ramp_v, 0, NULL)) ; #endif + LG_TRY (LAGraph_FastAssign_Semiring ( edges_per_deg, NULL, GrB_PLUS_INT64, deg_x, node_edges_x, ramp_v, LAGraph_plus_second_int64, NULL, msg @@ -348,7 +292,8 @@ int LAGraph_RichClubCoefficient &vpd_n, &vpd_size, &vpd_h, NULL)) ; LG_ASSERT (max_deg == vpd_n && max_deg == epd_n, GrB_INVALID_VALUE) ; - //run a cummulative sum (backwards) on vpd_arr + // TODO: should be a GraphBLAS cumulative-sum primitive + // run a cummulative sum (backwards) on vpd_arr for(int64_t i = max_deg - 1; i > 0; --i) { vpd_arr[i-1] += vpd_arr[i] ; @@ -360,7 +305,134 @@ int LAGraph_RichClubCoefficient GRB_TRY(GxB_Vector_load( verts_per_deg, (void **) &vpd_arr, vpd_type, vpd_n, vpd_size, vpd_h, NULL)) ; -#else + + // Computes the RCC of a matrix + GRB_TRY(GrB_eWiseMult( + *rccs, NULL, NULL, rcCalculation, edges_per_deg, verts_per_deg, NULL + )) ; + + LG_FREE_WORK ; + return (GrB_SUCCESS) ; +} +#endif + +//------------------------------------------------------------------------------ +// LG_RCC_compute_vanilla: vanilla GraphBLAS / pre-v10 implementation +//------------------------------------------------------------------------------ + +// Cannot use index binary ops, so materializes row degrees into A_deg = D * A +// (via GrB_Matrix_diag + mxm), where each entry A_deg(i,j) = out_degree(i)-1. +// Uses the value-only binary op iseq_2lt for mxv. The scatter of per-node +// edge counts into per-degree buckets uses manual C-array extraction followed +// by GrB_Vector_build. + +#undef LG_FREE_WORK +#define LG_FREE_WORK \ +{ \ + GrB_free(&D) ; \ + GrB_free(&A_deg) ; \ + GrB_free(°rees) ; \ + GrB_free(&node_edges) ; \ + GrB_free(&edges_per_deg) ; \ + GrB_free(&verts_per_deg) ; \ + GrB_free(&plus_2le) ; \ + GrB_free(&rcCalculation) ; \ + GrB_free(&node_compare) ; \ + LAGraph_Free(&a_space, NULL) ; \ +} + +#undef LG_FREE_ALL +#define LG_FREE_ALL \ +{ \ + LG_FREE_WORK ; \ + GrB_free (rccs) ; \ +} + +static int LG_RCC_compute_vanilla +( + GrB_Vector *rccs, + LAGraph_Graph G, + char *msg +) +{ + //-------------------------------------------------------------------------- + // Declarations + //-------------------------------------------------------------------------- + GrB_Matrix A = NULL ; // G->A, the adjacency matrix + + // n x n: adjacency matrix with A_deg(i,j) = out_degree(i) - 1 + GrB_Matrix A_deg = NULL ; + + // n x n diagonal matrix with D(i,i) = out_degree(i) - 1 + GrB_Matrix D = NULL ; + + // n degrees vector (degrees = out_degree - 1) + GrB_Vector degrees = NULL ; + + // n x 1: number of edges each node is "responsible" for + // (smallest-degree end * 2 + same-degree edges, to handle double-counting) + GrB_Vector node_edges = NULL ; + + // max_degree x 1: edge count / vertex count bucketed by degree + GrB_Vector edges_per_deg = NULL ; + GrB_Vector verts_per_deg = NULL ; + + // [+].[iseq_2lt]: compares values directly (row degree stored in A_deg) + GrB_Semiring plus_2le = NULL ; + + // 2E_K / (N_k (N_k - 1)) + GrB_BinaryOp rcCalculation = NULL ; + + // 2 * (x < y) + (x == y) + GrB_BinaryOp node_compare = NULL ; + + GrB_Index n ; + GrB_Index edge_vec_nvals ; + int64_t max_deg ; + bool iso = false ; + + void *a_space = NULL ; + + int64_t *node_edges_arr = NULL, *deg_arr = NULL, *epd_arr = NULL, + *ones = NULL, *vpd_arr = NULL ; + GrB_Index *epd_index = NULL, *vpd_index = NULL ; + + //-------------------------------------------------------------------------- + // Initializations + //-------------------------------------------------------------------------- + A = G->A ; + GRB_TRY (GrB_Matrix_nrows (&n, A)) ; + GRB_TRY (GrB_Matrix_new (&A_deg, GrB_INT64, n, n)) ; + + GRB_TRY (GrB_Vector_new (°rees, GrB_INT64, n)) ; + GRB_TRY (GrB_Vector_new (&node_edges, GrB_INT64, n)) ; + + // degrees = G->out_degree - 1 + // Fill out degree vector, to target col_scale mxm on graphs + // with singletons, scalar value irrelevant. + GRB_TRY (GrB_Vector_assign_INT64 ( + degrees, NULL, NULL, (int64_t) -1, GrB_ALL, 0, NULL)) ; + GRB_TRY (GrB_Vector_assign ( + degrees, NULL, GrB_PLUS_INT64, G->out_degree, GrB_ALL, 0, NULL)) ; + + GRB_TRY (GrB_BinaryOp_new( + &node_compare, (GxB_binary_function) (LG_RCC_iseq_2islt), + GrB_INT64, GrB_INT64, GrB_INT64)) ; + GRB_TRY (GrB_BinaryOp_new( + &rcCalculation, (GxB_binary_function) (LG_RCC_rich_club_formula), + GrB_FP64, GrB_INT64, GrB_INT64 )) ; + + GRB_TRY (GrB_Semiring_new (&plus_2le, GrB_PLUS_MONOID_INT64, node_compare)) ; + + //-------------------------------------------------------------------------- + // Calculations + //-------------------------------------------------------------------------- + GRB_TRY (GrB_Vector_reduce_INT64( + &max_deg, NULL, GrB_MAX_MONOID_INT64, G->out_degree, NULL)) ; + GRB_TRY (GrB_Vector_new (&edges_per_deg, GrB_INT64, max_deg)) ; + GRB_TRY (GrB_Vector_new (&verts_per_deg, GrB_INT64, max_deg)) ; + GRB_TRY (GrB_Vector_new (rccs, GrB_FP64, max_deg)) ; + // Each edge in the graph gets the value of the degree of its row node GRB_TRY (GrB_Matrix_diag (&D, degrees, 0)) ; GRB_TRY (GrB_mxm ( @@ -378,19 +450,19 @@ int LAGraph_RichClubCoefficient LG_TRY (LAGraph_Malloc( &a_space, edge_vec_nvals * 3 + max_deg * 4, sizeof(int64_t), NULL )) ; - int64_t *T = a_space; - deg_arr = T; T += edge_vec_nvals; - node_edges_arr = T; T += edge_vec_nvals; - ones = T; T += edge_vec_nvals; - epd_arr = T; T += max_deg; - vpd_arr = T; T += max_deg; - epd_index = T; T += max_deg; - vpd_index = T; T += max_deg; + int64_t *T = a_space ; + deg_arr = T ; T += edge_vec_nvals ; + node_edges_arr = T ; T += edge_vec_nvals ; + ones = T ; T += edge_vec_nvals ; + epd_arr = T ; T += max_deg ; + vpd_arr = T ; T += max_deg ; + epd_index = T ; T += max_deg ; + vpd_index = T ; T += max_deg ; #pragma omp parallel for schedule(static) for(uint64_t i = 0; i < edge_vec_nvals; ++i) { - ones[i] = 1ll; + ones[i] = 1ll ; } GRB_TRY (GrB_Vector_apply_BinaryOp2nd_INT64( degrees, NULL, NULL, GrB_MINUS_INT64, G->out_degree, 1, NULL)) ; @@ -404,15 +476,15 @@ int LAGraph_RichClubCoefficient // Build with degrees as indecies and handle duplicates via adition GRB_TRY (GrB_Vector_build_INT64 ( - edges_per_deg, deg_arr, node_edges_arr, edge_vec_nvals, + edges_per_deg, deg_arr, node_edges_arr, edge_vec_nvals, GrB_PLUS_INT64)) ; GRB_TRY (GrB_Vector_build_INT64 ( verts_per_deg, deg_arr, ones, edge_vec_nvals, GrB_PLUS_INT64)) ; GRB_TRY (GrB_Vector_assign_INT64( - edges_per_deg, edges_per_deg, NULL, (int64_t) 0, + edges_per_deg, edges_per_deg, NULL, (int64_t) 0, GrB_ALL, 0, GrB_DESC_SC)) ; GRB_TRY (GrB_Vector_assign_INT64( - verts_per_deg, verts_per_deg, NULL, (int64_t) 0, + verts_per_deg, verts_per_deg, NULL, (int64_t) 0, GrB_ALL, 0, GrB_DESC_SC)) ; // Extract into arrays @@ -422,7 +494,8 @@ int LAGraph_RichClubCoefficient GRB_TRY (GrB_Vector_extractTuples_INT64( vpd_index, vpd_arr, &max_deg, verts_per_deg )) ; - //run a cummulative sum (backwards) on vpd_arr + // TODO: should be a GraphBLAS cumulative-sum primitive + // run a cummulative sum (backwards) on vpd_arr for(int64_t i = max_deg - 1; i > 0; --i) { vpd_arr[i-1] += vpd_arr[i] ; @@ -438,26 +511,8 @@ int LAGraph_RichClubCoefficient )) ; T = deg_arr = node_edges_arr = ones = NULL ; epd_index = vpd_index = epd_arr = vpd_arr = NULL ; -#endif - /** - * Cumulative sum (TODO: should be a GBLAS method!) - * - * GrB_cumsum(GrB_Matrix C, const GrB_Matrix mask, const GrB_BinaryOp accum, - * const GrB_BinaryOp plus, GrB_Matrix A, const GrB_Descriptor desc) - * - * By default sums rows. Returns a nearly full matrix: - * [., ., 1, 1, 1, 1, ., ., 1] --> [., ., 1, 2, 3, 4, 4, 4, 5] - * Mask can be A, then returns a matrix with the same pattern. - * [., ., 1, 1, 1, 1, ., ., 1] --> [., ., 1, 2, 3, 4, ., ., 5] - * - * Should we be able to sum in the opposite direction? - * Yes since not all monoids have inverse operations. - * - * If plus biop is not a monoid, this method should still work? - */ - - //Computes the RCC of a matrix + // Computes the RCC of a matrix GRB_TRY(GrB_eWiseMult( *rccs, NULL, NULL, rcCalculation, edges_per_deg, verts_per_deg, NULL )) ; @@ -465,3 +520,44 @@ int LAGraph_RichClubCoefficient LG_FREE_WORK ; return (GrB_SUCCESS) ; } + +//------------------------------------------------------------------------------ +// LAGraph_RichClubCoefficient: public dispatcher +//------------------------------------------------------------------------------ + +#undef LG_FREE_WORK +#define LG_FREE_WORK { } + +#undef LG_FREE_ALL +#define LG_FREE_ALL { } + +int LAGraph_RichClubCoefficient +( + // output: + //rccs(i): rich club coefficent of i + GrB_Vector *rccs, + + // input: + LAGraph_Graph G, //input graph + char *msg +) +{ + LG_CLEAR_MSG ; + + LG_TRY (LAGraph_CheckGraph (G, msg)) ; + LG_ASSERT (rccs != NULL, GrB_NULL_POINTER) ; + + LG_ASSERT_MSG( + G->kind == LAGraph_ADJACENCY_UNDIRECTED, GrB_INVALID_VALUE, + "G->A must be symmetric") ; + LG_ASSERT_MSG (G->out_degree != NULL, GrB_EMPTY_OBJECT, + "G->out_degree must be defined") ; + LG_ASSERT_MSG (G->nself_edges == 0, GrB_INVALID_VALUE, + "G->nself_edges must be zero") ; + +#if LG_SUITESPARSE_GRAPHBLAS_V10 + return LG_RCC_compute_v10 (rccs, G, msg) ; +#else + return LG_RCC_compute_vanilla(rccs, G, msg) ; +#endif +}