From 43ee0acef248d8db99dc49edebf0c919c3d95ccf Mon Sep 17 00:00:00 2001 From: Michel Pelletier Date: Thu, 11 Jun 2026 18:24:27 -0700 Subject: [PATCH 1/3] Add LAGraph_Matrix_Sum utility LAGraph_Matrix_Sum combines an array of GrB_Matrix objects into a single matrix using a binary operator to resolve duplicate entries. It computes the total number of entries across all inputs, allocates one shared tuple buffer (I, J, X), extracts the tuples of each input matrix into the buffer, and calls GrB_Matrix_build with the dup operator to combine duplicate (i,j) coordinates. With dup = GrB_PLUS_FP64 (for example) this computes the element-wise sum of all the matrices. All input matrices must have identical dimensions and the same built-in type; user-defined types return GrB_NOT_IMPLEMENTED. Includes a test (src/test/test_Matrix_Sum.c) covering correctness, all built-in type branches, error handling, and a brutal malloc-failure variant. --- Config/LAGraph.h.in | 41 +++++ include/LAGraph.h | 41 +++++ src/test/test_Matrix_Sum.c | 274 +++++++++++++++++++++++++++++++ src/utility/LAGraph_Matrix_Sum.c | 170 +++++++++++++++++++ 4 files changed, 526 insertions(+) create mode 100644 src/test/test_Matrix_Sum.c create mode 100644 src/utility/LAGraph_Matrix_Sum.c diff --git a/Config/LAGraph.h.in b/Config/LAGraph.h.in index a46c06c8eb..77aed33386 100644 --- a/Config/LAGraph.h.in +++ b/Config/LAGraph.h.in @@ -1572,6 +1572,47 @@ int LAGraph_Matrix_Structure char *msg ) ; +//------------------------------------------------------------------------------ +// LAGraph_Matrix_Sum: sum an array of matrices with a binary operator +//------------------------------------------------------------------------------ + +/** LAGraph_Matrix_Sum: combines an array of matrices into a single matrix C. + * The tuples of all input matrices are concatenated into a single buffer and + * 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 the matrices. + * + * All input matrices must have identical dimensions and identical built-in + * type; C is created with that same type and dimensions. + * + * @param[out] C the resulting matrix. + * @param[in] Matrices array of nmatrices input matrices. + * @param[in] nmatrices number of matrices in Matrices (must be >= 1). + * @param[in] dup binary operator to combine duplicate (i,j) entries; + * if NULL, duplicates are handled per GrB_Matrix_build. + * @param[in,out] msg any error messages. + * + * @retval GrB_SUCCESS if successful. + * @retval GrB_NULL_POINTER if C, Matrices, or any Matrices[k] is NULL. + * @retval GrB_INVALID_VALUE if nmatrices is zero. + * @retval GrB_DIMENSION_MISMATCH if the inputs do not all share dimensions. + * @retval GrB_DOMAIN_MISMATCH if the inputs do not all share a type. + * @retval GrB_NOT_IMPLEMENTED if the matrix type is user-defined. + * @returns any GraphBLAS errors that may have been encountered. + */ + +LAGRAPH_PUBLIC +int LAGraph_Matrix_Sum +( + // output: + GrB_Matrix *C, // result = combination of all input matrices + // input: + GrB_Matrix *Matrices, // array of nmatrices input matrices + GrB_Index nmatrices, // number of matrices in the array (must be >= 1) + GrB_BinaryOp dup, // operator to combine duplicate (i,j) entries + char *msg +) ; + //------------------------------------------------------------------------------ // LAGraph_Vector_Structure: return the structure of a vector //------------------------------------------------------------------------------ diff --git a/include/LAGraph.h b/include/LAGraph.h index ea88c7c7dc..073f2bdf6c 100644 --- a/include/LAGraph.h +++ b/include/LAGraph.h @@ -1572,6 +1572,47 @@ int LAGraph_Matrix_Structure char *msg ) ; +//------------------------------------------------------------------------------ +// LAGraph_Matrix_Sum: sum an array of matrices with a binary operator +//------------------------------------------------------------------------------ + +/** LAGraph_Matrix_Sum: combines an array of matrices into a single matrix C. + * The tuples of all input matrices are concatenated into a single buffer and + * 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 the matrices. + * + * All input matrices must have identical dimensions and identical built-in + * type; C is created with that same type and dimensions. + * + * @param[out] C the resulting matrix. + * @param[in] Matrices array of nmatrices input matrices. + * @param[in] nmatrices number of matrices in Matrices (must be >= 1). + * @param[in] dup binary operator to combine duplicate (i,j) entries; + * if NULL, duplicates are handled per GrB_Matrix_build. + * @param[in,out] msg any error messages. + * + * @retval GrB_SUCCESS if successful. + * @retval GrB_NULL_POINTER if C, Matrices, or any Matrices[k] is NULL. + * @retval GrB_INVALID_VALUE if nmatrices is zero. + * @retval GrB_DIMENSION_MISMATCH if the inputs do not all share dimensions. + * @retval GrB_DOMAIN_MISMATCH if the inputs do not all share a type. + * @retval GrB_NOT_IMPLEMENTED if the matrix type is user-defined. + * @returns any GraphBLAS errors that may have been encountered. + */ + +LAGRAPH_PUBLIC +int LAGraph_Matrix_Sum +( + // output: + GrB_Matrix *C, // result = combination of all input matrices + // input: + GrB_Matrix *Matrices, // array of nmatrices input matrices + GrB_Index nmatrices, // number of matrices in the array (must be >= 1) + GrB_BinaryOp dup, // operator to combine duplicate (i,j) entries + char *msg +) ; + //------------------------------------------------------------------------------ // LAGraph_Vector_Structure: return the structure of a vector //------------------------------------------------------------------------------ diff --git a/src/test/test_Matrix_Sum.c b/src/test/test_Matrix_Sum.c new file mode 100644 index 0000000000..a46e78eb98 --- /dev/null +++ b/src/test/test_Matrix_Sum.c @@ -0,0 +1,274 @@ +//------------------------------------------------------------------------------ +// LAGraph/src/test/test_Matrix_Sum.c: test LAGraph_Matrix_Sum +//------------------------------------------------------------------------------ + +// 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 Michel Pelletier. + +//------------------------------------------------------------------------------ + +#include "LAGraph_test.h" +#include "LG_internal.h" + +//------------------------------------------------------------------------------ +// global variables +//------------------------------------------------------------------------------ + +char msg [LAGRAPH_MSG_LEN] ; +GrB_Matrix A = NULL, B = NULL, C = NULL, Expected = NULL ; + +//------------------------------------------------------------------------------ +// setup and teardown +//------------------------------------------------------------------------------ + +void setup (void) +{ + OK (LAGraph_Init (msg)) ; +} + +void teardown (void) +{ + OK (LAGraph_Finalize (msg)) ; +} + +//------------------------------------------------------------------------------ +// make_A, make_B: construct two 4x4 FP64 matrices with overlapping and +// distinct (i,j) coordinates +//------------------------------------------------------------------------------ + +static void make_AB (void) +{ + // A has entries at (0,0), (1,1), (2,2), (0,3) + GrB_Index Ai [ ] = { 0, 1, 2, 0 } ; + GrB_Index Aj [ ] = { 0, 1, 2, 3 } ; + double Ax [ ] = { 1, 2, 3, 4 } ; + OK (GrB_Matrix_new (&A, GrB_FP64, 4, 4)) ; + OK (GrB_Matrix_build_FP64 (A, Ai, Aj, Ax, 4, NULL)) ; + + // B has entries at (0,0), (1,1), (3,3), (2,0) + // shares (0,0) and (1,1) with A, distinct elsewhere + GrB_Index Bi [ ] = { 0, 1, 3, 2 } ; + GrB_Index Bj [ ] = { 0, 1, 3, 0 } ; + double Bx [ ] = { 10, 20, 30, 40 } ; + OK (GrB_Matrix_new (&B, GrB_FP64, 4, 4)) ; + OK (GrB_Matrix_build_FP64 (B, Bi, Bj, Bx, 4, NULL)) ; +} + +//------------------------------------------------------------------------------ +// test_Matrix_Sum: basic correctness +//------------------------------------------------------------------------------ + +void test_Matrix_Sum (void) +{ + setup ( ) ; + + make_AB ( ) ; + + // Expected = A + B (element-wise add, set union) + OK (GrB_Matrix_new (&Expected, GrB_FP64, 4, 4)) ; + OK (GrB_eWiseAdd (Expected, NULL, NULL, GrB_PLUS_FP64, A, B, NULL)) ; + + // C = sum of {A, B} using GrB_PLUS_FP64 for duplicates + GrB_Matrix Mats [2] = { A, B } ; + OK (LAGraph_Matrix_Sum (&C, Mats, 2, GrB_PLUS_FP64, msg)) ; + + bool ok ; + OK (LAGraph_Matrix_IsEqual (&ok, C, Expected, msg)) ; + TEST_CHECK (ok) ; + TEST_MSG ("sum of {A,B} did not equal A+B") ; + OK (GrB_free (&C)) ; + + // single-matrix array: result must be a copy of A + GrB_Matrix One [1] = { A } ; + OK (LAGraph_Matrix_Sum (&C, One, 1, GrB_PLUS_FP64, msg)) ; + OK (LAGraph_Matrix_IsEqual (&ok, C, A, msg)) ; + TEST_CHECK (ok) ; + TEST_MSG ("sum of {A} did not equal A") ; + OK (GrB_free (&C)) ; + + // an empty matrix in the array contributes nothing + GrB_Matrix Empty = NULL ; + OK (GrB_Matrix_new (&Empty, GrB_FP64, 4, 4)) ; + GrB_Matrix WithEmpty [3] = { A, Empty, B } ; + OK (LAGraph_Matrix_Sum (&C, WithEmpty, 3, GrB_PLUS_FP64, msg)) ; + OK (LAGraph_Matrix_IsEqual (&ok, C, Expected, msg)) ; + TEST_CHECK (ok) ; + TEST_MSG ("sum of {A,Empty,B} did not equal A+B") ; + OK (GrB_free (&Empty)) ; + OK (GrB_free (&C)) ; + + OK (GrB_free (&A)) ; + OK (GrB_free (&B)) ; + OK (GrB_free (&Expected)) ; + + teardown ( ) ; +} + +//------------------------------------------------------------------------------ +// test_Matrix_Sum_types: exercise every built-in type branch +//------------------------------------------------------------------------------ + +void test_Matrix_Sum_types (void) +{ + setup ( ) ; + + // (type, GrB_PLUS operator) pairs, one per built-in type branch + GrB_Type types [ ] = { GrB_BOOL, GrB_INT8, GrB_INT16, GrB_INT32, + GrB_INT64, GrB_UINT8, GrB_UINT16, GrB_UINT32, GrB_UINT64, + GrB_FP32, GrB_FP64 } ; + GrB_BinaryOp ops [ ] = { GrB_LOR, GrB_PLUS_INT8, GrB_PLUS_INT16, + GrB_PLUS_INT32, GrB_PLUS_INT64, GrB_PLUS_UINT8, GrB_PLUS_UINT16, + GrB_PLUS_UINT32, GrB_PLUS_UINT64, GrB_PLUS_FP32, GrB_PLUS_FP64 } ; + int ntypes = sizeof (types) / sizeof (types [0]) ; + + GrB_Index Ai [ ] = { 0, 1, 2 } ; + GrB_Index Aj [ ] = { 0, 1, 0 } ; + GrB_Index Bi [ ] = { 0, 1, 0 } ; + GrB_Index Bj [ ] = { 0, 1, 2 } ; + + for (int t = 0 ; t < ntypes ; t++) + { + OK (GrB_Matrix_new (&A, types [t], 3, 3)) ; + OK (GrB_Matrix_new (&B, types [t], 3, 3)) ; + // build with INT32 values; GraphBLAS typecasts to the matrix type + int32_t Ax [ ] = { 1, 1, 1 } ; + int32_t Bx [ ] = { 1, 1, 1 } ; + OK (GrB_Matrix_build_INT32 (A, Ai, Aj, Ax, 3, NULL)) ; + OK (GrB_Matrix_build_INT32 (B, Bi, Bj, Bx, 3, NULL)) ; + + OK (GrB_Matrix_new (&Expected, types [t], 3, 3)) ; + OK (GrB_eWiseAdd (Expected, NULL, NULL, ops [t], A, B, NULL)) ; + + GrB_Matrix Mats [2] = { A, B } ; + OK (LAGraph_Matrix_Sum (&C, Mats, 2, ops [t], msg)) ; + + bool ok ; + OK (LAGraph_Matrix_IsEqual (&ok, C, Expected, msg)) ; + TEST_CHECK (ok) ; + TEST_MSG ("type index %d failed", t) ; + + OK (GrB_free (&A)) ; + OK (GrB_free (&B)) ; + OK (GrB_free (&C)) ; + OK (GrB_free (&Expected)) ; + } + + teardown ( ) ; +} + +//------------------------------------------------------------------------------ +// test_Matrix_Sum_brutal +//------------------------------------------------------------------------------ + +#if LG_BRUTAL_TESTS +void test_Matrix_Sum_brutal (void) +{ + OK (LG_brutal_setup (msg)) ; + + make_AB ( ) ; + OK (GrB_Matrix_new (&Expected, GrB_FP64, 4, 4)) ; + OK (GrB_eWiseAdd (Expected, NULL, NULL, GrB_PLUS_FP64, A, B, NULL)) ; + + GrB_Matrix Mats [2] = { A, B } ; + LG_BRUTAL (LAGraph_Matrix_Sum (&C, Mats, 2, GrB_PLUS_FP64, msg)) ; + + bool ok ; + OK (LAGraph_Matrix_IsEqual (&ok, C, Expected, msg)) ; + TEST_CHECK (ok) ; + + OK (GrB_free (&A)) ; + OK (GrB_free (&B)) ; + OK (GrB_free (&C)) ; + OK (GrB_free (&Expected)) ; + + OK (LG_brutal_teardown (msg)) ; +} +#endif + +//------------------------------------------------------------------------------ +// test_Matrix_Sum_failures: test error handling +//------------------------------------------------------------------------------ + +void test_Matrix_Sum_failures (void) +{ + setup ( ) ; + + make_AB ( ) ; + GrB_Matrix Mats [2] = { A, B } ; + int result ; + + // NULL output + result = LAGraph_Matrix_Sum (NULL, Mats, 2, GrB_PLUS_FP64, msg) ; + TEST_CHECK (result == GrB_NULL_POINTER) ; + + // NULL Matrices array + C = NULL ; + result = LAGraph_Matrix_Sum (&C, NULL, 2, GrB_PLUS_FP64, msg) ; + TEST_CHECK (result == GrB_NULL_POINTER) ; + TEST_CHECK (C == NULL) ; + + // nmatrices == 0 + result = LAGraph_Matrix_Sum (&C, Mats, 0, GrB_PLUS_FP64, msg) ; + TEST_CHECK (result == GrB_INVALID_VALUE) ; + + // NULL entry in the array + GrB_Matrix WithNull [2] = { A, NULL } ; + result = LAGraph_Matrix_Sum (&C, WithNull, 2, GrB_PLUS_FP64, msg) ; + TEST_CHECK (result == GrB_NULL_POINTER) ; + + // dimension mismatch + GrB_Matrix D = NULL ; + OK (GrB_Matrix_new (&D, GrB_FP64, 5, 5)) ; + GrB_Matrix BadDim [2] = { A, D } ; + result = LAGraph_Matrix_Sum (&C, BadDim, 2, GrB_PLUS_FP64, msg) ; + TEST_CHECK (result == GrB_DIMENSION_MISMATCH) ; + OK (GrB_free (&D)) ; + + // type mismatch + GrB_Matrix E = NULL ; + OK (GrB_Matrix_new (&E, GrB_INT32, 4, 4)) ; + GrB_Matrix BadType [2] = { A, E } ; + result = LAGraph_Matrix_Sum (&C, BadType, 2, GrB_PLUS_FP64, msg) ; + TEST_CHECK (result == GrB_DOMAIN_MISMATCH) ; + OK (GrB_free (&E)) ; + + // user-defined type not supported + typedef struct { double x, y ; } udt_t ; + GrB_Type UDT = NULL ; + OK (GrB_Type_new (&UDT, sizeof (udt_t))) ; + GrB_Matrix U = NULL ; + OK (GrB_Matrix_new (&U, UDT, 4, 4)) ; + GrB_Matrix Udt [1] = { U } ; + result = LAGraph_Matrix_Sum (&C, Udt, 1, NULL, msg) ; + TEST_CHECK (result == GrB_NOT_IMPLEMENTED) ; + OK (GrB_free (&U)) ; + OK (GrB_free (&UDT)) ; + + OK (GrB_free (&A)) ; + OK (GrB_free (&B)) ; + + teardown ( ) ; +} + +//------------------------------------------------------------------------------ +// TEST_LIST: the list of tasks for this entire test +//------------------------------------------------------------------------------ + +TEST_LIST = +{ + { "Matrix_Sum", test_Matrix_Sum }, + { "Matrix_Sum_types", test_Matrix_Sum_types }, + { "Matrix_Sum_failures", test_Matrix_Sum_failures }, + #if LG_BRUTAL_TESTS + { "Matrix_Sum_brutal", test_Matrix_Sum_brutal }, + #endif + { NULL, NULL } +} ; diff --git a/src/utility/LAGraph_Matrix_Sum.c b/src/utility/LAGraph_Matrix_Sum.c new file mode 100644 index 0000000000..7532f8fb72 --- /dev/null +++ b/src/utility/LAGraph_Matrix_Sum.c @@ -0,0 +1,170 @@ +//------------------------------------------------------------------------------ +// LAGraph_Matrix_Sum: sum an array of matrices with a binary operator +//------------------------------------------------------------------------------ + +// 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 Michel Pelletier. + +//------------------------------------------------------------------------------ + +// 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 = +// GrB_PLUS_FP64 (for example) this computes the element-wise sum of all input +// matrices. + +// All input matrices must have identical dimensions and identical built-in +// type; C is created with that same type and dimensions. + +#define LG_FREE_WORK \ +{ \ + LAGraph_Free ((void **) &I, NULL) ; \ + LAGraph_Free ((void **) &J, NULL) ; \ + LAGraph_Free ((void **) &X, NULL) ; \ +} + +#define LG_FREE_ALL \ +{ \ + LG_FREE_WORK ; \ + GrB_free (C) ; \ +} + +#include "LG_internal.h" + +int LAGraph_Matrix_Sum +( + // output: + GrB_Matrix *C, // result = combination of all input matrices + // input: + GrB_Matrix *Matrices, // array of nmatrices input matrices + GrB_Index nmatrices, // number of matrices in the array (must be >= 1) + GrB_BinaryOp dup, // operator to combine duplicate (i,j) entries + char *msg +) +{ + + //-------------------------------------------------------------------------- + // check inputs + //-------------------------------------------------------------------------- + + LG_CLEAR_MSG ; + GrB_Index *I = NULL, *J = NULL ; + void *X = NULL ; + LG_ASSERT_MSG (C != NULL, GrB_NULL_POINTER, "&C != NULL") ; + LG_ASSERT (Matrices != NULL, GrB_NULL_POINTER) ; + (*C) = NULL ; + LG_ASSERT_MSG (nmatrices >= 1, GrB_INVALID_VALUE, + "nmatrices must be >= 1") ; + LG_ASSERT (Matrices [0] != NULL, GrB_NULL_POINTER) ; + + //-------------------------------------------------------------------------- + // determine the reference dimensions and type from the first matrix + //-------------------------------------------------------------------------- + + GrB_Index nrows, ncols ; + int32_t typecode ; + GRB_TRY (GrB_Matrix_nrows (&nrows, Matrices [0])) ; + GRB_TRY (GrB_Matrix_ncols (&ncols, Matrices [0])) ; + GRB_TRY (GrB_get (Matrices [0], &typecode, GrB_EL_TYPE_CODE)) ; + + //-------------------------------------------------------------------------- + // validate every matrix and accumulate the total number of entries + //-------------------------------------------------------------------------- + + GrB_Index total = 0 ; + for (GrB_Index k = 0 ; k < nmatrices ; k++) + { + GrB_Matrix Ak = Matrices [k] ; + LG_ASSERT (Ak != NULL, GrB_NULL_POINTER) ; + GrB_Index r, c, n ; + int32_t code ; + GRB_TRY (GrB_Matrix_nrows (&r, Ak)) ; + GRB_TRY (GrB_Matrix_ncols (&c, Ak)) ; + LG_ASSERT_MSG (r == nrows && c == ncols, GrB_DIMENSION_MISMATCH, + "all input matrices must have the same dimensions") ; + GRB_TRY (GrB_get (Ak, &code, GrB_EL_TYPE_CODE)) ; + 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 ; + } + + //-------------------------------------------------------------------------- + // allocate the shared row/column index buffers (guard against size 0) + //-------------------------------------------------------------------------- + + GrB_Index alloc = (total == 0) ? 1 : total ; + LG_TRY (LAGraph_Malloc ((void **) &I, alloc, sizeof (GrB_Index), msg)) ; + LG_TRY (LAGraph_Malloc ((void **) &J, alloc, sizeof (GrB_Index), msg)) ; + + //-------------------------------------------------------------------------- + // 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. + + #define LG_SUM_CASE(code, ctype, gtype, suffix) \ + case code : \ + { \ + ctype *Xt = NULL ; \ + 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++) \ + { \ + 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_TRY (GrB_Matrix_new (C, gtype, nrows, ncols)) ; \ + GRB_TRY (GrB_Matrix_build_ ## suffix (*C, I, J, Xt, total, \ + dup)) ; \ + } \ + break ; + + switch (typecode) + { + LG_SUM_CASE (GrB_BOOL_CODE, bool, GrB_BOOL, BOOL ) + LG_SUM_CASE (GrB_INT8_CODE, int8_t, GrB_INT8, INT8 ) + LG_SUM_CASE (GrB_INT16_CODE, int16_t, GrB_INT16, INT16 ) + LG_SUM_CASE (GrB_INT32_CODE, int32_t, GrB_INT32, INT32 ) + LG_SUM_CASE (GrB_INT64_CODE, int64_t, GrB_INT64, INT64 ) + LG_SUM_CASE (GrB_UINT8_CODE, uint8_t, GrB_UINT8, UINT8 ) + LG_SUM_CASE (GrB_UINT16_CODE, uint16_t, GrB_UINT16, UINT16) + LG_SUM_CASE (GrB_UINT32_CODE, uint32_t, GrB_UINT32, UINT32) + LG_SUM_CASE (GrB_UINT64_CODE, uint64_t, GrB_UINT64, UINT64) + LG_SUM_CASE (GrB_FP32_CODE, float, GrB_FP32, FP32 ) + LG_SUM_CASE (GrB_FP64_CODE, double, GrB_FP64, FP64 ) + default : + LG_ASSERT_MSG (false, GrB_NOT_IMPLEMENTED, + "user-defined types are not supported") ; + } + + #undef LG_SUM_CASE + + //-------------------------------------------------------------------------- + // free workspace and return result + //-------------------------------------------------------------------------- + + LG_FREE_WORK ; + return (GrB_SUCCESS) ; +} From 1895a885bea2cfc9cae03430bce42d565b7d8244 Mon Sep 17 00:00:00 2001 From: Michel Pelletier Date: Thu, 11 Jun 2026 21:35:34 -0700 Subject: [PATCH 2/3] 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. --- src/test/test_Matrix_Sum.c | 59 ++++++++++++++++++++++++++ src/utility/LAGraph_Matrix_Sum.c | 71 ++++++++++++++++++++++---------- 2 files changed, 109 insertions(+), 21 deletions(-) diff --git a/src/test/test_Matrix_Sum.c b/src/test/test_Matrix_Sum.c index a46e78eb98..3863ced258 100644 --- a/src/test/test_Matrix_Sum.c +++ b/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 }, diff --git a/src/utility/LAGraph_Matrix_Sum.c b/src/utility/LAGraph_Matrix_Sum.c index 7532f8fb72..26dbf5f773 100644 --- a/src/utility/LAGraph_Matrix_Sum.c +++ b/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)) ; \ From 6d3dc830f592bdae605a58cb59b5ffc57ee2d8d8 Mon Sep 17 00:00:00 2001 From: Michel Pelletier Date: Thu, 25 Jun 2026 11:09:56 -0700 Subject: [PATCH 3/3] Add LAGraph_Matrix_Binary_Sum utility Adds LAGraph_Matrix_Binary_Sum, which combines an array of matrices into a single matrix with the same interface and result as LAGraph_Matrix_Sum but a different technique: a pairwise binary reduction tree built from GrB_eWiseAdd (per the GraphChallenge "Anonymized Network Sensing" paper, "Binary Summation of Traffic Matrices"), rather than one large extract/build. At each level the matrices are summed in disjoint adjacent pairs using the dup operator; an unpaired (odd) trailing matrix is carried up to the next level unchanged, until a single matrix remains. The independent pair-sums within a level are parallelized across LG_nthreads_outer threads. Working on the smaller intermediate matrices of the tree keeps the active data in faster memory. All intermediate matrices are tracked in a single pre-NULLed Pool array with deterministic per-pair slots, so cleanup is a simple sweep that is correct on any failure path. Unlike LAGraph_Matrix_Sum, dup must be non-NULL since GrB_eWiseAdd requires a binary operator. Includes tests for correctness, odd counts (carry path), all built-in types, parallel reduction with multiple outer threads, error handling, and a brutal malloc-failure variant. --- Config/LAGraph.h.in | 51 +++ include/LAGraph.h | 51 +++ src/test/test_Matrix_Binary_Sum.c | 419 ++++++++++++++++++++++++ src/utility/LAGraph_Matrix_Binary_Sum.c | 272 +++++++++++++++ 4 files changed, 793 insertions(+) create mode 100644 src/test/test_Matrix_Binary_Sum.c create mode 100644 src/utility/LAGraph_Matrix_Binary_Sum.c diff --git a/Config/LAGraph.h.in b/Config/LAGraph.h.in index 77aed33386..b5e6320569 100644 --- a/Config/LAGraph.h.in +++ b/Config/LAGraph.h.in @@ -1613,6 +1613,57 @@ int LAGraph_Matrix_Sum char *msg ) ; +//------------------------------------------------------------------------------ +// LAGraph_Matrix_Binary_Sum: sum an array of matrices by binary reduction +//------------------------------------------------------------------------------ + +/** LAGraph_Matrix_Binary_Sum: combines an array of matrices into a single + * matrix C, computing the same result as LAGraph_Matrix_Sum but with a + * different technique: a pairwise binary reduction tree built from + * GrB_eWiseAdd (as in the GraphChallenge "Anonymized Network Sensing" paper, + * "Binary Summation of Traffic Matrices"). At each level the matrices are + * summed in disjoint adjacent pairs using the binary operator dup; an unpaired + * (odd) trailing matrix is carried up to the next level unchanged, until a + * single matrix remains. The dup operator has the same set-union semantics as + * the dup of GrB_Matrix_build (applied where both inputs have an entry, lone + * entries pass through), so the result equals that of LAGraph_Matrix_Sum. + * With dup = GrB_PLUS_FP64 (for example) this computes the element-wise sum of + * all the matrices. The independent pair-sums within a level are parallelized + * across LG_nthreads_outer threads. + * + * All input matrices must have identical dimensions and identical built-in + * type; C is created with that same type and dimensions. Unlike + * LAGraph_Matrix_Sum, dup must be non-NULL, since GrB_eWiseAdd requires a + * binary operator. + * + * @param[out] C the resulting matrix. + * @param[in] Matrices array of nmatrices input matrices. + * @param[in] nmatrices number of matrices in Matrices (must be >= 1). + * @param[in] dup binary operator to combine (i,j) entries; must be + * non-NULL. + * @param[in,out] msg any error messages. + * + * @retval GrB_SUCCESS if successful. + * @retval GrB_NULL_POINTER if C, Matrices, any Matrices[k], or dup is NULL. + * @retval GrB_INVALID_VALUE if nmatrices is zero. + * @retval GrB_DIMENSION_MISMATCH if the inputs do not all share dimensions. + * @retval GrB_DOMAIN_MISMATCH if the inputs do not all share a type. + * @retval GrB_NOT_IMPLEMENTED if the matrix type is user-defined. + * @returns any GraphBLAS errors that may have been encountered. + */ + +LAGRAPH_PUBLIC +int LAGraph_Matrix_Binary_Sum +( + // output: + GrB_Matrix *C, // result = combination of all input matrices + // input: + GrB_Matrix *Matrices, // array of nmatrices input matrices + GrB_Index nmatrices, // number of matrices in the array (must be >= 1) + GrB_BinaryOp dup, // operator to combine (i,j) entries (must be != NULL) + char *msg +) ; + //------------------------------------------------------------------------------ // LAGraph_Vector_Structure: return the structure of a vector //------------------------------------------------------------------------------ diff --git a/include/LAGraph.h b/include/LAGraph.h index 073f2bdf6c..f004194a88 100644 --- a/include/LAGraph.h +++ b/include/LAGraph.h @@ -1613,6 +1613,57 @@ int LAGraph_Matrix_Sum char *msg ) ; +//------------------------------------------------------------------------------ +// LAGraph_Matrix_Binary_Sum: sum an array of matrices by binary reduction +//------------------------------------------------------------------------------ + +/** LAGraph_Matrix_Binary_Sum: combines an array of matrices into a single + * matrix C, computing the same result as LAGraph_Matrix_Sum but with a + * different technique: a pairwise binary reduction tree built from + * GrB_eWiseAdd (as in the GraphChallenge "Anonymized Network Sensing" paper, + * "Binary Summation of Traffic Matrices"). At each level the matrices are + * summed in disjoint adjacent pairs using the binary operator dup; an unpaired + * (odd) trailing matrix is carried up to the next level unchanged, until a + * single matrix remains. The dup operator has the same set-union semantics as + * the dup of GrB_Matrix_build (applied where both inputs have an entry, lone + * entries pass through), so the result equals that of LAGraph_Matrix_Sum. + * With dup = GrB_PLUS_FP64 (for example) this computes the element-wise sum of + * all the matrices. The independent pair-sums within a level are parallelized + * across LG_nthreads_outer threads. + * + * All input matrices must have identical dimensions and identical built-in + * type; C is created with that same type and dimensions. Unlike + * LAGraph_Matrix_Sum, dup must be non-NULL, since GrB_eWiseAdd requires a + * binary operator. + * + * @param[out] C the resulting matrix. + * @param[in] Matrices array of nmatrices input matrices. + * @param[in] nmatrices number of matrices in Matrices (must be >= 1). + * @param[in] dup binary operator to combine (i,j) entries; must be + * non-NULL. + * @param[in,out] msg any error messages. + * + * @retval GrB_SUCCESS if successful. + * @retval GrB_NULL_POINTER if C, Matrices, any Matrices[k], or dup is NULL. + * @retval GrB_INVALID_VALUE if nmatrices is zero. + * @retval GrB_DIMENSION_MISMATCH if the inputs do not all share dimensions. + * @retval GrB_DOMAIN_MISMATCH if the inputs do not all share a type. + * @retval GrB_NOT_IMPLEMENTED if the matrix type is user-defined. + * @returns any GraphBLAS errors that may have been encountered. + */ + +LAGRAPH_PUBLIC +int LAGraph_Matrix_Binary_Sum +( + // output: + GrB_Matrix *C, // result = combination of all input matrices + // input: + GrB_Matrix *Matrices, // array of nmatrices input matrices + GrB_Index nmatrices, // number of matrices in the array (must be >= 1) + GrB_BinaryOp dup, // operator to combine (i,j) entries (must be != NULL) + char *msg +) ; + //------------------------------------------------------------------------------ // LAGraph_Vector_Structure: return the structure of a vector //------------------------------------------------------------------------------ diff --git a/src/test/test_Matrix_Binary_Sum.c b/src/test/test_Matrix_Binary_Sum.c new file mode 100644 index 0000000000..18063f0c6c --- /dev/null +++ b/src/test/test_Matrix_Binary_Sum.c @@ -0,0 +1,419 @@ +//------------------------------------------------------------------------------ +// LAGraph/src/test/test_Matrix_Binary_Sum.c: test LAGraph_Matrix_Binary_Sum +//------------------------------------------------------------------------------ + +// 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 Michel Pelletier. + +//------------------------------------------------------------------------------ + +#include "LAGraph_test.h" +#include "LG_internal.h" + +//------------------------------------------------------------------------------ +// global variables +//------------------------------------------------------------------------------ + +char msg [LAGRAPH_MSG_LEN] ; +GrB_Matrix A = NULL, B = NULL, C = NULL, Expected = NULL ; + +//------------------------------------------------------------------------------ +// setup and teardown +//------------------------------------------------------------------------------ + +void setup (void) +{ + OK (LAGraph_Init (msg)) ; +} + +void teardown (void) +{ + OK (LAGraph_Finalize (msg)) ; +} + +//------------------------------------------------------------------------------ +// make_AB: construct two 4x4 FP64 matrices with overlapping and distinct +// (i,j) coordinates +//------------------------------------------------------------------------------ + +static void make_AB (void) +{ + // A has entries at (0,0), (1,1), (2,2), (0,3) + GrB_Index Ai [ ] = { 0, 1, 2, 0 } ; + GrB_Index Aj [ ] = { 0, 1, 2, 3 } ; + double Ax [ ] = { 1, 2, 3, 4 } ; + OK (GrB_Matrix_new (&A, GrB_FP64, 4, 4)) ; + OK (GrB_Matrix_build_FP64 (A, Ai, Aj, Ax, 4, NULL)) ; + + // B has entries at (0,0), (1,1), (3,3), (2,0) + // shares (0,0) and (1,1) with A, distinct elsewhere + GrB_Index Bi [ ] = { 0, 1, 3, 2 } ; + GrB_Index Bj [ ] = { 0, 1, 3, 0 } ; + double Bx [ ] = { 10, 20, 30, 40 } ; + OK (GrB_Matrix_new (&B, GrB_FP64, 4, 4)) ; + OK (GrB_Matrix_build_FP64 (B, Bi, Bj, Bx, 4, NULL)) ; +} + +//------------------------------------------------------------------------------ +// test_Matrix_Binary_Sum: basic correctness +//------------------------------------------------------------------------------ + +void test_Matrix_Binary_Sum (void) +{ + setup ( ) ; + + make_AB ( ) ; + + // Expected = A + B (element-wise add, set union) + OK (GrB_Matrix_new (&Expected, GrB_FP64, 4, 4)) ; + OK (GrB_eWiseAdd (Expected, NULL, NULL, GrB_PLUS_FP64, A, B, NULL)) ; + + // C = binary sum of {A, B} using GrB_PLUS_FP64 to combine entries + GrB_Matrix Mats [2] = { A, B } ; + OK (LAGraph_Matrix_Binary_Sum (&C, Mats, 2, GrB_PLUS_FP64, msg)) ; + + bool ok ; + OK (LAGraph_Matrix_IsEqual (&ok, C, Expected, msg)) ; + TEST_CHECK (ok) ; + TEST_MSG ("binary sum of {A,B} did not equal A+B") ; + OK (GrB_free (&C)) ; + + // single-matrix array: result must be a copy of A + GrB_Matrix One [1] = { A } ; + OK (LAGraph_Matrix_Binary_Sum (&C, One, 1, GrB_PLUS_FP64, msg)) ; + OK (LAGraph_Matrix_IsEqual (&ok, C, A, msg)) ; + TEST_CHECK (ok) ; + TEST_MSG ("binary sum of {A} did not equal A") ; + OK (GrB_free (&C)) ; + + // an empty matrix in the array contributes nothing + GrB_Matrix Empty = NULL ; + OK (GrB_Matrix_new (&Empty, GrB_FP64, 4, 4)) ; + GrB_Matrix WithEmpty [3] = { A, Empty, B } ; + OK (LAGraph_Matrix_Binary_Sum (&C, WithEmpty, 3, GrB_PLUS_FP64, msg)) ; + OK (LAGraph_Matrix_IsEqual (&ok, C, Expected, msg)) ; + TEST_CHECK (ok) ; + TEST_MSG ("binary sum of {A,Empty,B} did not equal A+B") ; + OK (GrB_free (&Empty)) ; + OK (GrB_free (&C)) ; + + OK (GrB_free (&A)) ; + OK (GrB_free (&B)) ; + OK (GrB_free (&Expected)) ; + + teardown ( ) ; +} + +//------------------------------------------------------------------------------ +// test_Matrix_Binary_Sum_odd: exercise the odd-count carry path +//------------------------------------------------------------------------------ + +// Sum odd numbers of matrices (3, 5, 7, ...) so that at one or more levels an +// unpaired trailing matrix must be carried up to the next level unchanged. +// Compare against an independently accumulated expected. + +void test_Matrix_Binary_Sum_odd (void) +{ + setup ( ) ; + + GrB_Index counts [ ] = { 3, 5, 7, 9 } ; + int ncounts = sizeof (counts) / sizeof (counts [0]) ; + + for (int t = 0 ; t < ncounts ; t++) + { + GrB_Index nmat = counts [t] ; + GrB_Matrix *Mats = NULL ; + OK (LAGraph_Malloc ((void **) &Mats, nmat, sizeof (GrB_Matrix), msg)) ; + + OK (GrB_Matrix_new (&Expected, GrB_FP64, 8, 8)) ; + for (GrB_Index k = 0 ; k < nmat ; k++) + { + // 3 distinct entries per matrix; (5,5) is shared by all matrices, + // others overlap across matrices, so entries must be combined + GrB_Index Mi [ ] = { (GrB_Index) (k % 8), 1, 5 } ; + GrB_Index Mj [ ] = { 2, (GrB_Index) (k % 8), 5 } ; + double Mx [ ] = { (double) (k + 1), 1, 2 } ; + Mats [k] = NULL ; + OK (GrB_Matrix_new (&Mats [k], GrB_FP64, 8, 8)) ; + OK (GrB_Matrix_build_FP64 (Mats [k], Mi, Mj, Mx, 3, NULL)) ; + OK (GrB_eWiseAdd (Expected, NULL, NULL, GrB_PLUS_FP64, Expected, + Mats [k], NULL)) ; + } + + OK (LAGraph_Matrix_Binary_Sum (&C, Mats, nmat, GrB_PLUS_FP64, msg)) ; + + bool ok ; + OK (LAGraph_Matrix_IsEqual (&ok, C, Expected, msg)) ; + TEST_CHECK (ok) ; + TEST_MSG ("binary sum of %d matrices did not match expected", + (int) nmat) ; + + for (GrB_Index k = 0 ; k < nmat ; k++) + { + OK (GrB_free (&Mats [k])) ; + } + OK (LAGraph_Free ((void **) &Mats, msg)) ; + OK (GrB_free (&C)) ; + OK (GrB_free (&Expected)) ; + } + + teardown ( ) ; +} + +//------------------------------------------------------------------------------ +// test_Matrix_Binary_Sum_types: exercise every built-in type branch +//------------------------------------------------------------------------------ + +void test_Matrix_Binary_Sum_types (void) +{ + setup ( ) ; + + // (type, GrB_PLUS operator) pairs, one per built-in type branch + GrB_Type types [ ] = { GrB_BOOL, GrB_INT8, GrB_INT16, GrB_INT32, + GrB_INT64, GrB_UINT8, GrB_UINT16, GrB_UINT32, GrB_UINT64, + GrB_FP32, GrB_FP64 } ; + GrB_BinaryOp ops [ ] = { GrB_LOR, GrB_PLUS_INT8, GrB_PLUS_INT16, + GrB_PLUS_INT32, GrB_PLUS_INT64, GrB_PLUS_UINT8, GrB_PLUS_UINT16, + GrB_PLUS_UINT32, GrB_PLUS_UINT64, GrB_PLUS_FP32, GrB_PLUS_FP64 } ; + int ntypes = sizeof (types) / sizeof (types [0]) ; + + GrB_Index Ai [ ] = { 0, 1, 2 } ; + GrB_Index Aj [ ] = { 0, 1, 0 } ; + GrB_Index Bi [ ] = { 0, 1, 0 } ; + GrB_Index Bj [ ] = { 0, 1, 2 } ; + + for (int t = 0 ; t < ntypes ; t++) + { + OK (GrB_Matrix_new (&A, types [t], 3, 3)) ; + OK (GrB_Matrix_new (&B, types [t], 3, 3)) ; + // build with INT32 values; GraphBLAS typecasts to the matrix type + int32_t Ax [ ] = { 1, 1, 1 } ; + int32_t Bx [ ] = { 1, 1, 1 } ; + OK (GrB_Matrix_build_INT32 (A, Ai, Aj, Ax, 3, NULL)) ; + OK (GrB_Matrix_build_INT32 (B, Bi, Bj, Bx, 3, NULL)) ; + + OK (GrB_Matrix_new (&Expected, types [t], 3, 3)) ; + OK (GrB_eWiseAdd (Expected, NULL, NULL, ops [t], A, B, NULL)) ; + + GrB_Matrix Mats [2] = { A, B } ; + OK (LAGraph_Matrix_Binary_Sum (&C, Mats, 2, ops [t], msg)) ; + + bool ok ; + OK (LAGraph_Matrix_IsEqual (&ok, C, Expected, msg)) ; + TEST_CHECK (ok) ; + TEST_MSG ("type index %d failed", t) ; + + OK (GrB_free (&A)) ; + OK (GrB_free (&B)) ; + OK (GrB_free (&C)) ; + OK (GrB_free (&Expected)) ; + } + + teardown ( ) ; +} + +//------------------------------------------------------------------------------ +// test_Matrix_Binary_Sum_brutal +//------------------------------------------------------------------------------ + +#if LG_BRUTAL_TESTS +void test_Matrix_Binary_Sum_brutal (void) +{ + OK (LG_brutal_setup (msg)) ; + + // five matrices: an odd count, so the carry path is exercised under brutal + // malloc-failure testing as well, validating the Pool / W / Wnext cleanup + #define NBRUTAL 5 + GrB_Matrix Mats [NBRUTAL] ; + OK (GrB_Matrix_new (&Expected, GrB_FP64, 6, 6)) ; + for (int k = 0 ; k < NBRUTAL ; k++) + { + GrB_Index Mi [ ] = { (GrB_Index) (k % 6), 1, 4 } ; + GrB_Index Mj [ ] = { 2, (GrB_Index) (k % 6), 4 } ; + double Mx [ ] = { (double) (k + 1), 1, 2 } ; + Mats [k] = NULL ; + OK (GrB_Matrix_new (&Mats [k], GrB_FP64, 6, 6)) ; + OK (GrB_Matrix_build_FP64 (Mats [k], Mi, Mj, Mx, 3, NULL)) ; + OK (GrB_eWiseAdd (Expected, NULL, NULL, GrB_PLUS_FP64, Expected, + Mats [k], NULL)) ; + } + + LG_BRUTAL (LAGraph_Matrix_Binary_Sum (&C, Mats, NBRUTAL, GrB_PLUS_FP64, + msg)) ; + + bool ok ; + OK (LAGraph_Matrix_IsEqual (&ok, C, Expected, msg)) ; + TEST_CHECK (ok) ; + + for (int k = 0 ; k < NBRUTAL ; k++) + { + OK (GrB_free (&Mats [k])) ; + } + OK (GrB_free (&C)) ; + OK (GrB_free (&Expected)) ; + #undef NBRUTAL + + OK (LG_brutal_teardown (msg)) ; +} +#endif + +//------------------------------------------------------------------------------ +// test_Matrix_Binary_Sum_parallel: exercise the parallel pair-sum path +//------------------------------------------------------------------------------ + +// Sum many matrices with multiple outer threads and confirm the result matches +// an independently accumulated expected. This stresses the concurrent +// pair-sums of the binary reduction under real outer parallelism, including an +// odd count to also drive the carry path. + +void test_Matrix_Binary_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)) ; + + GrB_Index counts [ ] = { 16, 15 } ; // even and odd + int ncounts = sizeof (counts) / sizeof (counts [0]) ; + + for (int t = 0 ; t < ncounts ; t++) + { + GrB_Index nmat = counts [t] ; + GrB_Matrix *Mats = NULL ; + OK (LAGraph_Malloc ((void **) &Mats, nmat, sizeof (GrB_Matrix), msg)) ; + + OK (GrB_Matrix_new (&Expected, GrB_FP64, 10, 10)) ; + for (GrB_Index 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 combined when the matrices are summed + 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_Binary_Sum (&C, Mats, nmat, GrB_PLUS_FP64, msg)) ; + + bool ok ; + OK (LAGraph_Matrix_IsEqual (&ok, C, Expected, msg)) ; + TEST_CHECK (ok) ; + TEST_MSG ("parallel binary sum of %d matrices did not match expected", + (int) nmat) ; + + for (GrB_Index k = 0 ; k < nmat ; k++) + { + OK (GrB_free (&Mats [k])) ; + } + OK (LAGraph_Free ((void **) &Mats, msg)) ; + OK (GrB_free (&C)) ; + OK (GrB_free (&Expected)) ; + } + + // restore the original thread settings + OK (LAGraph_SetNumThreads (save_outer, save_inner, msg)) ; + + teardown ( ) ; +} + +//------------------------------------------------------------------------------ +// test_Matrix_Binary_Sum_failures: test error handling +//------------------------------------------------------------------------------ + +void test_Matrix_Binary_Sum_failures (void) +{ + setup ( ) ; + + make_AB ( ) ; + GrB_Matrix Mats [2] = { A, B } ; + int result ; + + // NULL output + result = LAGraph_Matrix_Binary_Sum (NULL, Mats, 2, GrB_PLUS_FP64, msg) ; + TEST_CHECK (result == GrB_NULL_POINTER) ; + + // NULL Matrices array + C = NULL ; + result = LAGraph_Matrix_Binary_Sum (&C, NULL, 2, GrB_PLUS_FP64, msg) ; + TEST_CHECK (result == GrB_NULL_POINTER) ; + TEST_CHECK (C == NULL) ; + + // nmatrices == 0 + result = LAGraph_Matrix_Binary_Sum (&C, Mats, 0, GrB_PLUS_FP64, msg) ; + TEST_CHECK (result == GrB_INVALID_VALUE) ; + + // NULL dup operator + result = LAGraph_Matrix_Binary_Sum (&C, Mats, 2, NULL, msg) ; + TEST_CHECK (result == GrB_NULL_POINTER) ; + + // NULL entry in the array + GrB_Matrix WithNull [2] = { A, NULL } ; + result = LAGraph_Matrix_Binary_Sum (&C, WithNull, 2, GrB_PLUS_FP64, msg) ; + TEST_CHECK (result == GrB_NULL_POINTER) ; + + // dimension mismatch + GrB_Matrix D = NULL ; + OK (GrB_Matrix_new (&D, GrB_FP64, 5, 5)) ; + GrB_Matrix BadDim [2] = { A, D } ; + result = LAGraph_Matrix_Binary_Sum (&C, BadDim, 2, GrB_PLUS_FP64, msg) ; + TEST_CHECK (result == GrB_DIMENSION_MISMATCH) ; + OK (GrB_free (&D)) ; + + // type mismatch + GrB_Matrix E = NULL ; + OK (GrB_Matrix_new (&E, GrB_INT32, 4, 4)) ; + GrB_Matrix BadType [2] = { A, E } ; + result = LAGraph_Matrix_Binary_Sum (&C, BadType, 2, GrB_PLUS_FP64, msg) ; + TEST_CHECK (result == GrB_DOMAIN_MISMATCH) ; + OK (GrB_free (&E)) ; + + // user-defined type not supported + typedef struct { double x, y ; } udt_t ; + GrB_Type UDT = NULL ; + OK (GrB_Type_new (&UDT, sizeof (udt_t))) ; + GrB_Matrix U = NULL ; + OK (GrB_Matrix_new (&U, UDT, 4, 4)) ; + GrB_Matrix Udt [1] = { U } ; + result = LAGraph_Matrix_Binary_Sum (&C, Udt, 1, GrB_PLUS_FP64, msg) ; + TEST_CHECK (result == GrB_NOT_IMPLEMENTED) ; + OK (GrB_free (&U)) ; + OK (GrB_free (&UDT)) ; + + OK (GrB_free (&A)) ; + OK (GrB_free (&B)) ; + + teardown ( ) ; +} + +//------------------------------------------------------------------------------ +// TEST_LIST: the list of tasks for this entire test +//------------------------------------------------------------------------------ + +TEST_LIST = +{ + { "Matrix_Binary_Sum", test_Matrix_Binary_Sum }, + { "Matrix_Binary_Sum_odd", test_Matrix_Binary_Sum_odd }, + { "Matrix_Binary_Sum_types", test_Matrix_Binary_Sum_types }, + { "Matrix_Binary_Sum_parallel", test_Matrix_Binary_Sum_parallel }, + { "Matrix_Binary_Sum_failures", test_Matrix_Binary_Sum_failures }, + #if LG_BRUTAL_TESTS + { "Matrix_Binary_Sum_brutal", test_Matrix_Binary_Sum_brutal }, + #endif + { NULL, NULL } +} ; diff --git a/src/utility/LAGraph_Matrix_Binary_Sum.c b/src/utility/LAGraph_Matrix_Binary_Sum.c new file mode 100644 index 0000000000..d89293f047 --- /dev/null +++ b/src/utility/LAGraph_Matrix_Binary_Sum.c @@ -0,0 +1,272 @@ +//------------------------------------------------------------------------------ +// LAGraph_Matrix_Binary_Sum: sum an array of matrices by binary reduction +//------------------------------------------------------------------------------ + +// 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 Michel Pelletier. + +//------------------------------------------------------------------------------ + +// LAGraph_Matrix_Binary_Sum combines an array of matrices into a single matrix +// C, computing the same result as LAGraph_Matrix_Sum but with a different +// technique: a pairwise binary reduction tree built from GrB_eWiseAdd, as +// described in the GraphChallenge "Anonymized Network Sensing" paper (Fig. 3, +// "Binary Summation of Traffic Matrices"). At each level of the tree the +// matrices are summed in disjoint adjacent pairs; an unpaired (odd) trailing +// matrix is carried up to the next level unchanged. The levels repeat until a +// single matrix remains. Each pair-sum uses the binary operator dup, which has +// the same set-union semantics as the dup operator of GrB_Matrix_build: dup is +// applied wherever both inputs have an entry, and lone entries pass through. +// With dup = GrB_PLUS_FP64 (for example) this computes the element-wise sum of +// all the input matrices. + +// The independent pair-sums within a level are issued concurrently across +// LG_nthreads_outer threads with OpenMP; each such outer thread calls into +// GraphBLAS, which parallelizes each GrB_eWiseAdd internally with +// LG_nthreads_inner threads (the documented two-level model). Working on the +// smaller intermediate matrices of the tree, rather than concatenating every +// tuple at once, keeps the active data in faster memory: adding two matrices +// each with N entries yields a matrix with fewer than 2N entries, so the +// relative work shrinks as the matrices grow. + +// All input matrices must have identical dimensions and identical built-in +// type; C is created with that same type and dimensions. Unlike +// LAGraph_Matrix_Sum, the dup operator must be non-NULL, since GrB_eWiseAdd +// requires a binary operator. + +// Every intermediate matrix created by the reduction is held in a single Pool +// array (a binary reduction of n leaves has exactly n-1 internal sums), so the +// free-path is a simple sweep over Pool that is correct at any point, including +// partial failure inside a level (GrB_free of a NULL handle is a no-op). The W +// and Wnext arrays hold only pointers (to original inputs or to Pool entries) +// and are therefore never themselves freed as matrices. + +#define LG_FREE_WORK \ +{ \ + if (Pool != NULL) \ + { \ + for (GrB_Index k = 0 ; k < npool ; k++) \ + { \ + GrB_free (&Pool [k]) ; \ + } \ + } \ + LAGraph_Free ((void **) &Pool, NULL) ; \ + LAGraph_Free ((void **) &W, NULL) ; \ + LAGraph_Free ((void **) &Wnext, NULL) ; \ +} + +#define LG_FREE_ALL \ +{ \ + LG_FREE_WORK ; \ + GrB_free (C) ; \ +} + +#include "LG_internal.h" + +int LAGraph_Matrix_Binary_Sum +( + // output: + GrB_Matrix *C, // result = combination of all input matrices + // input: + GrB_Matrix *Matrices, // array of nmatrices input matrices + GrB_Index nmatrices, // number of matrices in the array (must be >= 1) + GrB_BinaryOp dup, // operator to combine (i,j) entries (must be != NULL) + char *msg +) +{ + + //-------------------------------------------------------------------------- + // check inputs + //-------------------------------------------------------------------------- + + LG_CLEAR_MSG ; + GrB_Matrix *W = NULL, *Wnext = NULL, *Pool = NULL ; + GrB_Index npool = 0 ; + LG_ASSERT_MSG (C != NULL, GrB_NULL_POINTER, "&C != NULL") ; + LG_ASSERT (Matrices != NULL, GrB_NULL_POINTER) ; + (*C) = NULL ; + LG_ASSERT_MSG (nmatrices >= 1, GrB_INVALID_VALUE, + "nmatrices must be >= 1") ; + LG_ASSERT_MSG (dup != NULL, GrB_NULL_POINTER, + "dup operator must be non-NULL") ; + LG_ASSERT (Matrices [0] != NULL, GrB_NULL_POINTER) ; + + //-------------------------------------------------------------------------- + // determine the reference dimensions and type from the first matrix + //-------------------------------------------------------------------------- + + GrB_Index nrows, ncols ; + int32_t typecode ; + GRB_TRY (GrB_Matrix_nrows (&nrows, Matrices [0])) ; + GRB_TRY (GrB_Matrix_ncols (&ncols, Matrices [0])) ; + GRB_TRY (GrB_get (Matrices [0], &typecode, GrB_EL_TYPE_CODE)) ; + + // map the type code to a built-in GrB_Type for the intermediate matrices; + // GrB_eWiseAdd itself is polymorphic, so only GrB_Matrix_new needs the type + GrB_Type gtype = NULL ; + switch (typecode) + { + case GrB_BOOL_CODE : gtype = GrB_BOOL ; break ; + case GrB_INT8_CODE : gtype = GrB_INT8 ; break ; + case GrB_INT16_CODE : gtype = GrB_INT16 ; break ; + case GrB_INT32_CODE : gtype = GrB_INT32 ; break ; + case GrB_INT64_CODE : gtype = GrB_INT64 ; break ; + case GrB_UINT8_CODE : gtype = GrB_UINT8 ; break ; + case GrB_UINT16_CODE : gtype = GrB_UINT16 ; break ; + case GrB_UINT32_CODE : gtype = GrB_UINT32 ; break ; + case GrB_UINT64_CODE : gtype = GrB_UINT64 ; break ; + case GrB_FP32_CODE : gtype = GrB_FP32 ; break ; + case GrB_FP64_CODE : gtype = GrB_FP64 ; break ; + default : + LG_ASSERT_MSG (false, GrB_NOT_IMPLEMENTED, + "user-defined types are not supported") ; + } + + //-------------------------------------------------------------------------- + // validate every matrix has the same dimensions and type + //-------------------------------------------------------------------------- + + for (GrB_Index k = 0 ; k < nmatrices ; k++) + { + GrB_Matrix Ak = Matrices [k] ; + LG_ASSERT (Ak != NULL, GrB_NULL_POINTER) ; + GrB_Index r, c ; + int32_t code ; + GRB_TRY (GrB_Matrix_nrows (&r, Ak)) ; + GRB_TRY (GrB_Matrix_ncols (&c, Ak)) ; + LG_ASSERT_MSG (r == nrows && c == ncols, GrB_DIMENSION_MISMATCH, + "all input matrices must have the same dimensions") ; + GRB_TRY (GrB_get (Ak, &code, GrB_EL_TYPE_CODE)) ; + LG_ASSERT_MSG (code == typecode, GrB_DOMAIN_MISMATCH, + "all input matrices must have the same type") ; + } + + //-------------------------------------------------------------------------- + // handle the single-matrix case: C is an independent copy of Matrices [0] + //-------------------------------------------------------------------------- + + if (nmatrices == 1) + { + GRB_TRY (GrB_Matrix_dup (C, Matrices [0])) ; + return (GrB_SUCCESS) ; + } + + //-------------------------------------------------------------------------- + // allocate the working pointer arrays and the pool of intermediate matrices + //-------------------------------------------------------------------------- + + // A binary reduction of nmatrices leaves creates exactly nmatrices-1 sums. + // Pool holds every one of them and is pre-filled with NULL, so LG_FREE_WORK + // can sweep it at any time (including on partial failure within a level). + // Each pair at a level writes a deterministic Pool slot (base + p), so no + // shared counter is touched by the parallel loop. + + npool = nmatrices - 1 ; + GrB_Index pool_alloc = (npool == 0) ? 1 : npool ; + LG_TRY (LAGraph_Malloc ((void **) &Pool, pool_alloc, sizeof (GrB_Matrix), + msg)) ; + for (GrB_Index k = 0 ; k < pool_alloc ; k++) + { + Pool [k] = NULL ; + } + LG_TRY (LAGraph_Malloc ((void **) &W, nmatrices, sizeof (GrB_Matrix), + msg)) ; + LG_TRY (LAGraph_Malloc ((void **) &Wnext, nmatrices, sizeof (GrB_Matrix), + msg)) ; + + // level 0 references the original input matrices (which we never free) + for (GrB_Index k = 0 ; k < nmatrices ; k++) + { + W [k] = Matrices [k] ; + } + + //-------------------------------------------------------------------------- + // pairwise binary reduction + //-------------------------------------------------------------------------- + + GrB_Index count = nmatrices ; // number of matrices at the current level + GrB_Index base = 0 ; // first Pool slot used by the current level + + while (count > 1) + { + GrB_Index npairs = count / 2 ; // number of pair-sums + int64_t p ; + + // outer threads for the independent pair-sums at this level; each calls + // GraphBLAS, which nests inner threads underneath (two-level model) + int nthreads = LG_nthreads_outer ; + nthreads = LAGRAPH_MIN (nthreads, (int) npairs) ; + nthreads = LAGRAPH_MAX (nthreads, 1) ; + + // GRB_TRY cannot be used inside the 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. + int sum_status = GrB_SUCCESS ; + + #pragma omp parallel for num_threads(nthreads) schedule(dynamic,1) + for (p = 0 ; p < (int64_t) npairs ; p++) + { + GrB_Matrix R = NULL ; + GrB_Info info = GrB_Matrix_new (&R, gtype, nrows, ncols) ; + if (info >= GrB_SUCCESS) + { + info = GrB_eWiseAdd (R, NULL, NULL, dup, + W [2*p], W [2*p+1], NULL) ; + } + // each p writes a distinct Pool slot and a distinct Wnext slot + Pool [base + p] = R ; + Wnext [p] = R ; + if (info < GrB_SUCCESS) + { + #pragma omp critical + { + if (sum_status >= GrB_SUCCESS) sum_status = info ; + } + } + } + GRB_TRY (sum_status) ; + + // carry an unpaired (odd) trailing matrix up to the next level + if (count % 2 == 1) + { + Wnext [npairs] = W [count-1] ; + } + + // advance to the next level: swap W and Wnext, update count and base + GrB_Matrix *tmp = W ; W = Wnext ; Wnext = tmp ; + base += npairs ; + count = (count + 1) / 2 ; + } + + //-------------------------------------------------------------------------- + // hand the root of the tree to the caller + //-------------------------------------------------------------------------- + + // With nmatrices >= 2 the final root is always a genuine sum (a Pool entry): + // count==1 is only ever reached from a count==2 level, which has no carry. + // Transfer ownership by clearing its Pool slot so LG_FREE_WORK won't free it. + (*C) = W [0] ; + for (GrB_Index k = 0 ; k < npool ; k++) + { + if (Pool [k] == (*C)) + { + Pool [k] = NULL ; + break ; + } + } + + //-------------------------------------------------------------------------- + // free workspace and return result + //-------------------------------------------------------------------------- + + LG_FREE_WORK ; + return (GrB_SUCCESS) ; +}