diff --git a/experimental/algorithm/LAGraph_CFL_optimized_matrix_opt.c b/experimental/algorithm/LAGraph_CFL_optimized_matrix_opt.c new file mode 100644 index 0000000000..24ddded6d2 --- /dev/null +++ b/experimental/algorithm/LAGraph_CFL_optimized_matrix_opt.c @@ -0,0 +1,1074 @@ +//------------------------------------------------------------------------------ +// LAGraph_CFL_optimized_matrix_opt.c: Implementation of operations for +// Optimized Context-Free Language Reachability Matrix-Based Algorithm +//------------------------------------------------------------------------------ +// +// LAGraph, (c) 2019-2025 by The LAGraph Contributors, All Rights Reserved. +// SPDX-License-Identifier: BSD-2-Clause + +// Contributed by Vlasenco Daniel, Ilhom Kombaev, Semyon Grigoriev, St. Petersburg State +// University. + +//------------------------------------------------------------------------------ + +// Code is an implementation of optimized matrix operations for CFL_reachability +// algorithms, described in the following paper: +// * Ilia Muravev, "Optimization of the Context-Free Language Reachability Matrix-Based +// Algorithm" and based on the python implementation from: +// https://github.com/FormalLanguageConstrainedPathQuerying/CFPQ_PyAlgo/tree/murav/optimize-matrix + +#include "LG_internal.h" +#include + +#define OPT_EMPTY (1 << 0) +#define OPT_FORMAT (1 << 1) +#define OPT_LAZY (1 << 2) +#define OPT_BLOCK (1 << 3) + +typedef CFL_Matrix Matrix; +typedef enum CFL_Matrix_block Matrix_block; + +#define TO_COL(matrix) GrB_set(matrix, GrB_COLMAJOR, GrB_STORAGE_ORIENTATION_HINT) +#define TO_ROW(matrix) GrB_set(matrix, GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT) + +GrB_Info matrix_to_format(Matrix *matrix, int32_t format, bool is_both) { + // No changes required + if (matrix->format == format) { + return GrB_SUCCESS; + } + + // Matrix contains both formats, need to update required format and switch to it + if (matrix->is_both) { + if (matrix->format == GrB_ROWMAJOR && format == GrB_COLMAJOR) { + GrB_Matrix_assign(matrix->base_col, GrB_NULL, GrB_NULL, matrix->base, GrB_ALL, + matrix->nrows, GrB_ALL, matrix->ncols, GrB_NULL); + TO_COL(matrix->base_col); + matrix->base = matrix->base_col; + } else if (matrix->format == GrB_COLMAJOR && format == GrB_ROWMAJOR) { + GrB_Matrix_assign(matrix->base_row, GrB_NULL, GrB_NULL, matrix->base, GrB_ALL, + matrix->nrows, GrB_ALL, matrix->ncols, GrB_NULL); + TO_ROW(matrix->base_row); + matrix->base = matrix->base_row; + } + matrix->format = format; + return GrB_SUCCESS; + } + + // Matrix contains just one matrix and format is not same + GrB_Matrix *new_matrix = + matrix->format == GrB_ROWMAJOR ? &matrix->base_col : &matrix->base_row; + GrB_Matrix *old_matrix = + matrix->format == GrB_ROWMAJOR ? &matrix->base_row : &matrix->base_col; + + if (is_both) { + GrB_Matrix_new(new_matrix, GrB_BOOL, matrix->nrows, matrix->ncols); + GrB_Matrix_assign(*new_matrix, GrB_NULL, GrB_NULL, *old_matrix, GrB_ALL, + matrix->nrows, GrB_ALL, matrix->ncols, GrB_NULL); + matrix->is_both = true; + } else { + *new_matrix = *old_matrix; + *old_matrix = NULL; + } + + matrix->base = format == GrB_ROWMAJOR ? matrix->base_row : matrix->base_col; + format == GrB_ROWMAJOR ? TO_ROW(matrix->base) : TO_COL(matrix->base); + matrix->format = format; + return GrB_SUCCESS; +} + +// clear methods + +GrB_Info matrix_clear(Matrix *A) { + GrB_Info result = GrB_Matrix_clear(A->base); + result; + CFL_matrix_update(A); + return result; +} + +GrB_Info matrix_clear_format(Matrix *A, int8_t optimizations) { + if (!(optimizations & OPT_FORMAT)) { + return matrix_clear(A); + } + + if (!A->is_both) { + return matrix_clear(A); + } + + matrix_to_format(A, GrB_ROWMAJOR, false); + GrB_Info result = matrix_clear(A); + + if (result < GrB_SUCCESS) { + return result; + } + + matrix_to_format(A, GrB_COLMAJOR, false); + return matrix_clear(A); +} + +GrB_Info matrix_clear_empty(Matrix *A, int8_t optimizations) { + if (!(optimizations & OPT_EMPTY)) { + return matrix_clear_format(A, optimizations); + } + + CFL_matrix_update(A); + if (A->nvals == 0) { + return GrB_SUCCESS; + } + + return matrix_clear_format(A, optimizations); +} + +// duplicate methods + +GrB_Info matrix_dup(Matrix *output, Matrix *input) { + if (output == input) { + return GrB_SUCCESS; + } + + GrB_Info result = GrB_Matrix_apply(output->base, GrB_NULL, GrB_NULL, + GrB_IDENTITY_BOOL, input->base, GrB_NULL); + CFL_matrix_update(output); + + return result; +} + +GrB_Info matrix_dup_format(Matrix *output, Matrix *input, int8_t optimizations) { + if (!(optimizations & OPT_FORMAT)) { + return matrix_dup(output, input); + } + + if (!output->is_both) { + CFL_matrix_update(output); + CFL_matrix_update(input); + Matrix *larger = output->nvals > input->nvals ? output : input; + + matrix_to_format(output, larger->format, false); + matrix_to_format(input, larger->format, false); + + return matrix_dup(output, input); + } + + matrix_to_format(output, GrB_ROWMAJOR, false); + GrB_Info result = matrix_dup(output, input); + + if (result < GrB_SUCCESS) { + return result; + } + + matrix_to_format(output, GrB_COLMAJOR, false); + return matrix_dup(output, input); +} + +GrB_Info matrix_dup_empty(Matrix *output, Matrix *input, int8_t optimizations) { + if (!(optimizations & OPT_EMPTY)) { + return matrix_dup_format(output, input, optimizations); + } + + CFL_matrix_update(input); + if (input->nvals == 0) { + return matrix_clear_empty(output, optimizations); + } + + return matrix_dup_format(output, input, optimizations); +} + +GrB_Info matrix_dup_lazy(Matrix *output, Matrix *input, int8_t optimizations) { + if (!(optimizations & OPT_LAZY)) { + return matrix_dup_empty(output, input, optimizations); + } + + if (!input->is_lazy) { + return matrix_dup_empty(output, input, optimizations); + } + + bool result = true; + for (size_t i = 0; i < input->base_matrices_count; i++) { + result &= matrix_dup_empty(&output->base_matrices[i], &input->base_matrices[i], + optimizations); + } + output->base_matrices_count = input->base_matrices_count; + + return result; +} + +GrB_Info block_matrix_hyper_rotate_i(Matrix *matrix, enum CFL_Matrix_block format); + +GrB_Info matrix_dup_block(Matrix *output, Matrix *input, int8_t optimizations) { + if (!(optimizations & OPT_BLOCK)) { + return matrix_dup_lazy(output, input, optimizations); + } + + if (output->block_type == CELL && input->block_type == CELL) { + return matrix_dup_lazy(output, input, optimizations); + } + + block_matrix_hyper_rotate_i(input, output->block_type); + return matrix_dup_lazy(output, input, optimizations); +} + +// block optimization specific methods + +GrB_Info block_matrix_hyper_rotate_i(Matrix *matrix, enum CFL_Matrix_block format) { + if (matrix->is_lazy) { + for (size_t i = 0; i < matrix->base_matrices_count; i++) { + block_matrix_hyper_rotate_i(&matrix->base_matrices[i], format); + } + + CFL_matrix_update(matrix); + return GrB_SUCCESS; + } + + if (matrix->block_type == CELL) { + return GrB_SUCCESS; + } + + if (matrix->block_type == format) { + return GrB_SUCCESS; + } + + GrB_Scalar scalar_true; + GrB_Scalar_new(&scalar_true, GrB_BOOL); + GrB_Scalar_setElement_BOOL(scalar_true, true); + + CFL_matrix_update(matrix); + + if (matrix->block_type == VEC_VERT) { + // fix: change to lagraph malloc + GrB_Index *nrows = malloc(matrix->nvals * sizeof(GrB_Index)); + GrB_Index *ncols = malloc(matrix->nvals * sizeof(GrB_Index)); + + GrB_Matrix_extractTuples_BOOL(nrows, ncols, NULL, &matrix->nvals, matrix->base); + + for (size_t i = 0; i < matrix->nvals; i++) { + ncols[i] = ncols[i] + nrows[i] / matrix->ncols * matrix->ncols; + nrows[i] = nrows[i] % matrix->ncols; + } + + GrB_Matrix new; + GrB_Matrix_new(&new, GrB_BOOL, matrix->ncols, matrix->nrows); + GxB_Matrix_build_Scalar(new, nrows, ncols, scalar_true, matrix->nvals); + CFL_matrix_free(matrix); + *matrix = + matrix->is_lazy ? CFL_matrix_from_base_lazy(new) : CFL_matrix_from_base(new); + free(nrows); + free(ncols); + GrB_free(&scalar_true); + return GrB_SUCCESS; + } + + GrB_Index *nrows = malloc(matrix->nvals * sizeof(GrB_Index)); + GrB_Index *ncols = malloc(matrix->nvals * sizeof(GrB_Index)); + + GrB_Matrix_extractTuples_BOOL(nrows, ncols, NULL, &matrix->nvals, matrix->base); + + for (size_t i = 0; i < matrix->nvals; i++) { + nrows[i] = nrows[i] + ncols[i] / matrix->nrows * matrix->nrows; + ncols[i] = ncols[i] % matrix->nrows; + } + + GrB_Matrix new; + GrB_Matrix_new(&new, GrB_BOOL, matrix->ncols, matrix->nrows); + GxB_Matrix_build_Scalar(new, nrows, ncols, scalar_true, matrix->nvals); + CFL_matrix_free(matrix); + *matrix = + matrix->is_lazy ? CFL_matrix_from_base_lazy(new) : CFL_matrix_from_base(new); + free(nrows); + free(ncols); + GrB_free(&scalar_true); + return GrB_SUCCESS; +} + +void block_matrix_to_diag(Matrix *diag, Matrix *input) { + if (input->block_type == CELL) { + return; + } + + GrB_Scalar scalar_true; + GrB_Scalar_new(&scalar_true, GrB_BOOL); + GrB_Scalar_setElement_BOOL(scalar_true, true); + + CFL_matrix_update(input); + + GrB_Index *rows = malloc(input->nvals * sizeof(GrB_Index)); + GrB_Index *cols = malloc(input->nvals * sizeof(GrB_Index)); + GrB_Matrix_extractTuples_BOOL(rows, cols, NULL, &input->nvals, input->base); + + if (input->block_type == VEC_HORIZ) { + for (size_t i = 0; i < input->nvals; i++) { + rows[i] = rows[i] + cols[i] / input->nrows * input->nrows; + } + } + + if (input->block_type == VEC_VERT) { + for (size_t i = 0; i < input->nvals; i++) { + cols[i] = cols[i] + rows[i] / input->ncols * input->ncols; + } + } + + GxB_Matrix_build_Scalar(diag->base, rows, cols, scalar_true, input->nvals); + + free(rows); + free(cols); + GrB_free(&scalar_true); +} + +void block_matrix_reduce(Matrix *matrix, Matrix *input, int8_t optimizations) { + if (input->block_type == CELL) { + matrix_dup_block(matrix, input, optimizations); + } + + GrB_Scalar scalar_true; + GrB_Scalar_new(&scalar_true, GrB_BOOL); + GrB_Scalar_setElement_BOOL(scalar_true, true); + + CFL_matrix_update(input); + + GrB_Index *rows = malloc(input->nvals * sizeof(GrB_Index)); + GrB_Index *cols = malloc(input->nvals * sizeof(GrB_Index)); + GrB_Matrix_extractTuples_BOOL(rows, cols, NULL, &input->nvals, input->base); + + if (input->block_type == VEC_VERT) { + for (size_t i = 0; i < input->nvals; i++) { + rows[i] = rows[i] % input->ncols; + } + } + + if (input->block_type == VEC_HORIZ) { + for (size_t i = 0; i < input->nvals; i++) { + cols[i] = cols[i] % input->nrows; + } + } + + GxB_Matrix_build_Scalar(matrix->base, rows, cols, scalar_true, input->nvals); + CFL_matrix_update(matrix); + + free(rows); + free(cols); + GrB_free(&scalar_true); +} + +void block_matrix_repeat_into_vector(Matrix *matrix, Matrix *input, + GrB_Index block_count) { + GrB_Matrix *tiles = malloc(block_count * sizeof(GrB_Matrix)); + for (size_t i = 0; i < block_count; i++) { + tiles[i] = input->base; + } + + GxB_Matrix_concat(matrix->base, tiles, block_count, 1, GrB_NULL); + CFL_matrix_update(matrix); + free(tiles); +} + +// lazy optimization specific methods + +GrB_Info matrix_wise_empty(Matrix *output, Matrix *first, Matrix *second, bool accum, + int8_t optimizations); + +GrB_Info matrix_sort_lazy(Matrix *A, bool reverse) { + for (size_t i = 0; i < A->base_matrices_count; i++) { + for (size_t j = i + 1; j < A->base_matrices_count; j++) { + Matrix *first = reverse ? &A->base_matrices[i] : &A->base_matrices[j]; + Matrix *second = reverse ? &A->base_matrices[j] : &A->base_matrices[i]; + CFL_matrix_update(first); + CFL_matrix_update(second); + + if (first->nvals < second->nvals) { + Matrix temp = A->base_matrices[i]; + A->base_matrices[i] = A->base_matrices[j]; + A->base_matrices[j] = temp; + } + } + } + + return GrB_SUCCESS; +} + +Matrix CFL_matrix_to_base(Matrix *input, int8_t optimizations) { + if (!input->is_lazy) { + Matrix result = CFL_matrix_create(input->nrows, input->ncols); + CFL_dup(&result, input, optimizations); + return result; + } + + Matrix acc = CFL_matrix_create(input->nrows, input->ncols); + + GrB_Matrix _matrix; + GrB_Matrix_new(&_matrix, GrB_BOOL, input->nrows, input->ncols); + Matrix matrix = CFL_matrix_from_base_lazy(_matrix); + GrB_free(&_matrix); + matrix.base = NULL; + + for (size_t i = 0; i < input->base_matrices_count; i++) { + Matrix base = CFL_matrix_create(input->nrows, input->ncols); + matrix.base_matrices[i] = base; + } + + CFL_dup(&matrix, input, optimizations); + + matrix_sort_lazy(&matrix, false); + for (size_t j = 0; j < matrix.base_matrices_count; j++) { + CFL_wise(&acc, &acc, &matrix.base_matrices[j], false, optimizations); + } + + CFL_matrix_free(&matrix); + + return acc; +} + +GrB_Info matrix_combine_lazy(Matrix *A, size_t threshold, int8_t optimizations) { + Matrix *new_matrices = malloc(sizeof(Matrix) * 50); + size_t new_size = 0; + + matrix_sort_lazy(A, false); + + for (size_t i = 0; i < A->base_matrices_count; i++) { + CFL_matrix_update(&A->base_matrices[i]); + if (new_size == 0 || A->base_matrices[i].nvals > threshold) { + new_matrices[new_size++] = A->base_matrices[i]; + continue; + } + + matrix_wise_empty(&new_matrices[new_size - 1], &new_matrices[new_size - 1], + &A->base_matrices[i], false, optimizations); + GrB_free(&A->base_matrices[i].base); + } + + free(A->base_matrices); + A->base_matrices = new_matrices; + A->base_matrices_count = new_size; + CFL_matrix_update(A); + + return GrB_SUCCESS; +} + +// create and update methods +void CFL_matrix_update(Matrix *matrix) { + if (!matrix->is_lazy) { + GrB_Matrix_nvals(&matrix->nvals, matrix->base); + } else { + size_t new_nnz = 0; + for (size_t i = 0; i < matrix->base_matrices_count; i++) { + GrB_Matrix_nvals(&matrix->base_matrices[i].nvals, + matrix->base_matrices[i].base); + new_nnz += matrix->base_matrices[i].nvals; + } + + matrix->nvals = new_nnz; + } + + if (!matrix->is_lazy) { + GrB_Matrix_nrows(&matrix->nrows, matrix->base); + GrB_Matrix_ncols(&matrix->ncols, matrix->base); + } else { + GrB_Matrix_nrows(&matrix->nrows, matrix->base_matrices[0].base); + GrB_Matrix_ncols(&matrix->ncols, matrix->base_matrices[0].base); + } + + if (matrix->nrows == matrix->ncols) + matrix->block_type = CELL; + else + matrix->block_type = matrix->nrows > matrix->ncols ? VEC_VERT : VEC_HORIZ; +} + +// Create optimizied matrix from base GrB_Matrix +Matrix CFL_matrix_from_base(GrB_Matrix matrix) { + Matrix result; + + result.base = matrix; + result.nvals = 0; // We will get actual info in update function + result.nrows = 0; + result.ncols = 0; + + // Format optimization fields + result.base_row = matrix; + result.base_col = NULL; + result.format = GrB_ROWMAJOR; + result.is_both = false; + + // Lazy addition optimization fields + result.is_lazy = false; + result.base_matrices = NULL; + result.base_matrices_count = 0; + + // Block optimization fields + result.block_type = CELL; + + CFL_matrix_update(&result); + return result; +} + +Matrix CFL_matrix_from_base_lazy(GrB_Matrix matrix) { + Matrix result = CFL_matrix_from_base(matrix); + + Matrix lazy_result = CFL_matrix_from_base(matrix); + lazy_result.is_lazy = true; + lazy_result.base_matrices = malloc(sizeof(CFL_Matrix) * 40); // TODO: dynamic size + lazy_result.base_matrices[0] = result; + lazy_result.base_matrices_count = 1; + CFL_matrix_update(&lazy_result); + + return lazy_result; +} + +Matrix CFL_matrix_create(GrB_Index nrows, GrB_Index ncols) { + GrB_Matrix _result; + GrB_Matrix_new(&_result, GrB_BOOL, nrows, ncols); + Matrix result = CFL_matrix_from_base(_result); + + return result; +} + +// TODO: free all base_matrices, free format matrices +void CFL_matrix_free(Matrix *matrix) { + if (matrix->is_lazy) { + for (size_t i = 0; i < matrix->base_matrices_count; i++) { + CFL_matrix_free(&matrix->base_matrices[i]); + } + + free(matrix->base_matrices); + matrix->base_matrices = NULL; + return; + } + + if (matrix->is_both) { + // matrix->base is base_row or base_col + GrB_Matrix_free(&matrix->base_row); + GrB_Matrix_free(&matrix->base_col); + return; + } + + GrB_Matrix_free(&matrix->base); +} + +// mxm operations + +GrB_Info matrix_mxm(Matrix *output, Matrix *first, Matrix *second, bool accum, + bool swap) { + Matrix *left = swap ? second : first; + Matrix *right = swap ? first : second; + + GrB_Info result = GrB_mxm(output->base, GrB_NULL, accum ? GxB_ANY_BOOL : GrB_NULL, + GxB_ANY_PAIR_BOOL, left->base, right->base, GrB_NULL); + // IS_ISO(output->base, "MXM output"); + CFL_matrix_update(output); + return result; +} + +GrB_Info matrix_mxm_format(Matrix *output, Matrix *first, Matrix *second, bool accum, + bool swap, int8_t optimizations) { + if (!(optimizations & OPT_FORMAT)) { + return matrix_mxm(output, first, second, accum, swap); + } + + CFL_matrix_update(first); + CFL_matrix_update(second); + GrB_Index left_nvals = swap ? second->nvals : first->nvals; + GrB_Index right_nvals = swap ? first->nvals : second->nvals; + + int32_t desired_orientation = left_nvals <= right_nvals ? GrB_ROWMAJOR : GrB_COLMAJOR; + + if (!first->is_both && first->format != desired_orientation && + !(first->nvals > second->nvals / 3.0)) { + GrB_Info result = matrix_mxm(output, first, second, accum, swap); + return result; + } + + matrix_to_format(first, desired_orientation, true); + matrix_to_format(second, desired_orientation, false); + matrix_to_format(output, desired_orientation, false); + GrB_Info result = matrix_mxm(output, first, second, accum, swap); + return result; +} + +GrB_Info matrix_mxm_empty(Matrix *output, Matrix *first, Matrix *second, bool accum, + bool swap, int8_t optimizations) { + if (!(optimizations & OPT_EMPTY)) { + return matrix_mxm_format(output, first, second, accum, swap, optimizations); + } + + CFL_matrix_update(first); + CFL_matrix_update(second); + + if (first->nvals == 0 || second->nvals == 0) { + if (accum) { + return GrB_SUCCESS; + } + + CFL_matrix_update(output); + if (output->nvals == 0) { + return GrB_SUCCESS; + } + + matrix_clear_empty(output, optimizations); + } + + return matrix_mxm_format(output, first, second, accum, swap, optimizations); +} + +GrB_Info matrix_mxm_lazy(Matrix *output, Matrix *first, Matrix *second, bool accum, + bool swap, int8_t optimizations) { + if (!(optimizations & OPT_LAZY)) { + return matrix_mxm_empty(output, first, second, accum, swap, optimizations); + } + + if (!first->is_lazy) { + return matrix_mxm_empty(output, first, second, accum, swap, optimizations); + } + + CFL_matrix_update(second); + + matrix_combine_lazy(first, second->nvals, optimizations); + matrix_sort_lazy(first, false); + + GrB_Matrix *accs = malloc(sizeof(GrB_Matrix) * first->base_matrices_count); + Matrix *acc_matrices = malloc(sizeof(Matrix) * first->base_matrices_count); + for (size_t i = 0; i < first->base_matrices_count; i++) { + GrB_Matrix_new(&accs[i], GrB_BOOL, swap ? second->nrows : first->nrows, + swap ? first->ncols : second->ncols); + acc_matrices[i] = CFL_matrix_from_base(accs[i]); + } + + for (size_t i = 0; i < first->base_matrices_count; i++) { + matrix_mxm_empty(&acc_matrices[i], &first->base_matrices[i], second, false, swap, + optimizations); + } + + GrB_Matrix acc; + GrB_Matrix_new(&acc, GrB_BOOL, swap ? second->nrows : first->nrows, + swap ? first->ncols : second->ncols); + Matrix acc_matrix = CFL_matrix_from_base(acc); + + for (size_t i = 0; i < first->base_matrices_count; i++) { + matrix_wise_empty(&acc_matrix, &acc_matrix, &acc_matrices[i], false, + optimizations); + GrB_free(&acc_matrices[i].base); + } + free(accs); + free(acc_matrices); + + if (accum) { + GrB_Info result = + matrix_wise_empty(output, output, &acc_matrix, false, optimizations); + CFL_matrix_free(&acc_matrix); + return result; + } + + GrB_Info result = matrix_dup_block(output, &acc_matrix, optimizations); + CFL_matrix_free(&acc_matrix); + + return result; +} + +GrB_Info matrix_wise_block(Matrix *output, Matrix *first, Matrix *second, bool accum, + int8_t optimizations); + +GrB_Info matrix_mxm_block(Matrix *output, Matrix *first, Matrix *second, bool accum, + bool swap, int8_t optimizations) { + if (!(optimizations & OPT_BLOCK)) { + return matrix_mxm_lazy(output, first, second, accum, swap, optimizations); + } + + if (first->block_type == CELL && second->block_type == CELL) { + return matrix_mxm_lazy(output, first, second, accum, swap, optimizations); + } + + if (first->block_type == CELL) { + block_matrix_hyper_rotate_i(second, swap ? VEC_VERT : VEC_HORIZ); + block_matrix_hyper_rotate_i(output, swap ? VEC_VERT : VEC_HORIZ); + + Matrix temp = CFL_matrix_create(swap ? second->nrows : first->nrows, + swap ? first->ncols : second->ncols); + matrix_mxm_lazy(&temp, first, second, accum, swap, optimizations); + matrix_wise_block(output, output, &temp, false, optimizations); + CFL_matrix_free(&temp); + + return GrB_SUCCESS; + } + + if (second->block_type == CELL) { + block_matrix_hyper_rotate_i(first, swap ? VEC_HORIZ : VEC_VERT); + block_matrix_hyper_rotate_i(output, swap ? VEC_HORIZ : VEC_VERT); + + Matrix temp = CFL_matrix_create(swap ? second->nrows : first->nrows, + swap ? first->ncols : first->ncols); + matrix_mxm_lazy(&temp, first, second, accum, swap, optimizations); + matrix_wise_block(output, output, &temp, false, optimizations); + CFL_matrix_free(&temp); + + return GrB_SUCCESS; + } + + GrB_Index size = first->nrows > first->ncols ? first->nrows : first->ncols; + GrB_Matrix _diag; + GrB_Matrix_new(&_diag, GrB_BOOL, size, size); + Matrix diag = CFL_matrix_from_base(_diag); + block_matrix_to_diag(&diag, second); + CFL_matrix_update(&diag); + + block_matrix_hyper_rotate_i(first, swap ? VEC_VERT : VEC_HORIZ); + block_matrix_hyper_rotate_i(output, swap ? VEC_VERT : VEC_HORIZ); + + Matrix temp = CFL_matrix_create(first->nrows, diag.ncols); + matrix_mxm_lazy(&temp, first, &diag, false, swap, optimizations); + GrB_Info result = matrix_wise_block(output, output, &temp, false, optimizations); + CFL_matrix_free(&temp); + CFL_matrix_free(&diag); + return result; +} + +// wise operations + +GrB_Info matrix_wise(Matrix *output, Matrix *first, Matrix *second, bool accum) { + GrB_BinaryOp accum_op = accum ? GxB_ANY_BOOL : GrB_NULL; + if (output == first) + accum_op = GrB_NULL; + + GrB_Info result = GrB_eWiseAdd(output->base, GrB_NULL, accum_op, GxB_ANY_BOOL, + first->base, second->base, GrB_NULL); + + CFL_matrix_update(output); + return result; +} + +GrB_Info matrix_wise_format(Matrix *output, Matrix *first, Matrix *second, bool accum, + int8_t optimizations) { + if (!(optimizations & OPT_FORMAT)) { + return matrix_wise(output, first, second, accum); + } + + if (!output->is_both) { + CFL_matrix_update(output); + CFL_matrix_update(first); + CFL_matrix_update(second); + + Matrix *larger = output->nvals > first->nvals ? output : first; + larger = larger->nvals > second->nvals ? larger : second; + + matrix_to_format(output, larger->format, false); + matrix_to_format(first, larger->format, false); + matrix_to_format(second, larger->format, false); + + return matrix_wise(output, first, second, accum); + } + + matrix_to_format(output, GrB_ROWMAJOR, false); + matrix_to_format(first, output->format, false); + matrix_to_format(second, output->format, false); + GrB_Info result = matrix_wise(output, first, second, accum); + + if (result < GrB_SUCCESS) { + return result; + } + + matrix_to_format(output, GrB_COLMAJOR, false); + matrix_to_format(first, output->format, false); + matrix_to_format(second, output->format, false); + return matrix_wise(output, first, second, accum); +} + +GrB_Info matrix_wise_empty(Matrix *output, Matrix *first, Matrix *second, bool accum, + int8_t optimizations) { + if (!(optimizations & OPT_EMPTY)) { + return matrix_wise_format(output, first, second, accum, optimizations); + } + + CFL_matrix_update(first); + CFL_matrix_update(second); + CFL_matrix_update(output); + + if (output == first) { + if (first->nvals == 0) { + return matrix_dup_empty(first, second, optimizations); + } + + if (second->nvals == 0) { + return GrB_SUCCESS; + } + + return matrix_wise_format(output, first, second, accum, optimizations); + } + + if (first->nvals == 0 && second->nvals == 0) { + if (accum || output->nvals == 0) { + return GrB_SUCCESS; + } + + return matrix_clear_empty(output, optimizations); + } + + if (first->nvals == 0) { + if (accum) { + return matrix_wise_empty(output, output, second, false, optimizations); + } + + return matrix_dup_empty(output, second, optimizations); + } + + if (second->nvals == 0) { + if (accum) { + return matrix_wise_empty(output, output, first, false, optimizations); + } + + return matrix_dup_empty(output, first, optimizations); + } + + return matrix_wise_format(output, first, second, accum, optimizations); +} + +GrB_Info matrix_wise_lazy(Matrix *output, Matrix *first, Matrix *second, bool accum, + int8_t optimizations) { + if (!(optimizations & OPT_LAZY)) { + return matrix_wise_empty(output, first, second, accum, optimizations); + } + + if (!first->is_lazy && !second->is_lazy) { + return matrix_wise_empty(output, first, second, accum, optimizations); + } + + if (!first->is_lazy && second->is_lazy) { + for (size_t i = 0; i < second->base_matrices_count; i++) { + matrix_wise_empty(output, first, &second->base_matrices[i], true, + optimizations); + } + + return GrB_SUCCESS; + } + + GrB_Matrix _other; + GrB_Matrix_new(&_other, GrB_BOOL, output->nrows, output->ncols); + Matrix other = CFL_matrix_from_base(_other); + matrix_dup_empty(&other, second, optimizations); + + size_t other_nvals = other.nvals >= 10 ? other.nvals : 10; + + while (true) { + bool found = false; + + for (size_t i = 0; i < first->base_matrices_count; i++) { + CFL_matrix_update(&first->base_matrices[i]); + size_t self_nvals = + first->base_matrices[i].nvals >= 10 ? first->base_matrices[i].nvals : 10; + + if (other_nvals / 10 <= self_nvals && self_nvals <= other_nvals * 10) { + matrix_wise_empty(&other, &other, &first->base_matrices[i], accum, + optimizations); + GrB_free(&first->base_matrices[i].base); + for (size_t j = i + 1; j < first->base_matrices_count; j++) { + first->base_matrices[j - 1] = first->base_matrices[j]; + } + first->base_matrices_count--; + found = true; + break; + } + } + + if (found) { + continue; + } + + first->base_matrices[first->base_matrices_count++] = other; + break; + } + + CFL_matrix_update(output); + + return GrB_SUCCESS; +} + +// - Any operation on two hyper vectors is performed block-wise +// - When hyper vector is added in-place to a cell, then sum of hyper vector's blocks is +// added to a cell +// - When cell is added in-place to a hyper vector, then cell is added to each of +// hyper vector's blocks +GrB_Info matrix_wise_block(Matrix *output, Matrix *first, Matrix *second, bool accum, + int8_t optimizations) { + if (!(optimizations & OPT_BLOCK)) { + return matrix_wise_lazy(output, first, second, accum, optimizations); + } + + if (output != first) { + // fprintf(stderr, "Matrix wise currently support only iadd operation"); + return GrB_INVALID_VALUE; + } + + if (first->block_type == CELL && second->block_type == CELL) { + return matrix_wise_lazy(output, first, second, accum, optimizations); + } + + // second is vector + if (first->block_type == CELL) { + Matrix temp_reduced = CFL_matrix_create(first->nrows, first->ncols); + block_matrix_reduce(&temp_reduced, second, optimizations); + + GrB_Info info = + matrix_wise_lazy(output, first, &temp_reduced, accum, optimizations); + CFL_matrix_free(&temp_reduced); + return info; + } + + // first is vector + if (second->block_type == CELL) { + // LG_SET_BURBLE(true); + Matrix temp_vector = CFL_matrix_create(first->nrows, first->ncols); + GrB_Index block_count = first->nrows > first->ncols ? first->nrows : first->ncols; + block_matrix_repeat_into_vector(&temp_vector, second, block_count); + + block_matrix_hyper_rotate_i(&temp_vector, first->block_type); + GrB_Info info = + matrix_wise_lazy(output, first, &temp_vector, accum, optimizations); + CFL_matrix_free(&temp_vector); + return info; + } + + // both are vector + block_matrix_hyper_rotate_i(second, first->block_type); + return matrix_wise_lazy(output, first, second, accum, optimizations); +} + +// rsub methods + +GrB_Info matrix_rsub(Matrix *output, Matrix *mask) { + GrB_Info result = GrB_eWiseAdd(output->base, mask->base, GrB_NULL, GxB_ANY_BOOL, + output->base, output->base, GrB_DESC_RSC); + + CFL_matrix_update(output); + return result; +} + +GrB_Info matrix_rsub_format(Matrix *output, Matrix *mask, int8_t optimizations) { + if (!(optimizations & OPT_FORMAT)) { + return matrix_rsub(output, mask); + } + + CFL_matrix_update(output); + CFL_matrix_update(mask); + + Matrix *larger_matrix = output->nvals > mask->nvals ? output : mask; + matrix_to_format(output, larger_matrix->format, false); + matrix_to_format(mask, larger_matrix->format, false); + + return matrix_rsub(output, mask); +} + +GrB_Info matrix_rsub_empty(Matrix *output, Matrix *mask, int8_t optimizations) { + if (!(optimizations & OPT_EMPTY)) { + return matrix_rsub_format(output, mask, optimizations); + } + + CFL_matrix_update(mask); + CFL_matrix_update(output); + if (mask->nvals == 0 || output->nvals == 0) { + return GrB_SUCCESS; + } + + return matrix_rsub_format(output, mask, optimizations); +} + +GrB_Info matrix_rsub_lazy(Matrix *output, Matrix *mask, int8_t optimizations) { + if (!(optimizations & OPT_LAZY)) { + return matrix_rsub_empty(output, mask, optimizations); + } + + if (!mask->is_lazy) { + return matrix_rsub_empty(output, mask, optimizations); + } + + CFL_matrix_update(output); + matrix_combine_lazy(mask, output->nvals, optimizations); + matrix_sort_lazy(mask, true); + + for (size_t i = 0; i < mask->base_matrices_count; i++) { + matrix_rsub_empty(output, &mask->base_matrices[i], optimizations); + } + + return GrB_SUCCESS; +} + +GrB_Info matrix_rsub_block(Matrix *output, Matrix *mask, int8_t optimizations) { + if (!(optimizations & OPT_BLOCK)) { + return matrix_rsub_lazy(output, mask, optimizations); + } + + if ((output->block_type == CELL && mask->block_type != CELL) || + (output->block_type != CELL && mask->block_type == CELL)) { + // fprintf(stderr, "Don't support rsub operation between cell and vector"); + return GrB_INVALID_VALUE; + } + + if (output->block_type == CELL) { + return matrix_rsub_lazy(output, mask, optimizations); + } + + block_matrix_hyper_rotate_i(output, mask->block_type); + return matrix_rsub_lazy(output, mask, optimizations); +} + +// utility methods + +// void matrix_print_lazy(Matrix *A) { +// return; +// GxB_Print_Level pr = 1; + +// if (!A->is_lazy) { +// CFL_matrix_update(A); +// GxB_print(A->base, pr); +// // printf("nnz: %ld\n", A->nvals); +// return; +// } + +// if (A->base_matrices_count == 1) { +// CFL_matrix_update(A); +// // printf("nnz: %ld\n", A->nvals); +// GxB_print(A->base_matrices[0].base, pr); +// return; +// } + +// GrB_Matrix _temp; +// GrB_Matrix_new(&_temp, GrB_BOOL, A->nrows, A->ncols); +// Matrix temp = CFL_matrix_from_base(_temp); +// for (size_t i = 0; i < A->base_matrices_count; i++) { +// matrix_wise_empty(&temp, &temp, &A->base_matrices[i], false, optimizations); +// } + +// A = &temp; +// GxB_print(A->base, pr); +// CFL_matrix_update(A); +// // printf("nnz: %ld\n", A->nvals); +// GrB_free(&_temp); +// } + +// void print_graph_info(Matrix *matrices, size_t count) { +// GrB_Index nnz = 0; + +// for (size_t i = 0; i < count; i++) { +// Matrix *A = &matrices[i]; +// CFL_matrix_update(A); +// nnz += A->nvals; +// } + +// printf("NNZ: %ld\n", nnz); +// } + +// order of optimizations: block -> lazy -> empty -> format + +GrB_Info CFL_mxm(Matrix *output, Matrix *first, Matrix *second, bool accum, bool swap, + int8_t optimizations) { + return matrix_mxm_block(output, first, second, accum, swap, optimizations); +} + +GrB_Info CFL_wise(Matrix *output, Matrix *first, Matrix *second, bool accum, + int8_t optimizations) { + return matrix_wise_block(output, first, second, accum, optimizations); +} + +GrB_Info CFL_rsub(Matrix *output, Matrix *mask, int8_t optimizations) { + return matrix_rsub_block(output, mask, optimizations); +} + +GrB_Info CFL_dup(Matrix *output, Matrix *input, int8_t optimizations) { + return matrix_dup_block(output, input, optimizations); +} + +GrB_Info CFL_clear(Matrix *A, int8_t optimizations) { + return matrix_clear_empty(A, optimizations); +} \ No newline at end of file diff --git a/experimental/algorithm/LAGraph_CFL_reachability_multsrc.c b/experimental/algorithm/LAGraph_CFL_reachability_multsrc.c new file mode 100644 index 0000000000..0b5f917436 --- /dev/null +++ b/experimental/algorithm/LAGraph_CFL_reachability_multsrc.c @@ -0,0 +1,429 @@ +//------------------------------------------------------------------------------ +// LAGraph_CFL_reachability_multsrc.c: Multiple-Source Context-Free Language Reachability +// Matrix-Based Algorithm +//------------------------------------------------------------------------------ +// +// LAGraph, (c) 2019-2025 by The LAGraph Contributors, All Rights Reserved. +// SPDX-License-Identifier: BSD-2-Clause + +// Contributed by Ilhom Kombaev, Vlasenco Daniel, Semyon Grigoriev, St. Petersburg State University. + +//------------------------------------------------------------------------------ + +// Code is based on the "multiple-source CFPQ", described in the following paper: +// * Arseny Terekhov et al., "Multiple-Source Context-Free Path Querying in Terms of Linear Algebra" +// * URL: https://openproceedings.org/2021/conf/edbt/p48.pdf + +// #define DEBUG_CFL_REACHABILITY + +#define LG_FREE_WORK \ + { \ + LAGraph_Free((void **) &nnzs_T, msg); \ + LAGraph_Free((void **) &nnzs_TSrc_B, msg); \ + LAGraph_Free((void **) &nnzs_TSrc_C, msg); \ + LAGraph_Free((void **) &ones_vec, msg); \ + LAGraph_Free((void **) &T, msg); \ + LAGraph_Free((void **) &TSrc, msg); \ + LAGraph_Free((void **) &MSrc, msg); \ + LAGraph_Free((void **) &identity_matrix, msg); \ + LAGraph_Free ((void **) &M, msg); \ + LAGraph_Free ((void **) &A, msg); \ + LAGraph_Free ((void **) &a, msg); \ + GrB_free(&true_scalar); \ + } + +#define LG_FREE_ALL \ + { \ + for (size_t i = 0; i < nonterms_count; i++) { \ + GrB_free(&T[i]); \ + GrB_free(&TSrc[i]); \ + } \ + \ + LG_FREE_WORK; \ + } + +#include "LG_internal.h" +#include + +#define ERROR_RULE(msg) \ + { \ + LG_ASSERT_MSGF(false, GrB_INVALID_VALUE, "Rule with index %ld is invalid. " msg, \ + i); \ + } + +#define ADD_TO_MSG(...) \ + { \ + if (msg_len == 0) { \ + msg_len += \ + snprintf(msg, LAGRAPH_MSG_LEN, \ + "LAGraph failure (file %s, line %d): ", __FILE__, __LINE__); \ + } \ + if (msg_len < LAGRAPH_MSG_LEN) { \ + msg_len += snprintf(msg + msg_len, LAGRAPH_MSG_LEN - msg_len, __VA_ARGS__); \ + } \ + } + +#define ADD_INDEX_TO_ERROR_RULE(rule, i) \ + { \ + rule.len_indices_str += snprintf(rule.indices_str + rule.len_indices_str, \ + LAGRAPH_MSG_LEN - rule.len_indices_str, \ + rule.count == 0 ? "%ld" : ", %ld", i); \ + rule.count++; \ + } + +#define PRINT_MATRIX(_m) { \ + for (size_t _i = 0; _i < n; _i++) { \ + for (size_t _j = 0; _j < n; _j++) { \ + if (GxB_Matrix_isStoredElement(_m, _i, _j) == GrB_SUCCESS) { \ + printf("1 "); \ + } \ + else { \ + printf("0 "); \ + } \ + } \ + printf("\n"); \ + } \ +} + +#define PRINT_VECTOR(_v) { \ + for (size_t _i = 0; _i < n; _i++) { \ + if (GxB_Vector_isStoredElement(_v, _i) == GrB_SUCCESS) { \ + printf("1 "); \ + } \ + else { \ + printf("0 "); \ + } \ + } \ + printf("\n"); \ +} + +#define PRINT_RULE(_r) { \ + printf("%c -> ", _r.nonterm + 'A' - 1); \ + if (_r.prod_A == -1 && _r.prod_B == -1) { \ + printf("eps\n"); \ + } else if (_r.prod_A != -1 && _r.prod_B == -1) { \ + printf("%c\n", _r.prod_A + 'a'); \ + } else { \ + printf("%c %c\n", _r.prod_A + 'A'- 1 ,_r.prod_B + 'A' - 1); \ + } \ +} + +// LAGraph_CFL_reachability_multsrc: Multiple-Source Context-Free Language Reachability +// Matrix-Based Algorithm +// +// This function determines the set of vertex pairs (u, v) in a graph (represented by +// adjacency matrices) such that u is in the given set of source vertices `src`, and +// there is a path from u to v, where the edge labels form a word from the language +// generated by the context-free grammar (represented by `rules`). + +GrB_Info LAGraph_CFL_reachability_multsrc +( + // Output + GrB_Matrix *output, // A handle for the matrix containing results. + // + // output: (i, j) = true if and only if there is a path + // from node i to node j whose edge labels form a word + // derivable from the specified CFG. + // Input + const GrB_Matrix *adj_matrices, // Array of adjacency matrices representing the graph. + // The length of this array is equal to the count of + // terminals (terms_count). + // + // adj_matrices[t]: (i, j) == 1 if and only if there + // is an edge between nodes i and j with the label of + // the terminal corresponding to index 't' (where t is + // in the range [0, terms_count - 1]). + GrB_Index *src, // Array of source vertices + // These parameters will become unsigned + int32_t src_count, // The total number of source vertices + int32_t terms_count, // The total number of terminal symbols in the CFG. + int32_t nonterms_count, // The total number of non-terminal symbols in the CFG. + const LAGraph_rule_WCNF *rules, // The rules of the CFG. + size_t rules_count, // The total number of rules in the CFG. + char *msg // Message string for error reporting. +) +{ + // Declare workspace and clear the msg string, if not NULL + GrB_Matrix *T; + GrB_Matrix *TSrc; + GrB_Matrix MSrc; + GrB_Matrix M; + GrB_Matrix A; + GrB_Matrix B; + GrB_Vector a; + GrB_Index n; // number of vertices in the graph + GrB_Matrix identity_matrix = NULL; + GrB_Index *nnzs_T = NULL; + GrB_Index *nnzs_TSrc_B = NULL; + GrB_Index *nnzs_TSrc_C = NULL; + LG_CLEAR_MSG; + size_t msg_len = 0; // For error formatting + bool iso_flag = false; + GrB_Scalar true_scalar; + GrB_Vector ones_vec; + + // Will change the interface and omit this check in the future + if (nonterms_count < 0 || terms_count < 0) + return GrB_INVALID_VALUE; + + if (!nonterms_count || !rules_count) + return GrB_INVALID_VALUE; + + bool t_empty_flags[nonterms_count]; // t_empty_flags[i] == true <=> T[i] is empty + bool t_src_empty_flags[nonterms_count]; // t_src_empty_flags[i] == true <=> TSrc[i] is empty + + if (!output || !rules || !adj_matrices || !src) + return GrB_NULL_POINTER; + if (src_count <= 0) + return GrB_INVALID_VALUE; + + // Find null adjacency matrices + bool found_null = false; + for (int32_t i = 0; i < terms_count; i++) { + if (adj_matrices[i] != NULL) + continue; + + if (!found_null) { + ADD_TO_MSG("Adjacency matrices with these indices are null: "); + ADD_TO_MSG("%d", i); + } else { + ADD_TO_MSG(", %d", i); + } + + found_null = true; + } + if (found_null) { + return GrB_NULL_POINTER; + } + + GRB_TRY(GrB_Matrix_ncols(&n, adj_matrices[0])); + if (n < src_count) return GrB_INVALID_VALUE; + + GrB_Scalar_new(&true_scalar, GrB_BOOL); + GrB_Scalar_setElement_BOOL(true_scalar, true); + + LG_TRY(LAGraph_Calloc((void **) &T, nonterms_count, sizeof(GrB_Matrix), msg)); + LG_TRY(LAGraph_Calloc((void **) &TSrc, nonterms_count, sizeof(GrB_Matrix), msg)); + + GRB_TRY(GrB_Vector_new(&ones_vec, GrB_BOOL, n)); + GRB_TRY(GrB_Vector_assign_BOOL(ones_vec, GrB_NULL, GrB_NULL, true, GrB_ALL, n, NULL)); + GRB_TRY(GrB_Matrix_diag(&identity_matrix, ones_vec, 0)); + + LG_TRY(LAGraph_Calloc((void **) &nnzs_T, nonterms_count, sizeof(GrB_Index), msg)); + LG_TRY(LAGraph_Calloc((void **) &nnzs_TSrc_B, nonterms_count, sizeof(GrB_Index), msg)); + LG_TRY(LAGraph_Calloc((void **) &nnzs_TSrc_C, nonterms_count, sizeof(GrB_Index), msg)); + + // Create nonterms matrices + for (int32_t i = 0; i < nonterms_count; i++) { + GRB_TRY(GrB_Matrix_new(&T[i], GrB_BOOL, n, n)); + GRB_TRY(GrB_Matrix_new(&TSrc[i], GrB_BOOL, n, n)); + t_empty_flags[i] = true; + t_src_empty_flags[i] = true; + } + + for (int32_t i = 0; i < src_count; i++) { + GrB_Matrix_setElement(TSrc[0], true, src[i], src[i]); + } + + t_src_empty_flags[0] = false; + + GRB_TRY(GrB_Matrix_dup(&MSrc, TSrc[0])); + GRB_TRY(GrB_Matrix_new(&A, GrB_BOOL, n, n)); + GRB_TRY(GrB_Matrix_new(&B, GrB_BOOL, n, n)); + GRB_TRY(GrB_Vector_new(&a, GrB_BOOL, n)); + GRB_TRY(GrB_Matrix_new(&M, GrB_BOOL, n, n)); + + // #ifdef DEBUG_CFL_REACHABILITY + // printf("MSrc:\n"); + // PRINT_MATRIX(MSrc); + // #endif + + // Arrays for processing rules + size_t eps_rules[rules_count], eps_rules_count = 0; // [Variable -> eps] + size_t term_rules[rules_count], term_rules_count = 0; // [Variable -> term] + size_t bin_rules[rules_count], bin_rules_count = 0; // [Variable -> AB] + + // Process rules + typedef struct { + size_t count; + size_t len_indices_str; + char indices_str[LAGRAPH_MSG_LEN]; + } rule_error_s; + + rule_error_s term_err = {0}; + rule_error_s nonterm_err = {0}; + rule_error_s invalid_err = {0}; + + for (size_t i = 0; i < rules_count; i++) { + LAGraph_rule_WCNF rule = rules[i]; + + bool is_rule_eps = rule.prod_A == -1 && rule.prod_B == -1; + bool is_rule_term = rule.prod_A != -1 && rule.prod_B == -1; + bool is_rule_bin = rule.prod_A != -1 && rule.prod_B != -1; + + // Check that all rules are well-formed + if (rule.nonterm < 0 || rule.nonterm >= nonterms_count) { + ADD_INDEX_TO_ERROR_RULE(nonterm_err, i); + } + + // [Variable -> eps] + if (is_rule_eps) { + eps_rules[eps_rules_count++] = i; + + continue; + } + + // [Variable -> term] + if (is_rule_term) { + term_rules[term_rules_count++] = i; + + if (rule.prod_A < -1 || rule.prod_A >= terms_count) { + ADD_INDEX_TO_ERROR_RULE(term_err, i); + } + + continue; + } + + // [Variable -> A B] + if (is_rule_bin) { + bin_rules[bin_rules_count++] = i; + + if (rule.prod_A < -1 || rule.prod_A >= nonterms_count || rule.prod_B < -1 || + rule.prod_B >= nonterms_count) { + ADD_INDEX_TO_ERROR_RULE(nonterm_err, i); + } + + continue; + } + + // [Variable -> _ B] + ADD_INDEX_TO_ERROR_RULE(invalid_err, i); + } + + if (term_err.count + nonterm_err.count + invalid_err.count > 0) { + ADD_TO_MSG("Count of invalid rules: %ld.\n", + term_err.count + nonterm_err.count + invalid_err.count); + + if (nonterm_err.count > 0) { + ADD_TO_MSG("Non-terminals must be in range [0, nonterms_count). "); + ADD_TO_MSG("indices of invalid rules: %s\n", nonterm_err.indices_str) + } + if (term_err.count > 0) { + ADD_TO_MSG("Terminals must be in range [-1, nonterms_count). "); + ADD_TO_MSG("indices of invalid rules: %s\n", term_err.indices_str) + } + if (invalid_err.count > 0) { + ADD_TO_MSG("[Variable -> _ B] type of rule is not acceptable. "); + ADD_TO_MSG("indices of invalid rules: %.120s\n", invalid_err.indices_str) + } + + LG_FREE_ALL; + return GrB_INVALID_VALUE; + } + + // Rule [Variable -> term] + for (size_t i = 0; i < term_rules_count; i++) { + LAGraph_rule_WCNF term_rule = rules[term_rules[i]]; + + GxB_eWiseUnion( + T[term_rule.nonterm], GrB_NULL, GrB_NULL, GxB_PAIR_BOOL, + T[term_rule.nonterm], true_scalar, adj_matrices[term_rule.prod_A], true_scalar, GrB_NULL + ); + + t_empty_flags[term_rule.nonterm] = false; + + // #ifdef DEBUG_CFL_REACHABILITY + // GxB_Matrix_iso(&iso_flag, T[term_rule.nonterm]); + // printf("[TERM] eWiseUnion: NONTERM: %d (ISO: %d)\n", term_rule.nonterm, iso_flag); + // #endif + } + + // Rule [Variable -> eps] + for (size_t i = 0; i < eps_rules_count; i++) { + LAGraph_rule_WCNF eps_rule = rules[eps_rules[i]]; + + GxB_eWiseUnion ( + T[eps_rule.nonterm],GrB_NULL,GxB_PAIR_BOOL,GxB_PAIR_BOOL, + T[eps_rule.nonterm],true_scalar,identity_matrix,true_scalar,GrB_NULL + ); + + t_empty_flags[eps_rule.nonterm] = false; + + // #ifdef DEBUG_CFL_REACHABILITY + // GxB_Matrix_iso(&iso_flag, T[eps_rule.nonterm]); + // printf("[EPS] eWiseUnion: NONTERM: %d (ISO: %d)\n", + // eps_rule.nonterm, iso_flag); + // #endif + } + + // Rule [Variable -> Variable1 Variable2] + bool changed = true; + while (changed) { + changed = false; + for (size_t i = 0; i < bin_rules_count; i++) { + LAGraph_rule_WCNF bin_rule = rules[bin_rules[i]]; + + GRB_TRY(GrB_mxm(M, GrB_NULL, GrB_NULL, GxB_ANY_PAIR_BOOL, + TSrc[bin_rule.nonterm], T[bin_rule.prod_A], GrB_NULL)); + + + GRB_TRY(GrB_mxm(B, GrB_NULL, GrB_NULL, GxB_ANY_PAIR_BOOL, M, T[bin_rule.prod_B], GrB_NULL)); + + GRB_TRY(GrB_eWiseAdd(T[bin_rule.nonterm], GrB_NULL, GrB_NULL, GxB_ANY_BOOL, T[bin_rule.nonterm], B, GrB_NULL)); + + + GRB_TRY(GrB_eWiseAdd(TSrc[bin_rule.prod_A], GrB_NULL, GrB_NULL, GxB_ANY_BOOL, + TSrc[bin_rule.prod_A], TSrc[bin_rule.nonterm], GrB_NULL)); + + // Update source vertices matrix to find appropriate paths only + // M[i, j] == 1 => A[j, j] == 1 + GRB_TRY(GrB_vxm(a, GrB_NULL, GrB_NULL, GxB_ANY_PAIR_BOOL, ones_vec, M, GrB_NULL)); + GrB_Matrix_free(&A); + GRB_TRY(GrB_Matrix_diag(&A, a, 0); + GRB_TRY(GrB_eWiseAdd(TSrc[bin_rule.prod_B], GrB_NULL, GrB_NULL, GxB_ANY_BOOL, + TSrc[bin_rule.prod_B], A, GrB_NULL))); + + // Check if any of the matrices changed. If not, job is done. + GrB_Index nnz_T, nnz_TSrc_B, nnz_TSrc_C; + + GRB_TRY(GrB_Matrix_nvals(&nnz_T, T[bin_rule.nonterm])); + GRB_TRY(GrB_Matrix_nvals(&nnz_TSrc_B, TSrc[bin_rule.prod_A])); + GRB_TRY(GrB_Matrix_nvals(&nnz_TSrc_C, TSrc[bin_rule.prod_B])); + + if (nnz_T != 0) t_empty_flags[bin_rule.nonterm] = false; + if (nnz_TSrc_B != 0) t_src_empty_flags[bin_rule.prod_A] = false; + if (nnz_TSrc_C != 0) t_src_empty_flags[bin_rule.prod_B] = false; + + changed = changed || (nnzs_T[bin_rule.nonterm] != nnz_T); + changed = changed || (nnzs_TSrc_B[bin_rule.prod_A] != nnz_TSrc_B); + changed = changed || (nnzs_TSrc_C[bin_rule.prod_B] != nnz_TSrc_C); + + nnzs_T[bin_rule.nonterm] = nnz_T; + nnzs_TSrc_B[bin_rule.prod_A] = nnz_TSrc_B; + nnzs_TSrc_C[bin_rule.prod_B] = nnz_TSrc_C; + } + } + + // #ifdef DEBUG_CFL_REACHABILITY + // printf("Before MSrc = MSrc * T^S\n"); + // printf("MSrc:\n"); + // PRINT_MATRIX(MSrc) + // printf("T^S:\n"); + // PRINT_MATRIX(T[0]); + // #endif + + GRB_TRY(GrB_mxm(MSrc, GrB_NULL, GrB_NULL, GxB_ANY_PAIR_BOOL, + MSrc, T[0], GrB_NULL)); + + // #ifdef DEBUG_CFL_REACHABILITY + // printf("After MSrc = MSrc * T^S\n"); + // printf("MSrc, output:\n"); + // PRINT_MATRIX(MSrc) + // #endif + + GRB_TRY(GrB_Matrix_dup(output, MSrc)); + + LG_FREE_ALL; + + return GrB_SUCCESS; +} diff --git a/experimental/algorithm/LAGraph_CFL_reachability_multsrc_fast.c b/experimental/algorithm/LAGraph_CFL_reachability_multsrc_fast.c new file mode 100644 index 0000000000..8fb957543f --- /dev/null +++ b/experimental/algorithm/LAGraph_CFL_reachability_multsrc_fast.c @@ -0,0 +1,564 @@ +//------------------------------------------------------------------------------ +// LAGraph_CFL_reachability_multsrc.c: Optimized Multiple-Source Context-Free +// Language Reachability Matrix-Based Algorithm +//------------------------------------------------------------------------------ +// +// LAGraph, (c) 2019-2025 by The LAGraph Contributors, All Rights Reserved. +// SPDX-License-Identifier: BSD-2-Clause + +// Contributed by Ilhom Kombaev, Vlasenco Daniel, Semyon Grigoriev, St. Petersburg State University. + +//------------------------------------------------------------------------------ + +// Code is based on the "multiple-source CFPQ", described in the following paper: +// * Arseny Terekhov et al., "Multiple-Source Context-Free Path Querying in Terms of Linear Algebra" +// * URL: https://openproceedings.org/2021/conf/edbt/p48.pdf + +// #define DEBUG_CFL_REACHABILITY + +#define LG_FREE_WORK \ + { \ + LAGraph_Free((void **) &nnzs_T, msg); \ + LAGraph_Free((void **) &nnzs_TSrc_B, msg); \ + LAGraph_Free((void **) &nnzs_TSrc_C, msg); \ + LAGraph_Free((void **) &ones_vec, msg); \ + LAGraph_Free((void **) &T, msg); \ + LAGraph_Free((void **) &TSrc, msg); \ + LAGraph_Free((void **) &MSrc, msg); \ + LAGraph_Free((void **) &identity_matrix, msg); \ + LAGraph_Free ((void **) &A, msg); \ + LAGraph_Free ((void **) &a, msg); \ + GrB_free(&true_scalar); \ + } + +#define LG_FREE_ALL \ + { \ + for (size_t i = 0; i < nonterms_count; i++) { \ + GrB_free(&T[i].base); \ + GrB_free(&dT[i].base); \ + GrB_free(&TSrc[i].base); \ + } \ + \ + LG_FREE_WORK; \ + } + +#include "LG_internal.h" +#include + +#define ERROR_RULE(msg) \ + { \ + LG_ASSERT_MSGF(false, GrB_INVALID_VALUE, "Rule with index %ld is invalid. " msg, \ + i); \ + } + +#define ADD_TO_MSG(...) \ + { \ + if (msg_len == 0) { \ + msg_len += \ + snprintf(msg, LAGRAPH_MSG_LEN, \ + "LAGraph failure (file %s, line %d): ", __FILE__, __LINE__); \ + } \ + if (msg_len < LAGRAPH_MSG_LEN) { \ + msg_len += snprintf(msg + msg_len, LAGRAPH_MSG_LEN - msg_len, __VA_ARGS__); \ + } \ + } + +#define ADD_INDEX_TO_ERROR_RULE(rule, i) \ + { \ + rule.len_indices_str += snprintf(rule.indices_str + rule.len_indices_str, \ + LAGRAPH_MSG_LEN - rule.len_indices_str, \ + rule.count == 0 ? "%ld" : ", %ld", i); \ + rule.count++; \ + } + +#define PRINT_MATRIX(_m) { \ + for (size_t _i = 0; _i < n; _i++) { \ + for (size_t _j = 0; _j < n; _j++) { \ + if (GxB_Matrix_isStoredElement(_m, _i, _j) == GrB_SUCCESS) { \ + printf("1 "); \ + } \ + else { \ + printf("0 "); \ + } \ + } \ + printf("\n"); \ + } \ +} + +#define PRINT_VECTOR(_v) { \ + for (size_t _i = 0; _i < n; _i++) { \ + if (GxB_Vector_isStoredElement(_v, _i) == GrB_SUCCESS) { \ + printf("1 "); \ + } \ + else { \ + printf("0 "); \ + } \ + } \ + printf("\n"); \ +} + +#define PRINT_RULE(_r) { \ + printf("%c -> ", _r.nonterm + 'A' - 1); \ + if (_r.prod_A == -1 && _r.prod_B == -1) { \ + printf("eps\n"); \ + } else if (_r.prod_A != -1 && _r.prod_B == -1) { \ + printf("%c\n", _r.prod_A + 'a'); \ + } else { \ + printf("%c %c\n", _r.prod_A + 'A'- 1 ,_r.prod_B + 'A' - 1); \ + } \ +} + +#define OPT_EMPTY (1 << 0) +#define OPT_FORMAT (1 << 1) +#define OPT_LAZY (1 << 2) +#define OPT_BLOCK (1 << 3) + +// LAGraph_CFL_reachability_multsrc: Optimized Multiple-Source Context-Free +// Language Reachability Matrix-Based Algorithm +// +// This function determines the set of vertex pairs (u, v) in a graph (represented by +// adjacency matrices) such that u is in the given set of source vertices `src`, and +// there is a path from u to v, where the edge labels form a word from the language +// generated by the context-free grammar (represented by `rules`). + +GrB_Info LAGraph_CFL_reachability_multsrc_fast +( + // Output + GrB_Matrix *output, // A handle for the matrix containing results. + // + // output: (i, j) = true if and only if there is a path + // from node i to node j whose edge labels form a word + // derivable from the specified CFG. + // Input + const GrB_Matrix *adj_matrices, // Array of adjacency matrices representing the graph. + // The length of this array is equal to the count of + // terminals (terms_count). + // + // adj_matrices[t]: (i, j) == 1 if and only if there + // is an edge between nodes i and j with the label of + // the terminal corresponding to index 't' (where t is + // in the range [0, terms_count - 1]). + GrB_Index *src, // Array of source vertices + // These parameters will become unsigned + int32_t src_count, // The total number of source vertices + int32_t terms_count, // The total number of terminal symbols in the CFG. + int32_t nonterms_count, // The total number of non-terminal symbols in the CFG. + const LAGraph_rule_WCNF *rules, // The rules of the CFG. + size_t rules_count, // The total number of rules in the CFG. + char *msg, // Message string for error reporting. + int8_t opt_mask // Optimizations mask +) +{ + // Declare workspace and clear the msg string, if not NULL + CFL_Matrix *T; + CFL_Matrix *dT; + CFL_Matrix *TSrc; + CFL_Matrix *Adj; + CFL_Matrix MSrc; + CFL_Matrix M1; + CFL_Matrix M2; + CFL_Matrix Temp1; + CFL_Matrix Temp2; + CFL_Matrix A; + GrB_Vector a; + GrB_Index n; // number of vertices in the graph + CFL_Matrix iden; + GrB_Matrix identity_matrix = NULL; + GrB_Index *nnzs_T = NULL; + GrB_Index *nnzs_TSrc_B = NULL; + GrB_Index *nnzs_TSrc_C = NULL; + LG_CLEAR_MSG; + size_t msg_len = 0; // For error formatting + bool iso_flag = false; + GrB_Scalar true_scalar; + GrB_Vector ones_vec; + + // Will change the interface and omit this check in the future + if (nonterms_count < 0 || terms_count < 0) + return GrB_INVALID_VALUE; + + if (!nonterms_count || !rules_count) + return GrB_INVALID_VALUE; + + if (!output || !rules || !adj_matrices || !src) + return GrB_NULL_POINTER; + if (src_count <= 0) + return GrB_INVALID_VALUE; + + // Find null adjacency matrices + bool found_null = false; + for (int32_t i = 0; i < terms_count; i++) { + if (adj_matrices[i] != NULL) + continue; + + if (!found_null) { + ADD_TO_MSG("Adjacency matrices with these indices are null: "); + ADD_TO_MSG("%d", i); + } else { + ADD_TO_MSG(", %d", i); + } + + found_null = true; + } + if (found_null) { + return GrB_NULL_POINTER; + } + + GRB_TRY(GrB_Matrix_ncols(&n, adj_matrices[0])); + if (n < src_count) return GrB_INVALID_VALUE; + + GrB_Scalar_new(&true_scalar, GrB_BOOL); + GrB_Scalar_setElement_BOOL(true_scalar, true); + + LG_TRY(LAGraph_Calloc((void **) &T, nonterms_count, sizeof(CFL_Matrix), msg)); + LG_TRY(LAGraph_Calloc((void **) &dT, nonterms_count, sizeof(CFL_Matrix), msg)); + LG_TRY(LAGraph_Calloc((void **) &TSrc, nonterms_count, sizeof(CFL_Matrix), msg)); + + LG_TRY(LAGraph_Calloc((void **) &Adj, terms_count, sizeof(CFL_Matrix), msg)); + + GRB_TRY(GrB_Vector_new(&ones_vec, GrB_BOOL, n)); + GRB_TRY(GrB_Vector_assign_BOOL(ones_vec, GrB_NULL, GrB_NULL, true, GrB_ALL, n, NULL)); + GRB_TRY(GrB_Matrix_diag(&identity_matrix, ones_vec, 0)); + + LG_TRY(LAGraph_Calloc((void **) &nnzs_T, nonterms_count, sizeof(GrB_Index), msg)); + LG_TRY(LAGraph_Calloc((void **) &nnzs_TSrc_B, nonterms_count, sizeof(GrB_Index), msg)); + LG_TRY(LAGraph_Calloc((void **) &nnzs_TSrc_C, nonterms_count, sizeof(GrB_Index), msg)); + + for (int32_t i = 0; i < terms_count; i++) { + Adj[i] = ((opt_mask & OPT_LAZY) || (opt_mask & OPT_BLOCK)) + ? CFL_matrix_from_base_lazy(adj_matrices[i]) + : CFL_matrix_from_base(adj_matrices[i]); + } + + // Create nonterms matrices + for (int32_t i = 0; i < nonterms_count; i++) { + GrB_Matrix matrix; + + GRB_TRY(GrB_Matrix_new(&matrix, GrB_BOOL, n, n)); + T[i] = ((opt_mask & OPT_LAZY) || (opt_mask & OPT_BLOCK)) + ? CFL_matrix_from_base_lazy(matrix) + : CFL_matrix_from_base(matrix); + + GRB_TRY(GrB_Matrix_new(&matrix, GrB_BOOL, n, n)); + dT[i] = ((opt_mask & OPT_LAZY) || (opt_mask & OPT_BLOCK)) + ? CFL_matrix_from_base_lazy(matrix) + : CFL_matrix_from_base(matrix); + + TSrc[i] = CFL_matrix_create(n, n); + } + + for (int32_t i = 0; i < src_count; i++) { + GrB_Matrix_setElement(TSrc[0].base, true, src[i], src[i]); + } + TSrc[0] = ((opt_mask & OPT_LAZY) || (opt_mask & OPT_BLOCK)) + ? CFL_matrix_from_base_lazy(TSrc[0].base) + : CFL_matrix_from_base(TSrc[0].base); + + MSrc = CFL_matrix_create(n, n); + CFL_dup(&MSrc, &TSrc[0], opt_mask); + + GRB_TRY(GrB_Vector_new(&a, GrB_BOOL, n)); + + M1 = CFL_matrix_create(n, n); + M2 = CFL_matrix_create(n, n); + Temp1 = CFL_matrix_create(n, n); + Temp2 = CFL_matrix_create(n, n); + A = CFL_matrix_create(n, n); + + // #ifdef DEBUG_CFL_REACHABILITY + // printf("MSrc:\n"); + // PRINT_MATRIX(MSrc); + // #endif + + // Arrays for processing rules + size_t eps_rules[rules_count], eps_rules_count = 0; // [Variable -> eps] + size_t term_rules[rules_count], term_rules_count = 0; // [Variable -> term] + size_t bin_rules[rules_count], bin_rules_count = 0; // [Variable -> AB] + + // Process rules + typedef struct { + size_t count; + size_t len_indices_str; + char indices_str[LAGRAPH_MSG_LEN]; + } rule_error_s; + + rule_error_s term_err = {0}; + rule_error_s nonterm_err = {0}; + rule_error_s invalid_err = {0}; + + for (size_t i = 0; i < rules_count; i++) { + LAGraph_rule_WCNF rule = rules[i]; + + bool is_rule_eps = rule.prod_A == -1 && rule.prod_B == -1; + bool is_rule_term = rule.prod_A != -1 && rule.prod_B == -1; + bool is_rule_bin = rule.prod_A != -1 && rule.prod_B != -1; + + // Check that all rules are well-formed + if (rule.nonterm < 0 || rule.nonterm >= nonterms_count) { + ADD_INDEX_TO_ERROR_RULE(nonterm_err, i); + } + + // [Variable -> eps] + if (is_rule_eps) { + eps_rules[eps_rules_count++] = i; + + continue; + } + + // [Variable -> term] + if (is_rule_term) { + term_rules[term_rules_count++] = i; + + if (rule.prod_A < -1 || rule.prod_A >= terms_count) { + ADD_INDEX_TO_ERROR_RULE(term_err, i); + } + + continue; + } + + // [Variable -> A B] + if (is_rule_bin) { + bin_rules[bin_rules_count++] = i; + + if (rule.prod_A < -1 || rule.prod_A >= nonterms_count || rule.prod_B < -1 || + rule.prod_B >= nonterms_count) { + ADD_INDEX_TO_ERROR_RULE(nonterm_err, i); + } + + continue; + } + + // [Variable -> _ B] + ADD_INDEX_TO_ERROR_RULE(invalid_err, i); + } + + if (term_err.count + nonterm_err.count + invalid_err.count > 0) { + ADD_TO_MSG("Count of invalid rules: %ld.\n", + term_err.count + nonterm_err.count + invalid_err.count); + + if (nonterm_err.count > 0) { + ADD_TO_MSG("Non-terminals must be in range [0, nonterms_count). "); + ADD_TO_MSG("indices of invalid rules: %s\n", nonterm_err.indices_str) + } + if (term_err.count > 0) { + ADD_TO_MSG("Terminals must be in range [-1, nonterms_count). "); + ADD_TO_MSG("indices of invalid rules: %s\n", term_err.indices_str) + } + if (invalid_err.count > 0) { + ADD_TO_MSG("[Variable -> _ B] type of rule is not acceptable. "); + ADD_TO_MSG("indices of invalid rules: %.120s\n", invalid_err.indices_str) + } + + LG_FREE_ALL; + return GrB_INVALID_VALUE; + } + + // Rule [Variable -> term] + for (size_t i = 0; i < term_rules_count; i++) { + LAGraph_rule_WCNF term_rule = rules[term_rules[i]]; + + GRB_TRY(CFL_wise(&T[term_rule.nonterm], &T[term_rule.nonterm], &Adj[term_rule.prod_A], true, opt_mask)); + GRB_TRY(CFL_wise(&dT[term_rule.nonterm], &dT[term_rule.nonterm], &Adj[term_rule.prod_A], true, opt_mask)); + + // #ifdef DEBUG_CFL_REACHABILITY + // GxB_Matrix_iso(&iso_flag, T[term_rule.nonterm]); + // printf("[TERM] eWiseUnion: NONTERM: %d (ISO: %d)\n", term_rule.nonterm, iso_flag); + // #endif + } + + iden = CFL_matrix_from_base(identity_matrix); + + // Rule [Variable -> eps] + for (size_t i = 0; i < eps_rules_count; i++) { + LAGraph_rule_WCNF eps_rule = rules[eps_rules[i]]; + + GRB_TRY(CFL_wise(&T[eps_rule.nonterm], &T[eps_rule.nonterm], &iden, true, opt_mask)); + GRB_TRY(CFL_wise(&dT[eps_rule.nonterm], &dT[eps_rule.nonterm], &iden, true, opt_mask)); + + // #ifdef DEBUG_CFL_REACHABILITY + // GxB_Matrix_iso(&iso_flag, T[eps_rule.nonterm]); + // printf("[EPS] eWiseUnion: NONTERM: %d (ISO: %d)\n", + // eps_rule.nonterm, iso_flag); + // #endif + } + + // Rule [Variable -> Variable1 Variable2] + bool changed = true; + while (changed) { + changed = false; + for (size_t i = 0; i < bin_rules_count; i++) { + LAGraph_rule_WCNF bin_rule = rules[bin_rules[i]]; + // #ifdef DEBUG_CFL_REACHABILITY + // printf("Rule: "); + // PRINT_RULE(bin_rule) + // #endif + + // #ifdef DEBUG_CFL_REACHABILITY + // printf("Before M = TSrc^A * T^B:\n"); + // printf("TSrc^A:\n"); + // PRINT_MATRIX(TSrc[bin_rule.nonterm]); + // printf("T^B:\n"); + // PRINT_MATRIX(T[bin_rule.prod_A]); + // #endif + + + GRB_TRY(CFL_mxm(&M1, &TSrc[bin_rule.nonterm], &T[bin_rule.prod_A], false, false, opt_mask)); + + + GRB_TRY(CFL_wise(&T[bin_rule.nonterm], &T[bin_rule.nonterm], &dT[bin_rule.nonterm], false, opt_mask)); + + + // ? <- this means that without this step, unit tests pass, + // although algorithm requires it. + + // printf("TSrc:\n"); + // PRINT_MATRIX(TSrc[bin_rule.nonterm].base); + // printf("dT:\n"); + // PRINT_MATRIX(dT[bin_rule.prod_A].base); + + GRB_TRY(CFL_mxm(&M2, &TSrc[bin_rule.nonterm], &dT[bin_rule.prod_A], false, false, opt_mask)); + + // printf("row M2 = TSrc * dT:\n"); + // PRINT_MATRIX(M2.base_row); + // printf("col M2 = TSrc * dT:\n"); + // PRINT_MATRIX(M2.base_col); + + // #ifdef DEBUG_CFL_REACHABILITY + // printf("After M = TSrc^A * T^B:\n"); + // printf("M:\n"); + // PRINT_MATRIX(M); + // printf("Before T^A = T^A + M * T^C:\n"); + // printf("T^A:\n"); + // PRINT_MATRIX(T[bin_rule.nonterm]); + // printf("T^C:\n"); + // PRINT_MATRIX(T[bin_rule.prod_B]); + // #endif + + GRB_TRY(CFL_mxm(&Temp1, &M1, &dT[bin_rule.prod_B], false, false, opt_mask)); + + // ? + GRB_TRY(CFL_mxm(&Temp2, &M2, &T[bin_rule.prod_B], false, false, opt_mask)); + + GRB_TRY(CFL_dup(&dT[bin_rule.nonterm], &Temp1, opt_mask)); + + GRB_TRY(CFL_wise(&dT[bin_rule.nonterm], &dT[bin_rule.nonterm], &Temp2, false, opt_mask)); + + // ? + GRB_TRY(CFL_rsub(&dT[bin_rule.nonterm], &T[bin_rule.nonterm], opt_mask)); + + // #ifdef DEBUG_CFL_REACHABILITY + // printf("After T^A = T^A + M * T^C:\n"); + // printf("T^A:\n"); + // PRINT_MATRIX(T[bin_rule.nonterm]); + // #endif + + // #ifdef DEBUG_CFL_REACHABILITY + // printf("Before TSrc^B = TSrc^B + TSrc^A:\n"); + // printf("TSrc^A:\n"); + // PRINT_MATRIX(TSrc[bin_rule.nonterm]); + // printf("TSrc^B:\n"); + // PRINT_MATRIX(TSrc[bin_rule.prod_A]); + // #endif + + // #ifdef DEBUG_CFL_REACHABILITY + // printf("After TSrc^B = TSrc^B + TSrc^A:\n"); + // printf("TSrc^B:\n"); + // PRINT_MATRIX(TSrc[bin_rule.prod_A]); + // #endif + + // #ifdef DEBUG_CFL_REACHABILITY + // printf("Before A = dest(M)\n"); + // printf("M:\n"); + // PRINT_MATRIX(M); + // #endif + + // Update source vertices matrix to find appropriate paths only + // M1[i, j] == 1 => A[j, j] == 1 + GRB_TRY(GrB_vxm(a, GrB_NULL, GrB_NULL, GxB_ANY_PAIR_BOOL, ones_vec, M1.base, GrB_NULL)); + GRB_TRY(GrB_Matrix_free(&A.base)); + GRB_TRY(GrB_Matrix_diag(&A.base, a, 0)); + A = CFL_matrix_from_base(A.base); + + // #ifdef DEBUG_CFL_REACHABILITY + // printf("After A = dest(M)\n"); + // printf("a:\n"); + // PRINT_VECTOR(a); + // printf("A:\n"); + // PRINT_MATRIX(A); + // #endif + + // #ifdef DEBUG_CFL_REACHABILITY + // printf("Before TSrc^C = TSrc^c + A\n"); + // printf("TSrc^C:\n"); + // PRINT_MATRIX(TSrc[bin_rule.prod_B]) + // printf("A:\n"); + // PRINT_MATRIX(A) + // #endif + + GRB_TRY(CFL_wise(&TSrc[bin_rule.prod_A], &TSrc[bin_rule.prod_A], &TSrc[bin_rule.nonterm], false, opt_mask)); + + GRB_TRY(CFL_wise(&TSrc[bin_rule.prod_B], &TSrc[bin_rule.prod_B], &A, false, opt_mask)); + + // #ifdef DEBUG_CFL_REACHABILITY + // printf("After TSrc^C = TSrc^c + A\n"); + // printf("TSrc^C:\n"); + // PRINT_MATRIX(TSrc[bin_rule.prod_B]) + // #endif + + // #ifdef DEBUG_CFL_REACHABILITY + // GxB_Matrix_iso(&iso_flag, T[bin_rule.nonterm]); + // printf("[TERM1 TERM2] MULTIPLY, S: %d, A: %d, B: %d, " + // "I: %ld (ISO: %d)\n", + // bin_rule.nonterm, bin_rule.prod_A, bin_rule.prod_B, i, iso_flag); + // #endif + + // Check if any of the matrices changed. If not, job is done. + GrB_Index nnz_T, nnz_TSrc_B, nnz_TSrc_C; + + CFL_matrix_update(&T[bin_rule.nonterm]); + CFL_matrix_update(&TSrc[bin_rule.prod_A]); + CFL_matrix_update(&TSrc[bin_rule.prod_B]); + + nnz_T = T[bin_rule.nonterm].nvals; + nnz_TSrc_B = TSrc[bin_rule.prod_A].nvals; + nnz_TSrc_C = TSrc[bin_rule.prod_B].nvals; + + changed = changed || (nnzs_T[bin_rule.nonterm] != nnz_T); + changed = changed || (nnzs_TSrc_B[bin_rule.prod_A] != nnz_TSrc_B); + changed = changed || (nnzs_TSrc_C[bin_rule.prod_B] != nnz_TSrc_C); + + nnzs_T[bin_rule.nonterm] = nnz_T; + nnzs_TSrc_B[bin_rule.prod_A] = nnz_TSrc_B; + nnzs_TSrc_C[bin_rule.prod_B] = nnz_TSrc_C; + } + } + + // #ifdef DEBUG_CFL_REACHABILITY + // printf("Before MSrc = MSrc * T^S\n"); + // printf("MSrc:\n"); + // PRINT_MATRIX(MSrc) + // printf("T^S:\n"); + // PRINT_MATRIX(T[0]); + // #endif + + GRB_TRY(CFL_mxm(&MSrc, &MSrc, &T[0], false, false, opt_mask)); + // GRB_TRY(GrB_mxm(MSrc.base, GrB_NULL, GrB_NULL, GxB_ANY_PAIR_BOOL, + // MSrc.base, T[0].base, GrB_NULL)); + + // #ifdef DEBUG_CFL_REACHABILITY + // printf("After MSrc = MSrc * T^S\n"); + // printf("MSrc, output:\n"); + // PRINT_MATRIX(MSrc) + // #endif + + if (MSrc.base_matrices_count == 0) { + GRB_TRY(GrB_Matrix_dup(output, MSrc.base)); + } else { + CFL_Matrix res = CFL_matrix_to_base(&MSrc, opt_mask); + GRB_TRY(GrB_Matrix_dup(output, res.base)); + CFL_matrix_free(&res); + } + + LG_FREE_ALL; + return GrB_SUCCESS; +} diff --git a/experimental/test/test_CFL_optimized_matrix_opt.c b/experimental/test/test_CFL_optimized_matrix_opt.c new file mode 100644 index 0000000000..bad66e0640 --- /dev/null +++ b/experimental/test/test_CFL_optimized_matrix_opt.c @@ -0,0 +1,1279 @@ +//------------------------------------------------------------------------------ +// LAGraph/experimental/test/LAGraph_CFL_reachability.c: test cases for +// operations for Optimized Context-Free Language Reachability Matrix-Based +// Algorithm +//------------------------------------------------------------------------------ +// +// LAGraph, (c) 2019-2024 by The LAGraph Contributors, All Rights Reserved. +// SPDX-License-Identifier: BSD-2-Clause + +// Contributed by Ilhom Kombaev, Semyon Grigoriev, St. Petersburg State +// University. + +//------------------------------------------------------------------------------ + +#include +#include +#include +#include +#include +#include + +#define OPT_EMPTY (1 << 0) +#define OPT_FORMAT (1 << 1) +#define OPT_LAZY (1 << 2) +#define OPT_BLOCK (1 << 3) + +#define OPT_N_FLAGS 16 + +char msg[LAGRAPH_MSG_LEN]; + +typedef CFL_Matrix Matrix; +typedef enum CFL_Matrix_block Matrix_block; + +extern GrB_Info matrix_to_format(Matrix *matrix, int32_t format, bool is_bool); +extern GrB_Info matrix_clear_format(Matrix *A, int8_t optimizations); +extern GrB_Info matrix_dup_format(Matrix *output, Matrix *input, int8_t optimizations); +extern GrB_Info block_matrix_hyper_rotate_i(Matrix *matrix, enum CFL_Matrix_block format); +extern void block_matrix_to_diag(Matrix *diag, Matrix *input); +extern void block_matrix_reduce(Matrix *matrix, Matrix *input, int8_t optimizations); + +// extern GrB_Info matrix_clear_empty(Matrix *A, int8_t optimizations); +// extern GrB_Info matrix_mxm_empty(Matrix *output, Matrix *first, Matrix *second, +// bool accum, bool swap, int8_t optimizations); +// extern GrB_Info matrix_wise_empty(Matrix *output, Matrix *first, Matrix *second, +// bool accum, int8_t optimizations); + +//------------------------------------------------------------------------------ +// helpers +//------------------------------------------------------------------------------ + +static void setup(void) { OK(LAGraph_Init(msg)); } + +static void teardown(void) { OK(LAGraph_Finalize(msg)); } + +// 0 1 .. 0 +// 1 0 .. 0 +// .. .. .. .. +// 0 0 .. 0 +static Matrix make_simple_matrix(int n) { + GrB_Matrix A; + OK(GrB_Matrix_new(&A, GrB_BOOL, n, n)); + OK(GrB_Matrix_setElement_BOOL(A, true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(A, true, 1, 0)); + Matrix M = CFL_matrix_from_base(A); + return M; +} + +// 1 0 .. 0 +// 0 1 .. 0 +// .. .. .. .. +// 0 0 .. 0 +static Matrix make_simple_matrix_inverted(int n) { + GrB_Matrix A; + OK(GrB_Matrix_new(&A, GrB_BOOL, n, n)); + OK(GrB_Matrix_setElement_BOOL(A, true, 0, 0)); + OK(GrB_Matrix_setElement_BOOL(A, true, 1, 1)); + Matrix M = CFL_matrix_from_base(A); + return M; +} + +// 1 1 .. 1 +// 1 1 .. 1 +// .. .. .. .. +// 1 1 .. 1 +static Matrix make_ones_matrix(int n) { + GrB_Matrix A; + OK(GrB_Matrix_new(&A, GrB_BOOL, n, n)); + + for (size_t i = 0; i < n; i++) { + for (size_t j = 0; j < n; j++) { + OK(GrB_Matrix_setElement_BOOL(A, true, i, j)); + } + } + + Matrix M = CFL_matrix_from_base(A); + return M; +} + +static void free_matrix(Matrix *M) { CFL_matrix_free(M); } + +static bool compare_matrices(Matrix first, Matrix second) { + Matrix A, B, C; + A = CFL_matrix_to_base(&first, OPT_LAZY); + B = CFL_matrix_to_base(&second, OPT_LAZY); + C = CFL_matrix_create(A.nrows, A.ncols); + + if (A.nvals != B.nvals) { + CFL_matrix_free(&A); + CFL_matrix_free(&B); + CFL_matrix_free(&C); + + return false; + } + + GrB_eWiseMult(C.base, GrB_NULL, false, GrB_EQ_BOOL, A.base, B.base, GrB_NULL); + CFL_matrix_update(&C); + + bool result = C.nvals == A.nvals; + + CFL_matrix_free(&A); + CFL_matrix_free(&B); + CFL_matrix_free(&C); + + return result; +} + +static void test_compare_matrices_function(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + { + Matrix A = make_simple_matrix(5); + Matrix B = CFL_matrix_create(5, 5); + + TEST_CHECK(!compare_matrices(A, B)); + + CFL_matrix_free(&A); + CFL_matrix_free(&B); + } + + { + Matrix A = make_simple_matrix(5); + Matrix B = make_simple_matrix_inverted(5); + + TEST_CHECK(!compare_matrices(A, B)); + TEST_CHECK(compare_matrices(A, A)); + TEST_CHECK(compare_matrices(B, B)); + + CFL_matrix_free(&A); + CFL_matrix_free(&B); + } + + teardown(); +#endif +} + +//------------------------------------------------------------------------------ +// Сreation and Free +//------------------------------------------------------------------------------ + +static void test_CFL_create_free(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + CFL_Matrix A = CFL_matrix_create(4, 4); + TEST_CHECK(A.nrows == 4); + TEST_CHECK(A.ncols == 4); + TEST_CHECK(A.base != NULL); + + CFL_matrix_free(&A); + TEST_CHECK(A.base == NULL); + + teardown(); +#endif +} + +//------------------------------------------------------------------------------ +// Multiplication +//------------------------------------------------------------------------------ + +static void test_CFL_mxm(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (int mask = 0; mask < OPT_N_FLAGS; mask++) { + Matrix A = make_simple_matrix(2); + Matrix B = make_simple_matrix(2); + matrix_to_format(&A, GrB_ROWMAJOR, true); + Matrix C = CFL_matrix_create(2, 2); + + GrB_Info info = CFL_mxm(&C, &A, &B, false, false, mask); + TEST_CHECK(info == GrB_SUCCESS); + TEST_MSG("matrix_mxm_opt failed for mask=%d", mask); + + // Validate result (A*B should have two true element) + GrB_Index nvals = 0; + OK(GrB_Matrix_nvals(&nvals, C.base)); + TEST_CHECK(nvals == 2); + TEST_CHECK(C.nvals == 2); + TEST_MSG("Unexpected empty result, mask=%d", mask); + + free_matrix(&A); + free_matrix(&B); + free_matrix(&C); + } + + teardown(); +#endif +} + +//------------------------------------------------------------------------------ +// Wise +//------------------------------------------------------------------------------ + +static void test_CFL_wise(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (int mask = 0; mask < OPT_N_FLAGS; mask++) { + Matrix A = make_simple_matrix(2); + Matrix B = make_simple_matrix(2); + matrix_to_format(&A, GrB_ROWMAJOR, true); + + GrB_Info info = CFL_wise(&A, &A, &B, false, mask); + TEST_CHECK(info == GrB_SUCCESS); + TEST_MSG("matrix_wise failed for mask=%d", mask); + + // Validate result (A+B should have two true element) + GrB_Index nvals = 0; + OK(GrB_Matrix_nvals(&nvals, A.base)); + TEST_CHECK(nvals == 2); + TEST_CHECK(A.nvals == 2); + TEST_MSG("Unexpected empty result, mask=%d", mask); + + free_matrix(&A); + free_matrix(&B); + } + + teardown(); +#endif +} + +//------------------------------------------------------------------------------ +// Rsub +//------------------------------------------------------------------------------ + +static void test_CFL_rsub(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (int mask = 0; mask < OPT_N_FLAGS; mask++) { + Matrix A = make_simple_matrix(2); + Matrix B = make_simple_matrix(2); + matrix_to_format(&A, GrB_ROWMAJOR, true); + + GrB_Info info = CFL_rsub(&A, &B, mask); + TEST_CHECK(info == GrB_SUCCESS); + TEST_MSG("matrix_rsub failed for mask=%d", mask); + + // Validate result (A-B should have zero true element) + GrB_Index nvals = 0; + OK(GrB_Matrix_nvals(&nvals, A.base)); + TEST_CHECK(nvals == 0); + TEST_CHECK(A.nvals == 0); + TEST_MSG("Unexpected non-empty result, mask=%d", mask); + + free_matrix(&A); + free_matrix(&B); + } + + teardown(); +#endif +} + +//------------------------------------------------------------------------------ +// Clear +//------------------------------------------------------------------------------ + +static void test_CFL_clear(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (int mask = 0; mask < OPT_N_FLAGS; mask++) { + Matrix A = make_simple_matrix(3); + + // Clear nonempty matrix + CFL_clear(&A, mask); + GrB_Index nvals; + TEST_CHECK(A.nvals == 0); + + free_matrix(&A); + } + + teardown(); +#endif +} + +//------------------------------------------------------------------------------ +// Duplicate +//------------------------------------------------------------------------------ + +static void test_CFL_dup(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (int mask = 0; mask < OPT_N_FLAGS; mask++) { + CFL_Matrix A = make_simple_matrix(3); + CFL_Matrix B = CFL_matrix_create(3, 3); + + OK(CFL_dup(&A, &B, mask)); + + GrB_Index nvals_A, nvals_B; + OK(GrB_Matrix_nvals(&nvals_A, A.base)); + OK(GrB_Matrix_nvals(&nvals_B, B.base)); + + TEST_CHECK(nvals_A == nvals_B); + TEST_CHECK(A.nvals == B.nvals); + free_matrix(&A); + free_matrix(&B); + } + + teardown(); +#endif +} + +static void test_CFL_dup_same(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + CFL_Matrix A = make_simple_matrix(3); + OK(CFL_dup(&A, &A, 0)); + + free_matrix(&A); + + teardown(); +#endif +} + +//------------------------------------------------------------------------------ +// Format optimization +//------------------------------------------------------------------------------ + +static void test_CFL_format_create_rowmajor_matrix(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix A = make_simple_matrix(5); + TEST_CHECK(A.is_both == false); + TEST_CHECK(A.base_col == NULL); + TEST_CHECK(A.base_row != NULL); + TEST_CHECK(A.base == A.base_row); + TEST_CHECK(A.format == GrB_ROWMAJOR); + + CFL_matrix_free(&A); + + teardown(); +#endif +} + +static void test_CFL_format_to_format(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix A = make_simple_matrix(5); + + // same + OK(matrix_to_format(&A, GrB_ROWMAJOR, false)); + TEST_CHECK(A.is_both == false); + TEST_CHECK(A.base_col == NULL); + TEST_CHECK(A.base_row != NULL); + TEST_CHECK(A.base == A.base_row); + TEST_CHECK(A.format == GrB_ROWMAJOR); + + // switch to another + OK(matrix_to_format(&A, GrB_COLMAJOR, false)); + TEST_CHECK(A.is_both == false); + TEST_CHECK(A.base_col != NULL); + TEST_CHECK(A.base_row == NULL); + TEST_CHECK(A.base == A.base_col); + TEST_CHECK(A.format == GrB_COLMAJOR); + + // create both matrices, same format + OK(matrix_to_format(&A, GrB_COLMAJOR, true)); + TEST_CHECK(A.is_both == false); + TEST_CHECK(A.base_col != NULL); + TEST_CHECK(A.base_row == NULL); + TEST_CHECK(A.base == A.base_col); + TEST_CHECK(A.format == GrB_COLMAJOR); + + // switch to another, both + OK(matrix_to_format(&A, GrB_ROWMAJOR, true)); + TEST_CHECK(A.is_both == true); + TEST_CHECK(A.base_col != NULL); + TEST_CHECK(A.base_row != NULL); + TEST_CHECK(A.base == A.base_row); + TEST_CHECK(A.format == GrB_ROWMAJOR); + + // switch to another, when matrix already in both state + OK(matrix_to_format(&A, GrB_COLMAJOR, true)); + TEST_CHECK(A.is_both == true); + TEST_CHECK(A.base_col != NULL); + TEST_CHECK(A.base_row != NULL); + TEST_CHECK(A.base == A.base_col); + TEST_CHECK(A.format == GrB_COLMAJOR); + + CFL_matrix_free(&A); + + teardown(); +#endif +} + +static void test_CFL_format_clear_format(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix A = make_simple_matrix(5); + + // create matrix in both state + OK(matrix_to_format(&A, GrB_COLMAJOR, true)); + OK(matrix_clear_format(&A, OPT_FORMAT)); + + TEST_CHECK(A.nvals == 0); + + // NULL Matrix + GrB_Matrix old_base = A.base; + GrB_Matrix old_base_row = A.base_row; + GrB_Matrix old_base_col = A.base_col; + A.base = NULL; + A.base_row = NULL; + A.base_col = NULL; + GrB_Info result = matrix_clear_format(&A, OPT_FORMAT); + OK(!result); + + A.base = old_base; + A.base_row = old_base_row; + A.base_col = old_base_col; + CFL_matrix_free(&A); + + teardown(); +#endif +} + +static void test_CFL_format_dup_format(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix A = make_simple_matrix(5); + Matrix B = CFL_matrix_create(5, 5); + + // create matrix in both state + OK(matrix_to_format(&A, GrB_COLMAJOR, true)); + OK(matrix_dup_format(&A, &B, OPT_FORMAT)); + + TEST_CHECK(A.nvals == B.nvals); + TEST_CHECK(A.base != B.base); + + TEST_CHECK(compare_matrices(A, B)); + + // NULL Matrix + GrB_Matrix old_base = A.base; + GrB_Matrix old_base_row = A.base_row; + GrB_Matrix old_base_col = A.base_col; + A.base = NULL; + A.base_row = NULL; + A.base_col = NULL; + GrB_Info result = matrix_dup_format(&A, &B, OPT_FORMAT); + OK(!result); + + A.base = old_base; + A.base_row = old_base_row; + A.base_col = old_base_col; + CFL_matrix_free(&A); + CFL_matrix_free(&B); + + // Without optimization flag + A = make_simple_matrix(5); + B = CFL_matrix_create(5, 5); + OK(matrix_to_format(&A, GrB_COLMAJOR, true)); + OK(matrix_dup_format(&A, &B, 0)); + CFL_matrix_free(&A); + CFL_matrix_free(&B); + + teardown(); +#endif +} + +static void test_CFL_format_mxm_second_greather_then_k(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix A = CFL_matrix_create(5, 5); + Matrix B = make_ones_matrix(5); + Matrix C = CFL_matrix_create(5, 5); + + OK(matrix_to_format(&A, GrB_COLMAJOR, false)); + OK(CFL_mxm(&C, &A, &B, false, false, OPT_FORMAT)); + TEST_CHECK(C.nvals == 0); + + free_matrix(&A); + free_matrix(&B); + free_matrix(&C); + + teardown(); +#endif +} + +static void test_CFL_format_wise_when_both(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix A = CFL_matrix_create(5, 5); + Matrix B = make_ones_matrix(5); + + OK(matrix_to_format(&A, GrB_COLMAJOR, true)); + OK(CFL_wise(&A, &A, &B, false, OPT_FORMAT)); + TEST_CHECK(A.nvals == 25); + + // NULL Matrix + GrB_Matrix old_base = A.base; + GrB_Matrix old_base_row = A.base_row; + GrB_Matrix old_base_col = A.base_col; + A.base = NULL; + A.base_row = NULL; + A.base_col = NULL; + GrB_Info result = CFL_wise(&A, &A, &B, false, OPT_FORMAT); + A.base = old_base; + A.base_row = old_base_row; + A.base_col = old_base_col; + OK(!result); + + CFL_matrix_free(&A); + CFL_matrix_free(&B); + + teardown(); +#endif +} + +//------------------------------------------------------------------------------ +// Empty optimization +//------------------------------------------------------------------------------ + +static void test_CFL_empty_clear_empty_matrix(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix A = CFL_matrix_create(5, 5); + + OK(CFL_clear(&A, OPT_EMPTY)); + TEST_CHECK(A.nvals == 0); + + CFL_matrix_free(&A); + + teardown(); +#endif +} + +static void test_CFL_empty_dup(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix A = CFL_matrix_create(5, 5); + Matrix B = make_ones_matrix(5); + + OK(CFL_dup(&A, &B, OPT_EMPTY)); + TEST_CHECK(A.nvals == 25); + + CFL_matrix_free(&A); + CFL_matrix_free(&B); + + teardown(); +#endif +} + +static void test_CFL_empty_mxm_one_is_empty(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix A, B, C; + + // accum + A = CFL_matrix_create(5, 5); + B = make_ones_matrix(5); + C = make_ones_matrix(5); + + OK(CFL_mxm(&C, &A, &B, true, false, OPT_EMPTY)); + TEST_CHECK(C.nvals == 25); + + CFL_matrix_free(&A); + CFL_matrix_free(&B); + CFL_matrix_free(&C); + + // without accum + A = CFL_matrix_create(5, 5); + B = make_ones_matrix(5); + C = make_ones_matrix(5); + + OK(CFL_mxm(&C, &A, &B, false, false, OPT_EMPTY)); + TEST_CHECK(C.nvals == 0); + + CFL_matrix_free(&A); + CFL_matrix_free(&B); + CFL_matrix_free(&C); + + // output empty + A = CFL_matrix_create(5, 5); + B = make_ones_matrix(5); + C = CFL_matrix_create(5, 5); + + OK(CFL_mxm(&C, &A, &B, false, false, OPT_EMPTY)); + TEST_CHECK(C.nvals == 0); + + CFL_matrix_free(&A); + CFL_matrix_free(&B); + CFL_matrix_free(&C); + + teardown(); +#endif +} + +// wise(C, A, B, accum) +// matrix may be empty and full +// we have 16 states +static void test_CFL_empty_wise(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix A, B, C; + + for (size_t is_first_empty = 0; is_first_empty < 2; is_first_empty++) { + for (size_t is_second_empty = 0; is_second_empty < 2; is_second_empty++) { + for (size_t is_output_empty = 0; is_output_empty < 2; is_output_empty++) { + for (size_t is_accum = 0; is_accum < 2; is_accum++) { + A = is_first_empty ? CFL_matrix_create(5, 5) : make_ones_matrix(5); + B = is_second_empty ? CFL_matrix_create(5, 5) : make_ones_matrix(5); + C = is_output_empty ? CFL_matrix_create(5, 5) : make_ones_matrix(5); + + CFL_wise(&C, &A, &B, is_accum, OPT_EMPTY); + + TEST_CHECK(A.nvals == (is_first_empty ? 0 : 25)); + TEST_CHECK(B.nvals == (is_second_empty ? 0 : 25)); + + if (is_first_empty && is_second_empty) { + if (!is_accum) { + TEST_CHECK(C.nvals == 0); + } else { + TEST_CHECK(C.nvals == (is_output_empty ? 0 : 25)); + } + } else { + TEST_CHECK(C.nvals == 25); + } + + CFL_matrix_free(&A); + CFL_matrix_free(&B); + CFL_matrix_free(&C); + } + } + } + } + + // if output equal first argument (iadd) + for (size_t is_first_empty = 0; is_first_empty < 2; is_first_empty++) { + for (size_t is_second_empty = 0; is_second_empty < 2; is_second_empty++) { + for (size_t is_accum = 0; is_accum < 2; is_accum++) { + A = is_first_empty ? CFL_matrix_create(5, 5) : make_ones_matrix(5); + B = is_second_empty ? CFL_matrix_create(5, 5) : make_ones_matrix(5); + + CFL_wise(&A, &A, &B, is_accum, OPT_EMPTY); + + TEST_CHECK(B.nvals == (is_second_empty ? 0 : 25)); + + if (is_second_empty) { + TEST_CHECK(A.nvals == (is_first_empty ? 0 : 25)); + } else { + TEST_CHECK(A.nvals == 25); + } + + CFL_matrix_free(&A); + CFL_matrix_free(&B); + } + } + } + + teardown(); +#endif +} + +static void test_CFL_empty_rsub_both_empty(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix A, B; + A = CFL_matrix_create(5, 5); + B = CFL_matrix_create(5, 5); + + OK(CFL_rsub(&A, &B, OPT_EMPTY)); + TEST_CHECK(A.nvals == 0); + + free_matrix(&A); + free_matrix(&B); + + teardown(); +#endif +} + +//------------------------------------------------------------------------------ +// Lazy optimization +//------------------------------------------------------------------------------ + +static void test_CFL_lazy_create(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + GrB_Matrix A_base; + GrB_Matrix_new(&A_base, GrB_BOOL, 5, 5); + Matrix A = CFL_matrix_from_base_lazy(A_base); + TEST_CHECK(A.base_matrices != NULL); + TEST_CHECK(A.base_matrices_count == 1); + TEST_CHECK(A.is_lazy == true); + TEST_CHECK(A.nvals == 0); + + free_matrix(&A); + + teardown(); +#endif +} + +// Create lazy matrix +// 1 1 1 0 .. 0 (count of 1: 111) +// 1 0 0 0 .. 0 (count of 1: 1) +// 1 1 0 0 .. 0 (count of 1: 11) +// . . . . .. . +// 0 0 0 0 .. 0 +static Matrix make_lazy_matrix(size_t base_matrices_count) { + size_t n = 0; + Matrix base_matrices[base_matrices_count]; + // 1 -> 1, 2 -> 11, 3 -> 111, 4 -> 1111 + for (size_t i = 0; i < base_matrices_count; i++) { + n *= 10; + n++; + } + + size_t tmp_n = 0; + for (size_t i = 0; i < base_matrices_count; i++) { + tmp_n *= 10; + tmp_n++; + GrB_Matrix base_matrix; + GrB_Matrix_new(&base_matrix, GrB_BOOL, n, n); + for (size_t j = 0; j < tmp_n; j++) { + // i+1 for future testing of sorting matrices + GrB_Matrix_setElement_BOOL(base_matrix, true, (i + 1) % base_matrices_count, + j); + } + + base_matrices[(i + 1) % base_matrices_count] = CFL_matrix_from_base(base_matrix); + } + + GrB_Matrix _result; + GrB_Matrix_new(&_result, GrB_BOOL, n, n); + Matrix result = CFL_matrix_from_base_lazy(_result); + GrB_free(&_result); + result.base = NULL; + + result.is_lazy = true; + result.base_matrices_count = base_matrices_count; + for (size_t i = 0; i < base_matrices_count; i++) { + result.base_matrices[i] = base_matrices[i]; + } + + CFL_matrix_update(&result); + + return result; +} + +static void test_CFL_lazy_mxm(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t is_accum = 0; is_accum < 2; is_accum++) { + for (size_t is_reverse = 0; is_reverse < 2; is_reverse++) { + Matrix A = make_lazy_matrix(4); + Matrix B = make_ones_matrix(1111); + Matrix C = CFL_matrix_create(1111, 1111); + + Matrix A_base = CFL_matrix_to_base(&A, OPT_LAZY); + Matrix B_base = CFL_matrix_to_base(&B, OPT_LAZY); + Matrix C_base = CFL_matrix_to_base(&C, OPT_LAZY); + + CFL_mxm(&C, &A, &B, is_accum, is_reverse, OPT_LAZY); + CFL_mxm(&C_base, &A_base, &B_base, is_accum, is_reverse, 0); + + TEST_CHECK(compare_matrices(C, C_base)); + + free_matrix(&A); + free_matrix(&B); + free_matrix(&C); + free_matrix(&A_base); + free_matrix(&B_base); + free_matrix(&C_base); + } + } + + for (size_t is_accum = 0; is_accum < 2; is_accum++) { + for (size_t is_reverse = 0; is_reverse < 2; is_reverse++) { + Matrix A = make_lazy_matrix(4); + Matrix B = CFL_matrix_create(1111, 1111); + Matrix C = CFL_matrix_create(1111, 1111); + + Matrix A_base = CFL_matrix_to_base(&A, OPT_LAZY); + Matrix B_base = CFL_matrix_to_base(&B, OPT_LAZY); + Matrix C_base = CFL_matrix_to_base(&C, OPT_LAZY); + + CFL_mxm(&C, &A, &B, is_accum, is_reverse, OPT_LAZY); + CFL_mxm(&C_base, &A_base, &B_base, is_accum, is_reverse, 0); + + TEST_CHECK(compare_matrices(C, C_base)); + + free_matrix(&A); + free_matrix(&B); + free_matrix(&C); + free_matrix(&A_base); + free_matrix(&B_base); + free_matrix(&C_base); + } + } + + teardown(); +#endif +} + +static void test_CFL_lazy_wise(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t is_accum = 0; is_accum < 2; is_accum++) { + Matrix A = make_lazy_matrix(4); + Matrix B = CFL_matrix_create(1111, 1111); + Matrix C = CFL_matrix_create(1111, 1111); + + Matrix A_base = CFL_matrix_to_base(&A, OPT_LAZY); + Matrix B_base = CFL_matrix_to_base(&B, OPT_LAZY); + Matrix C_base = CFL_matrix_to_base(&C, OPT_LAZY); + + CFL_wise(&C, &A, &B, is_accum, OPT_LAZY); + CFL_wise(&C_base, &A_base, &B_base, is_accum, 0); + + TEST_CHECK(compare_matrices(A, C_base)); + + free_matrix(&A); + free_matrix(&B); + free_matrix(&C); + free_matrix(&A_base); + free_matrix(&B_base); + free_matrix(&C_base); + } + + for (size_t is_accum = 0; is_accum < 2; is_accum++) { + Matrix A = CFL_matrix_create(1111, 1111); + Matrix B = make_lazy_matrix(4); + Matrix C = CFL_matrix_create(1111, 1111); + + Matrix A_base = CFL_matrix_to_base(&A, OPT_LAZY); + Matrix B_base = CFL_matrix_to_base(&B, OPT_LAZY); + Matrix C_base = CFL_matrix_to_base(&C, OPT_LAZY); + + CFL_wise(&C, &A, &B, is_accum, OPT_LAZY); + CFL_wise(&C_base, &A_base, &B_base, is_accum, 0); + + TEST_CHECK(compare_matrices(C, C_base)); + + free_matrix(&A); + free_matrix(&B); + free_matrix(&C); + free_matrix(&A_base); + free_matrix(&B_base); + free_matrix(&C_base); + } + + teardown(); +#endif +} + +static void test_CFL_lazy_rsub(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + { + Matrix A = make_lazy_matrix(4); + Matrix B = CFL_matrix_create(1111, 1111); + + Matrix A_base = CFL_matrix_to_base(&A, OPT_LAZY); + Matrix B_base = CFL_matrix_to_base(&B, OPT_LAZY); + + CFL_rsub(&A, &B, OPT_LAZY); + CFL_rsub(&A_base, &B_base, 0); + + TEST_CHECK(compare_matrices(A, A_base)); + + free_matrix(&A); + free_matrix(&B); + free_matrix(&A_base); + free_matrix(&B_base); + } + + { + Matrix A = CFL_matrix_create(1111, 1111); + Matrix B = make_lazy_matrix(4); + + Matrix A_base = CFL_matrix_to_base(&A, OPT_LAZY); + Matrix B_base = CFL_matrix_to_base(&B, OPT_LAZY); + + CFL_rsub(&A, &B, OPT_LAZY); + CFL_rsub(&A_base, &B_base, 0); + + TEST_CHECK(compare_matrices(A, A_base)); + + free_matrix(&A); + free_matrix(&B); + free_matrix(&A_base); + free_matrix(&B_base); + } + + teardown(); +#endif +} + +//------------------------------------------------------------------------------ +// Block optimization +//------------------------------------------------------------------------------ + +// Create block vector +static Matrix make_block_vector(size_t n, size_t alpha, Matrix_block format) { + GrB_Matrix result; + switch (format) { + case CELL: + // 1 1 0 0 0 0 . + // 0 1 1 0 0 0 . + // 0 0 1 1 0 0 . + // 0 0 0 1 1 0 . + // 0 0 0 0 1 1 . + // 1 0 0 0 0 1 . + // 1 1 0 0 0 0 . + // . . . . . . . + GrB_Matrix_new(&result, GrB_BOOL, n, n); + for (size_t i = 0; i < n; i++) { + for (size_t j = 0; j < n; j++) { + bool value = j == (i % n) || j == ((i + 1) % n) ? true : false; + if (!value) + continue; + + GrB_Matrix_setElement_BOOL(result, value, i, j); + } + } + + return CFL_matrix_from_base(result); + + case VEC_VERT: + GrB_Matrix_new(&result, GrB_BOOL, n * alpha, n); + for (size_t i = 0; i < n * alpha; i++) { + for (size_t j = 0; j < n; j++) { + bool value = j == (i + 2 % n) || j == ((i + 3) % n) ? true : false; + if (!value) + continue; + + GrB_Matrix_setElement_BOOL(result, value, i, j); + } + } + + return CFL_matrix_from_base(result); + + default: + break; + } + + // VEC_HORIZ: + GrB_Matrix_new(&result, GrB_BOOL, n, n * alpha); + for (size_t i = 0; i < n * alpha; i++) { + for (size_t j = 0; j < n; j++) { + bool value = j == (i + 5 % n) || j == ((i + 6) % n) ? true : false; + if (!value) + continue; + + GrB_Matrix_setElement_BOOL(result, value, i, j); + } + } + + return CFL_matrix_from_base(result); +} + +// TODO: create equality check with basic mxm +static void test_CFL_block_mxm(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t is_accum = 0; is_accum < 2; is_accum++) { + for (size_t is_reverse = 0; is_reverse < 2; is_reverse++) { + for (size_t first_format = 0; first_format < 3; first_format++) { + for (size_t second_format = 0; second_format < 3; second_format++) { + Matrix A = make_block_vector(20, 20, first_format); + Matrix B = make_block_vector(20, 20, second_format); + + Matrix C; + if (first_format != CELL && second_format != CELL) { + size_t nrows = first_format == VEC_HORIZ ? 20 : 20 * 20; + size_t ncols = first_format == VEC_HORIZ ? 20 * 20 : 20; + C = CFL_matrix_create(nrows, ncols); + } else { + C = CFL_matrix_create(first_format == CELL ? 20 : 20 * 20, + second_format == CELL ? 20 : 20 * 20); + } + + OK(CFL_mxm(&C, &A, &B, is_accum, is_reverse, OPT_BLOCK)); + + free_matrix(&A); + free_matrix(&B); + free_matrix(&C); + } + } + } + } + + teardown(); +#endif +} + +static void test_CFL_block_dup(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + { + Matrix A = make_block_vector(20, 20, VEC_HORIZ); + Matrix B = CFL_matrix_create(20 * 20, 20); + + OK(CFL_dup(&B, &A, OPT_BLOCK)); + TEST_CHECK(B.nvals == A.nvals); + + free_matrix(&A); + free_matrix(&B); + } + + { + Matrix A = make_block_vector(20, 20, VEC_VERT); + Matrix B = CFL_matrix_create(20 * 20, 20); + + OK(CFL_dup(&B, &A, OPT_BLOCK)); + TEST_CHECK(B.nvals == A.nvals); + + free_matrix(&A); + free_matrix(&B); + } + + { + Matrix A = make_block_vector(20, 20, VEC_HORIZ); + Matrix B = CFL_matrix_create(20, 20 * 20); + + OK(CFL_dup(&B, &A, OPT_BLOCK)); + TEST_CHECK(B.nvals == A.nvals); + + free_matrix(&A); + free_matrix(&B); + } + + { + Matrix A = make_block_vector(20, 20, VEC_VERT); + Matrix B = CFL_matrix_create(20, 20 * 20); + + OK(CFL_dup(&B, &A, OPT_BLOCK)); + TEST_CHECK(B.nvals == A.nvals); + + free_matrix(&A); + free_matrix(&B); + } + + teardown(); +#endif +} + +static void test_CFL_block_wise(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + { + Matrix A = make_block_vector(20, 20, VEC_HORIZ); + Matrix B = make_block_vector(20, 20, VEC_HORIZ); + Matrix C = make_block_vector(20, 20, VEC_HORIZ); + + TEST_CHECK(CFL_wise(&C, &A, &B, false, OPT_BLOCK) == GrB_INVALID_VALUE); + + free_matrix(&A); + free_matrix(&B); + free_matrix(&C); + } + + for (size_t is_accum = 0; is_accum < 2; is_accum++) { + for (size_t is_reverse = 0; is_reverse < 2; is_reverse++) { + for (size_t first_format = 0; first_format < 3; first_format++) { + for (size_t second_format = 0; second_format < 3; second_format++) { + Matrix A = make_block_vector(20, 20, first_format); + Matrix B = make_block_vector(20, 20, second_format); + + OK(CFL_wise(&A, &A, &B, is_accum, OPT_BLOCK)); + + free_matrix(&A); + free_matrix(&B); + } + } + } + } + + teardown(); +#endif +} + +static void test_CFL_block_rsub(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t is_accum = 0; is_accum < 2; is_accum++) { + for (size_t is_reverse = 0; is_reverse < 2; is_reverse++) { + for (size_t first_format = 0; first_format < 3; first_format++) { + for (size_t second_format = 0; second_format < 3; second_format++) { + Matrix A = make_block_vector(20, 20, first_format); + Matrix B = make_block_vector(20, 20, second_format); + + if ((first_format != CELL && second_format == CELL) || + (first_format == CELL && second_format != CELL)) { + TEST_CHECK(CFL_rsub(&A, &B, OPT_BLOCK) == GrB_INVALID_VALUE); + } else { + OK(CFL_rsub(&A, &B, OPT_BLOCK)); + } + + free_matrix(&A); + free_matrix(&B); + } + } + } + } + + teardown(); +#endif +} + +static void test_CFL_block_hyper_rotate(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + // lazy + { + GrB_Matrix _A; + GrB_Matrix_new(&_A, GrB_BOOL, 20 * 20, 20); + Matrix A = CFL_matrix_from_base_lazy(_A); + GrB_free(&_A); + + GrB_Matrix _A0; + GrB_Matrix_new(&_A0, GrB_BOOL, 20 * 20, 20); + Matrix A0 = CFL_matrix_from_base(_A0); + GrB_Matrix _A1; + GrB_Matrix_new(&_A1, GrB_BOOL, 20 * 20, 20); + Matrix A1 = CFL_matrix_from_base(_A1); + GrB_Matrix _A2; + GrB_Matrix_new(&_A2, GrB_BOOL, 20 * 20, 20); + Matrix A2 = CFL_matrix_from_base(_A2); + + A.base_matrices[0] = A0; + A.base_matrices[1] = A1; + A.base_matrices[2] = A2; + A.base_matrices_count = 3; + + OK(block_matrix_hyper_rotate_i(&A, VEC_HORIZ)); + TEST_CHECK(A.nrows == 20); + TEST_CHECK(A.ncols == 20 * 20); + + free_matrix(&A); + // free_matrix(&A0); + // free_matrix(&A1); + // free_matrix(&A2); + } + + // cell + { + Matrix A = CFL_matrix_create(5, 5); + + OK(block_matrix_hyper_rotate_i(&A, VEC_HORIZ)); + + free_matrix(&A); + } + + // VERT + { + Matrix A = make_block_vector(5, 5, VEC_VERT); + TEST_CHECK(A.nvals != 0); + TEST_CHECK(A.is_lazy != true); + + OK(block_matrix_hyper_rotate_i(&A, VEC_HORIZ)); + + free_matrix(&A); + } + + teardown(); +#endif +} + +static void test_CFL_block_to_diag(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix A = CFL_matrix_create(5, 5); + Matrix B = CFL_matrix_create(5, 5); + + GrB_Index nvals_before; + GrB_Matrix_nvals(&nvals_before, A.base); + block_matrix_to_diag(&A, &B); + GrB_Index nvals_after; + GrB_Matrix_nvals(&nvals_after, A.base); + + // A didn't change + TEST_CHECK(nvals_after == nvals_before); + + free_matrix(&A); + free_matrix(&B); + + teardown(); +#endif +} + +static void test_CFL_block_to_reduce(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix A = CFL_matrix_create(5, 5); + Matrix B = make_ones_matrix(5); + + // cell matrix cause just dup + block_matrix_reduce(&A, &B, OPT_BLOCK); + TEST_CHECK(A.nvals == B.nvals); + + free_matrix(&A); + free_matrix(&B); + + teardown(); +#endif +} + +//------------------------------------------------------------------------------ +// TEST LIST +//------------------------------------------------------------------------------ + +TEST_LIST = { + {"test_compare_matrices_function", test_compare_matrices_function}, + {"test CFL create and free", test_CFL_create_free}, + {"test CFL mxm", test_CFL_mxm}, + {"test CFL clear", test_CFL_clear}, + {"test CFL dup", test_CFL_dup}, + {"test_CFL_dup_same", test_CFL_dup_same}, + {"test_CFL_wise", test_CFL_wise}, + {"test_CFL_rsub", test_CFL_rsub}, + {"test_CFL_format_create_rowmajor_matrix", test_CFL_format_create_rowmajor_matrix}, + {"test_CFL_format_to_format", test_CFL_format_to_format}, + {"test_CFL_format_clear_format", test_CFL_format_clear_format}, + {"test_CFL_format_dup_format", test_CFL_format_dup_format}, + {"test_CFL_format_wise_when_both", test_CFL_format_wise_when_both}, + {"test_CFL_format_mxm_second_greather_then_k", + test_CFL_format_mxm_second_greather_then_k}, + {"test_CFL_empty_clear_empty_matrix", test_CFL_empty_clear_empty_matrix}, + {"test_CFL_empty_dup", test_CFL_empty_dup}, + {"test_CFL_empty_mxm_one_is_empty", test_CFL_empty_mxm_one_is_empty}, + {"test_CFL_empty_wise", test_CFL_empty_wise}, + {"test_CFL_empty_rsub_both_empty", test_CFL_empty_rsub_both_empty}, + {"test_CFL_lazy_create", test_CFL_lazy_create}, + {"test_CFL_lazy_mxm", test_CFL_lazy_mxm}, + {"test_CFL_lazy_wise", test_CFL_lazy_wise}, + {"test_CFL_lazy_rsub", test_CFL_lazy_rsub}, + {"test_CFL_block_mxm", test_CFL_block_mxm}, + {"test_CFL_block_wise", test_CFL_block_wise}, + {"test_CFL_block_rsub", test_CFL_block_rsub}, + {"test_CFL_block_dup", test_CFL_block_dup}, + {"test_CFL_block_hyper_rotate", test_CFL_block_hyper_rotate}, + {"test_CFL_block_to_diag", test_CFL_block_to_diag}, + {"test_CFL_block_to_reduce", test_CFL_block_to_reduce}, + {NULL, NULL}}; \ No newline at end of file diff --git a/experimental/test/test_CFL_reachability_multsrc.c b/experimental/test/test_CFL_reachability_multsrc.c new file mode 100644 index 0000000000..e62b989226 --- /dev/null +++ b/experimental/test/test_CFL_reachability_multsrc.c @@ -0,0 +1,839 @@ +//------------------------------------------------------------------------------ +// LAGraph/experimental/test/test_CFL_reachability_multsrc.c: test cases for +// Multiple-Source Context-Free Language Reachability Matrix-Based Algorithm +//------------------------------------------------------------------------------ +// +// LAGraph, (c) 2019-2025 by The LAGraph Contributors, All Rights Reserved. +// SPDX-License-Identifier: BSD-2-Clause + +// Contributed by Vlasenco Daniel, Ilhom Kombaev, Semyon Grigoriev, St. Petersburg State University. + +//------------------------------------------------------------------------------ + +/* + * Note: + * Tests than end with *_allsrc are the same as tests for CFL_reachability, with + * all vertices selected as sources. Their behaviour should be identical to CFL_reachability. +*/ + +#include +#include +#include +#include +#include +#include + +#define run_algorithm() \ + LAGraph_CFL_reachability_multsrc(&output, adj_matrices, src, src_count, \ + grammar.terms_count, grammar.nonterms_count, grammar.rules, grammar.rules_count, \ + msg) + +#define check_error(error) \ + { \ + retval = run_algorithm(); \ + TEST_CHECK(retval == error); \ + TEST_MSG("%d (retval) != %d (error) (%s)", retval, error, msg); \ + } + +#define check_result(result) \ + { \ + char *expected = output_to_str(0); \ + TEST_CHECK(strcmp(result, expected) == 0); \ + TEST_MSG("Wrong result. Actual: %s", expected); \ + } + +typedef struct { + size_t nonterms_count; + size_t terms_count; + size_t rules_count; + LAGraph_rule_WCNF *rules; +} grammar_t; + +GrB_Matrix *adj_matrices = NULL; +GrB_Matrix output = NULL; +grammar_t grammar = {0, 0, 0, NULL}; +char msg[LAGRAPH_MSG_LEN]; + +void setup() { LAGraph_Init(msg); } + +void teardown(void) { LAGraph_Finalize(msg); } + +char *output_to_str(size_t nonterm) { + GrB_Index nnz = 0; + OK(GrB_Matrix_nvals(&nnz, output)); + GrB_Index *row = NULL ; + GrB_Index *col = NULL ; + bool *val = NULL ; + LAGraph_Malloc ((void **) &row, nnz, sizeof (GrB_Index), msg) ; + LAGraph_Malloc ((void **) &col, nnz, sizeof (GrB_Index), msg) ; + LAGraph_Malloc ((void **) &val, nnz, sizeof (GrB_Index), msg) ; + + OK(GrB_Matrix_extractTuples(row, col, val, &nnz, output)); + + // 11 - size of " (%ld, %ld)" +// char *result_str = malloc(11 * nnz * sizeof(char)); + char *result_str = NULL ; + LAGraph_Malloc ((void **) &result_str, 11*nnz, sizeof (char), msg) ; + + result_str[0] = '\0'; + for (size_t i = 0; i < nnz; i++) { + sprintf(result_str + strlen(result_str), i == 0 ? + "(%" PRIu64 ", %" PRIu64 ")" : " (%" PRIu64 ", %" PRIu64 ")", + row[i], col[i]); + } + + LAGraph_Free ((void **) &row, msg); + LAGraph_Free ((void **) &col, msg); + LAGraph_Free ((void **) &val, msg); + + return result_str; +} + +void free_workspace() { + for (size_t i = 0; i < grammar.terms_count; i++) { + GrB_free(&adj_matrices[i]); + } + LAGraph_Free ((void **) &adj_matrices, msg); + + GrB_free(&output); + + LAGraph_Free ((void **) &grammar.rules, msg); + grammar = (grammar_t){0, 0, 0, NULL}; +} + +//================================ +// Grammars +//================================ + +// S -> aSb | ab in WCNF +// +// Terms: [0 a] [1 b] +// Nonterms: [0 S] [1 A] [2 B] [3 C] +// S -> AB [0 1 2 0] +// S -> AC [0 1 3 0] +// C -> SB [3 0 2 0] +// A -> a [1 0 -1 0] +// B -> b [2 1 -1 0] +void init_grammar_aSb() { +// LAGraph_rule_WCNF *rules = calloc(5, sizeof(LAGraph_rule_WCNF)); + LAGraph_rule_WCNF *rules = NULL ; + LAGraph_Calloc ((void **) &rules, 5, sizeof(LAGraph_rule_WCNF), msg); + + rules[0] = (LAGraph_rule_WCNF){0, 1, 2, 0}; + rules[1] = (LAGraph_rule_WCNF){0, 1, 3, 0}; + rules[2] = (LAGraph_rule_WCNF){3, 0, 2, 0}; + rules[3] = (LAGraph_rule_WCNF){1, 0, -1, 0}; + rules[4] = (LAGraph_rule_WCNF){2, 1, -1, 0}; + + grammar = (grammar_t){ + .nonterms_count = 4, .terms_count = 2, .rules_count = 5, .rules = rules}; +} + +// S -> aSb | ab | eps in WCNF +// +// Terms: [0 a] [1 b] +// Nonterms: [0 S] [1 A] [2 B] [3 C] +// S -> AB [0 1 2 0] +// S -> AC [0 1 3 0] +// S -> eps [0 -1 -1 0] +// C -> SB [3 0 2 0] +// A -> a [1 0 -1 0] +// B -> b [2 1 -1 0] +void init_grammar_aSb_eps() { +// LAGraph_rule_WCNF *rules = calloc(5, sizeof(LAGraph_rule_WCNF)); + LAGraph_rule_WCNF *rules = NULL ; + LAGraph_Calloc ((void **) &rules, 6, sizeof(LAGraph_rule_WCNF), msg); + + rules[0] = (LAGraph_rule_WCNF){0, 1, 2, 0}; + rules[1] = (LAGraph_rule_WCNF){0, 1, 3, 0}; + rules[2] = (LAGraph_rule_WCNF){3, 0, 2, 0}; + rules[3] = (LAGraph_rule_WCNF){1, 0, -1, 0}; + rules[4] = (LAGraph_rule_WCNF){2, 1, -1, 0}; + rules[5] = (LAGraph_rule_WCNF){0, -1, -1, 0}; + + + grammar = (grammar_t){ + .nonterms_count = 4, .terms_count = 2, .rules_count = 6, .rules = rules}; +} + +// S -> aS | a in WCNF +// +// Terms: [0 a] +// Nonterms: [0 S] +// S -> SS [0 0 0 0] +// S -> a [0 0 -1 0] +void init_grammar_aS() { +// LAGraph_rule_WCNF *rules = calloc(2, sizeof(LAGraph_rule_WCNF)); + LAGraph_rule_WCNF *rules = NULL ; + LAGraph_Calloc ((void **) &rules, 2, sizeof(LAGraph_rule_WCNF), msg); + + rules[0] = (LAGraph_rule_WCNF){0, 0, 0, 0}; + rules[1] = (LAGraph_rule_WCNF){0, 0, -1, 0}; + + grammar = (grammar_t){ + .nonterms_count = 1, .terms_count = 1, .rules_count = 2, .rules = rules}; +} + +// Complex grammar +// aaaabbbb or aaabbb +// +// Terms: [0 a] [1 b] +// Nonterms: [0 S] [n Sn] +// S -> S1 S2 [0 1 2 0] +// S -> S15 S16 [0 15 16 0] +// S1 -> S3 S4 [1 3 4 0] +// S2 -> S5 S6 [2 5 6 0] +// S3 -> S7 S8 [3 7 8 0] +// S4 -> S9 S10 [4 9 10 0] +// S5 -> S11 S12 [5 11 12 0] +// S6 -> S13 S14 [6 13 14 0] +// S16 -> S17 S18 [16 17 18 0] +// S17 -> S19 S20 [17 19 20 0] +// S18 -> S21 S22 [18 21 22 0] +// S22 -> S23 S24 [22 23 24 0] +// S7 -> a [7 0 -1 0] +// S8 -> a [8 0 -1 0] +// S9 -> a [9 0 -1 0] +// S10 -> a [10 0 -1 0] +// S11 -> b [11 1 -1 0] +// S12 -> b [12 1 -1 0] +// S13 -> b [13 1 -1 0] +// S14 -> b [14 1 -1 0] +// S15 -> a [15 0 -1 0] +// S19 -> a [19 0 -1 0] +// S20 -> a [20 0 -1 0] +// S21 -> b [21 1 -1 0] +// S23 -> b [23 1 -1 0] +// S24 -> b [24 1 -1 0] +void init_grammar_complex() { +// LAGraph_rule_WCNF *rules = calloc(26, sizeof(LAGraph_rule_WCNF)); + LAGraph_rule_WCNF *rules = NULL ; + LAGraph_Calloc ((void **) &rules, 26, sizeof(LAGraph_rule_WCNF), msg); + + rules[0] = (LAGraph_rule_WCNF){0, 1, 2, 0}; + rules[1] = (LAGraph_rule_WCNF){0, 15, 16, 0}; + rules[2] = (LAGraph_rule_WCNF){1, 3, 4, 0}; + rules[3] = (LAGraph_rule_WCNF){2, 5, 6, 0}; + rules[4] = (LAGraph_rule_WCNF){3, 7, 8, 0}; + rules[5] = (LAGraph_rule_WCNF){4, 9, 10, 0}; + rules[6] = (LAGraph_rule_WCNF){5, 11, 12, 0}; + rules[7] = (LAGraph_rule_WCNF){6, 13, 14, 0}; + rules[8] = (LAGraph_rule_WCNF){16, 17, 18, 0}; + rules[9] = (LAGraph_rule_WCNF){17, 19, 20, 0}; + rules[10] = (LAGraph_rule_WCNF){18, 21, 22, 0}; + rules[11] = (LAGraph_rule_WCNF){22, 23, 24, 0}; + rules[12] = (LAGraph_rule_WCNF){7, 0, -1, 0}; + rules[13] = (LAGraph_rule_WCNF){8, 0, -1, 0}; + rules[14] = (LAGraph_rule_WCNF){9, 0, -1, 0}; + rules[15] = (LAGraph_rule_WCNF){10, 0, -1, 0}; + rules[16] = (LAGraph_rule_WCNF){11, 1, -1, 0}; + rules[17] = (LAGraph_rule_WCNF){12, 1, -1, 0}; + rules[18] = (LAGraph_rule_WCNF){13, 1, -1, 0}; + rules[19] = (LAGraph_rule_WCNF){14, 1, -1, 0}; + rules[20] = (LAGraph_rule_WCNF){15, 0, -1, 0}; + rules[21] = (LAGraph_rule_WCNF){19, 0, -1, 0}; + rules[22] = (LAGraph_rule_WCNF){20, 0, -1, 0}; + rules[23] = (LAGraph_rule_WCNF){21, 1, -1, 0}; + rules[24] = (LAGraph_rule_WCNF){23, 1, -1, 0}; + rules[25] = (LAGraph_rule_WCNF){24, 1, -1, 0}; + + grammar = (grammar_t){ + .nonterms_count = 25, .terms_count = 2, .rules_count = 26, .rules = rules}; +} + +//================================ +// Graphs without vertex labels +//================================ + +// Graph: +// +// 0 -a-> 1 +// 1 -a-> 2 +// 2 -a-> 0 +// 0 -b-> 3 +// 3 -b-> 0 +void init_graph_double_cycle() { +// adj_matrices = calloc(2, sizeof(GrB_Matrix)); + LAGraph_Calloc ((void **) &adj_matrices, 2, sizeof (GrB_Matrix), msg) ; + + GrB_Matrix adj_matrix_a, adj_matrix_b; + OK(GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, 4, 4)); + OK(GrB_Matrix_new(&adj_matrix_b, GrB_BOOL, 4, 4)); + + OK(GrB_Matrix_setElement(adj_matrix_a, true, 0, 1)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 1, 2)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 2, 0)); + + OK(GrB_Matrix_setElement(adj_matrix_b, true, 0, 3)); + OK(GrB_Matrix_setElement(adj_matrix_b, true, 3, 0)); + + adj_matrices[0] = adj_matrix_a; + adj_matrices[1] = adj_matrix_b; +} + +// Graph: +// +// 0 -a-> 1 +// 1 -a-> 2 +// 2 -a-> 3 +// 3 -a-> 4 +// 3 -b-> 5 +// 4 -b-> 3 +// 5 -b-> 6 +// 6 -b-> 7 +void init_graph_1() { +// adj_matrices = calloc(2, sizeof(GrB_Matrix)); + LAGraph_Calloc ((void **) &adj_matrices, 2, sizeof (GrB_Matrix), msg) ; + + GrB_Matrix adj_matrix_a, adj_matrix_b; + OK(GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, 8, 8)); + OK(GrB_Matrix_new(&adj_matrix_b, GrB_BOOL, 8, 8)); + + OK(GrB_Matrix_setElement(adj_matrix_a, true, 0, 1)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 1, 2)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 2, 3)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 3, 4)); + + OK(GrB_Matrix_setElement(adj_matrix_b, true, 3, 5)); + OK(GrB_Matrix_setElement(adj_matrix_b, true, 4, 3)); + OK(GrB_Matrix_setElement(adj_matrix_b, true, 5, 6)); + OK(GrB_Matrix_setElement(adj_matrix_b, true, 6, 7)); + + adj_matrices[0] = adj_matrix_a; + adj_matrices[1] = adj_matrix_b; +} + +// Graph: +// +// 0 -a-> 2 +// 1 -a-> 2 +// 3 -a-> 5 +// 4 -a-> 5 +// 2 -a-> 6 +// 5 -a-> 6 +// 2 -b-> 0 +// 2 -b-> 1 +// 5 -b-> 3 +// 5 -b-> 4 +// 6 -b-> 2 +// 6 -b-> 5 +void init_graph_tree() { +// adj_matrices = calloc(2, sizeof(GrB_Matrix)); + LAGraph_Calloc ((void **) &adj_matrices, 2, sizeof (GrB_Matrix), msg) ; + + GrB_Matrix adj_matrix_a, adj_matrix_b; + OK(GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, 7, 7)); + OK(GrB_Matrix_new(&adj_matrix_b, GrB_BOOL, 7, 7)); + + OK(GrB_Matrix_setElement(adj_matrix_a, true, 0, 2)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 1, 2)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 3, 5)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 4, 5)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 2, 6)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 5, 6)); + + OK(GrB_Matrix_setElement(adj_matrix_b, true, 2, 0)); + OK(GrB_Matrix_setElement(adj_matrix_b, true, 2, 1)); + OK(GrB_Matrix_setElement(adj_matrix_b, true, 5, 3)); + OK(GrB_Matrix_setElement(adj_matrix_b, true, 5, 4)); + OK(GrB_Matrix_setElement(adj_matrix_b, true, 6, 2)); + OK(GrB_Matrix_setElement(adj_matrix_b, true, 6, 5)); + + adj_matrices[0] = adj_matrix_a; + adj_matrices[1] = adj_matrix_b; +} + +// Graph: +// +// 0 -a-> 1 +// 1 -a-> 2 +// 2 -a-> 0 +void init_graph_one_cycle() { +// adj_matrices = calloc(1, sizeof(GrB_Matrix)); + LAGraph_Calloc ((void **) &adj_matrices, 1, sizeof (GrB_Matrix), msg) ; + + GrB_Matrix adj_matrix_a; + GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, 3, 3); + + OK(GrB_Matrix_setElement(adj_matrix_a, true, 0, 1)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 1, 2)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 2, 0)); + + adj_matrices[0] = adj_matrix_a; +} + +// Graph: + +// 0 -a-> 1 +// 1 -a-> 2 +// 2 -b-> 3 +// 3 -b-> 4 +void init_graph_line() { +// adj_matrices = calloc(2, sizeof(GrB_Matrix)); + LAGraph_Calloc ((void **) &adj_matrices, 2, sizeof (GrB_Matrix), msg) ; + + GrB_Matrix adj_matrix_a, adj_matrix_b; + GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, 5, 5); + GrB_Matrix_new(&adj_matrix_b, GrB_BOOL, 5, 5); + + OK(GrB_Matrix_setElement(adj_matrix_a, true, 0, 1)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 1, 2)); + + OK(GrB_Matrix_setElement(adj_matrix_b, true, 2, 3)); + OK(GrB_Matrix_setElement(adj_matrix_b, true, 3, 4)); + + adj_matrices[0] = adj_matrix_a; + adj_matrices[1] = adj_matrix_b; +} + +// Graph: + +// 0 -a-> 0 +// 0 -b-> 1 +// 1 -c-> 2 +void init_graph_2() { +// adj_matrices = calloc(3, sizeof(GrB_Matrix)); + LAGraph_Calloc ((void **) &adj_matrices, 3, sizeof (GrB_Matrix), msg) ; + + GrB_Matrix adj_matrix_a, adj_matrix_b, adj_matrix_c; + GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, 3, 3); + GrB_Matrix_new(&adj_matrix_b, GrB_BOOL, 3, 3); + GrB_Matrix_new(&adj_matrix_c, GrB_BOOL, 3, 3); + + OK(GrB_Matrix_setElement(adj_matrix_a, true, 0, 0)); + OK(GrB_Matrix_setElement(adj_matrix_b, true, 0, 1)); + OK(GrB_Matrix_setElement(adj_matrix_c, true, 1, 2)); + + adj_matrices[0] = adj_matrix_a; + adj_matrices[1] = adj_matrix_b; + adj_matrices[2] = adj_matrix_c; +} + +// Graph: + +// 0 -a-> 1 +// 1 -a-> 0 +// 0 -b-> 0 +void init_graph_3() { +// adj_matrices = calloc(2, sizeof(GrB_Matrix)); + LAGraph_Calloc ((void **) &adj_matrices, 2, sizeof (GrB_Matrix), msg) ; + + GrB_Matrix adj_matrix_a, adj_matrix_b; + GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, 2, 2); + GrB_Matrix_new(&adj_matrix_b, GrB_BOOL, 2, 2); + + OK(GrB_Matrix_setElement(adj_matrix_a, true, 0, 1)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 1, 0)); + OK(GrB_Matrix_setElement(adj_matrix_b, true, 0, 0)); + + adj_matrices[0] = adj_matrix_a; + adj_matrices[1] = adj_matrix_b; +} + +// Graph: +// +// 1 -a-> 0 +// 2 -a-> 0 +// 3 -a-> 0 +// 4 -a-> 0 +// 5 -a-> 0 +// 1 -a-> 2 +// 2 -a-> 3 +// 3 -a-> 4 +// 4 -a-> 5 +// 5 -a-> 1 +void init_graph_whirlpool() { +// adj_matrices = calloc(2, sizeof(GrB_Matrix)); + LAGraph_Calloc ((void **) &adj_matrices, 1, sizeof (GrB_Matrix), msg) ; + + GrB_Matrix adj_matrix_a; + OK(GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, 6, 6)); + + OK(GrB_Matrix_setElement(adj_matrix_a, true, 1, 0)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 2, 0)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 3, 0)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 4, 0)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 5, 0)); + + OK(GrB_Matrix_setElement(adj_matrix_a, true, 1, 2)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 2, 3)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 3, 4)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 4, 5)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 5, 1)); + + adj_matrices[0] = adj_matrix_a; +} + +// Graph: +// +// 0 -a-> 1 +// 0 -a-> 2 +// 0 -a-> 3 +// 0 -a-> 4 +// 0 -b-> 0 +void init_graph_allout() { +// adj_matrices = calloc(2, sizeof(GrB_Matrix)); + LAGraph_Calloc ((void **) &adj_matrices, 1, sizeof (GrB_Matrix), msg) ; + + GrB_Matrix adj_matrix_a, adj_matrix_b; + OK(GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, 5, 5)); + + OK(GrB_Matrix_setElement(adj_matrix_a, true, 0, 1)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 0, 2)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 0, 3)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 0, 2)); + OK(GrB_Matrix_setElement(adj_matrix_a, true, 0, 4)); + + adj_matrices[0] = adj_matrix_a; +} + +//================================ +// Tests with all vertices as sources +//================================ + +void test_CFL_reachability_cycle_allsrc(void) { + setup(); + GrB_Info retval; + GrB_Index src[] = { 0, 1, 2 }; + int32_t src_count = sizeof(src) / sizeof(GrB_Index); + + init_grammar_aS(); + init_graph_one_cycle(); + + OK(run_algorithm()); + check_result("(0, 0) (0, 1) (0, 2) (1, 0) (1, 1) (1, 2) (2, 0) (2, 1) (2, 2)"); + + free_workspace(); + teardown(); +} + +void test_CFL_reachability_two_cycle_allsrc(void) { + setup(); + GrB_Info retval; + GrB_Index src[] = { 0, 1, 2, 3 }; + int32_t src_count = sizeof(src) / sizeof(GrB_Index); + + init_grammar_aSb(); + init_graph_double_cycle(); + + OK(run_algorithm()); + check_result("(0, 0) (0, 3) (1, 0) (1, 3) (2, 0) (2, 3)"); + + free_workspace(); + teardown(); +} + +void test_CFL_reachability_labels_more_than_nonterms_allsrc(void) { + setup(); + GrB_Info retval; + GrB_Index src[] = { 0, 1, 2 }; + int32_t src_count = sizeof(src) / sizeof(GrB_Index); + + init_grammar_aSb(); + init_graph_2(); + + OK(run_algorithm()); + check_result("(0, 1)"); + + free_workspace(); + teardown(); +} + +void test_CFL_reachability_complex_grammar_allsrc(void) { + setup(); + GrB_Info retval; + GrB_Index src[] = { 0, 1, 2, 3, 4, 5, 6, 7 }; + int32_t src_count = sizeof(src) / sizeof(GrB_Index); + + init_grammar_complex(); + init_graph_1(); + + OK(run_algorithm()); + check_result("(0, 7) (1, 6)"); + + free_workspace(); + teardown(); +} + +void test_CFL_reachability_tree_allsrc(void) { + setup(); + GrB_Info retval; + GrB_Index src[] = { 0, 1, 2, 3, 4, 5, 6 }; + int32_t src_count = sizeof(src) / sizeof(GrB_Index); + + init_grammar_aSb(); + init_graph_tree(); + + OK(run_algorithm()); + check_result("(0, 0) (0, 1) (0, 3) (0, 4) (1, 0) (1, 1) (1, 3) (1, 4) (2, 2) (2, 5) " + "(3, 0) (3, 1) (3, 3) (3, 4) (4, 0) (4, 1) (4, 3) (4, 4) (5, 2) (5, 5)"); + + free_workspace(); + teardown(); +} + +void test_CFL_reachability_line_allsrc(void) { + setup(); + GrB_Info retval; + GrB_Index src[] = { 0, 1, 2, 3, 4 }; + int32_t src_count = sizeof(src) / sizeof(GrB_Index); + + init_grammar_aSb(); + init_graph_line(); + + OK(run_algorithm()); + check_result("(0, 4) (1, 3)"); + + free_workspace(); + teardown(); +} + +void test_CFL_reachability_two_nodes_cycle_allsrc(void) { + setup(); + GrB_Info retval; + GrB_Index src[] = { 0, 1 }; + int32_t src_count = sizeof(src) / sizeof(GrB_Index); + + init_grammar_aSb(); + init_graph_3(); + + OK(run_algorithm()); + check_result("(0, 0) (1, 0)"); + + free_workspace(); + teardown(); +} + +//================================ +// Tests with some vertices as sources +//================================ + +void test_CFL_reachability_tree_msrc(void) { + setup(); + GrB_Info retval; + GrB_Index src[] = { 2, 3 }; + int32_t src_count = sizeof(src) / sizeof(GrB_Index); + + init_grammar_aSb(); + init_graph_tree(); + + OK(run_algorithm()); + check_result("(2, 2) (2, 5) (3, 0) (3, 1) (3, 3) (3, 4)"); + + free_workspace(); + teardown(); +} + +void test_CFL_reachability_allin_1_4(void) { + setup(); + GrB_Info retval; + GrB_Index src[] = { 1, 4 }; + int32_t src_count = sizeof(src) / sizeof(GrB_Index); + + init_grammar_aS(); + init_graph_whirlpool(); + + OK(run_algorithm()); + check_result( + "(1, 0) (1, 1) (1, 2) (1, 3) (1, 4) (1, 5) (4, 0) (4, 1) (4, 2) (4, 3) (4, 4) (4, 5)"); + + free_workspace(); + teardown(); +} + +//================================ +// Tests with one source vertex +//================================ + +void test_CFL_reachability_cycle_onesrc(void) { + setup(); + GrB_Info retval; + GrB_Index src[] = { 0 }; + int32_t src_count = sizeof(src) / sizeof(GrB_Index); + + init_grammar_aS(); + init_graph_one_cycle(); + + OK(run_algorithm()); + check_result("(0, 0) (0, 1) (0, 2)"); + + free_workspace(); + teardown(); +} + +void test_CFL_reachability_allin_1(void) { + setup(); + GrB_Info retval; + GrB_Index src[] = { 0 }; + int32_t src_count = sizeof(src) / sizeof(GrB_Index); + + init_grammar_aS(); + init_graph_whirlpool(); + + OK(run_algorithm()); + check_result(""); + + free_workspace(); + teardown(); +} + +void test_CFL_reachability_allout_0(void) { + setup(); + GrB_Info retval; + GrB_Index src[] = { 0 }; + int32_t src_count = sizeof(src) / sizeof(GrB_Index); + + init_grammar_aS(); + init_graph_allout(); + + OK(run_algorithm()); + check_result("(0, 1) (0, 2) (0, 3) (0, 4)"); + + free_workspace(); + teardown(); +} + +void test_CFL_reachability_allout_1(void) { + setup(); + GrB_Info retval; + GrB_Index src[] = { 1 }; + int32_t src_count = sizeof(src) / sizeof(GrB_Index); + + init_grammar_aS(); + init_graph_allout(); + + OK(run_algorithm()); + check_result(""); + + free_workspace(); + teardown(); +} + +//==================== +// Tests with invalid parameters +//==================== + +void test_CFL_reachability_invalid_rules(void) { + setup(); + GrB_Info retval; + + GrB_Index src[] = { 0 }; + int32_t src_count = sizeof(src) / sizeof(GrB_Index); + + init_grammar_aSb_eps(); + init_graph_double_cycle(); + + // Rule [Variable -> _ B] + grammar.rules[0] = + (LAGraph_rule_WCNF){.nonterm = 0, .prod_A = -1, .prod_B = 1, .index = 0}; + check_error(GrB_INVALID_VALUE); + + // Rule [_ -> A B] + grammar.rules[0] = + (LAGraph_rule_WCNF){.nonterm = -1, .prod_A = 1, .prod_B = 2, .index = 0}; + check_error(GrB_INVALID_VALUE); + + // Rule [C -> A B], where C >= nonterms_count + grammar.rules[0] = + (LAGraph_rule_WCNF){.nonterm = 10, .prod_A = 1, .prod_B = 2, .index = 0}; + check_error(GrB_INVALID_VALUE); + + // Rule [C -> t], where t >= terms_count + grammar.rules[0] = + (LAGraph_rule_WCNF){.nonterm = 0, .prod_A = 10, .prod_B = -1, .index = 0}; + check_error(GrB_INVALID_VALUE); + + free_workspace(); + teardown(); +} + +// At least one rule of form S -> _ and +// at least one source node should be present. +// nonterms_count, rules_count, src_count > 0 +void test_CFL_reachability_invalid_params(void) { + setup(); + GrB_Info retval; + + GrB_Index src[] = { 0, 1, 2 }; + + init_grammar_aSb(); + init_graph_double_cycle(); + + int32_t src_count = 0; + check_error(GrB_INVALID_VALUE); + src_count = sizeof(src) / sizeof(GrB_Index); + + int32_t temp = grammar.nonterms_count; + grammar.nonterms_count = 0; + check_error(GrB_INVALID_VALUE); + grammar.nonterms_count = temp; + + temp = grammar.rules_count; + grammar.rules_count = 0; + check_error(GrB_INVALID_VALUE); + grammar.rules_count = temp; + + free_workspace(); + teardown(); +} + +// output, rules, adj_matrices, src should not be NULL +void test_CFL_reachability_null_params(void) { + setup(); + GrB_Info retval; + + GrB_Index _src[] = { 0, 1, 2 }; + GrB_Index* src = _src; + int32_t src_count = sizeof(_src) / sizeof(GrB_Index); + + init_grammar_aSb(); + init_graph_double_cycle(); + + GrB_Matrix *temp_m = adj_matrices; + adj_matrices = NULL; + check_error(GrB_NULL_POINTER); + adj_matrices = temp_m; + + GrB_Index *temp_src = src; + src = NULL; + check_error(GrB_NULL_POINTER); + src = temp_src; + + LAGraph_rule_WCNF *temp_rules = grammar.rules; + grammar.rules = NULL; + check_error(GrB_NULL_POINTER); + grammar.rules = temp_rules; + + #define run_algorithm() \ + LAGraph_CFL_reachability_multsrc(NULL, adj_matrices, src, src_count, \ + grammar.terms_count, grammar.nonterms_count, grammar.rules, grammar.rules_count, \ + msg) + check_error(GrB_NULL_POINTER); + + free_workspace(); + teardown(); +} + + +TEST_LIST = { + // Tests from LAGraph_CFL_reachability as special cases + {"test_CFL_reachability_cycle_allsrc", test_CFL_reachability_cycle_allsrc}, + {"CFL_reachability_two_cycle_allsrc", test_CFL_reachability_two_cycle_allsrc}, + {"CFL_reachability_labels_more_than_nonterms_allsrc", test_CFL_reachability_labels_more_than_nonterms_allsrc}, + {"CFL_reachability_complex_grammar_allsrc", test_CFL_reachability_complex_grammar_allsrc}, + {"CFL_reachability_tree_allsrc", test_CFL_reachability_tree_allsrc}, + {"CFL_reachability_line_allsrc", test_CFL_reachability_line_allsrc}, + {"CFL_reachability_two_nodes_cycle_allsrc", test_CFL_reachability_two_nodes_cycle_allsrc}, + // Tests for some source vertices. + {"CFL_reachability_cycle_onesrc", test_CFL_reachability_cycle_onesrc}, + {"CFL_reachability_tree_msrc", test_CFL_reachability_tree_msrc}, + {"CFL_reachability_allin_1", test_CFL_reachability_allin_1}, + {"CFL_reachability_allin_1_4", test_CFL_reachability_allin_1_4}, + {"CFL_reachability_allout_0", test_CFL_reachability_allout_0}, + {"CFL_reachability_allout_1", test_CFL_reachability_allout_1}, + // Tests with invalid parameters + {"CFL_reachability_invalid_rules", test_CFL_reachability_invalid_rules}, + {"CFL_reachability_invalid_params", test_CFL_reachability_invalid_params}, + {"CFL_reachability_null_params", test_CFL_reachability_null_params}, + {NULL, NULL} + }; diff --git a/include/LAGraphX.h b/include/LAGraphX.h index 3355addb2e..223ddeb349 100644 --- a/include/LAGraphX.h +++ b/include/LAGraphX.h @@ -1096,6 +1096,124 @@ GrB_Info LAGraph_CFL_reachability char *msg // Message string for error reporting. ) ; +// LAGraph_CFL_reachability_multsrc: Multiple-Source Context-Free Language Reachability +// Matrix-Based Algorithm +// +// This function determines the set of vertex pairs (u, v) in a graph (represented by +// adjacency matrices) such that u is from the given set of source vertices, and +// there is a path from u to v, where the edge labels form a word from the language +// generated by the context-free grammar (represented by `rules`). +GrB_Info LAGraph_CFL_reachability_multsrc +( + // Output + GrB_Matrix *output, // A handle for the matrix containing results. + // + // output: (i, j) = true if and only if there is a path + // from node i to node j whose edge labels form a word + // derivable from the specified CFG. + // Input + const GrB_Matrix *adj_matrices, // Array of adjacency matrices representing the graph. + // The length of this array is equal to the count of + // terminals (terms_count). + // + // adj_matrices[t]: (i, j) == 1 if and only if there + // is an edge between nodes i and j with the label of + // the terminal corresponding to index 't' (where t is + // in the range [0, terms_count - 1]). + GrB_Index *src, // Array of source vertices + int32_t src_count, // The total number of source vertices + int32_t terms_count, // The total number of terminal symbols in the CFG. + int32_t nonterms_count, // The total number of non-terminal symbols in the CFG. + const LAGraph_rule_WCNF *rules, // The rules of the CFG. + size_t rules_count, // The total number of rules in the CFG. + char *msg // Message string for error reporting. +); + +// LAGraph_CFL_reachability_multsrc: Optimized Multiple-Source Context-Free +// Language Reachability Matrix-Based Algorithm +// +GrB_Info LAGraph_CFL_reachability_multsrc_fast +( + // Output + GrB_Matrix *output, // A handle for the matrix containing results. + // + // output: (i, j) = true if and only if there is a path + // from node i to node j whose edge labels form a word + // derivable from the specified CFG. + // Input + const GrB_Matrix *adj_matrices, // Array of adjacency matrices representing the graph. + // The length of this array is equal to the count of + // terminals (terms_count). + // + // adj_matrices[t]: (i, j) == 1 if and only if there + // is an edge between nodes i and j with the label of + // the terminal corresponding to index 't' (where t is + // in the range [0, terms_count - 1]). + GrB_Index *src, // Array of source vertices + // These parameters will become unsigned + int32_t src_count, // The total number of source vertices + int32_t terms_count, // The total number of terminal symbols in the CFG. + int32_t nonterms_count, // The total number of non-terminal symbols in the CFG. + const LAGraph_rule_WCNF *rules, // The rules of the CFG. + size_t rules_count, // The total number of rules in the CFG. + char *msg, // Message string for error reporting. + int8_t opt_mask // Optimizations mask +); + +enum CFL_Matrix_block { CELL, VEC_HORIZ, VEC_VERT }; + +// Struct for matrix for CFL algorithms with optimizations +typedef struct CFL_Matrix { + GrB_Matrix base; // Base GrB_Matrix from we create CFL_Matrix + int8_t optimizations; // Optimizations flags + // Fields of base matrix + GrB_Index nvals; + GrB_Index nrows; + GrB_Index ncols; + // Fields of format optimization + GrB_Matrix base_row; + GrB_Matrix base_col; + int32_t format; + bool is_both; + // Fields of lazy addition optimization + struct CFL_Matrix *base_matrices; + size_t base_matrices_count; + bool is_lazy; + // Fields of block optimization + enum CFL_Matrix_block block_type; +} CFL_Matrix; + +void CFL_matrix_update(CFL_Matrix *matrix); +CFL_Matrix CFL_matrix_from_base(GrB_Matrix matrix); +CFL_Matrix CFL_matrix_from_base_lazy(GrB_Matrix matrix); +CFL_Matrix CFL_matrix_create(GrB_Index nrows, GrB_Index ncols); +void CFL_matrix_free(CFL_Matrix *matrix); + +GrB_Info CFL_mxm(CFL_Matrix *output, CFL_Matrix *first, CFL_Matrix *second, bool accum, + bool swap, int8_t optimizations); +GrB_Info CFL_wise(CFL_Matrix *output, CFL_Matrix *first, CFL_Matrix *second, bool accum, + int8_t optimizations); +GrB_Info CFL_rsub(CFL_Matrix *output, CFL_Matrix *mask, int8_t optimizations); +GrB_Info CFL_dup(CFL_Matrix *output, CFL_Matrix *input, int8_t optimizations); + +// Converts a matrix (lazy or regular) into its evaluated base form by merging +// all underlying base matrices. If the input matrix is not lazy, returns its copy. +// +// Parameters: +// input - [in] Pointer to the source matrix (may be lazy). +// optimizations - [in] Bitmask specifying enabled optimizations. +// +// Returns: +// A new matrix representing the fully evaluated base form. +CFL_Matrix CFL_matrix_to_base +( + // input + CFL_Matrix *matrix, + int8_t optimizations +) ; + +GrB_Info CFL_clear(CFL_Matrix *A, int8_t optimizations); + //------------------------------------------------------------------------------ // a simple example of an algorithm //------------------------------------------------------------------------------ @@ -1520,4 +1638,4 @@ int LAGr_MaxFlow( } #endif -#endif +#endif \ No newline at end of file