diff --git a/data/test_zero_cap.mtx b/data/test_zero_cap.mtx new file mode 100644 index 0000000000..8137244111 --- /dev/null +++ b/data/test_zero_cap.mtx @@ -0,0 +1,13 @@ +%%MatrixMarket matrix coordinate real general +% Example weighted directed graph adjacency (1-indexed) +% All weights are positive; some edges explicitly have zero weight +5 5 9 +1 2 0.500000 +1 4 0.000000 +2 3 1.250000 +2 5 2.000000 +3 1 0.000000 +3 4 3.141590 +4 2 0.750000 +5 4 0.125000 +5 5 1.000000 diff --git a/experimental/algorithm/LAGr_MaxFlow.c b/experimental/algorithm/LAGr_MaxFlow.c index 87cc1533b9..12deb4b517 100644 --- a/experimental/algorithm/LAGr_MaxFlow.c +++ b/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 the residual matrix. Used in the min-cut // input: LAGraph_Graph G, // graph to compute maxflow on GrB_Index src, // source node @@ -711,6 +712,10 @@ int LAGr_MaxFlow { (*flow_mtx) = NULL ; } + if (res_mtx != NULL) + { + (*res_mtx) = NULL ; + } LG_TRY(LAGraph_CheckGraph(G, msg)); LG_ASSERT (f != NULL, GrB_NULL_POINTER) ; (*f) = 0; @@ -720,6 +725,7 @@ int LAGr_MaxFlow LG_ASSERT_MSG(nrows == n, GrB_INVALID_VALUE, "Matrix must be square"); LG_ASSERT_MSG(src < n && src >= 0 && sink < n && sink >= 0, GrB_INVALID_VALUE, "src and sink must be a value between [0, n)"); + LG_ASSERT(G->emin_state != LAGraph_BOOLEAN_UNKNOWN, GrB_UNINITIALIZED_OBJECT) ; LG_ASSERT_MSG(G->emin > 0, GrB_INVALID_VALUE, "the edge weights (capacities) must be greater than 0"); @@ -1160,6 +1166,13 @@ 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_Matrix_new(res_mtx, GrB_FP64, n, n)); + 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 //---------------------------------------------------------------------------- diff --git a/experimental/algorithm/LAGraph_MinCut.c b/experimental/algorithm/LAGraph_MinCut.c new file mode 100644 index 0000000000..bc60f8de0d --- /dev/null +++ b/experimental/algorithm/LAGraph_MinCut.c @@ -0,0 +1,69 @@ +#include +#include "LG_internal.h" +#include + +#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; + GrB_Index n = 0; + GrB_Matrix_nrows(&n, R); + + LG_TRY(LAGraph_CheckGraph(G_origin, msg)); + LG_ASSERT (S != NULL, GrB_NULL_POINTER) ; + LG_ASSERT (S_bar != NULL, GrB_NULL_POINTER) ; + LG_ASSERT (cut_set != NULL, GrB_NULL_POINTER) ; + LG_ASSERT (s >= 0 && t >= 0, GrB_INVALID_VALUE) ; + + // printf("hit\n"); + GrB_Matrix A = G_origin->A; + + LG_TRY(GrB_Matrix_new(cut_set, GrB_FP64, n, n)); + LG_TRY(GrB_Vector_new(S_bar, GrB_INT64, n)); + //S is allocated during the bfs + + LG_TRY(LAGraph_New(&G, &R, LAGraph_ADJACENCY_DIRECTED, msg)); + LG_TRY(LAGr_BreadthFirstSearch(S, NULL, G, s, msg)); + + + 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)); + + 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, A, S_bar_diag, NULL)); + GRB_TRY(GrB_mxm(*cut_set, NULL, NULL, GxB_PLUS_TIMES_INT64, S_diag, *cut_set, NULL)); + + + LG_FREE_ALL; + return (GrB_SUCCESS); +} diff --git a/experimental/benchmark/maxflow_demo.c b/experimental/benchmark/maxflow_demo.c index 75fb56905c..53e6d04106 100644 --- a/experimental/benchmark/maxflow_demo.c +++ b/experimental/benchmark/maxflow_demo.c @@ -58,7 +58,7 @@ 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); @@ -66,7 +66,7 @@ int main (int argc, char ** argv){ 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); diff --git a/experimental/test/test_MaxFlow.c b/experimental/test/test_MaxFlow.c index 902e605ad4..e8c56c2384 100644 --- a/experimental/test/test_MaxFlow.c +++ b/experimental/test/test_MaxFlow.c @@ -30,7 +30,7 @@ GrB_Matrix A = NULL; #ifdef GRAPHBLAS_HAS_CUDA #define NTESTS 4 #else -#define NTESTS 7 +#define NTESTS 8 #endif char filename[LEN + 1]; @@ -47,6 +47,7 @@ test_info tests[] = { {"matrix_random_flow.mtx", 0,9, 22, LAGraph_ADJACENCY_DIRECTED}, {"rand.mtx", 0, 19, 37, LAGraph_ADJACENCY_DIRECTED}, {"mcl.mtx", 0, 9, 0, LAGraph_ADJACENCY_DIRECTED}, + {"test_zero_cap.mtx", 0, 4, 0.5, LAGraph_ADJACENCY_DIRECTED}, #ifndef GRAPHBLAS_HAS_CUDA // FIXME: the CUDA cases are currently very slow for these matrices, // when the GPU is hacked to always be used regardless of problem size: @@ -84,14 +85,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)); @@ -107,7 +108,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)); } } } @@ -152,7 +153,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); @@ -165,7 +166,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) ; diff --git a/experimental/test/test_MinCut.c b/experimental/test/test_MinCut.c new file mode 100644 index 0000000000..1ea1fad389 --- /dev/null +++ b/experimental/test/test_MinCut.c @@ -0,0 +1,93 @@ +#include +#include +#include +#include +#include + + +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); +#if LG_SUITESPARSE_GRAPHBLAS_V10 + 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(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); + + OK(GrB_reduce(&min_cut, NULL, GrB_PLUS_MONOID_FP64, 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); + } + +#endif + + + LAGraph_Finalize(msg); + +} + +TEST_LIST = {{"MinCut", test_MinCut}, {NULL, NULL}}; diff --git a/include/LAGraphX.h b/include/LAGraphX.h index 3355addb2e..6182f54805 100644 --- a/include/LAGraphX.h +++ b/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 ) }