Parallelize tuple extraction in LAGraph_Matrix_Sum

The per-matrix extraction loop precomputes an offset (prefix-sum) array so
each matrix's tuples occupy a disjoint region of the shared (I, J, X) buffer.
That removes the loop-carried offset dependency and lets the extraction run
across LG_nthreads_outer threads with OpenMP; each GrB_Matrix_extractTuples
is still parallelized internally by GraphBLAS with LG_nthreads_inner threads
(the two-level nested model).  The public signature is unchanged: the thread
count follows the usual LAGraph convention via LAGraph_SetNumThreads.

Because GRB_TRY cannot return out of an OpenMP region, the first extraction
error is captured under a critical section and checked after the loop.

Adds test_Matrix_Sum_parallel, which sums many overlapping matrices with
multiple outer threads and compares against an independently accumulated
result.

Author: michelp

src/test/test_Matrix_Sum.c

@@ -193,6 +193,64 @@ void test_Matrix_Sum_brutal (void)
 }
 #endif
 
+//------------------------------------------------------------------------------
+// test_Matrix_Sum_parallel: exercise the parallel extraction path
+//------------------------------------------------------------------------------
+
+// Sum many matrices with multiple outer threads and confirm the result matches
+// an independently accumulated expected.  This stresses the disjoint-offset
+// extraction under real outer parallelism.
+
+void test_Matrix_Sum_parallel (void)
+{
+    setup ( ) ;
+
+    // request 4 outer threads (saving and restoring the prior settings)
+    int save_outer, save_inner ;
+    OK (LAGraph_GetNumThreads (&save_outer, &save_inner, msg)) ;
+    OK (LAGraph_SetNumThreads (4, save_inner, msg)) ;
+
+    #define NMAT 16
+    GrB_Matrix Mats [NMAT] ;
+    OK (GrB_Matrix_new (&Expected, GrB_FP64, 10, 10)) ;
+
+    for (int k = 0 ; k < NMAT ; k++)
+    {
+        // each matrix has 3 distinct (i,j) entries; the (7,7) entry is shared
+        // by every matrix and others overlap across matrices, so duplicates
+        // must be summed when the matrices are combined
+        GrB_Index Mi [ ] = { (GrB_Index) (k % 10), 2, 7 } ;
+        GrB_Index Mj [ ] = { 3, (GrB_Index) (k % 10), 7 } ;
+        double    Mx [ ] = { (double) (k + 1), 1, 2 } ;
+        Mats [k] = NULL ;
+        OK (GrB_Matrix_new (&Mats [k], GrB_FP64, 10, 10)) ;
+        OK (GrB_Matrix_build_FP64 (Mats [k], Mi, Mj, Mx, 3, NULL)) ;
+        // accumulate into Expected independently
+        OK (GrB_eWiseAdd (Expected, NULL, NULL, GrB_PLUS_FP64, Expected,
+            Mats [k], NULL)) ;
+    }
+
+    OK (LAGraph_Matrix_Sum (&C, Mats, NMAT, GrB_PLUS_FP64, msg)) ;
+
+    bool ok ;
+    OK (LAGraph_Matrix_IsEqual (&ok, C, Expected, msg)) ;
+    TEST_CHECK (ok) ;
+    TEST_MSG ("parallel sum of %d matrices did not match expected", NMAT) ;
+
+    for (int k = 0 ; k < NMAT ; k++)
+    {
+        OK (GrB_free (&Mats [k])) ;
+    }
+    OK (GrB_free (&C)) ;
+    OK (GrB_free (&Expected)) ;
+    #undef NMAT
+
+    // restore the original thread settings
+    OK (LAGraph_SetNumThreads (save_outer, save_inner, msg)) ;
+
+    teardown ( ) ;
+}
+
 //------------------------------------------------------------------------------
 // test_Matrix_Sum_failures: test error handling
 //------------------------------------------------------------------------------
