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..553303fe2b --- /dev/null +++ b/experimental/algorithm/LAGraph_CFL_optimized_matrix_opt.c @@ -0,0 +1,1413 @@ +//------------------------------------------------------------------------------ +// 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 "LAGraph_CFL_optimized_matrix_opt.h" +#include "LG_internal.h" +#include + +#define BENCH_CFL_REACHBILITY false + +#if BENCH_CFL_REACHBILITY + #define IS_ISO(matrix, str) \ + { \ + bool iso_flag; \ + GrB_Index nnz; \ + TRY(GxB_Matrix_iso(&iso_flag, matrix)); \ + TRY(GrB_Matrix_nvals(&nnz, matrix)); \ + if (!iso_flag && nnz) { \ + printf("-----ISO ALERT----- (%s)\n", str); \ + TRY(GxB_print(matrix, 1)); \ + printf("-------------------\n"); \ + TRY(-122); \ + } \ + } + + #define IS_ROW(matrix, str) \ + { \ + if (matrix != NULL) { \ + int32_t orientation; \ + TRY(GrB_get(matrix, &orientation, GrB_STORAGE_ORIENTATION_HINT)); \ + if (orientation != GrB_ROWMAJOR) { \ + printf("-----NOT A ROW----- (%s)\n", str); \ + TRY(GxB_print(matrix, 1)); \ + printf("-------------------\n"); \ + TRY(-122); \ + } \ + } \ + } + + #define IS_COL(matrix, str) \ + { \ + if (matrix != NULL) { \ + int32_t orientation; \ + TRY(GrB_get(matrix, &orientation, GrB_STORAGE_ORIENTATION_HINT)); \ + if (orientation != GrB_COLMAJOR) { \ + printf("-----NOT A COL----- (%s)\n", str); \ + TRY(GxB_print(matrix, 1)); \ + printf("-------------------\n"); \ + TRY(-122); \ + } \ + } \ + } +#else + #define IS_ISO(matrix, str) + #define IS_ROW(matrix, str) + #define IS_COL(matrix, str) +#endif + +#define TRY(GrB_method) \ + { \ + GrB_Info LG_GrB_Info = GrB_method; \ + if (LG_GrB_Info < GrB_SUCCESS) { \ + fprintf(stderr, "LAGraph failure (file %s, line %d): (%d) \n", __FILE__, \ + __LINE__, LG_GrB_Info); \ + return (LG_GrB_Info); \ + } \ + } + +#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_print_lazy(Matrix *A, int8_t optimizations); + +GrB_Info matrix_to_format(Matrix *matrix, int32_t format, bool is_both) { + TRY(CFL_matrix_update(matrix)); + + // Matrix contain both formats so just switch base matrix + if (matrix->is_both) { + matrix->base = format == GrB_ROWMAJOR ? matrix->base_row : matrix->base_col; + matrix->format = format; + return GrB_SUCCESS; + } + + // No changes required + if (matrix->format == format) { + return GrB_SUCCESS; + } + + // Matrix contain just one matrix and format is not same + GrB_Matrix new_matrix, old_matrix; + if (matrix->format == GrB_ROWMAJOR) { + new_matrix = matrix->base_col; + old_matrix = matrix->base_row; + } else { + new_matrix = matrix->base_row; + old_matrix = matrix->base_col; + } + + if (is_both) { + TRY(GrB_Matrix_new(&new_matrix, GrB_BOOL, matrix->nrows, matrix->ncols)); + TRY(GrB_set(new_matrix, format, GrB_STORAGE_ORIENTATION_HINT)) + TRY(GrB_Matrix_assign(new_matrix, GrB_NULL, GrB_NULL, old_matrix, GrB_ALL, + matrix->nrows, GrB_ALL, matrix->ncols, GrB_NULL)); + if (format == GrB_ROWMAJOR) { + matrix->base_row = new_matrix; + matrix->base_col = matrix->base; + } else { + matrix->base_row = matrix->base; + matrix->base_col = new_matrix; + } + matrix->is_both = true; + matrix->base = new_matrix; + } else { + if (format == GrB_ROWMAJOR) { + matrix->base_row = matrix->base; + matrix->base_col = NULL; + TRY(TO_ROW(matrix->base)); + } else { + matrix->base_row = NULL; + matrix->base_col = matrix->base; + TRY(TO_COL(matrix->base)); + } + } + + matrix->format = format; + + TRY(CFL_matrix_update(matrix)); + + return GrB_SUCCESS; +} + +// clear methods + +GrB_Info matrix_clear(Matrix *A) { + TRY(GrB_Matrix_clear(A->base)); + TRY(CFL_matrix_update(A)); + return GrB_SUCCESS; +} + +GrB_Info matrix_clear_format(Matrix *A, int8_t optimizations) { + if (!(optimizations & OPT_FORMAT)) { + TRY(matrix_clear(A)); + return GrB_SUCCESS; + } + + if (!A->is_both) { + TRY(matrix_clear(A)); + return GrB_SUCCESS; + } + + TRY(matrix_to_format(A, GrB_ROWMAJOR, false)); + TRY(matrix_clear(A)); + + TRY(matrix_to_format(A, GrB_COLMAJOR, false)); + TRY(matrix_clear(A)); + + return GrB_SUCCESS; +} + +GrB_Info matrix_clear_empty(Matrix *A, int8_t optimizations) { + if (!(optimizations & OPT_EMPTY)) { + TRY(matrix_clear_format(A, optimizations)); + return GrB_SUCCESS; + } + + TRY(CFL_matrix_update(A)); + if (A->nvals == 0) { + return GrB_SUCCESS; + } + + TRY(matrix_clear_format(A, optimizations)); + return GrB_SUCCESS; +} + +// duplicate methods + +GrB_Info matrix_dup(Matrix *output, Matrix *input) { + if (output == input) { + return GrB_SUCCESS; + } + + TRY(GrB_Matrix_apply(output->base, GrB_NULL, GrB_NULL, GrB_IDENTITY_BOOL, input->base, + GrB_NULL)); + TRY(CFL_matrix_update(output)); + + return GrB_SUCCESS; +} + +GrB_Info matrix_dup_format(Matrix *output, Matrix *input, int8_t optimizations) { + if (!(optimizations & OPT_FORMAT)) { + TRY(matrix_dup(output, input)); + return GrB_SUCCESS; + } + + if (!output->is_both) { + TRY(CFL_matrix_update(output)); + TRY(CFL_matrix_update(input)); + Matrix *larger = output->nvals > input->nvals ? output : input; + + TRY(matrix_to_format(output, larger->format, false)); + TRY(matrix_to_format(input, larger->format, false)); + + TRY(matrix_dup(output, input)); + + return GrB_SUCCESS; + } + + TRY(matrix_to_format(output, GrB_ROWMAJOR, false)); + TRY(matrix_dup(output, input)); + TRY(matrix_to_format(output, GrB_COLMAJOR, false)); + TRY(matrix_dup(output, input)); + + return GrB_SUCCESS; +} + +GrB_Info matrix_dup_empty(Matrix *output, Matrix *input, int8_t optimizations) { + if (!(optimizations & OPT_EMPTY)) { + TRY(matrix_dup_format(output, input, optimizations)); + return GrB_SUCCESS; + } + + TRY(CFL_matrix_update(input)); + if (input->nvals == 0) { + TRY(matrix_clear_empty(output, optimizations)); + return GrB_SUCCESS; + } + + TRY(matrix_dup_format(output, input, optimizations)); + return GrB_SUCCESS; +} + +GrB_Info matrix_dup_lazy(Matrix *output, Matrix *input, int8_t optimizations) { + if (!(optimizations & OPT_LAZY)) { + TRY(matrix_dup_empty(output, input, optimizations)); + return GrB_SUCCESS; + } + + if (!input->is_lazy) { + TRY(matrix_dup_empty(output, input, optimizations)); + return GrB_SUCCESS; + } + + for (size_t i = 0; i < input->base_matrices_count; i++) { + TRY(matrix_dup_empty(output->base_matrices[i], input->base_matrices[i], + optimizations)); + } + output->base_matrices_count = input->base_matrices_count; + + return GrB_SUCCESS; +} + +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)) { + TRY(matrix_dup_lazy(output, input, optimizations)); + return GrB_SUCCESS; + } + + if (output->block_type == CELL && input->block_type == CELL) { + TRY(matrix_dup_lazy(output, input, optimizations)); + return GrB_SUCCESS; + } + + TRY(block_matrix_hyper_rotate_i(input, output->block_type)); + TRY(matrix_dup_lazy(output, input, optimizations)); + return GrB_SUCCESS; +} + +// 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++) { + TRY(block_matrix_hyper_rotate_i(matrix->base_matrices[i], format)); + } + + TRY(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; + TRY(GrB_Scalar_new(&scalar_true, GrB_BOOL)); + TRY(GrB_Scalar_setElement_BOOL(scalar_true, true)); + + TRY(CFL_matrix_update(matrix)); + + GrB_Index *nrows; + GrB_Index *ncols; + TRY(LAGraph_Calloc((void **)&nrows, matrix->nvals, sizeof(GrB_Index), NULL)); + TRY(LAGraph_Calloc((void **)&ncols, matrix->nvals, sizeof(GrB_Index), NULL)); + + TRY(GrB_Matrix_extractTuples_BOOL(nrows, ncols, NULL, &matrix->nvals, matrix->base)); + + if (matrix->block_type == VEC_VERT) { + 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; + } + } else { + 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; + TRY(GrB_Matrix_new(&new, GrB_BOOL, matrix->ncols, matrix->nrows)); + if (matrix->is_both) { + TRY(GxB_Matrix_build_Scalar(new, nrows, ncols, scalar_true, matrix->nvals)); + TRY(GrB_Matrix_free(&matrix->base_row)); + TRY(GrB_Matrix_free(&matrix->base_col)); + matrix->base = new; + matrix->base_row = new; + TRY(TO_ROW(matrix->base)); + TRY(GrB_Matrix_new(&new, GrB_BOOL, matrix->ncols, matrix->nrows)); + TRY(TO_COL(new)); + TRY(GxB_Matrix_build_Scalar(new, nrows, ncols, scalar_true, matrix->nvals)); + matrix->base_col = new; + matrix->format = GrB_ROWMAJOR; + } else { + TRY(GrB_Matrix_free(&matrix->base)); + matrix->base = new; + if (matrix->base_row != NULL) { + TRY(GxB_Matrix_build_Scalar(new, nrows, ncols, scalar_true, matrix->nvals)); + matrix->base_row = new; + } else { + TRY(TO_COL(new)); + TRY(GxB_Matrix_build_Scalar(new, nrows, ncols, scalar_true, matrix->nvals)); + matrix->base_col = new; + } + matrix->format = GrB_ROWMAJOR; + } + TRY(CFL_matrix_update(matrix)); + // TRY(CFL_matrix_free(matrix_p)); + // TRY(CFL_matrix_from_base(matrix_p, new)); + + TRY(LAGraph_Free((void **)&nrows, NULL)); + TRY(LAGraph_Free((void **)&ncols, NULL)); + TRY(GrB_free(&scalar_true)); + + return GrB_SUCCESS; +} + +GrB_Info block_matrix_to_diag(Matrix *diag, Matrix *input) { + if (input->block_type == CELL) { + return GrB_SUCCESS; + } + + GrB_Scalar scalar_true; + TRY(GrB_Scalar_new(&scalar_true, GrB_BOOL)); + TRY(GrB_Scalar_setElement_BOOL(scalar_true, true)); + + TRY(CFL_matrix_update(input)); + + GrB_Index *rows; + GrB_Index *cols; + TRY(LAGraph_Calloc((void **)&rows, input->nvals, sizeof(GrB_Index), NULL)); + TRY(LAGraph_Calloc((void **)&cols, input->nvals, sizeof(GrB_Index), NULL)); + TRY(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; + } + } + + TRY(GxB_Matrix_build_Scalar(diag->base, rows, cols, scalar_true, input->nvals)); + TRY(CFL_matrix_update(diag)); + + TRY(LAGraph_Free((void **)&rows, NULL)); + TRY(LAGraph_Free((void **)&cols, NULL)); + TRY(GrB_free(&scalar_true)); + + return GrB_SUCCESS; +} + +GrB_Info block_matrix_reduce(Matrix *matrix, Matrix *input, int8_t optimizations) { + if (input->block_type == CELL) { + TRY(matrix_dup_block(matrix, input, optimizations)); + return GrB_SUCCESS; + } + + GrB_Scalar scalar_true; + TRY(GrB_Scalar_new(&scalar_true, GrB_BOOL)); + TRY(GrB_Scalar_setElement_BOOL(scalar_true, true)); + + TRY(CFL_matrix_update(input)); + + GrB_Index *rows; + GrB_Index *cols; + TRY(LAGraph_Calloc((void **)&rows, input->nvals, sizeof(GrB_Index), NULL)); + TRY(LAGraph_Calloc((void **)&cols, input->nvals, sizeof(GrB_Index), NULL)); + TRY(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; + } + } + + TRY(GxB_Matrix_build_Scalar(matrix->base, rows, cols, scalar_true, input->nvals)); + TRY(CFL_matrix_update(matrix)); + + TRY(LAGraph_Free((void **)&rows, NULL)); + TRY(LAGraph_Free((void **)&cols, NULL)); + TRY(GrB_free(&scalar_true)); + + return GrB_SUCCESS; +} + +GrB_Info block_matrix_repeat_into_vector(Matrix *matrix, Matrix *input, + GrB_Index block_count) { + GrB_Matrix *tiles; + TRY(LAGraph_Calloc((void **)&tiles, block_count, sizeof(GrB_Matrix), NULL)); + for (size_t i = 0; i < block_count; i++) { + tiles[i] = input->base; + } + + if (matrix->block_type == VEC_VERT) { + TRY(GxB_Matrix_concat(matrix->base, tiles, block_count, 1, GrB_NULL)); + } else { + TRY(GxB_Matrix_concat(matrix->base, tiles, 1, block_count, GrB_NULL)); + } + TRY(CFL_matrix_update(matrix)); + TRY(LAGraph_Free((void **)&tiles, NULL)); + + return GrB_SUCCESS; +} + +// 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]; + TRY(CFL_matrix_update(first)); + TRY(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; +} + +GrB_Info CFL_matrix_to_base(Matrix **matrix_p, Matrix *input, int8_t optimizations) { + if (!input->is_lazy) { + Matrix *result; + TRY(CFL_matrix_create(&result, input->nrows, input->ncols)); + TRY(CFL_dup(result, input, optimizations)); + *matrix_p = result; + + return GrB_SUCCESS; + } + + Matrix *matrix; + TRY(CFL_matrix_create_lazy(&matrix, input->nrows, input->ncols)); + GrB_Info info = LAGraph_Realloc((void **)&matrix->base_matrices, + input->base_matrices_count, + matrix->base_matrices_count, sizeof(Matrix *), NULL); + if (info < GrB_SUCCESS) { + TRY(CFL_matrix_free(&matrix)); + return info; + } + TRY(CFL_matrix_free(&matrix->base_matrices[0])); + matrix->base_matrices_count = 0; + + for (size_t i = 0; i < input->base_matrices_count; i++) { + Matrix *base; + TRY(CFL_matrix_create(&base, input->nrows, input->ncols)); + matrix->base_matrices[i] = base; + matrix->base_matrices_count++; + } + + TRY(CFL_matrix_update(matrix)); + TRY(CFL_dup(matrix, input, optimizations)); + TRY(matrix_sort_lazy(matrix, false)); + + Matrix *acc; + TRY(CFL_matrix_create(&acc, input->nrows, input->ncols)); + for (size_t j = 0; j < matrix->base_matrices_count; j++) { + TRY(CFL_wise(acc, acc, matrix->base_matrices[j], false, optimizations)); + } + + TRY(CFL_matrix_free(&matrix)); + + *matrix_p = acc; + + return GrB_SUCCESS; +} + +GrB_Info CFL_matrix_extract_base(GrB_Matrix *base, Matrix *matrix, int8_t optimizations) { + if (matrix->is_lazy) { + // GraphBLAS never adds sparse matrices in-place: each addition builds a new + // matrix, so its cost depends on the size of both operands. Smaller base + // matrices are summed first (from the smallest one), and the largest base + // matrix is added only once at the end + TRY(matrix_sort_lazy(matrix, false)); + + size_t last = matrix->base_matrices_count - 1; + Matrix *largest = matrix->base_matrices[last]; + if (last > 0) { + Matrix *acc = matrix->base_matrices[0]; + for (size_t i = 1; i < last; i++) { + TRY(CFL_wise(acc, acc, matrix->base_matrices[i], false, optimizations)); + } + + TRY(CFL_wise(largest, largest, acc, false, optimizations)); + } + + TRY(CFL_matrix_extract_base(base, largest, optimizations)); + return GrB_SUCCESS; + } + + // Matrix contains both formats, so free the one that is not extracted + if (matrix->is_both) { + GrB_Matrix other = + matrix->base == matrix->base_row ? matrix->base_col : matrix->base_row; + TRY(GrB_Matrix_free(&other)); + } + + *base = matrix->base; + matrix->base = NULL; + matrix->base_row = NULL; + matrix->base_col = NULL; + matrix->is_both = false; + matrix->nvals = 0; + + return GrB_SUCCESS; +} + +GrB_Info matrix_combine_lazy(Matrix *A, size_t threshold, int8_t optimizations) { + Matrix **new_matrices; + TRY(LAGraph_Calloc((void **)&new_matrices, A->base_matrices_count, + sizeof(Matrix *), NULL)); + + size_t new_size = 0; + + TRY(matrix_sort_lazy(A, false)); + + for (size_t i = 0; i < A->base_matrices_count; i++) { + TRY(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; + } + + TRY(matrix_wise_empty(new_matrices[new_size - 1], new_matrices[new_size - 1], + A->base_matrices[i], false, optimizations)); + TRY(CFL_matrix_free(&A->base_matrices[i])); + } + + TRY(LAGraph_Free((void **)&A->base_matrices, NULL)); + A->base_matrices = new_matrices; + A->base_matrices_count = new_size; + TRY(CFL_matrix_update(A)); + + return GrB_SUCCESS; +} + +// create and update methods +GrB_Info CFL_matrix_update(Matrix *matrix) { + if (!matrix->is_lazy) { + TRY(GrB_Matrix_nvals(&matrix->nvals, matrix->base)); + } else { + size_t new_nnz = 0; + for (size_t i = 0; i < matrix->base_matrices_count; i++) { + TRY(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) { + TRY(GrB_Matrix_nrows(&matrix->nrows, matrix->base)); + TRY(GrB_Matrix_ncols(&matrix->ncols, matrix->base)); + } else { + TRY(GrB_Matrix_nrows(&matrix->nrows, matrix->base_matrices[0]->base)); + TRY(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; + +#ifdef BENCH_CFL_REACHBILITY + if (matrix->is_lazy) { + for (size_t i = 0; i < matrix->base_matrices_count; i++) { + IS_ISO(matrix->base_matrices[i]->base, "lazy matrix base"); + IS_ROW(matrix->base_matrices[i]->base_row, ""); + IS_COL(matrix->base_matrices[i]->base_col, ""); + } + } else { + IS_ISO(matrix->base, ""); + IS_ROW(matrix->base_row, ""); + IS_COL(matrix->base_col, ""); + } +#endif + + return GrB_SUCCESS; +} + +// Create optimizied matrix from base GrB_Matrix +// Input: base +// Output: matrix +GrB_Info CFL_matrix_from_base(Matrix **matrix, GrB_Matrix base) { + Matrix *result; + TRY(LAGraph_Calloc((void **)&result, 1, sizeof(Matrix), NULL)); + + result->base = base; + result->nvals = 0; // We will get actual info in update function + result->nrows = 0; + result->ncols = 0; + + // Format optimization fields + result->base_row = base; + 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; + + TRY(CFL_matrix_update(result)); + *matrix = result; + + return GrB_SUCCESS; +} + +// Create optimizied lazy matrix from base GrB_Matrix +// Input: base +// Output: matrix_p that holds new allocated matrix. +// +// Base matrix will be in array of lazy matrix +GrB_Info CFL_matrix_from_base_lazy(Matrix **matrix_p, GrB_Matrix base) { + Matrix *result; + TRY(CFL_matrix_from_base(&result, base)); + + TRY(CFL_matrix_from_base(matrix_p, base)); + + Matrix *matrix = *matrix_p; + matrix->is_lazy = true; + TRY(LAGraph_Calloc( + (void **)&matrix->base_matrices, 40, sizeof(CFL_Matrix *), + NULL)); // Initial storage; addition and copying resize it as needed. + matrix->base_matrices[0] = result; + matrix->base_matrices_count = 1; + matrix->base = NULL; + TRY(CFL_matrix_update(matrix)); + + return GrB_SUCCESS; +} + +GrB_Info CFL_matrix_create(Matrix **matrix, GrB_Index nrows, GrB_Index ncols) { + GrB_Matrix _result; + TRY(GrB_Matrix_new(&_result, GrB_BOOL, nrows, ncols)); + GrB_Info info = CFL_matrix_from_base(matrix, _result); + if (info < GrB_SUCCESS) { + GrB_Matrix_free(&_result); + } + + return info; +} + +GrB_Info CFL_matrix_create_lazy(Matrix **matrix, GrB_Index nrows, GrB_Index ncols) { + GrB_Matrix _result; + TRY(GrB_Matrix_new(&_result, GrB_BOOL, nrows, ncols)); + TRY(CFL_matrix_from_base_lazy(matrix, _result)) + + return GrB_SUCCESS; +} + +GrB_Info CFL_matrix_free(Matrix **matrix_p) { + if (*matrix_p == NULL) { + return GrB_SUCCESS; + } + + Matrix *matrix = *matrix_p; + if (matrix->is_both) { + TRY(GrB_Matrix_free(&(matrix->base_col))); + TRY(GrB_Matrix_free(&(matrix->base_row))); + } else { + TRY(GrB_Matrix_free(&(matrix->base))); + } + + for (size_t i = 0; i < matrix->base_matrices_count; i++) { + TRY(CFL_matrix_free(&(matrix->base_matrices[i]))); + } + TRY(LAGraph_Free((void **)&matrix->base_matrices, NULL)); + TRY(LAGraph_Free((void **)matrix_p, NULL)); + + return GrB_SUCCESS; + + 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; + + // CFL_matrix_update(first); + // CFL_matrix_update(second); + // CFL_matrix_update(output); + TRY(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"); + TRY(CFL_matrix_update(output)); + + return GrB_SUCCESS; +} + +GrB_Info matrix_mxm_format(Matrix *output, Matrix *first, Matrix *second, bool accum, + bool swap, int8_t optimizations) { + if (!(optimizations & OPT_FORMAT)) { + TRY(matrix_mxm(output, first, second, accum, swap)) + return GrB_SUCCESS; + } + + TRY(CFL_matrix_update(first)); + TRY(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)) { + TRY(matrix_mxm(output, first, second, accum, swap)) + return GrB_SUCCESS; + } + + TRY(matrix_to_format(first, desired_orientation, true)); + TRY(matrix_to_format(second, desired_orientation, false)); + TRY(matrix_to_format(output, desired_orientation, false)); + TRY(matrix_mxm(output, first, second, accum, swap)); + + return GrB_SUCCESS; +} + +GrB_Info matrix_mxm_empty(Matrix *output, Matrix *first, Matrix *second, bool accum, + bool swap, int8_t optimizations) { + if (!(optimizations & OPT_EMPTY)) { + TRY(matrix_mxm_format(output, first, second, accum, swap, optimizations)) + return GrB_SUCCESS; + } + + TRY(CFL_matrix_update(first)); + TRY(CFL_matrix_update(second)); + + if (first->nvals == 0 || second->nvals == 0) { + if (accum) { + return GrB_SUCCESS; + } + + TRY(CFL_matrix_update(output)); + if (output->nvals == 0) { + return GrB_SUCCESS; + } + + TRY(matrix_clear_empty(output, optimizations)); + return GrB_SUCCESS; + } + + TRY(matrix_mxm_format(output, first, second, accum, swap, optimizations)); + return GrB_SUCCESS; +} + +GrB_Info matrix_mxm_lazy(Matrix *output, Matrix *first, Matrix *second, bool accum, + bool swap, int8_t optimizations) { + if (!(optimizations & OPT_LAZY)) { + TRY(matrix_mxm_empty(output, first, second, accum, swap, optimizations)); + return GrB_SUCCESS; + } + + if (!first->is_lazy) { + TRY(matrix_mxm_empty(output, first, second, accum, swap, optimizations)); + return GrB_SUCCESS; + } + + TRY(CFL_matrix_update(second)); + + TRY(matrix_combine_lazy(first, second->nvals, optimizations)); + TRY(matrix_sort_lazy(first, false)); + + GrB_Matrix *accs; + TRY(LAGraph_Calloc((void **)&accs, first->base_matrices_count, sizeof(GrB_Matrix), + NULL)); + Matrix **acc_matrices; + TRY(LAGraph_Calloc((void **)&acc_matrices, first->base_matrices_count, + sizeof(Matrix *), NULL)); + for (size_t i = 0; i < first->base_matrices_count; i++) { + TRY(GrB_Matrix_new(&accs[i], GrB_BOOL, swap ? second->nrows : first->nrows, + swap ? first->ncols : second->ncols)); + TRY(CFL_matrix_from_base(&acc_matrices[i], accs[i])) + } + + for (size_t i = 0; i < first->base_matrices_count; i++) { + TRY(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_matrix, acc); + + for (size_t i = 0; i < first->base_matrices_count; i++) { + TRY(matrix_wise_empty(acc_matrix, acc_matrix, acc_matrices[i], false, + optimizations)); + TRY(CFL_matrix_free(&acc_matrices[i])); + } + LAGraph_Free((void **)&accs, NULL); + LAGraph_Free((void **)&acc_matrices, NULL); + + if (accum) { + TRY(matrix_wise_empty(output, output, acc_matrix, false, optimizations)); + } else { + TRY(matrix_dup_block(output, acc_matrix, optimizations)); + } + + TRY(CFL_matrix_free(&acc_matrix)); + + return GrB_SUCCESS; +} + +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)) { + TRY(matrix_mxm_lazy(output, first, second, accum, swap, optimizations)); + return GrB_SUCCESS; + } + + if (first->block_type == CELL && second->block_type == CELL) { + if (output->block_type == CELL) { + TRY(matrix_mxm_lazy(output, first, second, accum, swap, optimizations)); + return GrB_SUCCESS; + } + + // N_i -> A B: compute the square product and add it to every output block. + Matrix *temp; + TRY(CFL_matrix_create(&temp, swap ? second->nrows : first->nrows, + swap ? first->ncols : second->ncols)); + GrB_Info info = matrix_mxm_lazy(temp, first, second, false, swap, optimizations); + if (info >= GrB_SUCCESS && !accum) { + info = matrix_clear_empty(output, optimizations); + } + if (info >= GrB_SUCCESS) { + info = matrix_wise_block(output, output, temp, false, optimizations); + } + TRY(CFL_matrix_free(&temp)); + return info; + } + + if (first->block_type == CELL) { + TRY(block_matrix_hyper_rotate_i(second, swap ? VEC_VERT : VEC_HORIZ)); + TRY(block_matrix_hyper_rotate_i(output, swap ? VEC_VERT : VEC_HORIZ)); + + Matrix *temp; + TRY(CFL_matrix_create(&temp, swap ? second->nrows : first->nrows, + swap ? first->ncols : second->ncols)); + TRY(matrix_mxm_lazy(temp, first, second, accum, swap, optimizations)); + TRY(matrix_wise_block(output, output, temp, false, optimizations)); + TRY(CFL_matrix_free(&temp)); + + return GrB_SUCCESS; + } + + if (second->block_type == CELL) { + TRY(block_matrix_hyper_rotate_i(first, swap ? VEC_HORIZ : VEC_VERT)); + TRY(block_matrix_hyper_rotate_i(output, swap ? VEC_HORIZ : VEC_VERT)); + + Matrix *temp; + TRY(CFL_matrix_create(&temp, swap ? second->nrows : first->nrows, + swap ? first->ncols : first->ncols)); + TRY(matrix_mxm_lazy(temp, first, second, accum, swap, optimizations)); + TRY(matrix_wise_block(output, output, temp, false, optimizations)); + TRY(CFL_matrix_free(&temp)); + + return GrB_SUCCESS; + } + + GrB_Index size = first->nrows > first->ncols ? first->nrows : first->ncols; + + Matrix *diag; + TRY(CFL_matrix_create(&diag, size, size)); + TRY(block_matrix_to_diag(diag, second)); + TRY(CFL_matrix_update(diag)); + + TRY(block_matrix_hyper_rotate_i(first, swap ? VEC_VERT : VEC_HORIZ)); + TRY(block_matrix_hyper_rotate_i(output, swap ? VEC_VERT : VEC_HORIZ)); + + Matrix *temp; + TRY(CFL_matrix_create(&temp, first->nrows, first->ncols)); + TRY(matrix_mxm_lazy(temp, first, diag, false, swap, optimizations)); + TRY(matrix_wise_block(output, output, temp, false, optimizations)); + + TRY(CFL_matrix_free(&temp)); + TRY(CFL_matrix_free(&diag)); + + return GrB_SUCCESS; +} + +// 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; + + TRY(CFL_matrix_update(first)); + TRY(CFL_matrix_update(second)); + TRY(CFL_matrix_update(output)); + + TRY(GrB_Matrix_eWiseAdd_BinaryOp(output->base, GrB_NULL, accum_op, GxB_ANY_BOOL, + first->base, second->base, GrB_NULL)); + + TRY(CFL_matrix_update(output)); + + return GrB_SUCCESS; +} + +GrB_Info matrix_wise_format(Matrix *output, Matrix *first, Matrix *second, bool accum, + int8_t optimizations) { + if (!(optimizations & OPT_FORMAT)) { + TRY(matrix_wise(output, first, second, accum)) + return GrB_SUCCESS; + } + + if (!output->is_both) { + TRY(CFL_matrix_update(output)); + TRY(CFL_matrix_update(first)); + TRY(CFL_matrix_update(second)); + + Matrix *larger = output->nvals > first->nvals ? output : first; + larger = larger->nvals > second->nvals ? larger : second; + + TRY(matrix_to_format(output, larger->format, false)); + TRY(matrix_to_format(first, larger->format, false)); + TRY(matrix_to_format(second, larger->format, false)); + + TRY(matrix_wise(output, first, second, accum)); + + return GrB_SUCCESS; + } + + TRY(matrix_to_format(output, GrB_ROWMAJOR, false)); + TRY(matrix_to_format(first, output->format, false)); + TRY(matrix_to_format(second, output->format, false)); + + TRY(matrix_wise(output, first, second, accum)); + + TRY(matrix_to_format(output, GrB_COLMAJOR, false)); + TRY(matrix_to_format(first, output->format, false)); + TRY(matrix_to_format(second, output->format, false)); + + TRY(matrix_wise(output, first, second, accum)); + + return GrB_SUCCESS; +} + +GrB_Info matrix_wise_empty(Matrix *output, Matrix *first, Matrix *second, bool accum, + int8_t optimizations) { + if (!(optimizations & OPT_EMPTY)) { + TRY(matrix_wise_format(output, first, second, accum, optimizations)); + return GrB_SUCCESS; + } + + TRY(CFL_matrix_update(first)); + TRY(CFL_matrix_update(second)); + TRY(CFL_matrix_update(output)); + + if (output == first) { + if (first->nvals == 0) { + TRY(matrix_dup_empty(first, second, optimizations)); + return GrB_SUCCESS; + } + + if (second->nvals == 0) { + return GrB_SUCCESS; + } + + TRY(matrix_wise_format(output, first, second, accum, optimizations)) + return GrB_SUCCESS; + } + + if (first->nvals == 0 && second->nvals == 0) { + if (accum || output->nvals == 0) { + return GrB_SUCCESS; + } + + TRY(matrix_clear_empty(output, optimizations)); + return GrB_SUCCESS; + } + + if (first->nvals == 0) { + if (accum) { + TRY(matrix_wise_empty(output, output, second, false, optimizations)); + return GrB_SUCCESS; + } + + TRY(matrix_dup_empty(output, second, optimizations)); + return GrB_SUCCESS; + } + + if (second->nvals == 0) { + if (accum) { + TRY(matrix_wise_empty(output, output, first, false, optimizations)); + return GrB_SUCCESS; + } + + TRY(matrix_dup_empty(output, first, optimizations)); + return GrB_SUCCESS; + } + + TRY(matrix_wise_format(output, first, second, accum, optimizations)); + return GrB_SUCCESS; +} + +GrB_Info matrix_wise_lazy(Matrix *output, Matrix *first, Matrix *second, bool accum, + int8_t optimizations) { + if (!(optimizations & OPT_LAZY)) { + TRY(matrix_wise_empty(output, first, second, accum, optimizations)) + return GrB_SUCCESS; + } + + if (!first->is_lazy && !second->is_lazy) { + TRY(matrix_wise_empty(output, first, second, accum, optimizations)); + return GrB_SUCCESS; + } + + if (!first->is_lazy && second->is_lazy) { + for (size_t i = 0; i < second->base_matrices_count; i++) { + TRY(matrix_wise_empty(output, first, second->base_matrices[i], true, + optimizations)); + } + + return GrB_SUCCESS; + } + + Matrix *other; + TRY(CFL_matrix_create(&other, output->nrows, output->ncols)); + TRY(matrix_dup_empty(other, second, optimizations)); + + // Reserve before merging so allocation failure leaves the components intact. + GrB_Info info = LAGraph_Realloc((void **)&first->base_matrices, + first->base_matrices_count + 1, + first->base_matrices_count, sizeof(Matrix *), NULL); + if (info < GrB_SUCCESS) { + TRY(CFL_matrix_free(&other)); + return info; + } + + while (true) { + // A merge changes other's size and can make another component comparable. + size_t other_nvals = other->nvals >= 10 ? other->nvals : 10; + bool found = false; + + for (size_t i = 0; i < first->base_matrices_count; i++) { + TRY(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) { + TRY(matrix_wise_empty(other, other, first->base_matrices[i], accum, + optimizations)); + TRY(CFL_matrix_free(&first->base_matrices[i])); + 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; + } + + TRY(matrix_sort_lazy(first, false)); + TRY(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)) { + TRY(matrix_wise_lazy(output, first, second, accum, optimizations)); + return GrB_SUCCESS; + } + + 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) { + TRY(matrix_wise_lazy(output, first, second, accum, optimizations)); + return GrB_SUCCESS; + } + + // second is vector + if (first->block_type == CELL) { + Matrix *temp_reduced; + TRY(CFL_matrix_create(&temp_reduced, first->nrows, first->ncols)); + TRY(block_matrix_reduce(temp_reduced, second, optimizations)); + + TRY(matrix_wise_lazy(output, first, temp_reduced, accum, optimizations)); + TRY(CFL_matrix_free(&temp_reduced)); + + return GrB_SUCCESS; + } + + // first is vector + if (second->block_type == CELL) { + // LG_SET_BURBLE(true); + Matrix *temp_vector; + TRY(CFL_matrix_create(&temp_vector, first->nrows, first->ncols)); + + GrB_Index block_count = 0; + if (first->nrows > first->ncols) { + block_count = first->nrows / first->ncols; + } else { + block_count = first->ncols / first->nrows; + } + TRY(block_matrix_repeat_into_vector(temp_vector, second, block_count)); + TRY(block_matrix_hyper_rotate_i(temp_vector, first->block_type)); + TRY(matrix_wise_lazy(output, first, temp_vector, accum, optimizations)); + + TRY(CFL_matrix_free(&temp_vector)); + + return GrB_SUCCESS; + } + + // both are vector + TRY(block_matrix_hyper_rotate_i(second, first->block_type)); + TRY(matrix_wise_lazy(output, first, second, accum, optimizations)); + + return GrB_SUCCESS; +} + +// rsub methods + +GrB_Info matrix_rsub(Matrix *output, Matrix *mask) { + TRY(GrB_eWiseAdd(output->base, mask->base, GrB_NULL, GxB_ANY_BOOL, output->base, + output->base, GrB_DESC_RSC)); + + TRY(CFL_matrix_update(output)); + + return GrB_SUCCESS; +} + +GrB_Info matrix_rsub_format(Matrix *output, Matrix *mask, int8_t optimizations) { + if (!(optimizations & OPT_FORMAT)) { + TRY(matrix_rsub(output, mask)); + return GrB_SUCCESS; + } + + TRY(CFL_matrix_update(output)); + TRY(CFL_matrix_update(mask)); + + Matrix *larger_matrix = output->nvals > mask->nvals ? output : mask; + TRY(matrix_to_format(output, larger_matrix->format, false)); + TRY(matrix_to_format(mask, larger_matrix->format, false)); + + TRY(matrix_rsub(output, mask)); + + return GrB_SUCCESS; +} + +GrB_Info matrix_rsub_empty(Matrix *output, Matrix *mask, int8_t optimizations) { + if (!(optimizations & OPT_EMPTY)) { + TRY(matrix_rsub_format(output, mask, optimizations)); + return GrB_SUCCESS; + } + + TRY(CFL_matrix_update(mask)); + TRY(CFL_matrix_update(output)); + if (mask->nvals == 0 || output->nvals == 0) { + return GrB_SUCCESS; + } + + TRY(matrix_rsub_format(output, mask, optimizations)); + + return GrB_SUCCESS; +} + +GrB_Info matrix_rsub_lazy(Matrix *output, Matrix *mask, int8_t optimizations) { + if (!(optimizations & OPT_LAZY)) { + TRY(matrix_rsub_empty(output, mask, optimizations)); + return GrB_SUCCESS; + } + + if (!mask->is_lazy) { + TRY(matrix_rsub_empty(output, mask, optimizations)) + return GrB_SUCCESS; + } + + TRY(CFL_matrix_update(output)); + TRY(matrix_combine_lazy(mask, output->nvals, optimizations)); + TRY(matrix_sort_lazy(mask, true)); + + for (size_t i = 0; i < mask->base_matrices_count; i++) { + TRY(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)) { + TRY(matrix_rsub_lazy(output, mask, optimizations)) + return GrB_SUCCESS; + } + + 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) { + TRY(matrix_rsub_lazy(output, mask, optimizations)); + return GrB_SUCCESS; + } + + TRY(block_matrix_hyper_rotate_i(output, mask->block_type)); + TRY(matrix_rsub_lazy(output, mask, optimizations)); + + return GrB_SUCCESS; +} + +// utility methods + +// GrB_Info matrix_print_lazy(Matrix *A, int8_t optimizations) { +// // return; +// GxB_Print_Level pr = 1; + +// if (!A->is_lazy) { +// TRY(CFL_matrix_update(A)); +// TRY(GxB_print(A->base, pr)); +// // printf("nnz: %ld\n", A->nvals); +// return GrB_SUCCESS; +// } + +// if (A->base_matrices_count == 1) { +// TRY(CFL_matrix_update(A)); +// // printf("nnz: %ld\n", A->nvals); +// TRY(GxB_print(A->base_matrices[0]->base, pr)); +// return GrB_SUCCESS; +// } + +// Matrix *temp; +// TRY(CFL_matrix_create(&temp, A->nrows, A->ncols)); +// for (size_t i = 0; i < A->base_matrices_count; i++) { +// TRY(matrix_wise_empty(temp, temp, A->base_matrices[i], false, optimizations)); +// } + +// A = temp; +// TRY(GxB_print(A->base, pr)); +// TRY(CFL_matrix_update(A)); +// // printf("nnz: %ld\n", A->nvals); +// TRY(CFL_matrix_free(&temp)); + +// return GrB_SUCCESS; +// } + +// GrB_Info print_graph_info(Matrix **matrices, size_t count) { +// GrB_Index nnz = 0; + +// for (size_t i = 0; i < count; i++) { +// Matrix *A = matrices[i]; +// TRY(CFL_matrix_update(A)); +// nnz += A->nvals; +// } + +// printf("NNZ: %ld\n", nnz); +// return GrB_SUCCESS; +// } + +// order of optimizations: block -> lazy -> empty -> format + +GrB_Info CFL_mxm(Matrix *output, Matrix *first, Matrix *second, bool accum, bool swap, + int8_t optimizations) { + TRY(CFL_matrix_update(first)); + TRY(CFL_matrix_update(second)); + TRY(CFL_matrix_update(output)); + + TRY(matrix_mxm_block(output, first, second, accum, swap, optimizations)); + + TRY(CFL_matrix_update(output)); + return GrB_SUCCESS; +} + +GrB_Info CFL_wise(Matrix *output, Matrix *first, Matrix *second, bool accum, + int8_t optimizations) { + TRY(CFL_matrix_update(first)); + TRY(CFL_matrix_update(second)); + TRY(CFL_matrix_update(output)); + + TRY(matrix_wise_block(output, first, second, accum, optimizations)); + + TRY(CFL_matrix_update(output)); + + return GrB_SUCCESS; +} + +GrB_Info CFL_rsub(Matrix *output, Matrix *mask, int8_t optimizations) { + TRY(CFL_matrix_update(mask)); + TRY(CFL_matrix_update(output)); + + TRY(matrix_rsub_block(output, mask, optimizations)) + + TRY(CFL_matrix_update(output)); + return GrB_SUCCESS; +} + +GrB_Info CFL_dup(Matrix *output, Matrix *input, int8_t optimizations) { + TRY(CFL_matrix_update(output)); + TRY(CFL_matrix_update(input)); + + TRY(matrix_dup_block(output, input, optimizations)) + + TRY(CFL_matrix_update(output)); + return GrB_SUCCESS; +} + +GrB_Info CFL_clear(Matrix *A, int8_t optimizations) { + TRY(CFL_matrix_update(A)); + + TRY(matrix_clear_empty(A, optimizations)) + + TRY(CFL_matrix_update(A)); + return GrB_SUCCESS; +} diff --git a/experimental/algorithm/LAGraph_CFL_optimized_matrix_opt.h b/experimental/algorithm/LAGraph_CFL_optimized_matrix_opt.h new file mode 100644 index 0000000000..c80bc4aee9 --- /dev/null +++ b/experimental/algorithm/LAGraph_CFL_optimized_matrix_opt.h @@ -0,0 +1,121 @@ +//------------------------------------------------------------------------------ +// LAGraph_CFL_optimized_matrix_opt.h: Header with operations for +// Optimized Context-Free Language Reachability Matrix-Based Algorithm +//------------------------------------------------------------------------------ +// +// LAGraph, (c) 2019-2026 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 + +// Several square matrices can be concatenated into a horizontal or vertical vector of +// matrices for fast multiplication +// +// [CELL] - the matrix is square: n x n +// [VEC_HORIZ] - the matrix is a horizontal block vector: n x (n * k), +// where k is the number of concatenated matrices +// [VEC_VERT] - the matrix is a vertical block vector: (n * k) x n, +// where k is the number of concatenated matrices +enum CFL_Matrix_block { CELL, VEC_HORIZ, VEC_VERT }; + +// Matrix wrapper for CFL reachability algorithms with optional optimizations +typedef struct CFL_Matrix { + GrB_Matrix base; // Underlying GrB_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; + +// Creates a CFL matrix from intialized GrB_Matrix +// +// Creates a CFL_Matrix wrapping an already-initialized GrB_Matrix +// On success, *matrix is set to a newly allocated CFL_Matrix +GrB_Info CFL_matrix_from_base(CFL_Matrix **matrix, GrB_Matrix base); + +// Creates a lazy CFL matrix from intialized GrB_Matrix +// +// Creates a CFL_Matrix wrapping an already-initialized GrB_Matrix +// On success, *matrix is set to a newly allocated CFL_Matrix +GrB_Info CFL_matrix_from_base_lazy(CFL_Matrix **matrix, GrB_Matrix base); + +// Creates an empty CFL_Matrix of the given dimensions +// +// The underlying GrB_Matrix is allocated internally +// On success, *matrix is set to a newly allocated CFL_Matrix +GrB_Info CFL_matrix_create(CFL_Matrix **matrix, GrB_Index nrows, GrB_Index ncols); + +// Creates an empty lazy CFL_Matrix of the given dimensions +// +// The underlying GrB_Matrix is allocated internally +// On success, *matrix is set to a newly allocated CFL_Matrix with new CFL_Matrix +GrB_Info CFL_matrix_create_lazy(CFL_Matrix **matrix, GrB_Index nrows, GrB_Index ncols); + +// Free a CFL_Matrix and sets *matrix to NULL +GrB_Info CFL_matrix_free(CFL_Matrix **matrix); + +// Recomputes internal fields (nvals, nrows, ncols, format, block_type, etc...) +// from the current state of the underlying base matrix +GrB_Info CFL_matrix_update(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); + +// Computes C = B \ A, i.e. C(i, j) = true iff B(i, j) = true and A(i, j) = false. +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: +// matrix_p - [out] Pointer that will be holds created matrix +// input - [in] Pointer to the source matrix (may be lazy). +// optimizations - [in] Bitmask specifying enabled optimizations. +GrB_Info CFL_matrix_to_base( + // output + CFL_Matrix **matrix_p, + // input + CFL_Matrix *matrix, int8_t optimizations); + +// Moves the underlying base matrix out of a matrix (lazy or regular) without +// copying it. A lazy matrix is first merged into its largest base matrix. +// The source matrix doesn't own the extracted matrix anymore, but it still must be +// freed with CFL_matrix_free +// +// Parameters: +// base - [out] Extracted matrix. Caller must free +// matrix - [in] Pointer to the source matrix (may be lazy). +// optimizations - [in] Bitmask specifying enabled optimizations. +GrB_Info CFL_matrix_extract_base( + // output + GrB_Matrix *base, + // input + CFL_Matrix *matrix, int8_t optimizations); + +GrB_Info CFL_clear(CFL_Matrix *A, int8_t optimizations); \ No newline at end of file diff --git a/experimental/algorithm/LAGraph_CFL_reachability_advanced.c b/experimental/algorithm/LAGraph_CFL_reachability_advanced.c new file mode 100644 index 0000000000..44ddf0b394 --- /dev/null +++ b/experimental/algorithm/LAGraph_CFL_reachability_advanced.c @@ -0,0 +1,975 @@ +//------------------------------------------------------------------------------ +// LAGraph_CFL_reachability.c: 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. + +//------------------------------------------------------------------------------ + +// Code is based on the "A matrix-based CFPQ algorithm" described in the +// following paper: * Rustam Azimov, Semyon Grigorev, "Context-Free Path +// Querying Using Linear Algebra", URL: +// https://disser.spbu.ru/files/2022/disser_azimov.pdf + +#define LG_FREE_WORK \ + { \ + TRY_INNER(CFL_matrix_free(&iden)); \ + TRY_INNER(GrB_free(&identity_matrix)); \ + TRY_INNER(GrB_free(&v_diag)); \ + TRY_INNER(GrB_free(&new_adj_matrix)); \ + TRY_INNER(LAGraph_Free((void **)&to_new_symbols_map, NULL)); \ + TRY_INNER(LAGraph_Free((void **)&new_rules, NULL)); \ + TRY_INNER(LAGraph_Free((void **)&used_symbols, NULL)); \ + TRY_INNER(LAGraph_Free((void **)&nonterm_symbols, NULL)); \ + TRY_INNER(LAGraph_Free((void **)&active_symbols, NULL)); \ + for (size_t i = 0; i < new_symbols_amount; i++) { \ + if (temp_matrices != NULL) \ + TRY_INNER(CFL_matrix_free(&temp_matrices[i])); \ + if (delta_matrices != NULL) \ + TRY_INNER(CFL_matrix_free(&delta_matrices[i])); \ + if (matrices != NULL) \ + TRY_INNER(CFL_matrix_free(&matrices[i])); \ + if (new_adj_matrices != NULL && new_adj_matrices != adj_matrices) { \ + TRY_INNER(GrB_free(&new_adj_matrices[i])); \ + } \ + } \ + if (new_adj_matrices != adj_matrices) { \ + TRY_INNER(LAGraph_Free((void **)&new_adj_matrices, NULL)); \ + } \ + TRY_INNER(LAGraph_Free((void **)&delta_matrices, NULL)); \ + TRY_INNER(LAGraph_Free((void **)&matrices, NULL)); \ + TRY_INNER(LAGraph_Free((void **)&temp_matrices, NULL)); \ + } + +#define LG_FREE_ALL \ + { \ + LG_FREE_WORK; \ + if (outputs_initialized) { \ + for (size_t k = 0; k < symbols_amount; k++) { \ + TRY_INNER(GrB_free(&outputs[k])); \ + } \ + } \ + } + +#include "LAGraph_CFL_optimized_matrix_opt.h" +#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 != NULL && msg_len == 0) { \ + msg_len += \ + snprintf(msg, LAGRAPH_MSG_LEN, \ + "LAGraph failure (file %s, line %d): ", __FILE__, __LINE__); \ + } \ + if (msg != NULL && msg_len < LAGRAPH_MSG_LEN) { \ + msg_len += snprintf(msg + msg_len, LAGRAPH_MSG_LEN - msg_len, __VA_ARGS__); \ + } \ + } + +#define BENCH_CFL_REACHBILITY false + +#if BENCH_CFL_REACHBILITY + #define IS_ISO(matrix, str) \ + { \ + bool iso_flag; \ + GrB_Index nnz; \ + GxB_Matrix_iso(&iso_flag, matrix); \ + GrB_Matrix_nvals(&nnz, matrix); \ + if (!iso_flag && nnz) { \ + printf("-----ISO ALERT----- (%s)\n", str); \ + GxB_print(matrix, 1); \ + printf("-------------------\n"); \ + } \ + } + + #define TIMER_START() \ + { \ + start_time = LAGraph_WallClockTime(); \ + } + + #define TIMER_STOP(label, accumulator) \ + { \ + end_time = LAGraph_WallClockTime(); \ + printf("%s %.3fs\n", label, end_time - start_time); \ + if (accumulator != NULL) { \ + *(accumulator) += (end_time - start_time); \ + } \ + } + + #define IS_ROW(matrix, str) \ + { \ + int32_t orientation; \ + GrB_get(matrix, &orientation, GrB_STORAGE_ORIENTATION_HINT); \ + if (orientation != GrB_ROWMAJOR) { \ + printf("-----NOT A ROW----- (%s)\n", str); \ + GxB_print(matrix, 1); \ + printf("-------------------\n"); \ + } \ + } + + #define IS_COL(matrix, str) \ + { \ + int32_t orientation; \ + GrB_get(matrix, &orientation, GrB_STORAGE_ORIENTATION_HINT); \ + if (orientation != GrB_COLMAJOR) { \ + printf("-----NOT A COL----- (%s)\n", str); \ + GxB_print(matrix, 1); \ + printf("-------------------\n"); \ + } \ + } +#else + #define IS_ISO(matrix, str) + #define TIMER_START() + #define TIMER_STOP(label, accumulator) + #define IS_ROW(matrix, str) + #define IS_COL(matrix, str) +#endif +// clang-format on + +#define OPT_EMPTY (1 << 0) +#define OPT_FORMAT (1 << 1) +#define OPT_LAZY (1 << 2) +#define OPT_BLOCK (1 << 3) + +#define TRY(GrB_method) \ + { \ + GrB_Info LG_GrB_Info = GrB_method; \ + if (LG_GrB_Info < GrB_SUCCESS) { \ + fprintf(stderr, "LAGraph failure (file %s, line %d): \n", __FILE__, \ + __LINE__); \ + LG_FREE_ALL; \ + return (LG_GrB_Info); \ + } \ + } + +// Checks the return value of a GraphBLAS/LAGraph call inside a helper function. +// On failure, logs the error location to stderr, invokes FREE_INNER() to release +// any resources allocated within the current function, and returns the error code. +// +// FREE_INNER() must be defined by the caller before using this macro: +// +// #define FREE_INNER() \ +// { \ +// LAGraph_Free(&a) \ +// GrB_free(&M) \ +// } +// +// Note: use TRY() instead when the function is not a helper (i.e., it has +// its own top-level FREE_ALL cleanup). +#define TRY_INNER(GrB_method) \ + { \ + GrB_Info LG_GrB_Info = GrB_method; \ + if (LG_GrB_Info < GrB_SUCCESS) { \ + fprintf(stderr, "LAGraph failure (file %s, line %d): \n", __FILE__, \ + __LINE__); \ + FREE_INNER(); \ + return (LG_GrB_Info); \ + } \ + } + +#define TRY_I(GrB_method) \ + { \ + GrB_Info LG_GrB_Info = GrB_method; \ + if (LG_GrB_Info < GrB_SUCCESS) { \ + fprintf(stderr, \ + "LAGraph failure (file %s, line %d) (Iteration: %zu, i: %zu): \n", \ + __FILE__, __LINE__, iteration, i); \ + LG_FREE_ALL; \ + return (LG_GrB_Info); \ + } \ + } + +typedef struct { + int32_t index; + int32_t base_index; + int32_t count; +} CFL_Symbol; + +// When using the OPT_BLOCK optimization, indexed symbols must be grouped together. +// This produces a mapping: [old_index -> (new_index, base_index, indexed_count)] +// - new_index: the index of the symbol in the new numeration +// - base_index: the original index of the first symbol in the indexed group +// - indexed_count: the number of indexed symbols in the group (0 if not indexed) +// +// Example: +// (0) S -> (0, 0, 0) - non-indexed symbol +// (1) A_0 -> (1, 1, 3) - indexed group of 3, starting at old index 1 +// (2) A_1 -> } +// (3) A_2 -> } members of the A group +// (4) B_0 -> (2, 4, 2) - indexed group of 2, starting at old index 4 +// (5) B_1 -> } member of the B group +// (6) C -> (3, 6, 0) - non-indexed symbol +// (7) a -> (4, 7, 0) - non-indexed symbol +// +// This mapping is used to build a compact matrix array and to expand +// production rules, both with and without the OPT_BLOCK optimization. +// +// Output: CFL_Symbol **symbols and size_t *size +static GrB_Info get_new_symbols(const LAGraph_rule_EWCNF *rules, size_t rules_count, + size_t symbols_amount, CFL_Symbol **symbols, size_t *size, + char *msg) { + *symbols = NULL; + bool *checked = NULL; + +#undef FREE_INNER_WORK +#undef FREE_INNER + +#define FREE_INNER_WORK() \ + { \ + LAGraph_Free((void **)&checked, NULL); \ + }; + +#define FREE_INNER() \ + { \ + FREE_INNER_WORK(); \ + LAGraph_Free((void **)symbols, NULL); \ + } + + TRY_INNER(LAGraph_Calloc((void **)&checked, symbols_amount, sizeof(bool), msg)); + for (size_t i = 0; i < symbols_amount; i++) { + checked[i] = false; + } + + size_t capacity = 1; + *size = 0; + TRY_INNER(LAGraph_Calloc((void **)symbols, capacity, sizeof(CFL_Symbol), msg)); + + for (size_t i = 0; i < capacity; i++) { + CFL_Symbol sym = {0}; + (*symbols)[i] = sym; + } + + for (size_t i = 0; i < rules_count; i++) { + LAGraph_rule_EWCNF rule = rules[i]; + + int32_t prods[3] = {rule.nonterm, rule.prod_A, rule.prod_B}; + int bitmasks[3] = {LAGraph_EWNCF_INDEX_NONTERM, LAGraph_EWNCF_INDEX_PROD_A, + LAGraph_EWNCF_INDEX_PROD_B}; + + for (size_t j = 0; j < 3; j++) { + if (prods[j] == -1) { + continue; + } + + if (checked[prods[j]]) + continue; + + checked[prods[j]] = true; + + CFL_Symbol sym; + sym.base_index = prods[j]; + sym.index = *size; + sym.count = rule.indexed & bitmasks[j] ? rule.indexed_count : 0; + + for (size_t k = sym.base_index; k < sym.base_index + sym.count; k++) { + checked[k] = true; + } + + if (*size == capacity) { + capacity *= 2; + TRY_INNER(LAGraph_Realloc((void **)symbols, capacity, capacity / 2, + sizeof(CFL_Symbol), msg)); + } + + (*symbols)[(*size)++] = sym; + } + } + + for (size_t i = 0; i < symbols_amount; i++) { + if (checked[i]) + continue; + + if (*size == capacity) { + capacity *= 2; + TRY_INNER(LAGraph_Realloc((void **)symbols, capacity, capacity / 2, + sizeof(CFL_Symbol), msg)); + } + + CFL_Symbol sym; + sym.base_index = i; + sym.index = *size; + sym.count = 0; + + (*symbols)[(*size)++] = sym; + checked[i] = true; + // printf("Inserted (%ld, %ld, %ld)\n", sym.index, sym.base_index, sym.count); + } + + FREE_INNER_WORK(); + + return GrB_SUCCESS; +} + +// Expands indexed grammar rules into a set of concrete rules. +// +// Before: A_i -> B_i C (indexed_count = 3) +// D -> E F +// After: A_0 -> B_0 C +// A_1 -> B_1 C +// A_2 -> B_2 C +// D -> E F +// +// Output: LAGraph_rule_EWCNF **new_rules and size_t *new_rules_count +static GrB_Info explode_rules(const LAGraph_rule_EWCNF *rules, size_t rules_count, + LAGraph_rule_EWCNF **new_rules, size_t *new_rules_count, + char *msg) { + *new_rules = NULL; + +#undef FREE_INNER + +#define FREE_INNER() \ + { \ + LAGraph_Free((void **)new_rules, NULL); \ + } + + size_t new_rules_size = 0; + size_t new_rules_capacity = 1; + + TRY_INNER(LAGraph_Calloc((void **)new_rules, new_rules_capacity, + sizeof(LAGraph_rule_EWCNF), msg)); + + // explode rules + for (size_t i_rule = 0; i_rule < rules_count; i_rule++) { + LAGraph_rule_EWCNF rule = rules[i_rule]; + + if (new_rules_size == new_rules_capacity) { + TRY_INNER(LAGraph_Realloc((void **)new_rules, new_rules_capacity * 2, + new_rules_capacity, sizeof(LAGraph_rule_EWCNF), + msg)); + new_rules_capacity *= 2; + } + + if (rule.indexed_count == 0) { + (*new_rules)[new_rules_size++] = rule; + } else { + while (new_rules_size + rule.indexed_count >= new_rules_capacity) { + TRY_INNER(LAGraph_Realloc((void **)new_rules, new_rules_capacity * 2, + new_rules_capacity, sizeof(LAGraph_rule_EWCNF), + msg)); + new_rules_capacity *= 2; + } + + for (size_t rule_index = 0; rule_index < rule.indexed_count; rule_index++) { + LAGraph_rule_EWCNF new_rule = rule; + new_rule.indexed_count = 0; + new_rule.indexed = 0; + if (rule.nonterm != -1 && rule.indexed & LAGraph_EWNCF_INDEX_NONTERM) { + new_rule.nonterm = rule.nonterm + rule_index; + } + if (rule.prod_A != -1 && rule.indexed & LAGraph_EWNCF_INDEX_PROD_A) { + new_rule.prod_A = rule.prod_A + rule_index; + } + if (rule.prod_B != -1 && rule.indexed & LAGraph_EWNCF_INDEX_PROD_B) { + new_rule.prod_B = rule.prod_B + rule_index; + } + + (*new_rules)[new_rules_size++] = new_rule; + } + } + } + + *new_rules_count = new_rules_size; + + return GrB_SUCCESS; +} + +// Moves a CFL_Matrix into an array of GrB_Matrix matrices +// +// If the matrix is not a horizontal or vertical block vector (i.e., it is a square +// matrix), the underlying GrB_Matrix is moved into outputs[0] without copying +// +// Otherwise, the matrix split into square sub-matrices of size (graph_size x graph_size) +// Parameters: +// outputs - [out] Caller-allocated array of GrB_Matrix to write results into +// Must have enough space for all sub-matrices +// matrix - [in] Source CFL_Matrix to move. It doesn't own its base +// matrix after the call, but it still must be freed +static GrB_Info split_CFL_matrix(GrB_Matrix *outputs, CFL_Matrix *matrix, + int8_t optimizations) { + char msg[LAGRAPH_MSG_LEN]; + GrB_Matrix base = NULL; + GrB_Index *nrows = NULL, *ncols = NULL; + +#undef FREE_INNER + +#define FREE_INNER() \ + { \ + GrB_free(&base); \ + LAGraph_Free((void **)&nrows, NULL); \ + LAGraph_Free((void **)&ncols, NULL); \ + } + + TRY_INNER(CFL_matrix_extract_base(&base, matrix, optimizations)); + + GrB_Index base_nrows, base_ncols; + TRY_INNER(GrB_Matrix_nrows(&base_nrows, base)); + TRY_INNER(GrB_Matrix_ncols(&base_ncols, base)); + + if (base_nrows == base_ncols) { + *outputs = base; + return GrB_SUCCESS; + } + + // we can create nrows and ncols array with the same size, it will be mush easier than + // calculate size of each array :) + GrB_Index matrices_count = + base_nrows > base_ncols ? base_nrows / base_ncols : base_ncols / base_nrows; + + TRY_INNER(LAGraph_Calloc((void **)&nrows, matrices_count, sizeof(GrB_Index), msg)); + TRY_INNER(LAGraph_Calloc((void **)&ncols, matrices_count, sizeof(GrB_Index), msg)); + + GrB_Index graph_size = base_nrows < base_ncols ? base_nrows : base_ncols; + for (size_t i = 0; i < matrices_count; i++) { + nrows[i] = graph_size; + ncols[i] = graph_size; + } + + GrB_Index m = 0, n = 0; + if (base_nrows > base_ncols) { + m = matrices_count; + n = 1; + } else { + m = 1; + n = matrices_count; + } + + TRY_INNER(GxB_Matrix_split(outputs, m, n, nrows, ncols, base, GrB_NULL)); + + TRY_INNER(LAGraph_Free((void **)&nrows, NULL)); + TRY_INNER(LAGraph_Free((void **)&ncols, NULL)); + TRY_INNER(GrB_free(&base)); + + return GrB_SUCCESS; +} + +// Builds the mapping [old_index -> (new_index, base_index, indexed_count)] +// See get_new_symbols() for a detailed description of the mapping format +// +// Parameters: +// map - [out] Allocated array of CFL_Symbol. Caller must free +// size - [out] Number of entries in map +static GrB_Info get_new_symbols_map(const LAGraph_rule_EWCNF *rules, size_t rules_count, + size_t symbols_amount, CFL_Symbol **map, size_t *size, + char *msg, int8_t optimizations) { +#undef FREE_INNER + +#define FREE_INNER() \ + { \ + LAGraph_Free((void **)map, NULL); \ + } + + if (optimizations & OPT_BLOCK) { + TRY_INNER(get_new_symbols(rules, rules_count, symbols_amount, map, size, msg)); + } else { + *size = symbols_amount; + TRY_INNER(LAGraph_Calloc((void **)map, symbols_amount, sizeof(CFL_Symbol), msg)); + for (size_t i = 0; i < symbols_amount; i++) { + CFL_Symbol *sym = &(*map)[i]; + sym->index = i; + sym->base_index = i; + sym->count = 0; + } + } + + return GrB_SUCCESS; +} + +// Builds a new array of adjacency matrices according to the symbol mapping +// +// Parameters: +// new_adj_matrices_p - [out] Resulting matrix array. Caller must free +// (only if OPT_BLOCK is set). +static GrB_Info get_new_adj_matrices(const GrB_Matrix *adj_matrices, CFL_Symbol *map, + GrB_Index map_size, GrB_Matrix **new_adj_matrices_p, + char *msg, int8_t optimizations) { +#undef FREE_INNER + +#define FREE_INNER() \ + { \ + LAGraph_Free((void **)new_adj_matrices_p, NULL); \ + } + + GrB_Index n; + TRY_INNER(GrB_Matrix_ncols(&n, adj_matrices[0])); + + if (!(optimizations & OPT_BLOCK)) { + *new_adj_matrices_p = (GrB_Matrix *)adj_matrices; + return GrB_SUCCESS; + } + + TRY_INNER( + LAGraph_Calloc((void **)new_adj_matrices_p, map_size, sizeof(GrB_Matrix), msg)); + + for (size_t i = 0; i < map_size; i++) { + CFL_Symbol sym = map[i]; + + if (sym.count == 0) { + TRY_INNER( + GrB_Matrix_dup(&(*new_adj_matrices_p)[i], adj_matrices[sym.base_index])); + continue; + } + + GrB_Matrix new_col_matrix; + TRY_INNER(GrB_Matrix_new(&new_col_matrix, GrB_BOOL, n * sym.count, n)); + GrB_Matrix *Tiles = (GrB_Matrix *)adj_matrices + sym.base_index; + TRY_INNER(GxB_Matrix_concat(new_col_matrix, Tiles, sym.count, 1, GrB_NULL)); + (*new_adj_matrices_p)[i] = new_col_matrix; + } + + return GrB_SUCCESS; +} + +// Remaps rule symbol indices according to the symbol mapping +// +// Parameters: +// new_rules - [out] Allocated output rule array. Caller must free +// new_rules_count - [out] Number of rules written to new_rules +static GrB_Info get_new_rules(const LAGraph_rule_EWCNF *rules, size_t rules_count, + CFL_Symbol *map, size_t map_size, + LAGraph_rule_EWCNF **new_rules, size_t *new_rules_count, + char *msg, int8_t optimizations) { + *new_rules_count = 0; +#undef FREE_INNER + +#define FREE_INNER() \ + { \ + LAGraph_Free((void **)new_rules, NULL); \ + } + + if (!(optimizations & OPT_BLOCK)) { + TRY_INNER(explode_rules(rules, rules_count, new_rules, new_rules_count, msg)); + return GrB_SUCCESS; + } + + TRY_INNER( + LAGraph_Calloc((void **)new_rules, rules_count, sizeof(LAGraph_rule_EWCNF), msg)); + for (size_t i = 0; i < rules_count; i++) { + LAGraph_rule_EWCNF rule = rules[i]; + LAGraph_rule_EWCNF new_rule = rule; + + for (size_t i_sym = 0; i_sym < map_size; i_sym++) { + CFL_Symbol sym = map[i_sym]; + + if (rule.nonterm != -1 && rule.nonterm == sym.base_index) { + new_rule.nonterm = sym.index; + } + if (rule.prod_A != -1 && rule.prod_A == sym.base_index) { + new_rule.prod_A = sym.index; + } + if (rule.prod_B != -1 && rule.prod_B == sym.base_index) { + new_rule.prod_B = sym.index; + } + } + + (*new_rules)[(*new_rules_count)++] = new_rule; + } + + return GrB_SUCCESS; +} + +// LAGraph_CFL_reachability_adv: 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 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`). +// +// Terminals and non-terminals are enumerated by integers starting from zero. +// The start non-terminal is the non-terminal with index 0. +// +// Example: +// +// Graph: +// ┌───┐ ┌───┐ ┌───┐ ┌───┐ ┌───┐ +// │ 0 ├───► 1 ├───► 2 ├───► 3 ├───► 4 │ +// └───┘ a └─┬─┘ a └─▲─┘ b └───┘ b └───┘ +// │ │ +// │ ┌───┐ │ +// a└─► 5 ├─┘b +// └───┘ +// +// Grammar: S -> aSb | ab +// +// There are paths from node [1] to node [3] and from node [1] to node [2] that form +// the word "ab" ([1]-a->[2]-b->[3] and [1]-a->[5]-b->[2]). The word "ab" is in the +// language generated by our context-free grammar, so the pairs (1, 3) and (1, 2) will +// be included in the result. +// +// Note: It doesn't matter how many paths exist from node [A] to node [B] that form a +// word in the language. If at least one path exists, the pair ([A], [B]) will be +// included in the result. +// +// In contrast, the path from node [1] to node [4] forms the word "abb" +// ([1]-a->[2]-b->[3]-b->[4]) and the word "abbb" ([1]-a->[5]-b->[2]-b->[3]-b->[4]). +// The words "aab" and "abbb" are not in the language, so the pair (1, 4) will not be +// included in the result. +// +// With this graph and grammar, we obtain the following results: +// (0, 4) - because there exists a path (0-1-2-3-4) that forms the word "aabb" +// (1, 3) - because there exists a path (1-2-3) that forms "ab" +// (1, 2) - because there exists a path (1-5-2) that forms the word "ab" +// (0, 3) - because there exists a path (0-1-5-2-3) that forms the word "aabb" +GrB_Info LAGraph_CFL_reachability_adv( + // Output + GrB_Matrix *outputs, // Array of matrices containing results. + // The size of the array must be equal to the count + // of symbols (symbols_amount). + // + // Each matrix is square, with size equal to the number of + // vertices in the graph. Matrices are allocated by this + // method, the caller must free them. + // + // outputs[k]: (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 nonterminal 'k' of the specified CFG. + // + // Note: outputs[t] is NULL if the symbol 't' is not a + // nonterminal, i.e. it is not on the left side of any rule + // (terminals and symbols unused by rules). + // Input + const GrB_Matrix *adj_matrices, // Array of adjacency matrices representing the graph. + // The length of this array is equal to the count of + // symbols (symbols_amount). + // + // Each matrix is square, with size equal to the + // number of vertices in the graph. + // + // adj_matrices[i]: (i, j) == 1 if and only if there + // is an edge between nodes i and j with the label of + // the symbol corresponding to index 'i' (where i is + // in the range [0, symbols_amount - 1]). + // + // Note: adj_matrices[N], where N is the index of a + // nonterminal, must be initialized as an empty square + // matrix of size symbols_amount. + size_t symbols_amount, // Count of terminal and nonterminals + const LAGraph_rule_EWCNF *rules, // The rules of the CFG. + // Warning: now not ready for N -> A rules, where + // is N and A are nonterminals. + size_t rules_count, // The total number of rules in the CFG. + char *msg, // Message string for error reporting. + int8_t optimizations // Optimizations flags +) { +#undef FREE_INNER + +#define FREE_INNER() + + // Declare workspace and clear the msg string, if not NULL + CFL_Matrix **delta_matrices = NULL, **matrices = NULL, **temp_matrices = NULL; + CFL_Matrix *iden = NULL; + GrB_Matrix identity_matrix = NULL; + GrB_Vector v_diag = NULL; + GrB_Matrix new_adj_matrix = NULL; + + // for OPT_BLOCK optimization + size_t new_symbols_amount = 0; + GrB_Matrix *new_adj_matrices = NULL; + LAGraph_rule_EWCNF *new_rules = NULL; + size_t new_rules_count = 0; + CFL_Symbol *to_new_symbols_map = NULL; + + // Symbols that appear in no rule never change, so the main loop skips them + bool *used_symbols = NULL; + // Only nonterminals (symbols on the left side of a rule) are given to the caller + bool *nonterm_symbols = NULL; + bool outputs_initialized = false; + size_t *active_symbols = NULL; + size_t active_symbols_count = 0; + + LG_CLEAR_MSG; + size_t msg_len = 0; // For error formatting + + TRY(LAGraph_Calloc((void **)&delta_matrices, symbols_amount, sizeof(CFL_Matrix *), + msg)); + TRY(LAGraph_Calloc((void **)&matrices, symbols_amount, sizeof(CFL_Matrix *), msg)); + TRY(LAGraph_Calloc((void **)&temp_matrices, symbols_amount, sizeof(CFL_Matrix *), + msg)); + + LG_ASSERT_MSG(symbols_amount > 0, GrB_INVALID_VALUE, + "The number of symbols must be greater than zero."); + LG_ASSERT_MSG(rules_count > 0, GrB_INVALID_VALUE, + "The number of rules must be greater than zero."); + LG_ASSERT_MSG(outputs != NULL, GrB_NULL_POINTER, "The outputs array cannot be null."); + LG_ASSERT_MSG(rules != NULL, GrB_NULL_POINTER, "The rules array cannot be null."); + LG_ASSERT_MSG(adj_matrices != NULL, GrB_NULL_POINTER, + "The adjacency matrices array cannot be null."); + + for (size_t i = 0; i < symbols_amount; i++) { + outputs[i] = NULL; + } + outputs_initialized = true; + + // Find null adjacency matrices + bool found_null = false; + for (size_t i = 0; i < symbols_amount; i++) { + if (adj_matrices[i] != NULL) + continue; + + if (!found_null) { + ADD_TO_MSG("Adjacency matrices with these indexes are null: "); + ADD_TO_MSG("%ld", i); + } else { + ADD_TO_MSG(", %ld", i); + } + + found_null = true; + } + + if (found_null) { + LG_FREE_ALL; + return GrB_NULL_POINTER; + } + + // Validate original rules before grouping symbols or expanding indexed rules. + for (size_t i = 0; i < rules_count; i++) { + LAGraph_rule_EWCNF rule = rules[i]; + LG_ASSERT_MSGF(rule.prod_A != -1 || rule.prod_B == -1, GrB_INVALID_VALUE, + "Rule with index %zu has an invalid [Variable -> _ B] form.", i); + + int32_t prods[3] = {rule.nonterm, rule.prod_A, rule.prod_B}; + int bitmasks[3] = {LAGraph_EWNCF_INDEX_NONTERM, LAGraph_EWNCF_INDEX_PROD_A, + LAGraph_EWNCF_INDEX_PROD_B}; + + for (size_t j = 0; j < 3; j++) { + if (j != 0 && prods[j] == -1) { + continue; + } + + LG_ASSERT_MSGF(prods[j] >= 0 && (size_t)prods[j] < symbols_amount, + GrB_INVALID_VALUE, + "Rule with index %zu has an invalid symbol index %d.", + i, prods[j]); + + if (rule.indexed & bitmasks[j]) { + LG_ASSERT_MSGF(rule.indexed_count <= symbols_amount - (size_t)prods[j], + GrB_INVALID_VALUE, + "Rule with index %zu has an indexed group outside the " + "symbol range.", i); + } + } + } + + GrB_Index n; + TRY(GrB_Matrix_ncols(&n, adj_matrices[0])); + + TRY(get_new_symbols_map(rules, rules_count, symbols_amount, &to_new_symbols_map, + &new_symbols_amount, msg, optimizations)); + TRY(get_new_adj_matrices(adj_matrices, to_new_symbols_map, new_symbols_amount, + &new_adj_matrices, msg, optimizations)); + TRY(get_new_rules(rules, rules_count, to_new_symbols_map, new_symbols_amount, + &new_rules, &new_rules_count, msg, optimizations)); + + // Arrays for processing rules + size_t eps_rules[new_rules_count], eps_rules_count = 0; // [Variable -> eps] + size_t term_rules[new_rules_count], term_rules_count = 0; // [Variable -> term] + size_t bin_rules[new_rules_count], bin_rules_count = 0; // [Variable -> AB] + + // Process rules + for (size_t i = 0; i < new_rules_count; i++) { + LAGraph_rule_EWCNF rule = new_rules[i]; + + if (rule.prod_A == -1) { + // [Variable -> eps] + eps_rules[eps_rules_count++] = i; + } else if (rule.prod_B == -1) { + // [Variable -> term] + term_rules[term_rules_count++] = i; + } else { + // [Variable -> A B] + bin_rules[bin_rules_count++] = i; + } + } + + // Collect symbols used by rules. An unused symbol is never read or written + // by the algorithm, so its output is just a copy of its adjacency matrix. + TRY(LAGraph_Calloc((void **)&used_symbols, new_symbols_amount, sizeof(bool), msg)); + TRY(LAGraph_Calloc((void **)&active_symbols, new_symbols_amount, sizeof(size_t), + msg)); + TRY(LAGraph_Calloc((void **)&nonterm_symbols, new_symbols_amount, sizeof(bool), + msg)); + + for (size_t i = 0; i < new_rules_count; i++) { + LAGraph_rule_EWCNF rule = new_rules[i]; + used_symbols[rule.nonterm] = true; + nonterm_symbols[rule.nonterm] = true; + if (rule.prod_A != -1) + used_symbols[rule.prod_A] = true; + if (rule.prod_B != -1) + used_symbols[rule.prod_B] = true; + } + + for (size_t i = 0; i < new_symbols_amount; i++) { + if (used_symbols[i]) + active_symbols[active_symbols_count++] = i; + } + + // Create symbol matrices + for (size_t k = 0; k < active_symbols_count; k++) { + size_t i = active_symbols[k]; + GrB_Index nrows; + TRY(GrB_Matrix_nrows(&nrows, new_adj_matrices[i])); + GrB_Index ncols; + TRY(GrB_Matrix_ncols(&ncols, new_adj_matrices[i])); + + TRY(GrB_Matrix_dup(&new_adj_matrix, new_adj_matrices[i])); + TRY(CFL_matrix_from_base(&delta_matrices[i], new_adj_matrix)); + new_adj_matrix = NULL; // Ownership transferred to delta_matrices[i]. + + if (optimizations & OPT_LAZY) { + TRY(CFL_matrix_create_lazy(&matrices[i], nrows, ncols)); + } else { + TRY(CFL_matrix_create(&matrices[i], nrows, ncols)); + } + + TRY(CFL_matrix_create(&temp_matrices[i], nrows, ncols)); + } + + // Rule [Variable -> term] + for (size_t i = 0; i < term_rules_count; i++) { + LAGraph_rule_EWCNF term_rule = new_rules[term_rules[i]]; + CFL_Matrix *nonterm_matrix = delta_matrices[term_rule.nonterm]; + CFL_Matrix *term_matrix = delta_matrices[term_rule.prod_A]; + + TRY(CFL_wise(nonterm_matrix, nonterm_matrix, term_matrix, true, optimizations)); + } + + TRY(GrB_Vector_new(&v_diag, GrB_BOOL, n)); + TRY(GrB_Vector_assign_BOOL(v_diag, GrB_NULL, GrB_NULL, true, GrB_ALL, n, NULL)); + TRY(GrB_Matrix_diag(&identity_matrix, v_diag, 0)); + TRY(GrB_Vector_free(&v_diag)); + TRY(CFL_matrix_from_base(&iden, identity_matrix)); + identity_matrix = NULL; // Ownership transferred to iden. + + // Rule [Variable -> eps] + for (size_t i = 0; i < eps_rules_count; i++) { + LAGraph_rule_EWCNF eps_rule = new_rules[eps_rules[i]]; + CFL_Matrix *nonterm_matrix = delta_matrices[eps_rule.nonterm]; + + TRY(CFL_wise(nonterm_matrix, nonterm_matrix, iden, true, optimizations)); + } + + // Rule [Variable -> Variable1 Variable2] + double start_time, end_time; + bool changed = true; + size_t iteration = 0; + double mxm1 = 0.0; + double wise1 = 0.0; + double mxm2 = 0.0; + double wise2 = 0.0; + double rsubt = 0.0; + while (changed) { + iteration++; + changed = false; + +#if BENCH_CFL_REACHBILITY + printf("\n--- ITERATARION %ld ---\n", iteration); +#endif + + for (size_t k = 0; k < active_symbols_count; k++) { + size_t i = active_symbols[k]; + TRY_I(CFL_matrix_free(&temp_matrices[i])); + TRY_I(CFL_matrix_create(&temp_matrices[i], matrices[i]->nrows, + matrices[i]->ncols)); + } + + TIMER_START(); + for (size_t i = 0; i < bin_rules_count; i++) { + LAGraph_rule_EWCNF bin_rule = new_rules[bin_rules[i]]; + CFL_Matrix *A = matrices[bin_rule.prod_A]; + CFL_Matrix *B = delta_matrices[bin_rule.prod_B]; + CFL_Matrix *C = temp_matrices[bin_rule.nonterm]; + + TRY_I(CFL_mxm(C, A, B, true, false, optimizations)); + } + TIMER_STOP("MXM 1", &mxm1); + + TIMER_START() + for (size_t k = 0; k < active_symbols_count; k++) { + size_t i = active_symbols[k]; + CFL_Matrix *A = delta_matrices[i]; + CFL_Matrix *C = matrices[i]; + + TRY_I(CFL_wise(C, C, A, false, optimizations)); + } + TIMER_STOP("WISE 1", &wise1); + + TIMER_START() + for (size_t i = 0; i < bin_rules_count; i++) { + LAGraph_rule_EWCNF bin_rule = new_rules[bin_rules[i]]; + CFL_Matrix *A = matrices[bin_rule.prod_B]; + CFL_Matrix *B = delta_matrices[bin_rule.prod_A]; + CFL_Matrix *C = temp_matrices[bin_rule.nonterm]; + + TRY(CFL_mxm(C, A, B, true, true, optimizations)); + } + TIMER_STOP("MXM 2", &mxm2); + + // Rule [Variable -> term] + for (size_t i = 0; i < term_rules_count; i++) { + LAGraph_rule_EWCNF term_rule = new_rules[term_rules[i]]; + CFL_Matrix *A = temp_matrices[term_rule.nonterm]; + CFL_Matrix *B = delta_matrices[term_rule.prod_A]; + + TRY_I(CFL_wise(A, A, B, true, optimizations)); + } + + TIMER_START(); + for (size_t k = 0; k < active_symbols_count; k++) { + size_t i = active_symbols[k]; + TRY_I(CFL_dup(delta_matrices[i], temp_matrices[i], optimizations)); + } + TIMER_STOP("WISE 2 (copy)", &wise2); + + TIMER_START(); + for (size_t k = 0; k < active_symbols_count; k++) { + size_t i = active_symbols[k]; + CFL_Matrix *A = matrices[i]; + CFL_Matrix *C = delta_matrices[i]; + + TRY_I(CFL_rsub(C, A, optimizations)); + } + TIMER_STOP("WISE 3 (MASK)", &rsubt); + + size_t new_nnz = 0; + for (size_t k = 0; k < active_symbols_count; k++) { + size_t i = active_symbols[k]; + TRY(CFL_matrix_update(delta_matrices[i])); + new_nnz += delta_matrices[i]->nvals; + } + + if (new_nnz != 0) { + changed = true; + } + } + +#if BENCH_CFL_REACHBILITY + printf("MXM1: %.3f, wise1: %.3f, MXM2: %.3f, wise2: %.3f, rsub: %.3f", mxm1, wise1, + mxm2, wise2, rsubt); +#endif + + // get outputs matrices + // Nonterminal matrices are moved to the caller, outputs of other symbols stay NULL + for (size_t i = 0; i < new_symbols_amount; i++) { + if (!nonterm_symbols[i]) { + continue; + } + + CFL_Symbol sym = to_new_symbols_map[i]; + TRY(split_CFL_matrix(outputs + sym.base_index, matrices[i], optimizations)); + } + + LG_FREE_WORK; + 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..b9eb7da66d --- /dev/null +++ b/experimental/test/test_CFL_optimized_matrix_opt.c @@ -0,0 +1,1676 @@ +//------------------------------------------------------------------------------ +// 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 "../experimental/algorithm/LAGraph_CFL_optimized_matrix_opt.h" +#include +#include +#include +#include +#include +#include + +#define NOT_OK(method) TEST_CHECK(method != 0) + +#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_p, + enum CFL_Matrix_block format); +extern GrB_Info block_matrix_to_diag(Matrix *diag, Matrix *input); +extern GrB_Info block_matrix_reduce(Matrix *matrix, Matrix *input, int8_t optimizations); +extern GrB_Info matrix_combine_lazy(Matrix *matrix, size_t threshold, + 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)); + TEST_MSG(msg); +} + +static void teardown(void) { + OK(LAGraph_Finalize(msg)); + TEST_MSG(msg); +} + +// 0 1 .. 0 +// 1 0 .. 0 +// .. .. .. .. +// 0 0 .. 0 +static GrB_Info make_simple_matrix(Matrix **M, 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)); + + CFL_matrix_from_base(M, A); + return GrB_SUCCESS; +} + +// 1 0 .. 0 +// 0 1 .. 0 +// .. .. .. .. +// 0 0 .. 0 +static GrB_Info make_simple_matrix_inverted(Matrix **M, 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)); + + CFL_matrix_from_base(M, A); + return GrB_SUCCESS; +} + +// 1 1 .. 1 +// 1 1 .. 1 +// .. .. .. .. +// 1 1 .. 1 +static GrB_Info make_ones_matrix(Matrix **M, 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)); + } + } + + OK(CFL_matrix_from_base(M, A)); + return GrB_SUCCESS; +} + +// TODO: make free variadic function +static GrB_Info free_matrix(Matrix **M) { + OK(CFL_matrix_free(M)); + return GrB_SUCCESS; +} + +static bool compare_matrices(Matrix *first, Matrix *second) { + Matrix *A, *B, *C; + OK(CFL_matrix_to_base(&A, first, OPT_LAZY)); + OK(CFL_matrix_to_base(&B, second, OPT_LAZY)); + OK(CFL_matrix_create(&C, A->nrows, A->ncols)); + + if (A->nvals != B->nvals) { + OK(CFL_matrix_free(&A)); + OK(CFL_matrix_free(&B)); + OK(CFL_matrix_free(&C)); + + return false; + } + + OK(GrB_eWiseMult(C->base, GrB_NULL, false, GrB_EQ_BOOL, A->base, B->base, GrB_NULL)); + OK(CFL_matrix_update(C)); + + bool result = C->nvals == A->nvals; + + OK(CFL_matrix_free(&A)); + OK(CFL_matrix_free(&B)); + OK(CFL_matrix_free(&C)); + + return result; +} + +static void test_compare_matrices_function(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + { + Matrix *A, *B; + make_simple_matrix(&A, 5); + CFL_matrix_create(&B, 5, 5); + + TEST_CHECK(!compare_matrices(A, B)); + + CFL_matrix_free(&A); + CFL_matrix_free(&B); + } + + { + Matrix *A, *B; + make_simple_matrix(&A, 5); + make_simple_matrix_inverted(&B, 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(); + + { + Matrix *A; + OK(CFL_matrix_create(&A, 4, 8)); + TEST_CHECK(A->nvals == 0); + TEST_CHECK(A->nrows == 4); + TEST_CHECK(A->ncols == 8); + TEST_CHECK(A->base != NULL); + + OK(CFL_matrix_free(&A)); + TEST_CHECK(A == NULL); + } + + { + GrB_Matrix A_base = NULL; + OK(GrB_Matrix_new(&A_base, GrB_BOOL, 4, 8)); + OK(GrB_Matrix_setElement_BOOL(A_base, true, 1, 1)); + + Matrix *A; + OK(CFL_matrix_from_base(&A, A_base)); + TEST_CHECK(A->nvals == 1); + TEST_CHECK(A->nrows == 4); + TEST_CHECK(A->ncols == 8); + TEST_CHECK(A->base != NULL); + TEST_CHECK(A->base == A_base); + + OK(CFL_matrix_free(&A)); + TEST_CHECK(A == NULL); + } + + { + Matrix *A; + OK(CFL_matrix_create_lazy(&A, 4, 8)); + TEST_CHECK(A->nvals == 0); + TEST_CHECK(A->nrows == 4); + TEST_CHECK(A->ncols == 8); + TEST_CHECK(A->base == NULL); + TEST_CHECK(A->is_lazy == true); + TEST_CHECK(A->base_matrices_count != 0); + TEST_CHECK(A->base_matrices[0] != NULL); + TEST_CHECK(A->base_matrices[0]->nvals == 0); + TEST_CHECK(A->base_matrices[0]->nrows == 4); + TEST_CHECK(A->base_matrices[0]->ncols == 8); + + OK(CFL_matrix_free(&A)); + TEST_CHECK(A == NULL); + } + + { + GrB_Matrix A_base = NULL; + OK(GrB_Matrix_new(&A_base, GrB_BOOL, 4, 8)); + OK(GrB_Matrix_setElement_BOOL(A_base, true, 1, 1)); + + Matrix *A; + OK(CFL_matrix_from_base_lazy(&A, A_base)); + TEST_CHECK(A->nvals == 1); + TEST_CHECK(A->nrows == 4); + TEST_CHECK(A->ncols == 8); + TEST_CHECK(A->base == NULL); + TEST_CHECK(A->is_lazy == true); + TEST_CHECK(A->base_matrices_count != 0); + TEST_CHECK(A->base_matrices[0] != NULL); + TEST_CHECK(A->base_matrices[0]->nvals == 1); + TEST_CHECK(A->base_matrices[0]->nrows == 4); + TEST_CHECK(A->base_matrices[0]->ncols == 8); + + OK(CFL_matrix_free(&A)); + TEST_CHECK(A == 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, *B, *C; + make_simple_matrix(&A, 2); + make_simple_matrix(&B, 2); + matrix_to_format(A, GrB_ROWMAJOR, true); + CFL_matrix_create(&C, 2, 2); + + OK(CFL_mxm(C, A, B, false, false, mask)); + TEST_MSG("matrix_mxm_opt failed for mask=%d", mask); + + 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, *B; + OK(make_simple_matrix(&A, 2)); + OK(make_simple_matrix(&B, 2)); + OK(matrix_to_format(A, GrB_ROWMAJOR, true)); + + OK(CFL_wise(A, A, B, false, mask)); + TEST_MSG("matrix_wise failed for mask=%d", mask); + + 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); + + OK(free_matrix(&A)); + OK(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, *B; + OK(make_simple_matrix(&A, 2)); + OK(make_simple_matrix(&B, 2)); + OK(matrix_to_format(A, GrB_ROWMAJOR, true)); + + OK(CFL_rsub(A, B, mask)); + TEST_MSG("matrix_rsub failed for mask=%d", mask); + + 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); + + OK(free_matrix(&A)); + OK(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(&A, 3); + + OK(CFL_clear(A, mask)); + GrB_Index nvals; + TEST_CHECK(A->nvals == 0); + + OK(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++) { + Matrix *A, *B; + OK(make_simple_matrix(&A, 3)); + OK(CFL_matrix_create(&B, 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); + OK(free_matrix(&A)); + OK(free_matrix(&B)); + } + + teardown(); +#endif +} + +static void test_CFL_dup_same(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix *A; + OK(make_simple_matrix(&A, 3)); + OK(CFL_dup(A, A, 0)); + + OK(free_matrix(&A)); + + teardown(); +#endif +} + +//------------------------------------------------------------------------------ +// Format optimization +//------------------------------------------------------------------------------ + +static void test_CFL_format_create_rowmajor_matrix(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix *A; + OK(make_simple_matrix(&A, 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); + + OK(free_matrix(&A)); + + teardown(); +#endif +} + +static void test_CFL_format_to_format(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix *A; + OK(make_simple_matrix(&A, 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); + + OK(free_matrix(&A)); + + teardown(); +#endif +} + +static void test_CFL_format_clear_format(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix *A; + OK(make_simple_matrix(&A, 5)); + OK(matrix_to_format(A, GrB_ROWMAJOR, true)); + + // 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; + NOT_OK(matrix_clear_format(A, OPT_FORMAT)); + + A->base = old_base; + A->base_row = old_base_row; + A->base_col = old_base_col; + OK(CFL_matrix_free(&A)); + + teardown(); +#endif +} + +static void test_CFL_format_dup_format(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix *A, *B; + OK(make_simple_matrix(&A, 5)); + OK(CFL_matrix_create(&B, 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; + NOT_OK(matrix_dup_format(A, B, OPT_FORMAT)); + + A->base = old_base; + A->base_row = old_base_row; + A->base_col = old_base_col; + OK(CFL_matrix_free(&A)); + OK(CFL_matrix_free(&B)); + + // Without optimization flag + make_simple_matrix(&A, 5); + CFL_matrix_create(&B, 5, 5); + OK(matrix_to_format(A, GrB_COLMAJOR, true)); + OK(matrix_dup_format(A, B, 0)); + OK(CFL_matrix_free(&A)); + OK(CFL_matrix_free(&B)); + + teardown(); +#endif +} + +static void test_CFL_format_mxm_second_greather_then_k(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix *A, *B, *C; + CFL_matrix_create(&A, 5, 5); + make_ones_matrix(&B, 5); + CFL_matrix_create(&C, 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); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + OK(free_matrix(&C)); + + teardown(); +#endif +} + +static void test_CFL_format_wise_when_both(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix *A, *B; + + OK(CFL_matrix_create(&A, 5, 5)); + OK(make_ones_matrix(&B, 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; + NOT_OK(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(CFL_matrix_free(&A)); + OK(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(&A, 5, 5); + + OK(CFL_clear(A, OPT_EMPTY)); + TEST_CHECK(A->nvals == 0); + + OK(CFL_matrix_free(&A)); + + teardown(); +#endif +} + +static void test_CFL_empty_dup(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix *A, *B; + CFL_matrix_create(&A, 5, 5); + make_ones_matrix(&B, 5); + + OK(CFL_dup(A, B, OPT_EMPTY)); + TEST_CHECK(A->nvals == 25); + + OK(CFL_matrix_free(&A)); + OK(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 + CFL_matrix_create(&A, 5, 5); + make_ones_matrix(&B, 5); + make_ones_matrix(&C, 5); + + OK(CFL_mxm(C, A, B, true, false, OPT_EMPTY)); + TEST_CHECK(C->nvals == 25); + + OK(CFL_matrix_free(&A)); + OK(CFL_matrix_free(&B)); + OK(CFL_matrix_free(&C)); + + // without accum + CFL_matrix_create(&A, 5, 5); + make_ones_matrix(&B, 5); + make_ones_matrix(&C, 5); + + OK(CFL_mxm(C, A, B, false, false, OPT_EMPTY)); + TEST_CHECK(C->nvals == 0); + + OK(CFL_matrix_free(&A)); + OK(CFL_matrix_free(&B)); + OK(CFL_matrix_free(&C)); + + // output empty + CFL_matrix_create(&A, 5, 5); + make_ones_matrix(&B, 5); + CFL_matrix_create(&C, 5, 5); + + OK(CFL_mxm(C, A, B, false, false, OPT_EMPTY)); + TEST_CHECK(C->nvals == 0); + + OK(CFL_matrix_free(&A)); + OK(CFL_matrix_free(&B)); + OK(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++) { + + is_first_empty ? CFL_matrix_create(&A, 5, 5) + : make_ones_matrix(&A, 5); + is_second_empty ? CFL_matrix_create(&B, 5, 5) + : make_ones_matrix(&B, 5); + is_output_empty ? CFL_matrix_create(&C, 5, 5) + : make_ones_matrix(&C, 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); + } + + OK(CFL_matrix_free(&A)); + OK(CFL_matrix_free(&B)); + OK(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++) { + is_first_empty ? CFL_matrix_create(&A, 5, 5) : make_ones_matrix(&A, 5); + is_second_empty ? CFL_matrix_create(&B, 5, 5) : make_ones_matrix(&B, 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); + } + + OK(CFL_matrix_free(&A)); + OK(CFL_matrix_free(&B)); + } + } + } + + teardown(); +#endif +} + +static void test_CFL_empty_rsub_both_empty(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix *A, *B; + OK(CFL_matrix_create(&A, 5, 5)); + OK(CFL_matrix_create(&B, 5, 5)); + + OK(CFL_rsub(A, B, OPT_EMPTY)); + TEST_CHECK(A->nvals == 0); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + + teardown(); +#endif +} + +//------------------------------------------------------------------------------ +// Lazy optimization +//------------------------------------------------------------------------------ + +static void test_CFL_lazy_create(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix *A; + CFL_matrix_create_lazy(&A, 5, 5); + TEST_CHECK(A->base_matrices != NULL); + TEST_CHECK(A->base_matrices_count == 1); + TEST_CHECK(A->is_lazy == true); + TEST_CHECK(A->nvals == 0); + + OK(free_matrix(&A)); + + teardown(); +#endif +} + +// Create lazy matrix +// 1 0 0 0 .. 0 (count of 1: 1) +// 1 1 0 0 .. 0 (count of 1: 11) +// 1 1 1 0 .. 0 (count of 1: 111) +// . . . . .. . +// 0 0 0 0 .. 0 +static GrB_Info make_lazy_matrix(Matrix **M, 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; + OK(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 + OK(GrB_Matrix_setElement_BOOL(base_matrix, true, + (i + 1) % base_matrices_count, j)); + } + + OK(CFL_matrix_from_base(&base_matrices[(i + 1) % base_matrices_count], + base_matrix)); + } + + Matrix *result = *M; + OK(CFL_matrix_create_lazy(&result, n, n)); + + result->base_matrices_count = base_matrices_count; + for (size_t i = 0; i < base_matrices_count; i++) { + result->base_matrices[i] = base_matrices[i]; + } + + OK(CFL_matrix_update(result)); + *M = result; + + return GrB_SUCCESS; +} + +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, *B, *C; + OK(make_lazy_matrix(&A, 4)); + OK(make_ones_matrix(&B, 1111)); + OK(CFL_matrix_create(&C, 1111, 1111)); + + Matrix *A_base, *B_base, *C_base; + OK(CFL_matrix_to_base(&A_base, A, OPT_LAZY)); + OK(CFL_matrix_to_base(&B_base, B, OPT_LAZY)); + OK(CFL_matrix_to_base(&C_base, C, OPT_LAZY)); + + OK(CFL_mxm(C, A, B, is_accum, is_reverse, OPT_LAZY)); + OK(CFL_mxm(C_base, A_base, B_base, is_accum, is_reverse, 0)); + + TEST_CHECK(compare_matrices(C, C_base)); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + OK(free_matrix(&C)); + OK(free_matrix(&A_base)); + OK(free_matrix(&B_base)); + OK(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, *B, *C; + make_lazy_matrix(&A, 4); + CFL_matrix_create(&B, 1111, 1111); + CFL_matrix_create(&C, 1111, 1111); + + Matrix *A_base, *B_base, *C_base; + OK(CFL_matrix_to_base(&A_base, A, OPT_LAZY)); + OK(CFL_matrix_to_base(&B_base, B, OPT_LAZY)); + OK(CFL_matrix_to_base(&C_base, C, OPT_LAZY)); + + OK(CFL_mxm(C, A, B, is_accum, is_reverse, OPT_LAZY)); + OK(CFL_mxm(C_base, A_base, B_base, is_accum, is_reverse, 0)); + + TEST_CHECK(compare_matrices(C, C_base)); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + OK(free_matrix(&C)); + OK(free_matrix(&A_base)); + OK(free_matrix(&B_base)); + OK(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, *B, *C; + make_lazy_matrix(&A, 4); + CFL_matrix_create(&B, 1111, 1111); + CFL_matrix_create(&C, 1111, 1111); + + Matrix *A_base, *B_base, *C_base; + OK(CFL_matrix_to_base(&A_base, A, OPT_LAZY)); + OK(CFL_matrix_to_base(&B_base, B, OPT_LAZY)); + OK(CFL_matrix_to_base(&C_base, C, OPT_LAZY)); + + OK(CFL_wise(C, A, B, is_accum, OPT_LAZY)); + OK(CFL_wise(C_base, A_base, B_base, is_accum, 0)); + + TEST_CHECK(compare_matrices(A, C_base)); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + OK(free_matrix(&C)); + OK(free_matrix(&A_base)); + OK(free_matrix(&B_base)); + OK(free_matrix(&C_base)); + } + + for (size_t is_accum = 0; is_accum < 2; is_accum++) { + Matrix *A, *B, *C; + CFL_matrix_create(&A, 1111, 1111); + make_lazy_matrix(&B, 4); + CFL_matrix_create(&C, 1111, 1111); + + Matrix *A_base, *B_base, *C_base; + OK(CFL_matrix_to_base(&A_base, A, OPT_LAZY)); + OK(CFL_matrix_to_base(&B_base, B, OPT_LAZY)); + OK(CFL_matrix_to_base(&C_base, C, OPT_LAZY)); + + OK(CFL_wise(C, A, B, is_accum, OPT_LAZY)); + OK(CFL_wise(C_base, A_base, B_base, is_accum, 0)); + + TEST_CHECK(compare_matrices(C, C_base)); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + OK(free_matrix(&C)); + OK(free_matrix(&A_base)); + OK(free_matrix(&B_base)); + OK(free_matrix(&C_base)); + } + + teardown(); +#endif +} + +// Components occupy disjoint ranges of positions in a 64-by-64 matrix. +static void make_component(Matrix **matrix, size_t offset, size_t count) { + OK(CFL_matrix_create(matrix, 64, 64)); + for (size_t i = offset; i < offset + count; i++) { + OK(GrB_Matrix_setElement_BOOL((*matrix)->base, true, i / 64, i % 64)); + } + OK(CFL_matrix_update(*matrix)); +} + +static void make_lazy_components(Matrix **matrix, const size_t *counts, size_t count) { + OK(CFL_matrix_create_lazy(matrix, 64, 64)); + OK(CFL_matrix_free(&(*matrix)->base_matrices[0])); + (*matrix)->base_matrices_count = 0; + OK(LAGraph_Realloc((void **)&(*matrix)->base_matrices, count, 1, + sizeof(Matrix *), msg)); + size_t offset = 0; + for (size_t i = 0; i < count; i++) { + make_component(&(*matrix)->base_matrices[i], offset, counts[i]); + (*matrix)->base_matrices_count++; + offset += counts[i]; + } + OK(CFL_matrix_update(*matrix)); +} + +static void test_CFL_lazy_merge_chain(void) { +#if LAGRAPH_SUITESPARSE + setup(); + for (int8_t mask = OPT_LAZY; mask < OPT_N_FLAGS; mask++) { + if (!(mask & OPT_LAZY)) continue; + const size_t counts[] = {100, 1100}; + Matrix *A, *B, *expected; + make_lazy_components(&A, counts, 2); + make_component(&B, 1200, 90); + make_component(&expected, 0, 1290); + + OK(CFL_wise(A, A, B, false, mask)); + TEST_CHECK(A->base_matrices_count == 1); + TEST_CHECK(A->nvals == 1290); + TEST_CHECK(compare_matrices(A, expected)); + + OK(CFL_matrix_free(&A)); + OK(CFL_matrix_free(&B)); + OK(CFL_matrix_free(&expected)); + } + teardown(); +#endif +} + +static void test_CFL_lazy_many_components(void) { +#if LAGRAPH_SUITESPARSE + setup(); + for (int8_t mask = OPT_LAZY; mask < OPT_N_FLAGS; mask++) { + if (!(mask & OPT_LAZY)) continue; + size_t counts[60]; + for (size_t i = 0; i < 60; i++) counts[i] = 1; + Matrix *A, *B, *base, *expected; + make_lazy_components(&A, counts, 60); + + OK(matrix_combine_lazy(A, 0, mask)); + TEST_CHECK(A->base_matrices_count == 60); + OK(CFL_matrix_to_base(&base, A, mask)); + make_component(&expected, 0, 60); + TEST_CHECK(compare_matrices(base, expected)); + OK(CFL_matrix_free(&base)); + OK(CFL_matrix_free(&expected)); + + make_component(&B, 60, 1000); + OK(CFL_wise(A, A, B, false, mask)); + TEST_CHECK(A->base_matrices_count == 61); + TEST_CHECK(A->nvals == 1060); + make_component(&expected, 0, 1060); + TEST_CHECK(compare_matrices(A, expected)); + + OK(CFL_matrix_free(&A)); + OK(CFL_matrix_free(&B)); + OK(CFL_matrix_free(&expected)); + } + teardown(); +#endif +} + +// Checks that the extracted matrix is the same GrB_Matrix as the source one (it is +// moved, not copied) and that it has the expected entries +static void check_extracted_base(Matrix *A, GrB_Matrix source, size_t count, + int8_t mask) { + GrB_Matrix base = NULL; + OK(CFL_matrix_extract_base(&base, A, mask)); + TEST_CHECK(base == source); + TEST_MSG("Mask: %d", mask); + + Matrix *result, *expected; + OK(CFL_matrix_from_base(&result, base)); + make_component(&expected, 0, count); + TEST_CHECK(compare_matrices(result, expected)); + TEST_MSG("Mask: %d", mask); + + OK(CFL_matrix_free(&result)); + OK(CFL_matrix_free(&expected)); +} + +static void test_CFL_extract_base(void) { +#if LAGRAPH_SUITESPARSE + setup(); + for (int8_t mask = 0; mask < OPT_N_FLAGS; mask++) { + // regular matrix + { + Matrix *A; + make_component(&A, 0, 100); + check_extracted_base(A, A->base, 100, mask); + TEST_CHECK(A->base == NULL); + OK(CFL_matrix_free(&A)); + } + + // matrix with both formats + { + Matrix *A; + make_component(&A, 0, 100); + OK(matrix_to_format(A, GrB_COLMAJOR, true)); + TEST_CHECK(A->is_both); + check_extracted_base(A, A->base_col, 100, mask); + TEST_CHECK(A->base_row == NULL && A->base_col == NULL); + OK(CFL_matrix_free(&A)); + } + + // lazy matrix is merged into its largest base matrix + { + const size_t counts[] = {10, 1000, 100}; + Matrix *A; + make_lazy_components(&A, counts, 3); + check_extracted_base(A, A->base_matrices[1]->base, 1110, mask); + OK(CFL_matrix_free(&A)); + } + } + teardown(); +#endif +} + +static void *fail_lazy_realloc(void *p, size_t size) { + (void)p; + (void)size; + return NULL; +} + +static void test_CFL_lazy_growth_failure(void) { +#if LAGRAPH_SUITESPARSE + OK(LG_brutal_setup(msg)); + Matrix *A, *B; + OK(CFL_matrix_create_lazy(&A, 64, 64)); + make_component(&B, 0, 1); + // Warm up GraphBLAS before checking the allocation count. + OK(CFL_wise(A, A, B, false, OPT_LAZY)); + int64_t baseline = LG_nmalloc; + LAGraph_Realloc_function = fail_lazy_realloc; + GrB_Info info = CFL_wise(A, A, B, false, OPT_LAZY); + LAGraph_Realloc_function = LG_brutal_realloc; + + TEST_CHECK(info == GrB_OUT_OF_MEMORY); + TEST_CHECK(LG_nmalloc == baseline); + TEST_CHECK(A->base_matrices_count == 1); + TEST_CHECK(A->nvals == 1); + + OK(CFL_matrix_free(&A)); + OK(CFL_matrix_free(&B)); + + const size_t counts[] = {1, 1000}; + make_lazy_components(&A, counts, 2); + Matrix *base = NULL; + OK(CFL_matrix_to_base(&base, A, OPT_LAZY)); + OK(CFL_matrix_free(&base)); + baseline = LG_nmalloc; + LAGraph_Realloc_function = fail_lazy_realloc; + info = CFL_matrix_to_base(&base, A, OPT_LAZY); + LAGraph_Realloc_function = LG_brutal_realloc; + + TEST_CHECK(info == GrB_OUT_OF_MEMORY); + TEST_CHECK(base == NULL); + TEST_CHECK(LG_nmalloc == baseline); + TEST_CHECK(A->base_matrices_count == 2); + TEST_CHECK(A->nvals == 1001); + OK(CFL_matrix_free(&A)); + OK(LG_brutal_teardown(msg)); +#endif +} + +static void test_CFL_lazy_rsub(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + // { + // Matrix *A, *B; + // OK(make_lazy_matrix(&A, 4)); + // OK(CFL_matrix_create(&B, 1111, 1111)); + + // Matrix *A_base, *B_base; + // OK(CFL_matrix_to_base(&A_base, A, OPT_LAZY)); + // OK(CFL_matrix_to_base(&B_base, B, OPT_LAZY)); + + // OK(CFL_rsub(A, B, OPT_LAZY)); + // OK(CFL_rsub(A_base, B_base, 0)); + + // TEST_CHECK(compare_matrices(A, A_base)); + + // OK(free_matrix(&A)); + // OK(free_matrix(&B)); + // OK(free_matrix(&A_base)); + // OK(free_matrix(&B_base)); + // } + + { + Matrix *A, *B; + OK(CFL_matrix_create(&A, 1111, 1111)); + OK(make_lazy_matrix(&B, 4)); + + Matrix *A_base, *B_base; + OK(CFL_matrix_to_base(&A_base, A, OPT_LAZY)); + OK(CFL_matrix_to_base(&B_base, B, OPT_LAZY)); + + OK(CFL_rsub(A, B, OPT_LAZY)); + OK(CFL_rsub(A_base, B_base, 0)); + + TEST_CHECK(compare_matrices(A, A_base)); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + OK(free_matrix(&A_base)); + OK(free_matrix(&B_base)); + } + + teardown(); +#endif +} + +//------------------------------------------------------------------------------ +// Block optimization +//------------------------------------------------------------------------------ + +// Create block vector +// 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 . +// . . . . . . . +static GrB_Info make_block_vector(Matrix **M, size_t n, size_t alpha, + Matrix_block format) { + GrB_Matrix result; + switch (format) { + case CELL: + OK(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; + + OK(GrB_Matrix_setElement_BOOL(result, value, i, j)); + } + } + + OK(CFL_matrix_from_base(M, result)); + break; + case VEC_VERT: + OK(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; + + OK(GrB_Matrix_setElement_BOOL(result, value, i, j)); + } + } + + OK(CFL_matrix_from_base(M, result)); + break; + case VEC_HORIZ: + OK(GrB_Matrix_new(&result, GrB_BOOL, n, n * alpha)); + for (size_t i = 0; i < n; i++) { + for (size_t j = 0; j < n * alpha; j++) { + bool value = j == (i + 5 % n) || j == ((i + 6) % n) ? true : false; + if (!value) + continue; + + OK(GrB_Matrix_setElement_BOOL(result, value, i, j)); + TEST_MSG("%d %d %d %d %d", value, i, j, n, n * alpha); + } + } + + OK(CFL_matrix_from_base(M, result)); + break; + } + + return GrB_SUCCESS; +} + +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, *B; + OK(make_block_vector(&A, 20, 20, first_format)); + OK(make_block_vector(&B, 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; + CFL_matrix_create(&C, nrows, ncols); + } else { + CFL_matrix_create(&C, first_format == CELL ? 20 : 20 * 20, + second_format == CELL ? 20 : 20 * 20); + } + + OK(CFL_mxm(C, A, B, is_accum, is_reverse, OPT_BLOCK)); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + OK(free_matrix(&C)); + } + } + } + } + + teardown(); +#endif +} + +static void test_CFL_block_mxm_cells_into_vector(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = OPT_BLOCK; mask < OPT_N_FLAGS; mask++) { + for (size_t format = VEC_HORIZ; format <= VEC_VERT; format++) { + 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 is_empty = 0; is_empty < 2; is_empty++) { + Matrix *A, *B, *C, *expected; + if (mask & OPT_LAZY) { + OK(CFL_matrix_create_lazy(&A, 3, 3)); + } else { + OK(CFL_matrix_create(&A, 3, 3)); + } + GrB_Matrix base_A = A->is_lazy ? A->base_matrices[0]->base : A->base; + if (!is_empty) { + OK(GrB_Matrix_setElement_BOOL(base_A, true, 0, 1)); + } + OK(CFL_matrix_create(&B, 3, 3)); + OK(GrB_Matrix_setElement_BOOL(B->base, true, 1, 2)); + OK(GrB_Matrix_setElement_BOOL(B->base, true, 2, 0)); + + GrB_Index nrows = format == VEC_VERT ? 9 : 3; + GrB_Index ncols = format == VEC_HORIZ ? 9 : 3; + OK(CFL_matrix_create(&C, nrows, ncols)); + OK(CFL_matrix_create(&expected, nrows, ncols)); + for (size_t i = 0; i < 3; i++) { + GrB_Index row_offset = format == VEC_VERT ? 3 * i : 0; + GrB_Index col_offset = format == VEC_HORIZ ? 3 * i : 0; + OK(GrB_Matrix_setElement_BOOL(C->base, true, row_offset, col_offset)); + if (is_accum) { + OK(GrB_Matrix_setElement_BOOL(expected->base, true, + row_offset, col_offset)); + } + if (!is_empty) { + OK(GrB_Matrix_setElement_BOOL(expected->base, true, + row_offset + (is_reverse ? 2 : 0), + col_offset + (is_reverse ? 1 : 2))); + } + } + + GrB_Info info = CFL_mxm(C, A, B, is_accum, is_reverse, mask); + OK(info); + if (info == GrB_SUCCESS) { + TEST_CHECK(compare_matrices(C, expected)); + } + TEST_MSG("Mask: %zu, format: %zu, accum: %zu, reverse: %zu, empty: %zu", + mask, format, is_accum, is_reverse, is_empty); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + OK(free_matrix(&C)); + OK(free_matrix(&expected)); + } + } + } + } + } + + teardown(); +#endif +} + +static void test_CFL_block_dup(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix *A, *B; + + { + make_block_vector(&A, 20, 20, VEC_HORIZ); + CFL_matrix_create(&B, 20 * 20, 20); + + OK(CFL_dup(B, A, OPT_BLOCK)); + TEST_CHECK(B->nvals == A->nvals); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + } + + { + make_block_vector(&A, 20, 20, VEC_VERT); + CFL_matrix_create(&B, 20 * 20, 20); + + OK(CFL_dup(B, A, OPT_BLOCK)); + TEST_CHECK(B->nvals == A->nvals); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + } + + { + make_block_vector(&A, 20, 20, VEC_HORIZ); + CFL_matrix_create(&B, 20, 20 * 20); + + OK(CFL_dup(B, A, OPT_BLOCK)); + TEST_CHECK(B->nvals == A->nvals); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + } + + { + make_block_vector(&A, 20, 20, VEC_VERT); + CFL_matrix_create(&B, 20, 20 * 20); + + OK(CFL_dup(B, A, OPT_BLOCK)); + TEST_CHECK(B->nvals == A->nvals); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + } + + teardown(); +#endif +} + +static void test_CFL_block_wise(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + { + Matrix *A, *B, *C; + OK(make_block_vector(&A, 20, 20, VEC_HORIZ)); + OK(make_block_vector(&B, 20, 20, VEC_HORIZ)); + OK(make_block_vector(&C, 20, 20, VEC_HORIZ)); + + TEST_CHECK(CFL_wise(C, A, B, false, OPT_BLOCK) == GrB_INVALID_VALUE); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + OK(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, *B; + OK(make_block_vector(&A, 20, 20, first_format)); + OK(make_block_vector(&B, 20, 20, second_format)); + + OK(CFL_wise(A, A, B, is_accum, OPT_BLOCK)); + + OK(free_matrix(&A)); + OK(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, *B; + OK(make_block_vector(&A, 20, 20, first_format)); + OK(make_block_vector(&B, 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)); + } + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + } + } + } + } + + teardown(); +#endif +} + +static void test_CFL_block_hyper_rotate(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + // lazy + { + Matrix *A; + OK(CFL_matrix_create_lazy(&A, 20 * 20, 20)); + + OK(CFL_matrix_free(&A->base_matrices[0])); + A->base_matrices_count = 0; + + Matrix *A0, *A1, *A2; + CFL_matrix_create(&A0, 20 * 20, 20); + CFL_matrix_create(&A1, 20 * 20, 20); + CFL_matrix_create(&A2, 20 * 20, 20); + + A->base_matrices[0] = A0; + A->base_matrices[1] = A1; + A->base_matrices[2] = A2; + A->base_matrices_count = 3; + OK(CFL_matrix_update(A)); + + OK(block_matrix_hyper_rotate_i(A, VEC_HORIZ)); + TEST_CHECK(A->nrows == 20); + TEST_CHECK(A->ncols == 20 * 20); + + OK(free_matrix(&A)); + // free_matrix(&A0); + // free_matrix(&A1); + // free_matrix(&A2); + } + + // cell + { + Matrix *A; + CFL_matrix_create(&A, 5, 5); + + OK(block_matrix_hyper_rotate_i(A, VEC_HORIZ)); + + OK(free_matrix(&A)); + } + + // VERT + { + Matrix *A; + make_block_vector(&A, 5, 5, VEC_VERT); + TEST_CHECK(A->nvals != 0); + TEST_CHECK(A->is_lazy != true); + + OK(block_matrix_hyper_rotate_i(A, VEC_HORIZ)); + + OK(free_matrix(&A)); + } + + // Format + { + Matrix *A; + make_block_vector(&A, 5, 5, VEC_VERT); + matrix_to_format(A, GrB_COLMAJOR, true); + + OK(block_matrix_hyper_rotate_i(A, VEC_HORIZ)); + TEST_CHECK(A->nvals != 0); + TEST_CHECK(A->is_both = true); + + GrB_Index nrows_A, ncols_A, nrows_B, ncols_B; + OK(GrB_Matrix_nrows(&nrows_A, A->base_row)); + OK(GrB_Matrix_ncols(&ncols_A, A->base_row)); + OK(GrB_Matrix_nrows(&nrows_B, A->base_col)); + OK(GrB_Matrix_ncols(&ncols_B, A->base_col)); + + TEST_CHECK(nrows_A == nrows_B); + TEST_CHECK(ncols_A == ncols_B); + + OK(block_matrix_hyper_rotate_i(A, VEC_VERT)); + + OK(GrB_Matrix_nrows(&nrows_A, A->base_row)); + OK(GrB_Matrix_ncols(&ncols_A, A->base_row)); + OK(GrB_Matrix_nrows(&nrows_B, A->base_col)); + OK(GrB_Matrix_ncols(&ncols_B, A->base_col)); + + TEST_CHECK(nrows_A == nrows_B); + TEST_CHECK(ncols_A == ncols_B); + + OK(free_matrix(&A)); + } + + teardown(); +#endif +} + +static void test_CFL_block_to_diag(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + { + Matrix *A, *B; + OK(CFL_matrix_create(&A, 5, 5)); + OK(CFL_matrix_create(&B, 5, 5)); + + GrB_Index nvals_before, nvals_after; + + OK(GrB_Matrix_nvals(&nvals_before, A->base)); + OK(block_matrix_to_diag(A, B)); + OK(GrB_Matrix_nvals(&nvals_after, A->base)); + + // A didn't change + TEST_CHECK(nvals_after == nvals_before); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + } + + teardown(); +#endif +} + +static void test_CFL_block_to_reduce(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + Matrix *A, *B; + OK(CFL_matrix_create(&A, 5, 5)); + OK(make_ones_matrix(&B, 5)); + + // cell matrix cause just dup + OK(block_matrix_reduce(A, B, OPT_BLOCK)); + TEST_CHECK(A->nvals == B->nvals); + + OK(free_matrix(&A)); + OK(free_matrix(&B)); + + teardown(); +#endif +} + +//------------------------------------------------------------------------------ +// TEST LIST +//------------------------------------------------------------------------------ + +TEST_LIST = { + {"test_CFL_format_wise_when_both", test_CFL_format_wise_when_both}, + {"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_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_merge_chain", test_CFL_lazy_merge_chain}, + {"test_CFL_lazy_many_components", test_CFL_lazy_many_components}, + {"test_CFL_extract_base", test_CFL_extract_base}, +#if !defined(GRAPHBLAS_HAS_CUDA) + {"test_CFL_lazy_growth_failure", test_CFL_lazy_growth_failure}, +#endif + {"test_CFL_lazy_rsub", test_CFL_lazy_rsub}, + {"test_CFL_block_mxm", test_CFL_block_mxm}, + {"test_CFL_block_mxm_cells_into_vector", test_CFL_block_mxm_cells_into_vector}, + {"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}}; diff --git a/experimental/test/test_CFL_reachability_adv.c b/experimental/test/test_CFL_reachability_adv.c new file mode 100644 index 0000000000..6279f884cb --- /dev/null +++ b/experimental/test/test_CFL_reachability_adv.c @@ -0,0 +1,1286 @@ +//------------------------------------------------------------------------------ +// LAGraph/experimental/test/LAGraph_CFL_reachability.c: test cases for 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 run_algorithm(mask) \ + LAGraph_CFL_reachability_adv(outputs, adj_matrices, \ + grammar.terms_count + grammar.nonterms_count, \ + grammar.rules, grammar.rules_count, msg, mask) + +#define check_error(error, mask) \ + { \ + retval = run_algorithm(mask); \ + TEST_CHECK(retval == error); \ + TEST_MSG("retval = %d (%s)", retval, msg); \ + } + +#define check_result(result) \ + { \ + char *expected = output_to_str(0); \ + TEST_CHECK(strcmp(result, expected) == 0); \ + TEST_MSG("Wrong result. Mask: %x. Actual: %s", mask, expected); \ + LAGraph_Free((void **)&expected, msg); \ + } + +typedef struct { + size_t nonterms_count; + size_t terms_count; + size_t rules_count; + LAGraph_rule_EWCNF *rules; +} grammar_t; + +GrB_Matrix *adj_matrices = NULL; +int n_adj_matrices = 0; +GrB_Matrix *outputs = NULL; +grammar_t grammar = {0, 0, 0, NULL}; +char msg[LAGRAPH_MSG_LEN]; + +void setup() { LAGraph_Init(msg); } + +void teardown(void) { LAGraph_Finalize(msg); } + +void init_outputs() { + LAGraph_Calloc((void **)&outputs, grammar.nonterms_count + grammar.terms_count, + sizeof(GrB_Matrix), msg); +} + +char *output_to_str(size_t nonterm) { + GrB_Index nnz = 0; + GrB_Info info = GrB_Matrix_nvals(&nnz, outputs[nonterm]); + OK(info); + TEST_MSG("Error: %d\n", info); + 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); + + GrB_set(outputs[nonterm], GrB_ROWMAJOR, GrB_STORAGE_ORIENTATION_HINT); + OK(GrB_Matrix_extractTuples_BOOL(row, col, val, &nnz, outputs[nonterm])); + + // 11 - size of " (%ld, %ld)" + 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() { + if (adj_matrices != NULL) { + for (size_t i = 0; i < n_adj_matrices; i++) { + GrB_free(&adj_matrices[i]); + } + } + LAGraph_Free((void **)&adj_matrices, msg); + + if (outputs != NULL) { + for (size_t i = 0; i < grammar.nonterms_count + grammar.terms_count; i++) { + GrB_free(&outputs[i]); + } + } + LAGraph_Free((void **)&outputs, msg); + + LAGraph_Free((void **)&grammar.rules, msg); + grammar = (grammar_t){0, 0, 0, NULL}; +} + +//==================== +// Grammars +//==================== + +// S -> aSb | ab in WCNF +// +// Terms: [4 a] [5 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 4 -1 0] +// B -> b [2 5 -1 0] +void init_grammar_aSb() { + LAGraph_rule_EWCNF *rules = NULL; + LAGraph_Calloc((void **)&rules, 5, sizeof(LAGraph_rule_EWCNF), msg); + + rules[0] = (LAGraph_rule_EWCNF){0, 1, 2, 0, 0}; + rules[1] = (LAGraph_rule_EWCNF){0, 1, 3, 0, 0}; + rules[2] = (LAGraph_rule_EWCNF){3, 0, 2, 0, 0}; + rules[3] = (LAGraph_rule_EWCNF){1, 4, -1, 0, 0}; + rules[4] = (LAGraph_rule_EWCNF){2, 5, -1, 0, 0}; + + grammar = (grammar_t){ + .nonterms_count = 4, .terms_count = 2, .rules_count = 5, .rules = rules}; +} + +// Terms: [0 a] [1 b_0] [2 b_1] +// Nonterms: [3 S_0] [4 S_1] +// S_i -> a b_i [3 1 2 2 LAGraph_EWNCF_INDEX_NONTERM | LAGraph_EWNCF_INDEX_PROD_B] +void init_indexed_grammar() { + LAGraph_rule_EWCNF *rules = NULL; + LAGraph_Calloc((void **)&rules, 1, sizeof(LAGraph_rule_EWCNF), msg); + + rules[0] = (LAGraph_rule_EWCNF){ + 3, 0, 1, 2, LAGraph_EWNCF_INDEX_NONTERM | LAGraph_EWNCF_INDEX_PROD_B}; + + grammar = (grammar_t){ + .nonterms_count = 2, .terms_count = 3, .rules_count = 1, .rules = rules}; +} + +// Terms: [0 a_0] [1 a_1] [2 b_0] [3 b_1] +// Nonterms: [4 S_0] [5 S_1] +// S_i -> a_i b_i [4 0 2 2 LAGraph_EWNCF_INDEX_NONTERM | LAGraph_EWNCF_INDEX_PROD_A | +// LAGraph_EWNCF_INDEX_PROD_B] +void init_indexed_grammar_2() { + LAGraph_rule_EWCNF *rules = NULL; + LAGraph_Calloc((void **)&rules, 1, sizeof(LAGraph_rule_EWCNF), msg); + + rules[0] = + (LAGraph_rule_EWCNF){4, 0, 2, 2, + LAGraph_EWNCF_INDEX_NONTERM | LAGraph_EWNCF_INDEX_PROD_A | + LAGraph_EWNCF_INDEX_PROD_B}; + + grammar = (grammar_t){ + .nonterms_count = 2, .terms_count = 4, .rules_count = 1, .rules = rules}; +} + +// Terms: [0 a_0] [1 a_1] +// Nonterms: [2 S_0] [3 S_1] +// S_i -> a_i [2 0 -1 2 LAGraph_EWNCF_INDEX_NONTERM | LAGraph_EWNCF_INDEX_PROD_A] +void init_indexed_grammar_simple() { + LAGraph_rule_EWCNF *rules = NULL; + LAGraph_Calloc((void **)&rules, 1, sizeof(LAGraph_rule_EWCNF), msg); + + rules[0] = (LAGraph_rule_EWCNF){ + 2, 0, -1, 2, LAGraph_EWNCF_INDEX_NONTERM | LAGraph_EWNCF_INDEX_PROD_A}; + + grammar = (grammar_t){ + .nonterms_count = 2, .terms_count = 2, .rules_count = 1, .rules = rules}; +} + +// Terms: [0 a_0] [1 a_1] +// Nonterms: [2 S_0] [3 S_1] +// S_0 -> a_0 [2 0 -1 0 0] +// S_1 -> a_1 [3 1 -1 0 0] +void init_indexed_grammar_simple_exloded() { + LAGraph_rule_EWCNF *rules = NULL; + LAGraph_Calloc((void **)&rules, 2, sizeof(LAGraph_rule_EWCNF), msg); + + rules[0] = (LAGraph_rule_EWCNF){2, 0, -1, 0, 0}; + rules[1] = (LAGraph_rule_EWCNF){3, 1, -1, 0, 0}; + + grammar = (grammar_t){ + .nonterms_count = 2, .terms_count = 2, .rules_count = 2, .rules = rules}; +} + +// S -> aS | a | eps in WCNF +// +// Terms: [1 a] +// Nonterms: [0 S] +// S -> SS [0 0 0 0] +// S -> a [0 1 -1 0] +// S -> eps [0 -1 -1 0] +void init_grammar_aS() { + LAGraph_rule_EWCNF *rules = NULL; + LAGraph_Calloc((void **)&rules, 3, sizeof(LAGraph_rule_EWCNF), msg); + + rules[0] = (LAGraph_rule_EWCNF){0, 0, 0, 0}; + rules[1] = (LAGraph_rule_EWCNF){0, 1, -1, 0}; + rules[2] = (LAGraph_rule_EWCNF){0, -1, -1, 0}; + + grammar = (grammar_t){ + .nonterms_count = 1, .terms_count = 1, .rules_count = 3, .rules = rules}; +} + +// Complex grammar +// aaaabbbb or aaabbb +// +// Terms: [25 a] [26 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 25 -1 0] +// S8 -> a [8 25 -1 0] +// S9 -> a [9 25 -1 0] +// S10 -> a [10 25 -1 0] +// S11 -> b [11 26 -1 0] +// S12 -> b [12 26 -1 0] +// S13 -> b [13 26 -1 0] +// S14 -> b [14 26 -1 0] +// S15 -> a [15 25 -1 0] +// S19 -> a [19 25 -1 0] +// S20 -> a [20 25 -1 0] +// S21 -> b [21 26 -1 0] +// S23 -> b [23 26 -1 0] +// S24 -> b [24 26 -1 0] +void init_grammar_complex() { + LAGraph_rule_EWCNF *rules = NULL; + LAGraph_Calloc((void **)&rules, 26, sizeof(LAGraph_rule_EWCNF), msg); + + rules[0] = (LAGraph_rule_EWCNF){0, 1, 2, 0, 0}; + rules[1] = (LAGraph_rule_EWCNF){0, 15, 16, 0, 0}; + rules[2] = (LAGraph_rule_EWCNF){1, 3, 4, 0, 0}; + rules[3] = (LAGraph_rule_EWCNF){2, 5, 6, 0, 0}; + rules[4] = (LAGraph_rule_EWCNF){3, 7, 8, 0, 0}; + rules[5] = (LAGraph_rule_EWCNF){4, 9, 10, 0, 0}; + rules[6] = (LAGraph_rule_EWCNF){5, 11, 12, 0, 0}; + rules[7] = (LAGraph_rule_EWCNF){6, 13, 14, 0, 0}; + rules[8] = (LAGraph_rule_EWCNF){16, 17, 18, 0, 0}; + rules[9] = (LAGraph_rule_EWCNF){17, 19, 20, 0, 0}; + rules[10] = (LAGraph_rule_EWCNF){18, 21, 22, 0, 0}; + rules[11] = (LAGraph_rule_EWCNF){22, 23, 24, 0, 0}; + rules[12] = (LAGraph_rule_EWCNF){7, 25, -1, 0, 0}; + rules[13] = (LAGraph_rule_EWCNF){8, 25, -1, 0, 0}; + rules[14] = (LAGraph_rule_EWCNF){9, 25, -1, 0, 0}; + rules[15] = (LAGraph_rule_EWCNF){10, 25, -1, 0, 0}; + rules[16] = (LAGraph_rule_EWCNF){11, 26, -1, 0, 0}; + rules[17] = (LAGraph_rule_EWCNF){12, 26, -1, 0, 0}; + rules[18] = (LAGraph_rule_EWCNF){13, 26, -1, 0, 0}; + rules[19] = (LAGraph_rule_EWCNF){14, 26, -1, 0, 0}; + rules[20] = (LAGraph_rule_EWCNF){15, 25, -1, 0, 0}; + rules[21] = (LAGraph_rule_EWCNF){19, 25, -1, 0, 0}; + rules[22] = (LAGraph_rule_EWCNF){20, 25, -1, 0, 0}; + rules[23] = (LAGraph_rule_EWCNF){21, 26, -1, 0, 0}; + rules[24] = (LAGraph_rule_EWCNF){23, 26, -1, 0, 0}; + rules[25] = (LAGraph_rule_EWCNF){24, 26, -1, 0, 0}; + + grammar = (grammar_t){ + .nonterms_count = 25, .terms_count = 2, .rules_count = 26, .rules = rules}; +} + +//==================== +// Graphs +//==================== + +// Graph: +// +// 0 -a-> 1 +// 1 -a-> 2 +// 2 -a-> 0 +// 0 -b-> 3 +// 3 -b-> 0 +void init_graph_double_cycle() { + n_adj_matrices = grammar.nonterms_count + grammar.terms_count; + LAGraph_Calloc((void **)&adj_matrices, n_adj_matrices, sizeof(GrB_Matrix), msg); + size_t N = 4; + + GrB_Matrix adj_matrix_a, adj_matrix_b; + OK(GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, N, N)); + OK(GrB_Matrix_new(&adj_matrix_b, GrB_BOOL, N, N)); + + 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[grammar.nonterms_count + 0] = adj_matrix_a; + adj_matrices[grammar.nonterms_count + 1] = adj_matrix_b; + + for (size_t i = 0; i < grammar.nonterms_count; i++) { + GrB_Matrix_new(&adj_matrices[i], GrB_BOOL, N, N); + } +} + +// 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() { + n_adj_matrices = grammar.nonterms_count + grammar.terms_count; + LAGraph_Calloc((void **)&adj_matrices, n_adj_matrices, sizeof(GrB_Matrix), msg); + size_t N = 8; + + GrB_Matrix adj_matrix_a, adj_matrix_b; + OK(GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, N, N)); + OK(GrB_Matrix_new(&adj_matrix_b, GrB_BOOL, N, N)); + + 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[grammar.nonterms_count + 0] = adj_matrix_a; + adj_matrices[grammar.nonterms_count + 1] = adj_matrix_b; + + for (size_t i = 0; i < grammar.nonterms_count; i++) { + GrB_Matrix_new(&adj_matrices[i], GrB_BOOL, N, N); + } +} + +// 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() { + n_adj_matrices = grammar.nonterms_count + grammar.terms_count; + LAGraph_Calloc((void **)&adj_matrices, n_adj_matrices, sizeof(GrB_Matrix), msg); + size_t N = 7; + + GrB_Matrix adj_matrix_a, adj_matrix_b; + OK(GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, N, N)); + OK(GrB_Matrix_new(&adj_matrix_b, GrB_BOOL, N, N)); + + 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[grammar.nonterms_count + 0] = adj_matrix_a; + adj_matrices[grammar.nonterms_count + 1] = adj_matrix_b; + + for (size_t i = 0; i < grammar.nonterms_count; i++) { + GrB_Matrix_new(&adj_matrices[i], GrB_BOOL, N, N); + } +} + +// Graph: +// +// 0 -a-> 1 +// 1 -a-> 2 +// 2 -a-> 0 +void init_graph_one_cycle() { + n_adj_matrices = grammar.nonterms_count + grammar.terms_count; + LAGraph_Calloc((void **)&adj_matrices, n_adj_matrices, sizeof(GrB_Matrix), msg); + size_t N = 3; + + GrB_Matrix adj_matrix_a; + GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, N, N); + + 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[grammar.nonterms_count + 0] = adj_matrix_a; + + for (size_t i = 0; i < grammar.nonterms_count; i++) { + GrB_Matrix_new(&adj_matrices[i], GrB_BOOL, N, N); + } +} + +// Graph: + +// 0 -a-> 1 +// 1 -a-> 2 +// 2 -b-> 3 +// 3 -b-> 4 +void init_graph_line() { + n_adj_matrices = grammar.nonterms_count + grammar.terms_count; + LAGraph_Calloc((void **)&adj_matrices, n_adj_matrices, sizeof(GrB_Matrix), msg); + size_t N = 5; + + GrB_Matrix adj_matrix_a, adj_matrix_b; + GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, N, N); + GrB_Matrix_new(&adj_matrix_b, GrB_BOOL, N, N); + + 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[grammar.nonterms_count + 0] = adj_matrix_a; + adj_matrices[grammar.nonterms_count + 1] = adj_matrix_b; + + for (size_t i = 0; i < grammar.nonterms_count; i++) { + GrB_Matrix_new(&adj_matrices[i], GrB_BOOL, N, N); + } +} + +// Graph: + +// 0 -a-> 0 +// 0 -b-> 1 +// 1 -c-> 2 +void init_graph_2() { + n_adj_matrices = grammar.nonterms_count + grammar.terms_count; + LAGraph_Calloc((void **)&adj_matrices, n_adj_matrices, sizeof(GrB_Matrix), msg); + size_t N = 3; + + GrB_Matrix adj_matrix_a, adj_matrix_b, adj_matrix_c; + GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, N, N); + GrB_Matrix_new(&adj_matrix_b, GrB_BOOL, N, N); + GrB_Matrix_new(&adj_matrix_c, GrB_BOOL, N, N); + + 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[grammar.nonterms_count + 0] = adj_matrix_a; + adj_matrices[grammar.nonterms_count + 1] = adj_matrix_b; + adj_matrices[grammar.nonterms_count + 2] = adj_matrix_c; + + for (size_t i = 0; i < grammar.nonterms_count; i++) { + GrB_Matrix_new(&adj_matrices[i], GrB_BOOL, N, N); + } +} + +// Graph: + +// 0 -a-> 1 +// 1 -a-> 0 +// 0 -b-> 0 +void init_graph_3() { + n_adj_matrices = grammar.nonterms_count + grammar.terms_count; + LAGraph_Calloc((void **)&adj_matrices, n_adj_matrices, sizeof(GrB_Matrix), msg); + size_t N = 2; + + GrB_Matrix adj_matrix_a, adj_matrix_b; + GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, N, N); + GrB_Matrix_new(&adj_matrix_b, GrB_BOOL, N, N); + + 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[grammar.nonterms_count + 0] = adj_matrix_a; + adj_matrices[grammar.nonterms_count + 1] = adj_matrix_b; + + for (size_t i = 0; i < grammar.nonterms_count; i++) { + GrB_Matrix_new(&adj_matrices[i], GrB_BOOL, N, N); + } +} + +// Graph: + +// 0 -b-> 1 +// 1 -b-> 0 +void init_graph_4() { + n_adj_matrices = grammar.nonterms_count + grammar.terms_count; + LAGraph_Calloc((void **)&adj_matrices, n_adj_matrices, sizeof(GrB_Matrix), msg); + size_t N = 2; + + GrB_Matrix adj_matrix_a, adj_matrix_b; + GrB_Matrix_new(&adj_matrix_a, GrB_BOOL, N, N); + GrB_Matrix_new(&adj_matrix_b, GrB_BOOL, N, N); + + OK(GrB_Matrix_setElement(adj_matrix_b, true, 0, 1)); + OK(GrB_Matrix_setElement(adj_matrix_b, true, 1, 0)); + + adj_matrices[grammar.nonterms_count + 0] = adj_matrix_a; + adj_matrices[grammar.nonterms_count + 1] = adj_matrix_b; + + for (size_t i = 0; i < grammar.nonterms_count; i++) { + GrB_Matrix_new(&adj_matrices[i], GrB_BOOL, N, N); + } +} + +// Graph: + +// 0 -a-> 1 +// 1 -b_1-> 2 +void init_graph_5() { + n_adj_matrices = grammar.nonterms_count + grammar.terms_count; + LAGraph_Calloc((void **)&adj_matrices, n_adj_matrices, sizeof(GrB_Matrix), msg); + size_t N = 3; + for (size_t i = 0; i < n_adj_matrices; i++) { + GrB_Matrix_new(&adj_matrices[i], GrB_BOOL, N, N); + } + + OK(GrB_Matrix_setElement(adj_matrices[0], true, 0, 1)); + OK(GrB_Matrix_setElement(adj_matrices[2], true, 1, 2)); +} + +// Graph: + +// 0 -a_1-> 1 +void init_graph_6() { + n_adj_matrices = grammar.nonterms_count + grammar.terms_count; + LAGraph_Calloc((void **)&adj_matrices, n_adj_matrices, sizeof(GrB_Matrix), msg); + size_t N = 3; + for (size_t i = 0; i < n_adj_matrices; i++) { + GrB_Matrix_new(&adj_matrices[i], GrB_BOOL, N, N); + } + + OK(GrB_Matrix_setElement(adj_matrices[1], true, 0, 1)); +} + +// 0 -a_0-> 1 +// 0 -a_0-> 2 +// 1 -a_0-> 2 +// 1 -b_0-> 0 +// 2 -b_0-> 1 +// +// Terms: [0 a_0] [2 b_0] +void init_graph_7() { + n_adj_matrices = grammar.nonterms_count + grammar.terms_count; + LAGraph_Calloc((void **)&adj_matrices, n_adj_matrices, sizeof(GrB_Matrix), msg); + size_t N = 3; + for (size_t i = 0; i < n_adj_matrices; i++) { + GrB_Matrix_new(&adj_matrices[i], GrB_BOOL, N, N); + } + + OK(GrB_Matrix_setElement(adj_matrices[0], true, 0, 1)); + OK(GrB_Matrix_setElement(adj_matrices[0], true, 0, 2)); + OK(GrB_Matrix_setElement(adj_matrices[0], true, 1, 2)); + OK(GrB_Matrix_setElement(adj_matrices[2], true, 1, 0)); + OK(GrB_Matrix_setElement(adj_matrices[2], true, 2, 1)); +} + +//==================== +// Tests with valid result +//==================== + +void test_CFL_indexed(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + GrB_Info retval; + + init_indexed_grammar(); + init_graph_5(); + init_outputs(); + + OK(run_algorithm(mask)); + TEST_MSG("MASK: %zx, ERROR: %s\n", mask, msg); + char *expected = output_to_str(4); + TEST_CHECK(strcmp("(0, 2)", expected) == 0); + TEST_MSG("Wrong result. Mask: %zx. Actual: %s", mask, expected); + LAGraph_Free((void **)&expected, msg); + + free_workspace(); + } + + teardown(); +#endif +} + +void test_CFL_indexed_2(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + GrB_Info retval; + + init_indexed_grammar_2(); + init_graph_7(); + init_outputs(); + + OK(run_algorithm(mask)); + TEST_MSG("MASK: %zx, ERROR: %s\n", mask, msg); + char *expected = output_to_str(4); + TEST_CHECK(strcmp("(0, 0) (0, 1) (1, 1)", expected) == 0); + TEST_MSG("Wrong result. Mask: %zx. Actual: %s", mask, expected); + LAGraph_Free((void **)&expected, msg); + + free_workspace(); + } + + teardown(); +#endif +} + +void test_CFL_indexed_nonterm_only(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + // Symbols: [0 S_0] [1 S_1] [2 a] [3 b]. Rules: S_i -> a b | b a. + OK(LAGraph_Calloc((void **)&grammar.rules, 2, sizeof(LAGraph_rule_EWCNF), msg)); + grammar.nonterms_count = 2; + grammar.terms_count = 2; + grammar.rules_count = 2; + grammar.rules[0] = (LAGraph_rule_EWCNF){0, 2, 3, 2, LAGraph_EWNCF_INDEX_NONTERM}; + grammar.rules[1] = (LAGraph_rule_EWCNF){0, 3, 2, 2, LAGraph_EWNCF_INDEX_NONTERM}; + + n_adj_matrices = 4; + OK(LAGraph_Calloc((void **)&adj_matrices, n_adj_matrices, sizeof(GrB_Matrix), msg)); + for (size_t i = 0; i < n_adj_matrices; i++) { + OK(GrB_Matrix_new(&adj_matrices[i], GrB_BOOL, 3, 3)); + } + OK(GrB_Matrix_setElement_BOOL(adj_matrices[2], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(adj_matrices[3], true, 1, 2)); + OK(GrB_Matrix_setElement_BOOL(adj_matrices[3], true, 2, 0)); + init_outputs(); + + GrB_Matrix expected = NULL; + OK(GrB_Matrix_new(&expected, GrB_BOOL, 3, 3)); + OK(GrB_Matrix_setElement_BOOL(expected, true, 0, 2)); + OK(GrB_Matrix_setElement_BOOL(expected, true, 2, 1)); + + GrB_Info info = run_algorithm(mask); + OK(info); + TEST_MSG("Mask: %zu", mask); + if (info == GrB_SUCCESS) { + for (size_t i = 0; i < n_adj_matrices; i++) { + if (i >= 2) { + TEST_CHECK(outputs[i] == NULL); + TEST_MSG("Mask: %zu, symbol: %zu", mask, i); + continue; + } + + bool equal = false; + OK(LAGraph_Matrix_IsEqual(&equal, outputs[i], expected, msg)); + TEST_CHECK(equal); + TEST_MSG("Mask: %zu, symbol: %zu", mask, i); + } + } + + OK(GrB_Matrix_free(&expected)); + free_workspace(); + } + + teardown(); +#endif +} + +void test_CFL_indexed_simple(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + GrB_Info retval; + + init_indexed_grammar_simple(); + init_graph_6(); + init_outputs(); + + OK(run_algorithm(mask)); + TEST_MSG("MASK: %zx, ERROR: %s\n", mask, msg); + + char *expected = output_to_str(3); + TEST_CHECK(strcmp("(0, 1)", expected) == 0); + TEST_MSG("Wrong result. Mask: %zx. Actual: %s", mask, expected); + LAGraph_Free((void **)&expected, msg); + + free_workspace(); + } + + teardown(); +#endif +} + +void test_CFL_indexed_simple_exploded(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + GrB_Info retval; + + init_indexed_grammar_simple_exloded(); + init_graph_6(); + init_outputs(); + + OK(run_algorithm(mask)); + TEST_MSG("MASK: %zx, ERROR: %s\n", mask, msg); + char *expected = output_to_str(3); + TEST_CHECK(strcmp("(0, 1)", expected) == 0); + TEST_MSG("Wrong result. Mask: %zx. Actual: %s", mask, expected); + LAGraph_Free((void **)&expected, msg); + + free_workspace(); + } + + teardown(); +#endif +} + +void test_CFL_reachability_cycle(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 15; mask == 15; mask++) { + GrB_Info retval; + + init_grammar_aS(); + init_graph_one_cycle(); + init_outputs(); + + OK(run_algorithm(mask)); + TEST_MSG("MASK: %x, ERROR: %s\n", mask, msg); + check_result("(0, 0) (0, 1) (0, 2) (1, 0) (1, 1) (1, 2) (2, 0) (2, 1) (2, 2)"); + + free_workspace(); + } + + teardown(); +#endif +} + +void test_CFL_reachability_two_cycle(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + GrB_Info retval; + + init_grammar_aSb(); + init_graph_double_cycle(); + init_outputs(); + + OK(run_algorithm(mask)); + check_result("(0, 0) (0, 3) (1, 0) (1, 3) (2, 0) (2, 3)"); + + free_workspace(); + } + + teardown(); +#endif +} + +void test_CFL_reachability_indexed_rules(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + GrB_Info retval; + + init_indexed_grammar(); + init_graph_5(); + init_outputs(); + + OK(run_algorithm(mask)); + // Terminals a, b_0 and b_1 are not nonterminals, so they have no outputs + for (size_t i = 0; i < 3; i++) { + TEST_CHECK(outputs[i] == NULL); + TEST_MSG("Mask: %zx, symbol: %zu", mask, i); + } + + char *expected = output_to_str(3); + TEST_CHECK(strcmp("", expected) == 0); + TEST_MSG("Wrong result. Mask: %zx. Actual: %s", mask, expected); + LAGraph_Free((void **)&expected, msg); + + free_workspace(); + } + + teardown(); +#endif +} + +void test_CFL_reachability_labels_more_than_nonterms(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + GrB_Info retval; + + init_grammar_aSb(); + grammar.terms_count = 3; + init_graph_2(); + init_outputs(); + + OK(run_algorithm(mask)); + check_result("(0, 1)"); + + free_workspace(); + } + + teardown(); +#endif +} + +void test_CFL_reachability_complex_grammar(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + GrB_Info retval; + + init_grammar_complex(); + init_graph_1(); + init_outputs(); + + OK(run_algorithm(mask)); + check_result("(0, 7) (1, 6)"); + + free_workspace(); + } + + teardown(); +#endif +} + +void test_CFL_reachability_tree(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + GrB_Info retval; + + init_grammar_aSb(); + init_graph_tree(); + init_outputs(); + + OK(run_algorithm(mask)); + 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(); +#endif +} + +void test_CFL_reachability_line(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + GrB_Info retval; + + init_grammar_aSb(); + init_graph_line(); + init_outputs(); + + OK(run_algorithm(mask)); + check_result("(0, 4) (1, 3)"); + + free_workspace(); + } + teardown(); +#endif +} + +void test_CFL_reachability_two_nodes_cycle(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + GrB_Info retval; + + init_grammar_aSb(); + init_graph_3(); + init_outputs(); + + OK(run_algorithm(mask)); + check_result("(0, 0) (1, 0)"); + + free_workspace(); + } + + teardown(); +#endif +} + +void test_CFL_reachability_with_empty_adj_matrix(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + GrB_Info retval; + + init_grammar_aS(); + grammar.terms_count = 2; + init_graph_4(); + init_outputs(); + + OK(run_algorithm(mask)); + check_result("(0, 0) (1, 1)"); + + free_workspace(); + } + teardown(); +#endif +} + +// Only nonterminals (symbols on the left side of a rule) have outputs. Outputs of +// terminals and symbols that appear in no rule must be NULL. +void check_outputs(GrB_Matrix *expected, size_t mask) { + for (size_t i = 0; i < n_adj_matrices; i++) { + if (expected[i] == NULL) { + TEST_CHECK(outputs[i] == NULL); + TEST_MSG("Mask: %zx, symbol: %zu", mask, i); + continue; + } + + bool equal = false; + OK(LAGraph_Matrix_IsEqual(&equal, outputs[i], expected[i], msg)); + TEST_CHECK(equal); + TEST_MSG("Mask: %zx, symbol: %zu", mask, i); + } +} + +void init_adj_matrices(size_t count, GrB_Index n) { + n_adj_matrices = count; + OK(LAGraph_Calloc((void **)&adj_matrices, n_adj_matrices, sizeof(GrB_Matrix), msg)); + for (size_t i = 0; i < n_adj_matrices; i++) { + OK(GrB_Matrix_new(&adj_matrices[i], GrB_BOOL, n, n)); + } +} + +void test_CFL_reachability_unused_symbols(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + // Symbols: [0 S] [1 x] [2 a] [3 y] [4 b] [5 T] [6 z] [7 C] + // x, y, z appear in no rule; x and z have edges, y is empty. + // C appears only on the left side of a rule. + // S -> a b | a T, T -> S b, C -> b a + OK(LAGraph_Calloc((void **)&grammar.rules, 4, sizeof(LAGraph_rule_EWCNF), msg)); + grammar.nonterms_count = 3; + grammar.terms_count = 5; + grammar.rules_count = 4; + grammar.rules[0] = (LAGraph_rule_EWCNF){0, 2, 4, 0, 0}; + grammar.rules[1] = (LAGraph_rule_EWCNF){0, 2, 5, 0, 0}; + grammar.rules[2] = (LAGraph_rule_EWCNF){5, 0, 4, 0, 0}; + grammar.rules[3] = (LAGraph_rule_EWCNF){7, 4, 2, 0, 0}; + + init_adj_matrices(8, 4); + OK(GrB_Matrix_setElement_BOOL(adj_matrices[1], true, 0, 3)); + OK(GrB_Matrix_setElement_BOOL(adj_matrices[1], true, 3, 0)); + OK(GrB_Matrix_setElement_BOOL(adj_matrices[2], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(adj_matrices[2], true, 1, 2)); + OK(GrB_Matrix_setElement_BOOL(adj_matrices[4], true, 2, 3)); + OK(GrB_Matrix_setElement_BOOL(adj_matrices[4], true, 3, 0)); + OK(GrB_Matrix_setElement_BOOL(adj_matrices[6], true, 2, 2)); + init_outputs(); + + GrB_Matrix expected[8] = {NULL}; + OK(GrB_Matrix_new(&expected[0], GrB_BOOL, 4, 4)); + OK(GrB_Matrix_setElement_BOOL(expected[0], true, 1, 3)); + OK(GrB_Matrix_setElement_BOOL(expected[0], true, 0, 0)); + OK(GrB_Matrix_new(&expected[5], GrB_BOOL, 4, 4)); + OK(GrB_Matrix_setElement_BOOL(expected[5], true, 1, 0)); + OK(GrB_Matrix_new(&expected[7], GrB_BOOL, 4, 4)); + OK(GrB_Matrix_setElement_BOOL(expected[7], true, 3, 1)); + + GrB_Info info = run_algorithm(mask); + OK(info); + TEST_MSG("Mask: %zx, error: %s", mask, msg); + if (info == GrB_SUCCESS) { + check_outputs(expected, mask); + } + + for (size_t i = 0; i < 8; i++) { + GrB_free(&expected[i]); + } + free_workspace(); + } + + teardown(); +#endif +} + +void test_CFL_reachability_unused_symbols_indexed(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + // Symbols: [0 a] [1 b_0] [2 b_1] [3 S_0] [4 S_1] [5 u_0] [6 u_1] [7 u_2] [8 v] + // u_0..u_2 and v appear in no rule; u_1 and v are empty. + // S_i -> a b_i + OK(LAGraph_Calloc((void **)&grammar.rules, 1, sizeof(LAGraph_rule_EWCNF), msg)); + grammar.nonterms_count = 2; + grammar.terms_count = 7; + grammar.rules_count = 1; + grammar.rules[0] = (LAGraph_rule_EWCNF){ + 3, 0, 1, 2, LAGraph_EWNCF_INDEX_NONTERM | LAGraph_EWNCF_INDEX_PROD_B}; + + init_adj_matrices(9, 3); + OK(GrB_Matrix_setElement_BOOL(adj_matrices[0], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(adj_matrices[1], true, 1, 2)); + OK(GrB_Matrix_setElement_BOOL(adj_matrices[2], true, 1, 0)); + OK(GrB_Matrix_setElement_BOOL(adj_matrices[5], true, 0, 0)); + OK(GrB_Matrix_setElement_BOOL(adj_matrices[7], true, 2, 1)); + init_outputs(); + + GrB_Matrix expected[9] = {NULL}; + OK(GrB_Matrix_new(&expected[3], GrB_BOOL, 3, 3)); + OK(GrB_Matrix_setElement_BOOL(expected[3], true, 0, 2)); + OK(GrB_Matrix_new(&expected[4], GrB_BOOL, 3, 3)); + OK(GrB_Matrix_setElement_BOOL(expected[4], true, 0, 0)); + + GrB_Info info = run_algorithm(mask); + OK(info); + TEST_MSG("Mask: %zx, error: %s", mask, msg); + if (info == GrB_SUCCESS) { + check_outputs(expected, mask); + } + + for (size_t i = 0; i < 9; i++) { + GrB_free(&expected[i]); + } + free_workspace(); + } + + teardown(); +#endif +} + +//==================== +// Tests with invalid result +//==================== + +void test_CFL_reachability_invalid_rules(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + GrB_Info retval; + + init_grammar_aSb(); + init_graph_double_cycle(); + init_outputs(); + + const LAGraph_rule_EWCNF invalid_rules[] = { + {0, -1, 1, 0, 0}, // Missing first RHS symbol. + {-1, 1, 2, 0, 0}, + {-2, 1, 2, 0, 0}, + {6, 1, 2, 0, 0}, // Symbol indices must be less than six. + {0, -2, 2, 0, 0}, + {0, 6, 2, 0, 0}, + {0, 1, -2, 0, 0}, + {0, 1, 6, 0, 0}, + {0, 6, -1, 0, 0}, + // Valid base indices, but indexed groups extend beyond the symbol array. + {5, 1, 2, 2, LAGraph_EWNCF_INDEX_NONTERM}, + {0, 5, 2, 2, LAGraph_EWNCF_INDEX_PROD_A}, + {0, 1, 5, 2, LAGraph_EWNCF_INDEX_PROD_B}, + {0, 1, 2, UINT32_MAX, LAGraph_EWNCF_INDEX_PROD_A}, + }; + + for (size_t mask = 0; mask < 16; mask++) { + for (size_t i = 0; i < sizeof(invalid_rules) / sizeof(invalid_rules[0]); i++) { + grammar.rules[0] = invalid_rules[i]; + check_error(GrB_INVALID_VALUE, mask); + TEST_CHECK(strstr(msg, "Rule with index 0") != NULL); + TEST_MSG("Mask: %zu, invalid rule: %zu", mask, i); + } + } + + free_workspace(); + teardown(); +#endif +} + +// Fail one LAGraph allocation while the brutal allocator tracks both LAGraph +// and GraphBLAS resources. GraphBLAS allocations remain enabled for cleanup. +static size_t calloc_count, calloc_fail_at; + +static void *fail_calloc(size_t n, size_t size) { + if (++calloc_count == calloc_fail_at) { + return NULL; + } + return LG_brutal_calloc(n, size); +} + +void test_CFL_reachability_allocation_failure(void) { +#if LAGRAPH_SUITESPARSE + OK(LG_brutal_setup(msg)); + + GrB_Matrix adj[3] = {NULL}, result[3] = {NULL}; + for (size_t i = 0; i < 3; i++) { + OK(GrB_Matrix_new(&adj[i], GrB_BOOL, 3, 3)); + } + OK(GrB_Matrix_setElement_BOOL(adj[1], true, 0, 1)); + OK(GrB_Matrix_setElement_BOOL(adj[2], true, 1, 2)); + for (size_t i = 0; i < 3; i++) { + OK(GrB_Matrix_wait(adj[i], GrB_MATERIALIZE)); + } + + LAGraph_rule_EWCNF rule = {0, 1, 2, 0, 0}; + // Warm up GraphBLAS caches before taking the allocation baseline. + OK(LAGraph_CFL_reachability_adv(result, adj, 3, &rule, 1, msg, 0)); + for (size_t i = 0; i < 3; i++) { + OK(GrB_free(&result[i])); + } + // First three workspace arrays, symbol/rule setup, then a temporary matrix + // inside the first fixed-point iteration (allocation 16, through TRY_I). + const size_t failures[] = {1, 2, 3, 4, 5, 6, 16, 0}; + int64_t baseline = LG_nmalloc; + for (size_t i = 0; i < sizeof(failures) / sizeof(failures[0]); i++) { + calloc_count = 0; + calloc_fail_at = failures[i]; + LAGraph_Calloc_function = fail_calloc; + GrB_Info info = LAGraph_CFL_reachability_adv(result, adj, 3, &rule, 1, msg, 0); + LAGraph_Calloc_function = LG_brutal_calloc; + TEST_CHECK(info == (calloc_fail_at == 0 ? GrB_SUCCESS : GrB_OUT_OF_MEMORY)); + TEST_MSG("Failed allocation: %zu, retval: %d", calloc_fail_at, info); + for (size_t j = 0; j < 3; j++) { + OK(GrB_free(&result[j])); + } + TEST_CHECK(LG_nmalloc == baseline); + TEST_MSG("Failed allocation: %zu, live allocations: %" PRId64 + ", expected: %" PRId64, calloc_fail_at, LG_nmalloc, baseline); + } + + for (size_t i = 0; i < 3; i++) { + OK(GrB_free(&adj[i])); + } + OK(LG_brutal_teardown(msg)); +#endif +} + +void test_CFL_reachability_null_msg(void) { +#if LAGRAPH_SUITESPARSE + setup(); + + for (size_t mask = 0; mask < 16; mask++) { + init_grammar_aSb(); + init_graph_double_cycle(); + init_outputs(); + + OK(LAGraph_CFL_reachability_adv(outputs, adj_matrices, + grammar.terms_count + grammar.nonterms_count, + grammar.rules, grammar.rules_count, NULL, mask)); + free_workspace(); + + init_grammar_aSb(); + init_graph_double_cycle(); + init_outputs(); + GrB_free(&adj_matrices[0]); + GrB_free(&adj_matrices[1]); + + GrB_Info info = LAGraph_CFL_reachability_adv( + outputs, adj_matrices, grammar.terms_count + grammar.nonterms_count, + grammar.rules, grammar.rules_count, NULL, mask); + TEST_CHECK(info == GrB_NULL_POINTER); + TEST_MSG("Mask: %zu, retval: %d", mask, info); + free_workspace(); + } + + teardown(); +#endif +} + +void test_CFL_reachability_null_pointers(void) { +#if LAGRAPH_SUITESPARSE + + setup(); + + GrB_Info retval; + + init_grammar_aSb(); + init_graph_double_cycle(); + init_outputs(); + + // adj_matrices[0] = NULL; + // adj_matrices[1] = NULL; + GrB_free(&adj_matrices[0]); + GrB_free(&adj_matrices[1]); + + check_error(GrB_NULL_POINTER, 0); + TEST_CHECK(strstr(msg, "Adjacency matrices with these indexes are null: 0, 1") + != NULL); + + // adj_matrices = NULL; + for (size_t i = 0; i < n_adj_matrices; i++) { + GrB_free(&adj_matrices[i]); + } + + LAGraph_Free((void **)&adj_matrices, msg); + check_error(GrB_NULL_POINTER, 0); + TEST_CHECK(strstr(msg, "The adjacency matrices array cannot be null.") != NULL); + + free_workspace(); + init_grammar_aSb(); + init_graph_double_cycle(); + init_outputs(); + + // outputs = NULL; + LAGraph_Free((void **)&outputs, msg); + check_error(GrB_NULL_POINTER, 0); + TEST_CHECK(strstr(msg, "The outputs array cannot be null.") != NULL); + + free_workspace(); + init_grammar_aSb(); + init_graph_double_cycle(); + init_outputs(); + + // grammar.rules = NULL; + LAGraph_Free((void **)&grammar.rules, msg); + check_error(GrB_NULL_POINTER, 0); + TEST_CHECK(strstr(msg, "The rules array cannot be null.") != NULL); + + free_workspace(); + teardown(); +#endif +} + +TEST_LIST = { + {"CFG_reachability_indexed", test_CFL_indexed}, + {"CFG_reachability_indexed_2", test_CFL_indexed_2}, + {"CFG_reachability_indexed_nonterm_only", test_CFL_indexed_nonterm_only}, + {"CFL_reachability_complex_grammar", test_CFL_reachability_complex_grammar}, + {"CFG_reachability_indexed_simple", test_CFL_indexed_simple}, + {"CFG_reachability_indexed_simple_exploded", test_CFL_indexed_simple_exploded}, + {"CFL_reachability_cycle", test_CFL_reachability_cycle}, + {"CFG_reachability_indexed_rules", test_CFL_reachability_indexed_rules}, + {"CFL_reachability_two_cycle", test_CFL_reachability_two_cycle}, + {"CFL_reachability_labels_more_than_nonterms", + test_CFL_reachability_labels_more_than_nonterms}, + {"CFL_reachability_tree", test_CFL_reachability_tree}, + {"CFL_reachability_line", test_CFL_reachability_line}, + {"CFL_reachability_two_nodes_cycle", test_CFL_reachability_two_nodes_cycle}, + {"CFG_reach_basic_invalid_rules", test_CFL_reachability_invalid_rules}, + {"test_CFL_reachability_with_empty_adj_matrix", + test_CFL_reachability_with_empty_adj_matrix}, +#if !defined(GRAPHBLAS_HAS_CUDA) + {"CFL_reachability_unused_symbols", test_CFL_reachability_unused_symbols}, + {"CFL_reachability_unused_symbols_indexed", + test_CFL_reachability_unused_symbols_indexed}, + {"CFG_reachability_allocation_failure", test_CFL_reachability_allocation_failure}, + {"CFG_reachability_null_msg", test_CFL_reachability_null_msg}, + {"CFG_reachability_null_pointers", test_CFL_reachability_null_pointers}, +#endif + {NULL, NULL}}; diff --git a/include/LAGraphX.h b/include/LAGraphX.h index 3355addb2e..6b7c14704e 100644 --- a/include/LAGraphX.h +++ b/include/LAGraphX.h @@ -1028,6 +1028,95 @@ int LAGraph_SquareClustering int32_t index; // For rules that can be grouped by index } LAGraph_rule_WCNF; +// Flags for indecies of EWCNF rules +enum { + LAGraph_EWNCF_INDEX_NONTERM = 1 << 0, + LAGraph_EWNCF_INDEX_PROD_A = 1 << 1, + LAGraph_EWNCF_INDEX_PROD_B = 1 << 2, +}; + +// Production rule of Context-free grammar in Extended Weak Chomsky Normal Form +// +// All grammar symbols (terminal and nonterminals) are presented by natural numbers from range [0; symbols_amount). +// +// Rule without indecies defined by tuple of [NONTERM, PROD_A, PROD_B, INDEXED_COUNT, INDEXED] +// in Extended Weak Chomsky Normal Form where INDEXED_COUNT and INDEXED are zero: +// Variable -> eps: [NONTERM, -1, -1, 0, 0] +// Variable -> X: [NONTERM, PROD_X, -1, 0, 0] +// Variable -> X Y: [NONTERM, PROD_X, PROD_Y, 0, 0] +// +// Where X, Y are arbitrary grammar symbols (terminal or nonterminal) +// +// Example: +// Nonterms: [0 S] [1 A] [2 B] [3 C] +// Terms: [4 a] [5 b] +// S -> A B [0 1 2 0 0] +// S -> A b [0 1 5 0 0] +// S -> a B [0 4 2 0 0] +// S -> C [0 3 -1 0 0] +// A -> a [1 4 -1 0 0] +// B -> b [2 5 -1 0 0] +// C -> b [3 5 -1 0 0] +// S -> eps [0 -1 -1 0 0] +// +// Warning: +// Variable -> _ B: [NONTERM, -1, PROD_B, INDEXED_COUNT, INDEXED] is not valid rule and may causes errors +// +// --------------- Indexed rules --------------- +// Rule may represent a family of indexed productions of size `indexed_count`. +// +// Indexed symbols are specified by 'indexed' bit mask +// Posible indexed forms: +// +// 1) N -> A_i B [N, A, B, INDEXED_COUNT, INDEX_PROD_A] +// 2) N -> A B_i [N, A, B, INDEXED_COUNT, INDEX_PROD_A] +// 3) N -> A_i B_i [N, A, B, INDEXED_COUNT, INDEX_PROD_A | INDEX_PROD_B] +// 4) N_i -> A B [N, A, B, INDEXED_COUNT, INDEX_NONTERM] +// 5) N_i -> A_i B [N, A, B, INDEXED_COUNT, INDEX_NONTERM | INDEX_PROD_A] +// 6) N_i -> A B_i [N, A, B, INDEXED_COUNT, INDEX_NONTERM | INDEX_PROD_B] +// 7) N_i -> A_i B_i [N, A, B, INDEXED_COUNT, INDEX_NONTERM | INDEX_PROD_A | INDEX_PROD_B] +// +// There are constraints about indexed rules: +// +// 1) If symbol A is indexed in one rule, it must be indexed +// in all rules where it appears. +// 2) If multiple symbols are indexed in the same rule, +// their indexed_count must be equal. +// 3) Indexed symbols must be numbered consecutively: +// A_0, A_1, ..., A_{k-1} -> [base, base+1, ..., base+k-1] +// 4) For indexed rules, base symbol A must represent by A_0. +// +// Actual rules are expanded as: +// +// for i in range [0; indexed_count) +// N[_i] -> A[_i] B[_i] (only for symbols marked as indexed) +// +// Example: +// Z -> X_i Y (represented by [Z; X, Y, 3, INDEX_PROD_A]) expanded as: +// +// Z -> X_0 Y +// Z -> X_1 Y +// Z -> X_2 Y +// +// --------------- Difference from standard Weak Chomsky Normal Form (WCNF) --------------- +// +// Standard WCNF allows only: +// N -> A B (A,B are nonterminals) +// N -> a (a are terminal) +// N -> eps +// +// EWCNF allows: +// N -> X Y (X are terminal or nonterminal, Y are terminal or nonterminal) +// N -> X (X are terminal or nonterminal) +// N -> eps +// and indexed rules. + typedef struct { + int32_t nonterm; // LHS nonterminal symbol id + int32_t prod_A; // first RHS symbol id or -1 + int32_t prod_B; // second RHS symbol id or -1 + uint32_t indexed_count; // Number of indexed rules (zero if not indexed) + uint8_t indexed; // Bitmask of indexed symbols +} LAGraph_rule_EWCNF; // LAGraph_CFL_reachability: Context-Free Language Reachability Matrix-Based Algorithm // @@ -1096,6 +1185,95 @@ GrB_Info LAGraph_CFL_reachability char *msg // Message string for error reporting. ) ; + +// LAGraph_CFL_reachability_adv: 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 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`). +// +// Terminals and non-terminals are enumerated by integers starting from zero. +// The start non-terminal is the non-terminal with index 0. +// +// Example: +// +// Graph: +// ┌───┐ ┌───┐ ┌───┐ ┌───┐ ┌───┐ +// │ 0 ├───► 1 ├───► 2 ├───► 3 ├───► 4 │ +// └───┘ a └─┬─┘ a └─▲─┘ b └───┘ b └───┘ +// │ │ +// │ ┌───┐ │ +// a└─► 5 ├─┘b +// └───┘ +// +// Grammar: S -> aSb | ab +// +// There are paths from node [1] to node [3] and from node [1] to node [2] that form +// the word "ab" ([1]-a->[2]-b->[3] and [1]-a->[5]-b->[2]). The word "ab" is in the +// language generated by our context-free grammar, so the pairs (1, 3) and (1, 2) will +// be included in the result. +// +// Note: It doesn't matter how many paths exist from node [A] to node [B] that form a +// word in the language. If at least one path exists, the pair ([A], [B]) will be +// included in the result. +// +// In contrast, the path from node [1] to node [4] forms the word "abb" +// ([1]-a->[2]-b->[3]-b->[4]) and the word "abbb" ([1]-a->[5]-b->[2]-b->[3]-b->[4]). +// The words "aab" and "abbb" are not in the language, so the pair (1, 4) will not be +// included in the result. +// +// With this graph and grammar, we obtain the following results: +// (0, 4) - because there exists a path (0-1-2-3-4) that forms the word "aabb" +// (1, 3) - because there exists a path (1-2-3) that forms "ab" +// (1, 2) - because there exists a path (1-5-2) that forms the word "ab" +// (0, 3) - because there exists a path (0-1-5-2-3) that forms the word "aabb" +GrB_Info LAGraph_CFL_reachability_adv +( + // Output + GrB_Matrix *outputs, // Array of matrices containing results. + // The size of the array must be equal to the count + // of symbols (symbols_amount). + // + // Each matrix is square, with size equal to the number of + // vertices in the graph. Matrices are allocated by this + // method, the caller must free them. + // + // outputs[k]: (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 nonterminal 'k' of the specified CFG. + // + // Note: outputs[t] is NULL if the symbol 't' is not a + // nonterminal, i.e. it is not on the left side of any rule + // (terminals and symbols unused by rules). + // Input + const GrB_Matrix *adj_matrices, // Array of adjacency matrices representing the graph. + // The length of this array is equal to the count of + // symbols (symbols_amount). + // + // Each matrix is square, with size equal to the + // number of vertices in the graph. + // + // adj_matrices[i]: (i, j) == 1 if and only if there + // is an edge between nodes i and j with the label of + // the symbol corresponding to index 'i' (where i is + // in the range [0, symbols_amount - 1]). + // + // Note: adj_matrices[N], where N is the index of a + // nonterminal, must be initialized as an empty square + // matrix of size symbols_amount. + size_t symbols_amount, // Count of terminal and nonterminals + const LAGraph_rule_EWCNF *rules, // The rules of the CFG. + // Warning: now not ready for N -> A rules, where + // is N and A are nonterminals. + size_t rules_count, // The total number of rules in the CFG. + char *msg, // Message string for error reporting. + int8_t optimizations // Optimizations flags +) ; + + + //------------------------------------------------------------------------------ // a simple example of an algorithm //------------------------------------------------------------------------------