From a6e4dbe1a0658f7eba0a0a39d4f5b02c5b19482e Mon Sep 17 00:00:00 2001 From: Sterkachov Date: Wed, 3 Dec 2025 09:28:14 +0300 Subject: [PATCH 01/10] 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 02/10] 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 03/10] 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 From 767268b75eb8c0759fd476fc0db9590fdad2f76d Mon Sep 17 00:00:00 2001 From: Sterkachov Date: Tue, 3 Mar 2026 16:22:50 +0300 Subject: [PATCH 04/10] feat: add the algorithm's description --- .../algorithm/LAGraph_RSM_reachability.c | 112 ++++++++++++++++++ 1 file changed, 112 insertions(+) diff --git a/experimental/algorithm/LAGraph_RSM_reachability.c b/experimental/algorithm/LAGraph_RSM_reachability.c index 8972116a4c..cf6faf523d 100644 --- a/experimental/algorithm/LAGraph_RSM_reachability.c +++ b/experimental/algorithm/LAGraph_RSM_reachability.c @@ -2,6 +2,118 @@ // LAGraph_RSM_reachability.c //------------------------------------------------------------------------------ +/* +Let: + + Q = {0, …, Q-1} — set of RSM states + V = {0, …, V-1} — set of graph vertices + Σ — terminal alphabet + N — set of nonterminals + +Boolean matrices: + + G_a ∈ {0,1}^{V×V} — graph edges labeled by a ∈ Σ + N_a ∈ {0,1}^{Q×Q} — RSM transitions labeled by a ∈ Σ + Call_S ∈ {0,1}^{Q×Q} — call transitions for S ∈ N + Ret_S ∈ {0,1}^{Q×Q} — return transitions for S ∈ N + N_S ∈ {0,1}^{Q×Q} — internal nonterminal transitions + G_S ∈ {0,1}^{V×V} — derived graph edges for S + + M ∈ {0,1}^{Q×V} + P ∈ {0,1}^{Q×V} + + +------------------------------------------------------------ +Phase 1: Saturation (computation of reachable configurations) +------------------------------------------------------------ + +Initialization: + M[q₀, v₀] = 1 for all q₀ ∈ start_states + P[q₀, v₀] = 1 for all q₀ ∈ start_states + +while M != 0: + + Forward propagation: + M_new = 0 + for all a ∈ Σ: + M_new ← M_new ∪ (N_aᵀ · M · G_a) ∧ ~P + + For all S ∈ N: + M_new ← M_new ∪ (N_Sᵀ · M · G_S) ∧ ~P + M_new ← M_new ∪ (Call_Sᵀ · M) ∧ ~P + + Backward propagation (summary edge construction): + + + For all S ∈ N: + M_ret = Ret_Sᵀ · M + + if M_ret = 0: continue + + Restore paths to matching call positions by computing + the least fixed point of: + + X₀ = (Ret_S · M_ret) ∧ P + X_{i+1} = + ((⋁_{a∈Σ} N_a · X_i · G_aᵀ) + ∨ (⋁_{S1∈N} N_{S1} · X_i · G_{S1}ᵀ)) ∧ ~P + + After reaching fixed point X: + + X_call = (Call_S · X) ∧ P + + if X_call = 0: continue + + New summary edges: + + G_S_new = + (X_callᵀ · N_S · M_ret) ∧ ~G_S + + Update: + G_S ← G_S ∪ G_S_new + M_new = M_new ∪ M_ret + + M = M_new + P = P ∪ M + +Iterate until no new configurations are added to P. + +------------------------------------------------------------ +Phase 2: Reachability restricted to valid CF paths +------------------------------------------------------------ + +M = 0 + +Compute K ⊆ Q × V: + M[q₀, v₀] = 1 + K[q₀, v₀] = 1 for start configurations + +while M != 0: + + Forward propagation: + M_new = 0 + for all a ∈ Σ: + M_new ← M_new ∪ (N_aᵀ · M · G_a) + + For all S ∈ N: + M_new ← M_new ∪ (N_Sᵀ · M · G_S) + + M = M_new + P = P ∪ M + +until fixed point. + +------------------------------------------------------------ +Result +------------------------------------------------------------ + +Let F ∈ {0,1}^{1×Q} be final-state selector. + +Output: + result = F · K + +*/ + #include #include #include From 745f2157ab03d6dfadbc3380b869f6d08b8fbdfb Mon Sep 17 00:00:00 2001 From: Sterkachov Date: Wed, 8 Apr 2026 19:16:30 +0300 Subject: [PATCH 05/10] Implement second version of new algorithm --- .../algorithm/LAGraph_RSM_reachability.c | 141 ++++++------------ 1 file changed, 48 insertions(+), 93 deletions(-) diff --git a/experimental/algorithm/LAGraph_RSM_reachability.c b/experimental/algorithm/LAGraph_RSM_reachability.c index cf6faf523d..aa5336aa33 100644 --- a/experimental/algorithm/LAGraph_RSM_reachability.c +++ b/experimental/algorithm/LAGraph_RSM_reachability.c @@ -3,114 +3,69 @@ //------------------------------------------------------------------------------ /* -Let: +Initialize: + P := 0 + M := 0 + M[q_start, vs·V + vs] := 1 - Q = {0, …, Q-1} — set of RSM states - V = {0, …, V-1} — set of graph vertices - Σ — terminal alphabet - N — set of nonterminals + For each S: + Stack[S] := 0 // (q_ret, v_before, v_entry) + graph_nt[S] := 0 // (v_entry, v_final) -Boolean matrices: +while M != 0 do - G_a ∈ {0,1}^{V×V} — graph edges labeled by a ∈ Σ - N_a ∈ {0,1}^{Q×Q} — RSM transitions labeled by a ∈ Σ - Call_S ∈ {0,1}^{Q×Q} — call transitions for S ∈ N - Ret_S ∈ {0,1}^{Q×Q} — return transitions for S ∈ N - N_S ∈ {0,1}^{Q×Q} — internal nonterminal transitions - G_S ∈ {0,1}^{V×V} — derived graph edges for S + // -------------------------------------------------- + // 1. TERMINAL TRANSITIONS + // -------------------------------------------------- + M_term := 0 + for each label a do + for each q' do + M_term |= N_a^T @ M @ G_a - M ∈ {0,1}^{Q×V} - P ∈ {0,1}^{Q×V} + // -------------------------------------------------- + // 2. CALL TRANSITIONS + // -------------------------------------------------- + M_call := 0 + for each box S do + // select call states + mask_call := call[S]^T @ M // (q_start(S), v_start, v_entry) ------------------------------------------------------------- -Phase 1: Saturation (computation of reachable configurations) ------------------------------------------------------------- + if mask_call != 0: -Initialization: - M[q₀, v₀] = 1 for all q₀ ∈ start_states - P[q₀, v₀] = 1 for all q₀ ∈ start_states + // ---- push frames ---- + frame := rsm_nt[S]^T @ M // (q_ret, v_before, v_entry) + Stack[S] |= frame -while M != 0: + // ---- activate callee ---- + entry_vec := OR_rows(M[q_start]) // bcs exist only one q_start for each box + D := diag(entry_vec) + M_call[q_start] |= D - Forward propagation: - M_new = 0 - for all a ∈ Σ: - M_new ← M_new ∪ (N_aᵀ · M · G_a) ∧ ~P + // -------------------------------------------------- + // 3. RETURN TRANSITIONS + // -------------------------------------------------- + M_return := 0 + for each box S do - For all S ∈ N: - M_new ← M_new ∪ (N_Sᵀ · M · G_S) ∧ ~P - M_new ← M_new ∪ (Call_Sᵀ · M) ∧ ~P + M contains (q_final_S + other states, ..., ...) + Mask them with final_states and get (q_final_S, v_start_S, v_final_S) + derive (v_start_S, v_final_S) edges to graph_nt[S] + M_return |= Stack[S] @ graph_nt[S] // restore traviersion from point (q_return, v_before, +v_final_S) - Backward propagation (summary edge construction): + // -------------------------------------------------- + // 4. MERGE + // -------------------------------------------------- + M_new := (M_term | M_call | M_return) & ~P - For all S ∈ N: - M_ret = Ret_Sᵀ · M + P |= M_new + M := M_new - if M_ret = 0: continue +end while - Restore paths to matching call positions by computing - the least fixed point of: - - X₀ = (Ret_S · M_ret) ∧ P - X_{i+1} = - ((⋁_{a∈Σ} N_a · X_i · G_aᵀ) - ∨ (⋁_{S1∈N} N_{S1} · X_i · G_{S1}ᵀ)) ∧ ~P - - After reaching fixed point X: - - X_call = (Call_S · X) ∧ P - - if X_call = 0: continue - - New summary edges: - - G_S_new = - (X_callᵀ · N_S · M_ret) ∧ ~G_S - - Update: - G_S ← G_S ∪ G_S_new - M_new = M_new ∪ M_ret - - M = M_new - P = P ∪ M - -Iterate until no new configurations are added to P. - ------------------------------------------------------------- -Phase 2: Reachability restricted to valid CF paths ------------------------------------------------------------- - -M = 0 - -Compute K ⊆ Q × V: - M[q₀, v₀] = 1 - K[q₀, v₀] = 1 for start configurations - -while M != 0: - - Forward propagation: - M_new = 0 - for all a ∈ Σ: - M_new ← M_new ∪ (N_aᵀ · M · G_a) - - For all S ∈ N: - M_new ← M_new ∪ (N_Sᵀ · M · G_S) - - M = M_new - P = P ∪ M - -until fixed point. - ------------------------------------------------------------- -Result ------------------------------------------------------------- - -Let F ∈ {0,1}^{1×Q} be final-state selector. - -Output: - result = F · K +return P */ From f4581377a504bf3d04e2793720ddb3b003591e84 Mon Sep 17 00:00:00 2001 From: Sterkachov Date: Fri, 17 Apr 2026 05:20:15 +0300 Subject: [PATCH 06/10] CFPQ RSM algorithm v3 --- experimental/algorithm/LAGraph_CFPQ_RSM.c | 350 +++++++++ .../algorithm/LAGraph_RSM_reachability.c | 704 ----------------- experimental/test/test_CFPQ_RSM.c | 709 ++++++++++++++++++ experimental/test/test_RSM_reachability.c | 333 -------- include/LAGraphX.h | 45 +- 5 files changed, 1076 insertions(+), 1065 deletions(-) create mode 100644 experimental/algorithm/LAGraph_CFPQ_RSM.c delete mode 100644 experimental/algorithm/LAGraph_RSM_reachability.c create mode 100644 experimental/test/test_CFPQ_RSM.c delete mode 100644 experimental/test/test_RSM_reachability.c diff --git a/experimental/algorithm/LAGraph_CFPQ_RSM.c b/experimental/algorithm/LAGraph_CFPQ_RSM.c new file mode 100644 index 0000000000..5d107fc412 --- /dev/null +++ b/experimental/algorithm/LAGraph_CFPQ_RSM.c @@ -0,0 +1,350 @@ +#include "LAGraph.h" +#include "LG_internal.h" +#include + +#include +#include +#include + +#define LG_FREE_WORK \ + { \ + GrB_free(&M); \ + GrB_free(&P); \ + GrB_free(&M_term); \ + GrB_free(&M_nonterm); \ + GrB_free(&M_call); \ + GrB_free(&M_return); \ + GrB_free(&M_new); \ + GrB_free(&temp1); \ + GrB_free(&temp2); \ + GrB_free(&mask_call); \ + GrB_free(&frame); \ + GrB_free(&new_edges); \ + GrB_free(&diag_entry); \ + GrB_free(&outer_result); \ + GrB_free(&entry_vec); \ + if (Stack != NULL) \ + { \ + for (size_t _i = 0; _i < num_nonterminals; ++_i) \ + GrB_free(&Stack[_i]); \ + LAGraph_Free((void **)&Stack, GrB_NULL); \ + } \ + if (graph_nt != NULL) \ + { \ + for (size_t _i = 0; _i < num_nonterminals; ++_i) \ + GrB_free(&graph_nt[_i]); \ + LAGraph_Free((void **)&graph_nt, GrB_NULL); \ + } \ + } +#define LG_FREE_ALL \ + { \ + LG_FREE_WORK; \ + GrB_free(reachable); \ + } + +#define H_TRY(expr) \ + { \ + GrB_Info _hi = (expr); \ + if (_hi != GrB_SUCCESS) \ + return _hi; \ + } + +static GrB_Info s_mxm_chain(GrB_Matrix C, // output accumulator |Q|*|V*V| + GrB_Matrix A, // RSM matrix |Q|*|Q| + GrB_Matrix M, // frontier |Q|*|V*V| + GrB_Matrix B, // graph matrix |V|*|V| + GrB_Matrix temp1, // scratch |Q|*|V*V| + GrB_Matrix temp2, // scratch |Q|*|V*V| + GrB_Index Q, GrB_Index V) +{ + H_TRY(GrB_Matrix_clear(temp1)); + H_TRY(GrB_Matrix_clear(temp2)); + + // temp1 = A^T * M -> |Q|*|V*V| + H_TRY(GrB_mxm(temp1, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, A, M, GrB_DESC_T0)); + + // reshape: |Q| × |V*V| -> |Q*V|*|V| + H_TRY(GxB_Matrix_reshape(temp1, false, Q * V, V, GrB_NULL)); + H_TRY(GxB_Matrix_reshape(temp2, false, Q * V, V, GrB_NULL)); + + // temp2 = temp1 * B -> |Q*V|*|V| + H_TRY(GrB_mxm(temp2, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, temp1, B, GrB_NULL)); + + // reshape: |Q*V|*|V| -> |Q|*|V*V| + H_TRY(GxB_Matrix_reshape(temp1, false, Q, V * V, GrB_NULL)); + H_TRY(GxB_Matrix_reshape(temp2, false, Q, V * V, GrB_NULL)); + + // C |= temp2 + H_TRY(GrB_eWiseAdd(C, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, C, temp2, GrB_NULL)); + + return GrB_SUCCESS; +} +#undef H_TRY + +int LAGraph_CFPQ_RSM( + // output: + GrB_Vector *reachable, + + // RSM terminal input: + size_t num_terminals, GrB_Matrix *rsm_term, GrB_Matrix *graph_term, + + // RSM non-terminal input: + size_t num_nonterminals, GrB_Matrix *rsm_nonterm, GrB_Matrix *rsm_call, GrB_Vector *rsm_start, + GrB_Vector *rsm_final, + + size_t start_nonterm, GrB_Index start_vertex, + + GrB_Index Q, GrB_Index V, + + char *msg) +{ + LG_CLEAR_MSG; + + // working matrices / vectors + GrB_Matrix M = GrB_NULL; + GrB_Matrix P = GrB_NULL; + GrB_Matrix M_term = GrB_NULL; + GrB_Matrix M_nonterm = GrB_NULL; + GrB_Matrix M_call = GrB_NULL; + GrB_Matrix M_return = GrB_NULL; + GrB_Matrix M_new = GrB_NULL; + GrB_Matrix temp1 = GrB_NULL; + GrB_Matrix temp2 = GrB_NULL; + GrB_Matrix mask_call = GrB_NULL; + GrB_Matrix frame = GrB_NULL; + GrB_Matrix new_edges = GrB_NULL; + GrB_Matrix diag_entry = GrB_NULL; // add diagonal matrix for entry_vec + GrB_Matrix outer_result = GrB_NULL; // additional matrix + GrB_Vector entry_vec = GrB_NULL; + + GrB_Matrix *Stack = NULL; + GrB_Matrix *graph_nt = NULL; + + LG_ASSERT(reachable != NULL, GrB_NULL_POINTER); + LG_ASSERT(rsm_term != NULL, GrB_NULL_POINTER); + LG_ASSERT(rsm_nonterm != NULL, GrB_NULL_POINTER); + LG_ASSERT(rsm_call != NULL, GrB_NULL_POINTER); + LG_ASSERT(rsm_start != NULL, GrB_NULL_POINTER); + LG_ASSERT(rsm_final != NULL, GrB_NULL_POINTER); + LG_ASSERT(start_nonterm < num_nonterminals, GrB_INVALID_VALUE); + LG_ASSERT(start_vertex < V, GrB_INVALID_VALUE); + + *reachable = GrB_NULL; + + GrB_Index VV = V * V; + + // allocate per-nonterminal arrays + LG_TRY(LAGraph_Malloc((void **)&Stack, num_nonterminals, sizeof(GrB_Matrix), msg)); + LG_TRY(LAGraph_Malloc((void **)&graph_nt, num_nonterminals, sizeof(GrB_Matrix), msg)); + + for (size_t i = 0; i < num_nonterminals; ++i) + { + Stack[i] = GrB_NULL; + graph_nt[i] = GrB_NULL; + } + + for (size_t i = 0; i < num_nonterminals; ++i) + { + GRB_TRY(GrB_Matrix_new(&Stack[i], GrB_BOOL, Q, VV)); + GRB_TRY(GrB_Matrix_new(&graph_nt[i], GrB_BOOL, V, V)); + } + + // allocate working matrices + GRB_TRY(GrB_Matrix_new(&M, GrB_BOOL, Q, VV)); + GRB_TRY(GrB_Matrix_new(&P, GrB_BOOL, Q, VV)); + GRB_TRY(GrB_Matrix_new(&M_term, GrB_BOOL, Q, VV)); + GRB_TRY(GrB_Matrix_new(&M_nonterm, GrB_BOOL, Q, VV)); + GRB_TRY(GrB_Matrix_new(&M_call, GrB_BOOL, Q, VV)); + GRB_TRY(GrB_Matrix_new(&M_return, GrB_BOOL, Q, VV)); + GRB_TRY(GrB_Matrix_new(&M_new, GrB_BOOL, Q, VV)); + GRB_TRY(GrB_Matrix_new(&temp1, GrB_BOOL, Q, VV)); + GRB_TRY(GrB_Matrix_new(&temp2, GrB_BOOL, Q, VV)); + GRB_TRY(GrB_Matrix_new(&mask_call, GrB_BOOL, Q, VV)); + GRB_TRY(GrB_Matrix_new(&frame, GrB_BOOL, Q, VV)); + GRB_TRY(GrB_Matrix_new(&new_edges, GrB_BOOL, 1, VV)); + // GRB_TRY(GrB_Matrix_new(&diag_entry, GrB_BOOL, V, V)); + GRB_TRY(GrB_Matrix_new(&outer_result, GrB_BOOL, Q, VV)); + GRB_TRY(GrB_Vector_new(&entry_vec, GrB_BOOL, V)); + GRB_TRY(GrB_Vector_new(reachable, GrB_BOOL, V)); + + // seed initial frontier M and Stack + { + GrB_Index seed_col = start_vertex * V + start_vertex; + + GrB_Index nstart = 0; + GRB_TRY(GrB_Vector_nvals(&nstart, rsm_start[start_nonterm])); + + GrB_Index *start_states = NULL; + bool *start_vals = NULL; + LG_TRY(LAGraph_Malloc((void **)&start_states, nstart, sizeof(GrB_Index), msg)); + LG_TRY(LAGraph_Malloc((void **)&start_vals, nstart, sizeof(bool), msg)); + + GrB_Info info = GrB_Vector_extractTuples_BOOL(start_states, start_vals, &nstart, + rsm_start[start_nonterm]); + LAGraph_Free((void **)&start_vals, GrB_NULL); + + if (info != GrB_SUCCESS) + { + LAGraph_Free((void **)&start_states, GrB_NULL); + LG_FREE_ALL; + return info; + } + + for (GrB_Index k = 0; k < nstart; ++k) + { + GrB_Index q = start_states[k]; + GRB_TRY(GrB_Matrix_setElement_BOOL(M, true, q, seed_col)); + GRB_TRY(GrB_Matrix_setElement_BOOL(Stack[start_nonterm], true, q, seed_col)); + } + + LAGraph_Free((void **)&start_states, GrB_NULL); + } + + // main fixed-point loop + GrB_Index m_nvals = 0; + GRB_TRY(GrB_Matrix_nvals(&m_nvals, M)); + + while (m_nvals > 0) + { + + // PHASE 1: TERMINAL TRANSITIONS + GRB_TRY(GrB_Matrix_clear(M_term)); + for (size_t i = 0; i < num_terminals; ++i) + { + if (rsm_term[i] == GrB_NULL || graph_term[i] == GrB_NULL) + continue; + LG_TRY(s_mxm_chain(M_term, rsm_term[i], M, graph_term[i], temp1, temp2, Q, V)); + } + + // PHASE 2: NON-TERMINAL TRANSITIONS + GRB_TRY(GrB_Matrix_clear(M_nonterm)); + for (size_t i = 0; i < num_nonterminals; ++i) + { + if (rsm_nonterm[i] == GrB_NULL) + continue; + LG_TRY(s_mxm_chain(M_nonterm, rsm_nonterm[i], M, graph_nt[i], temp1, temp2, Q, V)); + } + + // PHASE 3: CALL TRANSITIONS + GRB_TRY(GrB_Matrix_clear(M_call)); + for (size_t i = 0; i < num_nonterminals; ++i) + { + if (rsm_call[i] == GrB_NULL || rsm_nonterm[i] == GrB_NULL) + continue; + + // mask_call = rsm_call[i]^T * M -> |Q|*|V*V| + GRB_TRY(GrB_Matrix_clear(mask_call)); + GRB_TRY(GrB_mxm(mask_call, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, rsm_call[i], + M, GrB_DESC_T0)); + + GrB_Index mc_nvals = 0; + GRB_TRY(GrB_Matrix_nvals(&mc_nvals, mask_call)); + if (mc_nvals == 0) + continue; + + // frame = rsm_nonterm[i]^T * M, accumulated into Stack[i] + GRB_TRY(GrB_Matrix_clear(frame)); + GRB_TRY(GrB_mxm(frame, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, rsm_nonterm[i], M, + GrB_DESC_T0)); + GRB_TRY(GrB_eWiseAdd(Stack[i], GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, Stack[i], + frame, GrB_NULL)); + + // entry_vec = column-OR of mask_call reshaped to |Q*V|*|V| + GRB_TRY(GrB_Matrix_clear(temp1)); + GRB_TRY(GrB_eWiseAdd(temp1, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, temp1, + mask_call, GrB_NULL)); + GRB_TRY(GxB_Matrix_reshape(temp1, false, Q * V, V, GrB_NULL)); + + GRB_TRY(GrB_Vector_clear(entry_vec)); + // column OR = row OR of transpose + GRB_TRY( + GrB_reduce(entry_vec, GrB_NULL, GrB_LOR, GrB_LOR_MONOID_BOOL, temp1, GrB_DESC_T0)); + + GRB_TRY(GxB_Matrix_reshape(temp1, false, Q, VV, GrB_NULL)); + + // Build diag(entry_vec) reshaped to 1*|V*V| + diag_entry = GrB_NULL; + GRB_TRY(GrB_Matrix_diag(&diag_entry, entry_vec, 0)); // |V|*|V| + GRB_TRY(GxB_Matrix_reshape(diag_entry, false, 1, VV, GrB_NULL)); + + GRB_TRY(GrB_Matrix_clear(outer_result)); + // |Q|*1 * 1*|V*V| -> |Q|*|V*V| + GRB_TRY(GrB_mxm(outer_result, GrB_NULL, GrB_NULL, GxB_LOR_LAND_BOOL, + (GrB_Matrix)rsm_start[i], diag_entry, GrB_NULL)); + + // M_call |= outer_result + GRB_TRY(GrB_eWiseAdd(M_call, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, M_call, + outer_result, GrB_NULL)); + + GRB_TRY(GrB_Matrix_free(&diag_entry)); + } + + // PHASE 4: RETURN TRANSITIONS + GRB_TRY(GrB_Matrix_clear(M_return)); + for (size_t i = 0; i < num_nonterminals; ++i) + { + if (rsm_final[i] == GrB_NULL) + continue; + + // 1*|Q| * |Q|*|V*V| -> 1*|V*V| + GRB_TRY(GrB_Matrix_clear(new_edges)); + GRB_TRY(GrB_mxm(new_edges, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, + (GrB_Matrix)rsm_final[i], M, GrB_DESC_T0)); + + GRB_TRY(GxB_Matrix_reshape(new_edges, false, V, V, GrB_NULL)); + + GRB_TRY(GrB_eWiseAdd(graph_nt[i], GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, + graph_nt[i], new_edges, GrB_NULL)); + + GRB_TRY(GrB_Matrix_clear(temp1)); + GRB_TRY(GrB_eWiseAdd(temp1, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, temp1, + Stack[i], GrB_NULL)); + GRB_TRY(GxB_Matrix_reshape(temp1, false, Q * V, V, GrB_NULL)); + + GRB_TRY(GrB_Matrix_clear(temp2)); + GRB_TRY(GxB_Matrix_reshape(temp2, false, Q * V, V, GrB_NULL)); + // |Q*V|*|V| * |V|*|V| -> |Q*V|*|V| + GRB_TRY(GrB_mxm(temp2, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, temp1, new_edges, + GrB_NULL)); + + GRB_TRY(GxB_Matrix_reshape(temp1, false, Q, VV, GrB_NULL)); + GRB_TRY(GxB_Matrix_reshape(temp2, false, Q, VV, GrB_NULL)); + + GRB_TRY(GrB_eWiseAdd(M_return, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, M_return, + temp2, GrB_NULL)); + GRB_TRY(GxB_Matrix_reshape(new_edges, false, 1, VV, GrB_NULL)); + } + + // PHASE 5: MERGE & MARK VISITED + // GrB_assign has big overhead + // M_new = (M_term | M_nonterm | M_call | M_return) & ~P + GRB_TRY(GrB_Matrix_clear(M_new)); + GRB_TRY(GrB_eWiseAdd(M_new, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, M_new, M_term, + GrB_NULL)); + GRB_TRY(GrB_eWiseAdd(M_new, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, M_new, M_nonterm, + GrB_NULL)); + GRB_TRY(GrB_eWiseAdd(M_new, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, M_new, M_call, + GrB_NULL)); + GRB_TRY(GrB_eWiseAdd(M_new, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, M_new, M_return, + GrB_NULL)); + + // M_new &= ~P + GRB_TRY(GrB_assign(M_new, P, GrB_NULL, M_new, GrB_ALL, Q, GrB_ALL, VV, GrB_DESC_RSC)); + + // P |= M_new + GRB_TRY(GrB_eWiseAdd(P, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, P, M_new, GrB_NULL)); + + GrB_Matrix swap = M; + M = M_new; + M_new = swap; + + GRB_TRY(GrB_Matrix_nvals(&m_nvals, M)); + } + + // GrB_DESC_T0 transposes graph_nt so Col_extract gives row start_vertex + GRB_TRY(GrB_Col_extract(*reachable, GrB_NULL, GrB_NULL, graph_nt[start_nonterm], GrB_ALL, V, + start_vertex, GrB_DESC_T0)); + + LG_FREE_WORK; + return GrB_SUCCESS; +} diff --git a/experimental/algorithm/LAGraph_RSM_reachability.c b/experimental/algorithm/LAGraph_RSM_reachability.c deleted file mode 100644 index aa5336aa33..0000000000 --- a/experimental/algorithm/LAGraph_RSM_reachability.c +++ /dev/null @@ -1,704 +0,0 @@ -//------------------------------------------------------------------------------ -// LAGraph_RSM_reachability.c -//------------------------------------------------------------------------------ - -/* -Initialize: - P := 0 - M := 0 - M[q_start, vs·V + vs] := 1 - - For each S: - Stack[S] := 0 // (q_ret, v_before, v_entry) - graph_nt[S] := 0 // (v_entry, v_final) - -while M != 0 do - - // -------------------------------------------------- - // 1. TERMINAL TRANSITIONS - // -------------------------------------------------- - M_term := 0 - for each label a do - for each q' do - M_term |= N_a^T @ M @ G_a - - // -------------------------------------------------- - // 2. CALL TRANSITIONS - // -------------------------------------------------- - M_call := 0 - for each box S do - - // select call states - mask_call := call[S]^T @ M // (q_start(S), v_start, v_entry) - - if mask_call != 0: - - // ---- push frames ---- - frame := rsm_nt[S]^T @ M // (q_ret, v_before, v_entry) - Stack[S] |= frame - - // ---- activate callee ---- - entry_vec := OR_rows(M[q_start]) // bcs exist only one q_start for each box - D := diag(entry_vec) - M_call[q_start] |= D - - // -------------------------------------------------- - // 3. RETURN TRANSITIONS - // -------------------------------------------------- - M_return := 0 - for each box S do - - M contains (q_final_S + other states, ..., ...) - Mask them with final_states and get (q_final_S, v_start_S, v_final_S) - derive (v_start_S, v_final_S) edges to graph_nt[S] - M_return |= Stack[S] @ graph_nt[S] // restore traviersion from point (q_return, v_before, -v_final_S) - - - // -------------------------------------------------- - // 4. MERGE - // -------------------------------------------------- - M_new := (M_term | M_call | M_return) & ~P - - P |= M_new - M := M_new - -end while - -return P - -*/ - -#include -#include -#include -#include -#include - -#include "LAGraph.h" -#include "LG_internal.h" - -#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, - - // 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; - - // 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_Index ng = 0; // # nodes in the graph - GrB_Index nr = 0; // # states in the RSM - GrB_Index states = ns; // # pairs in current frontier - - GrB_Index rows = 0, cols = 0, vals = 0; - - GrB_Matrix *A_a = GrB_NULL; // for G_a - GrB_Matrix *B_a = GrB_NULL; // for N_a - - 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 - - LG_ASSERT(reachable != GrB_NULL, GrB_NULL_POINTER); - - 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); - - (*reachable) = GrB_NULL; - - //-------------------------------------------------------------------------- - // 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)) - - 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 < 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_a[i] == GrB_NULL) - A_a[i] = GrB_NULL; - else - A_a[i] = G_a[i]->A; - } - LG_TRY(LAGraph_Malloc((void **)&B_a, num_terminals, sizeof(GrB_Matrix), msg)) - for (size_t i = 0; i < num_terminals; ++i) - { - if (N_a[i] == GrB_NULL) - B_a[i] = GrB_NULL; - else - B_a[i] = N_a[i]->A; - } - LG_TRY(LAGraph_Malloc((void **)&B_S, num_nonterminals, sizeof(GrB_Matrix), msg)) - for (size_t i = 0; i < num_nonterminals; ++i) - { - 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; - } - - //-------------------------------------------------------------------------- - // determine ng (graph size) and nr (RSM state count) - //-------------------------------------------------------------------------- - - for (size_t i = 0; i < num_terminals; ++i) - { - if (A_a[i] == GrB_NULL) - continue; - - 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; - - GRB_TRY(GrB_Matrix_nrows(&nr, B_a[i])); - break; - } - - //-------------------------------------------------------------------------- - // 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])); - - 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 (B_a[i] == GrB_NULL) - continue; - - GRB_TRY(GrB_Matrix_nrows(&rows, B_a[i])); - GRB_TRY(GrB_Matrix_ncols(&cols, B_a[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_S[i] == GrB_NULL) - continue; - - GRB_TRY(GrB_Matrix_nrows(&rows, B_S[i])); - GRB_TRY(GrB_Matrix_ncols(&cols, B_S[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 (B_call[i] == GrB_NULL) - continue; - - GRB_TRY(GrB_Matrix_nrows(&rows, B_call[i])); - GRB_TRY(GrB_Matrix_ncols(&cols, B_call[i])); - - 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") - } - for (size_t i = 0; i < num_nonterminals; ++i) - { - if (B_ret[i] == GrB_NULL) - continue; - - 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, - "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) - { - 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) - { - 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"); - } - - //-------------------------------------------------------------------------- - // 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)); - - // final_reducer[QF] = true - 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(&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; - - // Main loop - while (states != 0) - { - iteration++; - - GrB_Matrix old_frontier = frontier; - frontier = next_frontier; - next_frontier = old_frontier; - GRB_TRY(GrB_Matrix_clear(next_frontier)); - - //---------------------------------------------------------------------- - // FORWARD PHASE - //---------------------------------------------------------------------- - - // 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)); - - // 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)); - - // 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 & ~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_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)); - } - - // 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)); - - // 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)); - - // 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; - } - } - - 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; -} diff --git a/experimental/test/test_CFPQ_RSM.c b/experimental/test/test_CFPQ_RSM.c new file mode 100644 index 0000000000..8dd9953bd6 --- /dev/null +++ b/experimental/test/test_CFPQ_RSM.c @@ -0,0 +1,709 @@ +#include "LAGraph.h" +#include "LAGraphX.h" +#include "LAGraph_test.h" +#include + +char msg[LAGRAPH_MSG_LEN]; + +//============================================================================== +// test setup / teardown +//============================================================================== + +static void setup(void) +{ + OK(LAGraph_Init(msg)); +} + +static void teardown(void) +{ + OK(LAGraph_Finalize(msg)); +} + +static GrB_Info build_call_matrix(GrB_Matrix *call, GrB_Matrix rsm_nt, GrB_Vector rsm_start, + GrB_Index Q) +{ + GrB_Vector q_call_mask = GrB_NULL; + + GRB_TRY(GrB_Matrix_new(call, GrB_BOOL, Q, Q)); + GRB_TRY(GrB_Vector_new(&q_call_mask, GrB_BOOL, Q)); + + // [q_call, q_ret] -> [q_call] + GRB_TRY(GrB_reduce(q_call_mask, GrB_NULL, GrB_NULL, GrB_LOR_MONOID_BOOL, rsm_nt, GrB_NULL)); + // |Q|*1 * 1*|Q| -> |Q|*|Q| + GRB_TRY(GrB_mxm(*call, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, (GrB_Matrix)q_call_mask, + (GrB_Matrix)rsm_start, GrB_DESC_T1)); + + GRB_TRY(GrB_Vector_free(&q_call_mask)); + + return GrB_SUCCESS; +} + +static void check_reachable(GrB_Vector result, GrB_Index *expected, GrB_Index n_expected, + GrB_Index V, const char *test_name) +{ + GrB_Index nvals = 0; + OK(GrB_Vector_nvals(&nvals, result)); + + bool *got = (bool *)calloc(V, sizeof(bool)); + TEST_CHECK(got != NULL); + + GrB_Index *idx = GrB_NULL; + bool *val = GrB_NULL; + + if (nvals > 0) + { + idx = (GrB_Index *)malloc(nvals * sizeof(GrB_Index)); + val = (bool *)malloc(nvals * sizeof(bool)); + if (!TEST_CHECK(idx != NULL && val != NULL)) + { + free(got); + free(idx); + free(val); + return; + } + OK(GrB_Vector_extractTuples_BOOL(idx, val, &nvals, result)); + for (GrB_Index k = 0; k < nvals; ++k) + if (val[k]) + got[idx[k]] = true; + } + + // every expected vertex must be reachable + for (GrB_Index k = 0; k < n_expected; ++k) + { + TEST_CHECK(got[expected[k]]); + TEST_MSG("%s: vertex %llu should be reachable", test_name, (unsigned long long)expected[k]); + } + + GrB_Index actual_count = 0; + for (GrB_Index v = 0; v < V; ++v) + if (got[v]) + actual_count++; + + TEST_CHECK(actual_count == n_expected); + TEST_MSG("%s: expected %llu reachable vertices, got %llu", test_name, + (unsigned long long)n_expected, (unsigned long long)actual_count); + + free(got); + free(idx); + free(val); +} + +//============================================================================== +// RSM builders +//============================================================================== +// Grammar 1: S -> a S b | a b +static void rsm_aSb_ab(GrB_Matrix *term, GrB_Matrix *nt, GrB_Matrix *call, GrB_Vector *start, + GrB_Vector *final) +{ + GrB_Index Q = 4; + OK(GrB_Matrix_new(&term[0], GrB_BOOL, Q, Q)); + OK(GrB_Matrix_new(&term[1], GrB_BOOL, Q, Q)); + OK(GrB_Matrix_new(&nt[0], GrB_BOOL, Q, Q)); + OK(GrB_Vector_new(&start[0], GrB_BOOL, Q)); + OK(GrB_Vector_new(&final[0], GrB_BOOL, Q)); + + OK(GrB_Matrix_setElement_BOOL(term[0], true, 0, 1)); // 0 -a-> 1 + OK(GrB_Matrix_setElement_BOOL(term[1], true, 1, 3)); // 1 -b-> 3 + OK(GrB_Matrix_setElement_BOOL(term[1], true, 2, 3)); // 2 -b-> 3 + OK(GrB_Matrix_setElement_BOOL(nt[0], true, 1, 2)); // 1 -S-> 2 + OK(GrB_Vector_setElement_BOOL(start[0], true, 0)); + OK(GrB_Vector_setElement_BOOL(final[0], true, 3)); + + OK(build_call_matrix(call, *nt, *start, Q)); +} + +// Grammar 2: S -> a S b | c +static void rsm_aSb_c(GrB_Matrix *term, GrB_Matrix *nt, GrB_Matrix *call, GrB_Vector *start, + GrB_Vector *final) +{ + GrB_Index Q = 4; + OK(GrB_Matrix_new(&term[0], GrB_BOOL, Q, Q)); + OK(GrB_Matrix_new(&term[1], GrB_BOOL, Q, Q)); + OK(GrB_Matrix_new(&term[2], GrB_BOOL, Q, Q)); + OK(GrB_Matrix_new(&nt[0], GrB_BOOL, Q, Q)); + OK(GrB_Vector_new(&start[0], GrB_BOOL, Q)); + OK(GrB_Vector_new(&final[0], GrB_BOOL, Q)); + + OK(GrB_Matrix_setElement_BOOL(term[0], true, 0, 1)); // 0 -a-> 1 + OK(GrB_Matrix_setElement_BOOL(term[2], true, 0, 3)); // 1 -c-> 3 + OK(GrB_Matrix_setElement_BOOL(term[1], true, 2, 3)); // 2 -b-> 3 + OK(GrB_Matrix_setElement_BOOL(nt[0], true, 1, 2)); // 1 -S-> 2 + OK(GrB_Vector_setElement_BOOL(start[0], true, 0)); + OK(GrB_Vector_setElement_BOOL(final[0], true, 3)); + + OK(build_call_matrix(call, *nt, *start, Q)); +} + +// Grammar 3: S -> X | Y | S X | S Y +// X -> a S b | a b +// Y -> c S d | c d +static void rsm_complex(GrB_Matrix *term, GrB_Matrix *nt, GrB_Matrix *call, GrB_Vector *start, + GrB_Vector *final) +{ + GrB_Index Q = 11; + for (int t = 0; t < 4; ++t) + OK(GrB_Matrix_new(&term[t], GrB_BOOL, Q, Q)); + for (int n = 0; n < 3; ++n) + { + OK(GrB_Matrix_new(&nt[n], GrB_BOOL, Q, Q)); + OK(GrB_Vector_new(&start[n], GrB_BOOL, Q)); + OK(GrB_Vector_new(&final[n], GrB_BOOL, Q)); + } + + // --- S Component (States 0-2) --- + // S -> X | Y | S X | S Y + OK(GrB_Matrix_setElement_BOOL(nt[1], true, 0, 2)); // 0 -X-> 2 + OK(GrB_Matrix_setElement_BOOL(nt[2], true, 0, 2)); // 0 -Y-> 2 + OK(GrB_Matrix_setElement_BOOL(nt[0], true, 0, 1)); // 0 -S-> 1 + OK(GrB_Matrix_setElement_BOOL(nt[1], true, 1, 2)); // 1 -X-> 2 + OK(GrB_Matrix_setElement_BOOL(nt[2], true, 1, 2)); // 1 -Y-> 2 + OK(GrB_Vector_setElement_BOOL(start[0], true, 0)); + OK(GrB_Vector_setElement_BOOL(final[0], true, 2)); + + // --- X Component (States 3-6) --- + OK(GrB_Matrix_setElement_BOOL(term[0], true, 3, 4)); // 3 -a-> 4 + OK(GrB_Matrix_setElement_BOOL(term[1], true, 4, 6)); // 4 -b-> 6 + OK(GrB_Matrix_setElement_BOOL(term[1], true, 5, 6)); // 5 -b-> 6 + OK(GrB_Matrix_setElement_BOOL(nt[0], true, 4, 5)); // 4 -S-> 5 + OK(GrB_Vector_setElement_BOOL(start[1], true, 3)); + OK(GrB_Vector_setElement_BOOL(final[1], true, 6)); + + // --- Y Component (States 7-10) + OK(GrB_Matrix_setElement_BOOL(term[2], true, 7, 8)); // 7 -c-> 8 + OK(GrB_Matrix_setElement_BOOL(term[3], true, 8, 10)); // 8 -d-> 10 + OK(GrB_Matrix_setElement_BOOL(term[3], true, 9, 10)); // 9 -d-> 10 + OK(GrB_Matrix_setElement_BOOL(nt[0], true, 8, 9)); // 8 -S-> 9 + OK(GrB_Vector_setElement_BOOL(start[2], true, 7)); + OK(GrB_Vector_setElement_BOOL(final[2], true, 10)); + + for (int n = 0; n < 3; n++) + OK(build_call_matrix(&call[n], nt[n], start[n], Q)); +} + +//============================================================================== +// test cases +//============================================================================== + +//------------------------------------------------------------------------------ +// TC1: S -> a S b | a b +// Graph: 0-a->1; 1-a->2; 2-a->0; 0-b->3; 3-b->0 +// Source: 1 Expected: {0, 3} +//------------------------------------------------------------------------------ +void test_TC1_cyclic_ab(void) +{ + setup(); + + GrB_Matrix term[2], nt[1], call[1]; + GrB_Vector start[1], final[1]; + + term[0] = term[1] = GrB_NULL; + nt[0] = GrB_NULL; + call[0] = GrB_NULL; + start[0] = final[0] = GrB_NULL; + + rsm_aSb_ab(term, nt, call, start, final); + + GrB_Index V = 4; + GrB_Matrix gterm[2]; + gterm[0] = gterm[1] = GrB_NULL; + OK(GrB_Matrix_new(>erm[0], GrB_BOOL, V, V)); + OK(GrB_Matrix_new(>erm[1], GrB_BOOL, V, V)); + + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 1, 2)); + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 2, 0)); + OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 0, 3)); + OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 3, 0)); + + GrB_Vector result = GrB_NULL; + + GrB_Info info = + LAGraph_CFPQ_RSM(&result, 2, term, gterm, 1, nt, call, start, final, 0, 1, 4, V, msg); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {0, 3}; + check_reachable(result, expected, 2, V, "TC1"); + } + + GrB_free(&result); + for (int i = 0; i < 2; ++i) + { + GrB_free(&term[i]); + GrB_free(>erm[i]); + } + GrB_free(&nt[0]); + GrB_free(&call[0]); + GrB_free(&start[0]); + GrB_free(&final[0]); + + OK(info); + + teardown(); +} + +//------------------------------------------------------------------------------ +// TC2: S -> a S b | a b +// Graph: 0-a->0; 0-a->1; 1-b->1 +// Source: 0 Expected: {1} +//------------------------------------------------------------------------------ +void test_TC2_self_loops(void) +{ + setup(); + + GrB_Matrix term[2], nt[1], call[1]; + GrB_Vector start[1], final[1]; + + term[0] = term[1] = GrB_NULL; + nt[0] = GrB_NULL; + call[0] = GrB_NULL; + start[0] = final[0] = GrB_NULL; + + rsm_aSb_ab(term, nt, call, start, final); + + GrB_Index V = 2; + GrB_Matrix gterm[2]; + gterm[0] = gterm[1] = GrB_NULL; + OK(GrB_Matrix_new(>erm[0], GrB_BOOL, V, V)); + OK(GrB_Matrix_new(>erm[1], GrB_BOOL, V, V)); + + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 0)); + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 1, 1)); + + GrB_Vector result = GrB_NULL; + + GrB_Info info = + LAGraph_CFPQ_RSM(&result, 2, term, gterm, 1, nt, call, start, final, 0, 0, 4, V, msg); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {1}; + check_reachable(result, expected, 1, V, "TC2"); + } + + GrB_free(&result); + for (int i = 0; i < 2; ++i) + { + GrB_free(&term[i]); + GrB_free(>erm[i]); + } + GrB_free(&nt[0]); + GrB_free(&call[0]); + GrB_free(&start[0]); + GrB_free(&final[0]); + + OK(info); + + teardown(); +} + +//------------------------------------------------------------------------------ +// TC3: S -> a S b | a b +// Graph: 0-a->0; 0-b->0 +// Source: 0 Expected: {0} +//------------------------------------------------------------------------------ +void test_TC3_single_node(void) +{ + setup(); + + GrB_Matrix term[2], nt[1], call[1]; + GrB_Vector start[1], final[1]; + + term[0] = term[1] = GrB_NULL; + nt[0] = GrB_NULL; + call[0] = GrB_NULL; + start[0] = final[0] = GrB_NULL; + + rsm_aSb_ab(term, nt, call, start, final); + + GrB_Index V = 1; + GrB_Matrix gterm[2]; + gterm[0] = gterm[1] = GrB_NULL; + OK(GrB_Matrix_new(>erm[0], GrB_BOOL, V, V)); + OK(GrB_Matrix_new(>erm[1], GrB_BOOL, V, V)); + + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 0)); + OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 0, 0)); + + GrB_Vector result = GrB_NULL; + + GrB_Info info = + LAGraph_CFPQ_RSM(&result, 2, term, gterm, 1, nt, call, start, final, 0, 0, 4, V, msg); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {0}; + check_reachable(result, expected, 1, V, "TC3"); + } + + GrB_free(&result); + for (int i = 0; i < 2; ++i) + { + GrB_free(&term[i]); + GrB_free(>erm[i]); + } + GrB_free(&nt[0]); + GrB_free(&call[0]); + GrB_free(&start[0]); + GrB_free(&final[0]); + + OK(info); + + teardown(); +} + +//------------------------------------------------------------------------------ +// TC4: S -> a S b | c +// Graph: 0-a->0; 0-c->1; 1-b->1 +// Source: 0 Expected: {1} +//------------------------------------------------------------------------------ +void test_TC4_aSb_c_loops(void) +{ + setup(); + + GrB_Matrix term[3], nt[1], call[1]; + GrB_Vector start[1], final[1]; + + term[0] = term[1] = term[2] = GrB_NULL; + nt[0] = GrB_NULL; + call[0] = GrB_NULL; + start[0] = final[0] = GrB_NULL; + + rsm_aSb_c(term, nt, call, start, final); + + GrB_Index V = 2; + GrB_Matrix gterm[3]; + gterm[0] = gterm[1] = gterm[2] = GrB_NULL; + for (int i = 0; i < 3; ++i) + OK(GrB_Matrix_new(>erm[i], GrB_BOOL, V, V)); + + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 0)); + OK(GrB_Matrix_setElement_BOOL(gterm[2], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 1, 1)); + + GrB_Vector result = GrB_NULL; + + GrB_Info info = + LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, call, start, final, 0, 0, 4, V, msg); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {1}; + check_reachable(result, expected, 1, V, "TC4"); + } + + GrB_free(&result); + for (int i = 0; i < 3; ++i) + { + GrB_free(&term[i]); + GrB_free(>erm[i]); + } + GrB_free(&nt[0]); + GrB_free(&call[0]); + GrB_free(&start[0]); + GrB_free(&final[0]); + + OK(info); + + teardown(); +} + +//------------------------------------------------------------------------------ +// TC5: S -> a S b | c +// Graph: 0-a->1; 1-a->0; 0-c->2; 2-b->3; 3-b->2 +// Source: 1 Expected: {3} +//------------------------------------------------------------------------------ +void test_TC5_aSb_c_cycle1(void) +{ + setup(); + + GrB_Matrix term[3], nt[1], call[1]; + GrB_Vector start[1], final[1]; + + term[0] = term[1] = term[2] = GrB_NULL; + nt[0] = GrB_NULL; + call[0] = GrB_NULL; + start[0] = final[0] = GrB_NULL; + + rsm_aSb_c(term, nt, call, start, final); + + GrB_Index V = 4; + GrB_Matrix gterm[3]; + gterm[0] = gterm[1] = gterm[2] = GrB_NULL; + for (int i = 0; i < 3; ++i) + OK(GrB_Matrix_new(>erm[i], GrB_BOOL, V, V)); + + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 1, 0)); + OK(GrB_Matrix_setElement_BOOL(gterm[2], true, 0, 2)); + OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 2, 3)); + OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 3, 2)); + + GrB_Vector result = GrB_NULL; + + GrB_Info info = + LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, call, start, final, 0, 1, 4, V, msg); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {3}; + check_reachable(result, expected, 1, V, "TC5"); + } + + GrB_free(&result); + for (int i = 0; i < 3; ++i) + { + GrB_free(&term[i]); + GrB_free(>erm[i]); + } + GrB_free(&nt[0]); + GrB_free(&call[0]); + GrB_free(&start[0]); + GrB_free(&final[0]); + + OK(info); + + teardown(); +} +//------------------------------------------------------------------------------ +// TC6: S -> a S b | c +// Graph: 0-a->1; 1-a->0; 0-c->2; 2-b->2 +// Source: 1 Expected: {2} +//------------------------------------------------------------------------------ +void test_TC6_aSb_c_cycle2(void) +{ + setup(); + + GrB_Matrix term[3], nt[1], call[1]; + GrB_Vector start[1], final[1]; + + term[0] = term[1] = term[2] = GrB_NULL; + nt[0] = GrB_NULL; + call[0] = GrB_NULL; + start[0] = final[0] = GrB_NULL; + + rsm_aSb_c(term, nt, call, start, final); + + GrB_Index V = 3; + GrB_Matrix gterm[3]; + gterm[0] = gterm[1] = gterm[2] = GrB_NULL; + for (int i = 0; i < 3; ++i) + OK(GrB_Matrix_new(>erm[i], GrB_BOOL, V, V)); + + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 1, 0)); + OK(GrB_Matrix_setElement_BOOL(gterm[2], true, 0, 2)); + OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 2, 2)); + + GrB_Vector result = GrB_NULL; + + GrB_Info info = + LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, call, start, final, 0, 1, 4, V, msg); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {2}; + check_reachable(result, expected, 1, V, "TC6"); + } + + GrB_free(&result); + for (int i = 0; i < 3; ++i) + { + GrB_free(&term[i]); + GrB_free(>erm[i]); + } + GrB_free(&nt[0]); + GrB_free(&call[0]); + GrB_free(&start[0]); + GrB_free(&final[0]); + + OK(info); + + teardown(); +} + +// Graph: 0-a->1; 1-a->0; 0-c->2; 2-b->3; 3-b->4; 4-b->2 +static void graph_tc7_8(GrB_Matrix gterm[3]) +{ + GrB_Index V = 5; + for (int i = 0; i < 3; ++i) + OK(GrB_Matrix_new(>erm[i], GrB_BOOL, V, V)); + + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 1, 0)); + OK(GrB_Matrix_setElement_BOOL(gterm[2], true, 0, 2)); + OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 2, 3)); + OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 3, 4)); + OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 4, 2)); +} + +//------------------------------------------------------------------------------ +// TC7: S -> a S b | c +// same graph as TC7/8 above +// Source: 0 Expected: {2, 3, 4} +//------------------------------------------------------------------------------ +void test_TC7_aSb_c_cycle3(void) +{ + setup(); + + GrB_Matrix term[3], nt[1], call[1]; + GrB_Vector start[1], final[1]; + + term[0] = term[1] = term[2] = GrB_NULL; + nt[0] = GrB_NULL; + call[0] = GrB_NULL; + start[0] = final[0] = GrB_NULL; + + rsm_aSb_c(term, nt, call, start, final); + + GrB_Index V = 5; + GrB_Matrix gterm[3]; + gterm[0] = gterm[1] = gterm[2] = GrB_NULL; + graph_tc7_8(gterm); + + GrB_Vector result = GrB_NULL; + + GrB_Info info = + LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, call, start, final, 0, 0, 4, V, msg); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {2, 3, 4}; + check_reachable(result, expected, 3, V, "TC7"); + } + + GrB_free(&result); + for (int i = 0; i < 3; ++i) + { + GrB_free(&term[i]); + GrB_free(>erm[i]); + } + GrB_free(&nt[0]); + GrB_free(&call[0]); + GrB_free(&start[0]); + GrB_free(&final[0]); + + OK(info); + + teardown(); +} + +//------------------------------------------------------------------------------ +// TC8: S -> a S b | c +// same graph, different source vertex +// Source: 1 Expected: {2, 3, 4} +//------------------------------------------------------------------------------ +void test_TC8_aSb_c_cycle4(void) +{ + setup(); + + GrB_Matrix term[3], nt[1], call[1]; + GrB_Vector start[1], final[1]; + + term[0] = term[1] = term[2] = GrB_NULL; + nt[0] = GrB_NULL; + call[0] = GrB_NULL; + start[0] = final[0] = GrB_NULL; + + rsm_aSb_c(term, nt, call, start, final); + + GrB_Index V = 5; + GrB_Matrix gterm[3]; + gterm[0] = gterm[1] = gterm[2] = GrB_NULL; + graph_tc7_8(gterm); + + GrB_Vector result = GrB_NULL; + + GrB_Info info = + LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, call, start, final, 0, 1, 4, V, msg); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {2, 3, 4}; + check_reachable(result, expected, 3, V, "TC8"); + } + + GrB_free(&result); + for (int i = 0; i < 3; ++i) + { + GrB_free(&term[i]); + GrB_free(>erm[i]); + } + GrB_free(&nt[0]); + GrB_free(&call[0]); + GrB_free(&start[0]); + GrB_free(&final[0]); + + OK(info); + + teardown(); +} + +void test_TC9_mult_recursive_cycles(void) +{ + setup(); + + GrB_Matrix term[4], nt[3], call[3]; + GrB_Vector start[3], final[3]; + + term[0] = term[1] = term[2] = term[3] = GrB_NULL; + nt[0] = nt[1] = nt[2] = GrB_NULL; + call[0] = call[1] = call[2] = GrB_NULL; + start[0] = start[1] = start[2] = GrB_NULL; + final[0] = final[1] = final[2] = GrB_NULL; + + rsm_complex(term, nt, call, start, final); + + GrB_Index V = 6; + GrB_Matrix gterm[4]; + gterm[0] = gterm[1] = gterm[2] = gterm[3] = GrB_NULL; + for (int i = 0; i < 4; ++i) + OK(GrB_Matrix_new(>erm[i], GrB_BOOL, V, V)); + + // a-cycle (0,1,2) + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 1, 2)); + OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 2, 0)); + // b-cycle (2,3) + OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 2, 3)); + OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 3, 2)); + // c-cycle (3,4,5) + OK(GrB_Matrix_setElement_BOOL(gterm[2], true, 3, 4)); + OK(GrB_Matrix_setElement_BOOL(gterm[2], true, 4, 5)); + OK(GrB_Matrix_setElement_BOOL(gterm[2], true, 5, 3)); + // d-cycle + OK(GrB_Matrix_setElement_BOOL(gterm[3], true, 5, 2)); + OK(GrB_Matrix_setElement_BOOL(gterm[3], true, 2, 5)); + + GrB_Vector result = GrB_NULL; + + GrB_Info info = + LAGraph_CFPQ_RSM(&result, 4, term, gterm, 3, nt, call, start, final, 0, 0, 11, V, msg); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {2, 3, 5}; + check_reachable(result, expected, 3, V, "TC9"); + } + + GrB_free(&result); + for (int i = 0; i < 4; ++i) + { + GrB_free(&term[i]); + GrB_free(>erm[i]); + } + for (int i = 0; i < 3; ++i) + { + GrB_free(&nt[i]); + GrB_free(&call[i]); + GrB_free(&start[i]); + GrB_free(&final[i]); + } + + OK(info); + + teardown(); +} + +TEST_LIST = {{"TC1_cyclic_ab", test_TC1_cyclic_ab}, + {"TC2_self_loops", test_TC2_self_loops}, + {"TC3_single_node", test_TC3_single_node}, + {"TC4_aSb_c_loops", test_TC4_aSb_c_loops}, + {"TC5_aSb_c_cycle1", test_TC5_aSb_c_cycle1}, + {"TC4_aSb_c_loops", test_TC4_aSb_c_loops}, + {"TC7_aSb_c_cycle3", test_TC7_aSb_c_cycle3}, + {"TC8_aSb_c_cycle4", test_TC8_aSb_c_cycle4}, + {"TC9_mult_recursive_cycles", test_TC9_mult_recursive_cycles}, + {NULL, NULL}}; diff --git a/experimental/test/test_RSM_reachability.c b/experimental/test/test_RSM_reachability.c deleted file mode 100644 index bf10a389de..0000000000 --- a/experimental/test/test_RSM_reachability.c +++ /dev/null @@ -1,333 +0,0 @@ -#include -#include -#include -#include -#include - -#include "LAGraph.h" - -#include -#include -#include -#include - -#define LEN 512 -#define MAX_LABELS 3 -#define MAX_RESULTS 2000000 - -#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 msg[LAGRAPH_MSG_LEN]; - -static void setup(void) -{ - LAGraph_Init(msg); -} - -static void teardown(void) -{ - LAGraph_Finalize(msg); -} - -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 (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; -} - -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); -} - -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); - - for (GrB_Index v = 0; v < V; v++) - { - bool val = false; - GrB_Info info = GrB_Vector_extractElement_BOOL(&val, result, v); - - if (expected[v]) - { - 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); - } - } -} - -//============================================================================== -// TEST CASES -//============================================================================== - -void test_based(void) -{ - setup(); - - const GrB_Index Q = 4, V = 6; - - // 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}; - - // RSM: S: q0 -a-> q1 -S-> q2 -b-> q3 - - GrB_Index Na_I[] = {0}, Na_J[] = {1}; - GrB_Index Nb_I[] = {1, 2}, Nb_J[] = {3, 3}; // q1->q3, q2->q3 - - 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 - - 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); - - 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[] = {3}, 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(); -} - -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(); -} - -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(); -} - -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 = {{"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 daf5d5727d..ea81c2589e 100644 --- a/include/LAGraphX.h +++ b/include/LAGraphX.h @@ -856,34 +856,6 @@ int LAGraph_RegularPathQuery // nodes reachable from the starting by the char *msg // LAGraph output message ) ; -//**************************************************************************** -LAGRAPHX_PUBLIC -int LAGraph_RSM_reachability( - // output: - 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 int LAGraph_VertexCentrality_Triangle // vertex triangle-centrality @@ -1125,6 +1097,23 @@ GrB_Info LAGraph_CFL_reachability char *msg // Message string for error reporting. ) ; +int LAGraph_CFPQ_RSM( + // output: + GrB_Vector *reachable, + + // RSM terminal input: + size_t num_terminals, GrB_Matrix *rsm_term, GrB_Matrix *graph_term, + + // RSM non-terminal input: + size_t num_nonterminals, GrB_Matrix *rsm_nonterm, GrB_Matrix *rsm_call, GrB_Vector *rsm_start, + GrB_Vector *rsm_final, + + size_t start_nonterm, GrB_Index start_vertex, + + GrB_Index Q, GrB_Index V, + + char *msg); + //------------------------------------------------------------------------------ // a simple example of an algorithm //------------------------------------------------------------------------------ From 530be6c9b0775ea4ee361fe10bd9ced36f51da09 Mon Sep 17 00:00:00 2001 From: Sterkachov Date: Fri, 17 Apr 2026 22:48:24 +0300 Subject: [PATCH 07/10] refactor: built rsm_call matrix internally instead of accepting as parameter --- experimental/algorithm/LAGraph_CFPQ_RSM.c | 34 ++++++- experimental/test/test_CFPQ_RSM.c | 111 ++++++---------------- include/LAGraphX.h | 2 +- 3 files changed, 61 insertions(+), 86 deletions(-) diff --git a/experimental/algorithm/LAGraph_CFPQ_RSM.c b/experimental/algorithm/LAGraph_CFPQ_RSM.c index 5d107fc412..c04adc4050 100644 --- a/experimental/algorithm/LAGraph_CFPQ_RSM.c +++ b/experimental/algorithm/LAGraph_CFPQ_RSM.c @@ -35,6 +35,12 @@ GrB_free(&graph_nt[_i]); \ LAGraph_Free((void **)&graph_nt, GrB_NULL); \ } \ + if (rsm_call != NULL) \ + { \ + for (size_t _i = 0; _i < num_nonterminals; ++_i) \ + GrB_free(&rsm_call[_i]); \ + LAGraph_Free((void **)&rsm_call, GrB_NULL); \ + } \ } #define LG_FREE_ALL \ { \ @@ -79,6 +85,26 @@ static GrB_Info s_mxm_chain(GrB_Matrix C, // output accumulator |Q|*|V*V| return GrB_SUCCESS; } + +static GrB_Info build_call_matrix(GrB_Matrix *call, GrB_Matrix rsm_nt, GrB_Vector rsm_start, + GrB_Index Q) +{ + GrB_Vector q_call_mask = GrB_NULL; + + H_TRY(GrB_Matrix_new(call, GrB_BOOL, Q, Q)); + H_TRY(GrB_Vector_new(&q_call_mask, GrB_BOOL, Q)); + + // [q_call, q_ret] -> [q_call] + H_TRY(GrB_reduce(q_call_mask, GrB_NULL, GrB_NULL, GrB_LOR_MONOID_BOOL, rsm_nt, GrB_NULL)); + // |Q|*1 * 1*|Q| -> |Q|*|Q| + H_TRY(GrB_mxm(*call, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, (GrB_Matrix)q_call_mask, + (GrB_Matrix)rsm_start, GrB_DESC_T1)); + + H_TRY(GrB_Vector_free(&q_call_mask)); + + return GrB_SUCCESS; +} + #undef H_TRY int LAGraph_CFPQ_RSM( @@ -89,8 +115,7 @@ int LAGraph_CFPQ_RSM( size_t num_terminals, GrB_Matrix *rsm_term, GrB_Matrix *graph_term, // RSM non-terminal input: - size_t num_nonterminals, GrB_Matrix *rsm_nonterm, GrB_Matrix *rsm_call, GrB_Vector *rsm_start, - GrB_Vector *rsm_final, + size_t num_nonterminals, GrB_Matrix *rsm_nonterm, GrB_Vector *rsm_start, GrB_Vector *rsm_final, size_t start_nonterm, GrB_Index start_vertex, @@ -119,11 +144,11 @@ int LAGraph_CFPQ_RSM( GrB_Matrix *Stack = NULL; GrB_Matrix *graph_nt = NULL; + GrB_Matrix *rsm_call = NULL; LG_ASSERT(reachable != NULL, GrB_NULL_POINTER); LG_ASSERT(rsm_term != NULL, GrB_NULL_POINTER); LG_ASSERT(rsm_nonterm != NULL, GrB_NULL_POINTER); - LG_ASSERT(rsm_call != NULL, GrB_NULL_POINTER); LG_ASSERT(rsm_start != NULL, GrB_NULL_POINTER); LG_ASSERT(rsm_final != NULL, GrB_NULL_POINTER); LG_ASSERT(start_nonterm < num_nonterminals, GrB_INVALID_VALUE); @@ -136,17 +161,20 @@ int LAGraph_CFPQ_RSM( // allocate per-nonterminal arrays LG_TRY(LAGraph_Malloc((void **)&Stack, num_nonterminals, sizeof(GrB_Matrix), msg)); LG_TRY(LAGraph_Malloc((void **)&graph_nt, num_nonterminals, sizeof(GrB_Matrix), msg)); + LG_TRY(LAGraph_Malloc((void **)&rsm_call, num_nonterminals, sizeof(GrB_Matrix), msg)); for (size_t i = 0; i < num_nonterminals; ++i) { Stack[i] = GrB_NULL; graph_nt[i] = GrB_NULL; + rsm_call[i] = GrB_NULL; } for (size_t i = 0; i < num_nonterminals; ++i) { GRB_TRY(GrB_Matrix_new(&Stack[i], GrB_BOOL, Q, VV)); GRB_TRY(GrB_Matrix_new(&graph_nt[i], GrB_BOOL, V, V)); + GRB_TRY(build_call_matrix(&rsm_call[i], rsm_nonterm[i], rsm_start[i], Q)); } // allocate working matrices diff --git a/experimental/test/test_CFPQ_RSM.c b/experimental/test/test_CFPQ_RSM.c index 8dd9953bd6..763a0aad70 100644 --- a/experimental/test/test_CFPQ_RSM.c +++ b/experimental/test/test_CFPQ_RSM.c @@ -19,25 +19,6 @@ static void teardown(void) OK(LAGraph_Finalize(msg)); } -static GrB_Info build_call_matrix(GrB_Matrix *call, GrB_Matrix rsm_nt, GrB_Vector rsm_start, - GrB_Index Q) -{ - GrB_Vector q_call_mask = GrB_NULL; - - GRB_TRY(GrB_Matrix_new(call, GrB_BOOL, Q, Q)); - GRB_TRY(GrB_Vector_new(&q_call_mask, GrB_BOOL, Q)); - - // [q_call, q_ret] -> [q_call] - GRB_TRY(GrB_reduce(q_call_mask, GrB_NULL, GrB_NULL, GrB_LOR_MONOID_BOOL, rsm_nt, GrB_NULL)); - // |Q|*1 * 1*|Q| -> |Q|*|Q| - GRB_TRY(GrB_mxm(*call, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, (GrB_Matrix)q_call_mask, - (GrB_Matrix)rsm_start, GrB_DESC_T1)); - - GRB_TRY(GrB_Vector_free(&q_call_mask)); - - return GrB_SUCCESS; -} - static void check_reachable(GrB_Vector result, GrB_Index *expected, GrB_Index n_expected, GrB_Index V, const char *test_name) { @@ -92,8 +73,7 @@ static void check_reachable(GrB_Vector result, GrB_Index *expected, GrB_Index n_ // RSM builders //============================================================================== // Grammar 1: S -> a S b | a b -static void rsm_aSb_ab(GrB_Matrix *term, GrB_Matrix *nt, GrB_Matrix *call, GrB_Vector *start, - GrB_Vector *final) +static void rsm_aSb_ab(GrB_Matrix *term, GrB_Matrix *nt, GrB_Vector *start, GrB_Vector *final) { GrB_Index Q = 4; OK(GrB_Matrix_new(&term[0], GrB_BOOL, Q, Q)); @@ -108,13 +88,10 @@ static void rsm_aSb_ab(GrB_Matrix *term, GrB_Matrix *nt, GrB_Matrix *call, GrB_V OK(GrB_Matrix_setElement_BOOL(nt[0], true, 1, 2)); // 1 -S-> 2 OK(GrB_Vector_setElement_BOOL(start[0], true, 0)); OK(GrB_Vector_setElement_BOOL(final[0], true, 3)); - - OK(build_call_matrix(call, *nt, *start, Q)); } // Grammar 2: S -> a S b | c -static void rsm_aSb_c(GrB_Matrix *term, GrB_Matrix *nt, GrB_Matrix *call, GrB_Vector *start, - GrB_Vector *final) +static void rsm_aSb_c(GrB_Matrix *term, GrB_Matrix *nt, GrB_Vector *start, GrB_Vector *final) { GrB_Index Q = 4; OK(GrB_Matrix_new(&term[0], GrB_BOOL, Q, Q)); @@ -130,15 +107,12 @@ static void rsm_aSb_c(GrB_Matrix *term, GrB_Matrix *nt, GrB_Matrix *call, GrB_Ve OK(GrB_Matrix_setElement_BOOL(nt[0], true, 1, 2)); // 1 -S-> 2 OK(GrB_Vector_setElement_BOOL(start[0], true, 0)); OK(GrB_Vector_setElement_BOOL(final[0], true, 3)); - - OK(build_call_matrix(call, *nt, *start, Q)); } // Grammar 3: S -> X | Y | S X | S Y // X -> a S b | a b // Y -> c S d | c d -static void rsm_complex(GrB_Matrix *term, GrB_Matrix *nt, GrB_Matrix *call, GrB_Vector *start, - GrB_Vector *final) +static void rsm_complex(GrB_Matrix *term, GrB_Matrix *nt, GrB_Vector *start, GrB_Vector *final) { GrB_Index Q = 11; for (int t = 0; t < 4; ++t) @@ -175,9 +149,6 @@ static void rsm_complex(GrB_Matrix *term, GrB_Matrix *nt, GrB_Matrix *call, GrB_ OK(GrB_Matrix_setElement_BOOL(nt[0], true, 8, 9)); // 8 -S-> 9 OK(GrB_Vector_setElement_BOOL(start[2], true, 7)); OK(GrB_Vector_setElement_BOOL(final[2], true, 10)); - - for (int n = 0; n < 3; n++) - OK(build_call_matrix(&call[n], nt[n], start[n], Q)); } //============================================================================== @@ -193,15 +164,14 @@ void test_TC1_cyclic_ab(void) { setup(); - GrB_Matrix term[2], nt[1], call[1]; + GrB_Matrix term[2], nt[1]; GrB_Vector start[1], final[1]; term[0] = term[1] = GrB_NULL; nt[0] = GrB_NULL; - call[0] = GrB_NULL; start[0] = final[0] = GrB_NULL; - rsm_aSb_ab(term, nt, call, start, final); + rsm_aSb_ab(term, nt, start, final); GrB_Index V = 4; GrB_Matrix gterm[2]; @@ -217,8 +187,7 @@ void test_TC1_cyclic_ab(void) GrB_Vector result = GrB_NULL; - GrB_Info info = - LAGraph_CFPQ_RSM(&result, 2, term, gterm, 1, nt, call, start, final, 0, 1, 4, V, msg); + GrB_Info info = LAGraph_CFPQ_RSM(&result, 2, term, gterm, 1, nt, start, final, 0, 1, 4, V, msg); if (info == GrB_SUCCESS) { GrB_Index expected[] = {0, 3}; @@ -232,7 +201,6 @@ void test_TC1_cyclic_ab(void) GrB_free(>erm[i]); } GrB_free(&nt[0]); - GrB_free(&call[0]); GrB_free(&start[0]); GrB_free(&final[0]); @@ -250,15 +218,14 @@ void test_TC2_self_loops(void) { setup(); - GrB_Matrix term[2], nt[1], call[1]; + GrB_Matrix term[2], nt[1]; GrB_Vector start[1], final[1]; term[0] = term[1] = GrB_NULL; nt[0] = GrB_NULL; - call[0] = GrB_NULL; start[0] = final[0] = GrB_NULL; - rsm_aSb_ab(term, nt, call, start, final); + rsm_aSb_ab(term, nt, start, final); GrB_Index V = 2; GrB_Matrix gterm[2]; @@ -272,8 +239,7 @@ void test_TC2_self_loops(void) GrB_Vector result = GrB_NULL; - GrB_Info info = - LAGraph_CFPQ_RSM(&result, 2, term, gterm, 1, nt, call, start, final, 0, 0, 4, V, msg); + GrB_Info info = LAGraph_CFPQ_RSM(&result, 2, term, gterm, 1, nt, start, final, 0, 0, 4, V, msg); if (info == GrB_SUCCESS) { GrB_Index expected[] = {1}; @@ -287,7 +253,6 @@ void test_TC2_self_loops(void) GrB_free(>erm[i]); } GrB_free(&nt[0]); - GrB_free(&call[0]); GrB_free(&start[0]); GrB_free(&final[0]); @@ -305,15 +270,14 @@ void test_TC3_single_node(void) { setup(); - GrB_Matrix term[2], nt[1], call[1]; + GrB_Matrix term[2], nt[1]; GrB_Vector start[1], final[1]; term[0] = term[1] = GrB_NULL; nt[0] = GrB_NULL; - call[0] = GrB_NULL; start[0] = final[0] = GrB_NULL; - rsm_aSb_ab(term, nt, call, start, final); + rsm_aSb_ab(term, nt, start, final); GrB_Index V = 1; GrB_Matrix gterm[2]; @@ -326,8 +290,7 @@ void test_TC3_single_node(void) GrB_Vector result = GrB_NULL; - GrB_Info info = - LAGraph_CFPQ_RSM(&result, 2, term, gterm, 1, nt, call, start, final, 0, 0, 4, V, msg); + GrB_Info info = LAGraph_CFPQ_RSM(&result, 2, term, gterm, 1, nt, start, final, 0, 0, 4, V, msg); if (info == GrB_SUCCESS) { GrB_Index expected[] = {0}; @@ -341,7 +304,6 @@ void test_TC3_single_node(void) GrB_free(>erm[i]); } GrB_free(&nt[0]); - GrB_free(&call[0]); GrB_free(&start[0]); GrB_free(&final[0]); @@ -359,15 +321,14 @@ void test_TC4_aSb_c_loops(void) { setup(); - GrB_Matrix term[3], nt[1], call[1]; + GrB_Matrix term[3], nt[1]; GrB_Vector start[1], final[1]; term[0] = term[1] = term[2] = GrB_NULL; nt[0] = GrB_NULL; - call[0] = GrB_NULL; start[0] = final[0] = GrB_NULL; - rsm_aSb_c(term, nt, call, start, final); + rsm_aSb_c(term, nt, start, final); GrB_Index V = 2; GrB_Matrix gterm[3]; @@ -381,8 +342,7 @@ void test_TC4_aSb_c_loops(void) GrB_Vector result = GrB_NULL; - GrB_Info info = - LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, call, start, final, 0, 0, 4, V, msg); + GrB_Info info = LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, start, final, 0, 0, 4, V, msg); if (info == GrB_SUCCESS) { GrB_Index expected[] = {1}; @@ -396,7 +356,6 @@ void test_TC4_aSb_c_loops(void) GrB_free(>erm[i]); } GrB_free(&nt[0]); - GrB_free(&call[0]); GrB_free(&start[0]); GrB_free(&final[0]); @@ -414,15 +373,14 @@ void test_TC5_aSb_c_cycle1(void) { setup(); - GrB_Matrix term[3], nt[1], call[1]; + GrB_Matrix term[3], nt[1]; GrB_Vector start[1], final[1]; term[0] = term[1] = term[2] = GrB_NULL; nt[0] = GrB_NULL; - call[0] = GrB_NULL; start[0] = final[0] = GrB_NULL; - rsm_aSb_c(term, nt, call, start, final); + rsm_aSb_c(term, nt, start, final); GrB_Index V = 4; GrB_Matrix gterm[3]; @@ -438,8 +396,7 @@ void test_TC5_aSb_c_cycle1(void) GrB_Vector result = GrB_NULL; - GrB_Info info = - LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, call, start, final, 0, 1, 4, V, msg); + GrB_Info info = LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, start, final, 0, 1, 4, V, msg); if (info == GrB_SUCCESS) { GrB_Index expected[] = {3}; @@ -453,7 +410,6 @@ void test_TC5_aSb_c_cycle1(void) GrB_free(>erm[i]); } GrB_free(&nt[0]); - GrB_free(&call[0]); GrB_free(&start[0]); GrB_free(&final[0]); @@ -470,15 +426,14 @@ void test_TC6_aSb_c_cycle2(void) { setup(); - GrB_Matrix term[3], nt[1], call[1]; + GrB_Matrix term[3], nt[1]; GrB_Vector start[1], final[1]; term[0] = term[1] = term[2] = GrB_NULL; nt[0] = GrB_NULL; - call[0] = GrB_NULL; start[0] = final[0] = GrB_NULL; - rsm_aSb_c(term, nt, call, start, final); + rsm_aSb_c(term, nt, start, final); GrB_Index V = 3; GrB_Matrix gterm[3]; @@ -493,8 +448,7 @@ void test_TC6_aSb_c_cycle2(void) GrB_Vector result = GrB_NULL; - GrB_Info info = - LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, call, start, final, 0, 1, 4, V, msg); + GrB_Info info = LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, start, final, 0, 1, 4, V, msg); if (info == GrB_SUCCESS) { GrB_Index expected[] = {2}; @@ -508,7 +462,6 @@ void test_TC6_aSb_c_cycle2(void) GrB_free(>erm[i]); } GrB_free(&nt[0]); - GrB_free(&call[0]); GrB_free(&start[0]); GrB_free(&final[0]); @@ -541,15 +494,14 @@ void test_TC7_aSb_c_cycle3(void) { setup(); - GrB_Matrix term[3], nt[1], call[1]; + GrB_Matrix term[3], nt[1]; GrB_Vector start[1], final[1]; term[0] = term[1] = term[2] = GrB_NULL; nt[0] = GrB_NULL; - call[0] = GrB_NULL; start[0] = final[0] = GrB_NULL; - rsm_aSb_c(term, nt, call, start, final); + rsm_aSb_c(term, nt, start, final); GrB_Index V = 5; GrB_Matrix gterm[3]; @@ -558,8 +510,7 @@ void test_TC7_aSb_c_cycle3(void) GrB_Vector result = GrB_NULL; - GrB_Info info = - LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, call, start, final, 0, 0, 4, V, msg); + GrB_Info info = LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, start, final, 0, 0, 4, V, msg); if (info == GrB_SUCCESS) { GrB_Index expected[] = {2, 3, 4}; @@ -573,7 +524,6 @@ void test_TC7_aSb_c_cycle3(void) GrB_free(>erm[i]); } GrB_free(&nt[0]); - GrB_free(&call[0]); GrB_free(&start[0]); GrB_free(&final[0]); @@ -599,7 +549,7 @@ void test_TC8_aSb_c_cycle4(void) call[0] = GrB_NULL; start[0] = final[0] = GrB_NULL; - rsm_aSb_c(term, nt, call, start, final); + rsm_aSb_c(term, nt, start, final); GrB_Index V = 5; GrB_Matrix gterm[3]; @@ -608,8 +558,7 @@ void test_TC8_aSb_c_cycle4(void) GrB_Vector result = GrB_NULL; - GrB_Info info = - LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, call, start, final, 0, 1, 4, V, msg); + GrB_Info info = LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, start, final, 0, 1, 4, V, msg); if (info == GrB_SUCCESS) { GrB_Index expected[] = {2, 3, 4}; @@ -636,16 +585,15 @@ void test_TC9_mult_recursive_cycles(void) { setup(); - GrB_Matrix term[4], nt[3], call[3]; + GrB_Matrix term[4], nt[3]; GrB_Vector start[3], final[3]; term[0] = term[1] = term[2] = term[3] = GrB_NULL; nt[0] = nt[1] = nt[2] = GrB_NULL; - call[0] = call[1] = call[2] = GrB_NULL; start[0] = start[1] = start[2] = GrB_NULL; final[0] = final[1] = final[2] = GrB_NULL; - rsm_complex(term, nt, call, start, final); + rsm_complex(term, nt, start, final); GrB_Index V = 6; GrB_Matrix gterm[4]; @@ -671,7 +619,7 @@ void test_TC9_mult_recursive_cycles(void) GrB_Vector result = GrB_NULL; GrB_Info info = - LAGraph_CFPQ_RSM(&result, 4, term, gterm, 3, nt, call, start, final, 0, 0, 11, V, msg); + LAGraph_CFPQ_RSM(&result, 4, term, gterm, 3, nt, start, final, 0, 0, 11, V, msg); if (info == GrB_SUCCESS) { GrB_Index expected[] = {2, 3, 5}; @@ -687,7 +635,6 @@ void test_TC9_mult_recursive_cycles(void) for (int i = 0; i < 3; ++i) { GrB_free(&nt[i]); - GrB_free(&call[i]); GrB_free(&start[i]); GrB_free(&final[i]); } diff --git a/include/LAGraphX.h b/include/LAGraphX.h index ea81c2589e..2a5a7249a0 100644 --- a/include/LAGraphX.h +++ b/include/LAGraphX.h @@ -1105,7 +1105,7 @@ int LAGraph_CFPQ_RSM( size_t num_terminals, GrB_Matrix *rsm_term, GrB_Matrix *graph_term, // RSM non-terminal input: - size_t num_nonterminals, GrB_Matrix *rsm_nonterm, GrB_Matrix *rsm_call, GrB_Vector *rsm_start, + size_t num_nonterminals, GrB_Matrix *rsm_nonterm, GrB_Vector *rsm_start, GrB_Vector *rsm_final, size_t start_nonterm, GrB_Index start_vertex, From fca24cacdeef5efa01cb7bc857f9c1d328cd0802 Mon Sep 17 00:00:00 2001 From: Sterkachov Date: Sun, 19 Apr 2026 15:32:01 +0300 Subject: [PATCH 08/10] refactor: consolidate RSM parameters into struct and modularize tests --- experimental/algorithm/LAGraph_CFPQ_RSM.c | 107 ++-- experimental/test/test_CFPQ_RSM.c | 698 ++++++++++------------ include/LAGraphX.h | 30 +- 3 files changed, 386 insertions(+), 449 deletions(-) diff --git a/experimental/algorithm/LAGraph_CFPQ_RSM.c b/experimental/algorithm/LAGraph_CFPQ_RSM.c index c04adc4050..35ac23be58 100644 --- a/experimental/algorithm/LAGraph_CFPQ_RSM.c +++ b/experimental/algorithm/LAGraph_CFPQ_RSM.c @@ -25,21 +25,27 @@ GrB_free(&entry_vec); \ if (Stack != NULL) \ { \ - for (size_t _i = 0; _i < num_nonterminals; ++_i) \ + for (size_t _i = 0; _i < num_nonterm; ++_i) \ GrB_free(&Stack[_i]); \ - LAGraph_Free((void **)&Stack, GrB_NULL); \ + LAGraph_Free((void **)&Stack, msg); \ } \ if (graph_nt != NULL) \ { \ - for (size_t _i = 0; _i < num_nonterminals; ++_i) \ + for (size_t _i = 0; _i < num_nonterm; ++_i) \ GrB_free(&graph_nt[_i]); \ - LAGraph_Free((void **)&graph_nt, GrB_NULL); \ + LAGraph_Free((void **)&graph_nt, msg); \ + } \ + if (rsm_start != NULL) \ + { \ + for (size_t _i = 0; _i < num_nonterm; ++_i) \ + GrB_free(&rsm_start[_i]); \ + LAGraph_Free((void **)&rsm_start, msg); \ } \ if (rsm_call != NULL) \ { \ - for (size_t _i = 0; _i < num_nonterminals; ++_i) \ + for (size_t _i = 0; _i < num_nonterm; ++_i) \ GrB_free(&rsm_call[_i]); \ - LAGraph_Free((void **)&rsm_call, GrB_NULL); \ + LAGraph_Free((void **)&rsm_call, msg); \ } \ } #define LG_FREE_ALL \ @@ -89,42 +95,51 @@ static GrB_Info s_mxm_chain(GrB_Matrix C, // output accumulator |Q|*|V*V| static GrB_Info build_call_matrix(GrB_Matrix *call, GrB_Matrix rsm_nt, GrB_Vector rsm_start, GrB_Index Q) { - GrB_Vector q_call_mask = GrB_NULL; + GrB_Vector call_mask = GrB_NULL; H_TRY(GrB_Matrix_new(call, GrB_BOOL, Q, Q)); - H_TRY(GrB_Vector_new(&q_call_mask, GrB_BOOL, Q)); + H_TRY(GrB_Vector_new(&call_mask, GrB_BOOL, Q)); // [q_call, q_ret] -> [q_call] - H_TRY(GrB_reduce(q_call_mask, GrB_NULL, GrB_NULL, GrB_LOR_MONOID_BOOL, rsm_nt, GrB_NULL)); + H_TRY(GrB_reduce(call_mask, GrB_NULL, GrB_NULL, GrB_LOR_MONOID_BOOL, rsm_nt, GrB_NULL)); // |Q|*1 * 1*|Q| -> |Q|*|Q| - H_TRY(GrB_mxm(*call, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, (GrB_Matrix)q_call_mask, + H_TRY(GrB_mxm(*call, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, (GrB_Matrix)call_mask, (GrB_Matrix)rsm_start, GrB_DESC_T1)); - H_TRY(GrB_Vector_free(&q_call_mask)); + H_TRY(GrB_Vector_free(&call_mask)); return GrB_SUCCESS; } #undef H_TRY - -int LAGraph_CFPQ_RSM( - // output: - GrB_Vector *reachable, - - // RSM terminal input: - size_t num_terminals, GrB_Matrix *rsm_term, GrB_Matrix *graph_term, - - // RSM non-terminal input: - size_t num_nonterminals, GrB_Matrix *rsm_nonterm, GrB_Vector *rsm_start, GrB_Vector *rsm_final, - - size_t start_nonterm, GrB_Index start_vertex, - - GrB_Index Q, GrB_Index V, - - char *msg) +typedef struct +{ + GrB_Index state_count; // Q + GrB_Index terminal_count; // num_terminals + GrB_Index nonterminal_count; // num_nonterminals + GrB_Index start_nonterminal; // start_nonterm + GrB_Matrix *terminal_matrices; // rsm_term + GrB_Matrix *nonterminal_matrices; // rsm_nonterm + GrB_Index *start_states; // rsm_start + GrB_Vector *final_states; // rsm_final +} RSM; + +int LAGraph_CFPQ_RSM(GrB_Vector *reachable, RSM *rsm, GrB_Matrix *graph_term, + GrB_Index start_vertex, GrB_Index V, char *msg) { LG_CLEAR_MSG; + // shortcuts + GrB_Index VV = V * V; + GrB_Index Q = rsm->state_count; + GrB_Index num_term = rsm->terminal_count; + GrB_Index num_nonterm = rsm->nonterminal_count; + GrB_Index start_nonterm = rsm->start_nonterminal; + GrB_Matrix *rsm_term = rsm->terminal_matrices; + GrB_Matrix *rsm_nonterm = rsm->nonterminal_matrices; + GrB_Index *start_states = rsm->start_states; + GrB_Vector *final_states = rsm->final_states; + // working matrices / vectors GrB_Matrix M = GrB_NULL; GrB_Matrix P = GrB_NULL; @@ -145,35 +160,40 @@ int LAGraph_CFPQ_RSM( GrB_Matrix *Stack = NULL; GrB_Matrix *graph_nt = NULL; GrB_Matrix *rsm_call = NULL; + GrB_Vector *rsm_start = NULL; // rsm_start[S][i] = true, iff start_states[S] = i LG_ASSERT(reachable != NULL, GrB_NULL_POINTER); LG_ASSERT(rsm_term != NULL, GrB_NULL_POINTER); LG_ASSERT(rsm_nonterm != NULL, GrB_NULL_POINTER); - LG_ASSERT(rsm_start != NULL, GrB_NULL_POINTER); - LG_ASSERT(rsm_final != NULL, GrB_NULL_POINTER); - LG_ASSERT(start_nonterm < num_nonterminals, GrB_INVALID_VALUE); + LG_ASSERT(start_states != NULL, GrB_NULL_POINTER); + LG_ASSERT(final_states != NULL, GrB_NULL_POINTER); + LG_ASSERT(start_nonterm < num_nonterm, GrB_INVALID_VALUE); LG_ASSERT(start_vertex < V, GrB_INVALID_VALUE); *reachable = GrB_NULL; - GrB_Index VV = V * V; - // allocate per-nonterminal arrays - LG_TRY(LAGraph_Malloc((void **)&Stack, num_nonterminals, sizeof(GrB_Matrix), msg)); - LG_TRY(LAGraph_Malloc((void **)&graph_nt, num_nonterminals, sizeof(GrB_Matrix), msg)); - LG_TRY(LAGraph_Malloc((void **)&rsm_call, num_nonterminals, sizeof(GrB_Matrix), msg)); + LG_TRY(LAGraph_Malloc((void **)&Stack, rsm->nonterminal_count, sizeof(GrB_Matrix), msg)); + LG_TRY(LAGraph_Malloc((void **)&graph_nt, rsm->nonterminal_count, sizeof(GrB_Matrix), msg)); + LG_TRY(LAGraph_Malloc((void **)&rsm_start, rsm->nonterminal_count, sizeof(GrB_Vector), msg)); + LG_TRY(LAGraph_Malloc((void **)&rsm_call, rsm->nonterminal_count, sizeof(GrB_Matrix), msg)); - for (size_t i = 0; i < num_nonterminals; ++i) + for (size_t i = 0; i < num_nonterm; ++i) { Stack[i] = GrB_NULL; graph_nt[i] = GrB_NULL; + rsm_start[i] = GrB_NULL; rsm_call[i] = GrB_NULL; } - for (size_t i = 0; i < num_nonterminals; ++i) + for (size_t i = 0; i < num_nonterm; ++i) { GRB_TRY(GrB_Matrix_new(&Stack[i], GrB_BOOL, Q, VV)); GRB_TRY(GrB_Matrix_new(&graph_nt[i], GrB_BOOL, V, V)); + + GRB_TRY(GrB_Vector_new(&rsm_start[i], GrB_BOOL, Q)); + GrB_Vector_setElement_BOOL(rsm_start[i], true, start_states[i]); + GRB_TRY(build_call_matrix(&rsm_call[i], rsm_nonterm[i], rsm_start[i], Q)); } @@ -231,13 +251,12 @@ int LAGraph_CFPQ_RSM( // main fixed-point loop GrB_Index m_nvals = 0; GRB_TRY(GrB_Matrix_nvals(&m_nvals, M)); - while (m_nvals > 0) { // PHASE 1: TERMINAL TRANSITIONS GRB_TRY(GrB_Matrix_clear(M_term)); - for (size_t i = 0; i < num_terminals; ++i) + for (size_t i = 0; i < num_term; ++i) { if (rsm_term[i] == GrB_NULL || graph_term[i] == GrB_NULL) continue; @@ -246,7 +265,7 @@ int LAGraph_CFPQ_RSM( // PHASE 2: NON-TERMINAL TRANSITIONS GRB_TRY(GrB_Matrix_clear(M_nonterm)); - for (size_t i = 0; i < num_nonterminals; ++i) + for (size_t i = 0; i < num_nonterm; ++i) { if (rsm_nonterm[i] == GrB_NULL) continue; @@ -255,7 +274,7 @@ int LAGraph_CFPQ_RSM( // PHASE 3: CALL TRANSITIONS GRB_TRY(GrB_Matrix_clear(M_call)); - for (size_t i = 0; i < num_nonterminals; ++i) + for (size_t i = 0; i < num_nonterm; ++i) { if (rsm_call[i] == GrB_NULL || rsm_nonterm[i] == GrB_NULL) continue; @@ -309,15 +328,15 @@ int LAGraph_CFPQ_RSM( // PHASE 4: RETURN TRANSITIONS GRB_TRY(GrB_Matrix_clear(M_return)); - for (size_t i = 0; i < num_nonterminals; ++i) + for (size_t i = 0; i < num_nonterm; ++i) { - if (rsm_final[i] == GrB_NULL) + if (final_states[i] == GrB_NULL) continue; // 1*|Q| * |Q|*|V*V| -> 1*|V*V| GRB_TRY(GrB_Matrix_clear(new_edges)); GRB_TRY(GrB_mxm(new_edges, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, - (GrB_Matrix)rsm_final[i], M, GrB_DESC_T0)); + (GrB_Matrix)final_states[i], M, GrB_DESC_T0)); GRB_TRY(GxB_Matrix_reshape(new_edges, false, V, V, GrB_NULL)); diff --git a/experimental/test/test_CFPQ_RSM.c b/experimental/test/test_CFPQ_RSM.c index 763a0aad70..344c87ccc5 100644 --- a/experimental/test/test_CFPQ_RSM.c +++ b/experimental/test/test_CFPQ_RSM.c @@ -2,25 +2,101 @@ #include "LAGraphX.h" #include "LAGraph_test.h" #include +#include char msg[LAGRAPH_MSG_LEN]; +static RSM *rsm = NULL; +static size_t V = 0; +static size_t n_graph_term = 0; +static GrB_Matrix *graph_term = NULL; +static GrB_Vector result = GrB_NULL; + //============================================================================== // test setup / teardown //============================================================================== static void setup(void) { - OK(LAGraph_Init(msg)); + LAGraph_Init(msg); } static void teardown(void) { - OK(LAGraph_Finalize(msg)); + LAGraph_Finalize(msg); +} + +static RSM *alloc_rsm(GrB_Index state_count, GrB_Index term_count, GrB_Index nonterm_count, + GrB_Index start_nonterm) +{ + RSM *r = NULL; + LAGraph_Malloc((void **)&r, 1, sizeof(RSM), msg); + if (r == NULL) + return NULL; + + r->state_count = state_count; + r->terminal_count = term_count; + r->nonterminal_count = nonterm_count; + r->start_nonterminal = start_nonterm; + + LAGraph_Calloc((void **)&r->terminal_matrices, term_count, sizeof(GrB_Matrix), msg); + LAGraph_Calloc((void **)&r->nonterminal_matrices, nonterm_count, sizeof(GrB_Matrix), msg); + LAGraph_Calloc((void **)&r->start_states, nonterm_count, sizeof(GrB_Index), msg); + LAGraph_Calloc((void **)&r->final_states, nonterm_count, sizeof(GrB_Vector), msg); + + return r; +} + +static void free_rsm(void) +{ + if (rsm == NULL) + return; + + if (rsm->terminal_matrices != NULL) + { + for (size_t i = 0; i < rsm->terminal_count; ++i) + GrB_free(&rsm->terminal_matrices[i]); + LAGraph_Free((void **)&rsm->terminal_matrices, msg); + } + if (rsm->nonterminal_matrices != NULL) + { + for (size_t i = 0; i < rsm->nonterminal_count; ++i) + GrB_free(&rsm->nonterminal_matrices[i]); + LAGraph_Free((void **)&rsm->nonterminal_matrices, msg); + } + + if (rsm->nonterminal_matrices != NULL) + { + for (size_t i = 0; i < rsm->nonterminal_count; ++i) + { + GrB_free(&rsm->final_states[i]); + } + LAGraph_Free((void **)&rsm->final_states, msg); + } + + LAGraph_Free((void **)&rsm->start_states, msg); + LAGraph_Free((void **)&rsm, msg); +} + +static void free_graph(void) +{ + if (graph_term == NULL) + return; + for (size_t i = 0; i < n_graph_term; ++i) + GrB_free(&graph_term[i]); + LAGraph_Free((void **)&graph_term, msg); + n_graph_term = 0; + V = 0; +} + +static void free_workspace(void) +{ + GrB_free(&result); + free_rsm(); + free_graph(); } -static void check_reachable(GrB_Vector result, GrB_Index *expected, GrB_Index n_expected, - GrB_Index V, const char *test_name) +static void check_reachable(GrB_Index *expected, GrB_Index n_expected, const char *test_name) { GrB_Index nvals = 0; OK(GrB_Vector_nvals(&nvals, result)); @@ -46,6 +122,8 @@ static void check_reachable(GrB_Vector result, GrB_Index *expected, GrB_Index n_ for (GrB_Index k = 0; k < nvals; ++k) if (val[k]) got[idx[k]] = true; + free(idx); + free(val); } // every expected vertex must be reachable @@ -65,92 +143,210 @@ static void check_reachable(GrB_Vector result, GrB_Index *expected, GrB_Index n_ (unsigned long long)n_expected, (unsigned long long)actual_count); free(got); - free(idx); - free(val); } //============================================================================== // RSM builders //============================================================================== // Grammar 1: S -> a S b | a b -static void rsm_aSb_ab(GrB_Matrix *term, GrB_Matrix *nt, GrB_Vector *start, GrB_Vector *final) +static void init_rsm_aSb_ab(void) { GrB_Index Q = 4; - OK(GrB_Matrix_new(&term[0], GrB_BOOL, Q, Q)); - OK(GrB_Matrix_new(&term[1], GrB_BOOL, Q, Q)); - OK(GrB_Matrix_new(&nt[0], GrB_BOOL, Q, Q)); - OK(GrB_Vector_new(&start[0], GrB_BOOL, Q)); - OK(GrB_Vector_new(&final[0], GrB_BOOL, Q)); - - OK(GrB_Matrix_setElement_BOOL(term[0], true, 0, 1)); // 0 -a-> 1 - OK(GrB_Matrix_setElement_BOOL(term[1], true, 1, 3)); // 1 -b-> 3 - OK(GrB_Matrix_setElement_BOOL(term[1], true, 2, 3)); // 2 -b-> 3 - OK(GrB_Matrix_setElement_BOOL(nt[0], true, 1, 2)); // 1 -S-> 2 - OK(GrB_Vector_setElement_BOOL(start[0], true, 0)); - OK(GrB_Vector_setElement_BOOL(final[0], true, 3)); + rsm = alloc_rsm(Q, /*term_count=*/2, /*nonterm_count=*/1, /*start_nonterm=*/0); + + OK(GrB_Matrix_new(&rsm->terminal_matrices[0], GrB_BOOL, Q, Q)); + OK(GrB_Matrix_new(&rsm->terminal_matrices[1], GrB_BOOL, Q, Q)); + OK(GrB_Matrix_new(&rsm->nonterminal_matrices[0], GrB_BOOL, Q, Q)); + OK(GrB_Vector_new(&rsm->final_states[0], GrB_BOOL, Q)); + + OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[0], true, 0, 1)); // 0 -a-> 1 + OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[1], true, 1, 3)); // 1 -b-> 3 + OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[1], true, 2, 3)); // 2 -b-> 3 + OK(GrB_Matrix_setElement_BOOL(rsm->nonterminal_matrices[0], true, 1, 2)); // 1 -S-> 2 + rsm->start_states[0] = 0; + OK(GrB_Vector_setElement_BOOL(rsm->final_states[0], true, 3)); } // Grammar 2: S -> a S b | c -static void rsm_aSb_c(GrB_Matrix *term, GrB_Matrix *nt, GrB_Vector *start, GrB_Vector *final) +static void init_rsm_aSb_c(void) { GrB_Index Q = 4; - OK(GrB_Matrix_new(&term[0], GrB_BOOL, Q, Q)); - OK(GrB_Matrix_new(&term[1], GrB_BOOL, Q, Q)); - OK(GrB_Matrix_new(&term[2], GrB_BOOL, Q, Q)); - OK(GrB_Matrix_new(&nt[0], GrB_BOOL, Q, Q)); - OK(GrB_Vector_new(&start[0], GrB_BOOL, Q)); - OK(GrB_Vector_new(&final[0], GrB_BOOL, Q)); - - OK(GrB_Matrix_setElement_BOOL(term[0], true, 0, 1)); // 0 -a-> 1 - OK(GrB_Matrix_setElement_BOOL(term[2], true, 0, 3)); // 1 -c-> 3 - OK(GrB_Matrix_setElement_BOOL(term[1], true, 2, 3)); // 2 -b-> 3 - OK(GrB_Matrix_setElement_BOOL(nt[0], true, 1, 2)); // 1 -S-> 2 - OK(GrB_Vector_setElement_BOOL(start[0], true, 0)); - OK(GrB_Vector_setElement_BOOL(final[0], true, 3)); + rsm = alloc_rsm(Q, /*term_count=*/3, /*nonterm_count=*/1, /*start_nonterm=*/0); + + OK(GrB_Matrix_new(&rsm->terminal_matrices[0], GrB_BOOL, Q, Q)); + OK(GrB_Matrix_new(&rsm->terminal_matrices[1], GrB_BOOL, Q, Q)); + OK(GrB_Matrix_new(&rsm->terminal_matrices[2], GrB_BOOL, Q, Q)); + OK(GrB_Matrix_new(&rsm->nonterminal_matrices[0], GrB_BOOL, Q, Q)); + + OK(GrB_Vector_new(&rsm->final_states[0], GrB_BOOL, Q)); + + OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[0], true, 0, 1)); // 0 -a-> 1 + OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[2], true, 0, 3)); // 0 -c-> 3 + OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[1], true, 2, 3)); // 2 -b-> 3 + OK(GrB_Matrix_setElement_BOOL(rsm->nonterminal_matrices[0], true, 1, 2)); // 1 -S-> 2 + rsm->start_states[0] = 0; + OK(GrB_Vector_setElement_BOOL(rsm->final_states[0], true, 3)); } // Grammar 3: S -> X | Y | S X | S Y // X -> a S b | a b // Y -> c S d | c d -static void rsm_complex(GrB_Matrix *term, GrB_Matrix *nt, GrB_Vector *start, GrB_Vector *final) +static void init_rsm_complex(void) { GrB_Index Q = 11; + rsm = alloc_rsm(Q, /*term_count=*/4, /*nonterm_count=*/3, /*start_nonterm=*/0); for (int t = 0; t < 4; ++t) - OK(GrB_Matrix_new(&term[t], GrB_BOOL, Q, Q)); + OK(GrB_Matrix_new(&rsm->terminal_matrices[t], GrB_BOOL, Q, Q)); for (int n = 0; n < 3; ++n) { - OK(GrB_Matrix_new(&nt[n], GrB_BOOL, Q, Q)); - OK(GrB_Vector_new(&start[n], GrB_BOOL, Q)); - OK(GrB_Vector_new(&final[n], GrB_BOOL, Q)); + OK(GrB_Matrix_new(&rsm->nonterminal_matrices[n], GrB_BOOL, Q, Q)); + OK(GrB_Vector_new(&rsm->final_states[n], GrB_BOOL, Q)); } // --- S Component (States 0-2) --- // S -> X | Y | S X | S Y - OK(GrB_Matrix_setElement_BOOL(nt[1], true, 0, 2)); // 0 -X-> 2 - OK(GrB_Matrix_setElement_BOOL(nt[2], true, 0, 2)); // 0 -Y-> 2 - OK(GrB_Matrix_setElement_BOOL(nt[0], true, 0, 1)); // 0 -S-> 1 - OK(GrB_Matrix_setElement_BOOL(nt[1], true, 1, 2)); // 1 -X-> 2 - OK(GrB_Matrix_setElement_BOOL(nt[2], true, 1, 2)); // 1 -Y-> 2 - OK(GrB_Vector_setElement_BOOL(start[0], true, 0)); - OK(GrB_Vector_setElement_BOOL(final[0], true, 2)); + OK(GrB_Matrix_setElement_BOOL(rsm->nonterminal_matrices[1], true, 0, 2)); // 0 -X-> 2 + OK(GrB_Matrix_setElement_BOOL(rsm->nonterminal_matrices[2], true, 0, 2)); // 0 -Y-> 2 + OK(GrB_Matrix_setElement_BOOL(rsm->nonterminal_matrices[0], true, 0, 1)); // 0 -S-> 1 + OK(GrB_Matrix_setElement_BOOL(rsm->nonterminal_matrices[1], true, 1, 2)); // 1 -X-> 2 + OK(GrB_Matrix_setElement_BOOL(rsm->nonterminal_matrices[2], true, 1, 2)); // 1 -Y-> 2 + rsm->start_states[0] = 0; + OK(GrB_Vector_setElement_BOOL(rsm->final_states[0], true, 2)); // --- X Component (States 3-6) --- - OK(GrB_Matrix_setElement_BOOL(term[0], true, 3, 4)); // 3 -a-> 4 - OK(GrB_Matrix_setElement_BOOL(term[1], true, 4, 6)); // 4 -b-> 6 - OK(GrB_Matrix_setElement_BOOL(term[1], true, 5, 6)); // 5 -b-> 6 - OK(GrB_Matrix_setElement_BOOL(nt[0], true, 4, 5)); // 4 -S-> 5 - OK(GrB_Vector_setElement_BOOL(start[1], true, 3)); - OK(GrB_Vector_setElement_BOOL(final[1], true, 6)); + OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[0], true, 3, 4)); // 3 -a-> 4 + OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[1], true, 4, 6)); // 4 -b-> 6 + OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[1], true, 5, 6)); // 5 -b-> 6 + OK(GrB_Matrix_setElement_BOOL(rsm->nonterminal_matrices[0], true, 4, 5)); // 4 -S-> 5 + rsm->start_states[1] = 3; + OK(GrB_Vector_setElement_BOOL(rsm->final_states[1], true, 6)); // --- Y Component (States 7-10) - OK(GrB_Matrix_setElement_BOOL(term[2], true, 7, 8)); // 7 -c-> 8 - OK(GrB_Matrix_setElement_BOOL(term[3], true, 8, 10)); // 8 -d-> 10 - OK(GrB_Matrix_setElement_BOOL(term[3], true, 9, 10)); // 9 -d-> 10 - OK(GrB_Matrix_setElement_BOOL(nt[0], true, 8, 9)); // 8 -S-> 9 - OK(GrB_Vector_setElement_BOOL(start[2], true, 7)); - OK(GrB_Vector_setElement_BOOL(final[2], true, 10)); + OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[2], true, 7, 8)); // 7 -c-> 8 + OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[3], true, 8, 10)); // 8 -d-> 10 + OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[3], true, 9, 10)); // 9 -d-> 10 + OK(GrB_Matrix_setElement_BOOL(rsm->nonterminal_matrices[0], true, 8, 9)); // 8 -S-> 9 + rsm->start_states[2] = 7; + OK(GrB_Vector_setElement_BOOL(rsm->final_states[2], true, 10)); } +//============================================================================== +// Graph builders +//============================================================================== +// 0-a->1; 1-a->2; 2-a->0; 0-b->3; 3-b->0 (V=4, 2 labels: a=0, b=1) +static void init_graph_double_cycle_ab(void) +{ + V = 4; + n_graph_term = 2; + LAGraph_Calloc((void **)&graph_term, n_graph_term, sizeof(GrB_Matrix), msg); + for (size_t i = 0; i < n_graph_term; ++i) + OK(GrB_Matrix_new(&graph_term[i], GrB_BOOL, V, V)); + + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 1, 2)); + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 2, 0)); + OK(GrB_Matrix_setElement_BOOL(graph_term[1], true, 0, 3)); + OK(GrB_Matrix_setElement_BOOL(graph_term[1], true, 3, 0)); +} + +// 0-a->0; 0-a->1; 1-b->1 (V=2, 2 labels: a=0, b=1) +static void init_graph_self_loops(void) +{ + V = 2; + n_graph_term = 2; + LAGraph_Calloc((void **)&graph_term, n_graph_term, sizeof(GrB_Matrix), msg); + for (size_t i = 0; i < n_graph_term; ++i) + OK(GrB_Matrix_new(&graph_term[i], GrB_BOOL, V, V)); + + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 0, 0)); + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(graph_term[1], true, 1, 1)); +} + +// 0-a->0; 0-b->0 (V=2, 2 labels: a=0, b=1) +static void init_graph_single_node(void) +{ + V = 2; + n_graph_term = 2; + LAGraph_Calloc((void **)&graph_term, n_graph_term, sizeof(GrB_Matrix), msg); + for (size_t i = 0; i < n_graph_term; ++i) + OK(GrB_Matrix_new(&graph_term[i], GrB_BOOL, V, V)); + + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 0, 0)); + OK(GrB_Matrix_setElement_BOOL(graph_term[1], true, 0, 0)); +} + +// 0-a->0; 0-c->1; 1-b->1 (V=2, 3 labels: a=0, b=1, c=2) +static void init_graph_loops_abc(void) +{ + V = 2; + n_graph_term = 3; + LAGraph_Calloc((void **)&graph_term, n_graph_term, sizeof(GrB_Matrix), msg); + for (size_t i = 0; i < n_graph_term; ++i) + OK(GrB_Matrix_new(&graph_term[i], GrB_BOOL, V, V)); + + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 0, 0)); + OK(GrB_Matrix_setElement_BOOL(graph_term[2], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(graph_term[1], true, 0, 0)); +} + +// 0-a->1; 1-a->0; 0-c->2; 2-b->3; 3-b->2 (V=4, 3 labels: a=0, b=1, c=2) +static void init_graph_loops_abc_1(void) +{ + V = 4; + n_graph_term = 3; + LAGraph_Calloc((void **)&graph_term, n_graph_term, sizeof(GrB_Matrix), msg); + for (size_t i = 0; i < n_graph_term; ++i) + OK(GrB_Matrix_new(&graph_term[i], GrB_BOOL, V, V)); + + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 1, 0)); + OK(GrB_Matrix_setElement_BOOL(graph_term[2], true, 0, 2)); + OK(GrB_Matrix_setElement_BOOL(graph_term[1], true, 2, 3)); + OK(GrB_Matrix_setElement_BOOL(graph_term[1], true, 3, 2)); +} + +// 0-a->1; 1-a->0; 0-c->2; 2-b->3; 3-b->4; 4-b->2 (V=5, 3 labels: a=0, b=1, c=2) +static void init_graph_loops_abc_2(void) +{ + V = 5; + n_graph_term = 3; + LAGraph_Calloc((void **)&graph_term, n_graph_term, sizeof(GrB_Matrix), msg); + for (size_t i = 0; i < n_graph_term; ++i) + OK(GrB_Matrix_new(&graph_term[i], GrB_BOOL, V, V)); + + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 1, 0)); + OK(GrB_Matrix_setElement_BOOL(graph_term[2], true, 0, 2)); + OK(GrB_Matrix_setElement_BOOL(graph_term[1], true, 2, 3)); + OK(GrB_Matrix_setElement_BOOL(graph_term[1], true, 3, 4)); + OK(GrB_Matrix_setElement_BOOL(graph_term[1], true, 4, 2)); +} + +// a-cycle(0,1,2), b-cycle(2,3), c-cycle(3,4,5), d-edges(2<->5) (V=6, 4 labels) +static void init_graph_multi_cycle(void) +{ + V = 6; + n_graph_term = 4; + LAGraph_Calloc((void **)&graph_term, n_graph_term, sizeof(GrB_Matrix), msg); + for (size_t i = 0; i < n_graph_term; ++i) + OK(GrB_Matrix_new(&graph_term[i], GrB_BOOL, V, V)); + + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 1, 2)); + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 2, 0)); + OK(GrB_Matrix_setElement_BOOL(graph_term[1], true, 2, 3)); + OK(GrB_Matrix_setElement_BOOL(graph_term[1], true, 3, 2)); + OK(GrB_Matrix_setElement_BOOL(graph_term[2], true, 3, 4)); + OK(GrB_Matrix_setElement_BOOL(graph_term[2], true, 4, 5)); + OK(GrB_Matrix_setElement_BOOL(graph_term[2], true, 5, 3)); + OK(GrB_Matrix_setElement_BOOL(graph_term[3], true, 5, 2)); + OK(GrB_Matrix_setElement_BOOL(graph_term[3], true, 2, 5)); +} + +#define run_algorithm(start_vertex) \ + LAGraph_CFPQ_RSM(&result, rsm, graph_term, (start_vertex), V, msg); + //============================================================================== // test cases //============================================================================== @@ -164,48 +360,18 @@ void test_TC1_cyclic_ab(void) { setup(); - GrB_Matrix term[2], nt[1]; - GrB_Vector start[1], final[1]; - - term[0] = term[1] = GrB_NULL; - nt[0] = GrB_NULL; - start[0] = final[0] = GrB_NULL; - - rsm_aSb_ab(term, nt, start, final); - - GrB_Index V = 4; - GrB_Matrix gterm[2]; - gterm[0] = gterm[1] = GrB_NULL; - OK(GrB_Matrix_new(>erm[0], GrB_BOOL, V, V)); - OK(GrB_Matrix_new(>erm[1], GrB_BOOL, V, V)); - - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 1)); - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 1, 2)); - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 2, 0)); - OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 0, 3)); - OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 3, 0)); - - GrB_Vector result = GrB_NULL; + init_rsm_aSb_ab(); + init_graph_double_cycle_ab(); - GrB_Info info = LAGraph_CFPQ_RSM(&result, 2, term, gterm, 1, nt, start, final, 0, 1, 4, V, msg); + GrB_Info info = run_algorithm(/*start_vertex=*/1); if (info == GrB_SUCCESS) { GrB_Index expected[] = {0, 3}; - check_reachable(result, expected, 2, V, "TC1"); + check_reachable(expected, /*n_expected=*/2, "TC1"); } - GrB_free(&result); - for (int i = 0; i < 2; ++i) - { - GrB_free(&term[i]); - GrB_free(>erm[i]); - } - GrB_free(&nt[0]); - GrB_free(&start[0]); - GrB_free(&final[0]); - + free_workspace(); OK(info); - teardown(); } @@ -217,47 +383,18 @@ void test_TC1_cyclic_ab(void) void test_TC2_self_loops(void) { setup(); + init_rsm_aSb_ab(); + init_graph_self_loops(); - GrB_Matrix term[2], nt[1]; - GrB_Vector start[1], final[1]; - - term[0] = term[1] = GrB_NULL; - nt[0] = GrB_NULL; - start[0] = final[0] = GrB_NULL; - - rsm_aSb_ab(term, nt, start, final); - - GrB_Index V = 2; - GrB_Matrix gterm[2]; - gterm[0] = gterm[1] = GrB_NULL; - OK(GrB_Matrix_new(>erm[0], GrB_BOOL, V, V)); - OK(GrB_Matrix_new(>erm[1], GrB_BOOL, V, V)); - - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 0)); - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 1)); - OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 1, 1)); - - GrB_Vector result = GrB_NULL; - - GrB_Info info = LAGraph_CFPQ_RSM(&result, 2, term, gterm, 1, nt, start, final, 0, 0, 4, V, msg); + GrB_Info info = run_algorithm(/*start_vertex=*/0); if (info == GrB_SUCCESS) { GrB_Index expected[] = {1}; - check_reachable(result, expected, 1, V, "TC2"); + check_reachable(expected, /*n_expected=*/1, "TC2"); } - GrB_free(&result); - for (int i = 0; i < 2; ++i) - { - GrB_free(&term[i]); - GrB_free(>erm[i]); - } - GrB_free(&nt[0]); - GrB_free(&start[0]); - GrB_free(&final[0]); - + free_workspace(); OK(info); - teardown(); } @@ -269,46 +406,18 @@ void test_TC2_self_loops(void) void test_TC3_single_node(void) { setup(); + init_rsm_aSb_ab(); + init_graph_single_node(); - GrB_Matrix term[2], nt[1]; - GrB_Vector start[1], final[1]; - - term[0] = term[1] = GrB_NULL; - nt[0] = GrB_NULL; - start[0] = final[0] = GrB_NULL; - - rsm_aSb_ab(term, nt, start, final); - - GrB_Index V = 1; - GrB_Matrix gterm[2]; - gterm[0] = gterm[1] = GrB_NULL; - OK(GrB_Matrix_new(>erm[0], GrB_BOOL, V, V)); - OK(GrB_Matrix_new(>erm[1], GrB_BOOL, V, V)); - - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 0)); - OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 0, 0)); - - GrB_Vector result = GrB_NULL; - - GrB_Info info = LAGraph_CFPQ_RSM(&result, 2, term, gterm, 1, nt, start, final, 0, 0, 4, V, msg); + GrB_Info info = run_algorithm(/*start_vertex=*/0); if (info == GrB_SUCCESS) { GrB_Index expected[] = {0}; - check_reachable(result, expected, 1, V, "TC3"); - } - - GrB_free(&result); - for (int i = 0; i < 2; ++i) - { - GrB_free(&term[i]); - GrB_free(>erm[i]); + check_reachable(expected, /*n_expected=*/1, "TC3"); } - GrB_free(&nt[0]); - GrB_free(&start[0]); - GrB_free(&final[0]); + free_workspace(); OK(info); - teardown(); } @@ -320,47 +429,18 @@ void test_TC3_single_node(void) void test_TC4_aSb_c_loops(void) { setup(); + init_rsm_aSb_c(); + init_graph_loops_abc(); - GrB_Matrix term[3], nt[1]; - GrB_Vector start[1], final[1]; - - term[0] = term[1] = term[2] = GrB_NULL; - nt[0] = GrB_NULL; - start[0] = final[0] = GrB_NULL; - - rsm_aSb_c(term, nt, start, final); - - GrB_Index V = 2; - GrB_Matrix gterm[3]; - gterm[0] = gterm[1] = gterm[2] = GrB_NULL; - for (int i = 0; i < 3; ++i) - OK(GrB_Matrix_new(>erm[i], GrB_BOOL, V, V)); - - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 0)); - OK(GrB_Matrix_setElement_BOOL(gterm[2], true, 0, 1)); - OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 1, 1)); - - GrB_Vector result = GrB_NULL; - - GrB_Info info = LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, start, final, 0, 0, 4, V, msg); + GrB_Info info = run_algorithm(/*start_vertex=*/0); if (info == GrB_SUCCESS) { GrB_Index expected[] = {1}; - check_reachable(result, expected, 1, V, "TC4"); - } - - GrB_free(&result); - for (int i = 0; i < 3; ++i) - { - GrB_free(&term[i]); - GrB_free(>erm[i]); + check_reachable(expected, 1, "TC4"); } - GrB_free(&nt[0]); - GrB_free(&start[0]); - GrB_free(&final[0]); + free_workspace(); OK(info); - teardown(); } @@ -372,49 +452,18 @@ void test_TC4_aSb_c_loops(void) void test_TC5_aSb_c_cycle1(void) { setup(); + init_rsm_aSb_c(); + init_graph_loops_abc_1(); - GrB_Matrix term[3], nt[1]; - GrB_Vector start[1], final[1]; - - term[0] = term[1] = term[2] = GrB_NULL; - nt[0] = GrB_NULL; - start[0] = final[0] = GrB_NULL; - - rsm_aSb_c(term, nt, start, final); - - GrB_Index V = 4; - GrB_Matrix gterm[3]; - gterm[0] = gterm[1] = gterm[2] = GrB_NULL; - for (int i = 0; i < 3; ++i) - OK(GrB_Matrix_new(>erm[i], GrB_BOOL, V, V)); - - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 1)); - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 1, 0)); - OK(GrB_Matrix_setElement_BOOL(gterm[2], true, 0, 2)); - OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 2, 3)); - OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 3, 2)); - - GrB_Vector result = GrB_NULL; - - GrB_Info info = LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, start, final, 0, 1, 4, V, msg); + GrB_Info info = run_algorithm(/*start_vertex=*/1); if (info == GrB_SUCCESS) { GrB_Index expected[] = {3}; - check_reachable(result, expected, 1, V, "TC5"); - } - - GrB_free(&result); - for (int i = 0; i < 3; ++i) - { - GrB_free(&term[i]); - GrB_free(>erm[i]); + check_reachable(expected, 1, "TC5"); } - GrB_free(&nt[0]); - GrB_free(&start[0]); - GrB_free(&final[0]); + free_workspace(); OK(info); - teardown(); } //------------------------------------------------------------------------------ @@ -425,110 +474,51 @@ void test_TC5_aSb_c_cycle1(void) void test_TC6_aSb_c_cycle2(void) { setup(); + init_rsm_aSb_c(); - GrB_Matrix term[3], nt[1]; - GrB_Vector start[1], final[1]; - - term[0] = term[1] = term[2] = GrB_NULL; - nt[0] = GrB_NULL; - start[0] = final[0] = GrB_NULL; - - rsm_aSb_c(term, nt, start, final); + V = 3; + n_graph_term = 3; + LAGraph_Calloc((void **)&graph_term, n_graph_term, sizeof(GrB_Matrix), msg); + for (size_t i = 0; i < n_graph_term; ++i) + OK(GrB_Matrix_new(&graph_term[i], GrB_BOOL, V, V)); - GrB_Index V = 3; - GrB_Matrix gterm[3]; - gterm[0] = gterm[1] = gterm[2] = GrB_NULL; - for (int i = 0; i < 3; ++i) - OK(GrB_Matrix_new(>erm[i], GrB_BOOL, V, V)); + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(graph_term[0], true, 1, 0)); + OK(GrB_Matrix_setElement_BOOL(graph_term[2], true, 0, 2)); + OK(GrB_Matrix_setElement_BOOL(graph_term[1], true, 2, 2)); - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 1)); - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 1, 0)); - OK(GrB_Matrix_setElement_BOOL(gterm[2], true, 0, 2)); - OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 2, 2)); - - GrB_Vector result = GrB_NULL; - - GrB_Info info = LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, start, final, 0, 1, 4, V, msg); + GrB_Info info = run_algorithm(/*start_vertex=*/1); if (info == GrB_SUCCESS) { GrB_Index expected[] = {2}; - check_reachable(result, expected, 1, V, "TC6"); - } - - GrB_free(&result); - for (int i = 0; i < 3; ++i) - { - GrB_free(&term[i]); - GrB_free(>erm[i]); + check_reachable(expected, 1, "TC6"); } - GrB_free(&nt[0]); - GrB_free(&start[0]); - GrB_free(&final[0]); + free_workspace(); OK(info); - teardown(); } -// Graph: 0-a->1; 1-a->0; 0-c->2; 2-b->3; 3-b->4; 4-b->2 -static void graph_tc7_8(GrB_Matrix gterm[3]) -{ - GrB_Index V = 5; - for (int i = 0; i < 3; ++i) - OK(GrB_Matrix_new(>erm[i], GrB_BOOL, V, V)); - - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 1)); - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 1, 0)); - OK(GrB_Matrix_setElement_BOOL(gterm[2], true, 0, 2)); - OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 2, 3)); - OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 3, 4)); - OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 4, 2)); -} - //------------------------------------------------------------------------------ // TC7: S -> a S b | c -// same graph as TC7/8 above +// Graph: 0-a->1; 1-a->0; 0-c->2; 2-b->3; 3-b->4; 4-b->2 // Source: 0 Expected: {2, 3, 4} //------------------------------------------------------------------------------ void test_TC7_aSb_c_cycle3(void) { setup(); + init_rsm_aSb_c(); + init_graph_loops_abc_2(); - GrB_Matrix term[3], nt[1]; - GrB_Vector start[1], final[1]; - - term[0] = term[1] = term[2] = GrB_NULL; - nt[0] = GrB_NULL; - start[0] = final[0] = GrB_NULL; - - rsm_aSb_c(term, nt, start, final); - - GrB_Index V = 5; - GrB_Matrix gterm[3]; - gterm[0] = gterm[1] = gterm[2] = GrB_NULL; - graph_tc7_8(gterm); - - GrB_Vector result = GrB_NULL; - - GrB_Info info = LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, start, final, 0, 0, 4, V, msg); + GrB_Info info = run_algorithm(/*start_vertex=*/0); if (info == GrB_SUCCESS) { GrB_Index expected[] = {2, 3, 4}; - check_reachable(result, expected, 3, V, "TC7"); + check_reachable(expected, 3, "TC7"); } - GrB_free(&result); - for (int i = 0; i < 3; ++i) - { - GrB_free(&term[i]); - GrB_free(>erm[i]); - } - GrB_free(&nt[0]); - GrB_free(&start[0]); - GrB_free(&final[0]); - + free_workspace(); OK(info); - teardown(); } @@ -540,44 +530,18 @@ void test_TC7_aSb_c_cycle3(void) void test_TC8_aSb_c_cycle4(void) { setup(); + init_rsm_aSb_c(); + init_graph_loops_abc_2(); - GrB_Matrix term[3], nt[1], call[1]; - GrB_Vector start[1], final[1]; - - term[0] = term[1] = term[2] = GrB_NULL; - nt[0] = GrB_NULL; - call[0] = GrB_NULL; - start[0] = final[0] = GrB_NULL; - - rsm_aSb_c(term, nt, start, final); - - GrB_Index V = 5; - GrB_Matrix gterm[3]; - gterm[0] = gterm[1] = gterm[2] = GrB_NULL; - graph_tc7_8(gterm); - - GrB_Vector result = GrB_NULL; - - GrB_Info info = LAGraph_CFPQ_RSM(&result, 3, term, gterm, 1, nt, start, final, 0, 1, 4, V, msg); + GrB_Info info = run_algorithm(/*start_vertex=*/1); if (info == GrB_SUCCESS) { GrB_Index expected[] = {2, 3, 4}; - check_reachable(result, expected, 3, V, "TC8"); + check_reachable(expected, 3, "TC8"); } - GrB_free(&result); - for (int i = 0; i < 3; ++i) - { - GrB_free(&term[i]); - GrB_free(>erm[i]); - } - GrB_free(&nt[0]); - GrB_free(&call[0]); - GrB_free(&start[0]); - GrB_free(&final[0]); - + free_workspace(); OK(info); - teardown(); } @@ -585,62 +549,18 @@ void test_TC9_mult_recursive_cycles(void) { setup(); - GrB_Matrix term[4], nt[3]; - GrB_Vector start[3], final[3]; - - term[0] = term[1] = term[2] = term[3] = GrB_NULL; - nt[0] = nt[1] = nt[2] = GrB_NULL; - start[0] = start[1] = start[2] = GrB_NULL; - final[0] = final[1] = final[2] = GrB_NULL; - - rsm_complex(term, nt, start, final); - - GrB_Index V = 6; - GrB_Matrix gterm[4]; - gterm[0] = gterm[1] = gterm[2] = gterm[3] = GrB_NULL; - for (int i = 0; i < 4; ++i) - OK(GrB_Matrix_new(>erm[i], GrB_BOOL, V, V)); - - // a-cycle (0,1,2) - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 0, 1)); - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 1, 2)); - OK(GrB_Matrix_setElement_BOOL(gterm[0], true, 2, 0)); - // b-cycle (2,3) - OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 2, 3)); - OK(GrB_Matrix_setElement_BOOL(gterm[1], true, 3, 2)); - // c-cycle (3,4,5) - OK(GrB_Matrix_setElement_BOOL(gterm[2], true, 3, 4)); - OK(GrB_Matrix_setElement_BOOL(gterm[2], true, 4, 5)); - OK(GrB_Matrix_setElement_BOOL(gterm[2], true, 5, 3)); - // d-cycle - OK(GrB_Matrix_setElement_BOOL(gterm[3], true, 5, 2)); - OK(GrB_Matrix_setElement_BOOL(gterm[3], true, 2, 5)); - - GrB_Vector result = GrB_NULL; - - GrB_Info info = - LAGraph_CFPQ_RSM(&result, 4, term, gterm, 3, nt, start, final, 0, 0, 11, V, msg); + init_rsm_complex(); + init_graph_multi_cycle(); + + GrB_Info info = run_algorithm(/*start_vertex=*/0); if (info == GrB_SUCCESS) { GrB_Index expected[] = {2, 3, 5}; - check_reachable(result, expected, 3, V, "TC9"); - } - - GrB_free(&result); - for (int i = 0; i < 4; ++i) - { - GrB_free(&term[i]); - GrB_free(>erm[i]); - } - for (int i = 0; i < 3; ++i) - { - GrB_free(&nt[i]); - GrB_free(&start[i]); - GrB_free(&final[i]); + check_reachable(expected, 3, "TC9"); } + free_workspace(); OK(info); - teardown(); } @@ -649,7 +569,7 @@ TEST_LIST = {{"TC1_cyclic_ab", test_TC1_cyclic_ab}, {"TC3_single_node", test_TC3_single_node}, {"TC4_aSb_c_loops", test_TC4_aSb_c_loops}, {"TC5_aSb_c_cycle1", test_TC5_aSb_c_cycle1}, - {"TC4_aSb_c_loops", test_TC4_aSb_c_loops}, + {"TC6_aSb_c_cycle2", test_TC6_aSb_c_cycle2}, {"TC7_aSb_c_cycle3", test_TC7_aSb_c_cycle3}, {"TC8_aSb_c_cycle4", test_TC8_aSb_c_cycle4}, {"TC9_mult_recursive_cycles", test_TC9_mult_recursive_cycles}, diff --git a/include/LAGraphX.h b/include/LAGraphX.h index 2a5a7249a0..f0278fb704 100644 --- a/include/LAGraphX.h +++ b/include/LAGraphX.h @@ -1097,22 +1097,20 @@ GrB_Info LAGraph_CFL_reachability char *msg // Message string for error reporting. ) ; -int LAGraph_CFPQ_RSM( - // output: - GrB_Vector *reachable, - - // RSM terminal input: - size_t num_terminals, GrB_Matrix *rsm_term, GrB_Matrix *graph_term, - - // RSM non-terminal input: - size_t num_nonterminals, GrB_Matrix *rsm_nonterm, GrB_Vector *rsm_start, - GrB_Vector *rsm_final, - - size_t start_nonterm, GrB_Index start_vertex, - - GrB_Index Q, GrB_Index V, - - char *msg); +typedef struct +{ + GrB_Index state_count; + GrB_Index terminal_count; + GrB_Index nonterminal_count; + GrB_Index start_nonterminal; + GrB_Matrix *terminal_matrices; + GrB_Matrix *nonterminal_matrices; + GrB_Index *start_states; + GrB_Vector *final_states; +} RSM; + +int LAGraph_CFPQ_RSM(GrB_Vector *reachable, RSM *rsm, GrB_Matrix *graph_term, + GrB_Index start_vertex, GrB_Index V, char *msg); //------------------------------------------------------------------------------ // a simple example of an algorithm From c9efb510a07e599a4a1514de6056540760ff3c4e Mon Sep 17 00:00:00 2001 From: Sterkachov Date: Sun, 19 Apr 2026 22:08:48 +0300 Subject: [PATCH 09/10] feat: add multiple source support to CFPQ algorithm --- experimental/algorithm/LAGraph_CFPQ_RSM.c | 57 +++++++++++++++-------- experimental/test/test_CFPQ_RSM.c | 33 ++++++++----- include/LAGraphX.h | 4 +- 3 files changed, 61 insertions(+), 33 deletions(-) diff --git a/experimental/algorithm/LAGraph_CFPQ_RSM.c b/experimental/algorithm/LAGraph_CFPQ_RSM.c index 35ac23be58..5eb6317eba 100644 --- a/experimental/algorithm/LAGraph_CFPQ_RSM.c +++ b/experimental/algorithm/LAGraph_CFPQ_RSM.c @@ -23,6 +23,7 @@ GrB_free(&diag_entry); \ GrB_free(&outer_result); \ GrB_free(&entry_vec); \ + GrB_free(&source_mask); \ if (Stack != NULL) \ { \ for (size_t _i = 0; _i < num_nonterm; ++_i) \ @@ -124,8 +125,8 @@ typedef struct GrB_Vector *final_states; // rsm_final } RSM; -int LAGraph_CFPQ_RSM(GrB_Vector *reachable, RSM *rsm, GrB_Matrix *graph_term, - GrB_Index start_vertex, GrB_Index V, char *msg) +int LAGraph_CFPQ_RSM(GrB_Vector *reachable, const RSM *rsm, const GrB_Matrix *graph_term, + const GrB_Index *sources, size_t num_sources, GrB_Index V, char *msg) { LG_CLEAR_MSG; @@ -156,6 +157,7 @@ int LAGraph_CFPQ_RSM(GrB_Vector *reachable, RSM *rsm, GrB_Matrix *graph_term, GrB_Matrix diag_entry = GrB_NULL; // add diagonal matrix for entry_vec GrB_Matrix outer_result = GrB_NULL; // additional matrix GrB_Vector entry_vec = GrB_NULL; + GrB_Vector source_mask = GrB_NULL; GrB_Matrix *Stack = NULL; GrB_Matrix *graph_nt = NULL; @@ -168,7 +170,7 @@ int LAGraph_CFPQ_RSM(GrB_Vector *reachable, RSM *rsm, GrB_Matrix *graph_term, LG_ASSERT(start_states != NULL, GrB_NULL_POINTER); LG_ASSERT(final_states != NULL, GrB_NULL_POINTER); LG_ASSERT(start_nonterm < num_nonterm, GrB_INVALID_VALUE); - LG_ASSERT(start_vertex < V, GrB_INVALID_VALUE); + LG_ASSERT(num_sources <= V, GrB_INVALID_VALUE); *reachable = GrB_NULL; @@ -213,39 +215,43 @@ int LAGraph_CFPQ_RSM(GrB_Vector *reachable, RSM *rsm, GrB_Matrix *graph_term, // GRB_TRY(GrB_Matrix_new(&diag_entry, GrB_BOOL, V, V)); GRB_TRY(GrB_Matrix_new(&outer_result, GrB_BOOL, Q, VV)); GRB_TRY(GrB_Vector_new(&entry_vec, GrB_BOOL, V)); + GRB_TRY(GrB_Vector_new(&source_mask, GrB_BOOL, V)); GRB_TRY(GrB_Vector_new(reachable, GrB_BOOL, V)); // seed initial frontier M and Stack { - GrB_Index seed_col = start_vertex * V + start_vertex; GrB_Index nstart = 0; GRB_TRY(GrB_Vector_nvals(&nstart, rsm_start[start_nonterm])); - GrB_Index *start_states = NULL; - bool *start_vals = NULL; - LG_TRY(LAGraph_Malloc((void **)&start_states, nstart, sizeof(GrB_Index), msg)); - LG_TRY(LAGraph_Malloc((void **)&start_vals, nstart, sizeof(bool), msg)); + GrB_Index *seed_states = NULL; + bool *seed_vals = NULL; + LG_TRY(LAGraph_Malloc((void **)&seed_states, nstart, sizeof(GrB_Index), msg)); + LG_TRY(LAGraph_Malloc((void **)&seed_vals, nstart, sizeof(bool), msg)); - GrB_Info info = GrB_Vector_extractTuples_BOOL(start_states, start_vals, &nstart, + GrB_Info info = GrB_Vector_extractTuples_BOOL(seed_states, seed_vals, &nstart, rsm_start[start_nonterm]); - LAGraph_Free((void **)&start_vals, GrB_NULL); + LAGraph_Free((void **)&seed_vals, GrB_NULL); if (info != GrB_SUCCESS) { - LAGraph_Free((void **)&start_states, GrB_NULL); + LAGraph_Free((void **)&seed_states, GrB_NULL); LG_FREE_ALL; return info; } - for (GrB_Index k = 0; k < nstart; ++k) + for (size_t i = 0; i < num_sources; ++i) { - GrB_Index q = start_states[k]; - GRB_TRY(GrB_Matrix_setElement_BOOL(M, true, q, seed_col)); - GRB_TRY(GrB_Matrix_setElement_BOOL(Stack[start_nonterm], true, q, seed_col)); + GrB_Index seed_col = sources[i] * V + sources[i]; + for (GrB_Index k = 0; k < nstart; ++k) + { + GrB_Index q = seed_states[k]; + GRB_TRY(GrB_Matrix_setElement_BOOL(M, true, q, seed_col)); + GRB_TRY(GrB_Matrix_setElement_BOOL(Stack[start_nonterm], true, q, seed_col)); + } } - LAGraph_Free((void **)&start_states, GrB_NULL); + LAGraph_Free((void **)&seed_states, GrB_NULL); } // main fixed-point loop @@ -388,9 +394,22 @@ int LAGraph_CFPQ_RSM(GrB_Vector *reachable, RSM *rsm, GrB_Matrix *graph_term, GRB_TRY(GrB_Matrix_nvals(&m_nvals, M)); } - // GrB_DESC_T0 transposes graph_nt so Col_extract gives row start_vertex - GRB_TRY(GrB_Col_extract(*reachable, GrB_NULL, GrB_NULL, graph_nt[start_nonterm], GrB_ALL, V, - start_vertex, GrB_DESC_T0)); + // Extract all reachable vertices + + { + bool *vals = NULL; + LG_TRY(LAGraph_Malloc((void **)&vals, num_sources, sizeof(bool), msg)); + for (size_t i = 0; i < num_sources; ++i) + vals[i] = true; + + GrB_Info _bi = + GrB_Vector_build_BOOL(source_mask, sources, vals, num_sources, GrB_SECOND_BOOL); + LAGraph_Free((void **)&vals, GrB_NULL); + GRB_TRY(_bi); + } + + GRB_TRY(GrB_mxv(*reachable, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, + graph_nt[start_nonterm], source_mask, GrB_DESC_T0)); LG_FREE_WORK; return GrB_SUCCESS; diff --git a/experimental/test/test_CFPQ_RSM.c b/experimental/test/test_CFPQ_RSM.c index 344c87ccc5..2ac746ead4 100644 --- a/experimental/test/test_CFPQ_RSM.c +++ b/experimental/test/test_CFPQ_RSM.c @@ -65,7 +65,7 @@ static void free_rsm(void) LAGraph_Free((void **)&rsm->nonterminal_matrices, msg); } - if (rsm->nonterminal_matrices != NULL) + if (rsm->final_states != NULL) { for (size_t i = 0; i < rsm->nonterminal_count; ++i) { @@ -344,8 +344,8 @@ static void init_graph_multi_cycle(void) OK(GrB_Matrix_setElement_BOOL(graph_term[3], true, 2, 5)); } -#define run_algorithm(start_vertex) \ - LAGraph_CFPQ_RSM(&result, rsm, graph_term, (start_vertex), V, msg); +#define run_algorithm(sources, num_sources) \ + LAGraph_CFPQ_RSM(&result, rsm, graph_term, (sources), (num_sources), V, msg); //============================================================================== // test cases @@ -363,7 +363,8 @@ void test_TC1_cyclic_ab(void) init_rsm_aSb_ab(); init_graph_double_cycle_ab(); - GrB_Info info = run_algorithm(/*start_vertex=*/1); + GrB_Index sources = {1}; + GrB_Info info = run_algorithm(&sources, 1); if (info == GrB_SUCCESS) { GrB_Index expected[] = {0, 3}; @@ -386,7 +387,8 @@ void test_TC2_self_loops(void) init_rsm_aSb_ab(); init_graph_self_loops(); - GrB_Info info = run_algorithm(/*start_vertex=*/0); + GrB_Index sources = {0}; + GrB_Info info = run_algorithm(&sources, 1); if (info == GrB_SUCCESS) { GrB_Index expected[] = {1}; @@ -409,7 +411,8 @@ void test_TC3_single_node(void) init_rsm_aSb_ab(); init_graph_single_node(); - GrB_Info info = run_algorithm(/*start_vertex=*/0); + GrB_Index sources = {0}; + GrB_Info info = run_algorithm(&sources, 1); if (info == GrB_SUCCESS) { GrB_Index expected[] = {0}; @@ -432,7 +435,8 @@ void test_TC4_aSb_c_loops(void) init_rsm_aSb_c(); init_graph_loops_abc(); - GrB_Info info = run_algorithm(/*start_vertex=*/0); + GrB_Index sources = {0}; + GrB_Info info = run_algorithm(&sources, 1); if (info == GrB_SUCCESS) { GrB_Index expected[] = {1}; @@ -455,7 +459,8 @@ void test_TC5_aSb_c_cycle1(void) init_rsm_aSb_c(); init_graph_loops_abc_1(); - GrB_Info info = run_algorithm(/*start_vertex=*/1); + GrB_Index sources = {1}; + GrB_Info info = run_algorithm(&sources, 1); if (info == GrB_SUCCESS) { GrB_Index expected[] = {3}; @@ -487,7 +492,8 @@ void test_TC6_aSb_c_cycle2(void) OK(GrB_Matrix_setElement_BOOL(graph_term[2], true, 0, 2)); OK(GrB_Matrix_setElement_BOOL(graph_term[1], true, 2, 2)); - GrB_Info info = run_algorithm(/*start_vertex=*/1); + GrB_Index sources = {1}; + GrB_Info info = run_algorithm(&sources, 1); if (info == GrB_SUCCESS) { GrB_Index expected[] = {2}; @@ -510,7 +516,8 @@ void test_TC7_aSb_c_cycle3(void) init_rsm_aSb_c(); init_graph_loops_abc_2(); - GrB_Info info = run_algorithm(/*start_vertex=*/0); + GrB_Index sources = {0}; + GrB_Info info = run_algorithm(&sources, 1); if (info == GrB_SUCCESS) { GrB_Index expected[] = {2, 3, 4}; @@ -533,7 +540,8 @@ void test_TC8_aSb_c_cycle4(void) init_rsm_aSb_c(); init_graph_loops_abc_2(); - GrB_Info info = run_algorithm(/*start_vertex=*/1); + GrB_Index sources = {1}; + GrB_Info info = run_algorithm(&sources, 1); if (info == GrB_SUCCESS) { GrB_Index expected[] = {2, 3, 4}; @@ -552,7 +560,8 @@ void test_TC9_mult_recursive_cycles(void) init_rsm_complex(); init_graph_multi_cycle(); - GrB_Info info = run_algorithm(/*start_vertex=*/0); + GrB_Index sources = {1}; + GrB_Info info = run_algorithm(&sources, 1); if (info == GrB_SUCCESS) { GrB_Index expected[] = {2, 3, 5}; diff --git a/include/LAGraphX.h b/include/LAGraphX.h index f0278fb704..18171ed555 100644 --- a/include/LAGraphX.h +++ b/include/LAGraphX.h @@ -1109,8 +1109,8 @@ typedef struct GrB_Vector *final_states; } RSM; -int LAGraph_CFPQ_RSM(GrB_Vector *reachable, RSM *rsm, GrB_Matrix *graph_term, - GrB_Index start_vertex, GrB_Index V, char *msg); +int LAGraph_CFPQ_RSM(GrB_Vector *reachable, const RSM *rsm, const GrB_Matrix *graph_term, + const GrB_Index *sources, size_t num_sources, GrB_Index V, char *msg); //------------------------------------------------------------------------------ // a simple example of an algorithm From 35bee56067b09a3e0338723058246683afdefa74 Mon Sep 17 00:00:00 2001 From: Sterkachov Date: Thu, 23 Apr 2026 23:45:48 +0300 Subject: [PATCH 10/10] fix: Stack structure is redundant, use P instead --- experimental/algorithm/LAGraph_CFPQ_RSM.c | 75 +++++++++++------------ experimental/test/test_CFPQ_RSM.c | 28 +++------ 2 files changed, 44 insertions(+), 59 deletions(-) diff --git a/experimental/algorithm/LAGraph_CFPQ_RSM.c b/experimental/algorithm/LAGraph_CFPQ_RSM.c index 5eb6317eba..2f47f58caa 100644 --- a/experimental/algorithm/LAGraph_CFPQ_RSM.c +++ b/experimental/algorithm/LAGraph_CFPQ_RSM.c @@ -3,6 +3,7 @@ #include #include +#include #include #include @@ -247,6 +248,7 @@ int LAGraph_CFPQ_RSM(GrB_Vector *reachable, const RSM *rsm, const GrB_Matrix *gr { GrB_Index q = seed_states[k]; GRB_TRY(GrB_Matrix_setElement_BOOL(M, true, q, seed_col)); + GRB_TRY(GrB_Matrix_setElement_BOOL(P, true, q, seed_col)); GRB_TRY(GrB_Matrix_setElement_BOOL(Stack[start_nonterm], true, q, seed_col)); } } @@ -275,6 +277,12 @@ int LAGraph_CFPQ_RSM(GrB_Vector *reachable, const RSM *rsm, const GrB_Matrix *gr { if (rsm_nonterm[i] == GrB_NULL) continue; + + GrB_Index gnt_nvals = 0; + GRB_TRY(GrB_Matrix_nvals(&gnt_nvals, graph_nt[i])); + if (gnt_nvals == 0) + continue; + LG_TRY(s_mxm_chain(M_nonterm, rsm_nonterm[i], M, graph_nt[i], temp1, temp2, Q, V)); } @@ -287,7 +295,7 @@ int LAGraph_CFPQ_RSM(GrB_Vector *reachable, const RSM *rsm, const GrB_Matrix *gr // mask_call = rsm_call[i]^T * M -> |Q|*|V*V| GRB_TRY(GrB_Matrix_clear(mask_call)); - GRB_TRY(GrB_mxm(mask_call, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, rsm_call[i], + GRB_TRY(GrB_mxm(mask_call, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, rsm_call[i], M, GrB_DESC_T0)); GrB_Index mc_nvals = 0; @@ -295,27 +303,26 @@ int LAGraph_CFPQ_RSM(GrB_Vector *reachable, const RSM *rsm, const GrB_Matrix *gr if (mc_nvals == 0) continue; - // frame = rsm_nonterm[i]^T * M, accumulated into Stack[i] - GRB_TRY(GrB_Matrix_clear(frame)); - GRB_TRY(GrB_mxm(frame, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, rsm_nonterm[i], M, - GrB_DESC_T0)); - GRB_TRY(GrB_eWiseAdd(Stack[i], GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, Stack[i], - frame, GrB_NULL)); - // entry_vec = column-OR of mask_call reshaped to |Q*V|*|V| GRB_TRY(GrB_Matrix_clear(temp1)); - GRB_TRY(GrB_eWiseAdd(temp1, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, temp1, + GRB_TRY(GrB_eWiseAdd(temp1, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, temp1, mask_call, GrB_NULL)); GRB_TRY(GxB_Matrix_reshape(temp1, false, Q * V, V, GrB_NULL)); GRB_TRY(GrB_Vector_clear(entry_vec)); // column OR = row OR of transpose GRB_TRY( - GrB_reduce(entry_vec, GrB_NULL, GrB_LOR, GrB_LOR_MONOID_BOOL, temp1, GrB_DESC_T0)); + GrB_reduce(entry_vec, GrB_NULL, GrB_NULL, GrB_LOR_MONOID_BOOL, temp1, GrB_DESC_T0)); GRB_TRY(GxB_Matrix_reshape(temp1, false, Q, VV, GrB_NULL)); + GrB_Index ev_nvals = 0; + GRB_TRY(GrB_Vector_nvals(&ev_nvals, entry_vec)); + if (ev_nvals == 0) + continue; + // Build diag(entry_vec) reshaped to 1*|V*V| + GRB_TRY(GrB_Matrix_free(&diag_entry)); diag_entry = GrB_NULL; GRB_TRY(GrB_Matrix_diag(&diag_entry, entry_vec, 0)); // |V|*|V| GRB_TRY(GxB_Matrix_reshape(diag_entry, false, 1, VV, GrB_NULL)); @@ -326,17 +333,15 @@ int LAGraph_CFPQ_RSM(GrB_Vector *reachable, const RSM *rsm, const GrB_Matrix *gr (GrB_Matrix)rsm_start[i], diag_entry, GrB_NULL)); // M_call |= outer_result - GRB_TRY(GrB_eWiseAdd(M_call, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, M_call, - outer_result, GrB_NULL)); - - GRB_TRY(GrB_Matrix_free(&diag_entry)); + GRB_TRY( + GrB_eWiseAdd(M_call, GrB_NULL, GrB_LOR, GrB_LOR, M_call, outer_result, GrB_NULL)); } // PHASE 4: RETURN TRANSITIONS GRB_TRY(GrB_Matrix_clear(M_return)); for (size_t i = 0; i < num_nonterm; ++i) { - if (final_states[i] == GrB_NULL) + if (final_states[i] == GrB_NULL || rsm_nonterm[i] == GrB_NULL) continue; // 1*|Q| * |Q|*|V*V| -> 1*|V*V| @@ -346,25 +351,19 @@ int LAGraph_CFPQ_RSM(GrB_Vector *reachable, const RSM *rsm, const GrB_Matrix *gr GRB_TRY(GxB_Matrix_reshape(new_edges, false, V, V, GrB_NULL)); - GRB_TRY(GrB_eWiseAdd(graph_nt[i], GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, - graph_nt[i], new_edges, GrB_NULL)); - - GRB_TRY(GrB_Matrix_clear(temp1)); - GRB_TRY(GrB_eWiseAdd(temp1, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, temp1, - Stack[i], GrB_NULL)); - GRB_TRY(GxB_Matrix_reshape(temp1, false, Q * V, V, GrB_NULL)); + GrB_Index ne_vals = 0; + GRB_TRY(GrB_Matrix_nvals(&ne_vals, new_edges)); + if (ne_vals == 0) + { + GRB_TRY(GxB_Matrix_reshape(new_edges, false, 1, VV, GrB_NULL)); + continue; + } - GRB_TRY(GrB_Matrix_clear(temp2)); - GRB_TRY(GxB_Matrix_reshape(temp2, false, Q * V, V, GrB_NULL)); - // |Q*V|*|V| * |V|*|V| -> |Q*V|*|V| - GRB_TRY(GrB_mxm(temp2, GrB_NULL, GrB_NULL, GrB_LOR_LAND_SEMIRING_BOOL, temp1, new_edges, - GrB_NULL)); + GRB_TRY(GrB_eWiseAdd(graph_nt[i], GrB_NULL, GrB_LOR, GrB_LOR, graph_nt[i], new_edges, + GrB_NULL)); - GRB_TRY(GxB_Matrix_reshape(temp1, false, Q, VV, GrB_NULL)); - GRB_TRY(GxB_Matrix_reshape(temp2, false, Q, VV, GrB_NULL)); + LG_TRY(s_mxm_chain(M_return, rsm_nonterm[i], P, new_edges, temp1, temp2, Q, V)); - GRB_TRY(GrB_eWiseAdd(M_return, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, M_return, - temp2, GrB_NULL)); GRB_TRY(GxB_Matrix_reshape(new_edges, false, 1, VV, GrB_NULL)); } @@ -372,20 +371,16 @@ int LAGraph_CFPQ_RSM(GrB_Vector *reachable, const RSM *rsm, const GrB_Matrix *gr // GrB_assign has big overhead // M_new = (M_term | M_nonterm | M_call | M_return) & ~P GRB_TRY(GrB_Matrix_clear(M_new)); - GRB_TRY(GrB_eWiseAdd(M_new, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, M_new, M_term, - GrB_NULL)); - GRB_TRY(GrB_eWiseAdd(M_new, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, M_new, M_nonterm, - GrB_NULL)); - GRB_TRY(GrB_eWiseAdd(M_new, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, M_new, M_call, - GrB_NULL)); - GRB_TRY(GrB_eWiseAdd(M_new, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, M_new, M_return, - GrB_NULL)); + GRB_TRY(GrB_eWiseAdd(M_new, GrB_NULL, GrB_LOR, GrB_LOR, M_new, M_term, GrB_NULL)); + GRB_TRY(GrB_eWiseAdd(M_new, GrB_NULL, GrB_LOR, GrB_LOR, M_new, M_nonterm, GrB_NULL)); + GRB_TRY(GrB_eWiseAdd(M_new, GrB_NULL, GrB_LOR, GrB_LOR, M_new, M_call, GrB_NULL)); + GRB_TRY(GrB_eWiseAdd(M_new, GrB_NULL, GrB_LOR, GrB_LOR, M_new, M_return, GrB_NULL)); // M_new &= ~P GRB_TRY(GrB_assign(M_new, P, GrB_NULL, M_new, GrB_ALL, Q, GrB_ALL, VV, GrB_DESC_RSC)); // P |= M_new - GRB_TRY(GrB_eWiseAdd(P, GrB_NULL, GrB_LOR, GrB_LOR_LAND_SEMIRING_BOOL, P, M_new, GrB_NULL)); + GRB_TRY(GrB_eWiseAdd(P, GrB_NULL, GrB_LOR, GrB_LOR, P, M_new, GrB_NULL)); GrB_Matrix swap = M; M = M_new; diff --git a/experimental/test/test_CFPQ_RSM.c b/experimental/test/test_CFPQ_RSM.c index 2ac746ead4..6a37fff856 100644 --- a/experimental/test/test_CFPQ_RSM.c +++ b/experimental/test/test_CFPQ_RSM.c @@ -44,6 +44,15 @@ static RSM *alloc_rsm(GrB_Index state_count, GrB_Index term_count, GrB_Index non LAGraph_Calloc((void **)&r->start_states, nonterm_count, sizeof(GrB_Index), msg); LAGraph_Calloc((void **)&r->final_states, nonterm_count, sizeof(GrB_Vector), msg); + for (size_t i = 0; i < term_count; ++i) + GrB_Matrix_new(&r->terminal_matrices[i], GrB_BOOL, state_count, state_count); + + for (size_t i = 0; i < nonterm_count; ++i) + { + GrB_Matrix_new(&r->nonterminal_matrices[i], GrB_BOOL, state_count, state_count); + GrB_Vector_new(&r->final_states[i], GrB_BOOL, state_count); + } + return r; } @@ -154,11 +163,6 @@ static void init_rsm_aSb_ab(void) GrB_Index Q = 4; rsm = alloc_rsm(Q, /*term_count=*/2, /*nonterm_count=*/1, /*start_nonterm=*/0); - OK(GrB_Matrix_new(&rsm->terminal_matrices[0], GrB_BOOL, Q, Q)); - OK(GrB_Matrix_new(&rsm->terminal_matrices[1], GrB_BOOL, Q, Q)); - OK(GrB_Matrix_new(&rsm->nonterminal_matrices[0], GrB_BOOL, Q, Q)); - OK(GrB_Vector_new(&rsm->final_states[0], GrB_BOOL, Q)); - OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[0], true, 0, 1)); // 0 -a-> 1 OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[1], true, 1, 3)); // 1 -b-> 3 OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[1], true, 2, 3)); // 2 -b-> 3 @@ -173,13 +177,6 @@ static void init_rsm_aSb_c(void) GrB_Index Q = 4; rsm = alloc_rsm(Q, /*term_count=*/3, /*nonterm_count=*/1, /*start_nonterm=*/0); - OK(GrB_Matrix_new(&rsm->terminal_matrices[0], GrB_BOOL, Q, Q)); - OK(GrB_Matrix_new(&rsm->terminal_matrices[1], GrB_BOOL, Q, Q)); - OK(GrB_Matrix_new(&rsm->terminal_matrices[2], GrB_BOOL, Q, Q)); - OK(GrB_Matrix_new(&rsm->nonterminal_matrices[0], GrB_BOOL, Q, Q)); - - OK(GrB_Vector_new(&rsm->final_states[0], GrB_BOOL, Q)); - OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[0], true, 0, 1)); // 0 -a-> 1 OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[2], true, 0, 3)); // 0 -c-> 3 OK(GrB_Matrix_setElement_BOOL(rsm->terminal_matrices[1], true, 2, 3)); // 2 -b-> 3 @@ -195,13 +192,6 @@ static void init_rsm_complex(void) { GrB_Index Q = 11; rsm = alloc_rsm(Q, /*term_count=*/4, /*nonterm_count=*/3, /*start_nonterm=*/0); - for (int t = 0; t < 4; ++t) - OK(GrB_Matrix_new(&rsm->terminal_matrices[t], GrB_BOOL, Q, Q)); - for (int n = 0; n < 3; ++n) - { - OK(GrB_Matrix_new(&rsm->nonterminal_matrices[n], GrB_BOOL, Q, Q)); - OK(GrB_Vector_new(&rsm->final_states[n], GrB_BOOL, Q)); - } // --- S Component (States 0-2) --- // S -> X | Y | S X | S Y