Skip to content

Implement RSMs reachability algorithm - #7

Open
Zestria wants to merge 25 commits into
SparseLinearAlgebra:stablefrom
Zestria:feat/rsm-reach
Open

Zestria wants to merge 25 commits into
SparseLinearAlgebra:stablefrom
Zestria:feat/rsm-reach

Conversation

@Zestria

@Zestria Zestria commented Dec 3, 2025

Copy link
Copy Markdown

This patch implements a context-free path query algorithm using recursive state machines (RSMs) over edge-labeled directed graphs. It extends the approach of the previously unified regular path search algorithm by using linear algebra operations on adjacency matrices, as presented in https://arxiv.org/abs/2412.10287.

Comment thread experimental/algorithm/LAGraph_RSM_reachability.c Outdated
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)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree that this line of code at least deserves a named #defined or a constant MAX_ITERATIONS at the beginning of the file :)

Comment thread experimental/test/test_RSM_reachability.c Outdated
@Zestria

This comment was marked as outdated.

Zestria added 6 commits April 17, 2026 22:48
Matrix reshape from |V|*|V| to 1*|V*V| triggered 32-bit indexing
overflow when V*V > 2^32. The library correctly detects potential
overflow (Cj_is_32=false), but the in-place reshape operation
(in_place=true) inherits T->j_is_32=true from the original matrix,
causing J_work_is_32 to become true and overflow to occur.

Solution: use matrix_reshape_dup instead of matrix_reshape to create
a new matrix (in_place=false), so J_work_is_32 correctly evaluates
to false and 64-bit indexing is used.
@Zestria

This comment was marked as outdated.

@Zestria

Zestria commented May 11, 2026 •

Copy link
Copy Markdown
Author

Algorithm:
algorithm.pdf

@Zestria

This comment was marked as outdated.

@Zestria

This comment was marked as outdated.

@Zestria

Zestria commented Sep 15, 2026 •

Copy link
Copy Markdown
Author

Experiment's design and result:
experiments.pdf

Accumulate W = R_A^T * P per nonterminal instead of recomputing it in Phase 4.
…matrix

Replace M_term/M_nonterm/M_return with a single |Q*V|*|V| accumulator, reshaped once per iteration instead of per product.
Comment on lines +239 to +298
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;
}

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nit: that's a huge amount of work for (almost) the same code! Can we add a #define macro that would, basically, 1) verify the graph, 2) extract matrices out of it, 70 lines of code is quite a few.

if (B_ret[i] == GrB_NULL || B_call[i] == GrB_NULL || B_S[i] == GrB_NULL)
continue;

GRB_TRY(GrB_Matrix_clear(M_ret));

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No need for a separate Matrix_clear call: just use GrB_DESC_RT0 (add R) on line 536. R means "replace".

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same for the next line, use GrB_DESC_RS on line 545 instead.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, I think it is a general recommendation for all the clear()'s in the main procedure.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants