diff --git a/experimental/algorithm/LAGraph_RSM_reachability.c b/experimental/algorithm/LAGraph_RSM_reachability.c new file mode 100644 index 0000000000..8972116a4c --- /dev/null +++ b/experimental/algorithm/LAGraph_RSM_reachability.c @@ -0,0 +1,637 @@ +//------------------------------------------------------------------------------ +// LAGraph_RSM_reachability.c +//------------------------------------------------------------------------------ + +#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_RSM_reachability.c b/experimental/test/test_RSM_reachability.c new file mode 100644 index 0000000000..bf10a389de --- /dev/null +++ b/experimental/test/test_RSM_reachability.c @@ -0,0 +1,333 @@ +#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 3355addb2e..daf5d5727d 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, + + // 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