diff --git a/Config/LAGraph.h.in b/Config/LAGraph.h.in index a46c06c8eb..b5e6320569 100644 --- a/Config/LAGraph.h.in +++ b/Config/LAGraph.h.in @@ -1572,6 +1572,98 @@ 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_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 ea88c7c7dc..f004194a88 100644 --- a/include/LAGraph.h +++ b/include/LAGraph.h @@ -1572,6 +1572,98 @@ 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_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/test/test_Matrix_Sum.c b/src/test/test_Matrix_Sum.c new file mode 100644 index 0000000000..3863ced258 --- /dev/null +++ b/src/test/test_Matrix_Sum.c @@ -0,0 +1,333 @@ +//------------------------------------------------------------------------------ +// 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_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 +//------------------------------------------------------------------------------ + +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_parallel", test_Matrix_Sum_parallel }, + { "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_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) ; +} diff --git a/src/utility/LAGraph_Matrix_Sum.c b/src/utility/LAGraph_Matrix_Sum.c new file mode 100644 index 0000000000..26dbf5f773 --- /dev/null +++ b/src/utility/LAGraph_Matrix_Sum.c @@ -0,0 +1,199 @@ +//------------------------------------------------------------------------------ +// 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 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. + +// 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) ; \ + LAGraph_Free ((void **) &Offsets, 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, *Offsets = 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 compute where its tuples begin in the buffer + //-------------------------------------------------------------------------- + + // 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] ; + 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)) ; + Offsets [k+1] = Offsets [k] + n ; + } + GrB_Index total = Offsets [nmatrices] ; + + //-------------------------------------------------------------------------- + // 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)) ; + + //-------------------------------------------------------------------------- + // 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 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 : \ + { \ + ctype *Xt = NULL ; \ + LG_TRY (LAGraph_Malloc ((void **) &Xt, alloc, sizeof (ctype), \ + msg)) ; \ + X = (void *) Xt ; \ + 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 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)) ; \ + } \ + 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) ; +}