Min cut dev

experimental/algorithm/LAGr_MaxFlow.c

@@ -612,6 +612,7 @@ int LAGr_MaxFlow
     // output:
     double *f,              // max flow from src node to sink node
     GrB_Matrix *flow_mtx,   // optional output flow matrix
+    GrB_Matrix *res_mtx,    // optional output for min cut
     // input:
     LAGraph_Graph G,        // graph to compute maxflow on
     GrB_Index src,          // source node
@@ -1160,6 +1161,12 @@ int LAGr_MaxFlow
     GRB_TRY(GrB_select(*flow_mtx, NULL, NULL, GrB_VALUEGT_FP64, *flow_mtx, 0, NULL));
   }
 
+  if(res_mtx != NULL){
+   GRB_TRY(GrB_apply(*res_mtx, NULL, NULL, GetResidual, R, NULL)) ;
+    // prune zeros and negative entries from R_hat
+   GRB_TRY(GrB_select(*res_mtx, NULL, NULL, GrB_VALUEGT_FP64, *res_mtx, 0, NULL)) ; 
+  }
+
   //----------------------------------------------------------------------------
   // for test coverage only
   //----------------------------------------------------------------------------

experimental/algorithm/LAGraph_MinCut.c

@@ -0,0 +1,63 @@
+#include <LAGraph.h>
+#include "LG_internal.h"
+#include <LAGraph.h>
+
+#undef LG_FREE_ALL
+#undef LG_FREE_WORK
+
+#define LG_FREE_WORK				\
+{						\
+ GrB_free(&S_diag);				\
+ GrB_free(&S_bar_diag);				\
+ LAGraph_Delete(&G, msg);			\
+}
+
+#define LG_FREE_ALL                             \
+  { LG_FREE_WORK }
+
+int LAGraph_MinCut
+(
+    // outputs
+    GrB_Vector* S,
+    GrB_Vector* S_bar,
+    GrB_Matrix* cut_set,
+    // inputs
+    LAGraph_Graph G_origin,
+    GrB_Matrix R,
+    GrB_Index s,
+    GrB_Index t,
+    char *msg
+)
+{
+  //do a bfs from the source to the sink, stop if the frontier is empty
+
+  LAGraph_Graph G = NULL;
+  GrB_Matrix S_diag=NULL, S_bar_diag=NULL, A = G_origin->A;
+  GrB_Index n = 0;
+  GrB_Matrix_nrows(&n, R);
+  if(cut_set == NULL){
+    GRB_TRY(GrB_Matrix_new(cut_set, GrB_INT64, n, n));
+  }
+
+  LG_TRY(LAGraph_New(&G, &R, LAGraph_ADJACENCY_DIRECTED, msg));
+  LG_TRY(LAGr_BreadthFirstSearch(S, NULL, G, s, msg));
+
+  //restore the saturated edges capacities to get the set of cut edges
+  //LG_TRY(GrB_assign(R, R, NULL, A, GrB_ALL, n, GrB_ALL, n, GrB_DESC_SC)); //also fails
+
+  LG_TRY(GrB_assign(*S_bar, *S, NULL, 1, GrB_ALL, n, GrB_DESC_SC));
+  LG_TRY(GrB_assign(*S, *S, NULL, 1, GrB_ALL, n, GrB_DESC_S));
+
+  //GxB_print(G_origin->A, 5);
+  GRB_TRY(GrB_Matrix_diag(&S_diag, *S, 0));
+  GRB_TRY(GrB_Matrix_diag(&S_bar_diag, *S_bar, 0));
+
+
+  GRB_TRY(GrB_mxm(*cut_set, NULL, NULL, GxB_PLUS_TIMES_INT64, G_origin->A, S_bar_diag, NULL));
+  //GxB_print(*cut_set, 5);
+  GRB_TRY(GrB_mxm(*cut_set, NULL, NULL, GxB_PLUS_TIMES_INT64, S_diag, *cut_set, NULL));
+
+
+  LG_FREE_ALL;
+  return (GrB_SUCCESS);
+}

experimental/benchmark/maxflow_demo.c