@@ -266,6 +324,7 @@ TEST_LIST =
 {
     { "Matrix_Sum", test_Matrix_Sum },
     { "Matrix_Sum_types", test_Matrix_Sum_types },
+    { "Matrix_Sum_parallel", test_Matrix_Sum_parallel },
     { "Matrix_Sum_failures", test_Matrix_Sum_failures },
     #if LG_BRUTAL_TESTS
     { "Matrix_Sum_brutal", test_Matrix_Sum_brutal },

src/utility/LAGraph_Matrix_Sum.c

@@ -16,10 +16,14 @@
 //------------------------------------------------------------------------------
 
 // LAGraph_Matrix_Sum combines an array of matrices into a single matrix C.  It
-// computes the total number of entries across all inputs, allocates a single
-// tuple buffer (I, J, X) large enough to hold every entry, extracts the tuples
-// of each input matrix into that buffer, and then calls GrB_Matrix_build with
-// the binary operator dup to combine any duplicate (i,j) entries.  With dup =
+// computes the total number of entries across all inputs and the offset at
+// which each matrix's tuples begin in a single shared tuple buffer (I, J, X)
+// large enough to hold every entry.  Because each matrix writes to a disjoint
+// region of that buffer, the per-matrix extraction is parallelized across
+// LG_nthreads_outer threads with OpenMP; SuiteSparse:GraphBLAS parallelizes
+// each GrB_Matrix_extractTuples internally with LG_nthreads_inner threads.  The
+// concatenated tuples are then passed to GrB_Matrix_build, using the binary
+// operator dup to combine any duplicate (i,j) entries.  With dup =
 // GrB_PLUS_FP64 (for example) this computes the element-wise sum of all input
 // matrices.
 
@@ -31,6 +35,7 @@
     LAGraph_Free ((void **) &I, NULL) ;         \
     LAGraph_Free ((void **) &J, NULL) ;         \
     LAGraph_Free ((void **) &X, NULL) ;         \
+    LAGraph_Free ((void **) &Offsets, NULL) ;   \
 }
 
 #define LG_FREE_ALL                             \
@@ -58,7 +63,7 @@ int LAGraph_Matrix_Sum
     //--------------------------------------------------------------------------
 
     LG_CLEAR_MSG ;
-    GrB_Index *I = NULL, *J = NULL ;
+    GrB_Index *I = NULL, *J = NULL, *Offsets = NULL ;
     void *X = NULL ;
     LG_ASSERT_MSG (C != NULL, GrB_NULL_POINTER, "&C != NULL") ;
     LG_ASSERT (Matrices != NULL, GrB_NULL_POINTER) ;
@@ -78,10 +83,17 @@ int LAGraph_Matrix_Sum
     GRB_TRY (GrB_get (Matrices [0], &typecode, GrB_EL_TYPE_CODE)) ;
 
     //--------------------------------------------------------------------------
-    // validate every matrix and accumulate the total number of entries
+    // validate every matrix and compute where its tuples begin in the buffer
     //--------------------------------------------------------------------------
 
-    GrB_Index total = 0 ;
+    // Offsets [k] is the position in (I, J, X) at which the tuples of matrix k
+    // begin; Offsets [k+1] - Offsets [k] is its number of entries.  This prefix
+    // sum gives each matrix a disjoint buffer region so the extraction below
+    // can run in parallel without any data races.
+
+    LG_TRY (LAGraph_Malloc ((void **) &Offsets, nmatrices + 1,
+        sizeof (GrB_Index), msg)) ;
+    Offsets [0] = 0 ;
     for (GrB_Index k = 0 ; k < nmatrices ; k++)
     {
         GrB_Matrix Ak = Matrices [k] ;
@@ -96,8 +108,9 @@ int LAGraph_Matrix_Sum
         LG_ASSERT_MSG (code == typecode, GrB_DOMAIN_MISMATCH,
             "all input matrices must have the same type") ;
         GRB_TRY (GrB_Matrix_nvals (&n, Ak)) ;
-        total += n ;
+        Offsets [k+1] = Offsets [k] + n ;
     }
