update: union, reduce and insert

and LAGraph_CFL_AllPaths_free_outputs

Author: DanyaLitva

experimental/algorithm/LAGraph_CFL_AllPaths.c

@@ -16,18 +16,12 @@
 #include <LAGraphX.h>
 
 //Merging two ordered arrays of internal vertices in the add function
-static GrB_Index* merge_all_paths(GrB_Index* n, const void* left, const GrB_Index na, const void* right, const GrB_Index nb){
-  GrB_Index* a = (GrB_Index*) left;
-  GrB_Index* b = (GrB_Index*) right;
-  
-  if (na == 0 && nb == 0) {
-    *n = 0;
-    return NULL;
-  }
-  
+static GrB_Index* merge_all_paths(GrB_Index* n, const GrB_Index* a, const size_t na, const GrB_Index* b, const size_t nb)
+{  
   GrB_Index *tmp = malloc((na + nb) * sizeof(GrB_Index));
 
-  GrB_Index ia = 0, ib = 0, outn = 0;
+  size_t ia = 0, ib = 0, outn = 0;
+  // Handle the first elements of both arrays to initialize tmp and avoid checking if tmp is empty in the loop
   if (na > 0 && nb > 0) {
     if (a[0] < b[0]) {
       tmp[outn++] = a[ia++];
@@ -73,49 +67,92 @@ static GrB_Index* merge_all_paths(GrB_Index* n, const void* left, const GrB_Inde
   return tmp;
 }
 
-void clear_elem_all_paths(AllPathsElem *z){
-  if(z->middle){
-    free(z->middle);
-    z->middle=NULL;
+static GrB_Index* insert_all_paths(size_t* n, const GrB_Index* arr, size_t len, GrB_Index value)
+{
+  // Array is ordered, so we can use binary search to find the position to insert value
+  size_t l = 0, r = len, m = 0;
+  while (l < r) {
+    m = l + (r - l) / 2;
+    if (arr[m] < value) l = m + 1;
+    else r = m;
   }
-  z->n = 0;
+
+    // If value is already in the array, return the original array
+    if (l < len && arr[l] == value) {
+        GrB_Index *tmp = malloc(len * sizeof(GrB_Index));
+        memcpy(tmp, arr, len * sizeof(GrB_Index));
+        *n = len;
+        return tmp;
+    }
+
+    GrB_Index *tmp = malloc((len + 1) * sizeof(GrB_Index));
+    memcpy(tmp, arr, l * sizeof(GrB_Index));
+    tmp[l] = value;
+    memcpy(tmp + l + 1, arr + l, (len - l) * sizeof(GrB_Index));
+
+    *n = len + 1;
+    return tmp;
 }
 
 void add_all_paths(AllPathsElem *z, AllPathsElem *x, AllPathsElem *y)
 {
-  //temp is needed to avoid freeing the memory of z in case z == x or z == y
-  AllPathsElem temp;
-  temp.middle = merge_all_paths(&temp.n, x->middle, x->n, y->middle, y->n);
-  clear_elem_all_paths(x);
-  clear_elem_all_paths(y);
-  z->middle = temp.middle;
-  z->n = temp.n;
+    // temp is needed to avoid freeing the memory of z in case z == x or z == y
+    AllPathsElem temp;
+
+    // x->n and y->n are not equal to zero
+    if (x->n == 1 && y->n == 1) {
+        if (x->data.single_elem == y->data.single_elem) {
+            temp.n = 1;
+            temp.data.single_elem = x->data.single_elem;
+        } else {
+            temp.n = 2;
+            temp.data.middle = malloc(2 * sizeof(GrB_Index));
+            if (x->data.single_elem < y->data.single_elem) {
+                temp.data.middle[0] = x->data.single_elem;
+                temp.data.middle[1] = y->data.single_elem;
+            } else {
+                temp.data.middle[0] = y->data.single_elem;
+                temp.data.middle[1] = x->data.single_elem;
+            }
+        }
+    } else if (x->n == 1) {
+        temp.data.middle = insert_all_paths(&temp.n, y->data.middle, y->n, x->data.single_elem);
+    } else if (y->n == 1) {
+        temp.data.middle = insert_all_paths(&temp.n, x->data.middle, x->n, y->data.single_elem);
+    } else {
+        temp.data.middle = merge_all_paths(&temp.n, x->data.middle, x->n, y->data.middle, y->n);
+    }
+
+    if(x->n > 1){
+      free(x->data.middle);
+    }
+    if(y->n > 1){
+      free(y->data.middle);
+    }
+
+    *z = temp;
 }
 
 void mult_all_paths(AllPathsElem *z,
                      const AllPathsElem *x, GrB_Index ix, GrB_Index jx,
                      const AllPathsElem *y, GrB_Index iy, GrB_Index jy,
                      const void *theta)
 {
-  z->middle = malloc(sizeof(GrB_Index));
-  z->middle[0] = jx;
+  z->data.single_elem = jx;
   z->n = 1;
 }
 
 void set_all_paths(AllPathsElem *z, const AllPathsElem *x, const bool *edge_exist)
 {
-  AllPathsElem temp;
-  temp.middle = NULL;
-  temp.n = 0;
-  
-  if (edge_exist && *edge_exist){
-    temp.middle = malloc(sizeof(GrB_Index));
-    temp.n = 1;
-    temp.middle[0] = GrB_INDEX_MAX; //GrB_INDEX_MAX - marker of A->eps and A->t
+  if (*edge_exist)
+  {
+    z->data.single_elem = GrB_INDEX_MAX; // A special value to indicate that this path corresponds to a terminal rule (A->t) or an epsilon rule (A->eps)
+    z->n = 1;
+  }
+  else{
+    z->n = x->n;
+    z->data.single_elem = x->data.single_elem;
   }
-  
-  z->middle = merge_all_paths(&z->n, x->middle, x->n, temp.middle, temp.n);
-  clear_elem_all_paths(&temp);
 }
 
 #define MULT_PATH_INDEX_DEFN                                                   \
@@ -124,26 +161,25 @@ void set_all_paths(AllPathsElem *z, const AllPathsElem *x, const bool *edge_exis
 "                     const AllPathsElem *y, GrB_Index iy, GrB_Index jy, \n"      \
 "                     const void *theta) \n"      \
 "{ \n"      \
-"  z->middle = malloc(sizeof(GrB_Index)); \n"      \
-"  z->middle[0] = jx; \n"      \
+"  z->data.single_elem = jx; \n"      \
 "  z->n = 1; \n"      \
 "}"
 
 //Adding the count of all internal vertices in a reduction
 void add_get_nvals_all_paths(AllPathsElem *z, const AllPathsElem *x, const AllPathsElem *y)
 {
   z->n = x->n + y->n;
-  z->middle = NULL;
 }
 
 //Made global so that the get_nvals_all_paths matches the GrB_Matrix_nvals signature
-static GrB_Type* AllPaths_type_get_nvals = NULL;
+static GrB_Type AllPaths_type = NULL;
 static GrB_Monoid AllPaths_monoid_get_nvals = NULL;
 
 //A function that replaces GrB_Matrix_nvals in Reachability to check if new vertices have been added to the matrix.
-GrB_Info get_nvals_all_paths(GrB_Index *nvals, const GrB_Matrix A){
+static GrB_Info get_nvals_all_paths(GrB_Index *nvals, const GrB_Matrix A)
+{
   GrB_Scalar s = NULL;
-  GrB_Scalar_new(&s, *AllPaths_type_get_nvals);
+  GrB_Scalar_new(&s, AllPaths_type);
   GrB_Info info = GrB_reduce(s, NULL, AllPaths_monoid_get_nvals, A, NULL);
 
   if (info != GrB_SUCCESS)
@@ -152,50 +188,23 @@ GrB_Info get_nvals_all_paths(GrB_Index *nvals, const GrB_Matrix A){
     return info;
   }
 
-  AllPathsElem result = {0, NULL};
-  GrB_Info extract_info = GrB_Scalar_extractElement_UDT(&result, s);
-    
-  if (extract_info == GrB_NO_VALUE)
+  AllPathsElem result = {0};
+  info = GrB_Scalar_extractElement_UDT(&result, s);
+  if (info == GrB_NO_VALUE)
   {
       result.n = 0;
   }
-  else if (extract_info != GrB_SUCCESS)
+  else if (info != GrB_SUCCESS)
   {
       GrB_free(&s);
-      return extract_info;
+      return info;
   }
 
   *nvals = result.n;
   GrB_free(&s);
   return GrB_SUCCESS;
 }
 
-//To test the non-reduction approach in the future
-GrB_Info get_nvals_all_paths2(GrB_Index *nvals, const GrB_Matrix A){
-  GrB_Index accum = 0;
-  GxB_Iterator iterator;
-  GxB_Iterator_new(&iterator);
-  GrB_Info info = GxB_Matrix_Iterator_attach(iterator, A, NULL);
-  info = GxB_Matrix_Iterator_seek(iterator, 0);
-  AllPathsElem val;
-  
-  while (info != GxB_EXHAUSTED)
-  {
-    GxB_Iterator_get_UDT(iterator, (void*) &val);
-    accum+=val.n;
-    info = GxB_Matrix_Iterator_next(iterator);
-  }
-  
-  GrB_free(&iterator);
-  *nvals = accum;
-  
-  return GrB_SUCCESS;
-}
-
-// all_paths_ptr_t is a pointer to the type of elements of the outputs matrices.
-// Use GrB_free(all_paths_ptr_t) after you finish working with the outputs matrices.
-// Important: Do not free all_paths_ptr_t until all work with the output matrices is complete.
-// Accessing matrices after freeing their type is undefined behavior.
 GrB_Info LAGraph_CFL_AllPaths(
     // Output
     GrB_Matrix *outputs, // Array of matrices containing results.
@@ -207,6 +216,8 @@ GrB_Info LAGraph_CFL_AllPaths(
                          // There are paths from I to M by nonterminal N1 and from M to J by nonterminal N2,
                          // and A->N1 N2 where outputs[k] corresponds to nonterminal A.
                          // GrB_INDEX_MAX in the array is a special value for A->eps and A->t.
+    // AllPaths type - elements of the output matrices.
+    GrB_Type *all_paths_ptr_t,      // Pass a pointer to GrB_Type.
     // Input
     const GrB_Matrix *adj_matrices, // Array of adjacency matrices representing the graph.
                                     // The length of this array is equal to the count of
@@ -216,9 +227,6 @@ GrB_Info LAGraph_CFL_AllPaths(
                                     // is an edge between nodes i and j with the label of
                                     // the terminal corresponding to index 't' (where t is
                                     // in the range [0, terms_count - 1]).
-    GrB_Type *all_paths_ptr_t,      // AllPaths type - elements of the output matrices.
-                                    // Pass a pointer to GrB_Type and
-                                    // free it after you finish working with outputs matrices.
     int64_t terms_count,            // The total number of terminal symbols in the CFG.
     int64_t nonterms_count,         // The total number of non-terminal symbols in the CFG.
     const LAGraph_rule_WCNF *rules, // The rules of the CFG.
@@ -230,7 +238,6 @@ GrB_Info LAGraph_CFL_AllPaths(
   return (GrB_NOT_IMPLEMENTED);
 #else
   // Semiring components
-  GrB_Type AllPaths_type = NULL;
   GrB_BinaryOp AllPaths_add = NULL;
   GrB_BinaryOp AllPaths_add_get_nvals = NULL;
   GrB_Monoid AllPaths_monoid = NULL;
@@ -244,12 +251,11 @@ GrB_Info LAGraph_CFL_AllPaths(
   GrB_free(all_paths_ptr_t);
   GRB_TRY(GrB_Type_new(all_paths_ptr_t, sizeof(AllPathsElem)));
   AllPaths_type = *all_paths_ptr_t;
-  AllPaths_type_get_nvals = all_paths_ptr_t;
 
   GRB_TRY(GrB_Scalar_new(&Theta, GrB_BOOL));
   GRB_TRY(GrB_Scalar_setElement_BOOL(Theta, false));
 
-  AllPathsElem bottom = {0, NULL};
+  AllPathsElem bottom = {0};
   GRB_TRY(GrB_Scalar_new(&bottom_scalar, AllPaths_type));
   GRB_TRY(GrB_Scalar_setElement_UDT(bottom_scalar, (void *)(&bottom)));
 
@@ -304,19 +310,62 @@ GrB_Info LAGraph_CFL_AllPaths(
                          AllPaths_add_get_nvals,
                          (void *)(&bottom)));
 
-  CFL_Semiring semiring = {.type = AllPaths_type,
+  CFL_Semiring semiring = {
+      .type = AllPaths_type,
       .semiring = AllPaths_semiring,
       .add = AllPaths_add,
       .mult = AllPaths_mult,
       .init_path = AllPaths_set,
       .bottom_scalar = bottom_scalar,
-    .get_nvals = get_nvals_all_paths};
+      .get_nvals = get_nvals_all_paths};
 
   LG_TRY(LAGraph_CFPQ_core(outputs, adj_matrices, terms_count, nonterms_count, rules, rules_count, &semiring, msg));
 
-  AllPaths_type = NULL;
-  AllPaths_type_get_nvals = NULL;
   LG_FREE_WORK;
   return GrB_SUCCESS;
 #endif
 }
+
+// Helper function to free the output matrix of LAGraph_CFL_AllPaths, which contains elements of type AllPathsElem with dynamically allocated arrays of intermediate vertices. 
+static void free_AllPaths_matrix(GrB_Matrix* ptr_output) 
+{
+  GxB_Iterator iterator;
+  GxB_Iterator_new(&iterator);
+  GrB_Info info = GxB_Matrix_Iterator_attach(iterator, *ptr_output, NULL);
+  info = GxB_Matrix_Iterator_seek(iterator, 0);
+  AllPathsElem val;
+
+  while (info != GxB_EXHAUSTED)
+  {
+    GxB_Iterator_get_UDT(iterator, (void*)&val);
+    if (val.n > 1 && val.data.middle != NULL) {
+      free(val.data.middle);
+    }
+    info = GxB_Matrix_Iterator_next(iterator);
+  }
+
+  GrB_free(&iterator);
+  GrB_free(ptr_output);
+}
+
+// Free outputs and all_paths_ptr_t after you have finished working with the output matrices from LAGraph_CFL_AllPaths.
+GrB_Info LAGraph_CFL_AllPaths_free_outputs(GrB_Matrix* outputs, int64_t nonterms_count, GrB_Type* all_paths_ptr_t)
+{
+#if GxB_IMPLEMENTATION < GxB_VERSION(9, 4, 5)
+  return (GrB_NOT_IMPLEMENTED);
+#else
+  if (outputs) {
+    for (size_t i = 0; i < nonterms_count; i++) {
+      if (outputs[i] == NULL)
+        continue;
+      free_AllPaths_matrix(&outputs[i]);
+      outputs[i] = NULL;
+    }
+    free(outputs);
+  }
+  outputs = NULL;
+  GrB_free(all_paths_ptr_t);
+  return GrB_SUCCESS;
+#endif
+}
+