performance bug fix for HPEC25 paper on maxflow

Author: DrTimothyAldenDavis

ChangeLog

@@ -1,3 +1,7 @@
+Sept 1, 2025: version 1.2.1
+
+    * MaxFlow: performance bug fix for results in HPEC'25 paper.
+
 July 25, 2025: version 1.2.0
 
     * v2.1 C API: using the new GrB_get/set methods in the C API;

experimental/algorithm/LAGr_MaxFlow.c

@@ -15,6 +15,8 @@
 
 //------------------------------------------------------------------------------
 
+// TODO: work in progress (flow_matrix and global relabel)
+
 // LAGr_MaxFlow is a GraphBLAS implementation of the push-relabel algorithm
 // of Baumstark et al. [1]
 //
@@ -27,6 +29,7 @@
 #include "LG_internal.h"
 #include <LAGraph.h>
 
+//#define DBG
 #if LG_SUITESPARSE_GRAPHBLAS_V10
 
 //------------------------------------------------------------------------------
@@ -90,6 +93,8 @@ static GrB_Info LG_global_relabel
     GrB_Index sink,             // sink node
     GrB_Vector src_and_sink,    // mask vector, with just [src sink]
     GrB_UnaryOp GetResidual,    // unary op to compute resid=capacity-flow
+    GrB_BinaryOp global_relabel_accum, // accum for the disconnected assign
+    GrB_Index relabel_value, // value to relabel the disconnected nodes
     // input/output:
     GrB_Vector d,       // d(i) = height/label of node i
     // outputs:
@@ -116,10 +121,10 @@ static GrB_Info LG_global_relabel
     LG_TRY(LAGraph_Cached_OutDegree(G2, msg));
     // compute lvl using bfs on G2, starting at sink node
     LG_TRY(LAGr_BreadthFirstSearch(lvl, NULL, G2, sink, msg));
-    // d<!struct([src,sink])> = lvl
-    GRB_TRY(GrB_assign(d, src_and_sink, NULL, *lvl, GrB_ALL, n, GrB_DESC_SC));
-    // d<!struct(lvl)> = n
-    GRB_TRY(GrB_assign(d, *lvl, NULL, n, GrB_ALL, n, GrB_DESC_SC));
+    // d<!struct([src,sink])> = max(d(i), lvl)
+    GRB_TRY(GrB_assign(d, src_and_sink, global_relabel_accum, *lvl, GrB_ALL, n, GrB_DESC_SC));
+    // d<!struct(lvl)> = max(d(i), relabel_value)
+    GRB_TRY(GrB_assign(d, *lvl, global_relabel_accum, relabel_value, GrB_ALL, n, GrB_DESC_SC));
     LG_FREE_WORK ;
     return (GrB_SUCCESS) ;
 }
