From a6e4dbe1a0658f7eba0a0a39d4f5b02c5b19482e Mon Sep 17 00:00:00 2001 From: Sterkachov Date: Wed, 3 Dec 2025 09:28:14 +0300 Subject: [PATCH 1/3] feat: add RSM_reachability algorithm --- .../algorithm/LAGraph_RSM_reachability.c | 326 ++++++++++++++++++ include/LAGraphX.h | 31 +- 2 files changed, 356 insertions(+), 1 deletion(-) create mode 100644 experimental/algorithm/LAGraph_RSM_reachability.c diff --git a/experimental/algorithm/LAGraph_RSM_reachability.c b/experimental/algorithm/LAGraph_RSM_reachability.c new file mode 100644 index 0000000000..afd3482a3c --- /dev/null +++ b/experimental/algorithm/LAGraph_RSM_reachability.c @@ -0,0 +1,326 @@ +//------------------------------------------------------------------------------ +// LAGraph_RSM_reachability.c +//------------------------------------------------------------------------------ + +// For and edge-labelled directed graph the algorithm computes the set of nodes +// for which these conditions are held: +// * The node is reachable by a path from one of the source nodes. +// * The concatenation of the labels over this path's edges is a word from the +// specified context-free language +// +// The context-free constraints are specified by a recursive state machine (RSM) +// over a subset of the graph edge labels. The algorithm assumes the labels are +// enumerated from 0 to some nl-1. The input graph and the RSM are defined by +// adjacency matrix decomposition. They are represented by arrays of graphs G and +// R both of length + +#include "GraphBLAS.h" +#define LG_FREE_WORK \ +{ \ + GrB_free (&frontier) ; \ + GrB_free (&next_frontier) ; \ + GrB_free (&symbol_frontier) ; \ + GrB_free (&visited) ; \ + GrB_free (&final_reducer) ; \ + LAGraph_Free ((void **) &A, GrB_NULL) ; \ + LAGraph_Free ((void **) &B, GrB_NULL) ; \ + LAGraph_Free ((void **) &BT, GrB_NULL) ; \ + LAGraph_Free ((void **) &B_call, GrB_NULL) ; \ + LAGraph_Free ((void **) &BT_call, GrB_NULL) ; \ +} + +#define LG_FREE_ALL \ +{ \ + LG_FREE_WORK ; \ + GrB_free (reachable) ; \ +} + +#include "LG_internal.h" +#include "LAGraphX.h" + +int LAGraph_RSM_reachability +( + // output: + GrB_Vector *reachable, // reachable(i) = true if node i is reachable + // from one of the starting nodes by a path + // satisfying regular constraints + + LAGraph_Graph *R, // input recursive state machine's + // adjacency matrix decomposition + size_t nl, // total label count, # of matrices graph and + // RSM adjacency matrix decomposition + const GrB_Index *QS, // starting states in RSM + size_t nqs, // number of starting states in RSM + const GrB_Index *QF, // final states in RSM + size_t nqf, // number of final states in RSM + + LAGraph_Graph *R_call, // recursive state machine's call/return nr*nr + // matrix, where (i, j) = 1 if there is transition + // from i state of one automata to j starting + // state of another automata + LAGraph_Graph *G, // input graph adjacency matrix decomposition + const GrB_Index *S, // source vertices to start searching paths + size_t ns, // number of source vertices + char *msg // LAGraph output message +) +{ + //-------------------------------------------------------------------------- + // check inputs + //-------------------------------------------------------------------------- + + LG_CLEAR_MSG ; + + GrB_Matrix frontier = GrB_NULL ; // traversal frontier representing + // correspondence between RSM states + // and graph vertices + GrB_Matrix symbol_frontier = GrB_NULL ; // part of the new frontier for the + // specific label + GrB_Matrix next_frontier = GrB_NULL ; // frontier value on the next + // traversal step + GrB_Matrix visited = GrB_NULL ; // visited pairs (state, vertex) + + GrB_Vector final_reducer = GrB_NULL ; // auxiliary vector for reducing the + // visited matrix to an answer + + GrB_Index ng = 0 ; // # nodes in the graph + GrB_Index nr = 0 ; // # states in the RSM + GrB_Index states = ns ; // # pairs in the current + // correspondence between the graph and + // the RSM + + GrB_Index rows = 0 ; // utility matrix row count + GrB_Index cols = 0 ; // utility matrix column count + GrB_Index vals = 0 ; // utility matrix valiue count + + GrB_Matrix *A = GrB_NULL ; + GrB_Matrix *B = GrB_NULL ; + GrB_Matrix *BT = GrB_NULL ; + + GrB_Matrix *B_call = GrB_NULL ; + GrB_Matrix *BT_call = GrB_NULL ; + + LG_ASSERT (reachable != GrB_NULL, GrB_NULL_POINTER) ; + LG_ASSERT (R != GrB_NULL, GrB_NULL_POINTER) ; + LG_ASSERT (G != GrB_NULL, GrB_NULL_POINTER) ; + LG_ASSERT (S != GrB_NULL, GrB_NULL_POINTER) ; + + LG_ASSERT (R_call != GrB_NULL, GrB_NULL_POINTER) ; + LG_ASSERT (*R_call != GrB_NULL, GrB_NULL_POINTER) ; + + (*reachable) = GrB_NULL ; + + for (size_t i = 0 ; i < nl ; ++i) + { + if (G[i] == GrB_NULL) continue ; + LG_TRY (LAGraph_CheckGraph (G[i], msg)) ; + } + + for (size_t i = 0 ; i < nl ; ++i) + { + if (R[i] == GrB_NULL) continue; + LG_TRY (LAGraph_CheckGraph (R[i], msg)) ; + } + + LG_TRY (LAGraph_Malloc ((void **) &A, nl, sizeof (GrB_Matrix), msg)) ; + + for (size_t i = 0 ; i < nl ; ++i) + { + if (G[i] == GrB_NULL) + { + A[i] = GrB_NULL ; + continue ; + } + + A[i] = G[i]->A ; + } + + LG_TRY (LAGraph_Malloc((void **) &B, nl, sizeof (GrB_Matrix), msg)) ; + LG_TRY (LAGraph_Malloc((void **) &BT, nl, sizeof (GrB_Matrix), msg)) ; + + for (size_t i = 0 ; i < nl ; ++i) + { + BT[i] = GrB_NULL ; + + if (R[i] == GrB_NULL) + { + B[i] = GrB_NULL ; + continue ; + } + + B[i] = R[i]->A ; + if (R[i]->is_symmetric_structure == LAGraph_TRUE) + { + BT[i] = B[i] ; + } else + { + // BT[i] could be NULL and the matrix will be transposed by a + // descriptor + BT[i] = R[i]->AT ; + } + } + + + + LG_TRY (LAGraph_Malloc((void **) &B_call, 1, sizeof (GrB_Matrix), msg)) ; + LG_TRY (LAGraph_Malloc((void **) &BT_call, 1, sizeof (GrB_Matrix), msg)) ; + + *B_call = (*R_call)->A ; + *BT_call = (*R_call)->AT ; + + + + for (size_t i = 0 ; i < nl ; ++i) + { + if (A[i] == GrB_NULL) continue ; + + GRB_TRY (GrB_Matrix_nrows (&ng, A[i])) ; + break ; + } + + for (size_t i = 0 ; i < nl ; ++i) + { + if (B[i] == GrB_NULL) continue ; + + GRB_TRY (GrB_Matrix_nrows (&nr, B[i])) ; + break ; + } + + // Check all the matrices in graph adjacency matrix decomposition are + // square and of the same dimensions + for (size_t i = 0 ; i < nl ; ++i) + { + if (A[i] == GrB_NULL) continue ; + + GRB_TRY (GrB_Matrix_nrows (&rows, A[i])) ; + GRB_TRY (GrB_Matrix_ncols (&cols, A[i])) ; + + LG_ASSERT_MSG (rows == ng && cols == ng, LAGRAPH_NOT_CACHED, + "all the matrices in the graph adjacency matrix decomposition " + "should have the same dimensions and be square") ; + } + + // Check all the matrices in RSM adjancency matrix decomposition are + // square and of the same dimensions + for (size_t i = 0 ; i < nl ; ++i) + { + if (B[i] == GrB_NULL) continue ; + + GRB_TRY (GrB_Matrix_nrows (&rows, B[i])) ; + GRB_TRY (GrB_Matrix_ncols (&cols, B[i])) ; + + LG_ASSERT_MSG (rows == nr && cols == nr, LAGRAPH_NOT_CACHED, + "all the matrices in the RSM adjacency matrix decomposition " + "should have the same dimensions and be square") ; + } + + // Check size of RSM call/return matrix + + GRB_TRY (GrB_Matrix_nrows (&rows, *B_call)) ; + GRB_TRY (GrB_Matrix_ncols (&cols, *B_call)) ; + + LG_ASSERT_MSG (rows == nr && cols == nr, LAGRAPH_NOT_CACHED, + "RSM's call/return matrix should have nr*nr size") ; + + // Check source nodes in the graph + for (size_t i = 0 ; i < ns ; ++i) + { + GrB_Index s = S[i] ; + LG_ASSERT_MSG(s < ng, GrB_INVALID_INDEX, "invalid graph source node") ; + } + + // Check starting states of the RSM + for (size_t i = 0 ; i < nqs ; ++i) + { + GrB_Index qs = QS[i] ; + LG_ASSERT_MSG (qs < nr, GrB_INVALID_INDEX, + "invalid RSM starting state") ; + } + + // Check final states of the RSM + for (size_t i = 0 ; i < nqf ; ++i) + { + GrB_Index qf = QF[i] ; + LG_ASSERT_MSG (qf < nr, GrB_INVALID_INDEX, "invalid RSM final state") ; + } + + // ------------------------------------------------------------------------- + // initialization + // ------------------------------------------------------------------------- + + GRB_TRY (GrB_Vector_new (reachable, GrB_BOOL, ng)) ; + GRB_TRY (GrB_Vector_new (&final_reducer, GrB_BOOL, nr)) ; + + // Initialize matrix for reducing the result + GrB_assign (final_reducer, GrB_NULL, GrB_NULL, true, QF, nqf, GrB_NULL) ; + + GRB_TRY (GrB_Matrix_new (&frontier, GrB_BOOL, nr, ng)) ; + GRB_TRY (GrB_Matrix_new (&next_frontier, GrB_BOOL, nr, ng)) ; + GRB_TRY (GrB_Matrix_new (&symbol_frontier, GrB_BOOL, nr, ng)) ; + GRB_TRY (GrB_Matrix_new (&visited, GrB_BOOL, nr, ng)) ; + + // Initialize frontier with the source nodes + GrB_assign (next_frontier, GrB_NULL, GrB_NULL, true, QS, nqs, S, ns, GrB_NULL) ; + GrB_assign (visited, GrB_NULL, GrB_NULL, true, QS, nqs, S, ns, GrB_NULL) ; + + // Main loop + while (states != 0) + { + GrB_Matrix old_frontier = frontier ; + frontier = next_frontier ; + next_frontier = old_frontier ; + + GRB_TRY (GrB_Matrix_clear (next_frontier)) ; + + // Obtain a new relation between the RSM states and the graph nodes + for (size_t i = 0 ; i < nl ; ++i) + { + if (A[i] == GrB_NULL || B[i] == GrB_NULL) continue ; + + // Traverse the RSM + // Try to use a provided transposed matrix or use the descriptor + + GRB_TRY (GrB_Matrix_clear (symbol_frontier)) ; + + if (BT[i] != GrB_NULL) + { + GRB_TRY (GrB_mxm (symbol_frontier, GrB_NULL, GrB_LOR, + GrB_LOR_LAND_SEMIRING_BOOL, BT[i], frontier, GrB_NULL)) ; + } + else + { + GRB_TRY (GrB_mxm (symbol_frontier, GrB_NULL, GrB_LOR, + GrB_LOR_LAND_SEMIRING_BOOL, B[i], frontier, GrB_DESC_T0)) ; + } + + if ((*BT_call) != GrB_NULL) + { + GRB_TRY (GrB_mxm (next_frontier, visited, GrB_LOR, + GrB_LOR_LAND_SEMIRING_BOOL, *BT_call, frontier, GrB_DESC_SC)) ; // не уверен по поводу дескриптора + } + else + { + GRB_TRY (GrB_mxm (next_frontier, visited, GrB_LOR, + GrB_LOR_LAND_SEMIRING_BOOL, *B_call, frontier, GrB_DESC_SCT0)) ; // не уверен по поводу дескриптора + } + + // я реально не понимаю как это работает, но окей + // Traverse the graph + GRB_TRY (GrB_mxm (next_frontier, visited, GrB_LOR, + GrB_LOR_LAND_SEMIRING_BOOL, symbol_frontier, A[i], GrB_DESC_SC)) ; + } + + // я реально не понимаю как это работает, но окей + // Accumulate the new state <-> node correspondence + GRB_TRY (GrB_assign (visited, visited, GrB_NULL, next_frontier, + GrB_ALL, nr, GrB_ALL, ng, GrB_DESC_SC)) ; + + GRB_TRY (GrB_Matrix_nvals (&states, next_frontier)) ; + } + + // Extract the nodes matching the final RSM states + GRB_TRY (GrB_vxm (*reachable, GrB_NULL, GrB_NULL, + GrB_LOR_LAND_SEMIRING_BOOL, final_reducer, visited, GrB_NULL)) ; + + LG_FREE_WORK ; + return (GrB_SUCCESS) ; +} diff --git a/include/LAGraphX.h b/include/LAGraphX.h index 3355addb2e..d7053d0244 100644 --- a/include/LAGraphX.h +++ b/include/LAGraphX.h @@ -854,7 +854,36 @@ int LAGraph_RegularPathQuery // nodes reachable from the starting by the const GrB_Index *S, // source vertices to start searching paths size_t ns, // number of source vertices char *msg // LAGraph output message -); +) ; + +//**************************************************************************** +LAGRAPHX_PUBLIC +int LAGraph_RSM_reachability +( + // output: + GrB_Vector *reachable, // reachable(i) = true if node i is reachable + // from one of the starting nodes by a path + // satisfying regular constraints + + LAGraph_Graph *R, // input recursive state machine's + // adjacency matrix decomposition + size_t nl, // total label count, # of matrices graph and + // RSM adjacency matrix decomposition + const GrB_Index *QS, // starting states in RSM + size_t nqs, // number of starting states in RSM + const GrB_Index *QF, // final states in RSM + size_t nqf, // number of final states in RSM + + LAGraph_Graph *R_call, // recursive state machine's call/return nr*nr + // matrix, where (i, j) = 1 if there is transition + // from i state of one automata to j starting + // state of another automata + LAGraph_Graph *G, // input graph adjacency matrix decomposition + const GrB_Index *S, // source vertices to start searching paths + size_t ns, // number of source vertices + char *msg // LAGraph output message +) ; + //**************************************************************************** LAGRAPHX_PUBLIC int LAGraph_VertexCentrality_Triangle // vertex triangle-centrality From 2dd9590f7ad0020feb4b9175bcb4d1df3dd2e237 Mon Sep 17 00:00:00 2001 From: Sterkachov Date: Wed, 3 Dec 2025 09:29:54 +0300 Subject: [PATCH 2/3] feat: add simple test for RSM_reachability algo --- data/rsm_data/1_a.mtx | 4 + data/rsm_data/1_call.mtx | 5 + data/rsm_data/1_meta.txt | 2 + data/rsm_data/1_sources.txt | 1 + data/rsm_data/a.mtx | 4 + data/rsm_data/b.mtx | 4 + experimental/test/test_RSM_reachability.c | 243 ++++++++++++++++++++++ 7 files changed, 263 insertions(+) create mode 100644 data/rsm_data/1_a.mtx create mode 100644 data/rsm_data/1_call.mtx create mode 100644 data/rsm_data/1_meta.txt create mode 100644 data/rsm_data/1_sources.txt create mode 100644 data/rsm_data/a.mtx create mode 100644 data/rsm_data/b.mtx create mode 100644 experimental/test/test_RSM_reachability.c diff --git a/data/rsm_data/1_a.mtx b/data/rsm_data/1_a.mtx new file mode 100644 index 0000000000..0430115a88 --- /dev/null +++ b/data/rsm_data/1_a.mtx @@ -0,0 +1,4 @@ +%%MatrixMarket matrix coordinate pattern general$ +%%GraphBLAS type bool$ +4 4 1 +3 4 diff --git a/data/rsm_data/1_call.mtx b/data/rsm_data/1_call.mtx new file mode 100644 index 0000000000..1340e468bd --- /dev/null +++ b/data/rsm_data/1_call.mtx @@ -0,0 +1,5 @@ +%%MatrixMarket matrix coordinate pattern general +%%GraphBLAS type bool +4 4 2 +1 3 +4 2 diff --git a/data/rsm_data/1_meta.txt b/data/rsm_data/1_meta.txt new file mode 100644 index 0000000000..1f25976c26 --- /dev/null +++ b/data/rsm_data/1_meta.txt @@ -0,0 +1,2 @@ +1 1 +1 2 diff --git a/data/rsm_data/1_sources.txt b/data/rsm_data/1_sources.txt new file mode 100644 index 0000000000..56a6051ca2 --- /dev/null +++ b/data/rsm_data/1_sources.txt @@ -0,0 +1 @@ +1 \ No newline at end of file diff --git a/data/rsm_data/a.mtx b/data/rsm_data/a.mtx new file mode 100644 index 0000000000..3095351475 --- /dev/null +++ b/data/rsm_data/a.mtx @@ -0,0 +1,4 @@ +%%MatrixMarket matrix coordinate pattern general$ +%%GraphBLAS type bool$ +3 3 1 +1 2 diff --git a/data/rsm_data/b.mtx b/data/rsm_data/b.mtx new file mode 100644 index 0000000000..170cfb40dd --- /dev/null +++ b/data/rsm_data/b.mtx @@ -0,0 +1,4 @@ +%%MatrixMarket matrix coordinate pattern general$ +%%GraphBLAS type bool$ +3 3 1 +1 3 diff --git a/experimental/test/test_RSM_reachability.c b/experimental/test/test_RSM_reachability.c new file mode 100644 index 0000000000..121dc511f7 --- /dev/null +++ b/experimental/test/test_RSM_reachability.c @@ -0,0 +1,243 @@ +#include "GraphBLAS.h" +#include "LAGraph.h" +#include +#include +#include +#include +#include +#include + +#define LEN 512 +#define MAX_LABELS 3 +#define MAX_RESULTS 2000000 + +LAGraph_Graph G[MAX_LABELS] ; +LAGraph_Graph R[MAX_LABELS] ; +LAGraph_Graph R_call; + +GrB_Matrix A ; + +char testcase_name[LEN+1] ; +char filename[LEN+1] ; +char msg[LAGRAPH_MSG_LEN] ; + +typedef struct +{ + const char* name ; + const char* graphs[MAX_LABELS] ; + const char* rsm[MAX_LABELS] ; + const char* rsm_call; + const char* rsm_meta ; + const char* sources ; + const GrB_Index expected[MAX_RESULTS] ; + const size_t expected_count ; +} +matrix_info ; + +const matrix_info files [ ] = +{ + {"simple 1", + {"rsm_data/a.mtx", "rsm_data/b.mtx", NULL}, + {"rsm_data/1_a.mtx", NULL}, + "rsm_data/1_call.mtx", + "rsm_data/1_meta.txt", + "rsm_data/1_sources.txt", + {2}, + 1}, + {NULL, NULL, NULL, NULL}, +} ; + + +void test_RSM_reachability (void) +{ + LAGraph_Init (msg) ; + + for (int k = 0 ; ; ++k) + { + if (files[k].sources == NULL) break ; + + snprintf (testcase_name, LEN, "basic context-free path query %s", files[k].name) ; + TEST_CASE (testcase_name) ; + + for (int check_symmetry = 0 ; check_symmetry < 2 ; ++check_symmetry) + { + // Load graph from MTX files represention its adjacency matrix + // decomposition + for (int i = 0 ; ; ++i) + { + const char *name = files[k].graphs[i] ; + + if (name == NULL) break ; + if (strlen(name) == 0) continue ; + + snprintf (filename, LEN, LG_DATA_DIR "%s", name) ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + + OK (LAGraph_New (&(G[i]), &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + + TEST_CHECK (A == NULL) ; + } + + // Load RSM from MTX files representing its adjacency matrix + // decomposition + for (int i = 0 ; ; ++i) + { + const char *name = files[k].rsm[i] ; + + if (name == NULL) break ; + if (strlen(name) == 0) continue ; + + snprintf (filename, LEN, LG_DATA_DIR "%s", name) ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + + OK (LAGraph_New (&(R[i]), &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + + if (check_symmetry) + { + // Check if the pattern is symmetric - if it isn't make it. + // Note this also computes R[i]->AT + OK (LAGraph_Cached_IsSymmetricStructure (R[i], msg)) ; + } + + TEST_CHECK (A == NULL) ; + } + + { + const char *name = files[k].rsm_call ; + + snprintf (filename, LEN, LG_DATA_DIR "%s", name) ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + + OK (LAGraph_New (&(R_call), &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + + if (check_symmetry) + { + OK (LAGraph_Cached_IsSymmetricStructure(R_call, msg)) ; + } + + TEST_CHECK (A == NULL) ; // зачем нам это проверять?, почему оно NULL + } + + } + + // Note the matrix rows/cols are enumerated from 0 to n-1. + // Meanwhile, in MTX format they are enumerated from 1 to n. Thus, + // when loading/comparing the results these values should be + // decremented/incremented correspondingly. + + // Load graph source nodes from the sources file + GrB_Index s ; + GrB_Index S[16] ; + size_t ns = 0 ; + + const char *name = files[k].sources ; + snprintf (filename, LEN, LG_DATA_DIR "%s", name) ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + + while (fscanf(f, "%" PRIu64, &s) != EOF) + { + S[ns++] = s - 1 ; + } + + OK (fclose(f)) ; + + // Load RSM starting states from the meta file + GrB_Index qs ; + GrB_Index QS[16] ; // почему именно 16, мб личный выбор + size_t nqs = 0 ; + + name = files[k].rsm_meta ; + snprintf (filename, LEN, LG_DATA_DIR "%s", name) ; + f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; ; + + uint64_t nqs64 = 0 ; + TEST_CHECK (fscanf (f, "%" PRIu64, &nqs64) != EOF) ; + nqs = (size_t) nqs64 ; + + for (size_t i = 0 ; i < nqs ; ++i) + { + TEST_CHECK(fscanf (f, "%" PRIu64, &qs) != EOF) ; + QS[i] = qs - 1 ; + } + + // Load RSM final states from the same file + uint64_t qf ; + uint64_t QF[16] ; + size_t nqf = 0 ; + uint64_t nqf64 = 0 ; + + TEST_CHECK (fscanf (f, "%" PRIu64, &nqf64) != EOF) ; + nqf = (size_t) nqf64 ; + + for (size_t i = 0 ; i < nqf ; ++i) + { + TEST_CHECK (fscanf (f, "%" PRIu64, &qf) != EOF) ; + QF[i] = qf - 1 ; + } + + OK (fclose (f)) ; + + // Evaluate the algorithm + GrB_Vector r = NULL ; + + OK (LAGraph_RSM_reachability (&r, R, MAX_LABELS, QS, nqs, + QF, nqf, &R_call, G, S, ns, msg)) ; + + // Extract results from the output vector + GrB_Index *reachable ; + bool *values ; + + GrB_Index nvals ; + GrB_Vector_nvals (&nvals, r) ; + + OK (LAGraph_Malloc ((void **) &reachable, MAX_RESULTS, sizeof (GrB_Index), msg)) ;; + OK (LAGraph_Malloc ((void **) &values, MAX_RESULTS, sizeof (bool), msg)) ; + + GrB_Vector_extractTuples (reachable, values, &nvals, r) ; + + TEST_MSG("returned %lu values:\n", nvals); + for (uint64_t i = 0; i < nvals; i++) + TEST_MSG(" %lu\n", reachable[i] + 1); + + // Compare the results with expected values + TEST_CHECK (nvals == files[k].expected_count) ; + for (uint64_t i = 0 ; i < nvals ; ++i) + TEST_CHECK (reachable[i] + 1 == files[k].expected[i]) ; + + // Cleanup + OK (LAGraph_Free ((void **) &values, NULL)) ; + OK (LAGraph_Free ((void **) &reachable, NULL)) ; + + OK (GrB_free (&r)) ; + + for (uint64_t i = 0 ; i < MAX_LABELS ; ++i) + { + if (G[i] == NULL) continue ; + OK (LAGraph_Delete (&(G[i]), msg)) ; + } + + for (uint64_t i = 0; i < MAX_LABELS ; ++i) + { + if (R[i] == NULL) continue ; + OK (LAGraph_Delete (&(R[i]), msg)) ; + } + } + + LAGraph_Finalize (msg) ; +} + +TEST_LIST = { + {"RSM_reachability", test_RSM_reachability}, + {NULL, NULL} +} ; From f079f768dab1841f37e1d7540365d028dc0c11e6 Mon Sep 17 00:00:00 2001 From: Sterkachov Date: Fri, 20 Feb 2026 22:37:03 +0300 Subject: [PATCH 3/3] feat: implement the algorithm of rsm reachability --- data/rsm_data/1_a.mtx | 4 - data/rsm_data/1_call.mtx | 5 - data/rsm_data/1_meta.txt | 2 - data/rsm_data/1_sources.txt | 1 - data/rsm_data/a.mtx | 4 - data/rsm_data/b.mtx | 4 - .../algorithm/LAGraph_RSM_reachability.c | 781 ++++++++++++------ experimental/test/test_RSM_reachability.c | 488 ++++++----- include/LAGraphX.h | 48 +- 9 files changed, 859 insertions(+), 478 deletions(-) delete mode 100644 data/rsm_data/1_a.mtx delete mode 100644 data/rsm_data/1_call.mtx delete mode 100644 data/rsm_data/1_meta.txt delete mode 100644 data/rsm_data/1_sources.txt delete mode 100644 data/rsm_data/a.mtx delete mode 100644 data/rsm_data/b.mtx diff --git a/data/rsm_data/1_a.mtx b/data/rsm_data/1_a.mtx deleted file mode 100644 index 0430115a88..0000000000 --- a/data/rsm_data/1_a.mtx +++ /dev/null @@ -1,4 +0,0 @@ -%%MatrixMarket matrix coordinate pattern general$ -%%GraphBLAS type bool$ -4 4 1 -3 4 diff --git a/data/rsm_data/1_call.mtx b/data/rsm_data/1_call.mtx deleted file mode 100644 index 1340e468bd..0000000000 --- a/data/rsm_data/1_call.mtx +++ /dev/null @@ -1,5 +0,0 @@ -%%MatrixMarket matrix coordinate pattern general -%%GraphBLAS type bool -4 4 2 -1 3 -4 2 diff --git a/data/rsm_data/1_meta.txt b/data/rsm_data/1_meta.txt deleted file mode 100644 index 1f25976c26..0000000000 --- a/data/rsm_data/1_meta.txt +++ /dev/null @@ -1,2 +0,0 @@ -1 1 -1 2 diff --git a/data/rsm_data/1_sources.txt b/data/rsm_data/1_sources.txt deleted file mode 100644 index 56a6051ca2..0000000000 --- a/data/rsm_data/1_sources.txt +++ /dev/null @@ -1 +0,0 @@ -1 \ No newline at end of file diff --git a/data/rsm_data/a.mtx b/data/rsm_data/a.mtx deleted file mode 100644 index 3095351475..0000000000 --- a/data/rsm_data/a.mtx +++ /dev/null @@ -1,4 +0,0 @@ -%%MatrixMarket matrix coordinate pattern general$ -%%GraphBLAS type bool$ -3 3 1 -1 2 diff --git a/data/rsm_data/b.mtx b/data/rsm_data/b.mtx deleted file mode 100644 index 170cfb40dd..0000000000 --- a/data/rsm_data/b.mtx +++ /dev/null @@ -1,4 +0,0 @@ -%%MatrixMarket matrix coordinate pattern general$ -%%GraphBLAS type bool$ -3 3 1 -1 3 diff --git a/experimental/algorithm/LAGraph_RSM_reachability.c b/experimental/algorithm/LAGraph_RSM_reachability.c index afd3482a3c..8972116a4c 100644 --- a/experimental/algorithm/LAGraph_RSM_reachability.c +++ b/experimental/algorithm/LAGraph_RSM_reachability.c @@ -2,325 +2,636 @@ // LAGraph_RSM_reachability.c //------------------------------------------------------------------------------ -// For and edge-labelled directed graph the algorithm computes the set of nodes -// for which these conditions are held: -// * The node is reachable by a path from one of the source nodes. -// * The concatenation of the labels over this path's edges is a word from the -// specified context-free language -// -// The context-free constraints are specified by a recursive state machine (RSM) -// over a subset of the graph edge labels. The algorithm assumes the labels are -// enumerated from 0 to some nl-1. The input graph and the RSM are defined by -// adjacency matrix decomposition. They are represented by arrays of graphs G and -// R both of length - -#include "GraphBLAS.h" -#define LG_FREE_WORK \ -{ \ - GrB_free (&frontier) ; \ - GrB_free (&next_frontier) ; \ - GrB_free (&symbol_frontier) ; \ - GrB_free (&visited) ; \ - GrB_free (&final_reducer) ; \ - LAGraph_Free ((void **) &A, GrB_NULL) ; \ - LAGraph_Free ((void **) &B, GrB_NULL) ; \ - LAGraph_Free ((void **) &BT, GrB_NULL) ; \ - LAGraph_Free ((void **) &B_call, GrB_NULL) ; \ - LAGraph_Free ((void **) &BT_call, GrB_NULL) ; \ -} - -#define LG_FREE_ALL \ -{ \ - LG_FREE_WORK ; \ - GrB_free (reachable) ; \ -} +#include +#include +#include +#include +#include +#include "LAGraph.h" #include "LG_internal.h" -#include "LAGraphX.h" -int LAGraph_RSM_reachability -( +#define LG_FREE_WORK \ + { \ + GrB_free(&frontier); \ + GrB_free(&next_frontier); \ + GrB_free(&front_temp1); \ + GrB_free(&front_temp2); \ + GrB_free(&P); \ + GrB_free(&K); \ + GrB_free(&final_reducer); \ + GrB_free(&M_ret); \ + GrB_free(&M_call); \ + GrB_free(&M_call_new); \ + GrB_free(&G_S_new); \ + GrB_free(&step1); \ + LAGraph_Free((void **)&A_a, GrB_NULL); \ + LAGraph_Free((void **)&B_a, GrB_NULL); \ + LAGraph_Free((void **)&B_S, GrB_NULL); \ + LAGraph_Free((void **)&B_call, GrB_NULL); \ + LAGraph_Free((void **)&B_ret, GrB_NULL); \ + LAGraph_Free((void **)&B_S, GrB_NULL); \ + if (A_S != GrB_NULL) \ + { \ + for (size_t _i = 0; _i < num_nonterminals; ++_i) \ + GrB_free(&A_S[_i]); \ + LAGraph_Free((void **)&A_S, GrB_NULL); \ + } \ + } + +#define LG_FREE_ALL \ + { \ + LG_FREE_WORK; \ + GrB_free(reachable); \ + } + +int LAGraph_RSM_reachability( // output: - GrB_Vector *reachable, // reachable(i) = true if node i is reachable - // from one of the starting nodes by a path - // satisfying regular constraints - - LAGraph_Graph *R, // input recursive state machine's - // adjacency matrix decomposition - size_t nl, // total label count, # of matrices graph and - // RSM adjacency matrix decomposition - const GrB_Index *QS, // starting states in RSM - size_t nqs, // number of starting states in RSM - const GrB_Index *QF, // final states in RSM - size_t nqf, // number of final states in RSM - - LAGraph_Graph *R_call, // recursive state machine's call/return nr*nr - // matrix, where (i, j) = 1 if there is transition - // from i state of one automata to j starting - // state of another automata - LAGraph_Graph *G, // input graph adjacency matrix decomposition - const GrB_Index *S, // source vertices to start searching paths - size_t ns, // number of source vertices - char *msg // LAGraph output message -) + GrB_Vector *reachable, + + // input: + size_t num_terminals, // # terminal labels + LAGraph_Graph *G_a, // terminal transitions + LAGraph_Graph *N_a, // terminal transitions + + size_t num_nonterminals, // # nonterminal labels + LAGraph_Graph *N_S, // nonterminal RSM transitions + LAGraph_Graph *Call_S, // call transition matrix + LAGraph_Graph *Ret_S, // return transition matrix + + const GrB_Index *QS, // staring states in RSM + size_t nqs, // # starting states + const GrB_Index *QF, // final states in RSM + size_t nqf, // # final states + + const GrB_Index *Source, // sources vertices + size_t ns, // # source vertices + + char *msg + +) { //-------------------------------------------------------------------------- // check inputs //-------------------------------------------------------------------------- - LG_CLEAR_MSG ; + LG_CLEAR_MSG; + + // Working matrices + GrB_Matrix frontier = GrB_NULL; // traversal frontier + GrB_Matrix next_frontier = GrB_NULL; + GrB_Matrix front_temp1 = GrB_NULL; + GrB_Matrix front_temp2 = GrB_NULL; + GrB_Matrix P = GrB_NULL; // visited pairs (state, vertex) + GrB_Matrix K = GrB_NULL; // visited pairs in second loop + GrB_Vector final_reducer = GrB_NULL; // vector for reducing the visited matrix + + // Backward-phase working matrices + GrB_Matrix M_ret = GrB_NULL; + GrB_Matrix M_call = GrB_NULL; + GrB_Matrix M_call_new = GrB_NULL; + GrB_Matrix G_S_new = GrB_NULL; + GrB_Matrix step1 = GrB_NULL; - GrB_Matrix frontier = GrB_NULL ; // traversal frontier representing - // correspondence between RSM states - // and graph vertices - GrB_Matrix symbol_frontier = GrB_NULL ; // part of the new frontier for the - // specific label - GrB_Matrix next_frontier = GrB_NULL ; // frontier value on the next - // traversal step - GrB_Matrix visited = GrB_NULL ; // visited pairs (state, vertex) + GrB_Index ng = 0; // # nodes in the graph + GrB_Index nr = 0; // # states in the RSM + GrB_Index states = ns; // # pairs in current frontier - GrB_Vector final_reducer = GrB_NULL ; // auxiliary vector for reducing the - // visited matrix to an answer + GrB_Index rows = 0, cols = 0, vals = 0; - GrB_Index ng = 0 ; // # nodes in the graph - GrB_Index nr = 0 ; // # states in the RSM - GrB_Index states = ns ; // # pairs in the current - // correspondence between the graph and - // the RSM + GrB_Matrix *A_a = GrB_NULL; // for G_a + GrB_Matrix *B_a = GrB_NULL; // for N_a - GrB_Index rows = 0 ; // utility matrix row count - GrB_Index cols = 0 ; // utility matrix column count - GrB_Index vals = 0 ; // utility matrix valiue count + GrB_Matrix *A_S = GrB_NULL; // Derived graph edge matrices + GrB_Matrix *B_S = GrB_NULL; // for N_S + GrB_Matrix *B_call = GrB_NULL; // for Call_S + GrB_Matrix *B_ret = GrB_NULL; // for Ret_S - GrB_Matrix *A = GrB_NULL ; - GrB_Matrix *B = GrB_NULL ; - GrB_Matrix *BT = GrB_NULL ; + LG_ASSERT(reachable != GrB_NULL, GrB_NULL_POINTER); - GrB_Matrix *B_call = GrB_NULL ; - GrB_Matrix *BT_call = GrB_NULL ; + LG_ASSERT(G_a != GrB_NULL, GrB_NULL_POINTER); + LG_ASSERT(N_a != GrB_NULL, GrB_NULL_POINTER); + LG_ASSERT(N_S != GrB_NULL, GrB_NULL_POINTER); + LG_ASSERT(Call_S != GrB_NULL, GrB_NULL_POINTER); + LG_ASSERT(Ret_S != GrB_NULL, GrB_NULL_POINTER); + LG_ASSERT(QS != GrB_NULL, GrB_NULL_POINTER); + LG_ASSERT(QF != GrB_NULL, GrB_NULL_POINTER); + LG_ASSERT(Source != GrB_NULL, GrB_NULL_POINTER); - LG_ASSERT (reachable != GrB_NULL, GrB_NULL_POINTER) ; - LG_ASSERT (R != GrB_NULL, GrB_NULL_POINTER) ; - LG_ASSERT (G != GrB_NULL, GrB_NULL_POINTER) ; - LG_ASSERT (S != GrB_NULL, GrB_NULL_POINTER) ; + (*reachable) = GrB_NULL; - LG_ASSERT (R_call != GrB_NULL, GrB_NULL_POINTER) ; - LG_ASSERT (*R_call != GrB_NULL, GrB_NULL_POINTER) ; + //-------------------------------------------------------------------------- + // validate input graphs + //-------------------------------------------------------------------------- + + for (size_t i = 0; i < num_terminals; ++i) + if (G_a[i] != GrB_NULL) + LG_TRY(LAGraph_CheckGraph(G_a[i], msg)) + + for (size_t i = 0; i < num_terminals; ++i) + if (N_a[i] != GrB_NULL) + LG_TRY(LAGraph_CheckGraph(N_a[i], msg)) - (*reachable) = GrB_NULL ; + for (size_t i = 0; i < num_nonterminals; ++i) + if (N_S[i] != GrB_NULL) + LG_TRY(LAGraph_CheckGraph(N_S[i], msg)) - for (size_t i = 0 ; i < nl ; ++i) + for (size_t i = 0; i < num_nonterminals; ++i) + if (Call_S[i] != GrB_NULL) + LG_TRY(LAGraph_CheckGraph(Call_S[i], msg)) + + for (size_t i = 0; i < num_nonterminals; ++i) + if (Ret_S[i] != GrB_NULL) + LG_TRY(LAGraph_CheckGraph(Ret_S[i], msg)) + + LG_TRY(LAGraph_Malloc((void **)&A_a, num_terminals, sizeof(GrB_Matrix), msg)) + for (size_t i = 0; i < num_terminals; ++i) { - if (G[i] == GrB_NULL) continue ; - LG_TRY (LAGraph_CheckGraph (G[i], msg)) ; + if (G_a[i] == GrB_NULL) + A_a[i] = GrB_NULL; + else + A_a[i] = G_a[i]->A; } - - for (size_t i = 0 ; i < nl ; ++i) + LG_TRY(LAGraph_Malloc((void **)&B_a, num_terminals, sizeof(GrB_Matrix), msg)) + for (size_t i = 0; i < num_terminals; ++i) { - if (R[i] == GrB_NULL) continue; - LG_TRY (LAGraph_CheckGraph (R[i], msg)) ; + if (N_a[i] == GrB_NULL) + B_a[i] = GrB_NULL; + else + B_a[i] = N_a[i]->A; } - - LG_TRY (LAGraph_Malloc ((void **) &A, nl, sizeof (GrB_Matrix), msg)) ; - - for (size_t i = 0 ; i < nl ; ++i) + LG_TRY(LAGraph_Malloc((void **)&B_S, num_nonterminals, sizeof(GrB_Matrix), msg)) + for (size_t i = 0; i < num_nonterminals; ++i) { - if (G[i] == GrB_NULL) - { - A[i] = GrB_NULL ; - continue ; - } - - A[i] = G[i]->A ; + if (N_S[i] == GrB_NULL) + B_S[i] = GrB_NULL; + else + B_S[i] = N_S[i]->A; + } + LG_TRY(LAGraph_Malloc((void **)&B_call, num_nonterminals, sizeof(GrB_Matrix), msg)) + for (size_t i = 0; i < num_nonterminals; ++i) + { + if (Call_S[i] == GrB_NULL) + B_call[i] = GrB_NULL; + else + B_call[i] = Call_S[i]->A; + } + LG_TRY(LAGraph_Malloc((void **)&B_ret, num_nonterminals, sizeof(GrB_Matrix), msg)) + for (size_t i = 0; i < num_nonterminals; ++i) + { + if (Ret_S[i] == GrB_NULL) + B_ret[i] = GrB_NULL; + else + B_ret[i] = Ret_S[i]->A; } - LG_TRY (LAGraph_Malloc((void **) &B, nl, sizeof (GrB_Matrix), msg)) ; - LG_TRY (LAGraph_Malloc((void **) &BT, nl, sizeof (GrB_Matrix), msg)) ; + //-------------------------------------------------------------------------- + // determine ng (graph size) and nr (RSM state count) + //-------------------------------------------------------------------------- - for (size_t i = 0 ; i < nl ; ++i) + for (size_t i = 0; i < num_terminals; ++i) { - BT[i] = GrB_NULL ; - - if (R[i] == GrB_NULL) - { - B[i] = GrB_NULL ; - continue ; - } + if (A_a[i] == GrB_NULL) + continue; - B[i] = R[i]->A ; - if (R[i]->is_symmetric_structure == LAGraph_TRUE) - { - BT[i] = B[i] ; - } else - { - // BT[i] could be NULL and the matrix will be transposed by a - // descriptor - BT[i] = R[i]->AT ; - } + GRB_TRY(GrB_Matrix_nrows(&ng, A_a[i])); + break; } - + for (size_t i = 0; i < num_terminals; ++i) + { + if (B_a[i] == GrB_NULL) + continue; - LG_TRY (LAGraph_Malloc((void **) &B_call, 1, sizeof (GrB_Matrix), msg)) ; - LG_TRY (LAGraph_Malloc((void **) &BT_call, 1, sizeof (GrB_Matrix), msg)) ; + GRB_TRY(GrB_Matrix_nrows(&nr, B_a[i])); + break; + } - *B_call = (*R_call)->A ; - *BT_call = (*R_call)->AT ; + //-------------------------------------------------------------------------- + // dimension checks + //-------------------------------------------------------------------------- + for (size_t i = 0; i < num_terminals; ++i) + { + if (A_a[i] == GrB_NULL) + continue; + GRB_TRY(GrB_Matrix_nrows(&rows, A_a[i])); + GRB_TRY(GrB_Matrix_ncols(&cols, A_a[i])); - for (size_t i = 0 ; i < nl ; ++i) + LG_ASSERT_MSG(rows == ng && cols == ng, LAGRAPH_NOT_CACHED, + "all the matrices in the graph adjacency matrix decomposition " + "should have the same dimensions and be square"); + } + for (size_t i = 0; i < num_terminals; ++i) { - if (A[i] == GrB_NULL) continue ; + if (B_a[i] == GrB_NULL) + continue; - GRB_TRY (GrB_Matrix_nrows (&ng, A[i])) ; - break ; - } + GRB_TRY(GrB_Matrix_nrows(&rows, B_a[i])); + GRB_TRY(GrB_Matrix_ncols(&cols, B_a[i])); - for (size_t i = 0 ; i < nl ; ++i) + LG_ASSERT_MSG(rows == nr && cols == nr, LAGRAPH_NOT_CACHED, + "all the matrices in the RSM adjacency matrix decomposition " + "should have the same dimensions and be square") + } + for (size_t i = 0; i < num_nonterminals; ++i) { - if (B[i] == GrB_NULL) continue ; + if (B_S[i] == GrB_NULL) + continue; - GRB_TRY (GrB_Matrix_nrows (&nr, B[i])) ; - break ; - } + GRB_TRY(GrB_Matrix_nrows(&rows, B_S[i])); + GRB_TRY(GrB_Matrix_ncols(&cols, B_S[i])); - // Check all the matrices in graph adjacency matrix decomposition are - // square and of the same dimensions - for (size_t i = 0 ; i < nl ; ++i) + LG_ASSERT_MSG(rows == nr && cols == nr, LAGRAPH_NOT_CACHED, + "all the matrices in the RSM nonterminal matrix decomposition " + "should have the same dimensions and be square") + } + for (size_t i = 0; i < num_nonterminals; ++i) { - if (A[i] == GrB_NULL) continue ; + if (B_call[i] == GrB_NULL) + continue; - GRB_TRY (GrB_Matrix_nrows (&rows, A[i])) ; - GRB_TRY (GrB_Matrix_ncols (&cols, A[i])) ; + GRB_TRY(GrB_Matrix_nrows(&rows, B_call[i])); + GRB_TRY(GrB_Matrix_ncols(&cols, B_call[i])); - LG_ASSERT_MSG (rows == ng && cols == ng, LAGRAPH_NOT_CACHED, - "all the matrices in the graph adjacency matrix decomposition " - "should have the same dimensions and be square") ; + LG_ASSERT_MSG(rows == nr && cols == nr, LAGRAPH_NOT_CACHED, + "all the matrices in the RSM call matrix decomposition " + "should have the same dimensions and be square") } - - // Check all the matrices in RSM adjancency matrix decomposition are - // square and of the same dimensions - for (size_t i = 0 ; i < nl ; ++i) + for (size_t i = 0; i < num_nonterminals; ++i) { - if (B[i] == GrB_NULL) continue ; - - GRB_TRY (GrB_Matrix_nrows (&rows, B[i])) ; - GRB_TRY (GrB_Matrix_ncols (&cols, B[i])) ; - - LG_ASSERT_MSG (rows == nr && cols == nr, LAGRAPH_NOT_CACHED, - "all the matrices in the RSM adjacency matrix decomposition " - "should have the same dimensions and be square") ; - } - - // Check size of RSM call/return matrix + if (B_ret[i] == GrB_NULL) + continue; - GRB_TRY (GrB_Matrix_nrows (&rows, *B_call)) ; - GRB_TRY (GrB_Matrix_ncols (&cols, *B_call)) ; + GRB_TRY(GrB_Matrix_nrows(&rows, B_ret[i])); + GRB_TRY(GrB_Matrix_ncols(&cols, B_ret[i])); - LG_ASSERT_MSG (rows == nr && cols == nr, LAGRAPH_NOT_CACHED, - "RSM's call/return matrix should have nr*nr size") ; + LG_ASSERT_MSG(rows == nr && cols == nr, LAGRAPH_NOT_CACHED, + "all the matrices in the RSM return matrix decomposition " + "should have the same dimensions and be square") + } // Check source nodes in the graph - for (size_t i = 0 ; i < ns ; ++i) + for (size_t i = 0; i < ns; ++i) { - GrB_Index s = S[i] ; - LG_ASSERT_MSG(s < ng, GrB_INVALID_INDEX, "invalid graph source node") ; + GrB_Index s = Source[i]; + LG_ASSERT_MSG(s < ng, GrB_INVALID_INDEX, "invalid graph source node"); } - // Check starting states of the RSM - for (size_t i = 0 ; i < nqs ; ++i) + for (size_t i = 0; i < nqs; ++i) { - GrB_Index qs = QS[i] ; - LG_ASSERT_MSG (qs < nr, GrB_INVALID_INDEX, - "invalid RSM starting state") ; + GrB_Index qs = QS[i]; + LG_ASSERT_MSG(qs < nr, GrB_INVALID_INDEX, "invalid RSM starting state"); } - // Check final states of the RSM - for (size_t i = 0 ; i < nqf ; ++i) + for (size_t i = 0; i < nqf; ++i) { - GrB_Index qf = QF[i] ; - LG_ASSERT_MSG (qf < nr, GrB_INVALID_INDEX, "invalid RSM final state") ; + GrB_Index qf = QF[i]; + LG_ASSERT_MSG(qf < nr, GrB_INVALID_INDEX, "invalid RSM final state"); } + //-------------------------------------------------------------------------- + // allocate derived graph-edge matrices A_S[s] (ng x ng, initially empty) + //-------------------------------------------------------------------------- + + LG_TRY(LAGraph_Malloc((void **)&A_S, num_nonterminals, sizeof(GrB_Matrix), msg)) + for (size_t i = 0; i < num_nonterminals; ++i) + GRB_TRY(GrB_Matrix_new(&A_S[i], GrB_BOOL, ng, ng)); + // ------------------------------------------------------------------------- // initialization // ------------------------------------------------------------------------- - GRB_TRY (GrB_Vector_new (reachable, GrB_BOOL, ng)) ; - GRB_TRY (GrB_Vector_new (&final_reducer, GrB_BOOL, nr)) ; + GRB_TRY(GrB_Vector_new(reachable, GrB_BOOL, ng)); + GRB_TRY(GrB_Vector_new(&final_reducer, GrB_BOOL, nr)); + + // final_reducer[QF] = true + GrB_assign(final_reducer, GrB_NULL, GrB_NULL, true, QF, nqf, GrB_NULL); - // Initialize matrix for reducing the result - GrB_assign (final_reducer, GrB_NULL, GrB_NULL, true, QF, nqf, GrB_NULL) ; + // frontier / visited matrices (nr x ng) + GRB_TRY(GrB_Matrix_new(&frontier, GrB_BOOL, nr, ng)); + GRB_TRY(GrB_Matrix_new(&next_frontier, GrB_BOOL, nr, ng)); + GRB_TRY(GrB_Matrix_new(&front_temp1, GrB_BOOL, nr, ng)); + GRB_TRY(GrB_Matrix_new(&front_temp2, GrB_BOOL, nr, ng)); + GRB_TRY(GrB_Matrix_new(&P, GrB_BOOL, nr, ng)); + GRB_TRY(GrB_Matrix_new(&K, GrB_BOOL, nr, ng)) - GRB_TRY (GrB_Matrix_new (&frontier, GrB_BOOL, nr, ng)) ; - GRB_TRY (GrB_Matrix_new (&next_frontier, GrB_BOOL, nr, ng)) ; - GRB_TRY (GrB_Matrix_new (&symbol_frontier, GrB_BOOL, nr, ng)) ; - GRB_TRY (GrB_Matrix_new (&visited, GrB_BOOL, nr, ng)) ; + GRB_TRY(GrB_Matrix_new(&M_ret, GrB_BOOL, nr, ng)); + GRB_TRY(GrB_Matrix_new(&M_call, GrB_BOOL, nr, ng)); + GRB_TRY(GrB_Matrix_new(&M_call_new, GrB_BOOL, nr, ng)); + GRB_TRY(GrB_Matrix_new(&G_S_new, GrB_BOOL, ng, ng)); + GRB_TRY(GrB_Matrix_new(&step1, GrB_BOOL, ng, nr)); + + // Seed next_frontier, P, K from (start-states × source-vertices) + GrB_assign(next_frontier, GrB_NULL, GrB_NULL, true, QS, nqs, Source, ns, GrB_NULL); + GrB_assign(P, GrB_NULL, GrB_NULL, true, QS, nqs, Source, ns, GrB_NULL); + GrB_assign(K, GrB_NULL, GrB_NULL, true, QS, nqs, Source, ns, GrB_NULL); + + GRB_TRY(GrB_Matrix_nvals(&states, next_frontier)); + + int iteration = 0; - // Initialize frontier with the source nodes - GrB_assign (next_frontier, GrB_NULL, GrB_NULL, true, QS, nqs, S, ns, GrB_NULL) ; - GrB_assign (visited, GrB_NULL, GrB_NULL, true, QS, nqs, S, ns, GrB_NULL) ; - // Main loop while (states != 0) { - GrB_Matrix old_frontier = frontier ; - frontier = next_frontier ; - next_frontier = old_frontier ; + iteration++; + + GrB_Matrix old_frontier = frontier; + frontier = next_frontier; + next_frontier = old_frontier; + GRB_TRY(GrB_Matrix_clear(next_frontier)); - GRB_TRY (GrB_Matrix_clear (next_frontier)) ; + //---------------------------------------------------------------------- + // FORWARD PHASE + //---------------------------------------------------------------------- - // Obtain a new relation between the RSM states and the graph nodes - for (size_t i = 0 ; i < nl ; ++i) + // Terminals + for (size_t i = 0; i < num_terminals; ++i) { - if (A[i] == GrB_NULL || B[i] == GrB_NULL) continue ; + if (A_a[i] == GrB_NULL || B_a[i] == GrB_NULL) + continue; - // Traverse the RSM - // Try to use a provided transposed matrix or use the descriptor + GRB_TRY(GrB_Matrix_clear(front_temp1)); + GRB_TRY(GrB_Matrix_clear(front_temp2)); - GRB_TRY (GrB_Matrix_clear (symbol_frontier)) ; + // front_temp1 = N_a[a]^T * M + GRB_TRY(GrB_mxm(front_temp1, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, B_a[i], + frontier, GrB_DESC_T0)); - if (BT[i] != GrB_NULL) - { - GRB_TRY (GrB_mxm (symbol_frontier, GrB_NULL, GrB_LOR, - GrB_LOR_LAND_SEMIRING_BOOL, BT[i], frontier, GrB_NULL)) ; - } - else - { - GRB_TRY (GrB_mxm (symbol_frontier, GrB_NULL, GrB_LOR, - GrB_LOR_LAND_SEMIRING_BOOL, B[i], frontier, GrB_DESC_T0)) ; - } + // front_temp2 = front_temp1 * G_a[a] + GRB_TRY(GrB_mxm(front_temp2, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, front_temp1, + A_a[i], GrB_NULL)); - if ((*BT_call) != GrB_NULL) - { - GRB_TRY (GrB_mxm (next_frontier, visited, GrB_LOR, - GrB_LOR_LAND_SEMIRING_BOOL, *BT_call, frontier, GrB_DESC_SC)) ; // не уверен по поводу дескриптора - } - else + // M_new |= front_temp2 & ~P + GRB_TRY(GrB_eWiseAdd(next_frontier, P, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, + next_frontier, front_temp2, GrB_DESC_SC)); + } + + // Nonterminals (derived A_S edges) + for (size_t i = 0; i < num_nonterminals; ++i) + { + if (B_S[i] == GrB_NULL) + continue; + + GRB_TRY(GrB_Matrix_clear(front_temp1)); + GRB_TRY(GrB_Matrix_clear(front_temp2)); + + // front_temp1 = N_S[i]^T * frontier + GRB_TRY(GrB_mxm(front_temp1, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, B_S[i], + frontier, GrB_DESC_T0)); + + // front_temp2 = front_temp1 * A_S[i] + GRB_TRY(GrB_mxm(front_temp2, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, + front_temp1, A_S[i], GrB_NULL)); + + // next_frontier |= symbol_frontier & ~P + GRB_TRY(GrB_eWiseAdd(next_frontier, P, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, + next_frontier, front_temp2, GrB_DESC_SC)); + } + + // Call edges + for (size_t i = 0; i < num_nonterminals; ++i) + { + if (B_call[i] == GrB_NULL) + continue; + + GRB_TRY(GrB_Matrix_clear(front_temp1)); + + // front_temp1 = B_call[s]^T * frontier + GRB_TRY(GrB_mxm(front_temp1, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, B_call[i], + frontier, GrB_DESC_T0)); + + // next_frontier |= front_temp1 & ~P + GRB_TRY(GrB_eWiseAdd(next_frontier, P, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, + next_frontier, front_temp1, GrB_DESC_SC)); + } + + //---------------------------------------------------------------------- + // BACKWARD PHASE – match return edges back to their call sites + //---------------------------------------------------------------------- + + for (size_t i = 0; i < num_nonterminals; ++i) + { + if (B_ret[i] == GrB_NULL || B_call[i] == GrB_NULL || B_S[i] == GrB_NULL) + continue; + + GRB_TRY(GrB_Matrix_clear(M_ret)); + GRB_TRY(GrB_Matrix_clear(M_call)); + + // M_ret = B_ret[s]^T * frontier (positions that hit a return state) + GRB_TRY(GrB_mxm(M_ret, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, B_ret[i], + frontier, GrB_DESC_T0)); + + GrB_Index ret_nvals; + GRB_TRY(GrB_Matrix_nvals(&ret_nvals, M_ret)); + if (ret_nvals == 0) + continue; + + // M_call = (B_ret[s] * M_ret) & frontier + GRB_TRY(GrB_mxm(M_call, frontier, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, B_ret[i], M_ret, + GrB_DESC_S)); + + // walk backwards through visited set P to find the + // matching call-site positions. + + GrB_Index call_nvals; + GRB_TRY(GrB_Matrix_nvals(&call_nvals, M_call)); + while (call_nvals > 0) { - GRB_TRY (GrB_mxm (next_frontier, visited, GrB_LOR, - GrB_LOR_LAND_SEMIRING_BOOL, *B_call, frontier, GrB_DESC_SCT0)) ; // не уверен по поводу дескриптора + GRB_TRY(GrB_Matrix_clear(M_call_new)); + + // Backward terminal step: M_call_new |= (N_a[i] * M_call * G_a[i]^T) & P + for (size_t i = 0; i < num_terminals; ++i) + { + if (A_a[i] == GrB_NULL || B_a[i] == GrB_NULL) + continue; + + GRB_TRY(GrB_Matrix_clear(front_temp1)); + GRB_TRY(GrB_Matrix_clear(front_temp2)); + + // front_temp1 = N_a[a] * M_call + GRB_TRY(GrB_mxm(front_temp1, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, + B_a[i], M_call, GrB_NULL)); + + // front_temp2 = front_temp1 * G_a[a]^T + GRB_TRY(GrB_mxm(front_temp2, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, + front_temp1, A_a[i], GrB_DESC_T1)); + + // M_call_new |= front_temp2 & P + GRB_TRY(GrB_eWiseAdd(M_call_new, P, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, + M_call_new, front_temp2, GrB_DESC_S)); + } + // Backward nonterminal step: M_call_new |= (N_S[i] * M_call * G_S[i]^T) & P + for (size_t i = 0; i < num_nonterminals; ++i) + { + if (A_S[i] == GrB_NULL) + continue; + + GRB_TRY(GrB_Matrix_clear(front_temp1)); + GRB_TRY(GrB_Matrix_clear(front_temp2)); + + // front_temp1 = N_S[S] * M_call + GRB_TRY(GrB_mxm(front_temp1, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, + B_S[i], M_call, GrB_NULL)); + + // front_temp2 = front_temp1 * G_S[S]^T + GRB_TRY(GrB_mxm(front_temp2, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, + front_temp1, A_S[i], GrB_DESC_T1)); + + // M_call_new |= front_temp2 & P + GRB_TRY(GrB_eWiseAdd(M_call_new, P, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, + M_call_new, front_temp2, GrB_DESC_S)); + } + + GRB_TRY(GrB_Matrix_nvals(&call_nvals, M_call_new)); + if (call_nvals == 0) + break; + + GRB_TRY(GrB_Matrix_clear(M_call)); + GRB_TRY(GrB_eWiseAdd(M_call, GrB_NULL, GrB_NULL, GrB_LOR, M_call, M_call_new, + GrB_NULL)); } - // я реально не понимаю как это работает, но окей - // Traverse the graph - GRB_TRY (GrB_mxm (next_frontier, visited, GrB_LOR, - GrB_LOR_LAND_SEMIRING_BOOL, symbol_frontier, A[i], GrB_DESC_SC)) ; + // need to check if we reach call pos + // M_call = (B_call[s] * M_call) & P + GRB_TRY(GrB_Matrix_clear(front_temp1)); + GRB_TRY(GrB_mxm(front_temp1, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, B_call[i], + M_call, GrB_NULL)); + GRB_TRY(GrB_Matrix_clear(M_call)); + GRB_TRY(GrB_eWiseAdd(M_call, P, GrB_NULL, GrB_LOR, M_call, front_temp1, GrB_DESC_S)); + + GRB_TRY(GrB_Matrix_nvals(&call_nvals, M_call)); + if (call_nvals == 0) + continue; + + // we found some call positions + // but still not sure if they are correct ones + // G_S_new = (M_call.T @ N_S[S] @ M_ret) & ~G_S[S] + GRB_TRY(GrB_Matrix_clear(step1)); + GRB_TRY(GrB_Matrix_clear(front_temp2)); + GRB_TRY(GrB_Matrix_clear(G_S_new)); + + GRB_TRY(GrB_mxm(step1, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, M_call, B_S[i], + GrB_DESC_T0)); // step1 is ng*nr + + GRB_TRY(GrB_Matrix_clear(G_S_new)); + GRB_TRY(GrB_mxm(G_S_new, A_S[i], GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, step1, M_ret, + GrB_DESC_SC)); + + GrB_Index new_edge_nvals; + GRB_TRY(GrB_Matrix_nvals(&new_edge_nvals, G_S_new)); + if (new_edge_nvals == 0) + continue; + + /* + call positions are correct + we derived some edges for graph + now we need to traverse this part again */ + + // G_S[S] |= G_S_new + GRB_TRY(GrB_eWiseAdd(A_S[i], GrB_NULL, GrB_NULL, GrB_LOR, A_S[i], G_S_new, GrB_NULL)); + + // P &= ~M_ret + GRB_TRY(GrB_Matrix_clear(front_temp1)); + GRB_TRY( + GrB_eWiseAdd(front_temp1, GrB_NULL, GrB_NULL, GrB_LOR, front_temp1, P, GrB_NULL)); + GRB_TRY(GrB_Matrix_clear(P)); + GRB_TRY(GrB_eWiseAdd(P, M_ret, GrB_NULL, GrB_LOR, P, front_temp1, GrB_DESC_SC)); + + // next_frontier |= M_call + GRB_TRY(GrB_eWiseAdd(next_frontier, GrB_NULL, GrB_NULL, GrB_LOR, next_frontier, M_call, + GrB_NULL)); + } + + //---------------------------------------------------------------------- + // ADVANCE: P |= next_frontier + //---------------------------------------------------------------------- + + GRB_TRY(GrB_eWiseAdd(P, GrB_NULL, GrB_NULL, GrB_LOR, P, next_frontier, GrB_NULL)); + GRB_TRY(GrB_Matrix_nvals(&states, next_frontier)); + + if (iteration > 3000) + { + printf("Warning: Maximum iterations reached\n"); + break; } + } + + GRB_TRY(GrB_Matrix_clear(next_frontier)); + GrB_assign(next_frontier, GrB_NULL, GrB_NULL, true, QS, nqs, Source, ns, GrB_NULL); + + GRB_TRY(GrB_Matrix_nvals(&states, next_frontier)); + // second loop + while (states != 0) + { + iteration++; + GrB_Matrix old_frontier = frontier; + frontier = next_frontier; + next_frontier = old_frontier; + GRB_TRY(GrB_Matrix_clear(next_frontier)); + + // Terminals + for (size_t i = 0; i < num_terminals; ++i) + { + if (A_a[i] == GrB_NULL || B_a[i] == GrB_NULL) + continue; + + GRB_TRY(GrB_Matrix_clear(front_temp1)); + GRB_TRY(GrB_Matrix_clear(front_temp2)); - // я реально не понимаю как это работает, но окей - // Accumulate the new state <-> node correspondence - GRB_TRY (GrB_assign (visited, visited, GrB_NULL, next_frontier, - GrB_ALL, nr, GrB_ALL, ng, GrB_DESC_SC)) ; + // front_temp1 = N_a[a]^T * M + GRB_TRY(GrB_mxm(front_temp1, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, B_a[i], + frontier, GrB_DESC_T0)); - GRB_TRY (GrB_Matrix_nvals (&states, next_frontier)) ; + // front_temp2 = front_temp1 * G_a[a] + GRB_TRY(GrB_mxm(front_temp2, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, front_temp1, + A_a[i], GrB_NULL)); + + // M_new |= front_temp2 & ~K + GRB_TRY(GrB_eWiseAdd(next_frontier, K, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, + next_frontier, front_temp2, GrB_DESC_SC)); + } + + // Nonterminals (derived A_S edges) + for (size_t i = 0; i < num_nonterminals; ++i) + { + if (B_S[i] == GrB_NULL) + continue; + + GRB_TRY(GrB_Matrix_clear(front_temp1)); + GRB_TRY(GrB_Matrix_clear(front_temp2)); + + // front_temp1 = N_S[i]^T * frontier + GRB_TRY(GrB_mxm(front_temp1, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, B_S[i], + frontier, GrB_DESC_T0)); + + // front_temp2 = front_temp1 * G_S[i] + GRB_TRY(GrB_mxm(front_temp2, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, + front_temp1, A_S[i], GrB_NULL)); + + // next_frontier |= symbol_frontier & ~K + GRB_TRY(GrB_eWiseAdd(next_frontier, K, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, + next_frontier, front_temp2, GrB_DESC_SC)); + } + + //---------------------------------------------------------------------- + // ADVANCE: P |= next_frontier + //---------------------------------------------------------------------- + + GRB_TRY(GrB_eWiseAdd(K, GrB_NULL, GrB_NULL, GrB_LOR, K, next_frontier, GrB_NULL)); + GRB_TRY(GrB_Matrix_nvals(&states, next_frontier)); + + if (iteration > 3000) + { + printf("Warning: Maximum iterations reached\n"); + break; + } } - // Extract the nodes matching the final RSM states - GRB_TRY (GrB_vxm (*reachable, GrB_NULL, GrB_NULL, - GrB_LOR_LAND_SEMIRING_BOOL, final_reducer, visited, GrB_NULL)) ; + GRB_TRY(GrB_vxm(*reachable, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, final_reducer, K, + GrB_NULL)) - LG_FREE_WORK ; - return (GrB_SUCCESS) ; + LG_FREE_WORK; + return GrB_SUCCESS; } diff --git a/experimental/test/test_RSM_reachability.c b/experimental/test/test_RSM_reachability.c index 121dc511f7..bf10a389de 100644 --- a/experimental/test/test_RSM_reachability.c +++ b/experimental/test/test_RSM_reachability.c @@ -1,243 +1,333 @@ -#include "GraphBLAS.h" -#include "LAGraph.h" +#include +#include #include -#include #include -#include +#include + +#include "LAGraph.h" + #include #include +#include +#include #define LEN 512 #define MAX_LABELS 3 #define MAX_RESULTS 2000000 -LAGraph_Graph G[MAX_LABELS] ; -LAGraph_Graph R[MAX_LABELS] ; -LAGraph_Graph R_call; - -GrB_Matrix A ; +#define CHECK(info) \ + do \ + { \ + GrB_Info _i = (info); \ + if (_i != GrB_SUCCESS && _i != GrB_NO_VALUE) \ + { \ + fprintf(stderr, "GraphBLAS error %d at %s:%d\n", _i, __FILE__, __LINE__); \ + exit(EXIT_FAILURE); \ + } \ + } while (0) -char testcase_name[LEN+1] ; -char filename[LEN+1] ; -char msg[LAGRAPH_MSG_LEN] ; +char msg[LAGRAPH_MSG_LEN]; -typedef struct +static void setup(void) { - const char* name ; - const char* graphs[MAX_LABELS] ; - const char* rsm[MAX_LABELS] ; - const char* rsm_call; - const char* rsm_meta ; - const char* sources ; - const GrB_Index expected[MAX_RESULTS] ; - const size_t expected_count ; + LAGraph_Init(msg); } -matrix_info ; -const matrix_info files [ ] = +static void teardown(void) { - {"simple 1", - {"rsm_data/a.mtx", "rsm_data/b.mtx", NULL}, - {"rsm_data/1_a.mtx", NULL}, - "rsm_data/1_call.mtx", - "rsm_data/1_meta.txt", - "rsm_data/1_sources.txt", - {2}, - 1}, - {NULL, NULL, NULL, NULL}, -} ; - - -void test_RSM_reachability (void) -{ - LAGraph_Init (msg) ; - - for (int k = 0 ; ; ++k) - { - if (files[k].sources == NULL) break ; + LAGraph_Finalize(msg); +} - snprintf (testcase_name, LEN, "basic context-free path query %s", files[k].name) ; - TEST_CASE (testcase_name) ; +static GrB_Matrix make_bool_matrix(GrB_Index rows, GrB_Index cols, const GrB_Index *I, + const GrB_Index *J, GrB_Index nvals) +{ + GrB_Matrix A = GrB_NULL; + TEST_CHECK(GrB_SUCCESS == GrB_Matrix_new(&A, GrB_BOOL, rows, cols)); - for (int check_symmetry = 0 ; check_symmetry < 2 ; ++check_symmetry) - { - // Load graph from MTX files represention its adjacency matrix - // decomposition - for (int i = 0 ; ; ++i) - { - const char *name = files[k].graphs[i] ; - - if (name == NULL) break ; - if (strlen(name) == 0) continue ; - - snprintf (filename, LEN, LG_DATA_DIR "%s", name) ; - FILE *f = fopen (filename, "r") ; - TEST_CHECK (f != NULL) ; - OK (LAGraph_MMRead (&A, f, msg)) ; - OK (fclose (f)) ; - - OK (LAGraph_New (&(G[i]), &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; - - TEST_CHECK (A == NULL) ; - } - - // Load RSM from MTX files representing its adjacency matrix - // decomposition - for (int i = 0 ; ; ++i) - { - const char *name = files[k].rsm[i] ; - - if (name == NULL) break ; - if (strlen(name) == 0) continue ; - - snprintf (filename, LEN, LG_DATA_DIR "%s", name) ; - FILE *f = fopen (filename, "r") ; - TEST_CHECK (f != NULL) ; - OK (LAGraph_MMRead (&A, f, msg)) ; - OK (fclose (f)) ; - - OK (LAGraph_New (&(R[i]), &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; - - if (check_symmetry) - { - // Check if the pattern is symmetric - if it isn't make it. - // Note this also computes R[i]->AT - OK (LAGraph_Cached_IsSymmetricStructure (R[i], msg)) ; - } - - TEST_CHECK (A == NULL) ; - } - - { - const char *name = files[k].rsm_call ; - - snprintf (filename, LEN, LG_DATA_DIR "%s", name) ; - FILE *f = fopen (filename, "r") ; - TEST_CHECK (f != NULL) ; - OK (LAGraph_MMRead (&A, f, msg)) ; - OK (fclose (f)) ; - - OK (LAGraph_New (&(R_call), &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; - - if (check_symmetry) - { - OK (LAGraph_Cached_IsSymmetricStructure(R_call, msg)) ; - } - - TEST_CHECK (A == NULL) ; // зачем нам это проверять?, почему оно NULL - } + for (GrB_Index k = 0; k < nvals; k++) + TEST_CHECK(GrB_SUCCESS == GrB_Matrix_setElement_BOOL(A, true, I[k], J[k])); + return A; +} - } +static LAGraph_Graph make_graph(GrB_Matrix *pA) +{ + LAGraph_Graph G = GrB_NULL; + int retval = LAGraph_New(&G, pA, LAGraph_ADJACENCY_DIRECTED, msg); + TEST_CHECK(retval == 0); + TEST_MSG("LAGraph_New failed: %s", msg); + return G; +} - // Note the matrix rows/cols are enumerated from 0 to n-1. - // Meanwhile, in MTX format they are enumerated from 1 to n. Thus, - // when loading/comparing the results these values should be - // decremented/incremented correspondingly. +static void free_graphs(LAGraph_Graph *arr, size_t n) +{ + for (size_t i = 0; i < n; ++i) + if (arr[i] != GrB_NULL) + LAGraph_Delete(&arr[i], msg); +} - // Load graph source nodes from the sources file - GrB_Index s ; - GrB_Index S[16] ; - size_t ns = 0 ; +static void check_result(GrB_Vector result, GrB_Index V, const bool *expected, + const char *test_name) +{ + GrB_Index n = 0; + TEST_CHECK(GrB_SUCCESS == GrB_Vector_size(&n, result)); + TEST_CHECK_(n == V, "%s: result vector size %llu != %llu", test_name, (unsigned long long)n, + (unsigned long long)V); - const char *name = files[k].sources ; - snprintf (filename, LEN, LG_DATA_DIR "%s", name) ; - FILE *f = fopen (filename, "r") ; - TEST_CHECK (f != NULL) ; + for (GrB_Index v = 0; v < V; v++) + { + bool val = false; + GrB_Info info = GrB_Vector_extractElement_BOOL(&val, result, v); - while (fscanf(f, "%" PRIu64, &s) != EOF) + if (expected[v]) { - S[ns++] = s - 1 ; + TEST_CHECK_(info == GrB_SUCCESS && val, + "%s: vertex %llu should be reachable but is not", test_name, + (unsigned long long)v); } + else + { + bool absent = (info == GrB_NO_VALUE) || (!val); + TEST_CHECK_(absent, "%s: vertex %llu should not be reachable but is", test_name, + (unsigned long long)v); + } + } +} - OK (fclose(f)) ; - - // Load RSM starting states from the meta file - GrB_Index qs ; - GrB_Index QS[16] ; // почему именно 16, мб личный выбор - size_t nqs = 0 ; +//============================================================================== +// TEST CASES +//============================================================================== - name = files[k].rsm_meta ; - snprintf (filename, LEN, LG_DATA_DIR "%s", name) ; - f = fopen (filename, "r") ; - TEST_CHECK (f != NULL) ; ; +void test_based(void) +{ + setup(); - uint64_t nqs64 = 0 ; - TEST_CHECK (fscanf (f, "%" PRIu64, &nqs64) != EOF) ; - nqs = (size_t) nqs64 ; + const GrB_Index Q = 4, V = 6; - for (size_t i = 0 ; i < nqs ; ++i) - { - TEST_CHECK(fscanf (f, "%" PRIu64, &qs) != EOF) ; - QS[i] = qs - 1 ; - } + // Graph: 0 -a-> 1 -a-> 2 -b-> 3 -b-> 4 -b-> 5 + GrB_Index Ga_I[] = {0, 1}, Ga_J[] = {1, 2}; + GrB_Index Gb_I[] = {2, 3, 4}, Gb_J[] = {3, 4, 5}; - // Load RSM final states from the same file - uint64_t qf ; - uint64_t QF[16] ; - size_t nqf = 0 ; - uint64_t nqf64 = 0 ; + // RSM: S: q0 -a-> q1 -S-> q2 -b-> q3 - TEST_CHECK (fscanf (f, "%" PRIu64, &nqf64) != EOF) ; - nqf = (size_t) nqf64 ; + GrB_Index Na_I[] = {0}, Na_J[] = {1}; + GrB_Index Nb_I[] = {1, 2}, Nb_J[] = {3, 3}; // q1->q3, q2->q3 - for (size_t i = 0 ; i < nqf ; ++i) - { - TEST_CHECK (fscanf (f, "%" PRIu64, &qf) != EOF) ; - QF[i] = qf - 1 ; - } + GrB_Index Call_I[] = {1}, Call_J[] = {0}; + GrB_Index Ret_I[] = {3}, Ret_J[] = {2}; + GrB_Index Ns_I[] = {1}, Ns_J[] = {2}; // N_S: q1->q2 - OK (fclose (f)) ; + GrB_Matrix Ga_raw = make_bool_matrix(V, V, Ga_I, Ga_J, 2); + GrB_Matrix Gb_raw = make_bool_matrix(V, V, Gb_I, Gb_J, 3); + GrB_Matrix Na_raw = make_bool_matrix(Q, Q, Na_I, Na_J, 1); + GrB_Matrix Nb_raw = make_bool_matrix(Q, Q, Nb_I, Nb_J, 2); + GrB_Matrix Cs_raw = make_bool_matrix(Q, Q, Call_I, Call_J, 1); + GrB_Matrix Rs_raw = make_bool_matrix(Q, Q, Ret_I, Ret_J, 1); + GrB_Matrix Ns_raw = make_bool_matrix(Q, Q, Ns_I, Ns_J, 1); - // Evaluate the algorithm - GrB_Vector r = NULL ; + LAGraph_Graph Ga[2], Na[2], Cs[1], Rs[1], Ns[1]; + Ga[0] = make_graph(&Ga_raw); + Ga[1] = make_graph(&Gb_raw); + Na[0] = make_graph(&Na_raw); + Na[1] = make_graph(&Nb_raw); + Cs[0] = make_graph(&Cs_raw); + Rs[0] = make_graph(&Rs_raw); + Ns[0] = make_graph(&Ns_raw); - OK (LAGraph_RSM_reachability (&r, R, MAX_LABELS, QS, nqs, - QF, nqf, &R_call, G, S, ns, msg)) ; - - // Extract results from the output vector - GrB_Index *reachable ; - bool *values ; - - GrB_Index nvals ; - GrB_Vector_nvals (&nvals, r) ; + GrB_Index QS[] = {0}, QF[] = {3}, Source[] = {0}; + GrB_Vector result = GrB_NULL; - OK (LAGraph_Malloc ((void **) &reachable, MAX_RESULTS, sizeof (GrB_Index), msg)) ;; - OK (LAGraph_Malloc ((void **) &values, MAX_RESULTS, sizeof (bool), msg)) ; + int retval = + LAGraph_RSM_reachability(&result, 2, Ga, Na, 1, Ns, Cs, Rs, QS, 1, QF, 1, Source, 1, msg); - GrB_Vector_extractTuples (reachable, values, &nvals, r) ; - - TEST_MSG("returned %lu values:\n", nvals); - for (uint64_t i = 0; i < nvals; i++) - TEST_MSG(" %lu\n", reachable[i] + 1); + TEST_CHECK(GrB_SUCCESS == retval); + TEST_MSG("test_based: retval = %d (%s)", retval, msg); - // Compare the results with expected values - TEST_CHECK (nvals == files[k].expected_count) ; - for (uint64_t i = 0 ; i < nvals ; ++i) - TEST_CHECK (reachable[i] + 1 == files[k].expected[i]) ; + const bool expected[6] = {false, false, false, false, true, false}; + check_result(result, V, expected, "test_based"); - // Cleanup - OK (LAGraph_Free ((void **) &values, NULL)) ; - OK (LAGraph_Free ((void **) &reachable, NULL)) ; + GrB_free(&result); + free_graphs(Ga, 2); + free_graphs(Na, 2); + free_graphs(Cs, 1); + free_graphs(Rs, 1); + free_graphs(Ns, 1); - OK (GrB_free (&r)) ; + teardown(); +} - for (uint64_t i = 0 ; i < MAX_LABELS ; ++i) - { - if (G[i] == NULL) continue ; - OK (LAGraph_Delete (&(G[i]), msg)) ; - } +void test_balanced_paren(void) +{ + setup(); + + const GrB_Index Q = 5, V = 6; + + // Graph: 0 -a-> 1 -a-> 2 -b-> 3 -b-> 4 -b-> 5 + GrB_Index Ga_I[] = {0, 1, 4}, Ga_J[] = {1, 2, 5}; + GrB_Index Gb_I[] = {2, 3}, Gb_J[] = {3, 4}; + + // RSM: S: q0 -(-> q1 -S-> q2 -)-> q3 -S-> q4 + GrB_Index Na_I[] = {0}, Na_J[] = {1}; + GrB_Index Nb_I[] = {1, 1, 2, 2}, Nb_J[] = {3, 4, 3, 4}; + + GrB_Index Call_I[] = {1, 3}, Call_J[] = {0, 0}; + GrB_Index Ret_I[] = {4, 4}, Ret_J[] = {2, 4}; + GrB_Index Ns_I[] = {1, 3}, Ns_J[] = {2, 4}; + + GrB_Matrix Ga_raw = make_bool_matrix(V, V, Ga_I, Ga_J, 3); + GrB_Matrix Gb_raw = make_bool_matrix(V, V, Gb_I, Gb_J, 2); + GrB_Matrix Na_raw = make_bool_matrix(Q, Q, Na_I, Na_J, 1); + GrB_Matrix Nb_raw = make_bool_matrix(Q, Q, Nb_I, Nb_J, 4); + GrB_Matrix Cs_raw = make_bool_matrix(Q, Q, Call_I, Call_J, 2); + GrB_Matrix Rs_raw = make_bool_matrix(Q, Q, Ret_I, Ret_J, 2); + GrB_Matrix Ns_raw = make_bool_matrix(Q, Q, Ns_I, Ns_J, 2); + + LAGraph_Graph Ga[2], Na[2], Cs[1], Rs[1], Ns[1]; + Ga[0] = make_graph(&Ga_raw); + Ga[1] = make_graph(&Gb_raw); + Na[0] = make_graph(&Na_raw); + Na[1] = make_graph(&Nb_raw); + Cs[0] = make_graph(&Cs_raw); + Rs[0] = make_graph(&Rs_raw); + Ns[0] = make_graph(&Ns_raw); + + GrB_Index QS[] = {0}, QF[] = {4}, Source[] = {0}; + GrB_Vector result = GrB_NULL; + + int retval = + LAGraph_RSM_reachability(&result, 2, Ga, Na, 1, Ns, Cs, Rs, QS, 1, QF, 1, Source, 1, msg); + + TEST_CHECK(GrB_SUCCESS == retval); + TEST_MSG("test_based: retval = %d (%s)", retval, msg); + + const bool expected[6] = {false, false, false, false, true, false}; + check_result(result, V, expected, "test_based"); + + GrB_free(&result); + free_graphs(Ga, 2); + free_graphs(Na, 2); + free_graphs(Cs, 1); + free_graphs(Rs, 1); + free_graphs(Ns, 1); + + teardown(); +} - for (uint64_t i = 0; i < MAX_LABELS ; ++i) - { - if (R[i] == NULL) continue ; - OK (LAGraph_Delete (&(R[i]), msg)) ; - } - } +void parens_lrlrlr(void) +{ + setup(); + + const GrB_Index Q = 5, V = 7; + + // Graph 0 -(-> 1 -)-> 2 -(-> 3 -)-> 4 -(-> 5 -)-> 6 + GrB_Index Ga_I[] = {0, 2, 4}, Ga_J[] = {1, 3, 5}; + GrB_Index Gb_I[] = {1, 3, 5}, Gb_J[] = {2, 4, 6}; + + // RSM: S: q0 -(-> q1 -S-> q2 -)-> q3 -S-> q4 + GrB_Index Na_I[] = {0}, Na_J[] = {1}; + GrB_Index Nb_I[] = {1, 1, 2, 2}, Nb_J[] = {3, 4, 3, 4}; + + GrB_Index Call_I[] = {1, 3}, Call_J[] = {0, 0}; + GrB_Index Ret_I[] = {4, 4}, Ret_J[] = {2, 4}; + GrB_Index Ns_I[] = {1, 3}, Ns_J[] = {2, 4}; + + GrB_Matrix Ga_raw = make_bool_matrix(V, V, Ga_I, Ga_J, 3); + GrB_Matrix Gb_raw = make_bool_matrix(V, V, Gb_I, Gb_J, 3); + GrB_Matrix Na_raw = make_bool_matrix(Q, Q, Na_I, Na_J, 1); + GrB_Matrix Nb_raw = make_bool_matrix(Q, Q, Nb_I, Nb_J, 4); + GrB_Matrix Cs_raw = make_bool_matrix(Q, Q, Call_I, Call_J, 2); + GrB_Matrix Rs_raw = make_bool_matrix(Q, Q, Ret_I, Ret_J, 2); + GrB_Matrix Ns_raw = make_bool_matrix(Q, Q, Ns_I, Ns_J, 2); + + LAGraph_Graph Ga[2], Na[2], Cs[1], Rs[1], Ns[1]; + Ga[0] = make_graph(&Ga_raw); + Ga[1] = make_graph(&Gb_raw); + Na[0] = make_graph(&Na_raw); + Na[1] = make_graph(&Nb_raw); + Cs[0] = make_graph(&Cs_raw); + Rs[0] = make_graph(&Rs_raw); + Ns[0] = make_graph(&Ns_raw); + + GrB_Index QS[] = {0}, QF[] = {4}, Source[] = {0}; + GrB_Vector result = GrB_NULL; + + int retval = + LAGraph_RSM_reachability(&result, 2, Ga, Na, 1, Ns, Cs, Rs, QS, 1, QF, 1, Source, 1, msg); + + TEST_CHECK(GrB_SUCCESS == retval); + TEST_MSG("test_based: retval = %d (%s)", retval, msg); + + const bool expected[7] = {0, 0, 1, 0, 1, 0, 1}; + check_result(result, V, expected, "test_based"); + + GrB_free(&result); + free_graphs(Ga, 2); + free_graphs(Na, 2); + free_graphs(Cs, 1); + free_graphs(Rs, 1); + free_graphs(Ns, 1); + + teardown(); +} - LAGraph_Finalize (msg) ; +void very_big_parens(void) +{ + setup(); + + const GrB_Index Q = 5, V = 13; + + // Graph (()(()()))() + GrB_Index Ga_I[] = {0, 1, 3, 4, 6, 10}, Ga_J[] = {1, 2, 4, 5, 7, 11}; + GrB_Index Gb_I[] = {2, 5, 7, 8, 9, 11}, Gb_J[] = {3, 6, 8, 9, 10, 12}; + + // RSM: S: q0 -(-> q1 -S-> q2 -)-> q3 -S-> q4 + GrB_Index Na_I[] = {0}, Na_J[] = {1}; + GrB_Index Nb_I[] = {1, 1, 2, 2}, Nb_J[] = {3, 4, 3, 4}; + + GrB_Index Call_I[] = {1, 3}, Call_J[] = {0, 0}; + GrB_Index Ret_I[] = {4, 4}, Ret_J[] = {2, 4}; + GrB_Index Ns_I[] = {1, 3}, Ns_J[] = {2, 4}; + + GrB_Matrix Ga_raw = make_bool_matrix(V, V, Ga_I, Ga_J, 6); + GrB_Matrix Gb_raw = make_bool_matrix(V, V, Gb_I, Gb_J, 6); + GrB_Matrix Na_raw = make_bool_matrix(Q, Q, Na_I, Na_J, 1); + GrB_Matrix Nb_raw = make_bool_matrix(Q, Q, Nb_I, Nb_J, 4); + GrB_Matrix Cs_raw = make_bool_matrix(Q, Q, Call_I, Call_J, 2); + GrB_Matrix Rs_raw = make_bool_matrix(Q, Q, Ret_I, Ret_J, 2); + GrB_Matrix Ns_raw = make_bool_matrix(Q, Q, Ns_I, Ns_J, 2); + + LAGraph_Graph Ga[2], Na[2], Cs[1], Rs[1], Ns[1]; + Ga[0] = make_graph(&Ga_raw); + Ga[1] = make_graph(&Gb_raw); + Na[0] = make_graph(&Na_raw); + Na[1] = make_graph(&Nb_raw); + Cs[0] = make_graph(&Cs_raw); + Rs[0] = make_graph(&Rs_raw); + Ns[0] = make_graph(&Ns_raw); + + GrB_Index QS[] = {0}, QF[] = {4}, Source[] = {0}; + GrB_Vector result = GrB_NULL; + + int retval = + LAGraph_RSM_reachability(&result, 2, Ga, Na, 1, Ns, Cs, Rs, QS, 1, QF, 1, Source, 1, msg); + + TEST_CHECK(GrB_SUCCESS == retval); + TEST_MSG("test_based: retval = %d (%s)", retval, msg); + + const bool expected[13] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1}; + check_result(result, V, expected, "test_based"); + + GrB_free(&result); + free_graphs(Ga, 2); + free_graphs(Na, 2); + free_graphs(Cs, 1); + free_graphs(Rs, 1); + free_graphs(Ns, 1); + + teardown(); } -TEST_LIST = { - {"RSM_reachability", test_RSM_reachability}, - {NULL, NULL} -} ; +TEST_LIST = {{"simple", test_based}, + {"test_balanced_paren", test_balanced_paren}, + {"parens_lrlrlr", parens_lrlrlr}, + {"very_big_parens", very_big_parens}, + {NULL, NULL}}; diff --git a/include/LAGraphX.h b/include/LAGraphX.h index d7053d0244..daf5d5727d 100644 --- a/include/LAGraphX.h +++ b/include/LAGraphX.h @@ -858,31 +858,31 @@ int LAGraph_RegularPathQuery // nodes reachable from the starting by the //**************************************************************************** LAGRAPHX_PUBLIC -int LAGraph_RSM_reachability -( +int LAGraph_RSM_reachability( // output: - GrB_Vector *reachable, // reachable(i) = true if node i is reachable - // from one of the starting nodes by a path - // satisfying regular constraints - - LAGraph_Graph *R, // input recursive state machine's - // adjacency matrix decomposition - size_t nl, // total label count, # of matrices graph and - // RSM adjacency matrix decomposition - const GrB_Index *QS, // starting states in RSM - size_t nqs, // number of starting states in RSM - const GrB_Index *QF, // final states in RSM - size_t nqf, // number of final states in RSM - - LAGraph_Graph *R_call, // recursive state machine's call/return nr*nr - // matrix, where (i, j) = 1 if there is transition - // from i state of one automata to j starting - // state of another automata - LAGraph_Graph *G, // input graph adjacency matrix decomposition - const GrB_Index *S, // source vertices to start searching paths - size_t ns, // number of source vertices - char *msg // LAGraph output message -) ; + GrB_Vector *reachable, + + // input: + size_t num_terminals, // # terminal labels + LAGraph_Graph *G_a, // terminal transitions + LAGraph_Graph *N_a, // terminal transitions + + size_t num_nonterminals, // # nonterminal labels + LAGraph_Graph *N_S, // nonterminal RSM transitions + LAGraph_Graph *Call_S, // call transition matrix + LAGraph_Graph *Ret_S, // return transition matrix + + const GrB_Index *QS, // staring states in RSM + size_t nqs, // # starting states + const GrB_Index *QF, // final states in RSM + size_t nqf, // # final states + + const GrB_Index *Source, // sources vertices + size_t ns, // # source vertices + + char *msg + +); //**************************************************************************** LAGRAPHX_PUBLIC