@@ -58,15 +58,15 @@ int main (int argc, char ** argv){
 
   // LG_SET_BURBLE(1);
   double time = LAGraph_WallClockTime();
-  LAGRAPH_TRY(LAGr_MaxFlow(&flow, NULL, G, S, T, msg));
+  LAGRAPH_TRY(LAGr_MaxFlow(&flow, NULL, 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 %" PRIu64 " to %" PRIu64
     ", with flow_matrix returned\n", S, T);
   time = LAGraph_WallClockTime();
-  LAGRAPH_TRY(LAGr_MaxFlow(&flow, &flow_matrix, G, S, T, msg));
+  LAGRAPH_TRY(LAGr_MaxFlow(&flow, &flow_matrix, NULL, 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);

experimental/test/test_MaxFlow.c

@@ -76,14 +76,14 @@ void test_MaxFlow(void) {
     // test with JIT
     OK(GxB_Global_Option_set(GxB_JIT_C_CONTROL, GxB_JIT_ON));
     double flow = 0;
-    OK(LAGr_MaxFlow(&flow, NULL, G, tests[test].S, tests[test].T, msg));
+    OK(LAGr_MaxFlow(&flow, NULL, NULL, G, tests[test].S, tests[test].T, msg));
     printf("%s\n", msg);
     printf("flow is: %lf\n", flow);
     TEST_CHECK(flow == tests[test].F);
 
     // test without JIT
     OK(GxB_Global_Option_set(GxB_JIT_C_CONTROL, GxB_JIT_OFF));
-    OK(LAGr_MaxFlow(&flow, NULL, G, tests[test].S, tests[test].T, msg));
+    OK(LAGr_MaxFlow(&flow, NULL, NULL, G, tests[test].S, tests[test].T, msg));
     TEST_CHECK(flow == tests[test].F);
     OK(GxB_Global_Option_set(GxB_JIT_C_CONTROL, GxB_JIT_ON));
 
@@ -99,7 +99,7 @@ void test_MaxFlow(void) {
             {
               printf("src: %d, dest: %d\n", (int) src, (int) dest);
                 if (src == dest) continue ;
-                OK(LAGr_MaxFlow(&flow, NULL, G, src, dest, msg));
+                OK(LAGr_MaxFlow(&flow, NULL, NULL, G, src, dest, msg));
             }
         }
     }
@@ -144,7 +144,7 @@ void test_MaxFlowMtx(void) {
     // test with JIT
     OK(GxB_Global_Option_set(GxB_JIT_C_CONTROL, GxB_JIT_ON));
     double flow = 0;
-    OK(LAGr_MaxFlow(&flow, &flow_mtx, G, tests[test].S, tests[test].T, msg));
+    OK(LAGr_MaxFlow(&flow, &flow_mtx, NULL, G, tests[test].S, tests[test].T, msg));
     TEST_CHECK (flow_mtx != NULL) ;
     GxB_print (flow_mtx, 2) ;
     int status = LG_check_flow(flow_mtx, msg);
@@ -157,7 +157,7 @@ void test_MaxFlowMtx(void) {
 
     // test without JIT
     OK(GxB_Global_Option_set(GxB_JIT_C_CONTROL, GxB_JIT_OFF));
-    OK(LAGr_MaxFlow(&flow, &flow_mtx, G, tests[test].S, tests[test].T, msg));
+    OK(LAGr_MaxFlow(&flow, &flow_mtx, NULL, G, tests[test].S, tests[test].T, msg));
     TEST_CHECK (flow_mtx != NULL) ;
     status = LG_check_flow(flow_mtx, msg);
     TEST_CHECK (status == GrB_SUCCESS) ;

experimental/test/test_MinCut.c

@@ -0,0 +1,101 @@
+#include <stdio.h>
+#include <acutest.h>
+#include <LAGraphX.h>
+#include <LAGraph_test.h>
+#include <LG_Xtest.h>
+
+
+char msg[LAGRAPH_MSG_LEN];
+LAGraph_Graph G = NULL;
+GrB_Matrix A = NULL;
+#define LEN 512
+#define NTESTS 7
+char filename[LEN + 1];
+
+
+typedef struct {
+  char* filename;
+  GrB_Index s;
+  GrB_Index t;
+  LAGraph_Kind kind;
+} test_case;
+
+test_case tests[] = {
+  {"wiki.mtx", 0, 5, LAGraph_ADJACENCY_DIRECTED},
+  {"matrix_random_flow.mtx", 0,9, LAGraph_ADJACENCY_DIRECTED},
+  {"rand.mtx", 0, 19, LAGraph_ADJACENCY_DIRECTED},
+  {"mcl.mtx", 0, 9, LAGraph_ADJACENCY_DIRECTED},
+  {"cycle_flow.mtx", 0, 89,  LAGraph_ADJACENCY_DIRECTED},
+  {"random_weighted_general2.mtx", 0, 299, LAGraph_ADJACENCY_UNDIRECTED},
+  {"random_weighted_general1.mtx", 0, 499, LAGraph_ADJACENCY_UNDIRECTED}
+};
+
+
+void test_MinCut() {
+
+  LAGraph_Init(msg);
+  //OK(LG_SET_BURBLE(1));
+  OK(LG_SET_BURBLE(0));
+
+  for(uint8_t test = 0; test < NTESTS; test++){
+    GrB_Matrix A=NULL, R=NULL, cut_set=NULL;
+    GrB_Vector S=NULL, S_bar=NULL;
+    double min_cut = 0;
+    GrB_Index n = 0;
+    printf ("\nMatrix: %s\n", tests[test].filename);
+    TEST_CASE(tests[test].filename);
+
+    snprintf(filename, LEN, LG_DATA_DIR "%s", tests[test].filename);
+    FILE* f = fopen(filename, "r");
+    TEST_CHECK(f != NULL);
+
+    OK(LAGraph_MMRead(&A, f, msg));
+    OK(GrB_Matrix_nrows(&n, A));
+    OK(GrB_Matrix_new(&R, GrB_INT64, n, n));
+    OK(GrB_Matrix_new(&cut_set, GrB_INT64, n, n));
+
+    OK(GrB_Vector_new(&S, GrB_INT64, n));
+    OK(GrB_Vector_new(&S_bar, GrB_INT64, n));
+
+    OK(fclose(f));
+    LAGraph_Kind kind = tests [test].kind ;
+    OK(LAGraph_New(&G, &A, kind, msg));
+    if (kind == LAGraph_ADJACENCY_DIRECTED)
+    {
+        OK(LAGraph_Cached_AT(G, msg));
+    }
+
+    OK(LAGraph_Cached_EMin(G, msg));
+
+    // test with JIT
+    OK(GxB_Global_Option_set(GxB_JIT_C_CONTROL, GxB_JIT_ON));
+    double flow = 0;
+    OK(LAGr_MaxFlow(&flow, NULL, &R, G, tests[test].s, tests[test].t, msg));
+    printf("%s\n", msg);
+    printf("flow is: %lf\n", flow);
+
+    OK(LAGraph_MinCut(&S, &S_bar, &cut_set, G, R, tests[test].s, tests[test].t, msg));
+    printf("%s\n", msg);
+
+    //GxB_print(S, 5);
+    //GxB_print(S_bar, 5);
+    //GxB_print(cut_set, 5);
+
+    OK(GrB_reduce(&min_cut, NULL, GrB_PLUS_MONOID_INT64, cut_set, NULL));
+
+    TEST_CHECK(flow == min_cut);
+    printf("The min cut: %lf\n", min_cut);
+
+    GrB_free(&A);
+    GrB_free(&R);
+    GrB_free(&S);
+    GrB_free(&S_bar);
+    GrB_free(&cut_set);
+    LAGraph_Delete(&G, msg);
+  }
+
+  LAGraph_Finalize(msg);
+
+}
+
+TEST_LIST = {{"MinCut", test_MinCut}, {NULL, NULL}};

include/LAGraphX.h

@@ -1507,6 +1507,7 @@ int LAGr_MaxFlow(
     //outputs
     double* f,
     GrB_Matrix* flow_mtx,
+    GrB_Matrix* res_mtx,
     //inputs
     LAGraph_Graph G,
     GrB_Index src, //source node index
@@ -1515,6 +1516,19 @@ int LAGr_MaxFlow(
     char* msg
 );
 
+LAGRAPHX_PUBLIC
+int LAGraph_MinCut(
+    //outputs
+    GrB_Vector* S,
+    GrB_Vector* S_bar,
+    GrB_Matrix* cut_set,
+    // inputs
+    LAGraph_Graph G_origin, //original graph with capacities
+    GrB_Matrix R, //residual graph
+    GrB_Index s, //source node index
+    GrB_Index t, //sink node index
+    char *msg
+);
 
 #if defined ( __cplusplus )
 }