1818
1919// TODO: ready for src, except need vanilla, and the method is a bit slow
2020// because of the internal kernels it uses in SuiteSparse:GraphBLAS need work.
21+ // FIXED in v10.2 but may need more work in GraphBLAS.
2122
2223// LAGr_EdgeBetweennessCentrality: Exact algorithm for computing
2324// betweeness centrality.
2728//------------------------------------------------------------------------------
2829
2930#define useAssign
30- // #define debug
3131
3232#define LG_FREE_WORK \
3333{ \
@@ -94,7 +94,9 @@ int LAGr_EdgeBetweennessCentrality
9494 GrB_Matrix * centrality , // centrality(i): betweeness centrality of i
9595 // input:
9696 LAGraph_Graph G , // input graph
97- GrB_Vector sources , // source vertices to compute shortest paths (if NULL or empty, use all vertices)
97+ GrB_Vector sources , // source vertices to compute shortest paths
98+ // (if NULL or empty, use all vertices)
99+
98100 char * msg
99101)
100102{
@@ -217,11 +219,14 @@ int LAGr_EdgeBetweennessCentrality
217219 GRB_TRY (GrB_Vector_new (& bc_vertex_flow , GrB_FP64 , n )) ;
218220
219221 // Initialize centrality matrix with zeros using A as structural mask
222+ // centrality<A,struct> = 0
220223 LG_TRY (GrB_Matrix_new (centrality , GrB_FP64 , n , n )) ;
221- GRB_TRY (GrB_assign (* centrality , A , NULL , 0.0 , GrB_ALL , n , GrB_ALL , n , GrB_DESC_S )) ;
224+ GRB_TRY (GrB_assign (* centrality , A , NULL , 0.0 , GrB_ALL , n , GrB_ALL , n ,
225+ GrB_DESC_S )) ;
222226
223227 // Allocate memory for the array of S vectors
224- LG_TRY (LAGraph_Calloc ((void * * ) & Search , n + 1 , sizeof (GrB_Vector ), msg )) ;
228+ LG_TRY (LAGraph_Calloc ((void * * ) & Search , n + 1 , sizeof (GrB_Vector ),
229+ msg )) ;
225230
226231 // =========================================================================
227232 // === Process source nodes ================================================
@@ -234,7 +239,8 @@ int LAGr_EdgeBetweennessCentrality
234239 GRB_TRY (GrB_Vector_new (& internal_sources , GrB_INT64 , n )) ;
235240
236241 // internal_sources (0:n-1) = 0
237- GRB_TRY (GrB_assign (internal_sources , NULL , NULL , 0 , GrB_ALL , n , NULL )) ;
242+ GRB_TRY (GrB_assign (internal_sources , NULL , NULL , 0 , GrB_ALL , n ,
243+ NULL )) ;
238244
239245 // internal_sources (0:n-1) = 0:n-1
240246 GRB_TRY (GrB_apply (internal_sources , NULL , NULL , GrB_ROWINDEX_INT64 ,
@@ -266,7 +272,7 @@ int LAGr_EdgeBetweennessCentrality
266272 GRB_TRY (GrB_Vector_new (& J_vec , GrB_FP64 , n )) ;
267273 GRB_TRY (GrB_Vector_new (& I_vec , GrB_FP64 , n )) ;
268274 GRB_TRY (GrB_Matrix_new (& Fd1A , GrB_FP64 , n , n )) ;
269- GRB_TRY (GrB_Vector_new (& temp_update , GrB_FP64 , n )) ; // Create a temporary vector
275+ GRB_TRY (GrB_Vector_new (& temp_update , GrB_FP64 , n )) ; // Create temp vector
270276
271277 GRB_TRY (GrB_Matrix_new (& HalfUpdate , GrB_FP64 , n , n )) ;
272278 GRB_TRY (GrB_Matrix_new (& HalfUpdateT , GrB_FP64 , n , n )) ;
@@ -298,51 +304,54 @@ int LAGr_EdgeBetweennessCentrality
298304 GrB_DESC_T0 )) ;
299305
300306 GRB_TRY (GrB_Vector_nvals (& frontier_size , frontier )) ;
301- GRB_TRY (GrB_assign (frontier , frontier , NULL , 1.0 , GrB_ALL , n , GrB_DESC_S )) ;
307+ GRB_TRY (GrB_assign (frontier , frontier , NULL , 1.0 , GrB_ALL , n ,
308+ GrB_DESC_S )) ;
302309
303310 while (frontier_size != 0 )
304311 {
305312 depth ++ ;
306313
307- //----------------------------------------------------------------------
314+ //------------------------------------------------------------------
308315 // paths += frontier
309316 // Accumulate path counts for vertices at current depth
310- //----------------------------------------------------------------------
317+ //------------------------------------------------------------------
311318
312- GRB_TRY (GrB_assign (paths , NULL , GrB_PLUS_FP64 , frontier , GrB_ALL , n ,
313- NULL )) ;
319+ GRB_TRY (GrB_assign (paths , NULL , GrB_PLUS_FP64 , frontier , GrB_ALL ,
320+ n , NULL )) ;
314321
315- //----------------------------------------------------------------------
322+ //------------------------------------------------------------------
316323 // Search[depth] = structure(frontier)
317324 // Record the frontier structure at current depth
318- //----------------------------------------------------------------------
325+ //------------------------------------------------------------------
319326
320327 GrB_free (& (Search [depth ])) ;
321- LG_TRY (LAGraph_Vector_Structure (& (Search [depth ]), frontier , msg )) ;
328+ LG_TRY (LAGraph_Vector_Structure (& (Search [depth ]), frontier ,
329+ msg )) ;
322330
323- //----------------------------------------------------------------------
331+ //------------------------------------------------------------------
324332 // frontier<!paths> = frontier * A
325- //----------------------------------------------------------------------
333+ //------------------------------------------------------------------
326334
327335 GRB_TRY (LG_SET_FORMAT_HINT (frontier , LG_SPARSE )) ;
328336 GRB_TRY (GrB_vxm (frontier , paths , NULL , LAGraph_plus_first_fp64 ,
329337 frontier , A , GrB_DESC_RSC )) ;
330338
331- //----------------------------------------------------------------------
339+ //------------------------------------------------------------------
332340 // Get size of current frontier: frontier_size = nvals(frontier)
333- //----------------------------------------------------------------------
341+ //------------------------------------------------------------------
334342
335343 last_frontier_size = frontier_size ;
336344 GRB_TRY (GrB_Vector_nvals (& frontier_size , frontier )) ;
337- }
338345
346+ LG_ASSERT (depth <= n , GrB_PANIC ) ;
347+ }
339348
340- // =========================================================================
341- // === Betweenness centrality computation phase ============================
342- // =========================================================================
349+ // =====================================================================
350+ // === Betweenness centrality computation phase ========================
351+ // =====================================================================
343352
344353 // bc_vertex_flow = ones (n, n) ; a full matrix (and stays full)
345- GRB_TRY (GrB_assign (bc_vertex_flow , NULL , NULL , 0.0 , GrB_ALL , n , NULL )) ;
354+ GRB_TRY (GrB_assign (bc_vertex_flow , NULL , NULL , 0.0 , GrB_ALL , n , NULL ));
346355
347356 GRB_TRY (GrB_Matrix_clear (HalfUpdate )) ;
348357 GRB_TRY (GrB_Matrix_clear (HalfUpdateT )) ;
@@ -352,36 +361,39 @@ int LAGr_EdgeBetweennessCentrality
352361 GRB_TRY (GrB_Vector_clear (I_vec )) ;
353362 GRB_TRY (GrB_Vector_clear (temp_update )) ;
354363
355- // Backtrack through the BFS and compute centrality updates for each vertex
356- // GrB_Index fd1_size;
364+ // Backtrack through the BFS and compute centrality updates for each
365+ // vertex
357366
358367 while (depth >= 1 )
359368 {
360369 GrB_Vector f_d = Search [depth ] ;
361370 GrB_Vector f_d1 = Search [depth - 1 ] ;
362371
363- //----------------------------------------------------------------------
372+ //------------------------------------------------------------------
364373 // j<S(depth, :)> = (1 + v) / p
365374 // J = diag(j)
366375 // Compute weighted contributions from current level
367- //----------------------------------------------------------------------
376+ //------------------------------------------------------------------
377+
378+ GRB_TRY (GrB_eWiseMult (J_vec , f_d , NULL , Add_One_Divide ,
379+ bc_vertex_flow , paths , GrB_DESC_RS )) ;
368380
369- GRB_TRY (GrB_eWiseMult (J_vec , f_d , NULL , Add_One_Divide , bc_vertex_flow , paths , GrB_DESC_RS )) ;
370381 GRB_TRY (GrB_Matrix_diag (& J_matrix , J_vec , 0 )) ;
371382
372- //----------------------------------------------------------------------
383+ //------------------------------------------------------------------
373384 // i<S(depth-1, :)> = p
374385 // I = diag(i)
375386 // Compute weighted contributions from previous level
376- //----------------------------------------------------------------------
387+ //------------------------------------------------------------------
377388
378- GRB_TRY (GrB_Vector_extract (I_vec , f_d1 , NULL , paths , GrB_ALL , n , GrB_DESC_RS )) ;
389+ GRB_TRY (GrB_Vector_extract (I_vec , f_d1 , NULL , paths , GrB_ALL , n ,
390+ GrB_DESC_RS )) ;
379391 GRB_TRY (GrB_Matrix_diag (& I_matrix , I_vec , 0 )) ;
380392
381- //----------------------------------------------------------------------
393+ //------------------------------------------------------------------
382394 // Update = I × A × J
383395 // Compute edge updates based on current level weights
384- //----------------------------------------------------------------------
396+ //------------------------------------------------------------------
385397
386398 double t1 = LAGraph_WallClockTime ();
387399 GRB_TRY (GrB_mxm (Fd1A , NULL , NULL , LAGraph_plus_first_fp64 ,
@@ -396,32 +408,62 @@ int LAGr_EdgeBetweennessCentrality
396408 t2_total += t2 ;
397409 GRB_TRY (GrB_free (& I_matrix )) ;
398410 GRB_TRY (GrB_free (& J_matrix )) ;
399- //----------------------------------------------------------------------
400- // centrality<A> += Update
411+
412+ //------------------------------------------------------------------
413+ // centrality<A,struct> += Update
401414 // Accumulate centrality values for edges
402- //----------------------------------------------------------------------
415+ //------------------------------------------------------------------
416+
417+ // The centrality and A matrices have the same pattern.
418+ // centrality<centrality,struct> += Update should be fast, since it
419+ // tells GraphBLAS that the pattern of centrality doesn't change.
420+ // No pending tuples or zombies will be added to the matrix.
421+ // However, it currently forces GraphBLAS to make copy of the
422+ // centrality matrix to use as a mask.
403423
404424 #ifdef useAssign
405- // centrality{A} += Update, using assign
425+ // centrality<A,struct> += Update, using assign
406426 double t3 = LAGraph_WallClockTime ();
407427
408428 if (G -> kind == LAGraph_ADJACENCY_UNDIRECTED ) {
409- // First divide the Update matrix by 2 for symmetric distribution
410- GrB_apply (HalfUpdate , NULL , NULL , GrB_DIV_FP64 , Update , 2.0 , NULL );
429+ // First divide the Update matrix by 2 for symmetric
430+ // distribution
431+ GrB_apply (HalfUpdate , NULL , NULL , GrB_DIV_FP64 , Update ,
432+ 2.0 , NULL );
411433
412434 // Create a transposed version of the update
413435 GrB_transpose (HalfUpdateT , NULL , NULL , HalfUpdate , NULL );
414436
415- // Add the original and transposed matrices to create a symmetric update
416- GrB_eWiseAdd (SymmetricUpdate , NULL , NULL , GrB_PLUS_FP64 , HalfUpdate , HalfUpdateT , NULL );
437+ // Add the original and transposed matrices to create a
438+ // symmetric update
439+ GrB_eWiseAdd (SymmetricUpdate , NULL , NULL , GrB_PLUS_FP64 ,
440+ HalfUpdate , HalfUpdateT , NULL );
417441
418442 // Apply the symmetric update to the centrality
419- GRB_TRY (GrB_assign (* centrality , A , GrB_PLUS_FP64 , SymmetricUpdate , GrB_ALL , n , GrB_ALL , n , GrB_DESC_S ));
443+ #if LG_SUITESPARSE_GRAPHBLAS_V10_2
444+ // This was slow prior to v10.2.0:
445+ // centrality<centrality,struct> += SymmetricUpdate
446+ GRB_TRY (GrB_assign (* centrality , * centrality , GrB_PLUS_FP64 ,
447+ SymmetricUpdate , GrB_ALL , n , GrB_ALL , n , GrB_DESC_S ));
448+ #else
449+ // centrality<A,struct> += SymmetricUpdate
450+ GRB_TRY (GrB_assign (* centrality , A , GrB_PLUS_FP64 ,
451+ SymmetricUpdate , GrB_ALL , n , GrB_ALL , n , GrB_DESC_S ));
452+ #endif
420453
421454 }
422- else {
423- GRB_TRY (GrB_assign (* centrality , A , GrB_PLUS_FP64 , Update , GrB_ALL , n , GrB_ALL , n ,
424- GrB_DESC_S ));
455+ else
456+ {
457+ #if LG_SUITESPARSE_GRAPHBLAS_V10_2
458+ // This was slow prior to v10.2.0:
459+ // centrality<centrality,struct> += Update
460+ GRB_TRY (GrB_assign (* centrality , * centrality , GrB_PLUS_FP64 ,
461+ Update , GrB_ALL , n , GrB_ALL , n , GrB_DESC_S ));
462+ #else
463+ // centrality<A,struct> += Update
464+ GRB_TRY (GrB_assign (* centrality , A , GrB_PLUS_FP64 , Update ,
465+ GrB_ALL , n , GrB_ALL , n , GrB_DESC_S ));
466+ #endif
425467 }
426468
427469 t3 = LAGraph_WallClockTime () - t3 ;
@@ -430,20 +472,22 @@ int LAGr_EdgeBetweennessCentrality
430472 // TODO: Approx update using ewise add not implemented
431473 // centrality = centrality + Update using eWiseAdd
432474 double t3 = LAGraph_WallClockTime ();
433- GRB_TRY (GrB_eWiseAdd (* centrality , NULL , NULL , GrB_PLUS_FP64 , * centrality , Update , NULL ));
475+ GRB_TRY (GrB_eWiseAdd (* centrality , NULL , NULL , GrB_PLUS_FP64 ,
476+ * centrality , Update , NULL ));
434477 t3 = LAGraph_WallClockTime () - t3 ;
435478 t3_total += t3 ;
436479 #endif
437480
438481
439- //----------------------------------------------------------------------
482+ //------------------------------------------------------------------
440483 // v = Update +.
441484 // Reduce update matrix to vector for next iteration
442- //----------------------------------------------------------------------
443-
444- GRB_TRY (GrB_reduce (temp_update , NULL , NULL , GrB_PLUS_MONOID_FP64 , Update , NULL )) ;
445- GRB_TRY (GrB_eWiseAdd (bc_vertex_flow , NULL , NULL , GrB_PLUS_FP64 , bc_vertex_flow , temp_update , NULL )) ;
485+ //------------------------------------------------------------------
446486
487+ GRB_TRY (GrB_reduce (temp_update , NULL , NULL , GrB_PLUS_MONOID_FP64 ,
488+ Update , NULL )) ;
489+ GRB_TRY (GrB_eWiseAdd (bc_vertex_flow , NULL , NULL , GrB_PLUS_FP64 ,
490+ bc_vertex_flow , temp_update , NULL )) ;
447491 // 24 d = d − 1
448492 depth -- ;
449493 }
@@ -452,17 +496,13 @@ int LAGr_EdgeBetweennessCentrality
452496
453497 #ifdef debug
454498 printf (" I*A time: %g\n" , t1_total );
455-
456499 printf (" (I*A)*J time: %g\n" , t2_total );
457-
458500 #ifdef useAssign
459501 printf (" Centrality update using assign time: %g\n" , t3_total );
460502 #else
461503 printf (" Centrality update using eWiseAdd time: %g\n" , t3_total );
462504 #endif
463-
464- GxB_print (* centrality , GxB_FULL ) ;
465-
505+ // GxB_print (*centrality, GxB_FULL) ;
466506 #endif
467507
468508
0 commit comments