diff --git a/experimental/algorithm/LAGraph_CFPQ_RSM.c b/experimental/algorithm/LAGraph_CFPQ_RSM.c new file mode 100644 index 0000000000..2f47f58caa --- /dev/null +++ b/experimental/algorithm/LAGraph_CFPQ_RSM.c @@ -0,0 +1,411 @@ +#include "LAGraph.h" +#include "LG_internal.h" +#include + +#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); \ + GrB_free(&source_mask); \ + if (Stack != NULL) \ + { \ + for (size_t _i = 0; _i < num_nonterm; ++_i) \ + GrB_free(&Stack[_i]); \ + LAGraph_Free((void **)&Stack, msg); \ + } \ + if (graph_nt != NULL) \ + { \ + for (size_t _i = 0; _i < num_nonterm; ++_i) \ + GrB_free(&graph_nt[_i]); \ + 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_nonterm; ++_i) \ + GrB_free(&rsm_call[_i]); \ + LAGraph_Free((void **)&rsm_call, msg); \ + } \ + } +#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; +} + +static GrB_Info build_call_matrix(GrB_Matrix *call, GrB_Matrix rsm_nt, GrB_Vector rsm_start, + GrB_Index Q) +{ + GrB_Vector call_mask = GrB_NULL; + + H_TRY(GrB_Matrix_new(call, GrB_BOOL, Q, Q)); + H_TRY(GrB_Vector_new(&call_mask, GrB_BOOL, Q)); + + // [q_call, q_ret] -> [q_call] + 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)call_mask, + (GrB_Matrix)rsm_start, GrB_DESC_T1)); + + H_TRY(GrB_Vector_free(&call_mask)); + + return GrB_SUCCESS; +} + +#undef H_TRY +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, const RSM *rsm, const GrB_Matrix *graph_term, + const GrB_Index *sources, size_t num_sources, 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; + 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_Vector source_mask = GrB_NULL; + + 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(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(num_sources <= V, GrB_INVALID_VALUE); + + *reachable = GrB_NULL; + + // allocate per-nonterminal arrays + 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_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_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)); + } + + // 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(&source_mask, GrB_BOOL, V)); + GRB_TRY(GrB_Vector_new(reachable, GrB_BOOL, V)); + + // seed initial frontier M and Stack + { + + GrB_Index nstart = 0; + GRB_TRY(GrB_Vector_nvals(&nstart, rsm_start[start_nonterm])); + + 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(seed_states, seed_vals, &nstart, + rsm_start[start_nonterm]); + LAGraph_Free((void **)&seed_vals, GrB_NULL); + + if (info != GrB_SUCCESS) + { + LAGraph_Free((void **)&seed_states, GrB_NULL); + LG_FREE_ALL; + return info; + } + + for (size_t i = 0; i < num_sources; ++i) + { + 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(P, true, q, seed_col)); + GRB_TRY(GrB_Matrix_setElement_BOOL(Stack[start_nonterm], true, q, seed_col)); + } + } + + LAGraph_Free((void **)&seed_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_term; ++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_nonterm; ++i) + { + 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)); + } + + // PHASE 3: CALL TRANSITIONS + GRB_TRY(GrB_Matrix_clear(M_call)); + for (size_t i = 0; i < num_nonterm; ++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_NULL, 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; + + // 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_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_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)); + + 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, 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 || rsm_nonterm[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)final_states[i], M, GrB_DESC_T0)); + + GRB_TRY(GxB_Matrix_reshape(new_edges, false, 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_eWiseAdd(graph_nt[i], GrB_NULL, GrB_LOR, GrB_LOR, graph_nt[i], new_edges, + GrB_NULL)); + + LG_TRY(s_mxm_chain(M_return, rsm_nonterm[i], P, new_edges, temp1, temp2, Q, V)); + + 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, 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, P, M_new, GrB_NULL)); + + GrB_Matrix swap = M; + M = M_new; + M_new = swap; + + GRB_TRY(GrB_Matrix_nvals(&m_nvals, M)); + } + + // 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 new file mode 100644 index 0000000000..6a37fff856 --- /dev/null +++ b/experimental/test/test_CFPQ_RSM.c @@ -0,0 +1,575 @@ +#include "LAGraph.h" +#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) +{ + LAGraph_Init(msg); +} + +static void teardown(void) +{ + 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); + + 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; +} + +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->final_states != 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_Index *expected, GrB_Index n_expected, 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; + free(idx); + free(val); + } + + // 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); +} + +//============================================================================== +// RSM builders +//============================================================================== +// Grammar 1: S -> a S b | a b +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_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 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_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 init_rsm_complex(void) +{ + GrB_Index Q = 11; + rsm = alloc_rsm(Q, /*term_count=*/4, /*nonterm_count=*/3, /*start_nonterm=*/0); + + // --- S Component (States 0-2) --- + // S -> X | Y | S X | S Y + 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(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(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(sources, num_sources) \ + LAGraph_CFPQ_RSM(&result, rsm, graph_term, (sources), (num_sources), V, msg); + +//============================================================================== +// 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(); + + init_rsm_aSb_ab(); + init_graph_double_cycle_ab(); + + GrB_Index sources = {1}; + GrB_Info info = run_algorithm(&sources, 1); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {0, 3}; + check_reachable(expected, /*n_expected=*/2, "TC1"); + } + + free_workspace(); + 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(); + init_rsm_aSb_ab(); + init_graph_self_loops(); + + GrB_Index sources = {0}; + GrB_Info info = run_algorithm(&sources, 1); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {1}; + check_reachable(expected, /*n_expected=*/1, "TC2"); + } + + free_workspace(); + 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(); + init_rsm_aSb_ab(); + init_graph_single_node(); + + GrB_Index sources = {0}; + GrB_Info info = run_algorithm(&sources, 1); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {0}; + check_reachable(expected, /*n_expected=*/1, "TC3"); + } + + free_workspace(); + 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(); + init_rsm_aSb_c(); + init_graph_loops_abc(); + + GrB_Index sources = {0}; + GrB_Info info = run_algorithm(&sources, 1); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {1}; + check_reachable(expected, 1, "TC4"); + } + + free_workspace(); + 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(); + init_rsm_aSb_c(); + init_graph_loops_abc_1(); + + GrB_Index sources = {1}; + GrB_Info info = run_algorithm(&sources, 1); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {3}; + check_reachable(expected, 1, "TC5"); + } + + free_workspace(); + 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(); + init_rsm_aSb_c(); + + 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)); + + 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)); + + GrB_Index sources = {1}; + GrB_Info info = run_algorithm(&sources, 1); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {2}; + check_reachable(expected, 1, "TC6"); + } + + free_workspace(); + OK(info); + teardown(); +} + +//------------------------------------------------------------------------------ +// TC7: S -> a S b | c +// 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_Index sources = {0}; + GrB_Info info = run_algorithm(&sources, 1); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {2, 3, 4}; + check_reachable(expected, 3, "TC7"); + } + + free_workspace(); + 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(); + init_rsm_aSb_c(); + init_graph_loops_abc_2(); + + GrB_Index sources = {1}; + GrB_Info info = run_algorithm(&sources, 1); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {2, 3, 4}; + check_reachable(expected, 3, "TC8"); + } + + free_workspace(); + OK(info); + teardown(); +} + +void test_TC9_mult_recursive_cycles(void) +{ + setup(); + + init_rsm_complex(); + init_graph_multi_cycle(); + + GrB_Index sources = {1}; + GrB_Info info = run_algorithm(&sources, 1); + if (info == GrB_SUCCESS) + { + GrB_Index expected[] = {2, 3, 5}; + check_reachable(expected, 3, "TC9"); + } + + free_workspace(); + 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}, + {"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}, + {NULL, NULL}}; diff --git a/include/LAGraphX.h b/include/LAGraphX.h index 3355addb2e..18171ed555 100644 --- a/include/LAGraphX.h +++ b/include/LAGraphX.h @@ -854,7 +854,8 @@ 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_VertexCentrality_Triangle // vertex triangle-centrality @@ -1096,6 +1097,21 @@ GrB_Info LAGraph_CFL_reachability char *msg // Message string for error reporting. ) ; +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, 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 //------------------------------------------------------------------------------