Merge branch 'v1.3.x' into Louvain

Author: DrTimothyAldenDavis

.gitignore

@@ -43,8 +43,6 @@ CMakeFiles*
 *~
 data/comments*.txt
 data/comments*.mtx
-
-*/build/*
 **/.venv/*
 
 # Do not ignore this file

Config/LAGraph.h.in

@@ -2410,11 +2410,11 @@ int LAGr_SingleSourceShortestPath
 ) ;
 
 //------------------------------------------------------------------------------
-// LAGr_Betweenness: betweeness centrality metric
+// LAGr_Betweenness: betweenness centrality metric
 //------------------------------------------------------------------------------
 
-/** LAGr_Betweenness: betweeness centrality metric.  This methods computes an
- * approximation of the betweeness-centrality metric of all nodes in the graph.
+/** LAGr_Betweenness: betweenness centrality metric.  This methods computes an
+ * approximation of the betweenness-centrality metric of all nodes in the graph.
  * Only a few given source nodes are used for the approximation.  This is an
  * Advanced algorithm (G->AT is required).
  *
@@ -2442,7 +2442,7 @@ LAGRAPH_PUBLIC
 int LAGr_Betweenness
 (
     // output:
-    GrB_Vector *centrality,     // centrality(i): betweeness centrality of i
+    GrB_Vector *centrality,     // centrality(i): betweenness centrality of i
     // input:
     const LAGraph_Graph G,      // input graph
     const GrB_Index *sources,   // source vertices to compute shortest paths

data/test_zero_cap.mtx

@@ -0,0 +1,13 @@
+%%MatrixMarket matrix coordinate real general
+% Example weighted directed graph adjacency (1-indexed)
+% All weights are positive; some edges explicitly have zero weight
+5 5 9
+1 2 0.500000
+1 4 0.000000
+2 3 1.250000
+2 5 2.000000
+3 1 0.000000
+3 4 3.141590
+4 2 0.750000
+5 4 0.125000
+5 5 1.000000

experimental/algorithm/LAGr_BreadthFirstSearch_Extended.c

@@ -27,6 +27,8 @@
 
 #include "LG_alg_internal.h"
 
+// FIXME: make src, max_level, dest, many_expected inputs GrB_Scalar
+
 int LAGr_BreadthFirstSearch_Extended
 (
     // output:

experimental/algorithm/LAGr_EdgeBetweennessCentrality.c

@@ -21,7 +21,7 @@
 // FIXED in v10.2 but may need more work in GraphBLAS.
 
 // LAGr_EdgeBetweennessCentrality: Exact algorithm for computing
-// betweeness centrality.
+// betweenness centrality.
 
 // This is an Advanced algorithm (no self edges allowed)
 
@@ -91,12 +91,11 @@ void LG_EBC_add_one_divide_function (double *z, const double *x, const double *y
 int LAGr_EdgeBetweennessCentrality
 (
     // output:
-    GrB_Matrix *centrality,     // centrality(i): betweeness centrality of i
+    GrB_Matrix *centrality,     // centrality(i): betweenness centrality of i
     // input:
     LAGraph_Graph G,            // input graph
     GrB_Vector sources,         // source vertices to compute shortest paths
                                 // (if NULL or empty, use all vertices)
-
     char *msg
 )
 {

experimental/algorithm/LAGr_Jaccard.c

@@ -0,0 +1,167 @@
+//------------------------------------------------------------------------------
+// LAGraph_Jaccard - parallel jaccard similarity
+//------------------------------------------------------------------------------
+
+// LAGraph, (c) 2019-2025 by The LAGraph Contributors, All Rights Reserved.
+// SPDX-License-Identifier: BSD-2-Clause
+//
+// For additional details (including references to third party source code and
+// other files) see the LICENSE file or contact permission@sei.cmu.edu. See
+// Contributors.txt for a full list of contributors. Created, in part, with
+// funding and support from the U.S. Government (see Acknowledgments.txt file).
+// DM22-0790
+
+// Contributed by Elaheh Hassani and Tim Davis, Texas A&M University
+
+//------------------------------------------------------------------------------
+// LAGr_Jaccard: compute Jaccard similarity (weight) coefficients for an undirected graph.
+//
+// Inputs:
+// G             - a valid LAGraph_Graph with:
+//                    * G->A structurally symmetric (undirected)
+//                    * no self-edges
+//                    * G->out_degree cached (no explicit zeros)
+// all_pairs     - if true, compute Jaccard for all pairs of nodes in the graph;
+//                 if false, compute only for neighboring nodes (current edges).
+//
+// Output:
+// JC            - an n-by-n matrix of Jaccard coefficients. JC is always
+//                 returned as an upper-triangular matrix. When all_pairs is
+//                 true, the diagonal entries are equal to 1. When all_pairs is
+//                 false, the diagonal is structurally empty.
+//
+// Jaccard similarity measures the overlap of neighbor sets.
+// For a pair of nodes i and j in the graph G:
+//         JC(i,j) = |N(i) ∩ N(j)| / |N(i) ∪ N(j)|
+// where N(i) is the set of neighbors of i in G.
+//
+// The computation proceeds in two stages.  First, an “intersection matrix” B
+// is formed where B(i,j) = |N(i) ∩ N(j)|.  
+// Second, a degree-based denominator matrix R is formed where
+// R(i,j) = deg(i) + deg(j) − B(i,j).
+// Lastly, the output matrix JC is computed as the element-wise ratio
+// JC = B ./ R.
+//
+
+//------------------------------------------------------------------------------
+// References:
+// (1) "Parallel Algorithms for Computing Jaccard Weights on Graphs using Linear Algebra," 
+//     in Proc. IEEE High Performance Extreme Computing Conference (HPEC), 2023.
+//     https://doi.org/10.1109/HPEC58863.2023.10363558
+//
+// (2) https://en.wikipedia.org/wiki/Jaccard_index
+
+
+#define LG_FREE_WORK                           \
+{                                              \
+    GrB_free(&R);                              \
+    GrB_free(&B);                              \
+    GrB_free(&D);                              \
+    GrB_free (&M) ;                            \
+}
+
+#define LG_FREE_ALL                            \
+{                                              \
+    LG_FREE_WORK ;                             \
+}
+
+#include "LG_internal.h"
+
+int LAGr_Jaccard
+(
+ // output
+ GrB_Matrix *JC,
+ // input:
+ LAGraph_Graph G,
+ bool all_pairs, 
+ char *msg
+ )
+{
+    GrB_Matrix B = NULL, R = NULL, D = NULL, M = NULL ;    
+    GrB_Index n;
+
+    //--------------------------------------------------------------------------
+    // check inputs
+    //--------------------------------------------------------------------------
+    LG_CLEAR_MSG ;
+    // error if G is directed (OK if G is directed but G->A is symmetric in structure
+    // error if G has self edges, or unknown
+    LG_ASSERT (JC != NULL, GrB_NULL_POINTER) ;
+    (*JC) = NULL ;
+    LG_TRY (LAGraph_CheckGraph (G, msg)) ;
+    LG_ASSERT_MSG (G->nself_edges == 0, LAGRAPH_NO_SELF_EDGES_ALLOWED, "G->nself_edges must be zero") ;
+    LG_ASSERT_MSG ((G->kind == LAGraph_ADJACENCY_UNDIRECTED || (G->kind == LAGraph_ADJACENCY_DIRECTED &&
+                    G->is_symmetric_structure == LAGraph_TRUE)),
+            LAGRAPH_SYMMETRIC_STRUCTURE_REQUIRED,
+            "G->A must be known to be symmetric") ;
+    //--------------------------------------------------------------------------
+    // degree vector deg
+    //--------------------------------------------------------------------------
+    GrB_Vector deg = G->out_degree ;
+    LG_ASSERT_MSG (deg != NULL, LAGRAPH_NOT_CACHED, "G->out_degree is required") ;
+    GrB_Matrix A = G->A ;
+    GRB_TRY( GrB_Matrix_nrows (&n, A) );
+    GrB_Type int_type  = (n > INT32_MAX) ? GrB_INT64 : GrB_INT32 ;
+
+    //--------------------------------------------------------------------------
+    // B is intersection matrix 
+    //--------------------------------------------------------------------------
+
+    // B(i,j) is the size of the intersection of the pattern of A(i,:) and
+    // A(:,j).  If all_pairs is true, B is computed for all entries in A^2.
+    // Otherwise, it is computed just for entries in triu(A).
+
+    // The final output matrix (B) is always upper triangular, in both cases.
+
+    GRB_TRY(GrB_Matrix_new(&B, GrB_FP64, n, n));
+
+    if (all_pairs)
+    {
+        // B = triu (A*A)
+        GRB_TRY(GrB_mxm(B, NULL, NULL, LAGraph_plus_one_fp64, A, A, NULL));
+        GRB_TRY( GrB_select(B, NULL, NULL, GrB_TRIU, B, (int64_t)0, NULL));
+    }
+    else
+    {
+        // B<triu(A)> = A*A'
+        GRB_TRY(GrB_Matrix_new(&M, GrB_BOOL, n, n));
+        GRB_TRY( GrB_select(M, NULL, NULL, GrB_TRIU, A, (int64_t)0, NULL));
+        GRB_TRY(GrB_mxm(B, M, NULL, LAGraph_plus_one_fp64, A, A, GrB_DESC_ST1));
+        GrB_free (&M) ;
+    }
+
+    //--------------------------------------------------------------------------
+    // R has summation of degree of corresponding nodes
+    // B is jaccard index B <- B / (R-B)
+    //--------------------------------------------------------------------------
+
+    // If deg vectors is sparse, make it dense
+    GrB_Index deg_nnz = 0;
+    GrB_Vector d = NULL;
+    GRB_TRY(GrB_Vector_nvals(&deg_nnz, deg));
+    if (deg_nnz < n ) {
+        GrB_Vector t = NULL;
+        GRB_TRY (GrB_Vector_dup(&t, deg));
+        GRB_TRY (GrB_assign (t, t, NULL, 0, GrB_ALL, n, GrB_DESC_SC)) ;
+        GRB_TRY (GrB_Vector_dup(&d, t));
+    }
+    else {
+        GRB_TRY (GrB_Vector_dup(&d, deg));
+    }
+
+    // D is degree matrix
+    GRB_TRY(GrB_Matrix_diag(&D, d, 0));   // use d here
+    GrB_Matrix_new(&R, int_type, n, n);
+
+    // R = B*D  -> R_ij = deg_j - b_ij
+    GRB_TRY(GrB_mxm(R, NULL, NULL, (int_type == GrB_INT64) ? GxB_PLUS_RMINUS_INT64 : GxB_PLUS_RMINUS_INT32, B, D, GrB_DESC_S));
+    // R = D*R  -> R_ij = deg_i + r_ij    
+    GRB_TRY(GrB_mxm(R, NULL, NULL, (int_type == GrB_INT64) ? GxB_PLUS_PLUS_INT64 : GxB_PLUS_PLUS_INT32, D, R, GrB_DESC_S));    
+
+    GRB_TRY(  GrB_eWiseMult(B, NULL, NULL, GrB_DIV_FP64, B, R, NULL) );
+    (*JC) = B;
+    B = NULL;           
+    LG_FREE_WORK;   
+
+    return (GrB_SUCCESS) ;
+}

experimental/algorithm/LAGr_MarkovClustering.c

@@ -40,6 +40,9 @@
 #include "LG_internal.h"
 #include <LAGraphX.h>
 
+// FIXME: make inputs e, i, pruning_threshold, convergence_threshold GrB_Scalar
+// and rename these parameters
+
 int LAGr_MarkovClustering(
     // output:
     GrB_Vector *c_f, // output cluster vector

experimental/algorithm/LAGr_MaxFlow.c

@@ -2,7 +2,7 @@
 // LAGr_MaxFlow: max flow
 //------------------------------------------------------------------------------
 
-// LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved.
+// LAGraph, (c) 2019-2026 by The LAGraph Contributors, All Rights Reserved.
 // SPDX-License-Identifier: BSD-2-Clause
 //
 // For additional details (including references to third party source code and
@@ -26,6 +26,27 @@
 // [2] D. Peries and T. Davis, "A parallel push-relabel maximum flow algorithm
 // in LAGraph and GraphBLAS", IEEE HPEC'25, Sept 2025.
 
+// FIXME: add something like the following:
+
+// LAGr_MaxFlow computes the maximum flow from the src to sink nodes, returning
+// the total maximum flow result as a single double scalar (f).  It can
+// optionally return the flows on each each, as the flow_mtx parameter.  This
+// matrix has no pre-flows left in it (which appear during the computation of a
+// push/relabel maxflow method).  However, computing the flow_mtx is very
+// expensive.  If not needed, pass in NULL instead of &flow_mtx (say), so it is
+// not computed.
+
+// FIXME ... ditto describe: res_mtx
+
+// G can be directed or undirected.  It must have edge weights of type
+// GrB_FP64, or types that can be typecasting to FP64.  This includes all
+// built-in types of GraphBLAS.  All computations are done in FP64.
+
+// FIXME: say something about convergence failure
+
+// See also LAGraph_MinCut to compute the minimum cut, using the results
+// from this method.
+
 // TODO: return the (optional) flow matrix can be costly in terms of run time.
 // The HPEC'25 results only benchmark the computation of the max flow, f.
 // Future work:  we plan on revising how the flow matrix is constructed.
@@ -607,11 +628,15 @@ void print_compareVec(const GrB_Vector vec) {
 // LAGraph_MaxFlow
 //------------------------------------------------------------------------------
 
+// FIXME: make output f a GrB_Scalar
+// FIXME: make input src,sink GrB_Scalar
+
 int LAGr_MaxFlow
 (
     // output:
     double *f,              // max flow from src node to sink node
     GrB_Matrix *flow_mtx,   // optional output flow matrix
+    GrB_Matrix *res_mtx,    // optional output for the residual matrix. Used in the min-cut
     // input:
     LAGraph_Graph G,        // graph to compute maxflow on
     GrB_Index src,          // source node
@@ -711,6 +736,10 @@ int LAGr_MaxFlow
   {
     (*flow_mtx) = NULL ;
   }
+  if (res_mtx != NULL)
+  {
+    (*res_mtx) = NULL ;
+  }
   LG_TRY(LAGraph_CheckGraph(G, msg));
   LG_ASSERT (f != NULL, GrB_NULL_POINTER) ;
   (*f) = 0;
@@ -720,6 +749,7 @@ int LAGr_MaxFlow
   LG_ASSERT_MSG(nrows == n, GrB_INVALID_VALUE, "Matrix must be square");
   LG_ASSERT_MSG(src < n && src >= 0 && sink < n && sink >= 0,
         GrB_INVALID_VALUE, "src and sink must be a value between [0, n)");
+  LG_ASSERT(G->emin_state != LAGraph_BOOLEAN_UNKNOWN, GrB_UNINITIALIZED_OBJECT) ;
   LG_ASSERT_MSG(G->emin > 0, GrB_INVALID_VALUE,
         "the edge weights (capacities) must be greater than 0");
 
@@ -1016,9 +1046,23 @@ int LAGr_MaxFlow
 
   for (int64_t iter = 0 ; n_active > 0 ; iter++)
   {
+
+    //--------------------------------------------------------------------------
+    // check for convergence failure
+    //--------------------------------------------------------------------------
+
+    // FIXME: test this case, and check if 2n is correct
+
+    // If the edge weights are vastly different in magnitude, it might be
+    // possible for this method to get stuck in an infinite loop.  The maximum
+    // number of iterations in theory is 2*n, so check if this is exceeded with
+    // 3*n as the upper bound just for good measure.
+
     #ifdef DBG
       printf ("iter: %ld, n_active %ld\n", iter, n_active) ;
     #endif
+    LG_ASSERT (iter < 3*n, LAGRAPH_CONVERGENCE_FAILURE) ;
+
     //--------------------------------------------------------------------------
     // Part 1: global relabeling
     //--------------------------------------------------------------------------
@@ -1160,6 +1204,13 @@ int LAGr_MaxFlow
     GRB_TRY(GrB_select(*flow_mtx, NULL, NULL, GrB_VALUEGT_FP64, *flow_mtx, 0, NULL));
   }
 
+  if(res_mtx != NULL){
+   GRB_TRY(GrB_Matrix_new(res_mtx, GrB_FP64, n, n)); 
+   GRB_TRY(GrB_apply(*res_mtx, NULL, NULL, GetResidual, R, NULL)) ;
+    // prune zeros and negative entries from R_hat
+   GRB_TRY(GrB_select(*res_mtx, NULL, NULL, GrB_VALUEGT_FP64, *res_mtx, 0, NULL)) ; 
+  }
+
   //----------------------------------------------------------------------------
   // for test coverage only
   //----------------------------------------------------------------------------

experimental/algorithm/LAGr_MaximumMatching.c

@@ -476,6 +476,9 @@ invert_2(GrB_Vector out,  // input/output
         GrB_free(&mateC);                                                      \
     }
 
+// FIXME: make col_init a GrB_Scalar?
+// FIXME need a bipartite LAGraph_Graph
+
 int LAGr_MaximumMatching(
     // outputs:
     GrB_Vector

experimental/algorithm/LAGr_Modularity.c

@@ -60,6 +60,9 @@
 #include "LG_internal.h"
 #include <LAGraphX.h>
 
+// FIXME: output modularity a GrB_Scalar
+// FIXME: input resolution a GrB_Scalar
+
 int LAGr_Modularity(
     // Outputs
     double *mod_handle, // output modularity