@@ -459,7 +464,7 @@ JIT_STR(void LG_MF_MxeMult32(LG_MF_resultTuple32 * z,
   bool j_active = ((*y) > 0) ;
   if ((x->di <  x->dj-1) /* case a */
   ||  (x->di == x->dj-1 && !j_active) /* case b */
-  ||  (x->di == x->dj   && (!j_active || (j_active && (i < j)))) /* case c */
+  ||  (x->di == x->dj && (!j_active || (j_active && (i < j)))) /* case c */
   ||  (x->di == x->dj+1))   /* case d */
   {
       z->residual = x->residual;
@@ -564,6 +569,28 @@ JIT_STR(void LG_MF_ExtractMatrixFlow(double* flow, const LG_MF_flowEdge* edge){*
 
 #endif
 
+#ifdef DBG
+void print_compareVec(const GrB_Vector vec) {
+  GxB_Iterator iter;
+  GxB_Iterator_new(&iter);
+  GrB_Info info = GxB_Vector_Iterator_attach(iter, vec, NULL);
+  if(info < 0){
+    printf("error with matrix passed in");
+  }
+  info = GxB_Vector_Iterator_seek(iter, 0);
+  while(info != GxB_EXHAUSTED){
+    GrB_Index i;
+    i = GxB_Vector_Iterator_getIndex(iter);
+    LG_MF_compareTuple32 e;
+    GxB_Iterator_get_UDT(iter, &e);
+    printf("(%ld, 0)         (di: %d, dj: %d, J: %d, residual: %lf) \n", i, e.di, e.dj, e.j, e.residual);
+    info = GxB_Vector_Iterator_next(iter);
+  }
+  GrB_free(&iter);
+}
+#endif
+
+
 //------------------------------------------------------------------------------
 // LAGraph_MaxFlow
 //------------------------------------------------------------------------------
@@ -766,6 +793,9 @@ int LAGr_MaxFlow
 
   GrB_Type Integer_Type = NULL ;
 
+  //accum operator for the global relabel
+  GrB_BinaryOp global_relabel_accum = NULL ;
+
   #ifdef COVERAGE
   // Just for test coverage, use 64-bit ints for n > 100.  Do not use this
   // rule in production!
@@ -782,6 +812,9 @@ int LAGr_MaxFlow
 
     Integer_Type = GrB_INT64 ;
 
+    // use the 64 bit max operator 
+    global_relabel_accum = GrB_MAX_INT64 ;
+
     // create types for computation
     GRB_TRY(GxB_Type_new(&ResultTuple, sizeof(LG_MF_resultTuple64),
         "LG_MF_resultTuple64", RESULTTUPLE_STR64));
@@ -852,6 +885,9 @@ int LAGr_MaxFlow
 
     Integer_Type = GrB_INT32 ;
 
+    // use 32 bit max op
+    global_relabel_accum = GrB_MAX_INT32 ;
+
     // create types for computation
     GRB_TRY(GxB_Type_new(&ResultTuple, sizeof(LG_MF_resultTuple32),
         "LG_MF_resultTuple32", RESULTTUPLE_STR32));
@@ -943,8 +979,14 @@ int LAGr_MaxFlow
   GRB_TRY(GrB_apply(R, A, NULL, ResidualBackward, AT, GrB_DESC_SC));
 
   // initial global relabeling
-  LG_TRY (LG_global_relabel (R, sink, src_and_sink, GetResidual, d, &lvl, msg)) ;
-
+  // relabel to 2*n to prevent any flow from going to the
+  // disconnected nodes.
+  GrB_Index relabel_value = 2*n ;
+  LG_TRY (LG_global_relabel (R, sink, src_and_sink, GetResidual, global_relabel_accum, relabel_value, d, &lvl, msg)) ;
+
+  // reset to n
+  relabel_value = n ;
+  
   // create excess vector e and initial flows from the src to its neighbors
   // e<struct(lvl)> = A (src,:)
   GRB_TRY(GrB_Vector_new(&e, GrB_FP64, n));
@@ -968,20 +1010,27 @@ int LAGr_MaxFlow
 
   for (int64_t iter = 0 ; n_active > 0 ; iter++)
   {
-
+    #ifdef DBG
+      printf ("iter: %ld, n_active %ld\n", iter, n_active) ;
+    #endif
     //--------------------------------------------------------------------------
     // Part 1: global relabeling
     //--------------------------------------------------------------------------
 
-    if ((iter > 0) && (flow_mtx != NULL) && (iter % 12 == 0))
+    if ((iter > 0)  && (iter % 12 == 0))
     {
-      LG_TRY (LG_global_relabel (R, sink, src_and_sink, GetResidual, d, &lvl, msg)) ;
-      // delete nodes in e that cannot be reached from the sink
-      // e<!struct(lvl)> = empty scalar
-      GrB_assign (e, lvl, NULL, empty, GrB_ALL, n, GrB_DESC_SC) ;
+      #ifdef DBG
+        printf ("relabel at : %ld\n", iter) ;
+      #endif
+      LG_TRY (LG_global_relabel (R, sink, src_and_sink, GetResidual, global_relabel_accum, relabel_value, d, &lvl, msg)) ;
+      if(flow_mtx == NULL){
+        // delete nodes in e that cannot be reached from the sink
+	//  e<!struct(lvl)> = empty scalar
+	GrB_assign (e, lvl, NULL, empty, GrB_ALL, n, GrB_DESC_SC) ;
+	GRB_TRY(GrB_Vector_nvals(&n_active, e));
+	if(n_active == 0) break;
+      }
       GrB_free(&lvl);
-      GRB_TRY(GrB_Vector_nvals(&n_active, e));
-      if(n_active == 0) break;
     }
 
     //--------------------------------------------------------------------------
@@ -1001,6 +1050,14 @@ int LAGr_MaxFlow
     // create Map matrix from pattern and values of yd
     // yd = CreateCompareVec (y,d) using eWiseMult
     GRB_TRY(GrB_eWiseMult(yd, NULL, NULL, CreateCompareVec, y,  d, NULL));
+
+    #ifdef DBG
+      print_compareVec(yd);
+      GxB_print(d, 3);
+      GxB_print(e, 3);
+    #endif
+    
+    
     // Jvec = ExtractJ (yd), where Jvec(i) = yd(i)->j
     GRB_TRY(GrB_apply(Jvec, NULL, NULL, ExtractJ, yd, NULL));
     GRB_TRY(GrB_Matrix_clear(Map));
@@ -1043,7 +1100,12 @@ int LAGr_MaxFlow
     // create the Delta matrix from delta_vec and y
     // note that delta_vec has the same structure as y
     // Jvec = ExtractYJ (y), where Jvec(i) = y(i)->j
+    // if Jvec has no values, then there is no possible
+    // candidates to push to, so the algorithm terminates
     GRB_TRY(GrB_apply(Jvec, NULL, NULL, ExtractYJ, y, NULL));
+    GrB_Index J_n;
+    GRB_TRY(GrB_Vector_nvals(&J_n, Jvec));
+    if(J_n == 0) break;
     GRB_TRY(GrB_Matrix_clear(Delta));
     GRB_TRY(GrB_Matrix_build(Delta, delta_vec, Jvec, delta_vec, GxB_IGNORE_DUP, desc));
 
@@ -1057,14 +1119,15 @@ int LAGr_MaxFlow
 
     // reduce Delta to delta_vec
     // delta_vec = sum (Delta), summing up each row of Delta
-    GRB_TRY(GrB_reduce(delta_vec, NULL, NULL, GrB_PLUS_FP64, Delta, GrB_DESC_T0));
+    GRB_TRY(GrB_reduce(delta_vec, NULL, NULL, GrB_PLUS_MONOID_FP64, Delta, GrB_DESC_T0));
 
     // add delta_vec to e
     // e<struct(delta_vec)> += delta_vec
     GRB_TRY(GrB_assign(e, delta_vec, GrB_PLUS_FP64, delta_vec, GrB_ALL, n, GrB_DESC_S));
 
     // augment maxflow for all active nodes
     LG_TRY (LG_augment_maxflow (f, e, sink, src_and_sink, &n_active, msg)) ;
+
   }
 
   //----------------------------------------------------------------------------

experimental/benchmark/do_maxflow

@@ -0,0 +1,14 @@
+#!/bin/bash
+
+
+MTX="/home/davis/matrices"
+GAP="/home/davis/GAP"
+
+../../build/experimental/benchmark/maxflow_demo $MTX/com-Youtube.mtx 4 275
+../../build/experimental/benchmark/maxflow_demo $MTX/com-LiveJournal.mtx 2 88
+../../build/experimental/benchmark/maxflow_demo $MTX/com-Orkut.mtx 43 75
+
+../../build/experimental/benchmark/maxflow_demo $GAP/GAP-web/GAP-web.grb 64712 500
+../../build/experimental/benchmark/maxflow_demo $GAP/GAP-road/GAP-road.grb 4 1075
+../../build/experimental/benchmark/maxflow_demo $GAP/GAP-twitter/GAP-twitter.grb 13 17
+

experimental/benchmark/maxflow_demo.c

@@ -7,6 +7,7 @@
 
 #define LAGRAPH_CATCH(info)                                                    \
   {                                                                            \
+    printf ("result: %d %s\n", info, msg) ;                                    \
     LAGraph_Delete(&G, msg);                                                   \
     return (info);                                                             \
   }
@@ -19,37 +20,58 @@ int main (int argc, char ** argv){
 
   char msg[LAGRAPH_MSG_LEN];
   LAGraph_Graph G = NULL;
-  
+
   double flow = 0;
-  GrB_Index T=0, S=0;
+  GrB_Index T=0, S=0, nflow ;
+  GrB_Matrix flow_matrix = NULL ;
 
   LAGRAPH_TRY(LAGraph_Init(msg));
-  
+
   //read in graph
   LAGRAPH_TRY(readproblem(&G, NULL, false, true, false, NULL, true, argc, argv));
+  printf ("read in the problem\n") ;
+  double t1 = LAGraph_WallClockTime();
   LAGRAPH_TRY(LAGraph_Cached_AT(G, msg));
+  t1 = LAGraph_WallClockTime() - t1 ;
+  printf ("cached AT, time %g\n", t1) ;
+  t1 = LAGraph_WallClockTime();
   LAGRAPH_TRY(LAGraph_Cached_EMin(G, msg));
+  t1 = LAGraph_WallClockTime() - t1 ;
+  printf ("cached EMin, time %g\n", t1) ;
 
   char* end1, *end2;
   S = strtoul(argv[2], &end1, 10);
   T = strtoul(argv[3], &end2, 10);
   if(argc > 4){
     int num_threads = atoi(argv[4]);
-    LAGRAPH_TRY(LAGraph_SetNumThreads(1, num_threads, msg)); 
+    LAGRAPH_TRY(LAGraph_SetNumThreads(1, num_threads, msg));
   }
 
+  int nthreads_outer, nthreads_inner ;
+  LAGRAPH_TRY(LAGraph_GetNumThreads(&nthreads_outer, &nthreads_inner, msg));
+  printf ("nthreads: %d %d\n", nthreads_outer, nthreads_inner) ;
+
   if(end1 == 0 || end2 == 0){
     printf("values for source and sink are incorrect.\n");
   }
-  printf("Starting max flow from %ld to %ld", S, T);
+  printf("Starting max flow from %ld to %ld\n", S, T);
 
-  //LG_SET_BURBLE(1);
+  // LG_SET_BURBLE(1);
   double time = LAGraph_WallClockTime();
   LAGRAPH_TRY(LAGr_MaxFlow(&flow, NULL, G, S, T, msg));
   time = LAGraph_WallClockTime() - time;
   printf("Time for LAGraph_MaxFlow: %g sec\n", time);
   printf("Max Flow is: %lf\n", flow);
- 
+
+  /* printf("Starting max flow from %ld to %ld, with flow_matrix returned\n", S, T); */
+  /* time = LAGraph_WallClockTime(); */
+  /* LAGRAPH_TRY(LAGr_MaxFlow(&flow, &flow_matrix, G, S, T, msg)); */
+  /* time = LAGraph_WallClockTime() - time; */
+  /* printf("Time for LAGraph_MaxFlow with flow matrix: %g sec\n", time); */
+  /* printf("Max Flow is: %lf\n", flow); */
+  /* GRB_TRY (GrB_Matrix_nvals (&nflow, flow_matrix)) ; */
+  /* printf("# of entries in flow matrix: %lu\n", nflow); */
+
   LAGraph_Delete(&G, msg);
   LAGRAPH_TRY(LAGraph_Finalize(msg));
 

experimental/test/test_MaxFlow.c

@@ -53,7 +53,7 @@ test_info tests[] = {
 void test_MaxFlow(void) {
 #if LG_SUITESPARSE_GRAPHBLAS_V10
   LAGraph_Init(msg);
-//OK(LG_SET_BURBLE(1));
+  //OK(LG_SET_BURBLE(1));
   OK(LG_SET_BURBLE(0));
   for(uint8_t test = 0; test < NTESTS; test++){
     GrB_Matrix A=NULL;
@@ -87,6 +87,23 @@ void test_MaxFlow(void) {
     TEST_CHECK(flow == tests[test].F);
     OK(GxB_Global_Option_set(GxB_JIT_C_CONTROL, GxB_JIT_ON));
 
+    // test all source/destination pairs for small problems
+    GrB_Index n ;
+    OK (GrB_Matrix_nrows (&n, G->A)) ;
+    printf ("n: %ld\n", (int64_t) n) ;
+    if (n < 100)
+    {
+        for (GrB_Index src = 0 ; src < n ; src++)
+        {
+            for (GrB_Index dest = 0 ; dest < n ; dest++)
+            {
+              printf("src: %ld, dest: %ld\n", src, dest);
+                if (src == dest) continue ;
+                OK(LAGr_MaxFlow(&flow, NULL, G, src, dest, msg));
+            }
+        }
+    }
+
     //free work
     OK(LAGraph_Delete(&G, msg));
   }