+    GrB_Index total = Offsets [nmatrices] ;
 
     //--------------------------------------------------------------------------
     // allocate the shared row/column index buffers (guard against size 0)
@@ -107,14 +120,25 @@ int LAGraph_Matrix_Sum
     LG_TRY (LAGraph_Malloc ((void **) &I, alloc, sizeof (GrB_Index), msg)) ;
     LG_TRY (LAGraph_Malloc ((void **) &J, alloc, sizeof (GrB_Index), msg)) ;
 
+    //--------------------------------------------------------------------------
+    // determine the number of threads for the outer extraction loop
+    //--------------------------------------------------------------------------
+
+    int nthreads = LG_nthreads_outer ;
+    nthreads = LAGRAPH_MIN (nthreads, (int) nmatrices) ;
+    nthreads = LAGRAPH_MAX (nthreads, 1) ;
+
     //--------------------------------------------------------------------------
     // extract tuples from every matrix, then build the result
     //--------------------------------------------------------------------------
 
     // For each built-in type: allocate the value buffer X with the correct
-    // element size, extract the tuples of every input matrix into the shared
-    // buffer at the running offset, create C, and build it with the dup
-    // operator to combine duplicate (i,j) entries.
+    // element size, extract the tuples of every input matrix into its disjoint
+    // region of the shared buffer (in parallel, since the regions never
+    // overlap), create C, and build it with the dup operator to combine
+    // duplicate (i,j) entries.  GRB_TRY cannot be used inside an OpenMP region
+    // (it returns from the function), so the first error is captured into
+    // sum_status under a critical section and checked after the region.
 
     #define LG_SUM_CASE(code, ctype, gtype, suffix)                          \
         case code :                                                          \
@@ -123,18 +147,23 @@ int LAGraph_Matrix_Sum
             LG_TRY (LAGraph_Malloc ((void **) &Xt, alloc, sizeof (ctype),    \
                 msg)) ;                                                      \
             X = (void *) Xt ;                                                \
-            GrB_Index offset = 0 ;                                           \
-            for (GrB_Index k = 0 ; k < nmatrices ; k++)                      \
+            int sum_status = GrB_SUCCESS ;                                   \
+            int64_t k ;                                                      \
+            _Pragma ("omp parallel for num_threads(nthreads) schedule(dynamic,1)") \
+            for (k = 0 ; k < (int64_t) nmatrices ; k++)                      \
             {                                                                \
-                GrB_Index n, got ;                                           \
-                GRB_TRY (GrB_Matrix_nvals (&n, Matrices [k])) ;              \
-                if (n == 0) continue ;                                       \
-                got = n ;                                                    \
-                GRB_TRY (GrB_Matrix_extractTuples_ ## suffix (               \
-                    I + offset, J + offset, Xt + offset, &got,               \
-                    Matrices [k])) ;                                         \
-                offset += n ;                                                \
+                GrB_Index off = Offsets [k] ;                                \
+                GrB_Index got = Offsets [k+1] - off ;                        \
+                if (got == 0) continue ;                                     \
+                GrB_Info info = GrB_Matrix_extractTuples_ ## suffix (        \
+                    I + off, J + off, Xt + off, &got, Matrices [k]) ;        \
+                if (info < GrB_SUCCESS)                                      \
+                {                                                            \
+                    _Pragma ("omp critical")                                 \
+                    { if (sum_status >= GrB_SUCCESS) sum_status = info ; }   \
+                }                                                            \
             }                                                                \
+            GRB_TRY (sum_status) ;                                           \
             GRB_TRY (GrB_Matrix_new (C, gtype, nrows, ncols)) ;              \
             GRB_TRY (GrB_Matrix_build_ ## suffix (*C, I, J, Xt, total,       \
                 dup)) ;                                                      \