From eff687b7bbb8a68239730401d3e9dccd14ce6d61 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Fri, 26 Jan 2024 13:47:44 +0100 Subject: [PATCH 01/51] Masked outer product implementation --- include/graphblas/reference/blas3.hpp | 200 ++++++++++++++++++++++++++ 1 file changed, 200 insertions(+) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 868547183..4baa2965a 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -916,6 +916,206 @@ namespace grb { return ret; } + + /** + * A masked outer product of two vectors. Assuming vectors \a u and \a v are oriented + * column-wise, the result matrix \a A will contain \f$ uv^T \f$, masked to non-zero values from \a mask. + */ + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType1, typename InputType2, typename OutputType, + typename Coords, + typename RIT, typename CIT, typename NIT + > + RC maskedOuter( + Matrix< OutputType, reference, RIT, CIT, NIT > &A, + const Matrix< OutputType, reference, RIT, CIT, NIT > &mask, + const Vector< InputType1, reference, Coords > &u, + const Vector< InputType2, reference, Coords > &v, + const Operator &mul = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + // static checks + NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || + std::is_same< typename Operator::D1, InputType1 >::value + ), "grb::maskedOuter", + "called with a prefactor vector that does not match the first domain " + "of the given multiplication operator" ); + NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || + std::is_same< typename Operator::D2, InputType2 >::value + ), "grb::maskedOuter", + "called with a postfactor vector that does not match the first domain " + "of the given multiplication operator" ); + NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || + std::is_same< typename Operator::D3, OutputType >::value + ), "grb::maskedOuter", + "called with an output matrix that does not match the output domain of " + "the given multiplication operator" ); + static_assert( ( !(descr & descriptors::invert_mask) + ), + "grb::maskedOuter: invert_mask descriptor cannot be used "); +#ifdef _DEBUG + std::cout << "In grb::maskedOuter (reference)\n"; +#endif + + const size_t nrows = size( u ); + const size_t ncols = size( v ); + + const size_t m = grb::nrows( mask ); + const size_t n = grb::ncols( mask ); + + assert( phase != TRY ); + if( nrows != grb::nrows( A ) || nrows != m ) { + return MISMATCH; + } + + if( ncols != grb::ncols( A ) || ncols != n ) { + return MISMATCH; + } + + + const auto &mask_raw = internal::getCRS( mask ); + + size_t nzc = 0; + for( size_t i = 0; i < nrows; ++i ) { + if( internal::getCoordinates( u ).assigned( i ) ) { + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + if( internal::getCoordinates( v ).assigned( k_col ) ) { + nzc++; + } + } + } + } + + if( phase == RESIZE ) { + return resize( A, nzc ); + } + + assert( phase == EXECUTE ); + if( capacity( A ) < nzc ) { +#ifdef _DEBUG + std::cout << "\t insufficient capacity to complete " + "requested masked outer-product computation\n"; +#endif + const RC clear_rc = clear( A ); + if( clear_rc != SUCCESS ) { + return PANIC; + } else { + return FAILED; + } + } + + grb::Monoid< + grb::operators::left_assign< OutputType >, + grb::identities::zero + > monoid; + + RC ret = SUCCESS; + if( phase == EXECUTE ) { + ret = grb::clear( A ); + } + assert( nnz( A ) == 0 ); + + auto &CRS_raw = internal::getCRS( A ); + auto &CCS_raw = internal::getCCS( A ); + + const InputType1 * __restrict__ const x = internal::getRaw( u ); + const InputType2 * __restrict__ const y = internal::getRaw( v ); + char * arr = nullptr; + char * buf = nullptr; + OutputType * valbuf = nullptr; + internal::getMatrixBuffers( arr, buf, valbuf, 1, A ); + config::NonzeroIndexType * A_col_index = internal::template + getReferenceBuffer< typename config::NonzeroIndexType >( ncols + 1 ); + + + internal::Coordinates< reference > coors; + coors.set( arr, false, buf, ncols ); + + CRS_raw.col_start[ 0 ] = 0; + for( size_t j = 0; j <= ncols; ++j ) { + CCS_raw.col_start[ j ] = 0; + } + + + nzc = 0; + for( size_t i = 0; i < nrows; ++i ) { + if( internal::getCoordinates( u ).assigned( i ) ) { + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + if( internal::getCoordinates( v ).assigned( k_col ) ) { + (void) ++nzc; + (void) ++CCS_raw.col_start[ k_col + 1 ]; + } + } + } + CRS_raw.col_start[ i + 1 ] = nzc; + } + + + A_col_index[ 0 ] = 0; + for( size_t j = 1; j < ncols; ++j ) { + CCS_raw.col_start[ j + 1 ] += CCS_raw.col_start[ j ]; + A_col_index[ j ] = 0; + } + + + const size_t old_nzc = nzc; + + // use previously computed CCS offset array to update CCS during the + // computational phase + nzc = 0; + for( size_t i = 0; i < nrows; ++i ) { + coors.clear(); + if( internal::getCoordinates( u ).assigned( i ) ) { + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + if( internal::getCoordinates( v ).assigned( k_col ) ) { + coors.assign( k_col ); + valbuf[ k_col ] = monoid.template getIdentity< OutputType >(); + (void) grb::apply( valbuf[ k_col ], + x[i], + y[k_col], + mul ); + } + } + } + for( size_t k = 0; k < coors.nonzeroes(); ++k ) { + assert( nzc < old_nzc ); + const size_t j = coors.index( k ); + // update CRS + CRS_raw.row_index[ nzc ] = j; + CRS_raw.setValue( nzc, valbuf[ j ] ); + // update CCS + const size_t CCS_index = A_col_index[ j ]++ + CCS_raw.col_start[ j ]; + CCS_raw.row_index[ CCS_index ] = i; + CCS_raw.setValue( CCS_index, valbuf[ j ] ); + // update count + (void) ++nzc; + } + CRS_raw.col_start[ i + 1 ] = nzc; + } + + for( size_t j = 0; j < ncols; ++j ) { + assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] == + A_col_index[ j ] ); + } + assert( nzc == old_nzc ); + + + internal::setCurrentNonzeroes( A, nzc ); + + return ret; + } + namespace internal { /** From 1a9663b49b3a2079c6ab9343a52dad76483dc83d Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Fri, 26 Jan 2024 14:56:01 +0100 Subject: [PATCH 02/51] Add unit test for the masked outer product --- tests/unit/CMakeLists.txt | 4 + tests/unit/maskedOuter.cpp | 226 +++++++++++++++++++++++++++++++++++++ tests/unit/unittests.sh | 6 + 3 files changed, 236 insertions(+) create mode 100644 tests/unit/maskedOuter.cpp diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 8efb144ba..7191679b2 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -246,6 +246,10 @@ add_grb_executables( outer outer.cpp BACKENDS reference reference_omp hyperdags nonblocking ) +add_grb_executables( maskedOuter maskedOuter.cpp + BACKENDS reference +) + # must generate the golden output for other tests force_add_grb_executable( mxv mxv.cpp BACKEND reference diff --git a/tests/unit/maskedOuter.cpp b/tests/unit/maskedOuter.cpp new file mode 100644 index 000000000..e0df449d2 --- /dev/null +++ b/tests/unit/maskedOuter.cpp @@ -0,0 +1,226 @@ + +/* + * Copyright 2021 Huawei Technologies Co., Ltd. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +#include "graphblas.hpp" + + +using namespace grb; + +// sample data +static const double vec1_vals[ 3 ] = { 1, 2, 3 }; +static const double vec2_vals[ 3 ] = { 4, 5, 6 }; + +static const double m1[ 3 ] = { 0.5, 3.4, 5 }; + +static const size_t I1[ 3 ] = { 0, 1, 2 }; +static const size_t J1[ 3 ] = { 0, 1, 2 }; + +static const double m2[ 6 ] = { 0.5, -1, 23, 34, -4.4, 3.14 }; + +static const size_t I2[ 6 ] = { 0, 0, 1, 1, 2, 2 }; +static const size_t J2[ 6 ] = { 1, 2, 0, 2, 0, 1 }; + +static const double mask_test1_in[ 3 ] = { 1, 1, 1 }; +static const double mask_test1_expect[ 3 ] = { 4, 10, 18 }; + +static const double mask_test2_in[ 3 ] = { 1, 1, 1 }; +static const double mask_test2_expect[ 3 ] = { 11, 20, 27 }; + +void grbProgram( const void *, const size_t in_size, int &error ) { + error = 0; + + if( in_size != 0 ) { + (void)fprintf( stderr, "Unit tests called with unexpected input\n" ); + error = 1; + return; + } + + // allocate + grb::Vector< double > u( 3 ); + grb::Vector< double > v( 3 ); + grb::Matrix< double > Result1( 3, 3 ); + grb::Matrix< double > Result2( 3, 3 ); + grb::Matrix< double > mask1( 3, 3 ); + grb::Matrix< double > mask2( 3, 3 ); + grb::Vector< double > test1( 3 ); + grb::Vector< double > out1( 3 ); + grb::Vector< double > test2( 3 ); + grb::Vector< double > out2( 3 ); + + // semiring + grb::Semiring< + grb::operators::add< double >, grb::operators::mul< double >, + grb::identities::zero, grb::identities::one + > ring; + + // initialise vec + const double * vec_iter = &(vec1_vals[ 0 ]); + grb::RC rc = grb::buildVector( u, vec_iter, vec_iter + 3, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t initial buildVector FAILED\n"; + error = 5; + } + + if( !error ) { + vec_iter = &(vec2_vals[ 0 ]); + rc = grb::buildVector( v, vec_iter, vec_iter + 3, SEQUENTIAL ); + } + if( rc != SUCCESS ) { + std::cerr << "\t initial buildVector FAILED\n"; + error = 10; + } + + if( !error ) { + rc = grb::buildMatrixUnique( mask1, &( I1[ 0 ] ), &( J1[ 0 ] ), m1, 3, SEQUENTIAL ); + } + + if( rc != SUCCESS ) { + std::cerr << "\t first buildMatrix FAILED\n"; + error = 15; + } + + if( !error ) { + rc = grb::buildMatrixUnique( mask2, &( I2[ 0 ] ), &( J2[ 0 ] ), m2, 6, SEQUENTIAL ); + } + + if( rc != SUCCESS ) { + std::cerr << "\t second buildMatrix FAILED\n"; + error = 20; + } + + if( !error ) { + rc = grb::maskedOuter( Result1, mask1, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::maskedOuter( Result1, mask1, u, v, ring.getMultiplicativeOperator() ); + } + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::maskedOuter: " + << toString( rc ) << ".\n"; + error = 25; + } + + if( !error && grb::nnz( Result1 ) != 3 ) { + std::cerr << "\t Unexpected number of nonzeroes in matrix: " + << grb::nnz( Result1 ) << ", expected 3.\n"; + error = 30; + } + + if( !error ) { + rc = grb::maskedOuter( Result2, mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::maskedOuter( Result2, mask2, u, v, ring.getMultiplicativeOperator() ); + } + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::maskedOuter: " + << toString( rc ) << ".\n"; + error = 35; + } + + if( !error && grb::nnz( Result2 ) != 6 ) { + std::cerr << "\t Unexpected number of nonzeroes in matrix: " + << grb::nnz( Result2 ) << ", expected 6.\n"; + error = 40; + } + + if( !error ) { + const double * const test1_iter = &( mask_test1_in[ 0 ] ); + rc = grb::buildVector( test1, test1_iter, test1_iter + 3, SEQUENTIAL ); + if( rc == grb::SUCCESS ) { + rc = grb::vxm( out1, test1, Result1, ring ); + } + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from premultiplying M by a vector (vxm): " + << toString( rc ) << ".\n"; + error = 45; + } + } + + if( !error ) { + if( grb::nnz( out1 ) != 3 ) { + std::cerr << "\t Unexpected number of nonzeroes (premultiply): " + << grb::nnz( out1 ) << ", expected 3\n"; + error = 50; + } + for( const auto &pair : out1 ) { + size_t i = pair.first; + if( pair.second != mask_test1_expect[ i ] ) { + std::cerr << "Premultiplying Result1 by a vector of all ones, " + << "unexpected value " << pair.second << " " + << "at coordinate " << i << ", expected " + << mask_test1_expect[ i ] << ".\n"; + error = 55; + break; + } + } + } + + if( !error ) { + const double * const test2_iter = &( mask_test2_in[ 0 ] ); + rc = grb::buildVector( test2, test2_iter, test2_iter + 3, SEQUENTIAL ); + if( rc == grb::SUCCESS ) { + rc = grb::vxm< grb::descriptors::transpose_matrix >( out2, test2, Result2, ring ); + } + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from postmultiplying M by a vector (vxm): " + << toString( rc ) << ".\n"; + error = 60; + } + } + + if( !error ) { + if( grb::nnz( out2 ) != 3 ) { + std::cerr << "\t Unexpected number of nonzeroes (postmultiply): " + << grb::nnz( out1 ) << ", expected 3\n"; + error = 65; + } + for( const auto &pair : out2 ) { + size_t i = pair.first; + if( pair.second != mask_test2_expect[ i ] ) { + std::cerr << "Postmultiplying M by a vector of all ones, " + << "unexpected value " << pair.second << " " + << "at coordinate " << i << ", " + << "expected " << mask_test2_expect[ i ] << ".\n"; + error = 70; + break; + } + } + } +} + +int main( int argc, char ** argv ) { + (void)argc; + std::cout << "Functional test executable: " << argv[ 0 ] << "\n"; + + int error; + grb::Launcher< AUTOMATIC > launcher; + + if( launcher.exec( &grbProgram, nullptr, 0, error ) != SUCCESS ) { + std::cerr << "Test failed to launch\n"; + error = 255; + } + if( error == 0 ) { + std::cout << "Test OK\n" << std::endl; + } else { + std::cerr << std::flush; + std::cout << "Test FAILED\n" << std::endl; + } + + // done + return error; +} + diff --git a/tests/unit/unittests.sh b/tests/unit/unittests.sh index ae167c0b9..29c4e4752 100755 --- a/tests/unit/unittests.sh +++ b/tests/unit/unittests.sh @@ -587,6 +587,12 @@ for MODE in ${MODES}; do grep 'Test OK' ${TEST_OUT_DIR}/outer_${MODE}_${BACKEND}_${P}_${T}.log || echo "Test FAILED" echo " " + echo ">>> [x] [ ] Testing grb::maskedOuter on a small matrix" + $runner ${TEST_BIN_DIR}/maskedOuter_${MODE}_${BACKEND} &> ${TEST_OUT_DIR}/maskedOuter_${MODE}_${BACKEND}_${P}_${T}.log + head -1 ${TEST_OUT_DIR}/maskedOuter_${MODE}_${BACKEND}_${P}_${T}.log + grep 'Test OK' ${TEST_OUT_DIR}/maskedOuter_${MODE}_${BACKEND}_${P}_${T}.log || echo "Test FAILED" + echo " " + echo ">>> [x] [ ] Testing vector times matrix using the normal (+,*)" echo " semiring over integers on a diagonal matrix" echo " " From de9a6ce0dd967a0af18ffacccb909ebcf9f9d608 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 13 Feb 2024 20:03:40 +0100 Subject: [PATCH 03/51] OpenMP support and code changes --- include/graphblas/reference/blas3.hpp | 56 +++++++++++++++++++-------- 1 file changed, 40 insertions(+), 16 deletions(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 4baa2965a..309f8ab49 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -971,6 +971,8 @@ namespace grb { const size_t m = grb::nrows( mask ); const size_t n = grb::ncols( mask ); + constexpr bool crs_only = descr & descriptors::force_row_major; + assert( phase != TRY ); if( nrows != grb::nrows( A ) || nrows != m ) { return MISMATCH; @@ -984,6 +986,9 @@ namespace grb { const auto &mask_raw = internal::getCRS( mask ); size_t nzc = 0; +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for reduction(+:nzc) +#endif for( size_t i = 0; i < nrows; ++i ) { if( internal::getCoordinates( u ).assigned( i ) ) { for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { @@ -1036,13 +1041,20 @@ namespace grb { config::NonzeroIndexType * A_col_index = internal::template getReferenceBuffer< typename config::NonzeroIndexType >( ncols + 1 ); - + internal::Coordinates< reference > coors; coors.set( arr, false, buf, ncols ); CRS_raw.col_start[ 0 ] = 0; - for( size_t j = 0; j <= ncols; ++j ) { - CCS_raw.col_start[ j ] = 0; + + if( !crs_only ) { + +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for simd +#endif + for( size_t j = 0; j <= ncols; ++j ) { + CCS_raw.col_start[ j ] = 0; + } } @@ -1053,18 +1065,27 @@ namespace grb { const size_t k_col = mask_raw.row_index[ k ]; if( internal::getCoordinates( v ).assigned( k_col ) ) { (void) ++nzc; - (void) ++CCS_raw.col_start[ k_col + 1 ]; + if( !crs_only ) { + (void) ++CCS_raw.col_start[ k_col + 1 ]; + } } } } CRS_raw.col_start[ i + 1 ] = nzc; } + if( !crs_only ) { + for( size_t j = 1; j < ncols; ++j ) { + CCS_raw.col_start[ j + 1 ] += CCS_raw.col_start[ j ]; + } - A_col_index[ 0 ] = 0; - for( size_t j = 1; j < ncols; ++j ) { - CCS_raw.col_start[ j + 1 ] += CCS_raw.col_start[ j ]; - A_col_index[ j ] = 0; + +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for simd +#endif + for( size_t j = 0; j < ncols; ++j ) { + A_col_index[ j ] = 0; + } } @@ -1074,8 +1095,8 @@ namespace grb { // computational phase nzc = 0; for( size_t i = 0; i < nrows; ++i ) { - coors.clear(); if( internal::getCoordinates( u ).assigned( i ) ) { + coors.clear(); for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { const size_t k_col = mask_raw.row_index[ k ]; if( internal::getCoordinates( v ).assigned( k_col ) ) { @@ -1095,18 +1116,21 @@ namespace grb { CRS_raw.row_index[ nzc ] = j; CRS_raw.setValue( nzc, valbuf[ j ] ); // update CCS - const size_t CCS_index = A_col_index[ j ]++ + CCS_raw.col_start[ j ]; - CCS_raw.row_index[ CCS_index ] = i; - CCS_raw.setValue( CCS_index, valbuf[ j ] ); + if( !crs_only ) { + const size_t CCS_index = A_col_index[ j ]++ + CCS_raw.col_start[ j ]; + CCS_raw.row_index[ CCS_index ] = i; + CCS_raw.setValue( CCS_index, valbuf[ j ] ); + } // update count (void) ++nzc; } CRS_raw.col_start[ i + 1 ] = nzc; } - - for( size_t j = 0; j < ncols; ++j ) { - assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] == - A_col_index[ j ] ); + if( !crs_only ) { + for( size_t j = 0; j < ncols; ++j ) { + assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] == + A_col_index[ j ] ); + } } assert( nzc == old_nzc ); From 99a1040ba3bc4c4b137b0a59a6be733208f80164 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 13 Feb 2024 20:04:16 +0100 Subject: [PATCH 04/51] test improvement --- tests/unit/CMakeLists.txt | 2 +- tests/unit/maskedOuter.cpp | 284 +++++++++++++++++++++++++++---------- 2 files changed, 213 insertions(+), 73 deletions(-) diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 7191679b2..31fe0e57b 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -247,7 +247,7 @@ add_grb_executables( outer outer.cpp ) add_grb_executables( maskedOuter maskedOuter.cpp - BACKENDS reference + BACKENDS reference reference_omp ) # must generate the golden output for other tests diff --git a/tests/unit/maskedOuter.cpp b/tests/unit/maskedOuter.cpp index e0df449d2..ba6d7f8bd 100644 --- a/tests/unit/maskedOuter.cpp +++ b/tests/unit/maskedOuter.cpp @@ -24,8 +24,6 @@ using namespace grb; // sample data -static const double vec1_vals[ 3 ] = { 1, 2, 3 }; -static const double vec2_vals[ 3 ] = { 4, 5, 6 }; static const double m1[ 3 ] = { 0.5, 3.4, 5 }; @@ -37,10 +35,7 @@ static const double m2[ 6 ] = { 0.5, -1, 23, 34, -4.4, 3.14 }; static const size_t I2[ 6 ] = { 0, 0, 1, 1, 2, 2 }; static const size_t J2[ 6 ] = { 1, 2, 0, 2, 0, 1 }; -static const double mask_test1_in[ 3 ] = { 1, 1, 1 }; static const double mask_test1_expect[ 3 ] = { 4, 10, 18 }; - -static const double mask_test2_in[ 3 ] = { 1, 1, 1 }; static const double mask_test2_expect[ 3 ] = { 11, 20, 27 }; void grbProgram( const void *, const size_t in_size, int &error ) { @@ -53,15 +48,15 @@ void grbProgram( const void *, const size_t in_size, int &error ) { } // allocate - grb::Vector< double > u( 3 ); - grb::Vector< double > v( 3 ); + grb::Vector< double > u = { 1, 2, 3 }; + grb::Vector< double > v = { 4, 5, 6 }; grb::Matrix< double > Result1( 3, 3 ); grb::Matrix< double > Result2( 3, 3 ); - grb::Matrix< double > mask1( 3, 3 ); - grb::Matrix< double > mask2( 3, 3 ); - grb::Vector< double > test1( 3 ); + grb::Matrix< double > Mask1( 3, 3 ); + grb::Matrix< double > Mask2( 3, 3 ); + grb::Vector< double > test1 = { 1, 1, 1 }; grb::Vector< double > out1( 3 ); - grb::Vector< double > test2( 3 ); + grb::Vector< double > test2 = { 1, 1, 1 }; grb::Vector< double > out2( 3 ); // semiring @@ -70,83 +65,63 @@ void grbProgram( const void *, const size_t in_size, int &error ) { grb::identities::zero, grb::identities::one > ring; - // initialise vec - const double * vec_iter = &(vec1_vals[ 0 ]); - grb::RC rc = grb::buildVector( u, vec_iter, vec_iter + 3, SEQUENTIAL ); - if( rc != SUCCESS ) { - std::cerr << "\t initial buildVector FAILED\n"; - error = 5; - } + grb::RC rc = grb::buildMatrixUnique( Mask1, &( I1[ 0 ] ), &( J1[ 0 ] ), m1, 3, SEQUENTIAL ); - if( !error ) { - vec_iter = &(vec2_vals[ 0 ]); - rc = grb::buildVector( v, vec_iter, vec_iter + 3, SEQUENTIAL ); - } if( rc != SUCCESS ) { - std::cerr << "\t initial buildVector FAILED\n"; - error = 10; + std::cerr << "\t first mask buildMatrix FAILED\n"; + error = 5; } if( !error ) { - rc = grb::buildMatrixUnique( mask1, &( I1[ 0 ] ), &( J1[ 0 ] ), m1, 3, SEQUENTIAL ); + rc = grb::buildMatrixUnique( Mask2, &( I2[ 0 ] ), &( J2[ 0 ] ), m2, 6, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t second mask buildMatrix FAILED\n"; + error = 10; + } } - if( rc != SUCCESS ) { - std::cerr << "\t first buildMatrix FAILED\n"; - error = 15; - } + if( !error ) { - rc = grb::buildMatrixUnique( mask2, &( I2[ 0 ] ), &( J2[ 0 ] ), m2, 6, SEQUENTIAL ); - } - - if( rc != SUCCESS ) { - std::cerr << "\t second buildMatrix FAILED\n"; - error = 20; + rc = grb::maskedOuter( Result1, Mask1, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::maskedOuter( Result1, Mask1, u, v, ring.getMultiplicativeOperator() ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::maskedOuter: " + << toString( rc ) << ".\n"; + error = 15; + } } - if( !error ) { - rc = grb::maskedOuter( Result1, mask1, u, v, ring.getMultiplicativeOperator(), RESIZE ); - rc = rc ? rc : grb::maskedOuter( Result1, mask1, u, v, ring.getMultiplicativeOperator() ); - } - if( rc != grb::SUCCESS ) { - std::cerr << "Unexpected return code from grb::maskedOuter: " - << toString( rc ) << ".\n"; - error = 25; - } if( !error && grb::nnz( Result1 ) != 3 ) { std::cerr << "\t Unexpected number of nonzeroes in matrix: " << grb::nnz( Result1 ) << ", expected 3.\n"; - error = 30; + error = 20; } if( !error ) { - rc = grb::maskedOuter( Result2, mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); - rc = rc ? rc : grb::maskedOuter( Result2, mask2, u, v, ring.getMultiplicativeOperator() ); - } - if( rc != grb::SUCCESS ) { - std::cerr << "Unexpected return code from grb::maskedOuter: " - << toString( rc ) << ".\n"; - error = 35; + rc = grb::maskedOuter< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::maskedOuter< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator() ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::maskedOuter: " + << toString( rc ) << ".\n"; + error = 25; + } } + if( !error && grb::nnz( Result2 ) != 6 ) { std::cerr << "\t Unexpected number of nonzeroes in matrix: " << grb::nnz( Result2 ) << ", expected 6.\n"; - error = 40; + error = 30; } if( !error ) { - const double * const test1_iter = &( mask_test1_in[ 0 ] ); - rc = grb::buildVector( test1, test1_iter, test1_iter + 3, SEQUENTIAL ); - if( rc == grb::SUCCESS ) { - rc = grb::vxm( out1, test1, Result1, ring ); - } + rc = grb::vxm( out1, test1, Result1, ring ); if( rc != grb::SUCCESS ) { - std::cerr << "Unexpected return code from premultiplying M by a vector (vxm): " + std::cerr << "Unexpected return code from premultiplying Result1 by a vector (vxm): " << toString( rc ) << ".\n"; - error = 45; + error = 35; } } @@ -154,7 +129,7 @@ void grbProgram( const void *, const size_t in_size, int &error ) { if( grb::nnz( out1 ) != 3 ) { std::cerr << "\t Unexpected number of nonzeroes (premultiply): " << grb::nnz( out1 ) << ", expected 3\n"; - error = 50; + error = 40; } for( const auto &pair : out1 ) { size_t i = pair.first; @@ -163,20 +138,16 @@ void grbProgram( const void *, const size_t in_size, int &error ) { << "unexpected value " << pair.second << " " << "at coordinate " << i << ", expected " << mask_test1_expect[ i ] << ".\n"; - error = 55; + error = 45; break; } } } if( !error ) { - const double * const test2_iter = &( mask_test2_in[ 0 ] ); - rc = grb::buildVector( test2, test2_iter, test2_iter + 3, SEQUENTIAL ); - if( rc == grb::SUCCESS ) { - rc = grb::vxm< grb::descriptors::transpose_matrix >( out2, test2, Result2, ring ); - } + rc = grb::vxm< grb::descriptors::transpose_matrix >( out2, test2, Result2, ring ); if( rc != grb::SUCCESS ) { - std::cerr << "Unexpected return code from postmultiplying M by a vector (vxm): " + std::cerr << "Unexpected return code from postmultiplying Result2 by a vector (vxm): " << toString( rc ) << ".\n"; error = 60; } @@ -191,7 +162,7 @@ void grbProgram( const void *, const size_t in_size, int &error ) { for( const auto &pair : out2 ) { size_t i = pair.first; if( pair.second != mask_test2_expect[ i ] ) { - std::cerr << "Postmultiplying M by a vector of all ones, " + std::cerr << "Postmultiplying Result2 by a vector of all ones, " << "unexpected value " << pair.second << " " << "at coordinate " << i << ", " << "expected " << mask_test2_expect[ i ] << ".\n"; @@ -202,6 +173,134 @@ void grbProgram( const void *, const size_t in_size, int &error ) { } } +void grb_program_custom_size( const size_t &n, int &error ) { + grb::Semiring< + grb::operators::add< double >, grb::operators::mul< double >, + grb::identities::zero, grb::identities::one + > ring; + + // initialize test + grb::Matrix< double > Result( n, n ); + grb::Matrix< double > Mask( n, n ); + size_t I[ 2 * n - 1 ], J[ 2 * n - 1 ]; + double mask_in [ 2 * n - 1 ]; + double u_in[ n ]; + double v_in[ n ]; + double test_in[ n ]; + double expect[ n ]; + grb::Vector< double > u( n ); + grb::Vector< double > v( n ); + grb::Vector< double > test( n ); + grb::Vector< double > out( n ); + for( size_t k = 0; k < n; ++k ) { + I[ k ] = J[ k ] = k; + mask_in[ k ] = 1; + u_in[ k ] = k + 1; + v_in[ k ] = k + 1; + test_in[ k ] = 1; + if( k < n - 1 ) { + I[ n + k ] = k; + J[ n + k ] = k + 1; + mask_in [ n + k ] = 1; + } + if( k == 0 ) { + expect [ k ] = 1; + } + else { + expect [ k ] = ( k + 1 ) * ( 2 * k + 1 ); + } + } + + const double * vec_iter = &(u_in[ 0 ]); + grb::RC rc = grb::buildVector( u, vec_iter, vec_iter + n, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildVector of u vector FAILED\n"; + error = 5; + } + + if( !error ) { + vec_iter = &(v_in[ 0 ]); + rc = grb::buildVector( v, vec_iter, vec_iter + n, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildVector of v vector FAILED\n"; + error = 10; + } + } + + if( !error ) { + vec_iter = &(test_in[ 0 ]); + rc = grb::buildVector( test, vec_iter, vec_iter + n, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildVector of test input vector FAILED\n"; + error = 15; + } + } + + if( !error ) { + rc = grb::resize( Mask, 2 * n - 1 ); + if( rc != SUCCESS ) { + std::cerr << "\t mask matrix resize FAILED\n"; + error = 20; + } + } + + if( !error ) { + rc = grb::buildMatrixUnique( Mask, I, J, mask_in, 2 * n - 1, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildMatrixUnique of mask matrix FAILED\n"; + error = 25; + } + } + + if( !error ) { + rc = grb::maskedOuter( Result, Mask, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::maskedOuter( Result, Mask, u, v, ring.getMultiplicativeOperator() ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::maskedOuter: " + << toString( rc ) << ".\n"; + error = 30; + } + } + + + if( !error && grb::nnz( Result ) != 2 * n -1 ) { + std::cerr << "\t Unexpected number of nonzeroes in matrix: " + << grb::nnz( Result ) << ", expected " + << 2 * n - 1 <<".\n"; + error = 35; + } + + if( !error ) { + rc = grb::vxm( out, test, Result, ring ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from premultiplying Result by a vector (vxm): " + << toString( rc ) << ".\n"; + error = 40; + } + } + + if( !error ) { + if( grb::nnz( out ) != n ) { + std::cerr << "\t Unexpected number of nonzeroes (premultiply): " + << grb::nnz( out ) << ", expected " + << n << ".\n"; + error = 45; + } + for( const auto &pair : out ) { + size_t i = pair.first; + if( pair.second != expect[ i ] ) { + std::cerr << "Premultiplying Result by a vector of all ones, " + << "unexpected value " << pair.second << " " + << "at coordinate " << i << ", expected " + << expect[ i ] << ".\n"; + error = 50; + break; + } + } + } + +} + int main( int argc, char ** argv ) { (void)argc; std::cout << "Functional test executable: " << argv[ 0 ] << "\n"; @@ -210,14 +309,55 @@ int main( int argc, char ** argv ) { grb::Launcher< AUTOMATIC > launcher; if( launcher.exec( &grbProgram, nullptr, 0, error ) != SUCCESS ) { - std::cerr << "Test failed to launch\n"; + std::cerr << "Test 1 failed to launch\n"; + error = 255; + } + if( error == 0 ) { + std::cout << "Test 1 OK\n" << std::endl; + } else { + std::cerr << std::flush; + std::cout << "Test 1 FAILED\n" << std::endl; + } + + bool printUsage = false; + size_t n = 100; + + // error checking + if( argc > 2 ) { + printUsage = true; + } + if( argc == 2 ) { + size_t read; + std::istringstream ss( argv[ 1 ] ); + if( ! ( ss >> read ) ) { + std::cerr << "Error parsing first argument\n"; + printUsage = true; + } else if( ! ss.eof() ) { + std::cerr << "Error parsing first argument\n"; + printUsage = true; + } else { + // all OK + n = read; + } + } + if( printUsage ) { + std::cerr << "Usage: " << argv[ 0 ] << " [n]\n"; + std::cerr << " -n (optional, default is 100): an integer, the " + "test size.\n"; + return 1; + } + + error = 0; + + if( launcher.exec( &grb_program_custom_size, n, error ) != SUCCESS ) { + std::cerr << "Launching test 2 FAILED\n"; error = 255; } if( error == 0 ) { - std::cout << "Test OK\n" << std::endl; + std::cout << "Test 2 OK\n" << std::endl; } else { std::cerr << std::flush; - std::cout << "Test FAILED\n" << std::endl; + std::cout << "Test 2 FAILED\n" << std::endl; } // done From f5562e107aaa871140e80881e5624be04ea188c2 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 13 Feb 2024 20:35:08 +0100 Subject: [PATCH 05/51] code cleanup --- include/graphblas/reference/blas3.hpp | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 309f8ab49..61d8dc42d 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -1018,11 +1018,6 @@ namespace grb { } } - grb::Monoid< - grb::operators::left_assign< OutputType >, - grb::identities::zero - > monoid; - RC ret = SUCCESS; if( phase == EXECUTE ) { ret = grb::clear( A ); @@ -1101,8 +1096,7 @@ namespace grb { const size_t k_col = mask_raw.row_index[ k ]; if( internal::getCoordinates( v ).assigned( k_col ) ) { coors.assign( k_col ); - valbuf[ k_col ] = monoid.template getIdentity< OutputType >(); - (void) grb::apply( valbuf[ k_col ], + grb::apply( valbuf[ k_col ], x[i], y[k_col], mul ); From f64dee8d2a1f1e479e3769eba5b5708bf0e2fbf2 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 13 Feb 2024 20:38:46 +0100 Subject: [PATCH 06/51] clear the matrix in special cases --- include/graphblas/reference/blas3.hpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 61d8dc42d..89f051f9f 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -982,6 +982,11 @@ namespace grb { return MISMATCH; } + if( nnz( mask ) == 0 || nnz( u ) == 0 || nnz( v ) == 0 ) { + clear( A ); + return SUCCESS; + } + const auto &mask_raw = internal::getCRS( mask ); From ef07fec6fe00b3090cab6938fdc1552d89973642 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 13 Feb 2024 20:45:59 +0100 Subject: [PATCH 07/51] special case of an empty mask --- include/graphblas/reference/blas3.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 89f051f9f..21b9125cb 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -982,11 +982,15 @@ namespace grb { return MISMATCH; } - if( nnz( mask ) == 0 || nnz( u ) == 0 || nnz( v ) == 0 ) { + if( nnz( u ) == 0 || nnz( v ) == 0 ) { clear( A ); return SUCCESS; } + if( nnz(mask) == 0 ) { + return outer( A, u, v, mul, phase ); + } + const auto &mask_raw = internal::getCRS( mask ); From b79e7ff1aa2c04ff1c418531dc7a55e02b9e7439 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Thu, 15 Feb 2024 22:27:32 +0100 Subject: [PATCH 08/51] empty mask special case --- include/graphblas/reference/blas3.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 21b9125cb..f9116368c 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -971,6 +971,10 @@ namespace grb { const size_t m = grb::nrows( mask ); const size_t n = grb::ncols( mask ); + if( m == 0 || n == 0 ) { + return outer< descr >( A, u, v, mul, phase ); + } + constexpr bool crs_only = descr & descriptors::force_row_major; assert( phase != TRY ); @@ -987,10 +991,6 @@ namespace grb { return SUCCESS; } - if( nnz(mask) == 0 ) { - return outer( A, u, v, mul, phase ); - } - const auto &mask_raw = internal::getCRS( mask ); From 7b3f795d4ee56d6dda227ecc9c16ae9d430b3cc4 Mon Sep 17 00:00:00 2001 From: Benjamin Lozes Date: Fri, 16 Feb 2024 10:38:26 +0100 Subject: [PATCH 09/51] nonblocking implementation --- include/graphblas/nonblocking/blas3.hpp | 44 +++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/include/graphblas/nonblocking/blas3.hpp b/include/graphblas/nonblocking/blas3.hpp index c52c9aaff..b8e49515b 100644 --- a/include/graphblas/nonblocking/blas3.hpp +++ b/include/graphblas/nonblocking/blas3.hpp @@ -416,6 +416,50 @@ namespace grb { ); } + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType1, typename InputType2, typename OutputType, + typename Coords, + typename RIT, typename CIT, typename NIT + > + RC maskedOuter( + Matrix< OutputType, nonblocking, RIT, CIT, NIT > &A, + const Matrix< OutputType, nonblocking, RIT, CIT, NIT > &mask, + const Vector< InputType1, nonblocking, Coords > &u, + const Vector< InputType2, nonblocking, Coords > &v, + const Operator &mul = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + if( internal::NONBLOCKING::warn_if_not_native && + config::PIPELINE::warn_if_not_native + ) { + std::cerr << "Warning: maskedOuter (nonblocking) currently delegates " + << "to a blocking implementation.\n" + << " Further similar such warnings will be suppressed.\n"; + internal::NONBLOCKING::warn_if_not_native = false; + } + + // nonblocking execution is not supported + // first, execute any computation that is not completed + internal::le.execution(); + + // second, delegate to the reference backend + return maskedOuter< descr, Operator >( + internal::getRefMatrix( A ), + internal::getRefMatrix( mask ), + internal::getRefVector( u ), + internal::getRefVector( v ), + mul, phase + ); + } + namespace internal { template< From 8dff749f2597fc52f074f6b553e7762225032e8a Mon Sep 17 00:00:00 2001 From: Benjamin Lozes Date: Fri, 16 Feb 2024 10:53:01 +0100 Subject: [PATCH 10/51] Implementation in base + documentation --- include/graphblas/base/blas3.hpp | 90 +++++++++++++++++++++++++ include/graphblas/nonblocking/blas3.hpp | 6 +- include/graphblas/reference/blas3.hpp | 15 +++-- 3 files changed, 102 insertions(+), 9 deletions(-) diff --git a/include/graphblas/base/blas3.hpp b/include/graphblas/base/blas3.hpp index c22bfba71..95b2b439b 100644 --- a/include/graphblas/base/blas3.hpp +++ b/include/graphblas/base/blas3.hpp @@ -442,6 +442,96 @@ namespace grb { return ret == SUCCESS ? UNSUPPORTED : ret; } + /** + * Masked outer product of two vectors. Assuming vectors \a u and \a v are + * oriented column-wise, the result matrix \a A will contain \f$ uv^T \f$, + * masked to non-zero values from \a mask.. This is an out-of-place + * function and will be updated soon to be in-place instead. + * + * \internal Implemented via mxm as a multiplication of a column vector with + * a row vector, while following the given mask. + * + * @tparam descr The descriptor to be used. Optional; the default is + * #grb::descriptors::no_operation. + * @tparam Operator The operator to use. + * @tparam InputType1 The value type of the left-hand vector. + * @tparam InputType2 The value type of the right-hand vector. + * @tparam MaskType The value type of the mask matrix. + * @tparam OutputType The value type of the ouput matrix + * + * @param[out] A The output matrix. + * @param[out] mask The mask matrix. + * @param[in] u The left-hand input vector. + * @param[in] v The right-hand input vector. + * @param[in] op The operator. + * @param[in] phase The #grb::Phase the call should execute. Optional; the + * default parameter is #grb::EXECUTE. + * + * @return #grb::SUCCESS On successful completion of this call. + * @return #grb::MISMATCH Whenever the dimensions of the inputs and/or outputs + * do not match. All input data containers are left + * untouched if this exit code is returned; it will be + * be as though this call was never made. + * @return #grb::FAILED If \a phase is #grb::EXECUTE, indicates that the + * capacity of \a A was insufficient. The output + * matrix \a A is cleared, and the call to this function + * has no further effects. + * @return #grb::OUTOFMEM If \a phase is #grb::RESIZE, indicates an + * out-of-memory exception. The call to this function + * shall have no other effects beyond returning this + * error code; the previous state of \a A is retained. + * @return #grb::PANIC A general unmitigable error has been encountered. If + * returned, ALP enters an undefined state and the user + * program is encouraged to exit as quickly as possible. + * + * \par Performance semantics + * + * Each backend must define performance semantics for this primitive. + * + * @see perfSemantics + */ + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType1, typename InputType2, + typename OutputType, typename MaskType, + typename Coords, + typename RIT, typename CIT, typename NIT, + Backend backend + > + RC maskedOuter( + Matrix< OutputType, backend, RIT, CIT, NIT > &A, + const Matrix< MaskType, backend, RIT, CIT, NIT > &mask, + const Vector< InputType1, backend, Coords > &u, + const Vector< InputType2, backend, Coords > &v, + const Operator &op = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + !grb::is_object< MaskType >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + (void) A; + (void) mask; + (void) u; + (void) v; + (void) op; + (void) phase; +#ifdef _DEBUG + std::cerr << "Selected backend does not implement grb::maskedOuter\n"; +#endif +#ifndef NDEBUG + const bool selected_backend_does_not_support_masked_outer = false; + assert( selected_backend_does_not_support_masked_outer ); +#endif + const RC ret = grb::clear( A ); + return ret == SUCCESS ? UNSUPPORTED : ret; + } + + /** * @} */ diff --git a/include/graphblas/nonblocking/blas3.hpp b/include/graphblas/nonblocking/blas3.hpp index b8e49515b..07e47d00f 100644 --- a/include/graphblas/nonblocking/blas3.hpp +++ b/include/graphblas/nonblocking/blas3.hpp @@ -419,13 +419,14 @@ namespace grb { template< Descriptor descr = descriptors::no_operation, class Operator, - typename InputType1, typename InputType2, typename OutputType, + typename InputType1, typename InputType2, + typename OutputType, typename MaskType, typename Coords, typename RIT, typename CIT, typename NIT > RC maskedOuter( Matrix< OutputType, nonblocking, RIT, CIT, NIT > &A, - const Matrix< OutputType, nonblocking, RIT, CIT, NIT > &mask, + const Matrix< MaskType, nonblocking, RIT, CIT, NIT > &mask, const Vector< InputType1, nonblocking, Coords > &u, const Vector< InputType2, nonblocking, Coords > &v, const Operator &mul = Operator(), @@ -434,6 +435,7 @@ namespace grb { grb::is_operator< Operator >::value && !grb::is_object< InputType1 >::value && !grb::is_object< InputType2 >::value && + !grb::is_object< MaskType >::value && !grb::is_object< OutputType >::value, void >::type * const = nullptr ) { diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index f9116368c..a514f8212 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -924,13 +924,14 @@ namespace grb { template< Descriptor descr = descriptors::no_operation, class Operator, - typename InputType1, typename InputType2, typename OutputType, + typename InputType1, typename InputType2, + typename MaskType, typename OutputType, typename Coords, typename RIT, typename CIT, typename NIT > RC maskedOuter( Matrix< OutputType, reference, RIT, CIT, NIT > &A, - const Matrix< OutputType, reference, RIT, CIT, NIT > &mask, + const Matrix< MaskType, reference, RIT, CIT, NIT > &mask, const Vector< InputType1, reference, Coords > &u, const Vector< InputType2, reference, Coords > &v, const Operator &mul = Operator(), @@ -939,6 +940,7 @@ namespace grb { grb::is_operator< Operator >::value && !grb::is_object< InputType1 >::value && !grb::is_object< InputType2 >::value && + !grb::is_object< MaskType >::value && !grb::is_object< OutputType >::value, void >::type * const = nullptr ) { @@ -958,8 +960,7 @@ namespace grb { ), "grb::maskedOuter", "called with an output matrix that does not match the output domain of " "the given multiplication operator" ); - static_assert( ( !(descr & descriptors::invert_mask) - ), + static_assert( !(descr & descriptors::invert_mask), "grb::maskedOuter: invert_mask descriptor cannot be used "); #ifdef _DEBUG std::cout << "In grb::maskedOuter (reference)\n"; @@ -972,6 +973,7 @@ namespace grb { const size_t n = grb::ncols( mask ); if( m == 0 || n == 0 ) { + // If the mask has a null size, it will be ignored return outer< descr >( A, u, v, mul, phase ); } @@ -1061,7 +1063,7 @@ namespace grb { } } - + nzc = 0; for( size_t i = 0; i < nrows; ++i ) { if( internal::getCoordinates( u ).assigned( i ) ) { @@ -1136,10 +1138,9 @@ namespace grb { } } assert( nzc == old_nzc ); - internal::setCurrentNonzeroes( A, nzc ); - + return ret; } From 0a89a08033a9b08743deeb70b4edc79b2b901ed7 Mon Sep 17 00:00:00 2001 From: Benjamin Lozes Date: Fri, 16 Feb 2024 10:58:37 +0100 Subject: [PATCH 11/51] hyperdags implementation --- include/graphblas/hyperdags/blas3.hpp | 53 +++++++++++++++++++++++ include/graphblas/hyperdags/hyperdags.hpp | 5 ++- src/graphblas/hyperdags/hyperdags.cpp | 5 ++- 3 files changed, 61 insertions(+), 2 deletions(-) diff --git a/include/graphblas/hyperdags/blas3.hpp b/include/graphblas/hyperdags/blas3.hpp index ee0c10f36..aa74dc6be 100644 --- a/include/graphblas/hyperdags/blas3.hpp +++ b/include/graphblas/hyperdags/blas3.hpp @@ -257,6 +257,59 @@ namespace grb { return ret; } + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType1, typename InputType2, + typename MaskType, typename OutputType, + typename Coords, + typename RIT, typename CIT, typename NIT + > + RC maskedOuter( + Matrix< OutputType, hyperdags, RIT, CIT, NIT > &A, + const Matrix< MaskType, hyperdags, RIT, CIT, NIT > &mask, + const Vector< InputType1, hyperdags, Coords > &u, + const Vector< InputType2, hyperdags, Coords > &v, + const Operator &mul = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + !grb::is_object< MaskType >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + const RC ret = maskedOuter< descr >( + internal::getMatrix( A ), + internal::getMatrix( mask ), + internal::getVector( u ), + internal::getVector( v ), + mul, + phase + ); + if( ret != SUCCESS ) { return ret; } + if( phase != EXECUTE ) { return ret; } + if( nrows( A ) == 0 || ncols( A ) == 0 ) { return ret; } + std::array< const void *, 0 > sourcesP{}; + std::array< uintptr_t, 4 > sourcesC{ + getID( internal::getVector(u) ), + getID( internal::getVector(v) ), + getID( internal::getMatrix(mask) ), + getID( internal::getMatrix(A) ) + }; + std::array< uintptr_t, 1 > destinations{ + getID( internal::getMatrix(A) ) + }; + internal::hyperdags::generator.addOperation( + internal::hyperdags::MASKED_OUTER, + sourcesP.begin(), sourcesP.end(), + sourcesC.begin(), sourcesC.end(), + destinations.begin(), destinations.end() + ); + return ret; + } + template< Descriptor descr = descriptors::no_operation, typename OutputType, typename InputType1, typename InputType2, diff --git a/include/graphblas/hyperdags/hyperdags.hpp b/include/graphblas/hyperdags/hyperdags.hpp index 4ef0e0059..db5d33203 100644 --- a/include/graphblas/hyperdags/hyperdags.hpp +++ b/include/graphblas/hyperdags/hyperdags.hpp @@ -386,6 +386,8 @@ namespace grb { OUTER, + MASKED_OUTER, + UNZIP_VECTOR_VECTOR_VECTOR, ZIP_MATRIX_VECTOR_VECTOR_VECTOR, @@ -493,7 +495,7 @@ namespace grb { }; /** \internal How many operation vertex types exist. */ - const constexpr size_t numOperationVertexTypes = 106; + const constexpr size_t numOperationVertexTypes = 107; /** \internal An array of all operation vertex types. */ const constexpr enum OperationVertexType @@ -553,6 +555,7 @@ namespace grb { MXM_MATRIX_MATRIX_MATRIX_SEMIRING, MXM_MATRIX_MATRIX_MATRIX_MONOID, OUTER, + MASKED_OUTER, UNZIP_VECTOR_VECTOR_VECTOR, ZIP_MATRIX_VECTOR_VECTOR_VECTOR, ZIP_MATRIX_VECTOR_VECTOR, diff --git a/src/graphblas/hyperdags/hyperdags.cpp b/src/graphblas/hyperdags/hyperdags.cpp index 6000f3af7..da1c62b19 100644 --- a/src/graphblas/hyperdags/hyperdags.cpp +++ b/src/graphblas/hyperdags/hyperdags.cpp @@ -219,7 +219,10 @@ std::string grb::internal::hyperdags::toString( return "mxm( matrix, matrix, matrix, monoid, scalar, scalar )"; case OUTER: - return "outer( matrix, vector, vector, scalar, scalar )"; + return "outer( matrix, vector, vector, operation )"; + + case MASKED_OUTER: + return "maskedOuter( matrix, matrix, vector, vector, operation )"; case MXV_VECTOR_VECTOR_MATRIX_VECTOR_VECTOR_R: return "mxv( vector, vector, matrix, vector, vector, ring )"; From 6843e53dce691bca404f425190182db237abef21 Mon Sep 17 00:00:00 2001 From: Benjamin Lozes Date: Fri, 16 Feb 2024 11:01:58 +0100 Subject: [PATCH 12/51] bsp1d+hybrid implementation --- include/graphblas/bsp1d/blas3.hpp | 44 +++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/include/graphblas/bsp1d/blas3.hpp b/include/graphblas/bsp1d/blas3.hpp index 386beb164..05f7d3c97 100644 --- a/include/graphblas/bsp1d/blas3.hpp +++ b/include/graphblas/bsp1d/blas3.hpp @@ -205,6 +205,50 @@ namespace grb { return internal::checkGlobalErrorStateOrClear( C, ret ); } + /** \internal Simply delegates to process-local backend */ + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType1, typename InputType2, + typename MaskType, typename OutputType, + typename Coords, + typename RIT, typename CIT, typename NIT + > + RC maskedOuter( + Matrix< OutputType, BSP1D, RIT, CIT, NIT > &A, + const Matrix< MaskType, BSP1D, RIT, CIT, NIT > &mask, + const Vector< InputType1, BSP1D, Coords > &u, + const Vector< InputType2, BSP1D, Coords > &v, + const Operator &mul = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + !grb::is_object< MaskType >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + assert( phase != TRY ); + RC ret = maskedOuter< descr, Operator >( + internal::getLocal( A ), + internal::getLocal( mask ), + internal::getLocal( u ), + internal::getLocal( v ), + mul, + phase + ); + if( phase == RESIZE ) { + if( collectives<>::allreduce( ret, operators::any_or< RC >() ) != SUCCESS ) { + return PANIC; + } else { + return SUCCESS; + } + } + assert( phase == EXECUTE ); + return internal::checkGlobalErrorStateOrClear( A, ret ); + } + } // namespace grb #endif From cbae6697baa5336c45087c0fe66a04bc7b209c21 Mon Sep 17 00:00:00 2001 From: Benjamin Lozes Date: Fri, 16 Feb 2024 11:02:13 +0100 Subject: [PATCH 13/51] Enable test for all backends --- tests/unit/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index a4b41fdb7..0d61e2521 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -260,7 +260,7 @@ add_grb_executables( outer outer.cpp ) add_grb_executables( maskedOuter maskedOuter.cpp - BACKENDS reference reference_omp + BACKENDS reference reference_omp bsp1d hybrid hyperdags nonblocking ) # must generate the golden output for other tests From 8b4813514f3b87c9a0fe12d0808fd7a6bef80618 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Fri, 16 Feb 2024 00:33:24 +0100 Subject: [PATCH 14/51] rename maskedOuter to outer --- include/graphblas/base/blas3.hpp | 4 ++-- include/graphblas/bsp1d/blas3.hpp | 4 ++-- include/graphblas/hyperdags/blas3.hpp | 4 ++-- include/graphblas/nonblocking/blas3.hpp | 6 +++--- include/graphblas/reference/blas3.hpp | 12 ++++++------ src/graphblas/hyperdags/hyperdags.cpp | 2 +- tests/unit/maskedOuter.cpp | 12 ++++++------ 7 files changed, 22 insertions(+), 22 deletions(-) diff --git a/include/graphblas/base/blas3.hpp b/include/graphblas/base/blas3.hpp index 95b2b439b..611c8c159 100644 --- a/include/graphblas/base/blas3.hpp +++ b/include/graphblas/base/blas3.hpp @@ -499,7 +499,7 @@ namespace grb { typename RIT, typename CIT, typename NIT, Backend backend > - RC maskedOuter( + RC outer( Matrix< OutputType, backend, RIT, CIT, NIT > &A, const Matrix< MaskType, backend, RIT, CIT, NIT > &mask, const Vector< InputType1, backend, Coords > &u, @@ -521,7 +521,7 @@ namespace grb { (void) op; (void) phase; #ifdef _DEBUG - std::cerr << "Selected backend does not implement grb::maskedOuter\n"; + std::cerr << "Selected backend does not implement grb::outer\n"; #endif #ifndef NDEBUG const bool selected_backend_does_not_support_masked_outer = false; diff --git a/include/graphblas/bsp1d/blas3.hpp b/include/graphblas/bsp1d/blas3.hpp index 05f7d3c97..6531dbeea 100644 --- a/include/graphblas/bsp1d/blas3.hpp +++ b/include/graphblas/bsp1d/blas3.hpp @@ -214,7 +214,7 @@ namespace grb { typename Coords, typename RIT, typename CIT, typename NIT > - RC maskedOuter( + RC outer( Matrix< OutputType, BSP1D, RIT, CIT, NIT > &A, const Matrix< MaskType, BSP1D, RIT, CIT, NIT > &mask, const Vector< InputType1, BSP1D, Coords > &u, @@ -230,7 +230,7 @@ namespace grb { void >::type * const = nullptr ) { assert( phase != TRY ); - RC ret = maskedOuter< descr, Operator >( + RC ret = outer< descr, Operator >( internal::getLocal( A ), internal::getLocal( mask ), internal::getLocal( u ), diff --git a/include/graphblas/hyperdags/blas3.hpp b/include/graphblas/hyperdags/blas3.hpp index aa74dc6be..16a52b9b9 100644 --- a/include/graphblas/hyperdags/blas3.hpp +++ b/include/graphblas/hyperdags/blas3.hpp @@ -265,7 +265,7 @@ namespace grb { typename Coords, typename RIT, typename CIT, typename NIT > - RC maskedOuter( + RC outer( Matrix< OutputType, hyperdags, RIT, CIT, NIT > &A, const Matrix< MaskType, hyperdags, RIT, CIT, NIT > &mask, const Vector< InputType1, hyperdags, Coords > &u, @@ -280,7 +280,7 @@ namespace grb { !grb::is_object< OutputType >::value, void >::type * const = nullptr ) { - const RC ret = maskedOuter< descr >( + const RC ret = outer< descr >( internal::getMatrix( A ), internal::getMatrix( mask ), internal::getVector( u ), diff --git a/include/graphblas/nonblocking/blas3.hpp b/include/graphblas/nonblocking/blas3.hpp index 07e47d00f..9cb925efd 100644 --- a/include/graphblas/nonblocking/blas3.hpp +++ b/include/graphblas/nonblocking/blas3.hpp @@ -424,7 +424,7 @@ namespace grb { typename Coords, typename RIT, typename CIT, typename NIT > - RC maskedOuter( + RC outer( Matrix< OutputType, nonblocking, RIT, CIT, NIT > &A, const Matrix< MaskType, nonblocking, RIT, CIT, NIT > &mask, const Vector< InputType1, nonblocking, Coords > &u, @@ -442,7 +442,7 @@ namespace grb { if( internal::NONBLOCKING::warn_if_not_native && config::PIPELINE::warn_if_not_native ) { - std::cerr << "Warning: maskedOuter (nonblocking) currently delegates " + std::cerr << "Warning: outer (nonblocking) currently delegates " << "to a blocking implementation.\n" << " Further similar such warnings will be suppressed.\n"; internal::NONBLOCKING::warn_if_not_native = false; @@ -453,7 +453,7 @@ namespace grb { internal::le.execution(); // second, delegate to the reference backend - return maskedOuter< descr, Operator >( + return outer< descr, Operator >( internal::getRefMatrix( A ), internal::getRefMatrix( mask ), internal::getRefVector( u ), diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index a514f8212..4540b7532 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -929,7 +929,7 @@ namespace grb { typename Coords, typename RIT, typename CIT, typename NIT > - RC maskedOuter( + RC outer( Matrix< OutputType, reference, RIT, CIT, NIT > &A, const Matrix< MaskType, reference, RIT, CIT, NIT > &mask, const Vector< InputType1, reference, Coords > &u, @@ -947,23 +947,23 @@ namespace grb { // static checks NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || std::is_same< typename Operator::D1, InputType1 >::value - ), "grb::maskedOuter", + ), "grb::outer", "called with a prefactor vector that does not match the first domain " "of the given multiplication operator" ); NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || std::is_same< typename Operator::D2, InputType2 >::value - ), "grb::maskedOuter", + ), "grb::outer", "called with a postfactor vector that does not match the first domain " "of the given multiplication operator" ); NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || std::is_same< typename Operator::D3, OutputType >::value - ), "grb::maskedOuter", + ), "grb::outer", "called with an output matrix that does not match the output domain of " "the given multiplication operator" ); static_assert( !(descr & descriptors::invert_mask), - "grb::maskedOuter: invert_mask descriptor cannot be used "); + "grb::outer: invert_mask descriptor cannot be used "); #ifdef _DEBUG - std::cout << "In grb::maskedOuter (reference)\n"; + std::cout << "In grb::outer (reference)\n"; #endif const size_t nrows = size( u ); diff --git a/src/graphblas/hyperdags/hyperdags.cpp b/src/graphblas/hyperdags/hyperdags.cpp index da1c62b19..011857b1f 100644 --- a/src/graphblas/hyperdags/hyperdags.cpp +++ b/src/graphblas/hyperdags/hyperdags.cpp @@ -222,7 +222,7 @@ std::string grb::internal::hyperdags::toString( return "outer( matrix, vector, vector, operation )"; case MASKED_OUTER: - return "maskedOuter( matrix, matrix, vector, vector, operation )"; + return "outer( matrix, matrix, vector, vector, operation )"; case MXV_VECTOR_VECTOR_MATRIX_VECTOR_VECTOR_R: return "mxv( vector, vector, matrix, vector, vector, ring )"; diff --git a/tests/unit/maskedOuter.cpp b/tests/unit/maskedOuter.cpp index ba6d7f8bd..1b4753666 100644 --- a/tests/unit/maskedOuter.cpp +++ b/tests/unit/maskedOuter.cpp @@ -83,8 +83,8 @@ void grbProgram( const void *, const size_t in_size, int &error ) { if( !error ) { - rc = grb::maskedOuter( Result1, Mask1, u, v, ring.getMultiplicativeOperator(), RESIZE ); - rc = rc ? rc : grb::maskedOuter( Result1, Mask1, u, v, ring.getMultiplicativeOperator() ); + rc = grb::outer( Result1, Mask1, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer( Result1, Mask1, u, v, ring.getMultiplicativeOperator() ); if( rc != grb::SUCCESS ) { std::cerr << "Unexpected return code from grb::maskedOuter: " << toString( rc ) << ".\n"; @@ -100,8 +100,8 @@ void grbProgram( const void *, const size_t in_size, int &error ) { } if( !error ) { - rc = grb::maskedOuter< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); - rc = rc ? rc : grb::maskedOuter< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator() ); + rc = grb::outer< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator() ); if( rc != grb::SUCCESS ) { std::cerr << "Unexpected return code from grb::maskedOuter: " << toString( rc ) << ".\n"; @@ -253,8 +253,8 @@ void grb_program_custom_size( const size_t &n, int &error ) { } if( !error ) { - rc = grb::maskedOuter( Result, Mask, u, v, ring.getMultiplicativeOperator(), RESIZE ); - rc = rc ? rc : grb::maskedOuter( Result, Mask, u, v, ring.getMultiplicativeOperator() ); + rc = grb::outer( Result, Mask, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer( Result, Mask, u, v, ring.getMultiplicativeOperator() ); if( rc != grb::SUCCESS ) { std::cerr << "Unexpected return code from grb::maskedOuter: " << toString( rc ) << ".\n"; From abcc35609927e1a0d463c7ae1fa347758c600833 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Fri, 26 Jan 2024 13:47:44 +0100 Subject: [PATCH 15/51] Masked outer product implementation --- include/graphblas/reference/blas3.hpp | 200 ++++++++++++++++++++++++++ 1 file changed, 200 insertions(+) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 868547183..4baa2965a 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -916,6 +916,206 @@ namespace grb { return ret; } + + /** + * A masked outer product of two vectors. Assuming vectors \a u and \a v are oriented + * column-wise, the result matrix \a A will contain \f$ uv^T \f$, masked to non-zero values from \a mask. + */ + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType1, typename InputType2, typename OutputType, + typename Coords, + typename RIT, typename CIT, typename NIT + > + RC maskedOuter( + Matrix< OutputType, reference, RIT, CIT, NIT > &A, + const Matrix< OutputType, reference, RIT, CIT, NIT > &mask, + const Vector< InputType1, reference, Coords > &u, + const Vector< InputType2, reference, Coords > &v, + const Operator &mul = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + // static checks + NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || + std::is_same< typename Operator::D1, InputType1 >::value + ), "grb::maskedOuter", + "called with a prefactor vector that does not match the first domain " + "of the given multiplication operator" ); + NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || + std::is_same< typename Operator::D2, InputType2 >::value + ), "grb::maskedOuter", + "called with a postfactor vector that does not match the first domain " + "of the given multiplication operator" ); + NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || + std::is_same< typename Operator::D3, OutputType >::value + ), "grb::maskedOuter", + "called with an output matrix that does not match the output domain of " + "the given multiplication operator" ); + static_assert( ( !(descr & descriptors::invert_mask) + ), + "grb::maskedOuter: invert_mask descriptor cannot be used "); +#ifdef _DEBUG + std::cout << "In grb::maskedOuter (reference)\n"; +#endif + + const size_t nrows = size( u ); + const size_t ncols = size( v ); + + const size_t m = grb::nrows( mask ); + const size_t n = grb::ncols( mask ); + + assert( phase != TRY ); + if( nrows != grb::nrows( A ) || nrows != m ) { + return MISMATCH; + } + + if( ncols != grb::ncols( A ) || ncols != n ) { + return MISMATCH; + } + + + const auto &mask_raw = internal::getCRS( mask ); + + size_t nzc = 0; + for( size_t i = 0; i < nrows; ++i ) { + if( internal::getCoordinates( u ).assigned( i ) ) { + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + if( internal::getCoordinates( v ).assigned( k_col ) ) { + nzc++; + } + } + } + } + + if( phase == RESIZE ) { + return resize( A, nzc ); + } + + assert( phase == EXECUTE ); + if( capacity( A ) < nzc ) { +#ifdef _DEBUG + std::cout << "\t insufficient capacity to complete " + "requested masked outer-product computation\n"; +#endif + const RC clear_rc = clear( A ); + if( clear_rc != SUCCESS ) { + return PANIC; + } else { + return FAILED; + } + } + + grb::Monoid< + grb::operators::left_assign< OutputType >, + grb::identities::zero + > monoid; + + RC ret = SUCCESS; + if( phase == EXECUTE ) { + ret = grb::clear( A ); + } + assert( nnz( A ) == 0 ); + + auto &CRS_raw = internal::getCRS( A ); + auto &CCS_raw = internal::getCCS( A ); + + const InputType1 * __restrict__ const x = internal::getRaw( u ); + const InputType2 * __restrict__ const y = internal::getRaw( v ); + char * arr = nullptr; + char * buf = nullptr; + OutputType * valbuf = nullptr; + internal::getMatrixBuffers( arr, buf, valbuf, 1, A ); + config::NonzeroIndexType * A_col_index = internal::template + getReferenceBuffer< typename config::NonzeroIndexType >( ncols + 1 ); + + + internal::Coordinates< reference > coors; + coors.set( arr, false, buf, ncols ); + + CRS_raw.col_start[ 0 ] = 0; + for( size_t j = 0; j <= ncols; ++j ) { + CCS_raw.col_start[ j ] = 0; + } + + + nzc = 0; + for( size_t i = 0; i < nrows; ++i ) { + if( internal::getCoordinates( u ).assigned( i ) ) { + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + if( internal::getCoordinates( v ).assigned( k_col ) ) { + (void) ++nzc; + (void) ++CCS_raw.col_start[ k_col + 1 ]; + } + } + } + CRS_raw.col_start[ i + 1 ] = nzc; + } + + + A_col_index[ 0 ] = 0; + for( size_t j = 1; j < ncols; ++j ) { + CCS_raw.col_start[ j + 1 ] += CCS_raw.col_start[ j ]; + A_col_index[ j ] = 0; + } + + + const size_t old_nzc = nzc; + + // use previously computed CCS offset array to update CCS during the + // computational phase + nzc = 0; + for( size_t i = 0; i < nrows; ++i ) { + coors.clear(); + if( internal::getCoordinates( u ).assigned( i ) ) { + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + if( internal::getCoordinates( v ).assigned( k_col ) ) { + coors.assign( k_col ); + valbuf[ k_col ] = monoid.template getIdentity< OutputType >(); + (void) grb::apply( valbuf[ k_col ], + x[i], + y[k_col], + mul ); + } + } + } + for( size_t k = 0; k < coors.nonzeroes(); ++k ) { + assert( nzc < old_nzc ); + const size_t j = coors.index( k ); + // update CRS + CRS_raw.row_index[ nzc ] = j; + CRS_raw.setValue( nzc, valbuf[ j ] ); + // update CCS + const size_t CCS_index = A_col_index[ j ]++ + CCS_raw.col_start[ j ]; + CCS_raw.row_index[ CCS_index ] = i; + CCS_raw.setValue( CCS_index, valbuf[ j ] ); + // update count + (void) ++nzc; + } + CRS_raw.col_start[ i + 1 ] = nzc; + } + + for( size_t j = 0; j < ncols; ++j ) { + assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] == + A_col_index[ j ] ); + } + assert( nzc == old_nzc ); + + + internal::setCurrentNonzeroes( A, nzc ); + + return ret; + } + namespace internal { /** From 89b118493abb80a8e023a1fbb723b51e2479d220 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Fri, 26 Jan 2024 14:56:01 +0100 Subject: [PATCH 16/51] Add unit test for the masked outer product --- tests/unit/CMakeLists.txt | 4 + tests/unit/maskedOuter.cpp | 226 +++++++++++++++++++++++++++++++++++++ tests/unit/unittests.sh | 6 + 3 files changed, 236 insertions(+) create mode 100644 tests/unit/maskedOuter.cpp diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index ab0dac498..caf3cf3e6 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -259,6 +259,10 @@ add_grb_executables( outer outer.cpp BACKENDS reference reference_omp hyperdags nonblocking ) +add_grb_executables( maskedOuter maskedOuter.cpp + BACKENDS reference +) + # must generate the golden output for other tests force_add_grb_executable( mxv mxv.cpp BACKEND reference diff --git a/tests/unit/maskedOuter.cpp b/tests/unit/maskedOuter.cpp new file mode 100644 index 000000000..e0df449d2 --- /dev/null +++ b/tests/unit/maskedOuter.cpp @@ -0,0 +1,226 @@ + +/* + * Copyright 2021 Huawei Technologies Co., Ltd. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +#include "graphblas.hpp" + + +using namespace grb; + +// sample data +static const double vec1_vals[ 3 ] = { 1, 2, 3 }; +static const double vec2_vals[ 3 ] = { 4, 5, 6 }; + +static const double m1[ 3 ] = { 0.5, 3.4, 5 }; + +static const size_t I1[ 3 ] = { 0, 1, 2 }; +static const size_t J1[ 3 ] = { 0, 1, 2 }; + +static const double m2[ 6 ] = { 0.5, -1, 23, 34, -4.4, 3.14 }; + +static const size_t I2[ 6 ] = { 0, 0, 1, 1, 2, 2 }; +static const size_t J2[ 6 ] = { 1, 2, 0, 2, 0, 1 }; + +static const double mask_test1_in[ 3 ] = { 1, 1, 1 }; +static const double mask_test1_expect[ 3 ] = { 4, 10, 18 }; + +static const double mask_test2_in[ 3 ] = { 1, 1, 1 }; +static const double mask_test2_expect[ 3 ] = { 11, 20, 27 }; + +void grbProgram( const void *, const size_t in_size, int &error ) { + error = 0; + + if( in_size != 0 ) { + (void)fprintf( stderr, "Unit tests called with unexpected input\n" ); + error = 1; + return; + } + + // allocate + grb::Vector< double > u( 3 ); + grb::Vector< double > v( 3 ); + grb::Matrix< double > Result1( 3, 3 ); + grb::Matrix< double > Result2( 3, 3 ); + grb::Matrix< double > mask1( 3, 3 ); + grb::Matrix< double > mask2( 3, 3 ); + grb::Vector< double > test1( 3 ); + grb::Vector< double > out1( 3 ); + grb::Vector< double > test2( 3 ); + grb::Vector< double > out2( 3 ); + + // semiring + grb::Semiring< + grb::operators::add< double >, grb::operators::mul< double >, + grb::identities::zero, grb::identities::one + > ring; + + // initialise vec + const double * vec_iter = &(vec1_vals[ 0 ]); + grb::RC rc = grb::buildVector( u, vec_iter, vec_iter + 3, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t initial buildVector FAILED\n"; + error = 5; + } + + if( !error ) { + vec_iter = &(vec2_vals[ 0 ]); + rc = grb::buildVector( v, vec_iter, vec_iter + 3, SEQUENTIAL ); + } + if( rc != SUCCESS ) { + std::cerr << "\t initial buildVector FAILED\n"; + error = 10; + } + + if( !error ) { + rc = grb::buildMatrixUnique( mask1, &( I1[ 0 ] ), &( J1[ 0 ] ), m1, 3, SEQUENTIAL ); + } + + if( rc != SUCCESS ) { + std::cerr << "\t first buildMatrix FAILED\n"; + error = 15; + } + + if( !error ) { + rc = grb::buildMatrixUnique( mask2, &( I2[ 0 ] ), &( J2[ 0 ] ), m2, 6, SEQUENTIAL ); + } + + if( rc != SUCCESS ) { + std::cerr << "\t second buildMatrix FAILED\n"; + error = 20; + } + + if( !error ) { + rc = grb::maskedOuter( Result1, mask1, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::maskedOuter( Result1, mask1, u, v, ring.getMultiplicativeOperator() ); + } + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::maskedOuter: " + << toString( rc ) << ".\n"; + error = 25; + } + + if( !error && grb::nnz( Result1 ) != 3 ) { + std::cerr << "\t Unexpected number of nonzeroes in matrix: " + << grb::nnz( Result1 ) << ", expected 3.\n"; + error = 30; + } + + if( !error ) { + rc = grb::maskedOuter( Result2, mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::maskedOuter( Result2, mask2, u, v, ring.getMultiplicativeOperator() ); + } + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::maskedOuter: " + << toString( rc ) << ".\n"; + error = 35; + } + + if( !error && grb::nnz( Result2 ) != 6 ) { + std::cerr << "\t Unexpected number of nonzeroes in matrix: " + << grb::nnz( Result2 ) << ", expected 6.\n"; + error = 40; + } + + if( !error ) { + const double * const test1_iter = &( mask_test1_in[ 0 ] ); + rc = grb::buildVector( test1, test1_iter, test1_iter + 3, SEQUENTIAL ); + if( rc == grb::SUCCESS ) { + rc = grb::vxm( out1, test1, Result1, ring ); + } + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from premultiplying M by a vector (vxm): " + << toString( rc ) << ".\n"; + error = 45; + } + } + + if( !error ) { + if( grb::nnz( out1 ) != 3 ) { + std::cerr << "\t Unexpected number of nonzeroes (premultiply): " + << grb::nnz( out1 ) << ", expected 3\n"; + error = 50; + } + for( const auto &pair : out1 ) { + size_t i = pair.first; + if( pair.second != mask_test1_expect[ i ] ) { + std::cerr << "Premultiplying Result1 by a vector of all ones, " + << "unexpected value " << pair.second << " " + << "at coordinate " << i << ", expected " + << mask_test1_expect[ i ] << ".\n"; + error = 55; + break; + } + } + } + + if( !error ) { + const double * const test2_iter = &( mask_test2_in[ 0 ] ); + rc = grb::buildVector( test2, test2_iter, test2_iter + 3, SEQUENTIAL ); + if( rc == grb::SUCCESS ) { + rc = grb::vxm< grb::descriptors::transpose_matrix >( out2, test2, Result2, ring ); + } + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from postmultiplying M by a vector (vxm): " + << toString( rc ) << ".\n"; + error = 60; + } + } + + if( !error ) { + if( grb::nnz( out2 ) != 3 ) { + std::cerr << "\t Unexpected number of nonzeroes (postmultiply): " + << grb::nnz( out1 ) << ", expected 3\n"; + error = 65; + } + for( const auto &pair : out2 ) { + size_t i = pair.first; + if( pair.second != mask_test2_expect[ i ] ) { + std::cerr << "Postmultiplying M by a vector of all ones, " + << "unexpected value " << pair.second << " " + << "at coordinate " << i << ", " + << "expected " << mask_test2_expect[ i ] << ".\n"; + error = 70; + break; + } + } + } +} + +int main( int argc, char ** argv ) { + (void)argc; + std::cout << "Functional test executable: " << argv[ 0 ] << "\n"; + + int error; + grb::Launcher< AUTOMATIC > launcher; + + if( launcher.exec( &grbProgram, nullptr, 0, error ) != SUCCESS ) { + std::cerr << "Test failed to launch\n"; + error = 255; + } + if( error == 0 ) { + std::cout << "Test OK\n" << std::endl; + } else { + std::cerr << std::flush; + std::cout << "Test FAILED\n" << std::endl; + } + + // done + return error; +} + diff --git a/tests/unit/unittests.sh b/tests/unit/unittests.sh index 5b7dfc16a..f3db533a0 100755 --- a/tests/unit/unittests.sh +++ b/tests/unit/unittests.sh @@ -643,6 +643,12 @@ for MODE in ${MODES}; do grep 'Test OK' ${TEST_OUT_DIR}/outer_${MODE}_${BACKEND}_${P}_${T}.log || echo "Test FAILED" echo " " + echo ">>> [x] [ ] Testing grb::maskedOuter on a small matrix" + $runner ${TEST_BIN_DIR}/maskedOuter_${MODE}_${BACKEND} &> ${TEST_OUT_DIR}/maskedOuter_${MODE}_${BACKEND}_${P}_${T}.log + head -1 ${TEST_OUT_DIR}/maskedOuter_${MODE}_${BACKEND}_${P}_${T}.log + grep 'Test OK' ${TEST_OUT_DIR}/maskedOuter_${MODE}_${BACKEND}_${P}_${T}.log || echo "Test FAILED" + echo " " + echo ">>> [x] [ ] Testing vector times matrix using the normal (+,*)" echo " semiring over integers on a diagonal matrix" echo " " From 4dc1e37a68a349c0e6862c9b2eeb62bd48cb906b Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 13 Feb 2024 20:03:40 +0100 Subject: [PATCH 17/51] OpenMP support and code changes --- include/graphblas/reference/blas3.hpp | 56 +++++++++++++++++++-------- 1 file changed, 40 insertions(+), 16 deletions(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 4baa2965a..309f8ab49 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -971,6 +971,8 @@ namespace grb { const size_t m = grb::nrows( mask ); const size_t n = grb::ncols( mask ); + constexpr bool crs_only = descr & descriptors::force_row_major; + assert( phase != TRY ); if( nrows != grb::nrows( A ) || nrows != m ) { return MISMATCH; @@ -984,6 +986,9 @@ namespace grb { const auto &mask_raw = internal::getCRS( mask ); size_t nzc = 0; +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for reduction(+:nzc) +#endif for( size_t i = 0; i < nrows; ++i ) { if( internal::getCoordinates( u ).assigned( i ) ) { for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { @@ -1036,13 +1041,20 @@ namespace grb { config::NonzeroIndexType * A_col_index = internal::template getReferenceBuffer< typename config::NonzeroIndexType >( ncols + 1 ); - + internal::Coordinates< reference > coors; coors.set( arr, false, buf, ncols ); CRS_raw.col_start[ 0 ] = 0; - for( size_t j = 0; j <= ncols; ++j ) { - CCS_raw.col_start[ j ] = 0; + + if( !crs_only ) { + +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for simd +#endif + for( size_t j = 0; j <= ncols; ++j ) { + CCS_raw.col_start[ j ] = 0; + } } @@ -1053,18 +1065,27 @@ namespace grb { const size_t k_col = mask_raw.row_index[ k ]; if( internal::getCoordinates( v ).assigned( k_col ) ) { (void) ++nzc; - (void) ++CCS_raw.col_start[ k_col + 1 ]; + if( !crs_only ) { + (void) ++CCS_raw.col_start[ k_col + 1 ]; + } } } } CRS_raw.col_start[ i + 1 ] = nzc; } + if( !crs_only ) { + for( size_t j = 1; j < ncols; ++j ) { + CCS_raw.col_start[ j + 1 ] += CCS_raw.col_start[ j ]; + } - A_col_index[ 0 ] = 0; - for( size_t j = 1; j < ncols; ++j ) { - CCS_raw.col_start[ j + 1 ] += CCS_raw.col_start[ j ]; - A_col_index[ j ] = 0; + +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for simd +#endif + for( size_t j = 0; j < ncols; ++j ) { + A_col_index[ j ] = 0; + } } @@ -1074,8 +1095,8 @@ namespace grb { // computational phase nzc = 0; for( size_t i = 0; i < nrows; ++i ) { - coors.clear(); if( internal::getCoordinates( u ).assigned( i ) ) { + coors.clear(); for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { const size_t k_col = mask_raw.row_index[ k ]; if( internal::getCoordinates( v ).assigned( k_col ) ) { @@ -1095,18 +1116,21 @@ namespace grb { CRS_raw.row_index[ nzc ] = j; CRS_raw.setValue( nzc, valbuf[ j ] ); // update CCS - const size_t CCS_index = A_col_index[ j ]++ + CCS_raw.col_start[ j ]; - CCS_raw.row_index[ CCS_index ] = i; - CCS_raw.setValue( CCS_index, valbuf[ j ] ); + if( !crs_only ) { + const size_t CCS_index = A_col_index[ j ]++ + CCS_raw.col_start[ j ]; + CCS_raw.row_index[ CCS_index ] = i; + CCS_raw.setValue( CCS_index, valbuf[ j ] ); + } // update count (void) ++nzc; } CRS_raw.col_start[ i + 1 ] = nzc; } - - for( size_t j = 0; j < ncols; ++j ) { - assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] == - A_col_index[ j ] ); + if( !crs_only ) { + for( size_t j = 0; j < ncols; ++j ) { + assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] == + A_col_index[ j ] ); + } } assert( nzc == old_nzc ); From dfafe0caddc3a81ac61f787b5d29c6a4b17b3a4d Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 13 Feb 2024 20:04:16 +0100 Subject: [PATCH 18/51] test improvement --- tests/unit/CMakeLists.txt | 2 +- tests/unit/maskedOuter.cpp | 284 +++++++++++++++++++++++++++---------- 2 files changed, 213 insertions(+), 73 deletions(-) diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index caf3cf3e6..a4b41fdb7 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -260,7 +260,7 @@ add_grb_executables( outer outer.cpp ) add_grb_executables( maskedOuter maskedOuter.cpp - BACKENDS reference + BACKENDS reference reference_omp ) # must generate the golden output for other tests diff --git a/tests/unit/maskedOuter.cpp b/tests/unit/maskedOuter.cpp index e0df449d2..ba6d7f8bd 100644 --- a/tests/unit/maskedOuter.cpp +++ b/tests/unit/maskedOuter.cpp @@ -24,8 +24,6 @@ using namespace grb; // sample data -static const double vec1_vals[ 3 ] = { 1, 2, 3 }; -static const double vec2_vals[ 3 ] = { 4, 5, 6 }; static const double m1[ 3 ] = { 0.5, 3.4, 5 }; @@ -37,10 +35,7 @@ static const double m2[ 6 ] = { 0.5, -1, 23, 34, -4.4, 3.14 }; static const size_t I2[ 6 ] = { 0, 0, 1, 1, 2, 2 }; static const size_t J2[ 6 ] = { 1, 2, 0, 2, 0, 1 }; -static const double mask_test1_in[ 3 ] = { 1, 1, 1 }; static const double mask_test1_expect[ 3 ] = { 4, 10, 18 }; - -static const double mask_test2_in[ 3 ] = { 1, 1, 1 }; static const double mask_test2_expect[ 3 ] = { 11, 20, 27 }; void grbProgram( const void *, const size_t in_size, int &error ) { @@ -53,15 +48,15 @@ void grbProgram( const void *, const size_t in_size, int &error ) { } // allocate - grb::Vector< double > u( 3 ); - grb::Vector< double > v( 3 ); + grb::Vector< double > u = { 1, 2, 3 }; + grb::Vector< double > v = { 4, 5, 6 }; grb::Matrix< double > Result1( 3, 3 ); grb::Matrix< double > Result2( 3, 3 ); - grb::Matrix< double > mask1( 3, 3 ); - grb::Matrix< double > mask2( 3, 3 ); - grb::Vector< double > test1( 3 ); + grb::Matrix< double > Mask1( 3, 3 ); + grb::Matrix< double > Mask2( 3, 3 ); + grb::Vector< double > test1 = { 1, 1, 1 }; grb::Vector< double > out1( 3 ); - grb::Vector< double > test2( 3 ); + grb::Vector< double > test2 = { 1, 1, 1 }; grb::Vector< double > out2( 3 ); // semiring @@ -70,83 +65,63 @@ void grbProgram( const void *, const size_t in_size, int &error ) { grb::identities::zero, grb::identities::one > ring; - // initialise vec - const double * vec_iter = &(vec1_vals[ 0 ]); - grb::RC rc = grb::buildVector( u, vec_iter, vec_iter + 3, SEQUENTIAL ); - if( rc != SUCCESS ) { - std::cerr << "\t initial buildVector FAILED\n"; - error = 5; - } + grb::RC rc = grb::buildMatrixUnique( Mask1, &( I1[ 0 ] ), &( J1[ 0 ] ), m1, 3, SEQUENTIAL ); - if( !error ) { - vec_iter = &(vec2_vals[ 0 ]); - rc = grb::buildVector( v, vec_iter, vec_iter + 3, SEQUENTIAL ); - } if( rc != SUCCESS ) { - std::cerr << "\t initial buildVector FAILED\n"; - error = 10; + std::cerr << "\t first mask buildMatrix FAILED\n"; + error = 5; } if( !error ) { - rc = grb::buildMatrixUnique( mask1, &( I1[ 0 ] ), &( J1[ 0 ] ), m1, 3, SEQUENTIAL ); + rc = grb::buildMatrixUnique( Mask2, &( I2[ 0 ] ), &( J2[ 0 ] ), m2, 6, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t second mask buildMatrix FAILED\n"; + error = 10; + } } - if( rc != SUCCESS ) { - std::cerr << "\t first buildMatrix FAILED\n"; - error = 15; - } + if( !error ) { - rc = grb::buildMatrixUnique( mask2, &( I2[ 0 ] ), &( J2[ 0 ] ), m2, 6, SEQUENTIAL ); - } - - if( rc != SUCCESS ) { - std::cerr << "\t second buildMatrix FAILED\n"; - error = 20; + rc = grb::maskedOuter( Result1, Mask1, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::maskedOuter( Result1, Mask1, u, v, ring.getMultiplicativeOperator() ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::maskedOuter: " + << toString( rc ) << ".\n"; + error = 15; + } } - if( !error ) { - rc = grb::maskedOuter( Result1, mask1, u, v, ring.getMultiplicativeOperator(), RESIZE ); - rc = rc ? rc : grb::maskedOuter( Result1, mask1, u, v, ring.getMultiplicativeOperator() ); - } - if( rc != grb::SUCCESS ) { - std::cerr << "Unexpected return code from grb::maskedOuter: " - << toString( rc ) << ".\n"; - error = 25; - } if( !error && grb::nnz( Result1 ) != 3 ) { std::cerr << "\t Unexpected number of nonzeroes in matrix: " << grb::nnz( Result1 ) << ", expected 3.\n"; - error = 30; + error = 20; } if( !error ) { - rc = grb::maskedOuter( Result2, mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); - rc = rc ? rc : grb::maskedOuter( Result2, mask2, u, v, ring.getMultiplicativeOperator() ); - } - if( rc != grb::SUCCESS ) { - std::cerr << "Unexpected return code from grb::maskedOuter: " - << toString( rc ) << ".\n"; - error = 35; + rc = grb::maskedOuter< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::maskedOuter< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator() ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::maskedOuter: " + << toString( rc ) << ".\n"; + error = 25; + } } + if( !error && grb::nnz( Result2 ) != 6 ) { std::cerr << "\t Unexpected number of nonzeroes in matrix: " << grb::nnz( Result2 ) << ", expected 6.\n"; - error = 40; + error = 30; } if( !error ) { - const double * const test1_iter = &( mask_test1_in[ 0 ] ); - rc = grb::buildVector( test1, test1_iter, test1_iter + 3, SEQUENTIAL ); - if( rc == grb::SUCCESS ) { - rc = grb::vxm( out1, test1, Result1, ring ); - } + rc = grb::vxm( out1, test1, Result1, ring ); if( rc != grb::SUCCESS ) { - std::cerr << "Unexpected return code from premultiplying M by a vector (vxm): " + std::cerr << "Unexpected return code from premultiplying Result1 by a vector (vxm): " << toString( rc ) << ".\n"; - error = 45; + error = 35; } } @@ -154,7 +129,7 @@ void grbProgram( const void *, const size_t in_size, int &error ) { if( grb::nnz( out1 ) != 3 ) { std::cerr << "\t Unexpected number of nonzeroes (premultiply): " << grb::nnz( out1 ) << ", expected 3\n"; - error = 50; + error = 40; } for( const auto &pair : out1 ) { size_t i = pair.first; @@ -163,20 +138,16 @@ void grbProgram( const void *, const size_t in_size, int &error ) { << "unexpected value " << pair.second << " " << "at coordinate " << i << ", expected " << mask_test1_expect[ i ] << ".\n"; - error = 55; + error = 45; break; } } } if( !error ) { - const double * const test2_iter = &( mask_test2_in[ 0 ] ); - rc = grb::buildVector( test2, test2_iter, test2_iter + 3, SEQUENTIAL ); - if( rc == grb::SUCCESS ) { - rc = grb::vxm< grb::descriptors::transpose_matrix >( out2, test2, Result2, ring ); - } + rc = grb::vxm< grb::descriptors::transpose_matrix >( out2, test2, Result2, ring ); if( rc != grb::SUCCESS ) { - std::cerr << "Unexpected return code from postmultiplying M by a vector (vxm): " + std::cerr << "Unexpected return code from postmultiplying Result2 by a vector (vxm): " << toString( rc ) << ".\n"; error = 60; } @@ -191,7 +162,7 @@ void grbProgram( const void *, const size_t in_size, int &error ) { for( const auto &pair : out2 ) { size_t i = pair.first; if( pair.second != mask_test2_expect[ i ] ) { - std::cerr << "Postmultiplying M by a vector of all ones, " + std::cerr << "Postmultiplying Result2 by a vector of all ones, " << "unexpected value " << pair.second << " " << "at coordinate " << i << ", " << "expected " << mask_test2_expect[ i ] << ".\n"; @@ -202,6 +173,134 @@ void grbProgram( const void *, const size_t in_size, int &error ) { } } +void grb_program_custom_size( const size_t &n, int &error ) { + grb::Semiring< + grb::operators::add< double >, grb::operators::mul< double >, + grb::identities::zero, grb::identities::one + > ring; + + // initialize test + grb::Matrix< double > Result( n, n ); + grb::Matrix< double > Mask( n, n ); + size_t I[ 2 * n - 1 ], J[ 2 * n - 1 ]; + double mask_in [ 2 * n - 1 ]; + double u_in[ n ]; + double v_in[ n ]; + double test_in[ n ]; + double expect[ n ]; + grb::Vector< double > u( n ); + grb::Vector< double > v( n ); + grb::Vector< double > test( n ); + grb::Vector< double > out( n ); + for( size_t k = 0; k < n; ++k ) { + I[ k ] = J[ k ] = k; + mask_in[ k ] = 1; + u_in[ k ] = k + 1; + v_in[ k ] = k + 1; + test_in[ k ] = 1; + if( k < n - 1 ) { + I[ n + k ] = k; + J[ n + k ] = k + 1; + mask_in [ n + k ] = 1; + } + if( k == 0 ) { + expect [ k ] = 1; + } + else { + expect [ k ] = ( k + 1 ) * ( 2 * k + 1 ); + } + } + + const double * vec_iter = &(u_in[ 0 ]); + grb::RC rc = grb::buildVector( u, vec_iter, vec_iter + n, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildVector of u vector FAILED\n"; + error = 5; + } + + if( !error ) { + vec_iter = &(v_in[ 0 ]); + rc = grb::buildVector( v, vec_iter, vec_iter + n, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildVector of v vector FAILED\n"; + error = 10; + } + } + + if( !error ) { + vec_iter = &(test_in[ 0 ]); + rc = grb::buildVector( test, vec_iter, vec_iter + n, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildVector of test input vector FAILED\n"; + error = 15; + } + } + + if( !error ) { + rc = grb::resize( Mask, 2 * n - 1 ); + if( rc != SUCCESS ) { + std::cerr << "\t mask matrix resize FAILED\n"; + error = 20; + } + } + + if( !error ) { + rc = grb::buildMatrixUnique( Mask, I, J, mask_in, 2 * n - 1, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildMatrixUnique of mask matrix FAILED\n"; + error = 25; + } + } + + if( !error ) { + rc = grb::maskedOuter( Result, Mask, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::maskedOuter( Result, Mask, u, v, ring.getMultiplicativeOperator() ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::maskedOuter: " + << toString( rc ) << ".\n"; + error = 30; + } + } + + + if( !error && grb::nnz( Result ) != 2 * n -1 ) { + std::cerr << "\t Unexpected number of nonzeroes in matrix: " + << grb::nnz( Result ) << ", expected " + << 2 * n - 1 <<".\n"; + error = 35; + } + + if( !error ) { + rc = grb::vxm( out, test, Result, ring ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from premultiplying Result by a vector (vxm): " + << toString( rc ) << ".\n"; + error = 40; + } + } + + if( !error ) { + if( grb::nnz( out ) != n ) { + std::cerr << "\t Unexpected number of nonzeroes (premultiply): " + << grb::nnz( out ) << ", expected " + << n << ".\n"; + error = 45; + } + for( const auto &pair : out ) { + size_t i = pair.first; + if( pair.second != expect[ i ] ) { + std::cerr << "Premultiplying Result by a vector of all ones, " + << "unexpected value " << pair.second << " " + << "at coordinate " << i << ", expected " + << expect[ i ] << ".\n"; + error = 50; + break; + } + } + } + +} + int main( int argc, char ** argv ) { (void)argc; std::cout << "Functional test executable: " << argv[ 0 ] << "\n"; @@ -210,14 +309,55 @@ int main( int argc, char ** argv ) { grb::Launcher< AUTOMATIC > launcher; if( launcher.exec( &grbProgram, nullptr, 0, error ) != SUCCESS ) { - std::cerr << "Test failed to launch\n"; + std::cerr << "Test 1 failed to launch\n"; + error = 255; + } + if( error == 0 ) { + std::cout << "Test 1 OK\n" << std::endl; + } else { + std::cerr << std::flush; + std::cout << "Test 1 FAILED\n" << std::endl; + } + + bool printUsage = false; + size_t n = 100; + + // error checking + if( argc > 2 ) { + printUsage = true; + } + if( argc == 2 ) { + size_t read; + std::istringstream ss( argv[ 1 ] ); + if( ! ( ss >> read ) ) { + std::cerr << "Error parsing first argument\n"; + printUsage = true; + } else if( ! ss.eof() ) { + std::cerr << "Error parsing first argument\n"; + printUsage = true; + } else { + // all OK + n = read; + } + } + if( printUsage ) { + std::cerr << "Usage: " << argv[ 0 ] << " [n]\n"; + std::cerr << " -n (optional, default is 100): an integer, the " + "test size.\n"; + return 1; + } + + error = 0; + + if( launcher.exec( &grb_program_custom_size, n, error ) != SUCCESS ) { + std::cerr << "Launching test 2 FAILED\n"; error = 255; } if( error == 0 ) { - std::cout << "Test OK\n" << std::endl; + std::cout << "Test 2 OK\n" << std::endl; } else { std::cerr << std::flush; - std::cout << "Test FAILED\n" << std::endl; + std::cout << "Test 2 FAILED\n" << std::endl; } // done From a73fcb19f29d83a775975c8e3a0f3a584a505832 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 13 Feb 2024 20:35:08 +0100 Subject: [PATCH 19/51] code cleanup --- include/graphblas/reference/blas3.hpp | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 309f8ab49..61d8dc42d 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -1018,11 +1018,6 @@ namespace grb { } } - grb::Monoid< - grb::operators::left_assign< OutputType >, - grb::identities::zero - > monoid; - RC ret = SUCCESS; if( phase == EXECUTE ) { ret = grb::clear( A ); @@ -1101,8 +1096,7 @@ namespace grb { const size_t k_col = mask_raw.row_index[ k ]; if( internal::getCoordinates( v ).assigned( k_col ) ) { coors.assign( k_col ); - valbuf[ k_col ] = monoid.template getIdentity< OutputType >(); - (void) grb::apply( valbuf[ k_col ], + grb::apply( valbuf[ k_col ], x[i], y[k_col], mul ); From 6b91b64bb95bbe09374e9e20db091c66b46560f4 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 13 Feb 2024 20:38:46 +0100 Subject: [PATCH 20/51] clear the matrix in special cases --- include/graphblas/reference/blas3.hpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 61d8dc42d..89f051f9f 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -982,6 +982,11 @@ namespace grb { return MISMATCH; } + if( nnz( mask ) == 0 || nnz( u ) == 0 || nnz( v ) == 0 ) { + clear( A ); + return SUCCESS; + } + const auto &mask_raw = internal::getCRS( mask ); From e91f9edaa35f43583817ab9bf0cc220676995b53 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 13 Feb 2024 20:45:59 +0100 Subject: [PATCH 21/51] special case of an empty mask --- include/graphblas/reference/blas3.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 89f051f9f..21b9125cb 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -982,11 +982,15 @@ namespace grb { return MISMATCH; } - if( nnz( mask ) == 0 || nnz( u ) == 0 || nnz( v ) == 0 ) { + if( nnz( u ) == 0 || nnz( v ) == 0 ) { clear( A ); return SUCCESS; } + if( nnz(mask) == 0 ) { + return outer( A, u, v, mul, phase ); + } + const auto &mask_raw = internal::getCRS( mask ); From 76bd1028b2282dcf91d670ad5a9e86b91e98330a Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Thu, 15 Feb 2024 22:27:32 +0100 Subject: [PATCH 22/51] empty mask special case --- include/graphblas/reference/blas3.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 21b9125cb..f9116368c 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -971,6 +971,10 @@ namespace grb { const size_t m = grb::nrows( mask ); const size_t n = grb::ncols( mask ); + if( m == 0 || n == 0 ) { + return outer< descr >( A, u, v, mul, phase ); + } + constexpr bool crs_only = descr & descriptors::force_row_major; assert( phase != TRY ); @@ -987,10 +991,6 @@ namespace grb { return SUCCESS; } - if( nnz(mask) == 0 ) { - return outer( A, u, v, mul, phase ); - } - const auto &mask_raw = internal::getCRS( mask ); From 279894c9a7576397b926e1333297f33c7415e10d Mon Sep 17 00:00:00 2001 From: Benjamin Lozes Date: Fri, 16 Feb 2024 10:38:26 +0100 Subject: [PATCH 23/51] nonblocking implementation --- include/graphblas/nonblocking/blas3.hpp | 44 +++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/include/graphblas/nonblocking/blas3.hpp b/include/graphblas/nonblocking/blas3.hpp index c52c9aaff..b8e49515b 100644 --- a/include/graphblas/nonblocking/blas3.hpp +++ b/include/graphblas/nonblocking/blas3.hpp @@ -416,6 +416,50 @@ namespace grb { ); } + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType1, typename InputType2, typename OutputType, + typename Coords, + typename RIT, typename CIT, typename NIT + > + RC maskedOuter( + Matrix< OutputType, nonblocking, RIT, CIT, NIT > &A, + const Matrix< OutputType, nonblocking, RIT, CIT, NIT > &mask, + const Vector< InputType1, nonblocking, Coords > &u, + const Vector< InputType2, nonblocking, Coords > &v, + const Operator &mul = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + if( internal::NONBLOCKING::warn_if_not_native && + config::PIPELINE::warn_if_not_native + ) { + std::cerr << "Warning: maskedOuter (nonblocking) currently delegates " + << "to a blocking implementation.\n" + << " Further similar such warnings will be suppressed.\n"; + internal::NONBLOCKING::warn_if_not_native = false; + } + + // nonblocking execution is not supported + // first, execute any computation that is not completed + internal::le.execution(); + + // second, delegate to the reference backend + return maskedOuter< descr, Operator >( + internal::getRefMatrix( A ), + internal::getRefMatrix( mask ), + internal::getRefVector( u ), + internal::getRefVector( v ), + mul, phase + ); + } + namespace internal { template< From eeed9107474f4fb61f589d7c430ee623df700ca7 Mon Sep 17 00:00:00 2001 From: Benjamin Lozes Date: Fri, 16 Feb 2024 10:53:01 +0100 Subject: [PATCH 24/51] Implementation in base + documentation --- include/graphblas/base/blas3.hpp | 90 +++++++++++++++++++++++++ include/graphblas/nonblocking/blas3.hpp | 6 +- include/graphblas/reference/blas3.hpp | 15 +++-- 3 files changed, 102 insertions(+), 9 deletions(-) diff --git a/include/graphblas/base/blas3.hpp b/include/graphblas/base/blas3.hpp index c22bfba71..95b2b439b 100644 --- a/include/graphblas/base/blas3.hpp +++ b/include/graphblas/base/blas3.hpp @@ -442,6 +442,96 @@ namespace grb { return ret == SUCCESS ? UNSUPPORTED : ret; } + /** + * Masked outer product of two vectors. Assuming vectors \a u and \a v are + * oriented column-wise, the result matrix \a A will contain \f$ uv^T \f$, + * masked to non-zero values from \a mask.. This is an out-of-place + * function and will be updated soon to be in-place instead. + * + * \internal Implemented via mxm as a multiplication of a column vector with + * a row vector, while following the given mask. + * + * @tparam descr The descriptor to be used. Optional; the default is + * #grb::descriptors::no_operation. + * @tparam Operator The operator to use. + * @tparam InputType1 The value type of the left-hand vector. + * @tparam InputType2 The value type of the right-hand vector. + * @tparam MaskType The value type of the mask matrix. + * @tparam OutputType The value type of the ouput matrix + * + * @param[out] A The output matrix. + * @param[out] mask The mask matrix. + * @param[in] u The left-hand input vector. + * @param[in] v The right-hand input vector. + * @param[in] op The operator. + * @param[in] phase The #grb::Phase the call should execute. Optional; the + * default parameter is #grb::EXECUTE. + * + * @return #grb::SUCCESS On successful completion of this call. + * @return #grb::MISMATCH Whenever the dimensions of the inputs and/or outputs + * do not match. All input data containers are left + * untouched if this exit code is returned; it will be + * be as though this call was never made. + * @return #grb::FAILED If \a phase is #grb::EXECUTE, indicates that the + * capacity of \a A was insufficient. The output + * matrix \a A is cleared, and the call to this function + * has no further effects. + * @return #grb::OUTOFMEM If \a phase is #grb::RESIZE, indicates an + * out-of-memory exception. The call to this function + * shall have no other effects beyond returning this + * error code; the previous state of \a A is retained. + * @return #grb::PANIC A general unmitigable error has been encountered. If + * returned, ALP enters an undefined state and the user + * program is encouraged to exit as quickly as possible. + * + * \par Performance semantics + * + * Each backend must define performance semantics for this primitive. + * + * @see perfSemantics + */ + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType1, typename InputType2, + typename OutputType, typename MaskType, + typename Coords, + typename RIT, typename CIT, typename NIT, + Backend backend + > + RC maskedOuter( + Matrix< OutputType, backend, RIT, CIT, NIT > &A, + const Matrix< MaskType, backend, RIT, CIT, NIT > &mask, + const Vector< InputType1, backend, Coords > &u, + const Vector< InputType2, backend, Coords > &v, + const Operator &op = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + !grb::is_object< MaskType >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + (void) A; + (void) mask; + (void) u; + (void) v; + (void) op; + (void) phase; +#ifdef _DEBUG + std::cerr << "Selected backend does not implement grb::maskedOuter\n"; +#endif +#ifndef NDEBUG + const bool selected_backend_does_not_support_masked_outer = false; + assert( selected_backend_does_not_support_masked_outer ); +#endif + const RC ret = grb::clear( A ); + return ret == SUCCESS ? UNSUPPORTED : ret; + } + + /** * @} */ diff --git a/include/graphblas/nonblocking/blas3.hpp b/include/graphblas/nonblocking/blas3.hpp index b8e49515b..07e47d00f 100644 --- a/include/graphblas/nonblocking/blas3.hpp +++ b/include/graphblas/nonblocking/blas3.hpp @@ -419,13 +419,14 @@ namespace grb { template< Descriptor descr = descriptors::no_operation, class Operator, - typename InputType1, typename InputType2, typename OutputType, + typename InputType1, typename InputType2, + typename OutputType, typename MaskType, typename Coords, typename RIT, typename CIT, typename NIT > RC maskedOuter( Matrix< OutputType, nonblocking, RIT, CIT, NIT > &A, - const Matrix< OutputType, nonblocking, RIT, CIT, NIT > &mask, + const Matrix< MaskType, nonblocking, RIT, CIT, NIT > &mask, const Vector< InputType1, nonblocking, Coords > &u, const Vector< InputType2, nonblocking, Coords > &v, const Operator &mul = Operator(), @@ -434,6 +435,7 @@ namespace grb { grb::is_operator< Operator >::value && !grb::is_object< InputType1 >::value && !grb::is_object< InputType2 >::value && + !grb::is_object< MaskType >::value && !grb::is_object< OutputType >::value, void >::type * const = nullptr ) { diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index f9116368c..a514f8212 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -924,13 +924,14 @@ namespace grb { template< Descriptor descr = descriptors::no_operation, class Operator, - typename InputType1, typename InputType2, typename OutputType, + typename InputType1, typename InputType2, + typename MaskType, typename OutputType, typename Coords, typename RIT, typename CIT, typename NIT > RC maskedOuter( Matrix< OutputType, reference, RIT, CIT, NIT > &A, - const Matrix< OutputType, reference, RIT, CIT, NIT > &mask, + const Matrix< MaskType, reference, RIT, CIT, NIT > &mask, const Vector< InputType1, reference, Coords > &u, const Vector< InputType2, reference, Coords > &v, const Operator &mul = Operator(), @@ -939,6 +940,7 @@ namespace grb { grb::is_operator< Operator >::value && !grb::is_object< InputType1 >::value && !grb::is_object< InputType2 >::value && + !grb::is_object< MaskType >::value && !grb::is_object< OutputType >::value, void >::type * const = nullptr ) { @@ -958,8 +960,7 @@ namespace grb { ), "grb::maskedOuter", "called with an output matrix that does not match the output domain of " "the given multiplication operator" ); - static_assert( ( !(descr & descriptors::invert_mask) - ), + static_assert( !(descr & descriptors::invert_mask), "grb::maskedOuter: invert_mask descriptor cannot be used "); #ifdef _DEBUG std::cout << "In grb::maskedOuter (reference)\n"; @@ -972,6 +973,7 @@ namespace grb { const size_t n = grb::ncols( mask ); if( m == 0 || n == 0 ) { + // If the mask has a null size, it will be ignored return outer< descr >( A, u, v, mul, phase ); } @@ -1061,7 +1063,7 @@ namespace grb { } } - + nzc = 0; for( size_t i = 0; i < nrows; ++i ) { if( internal::getCoordinates( u ).assigned( i ) ) { @@ -1136,10 +1138,9 @@ namespace grb { } } assert( nzc == old_nzc ); - internal::setCurrentNonzeroes( A, nzc ); - + return ret; } From 28281fd27cf687129c6f550171635d5efeaaa7c7 Mon Sep 17 00:00:00 2001 From: Benjamin Lozes Date: Fri, 16 Feb 2024 10:58:37 +0100 Subject: [PATCH 25/51] hyperdags implementation --- include/graphblas/hyperdags/blas3.hpp | 53 +++++++++++++++++++++++ include/graphblas/hyperdags/hyperdags.hpp | 5 ++- src/graphblas/hyperdags/hyperdags.cpp | 5 ++- 3 files changed, 61 insertions(+), 2 deletions(-) diff --git a/include/graphblas/hyperdags/blas3.hpp b/include/graphblas/hyperdags/blas3.hpp index ee0c10f36..aa74dc6be 100644 --- a/include/graphblas/hyperdags/blas3.hpp +++ b/include/graphblas/hyperdags/blas3.hpp @@ -257,6 +257,59 @@ namespace grb { return ret; } + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType1, typename InputType2, + typename MaskType, typename OutputType, + typename Coords, + typename RIT, typename CIT, typename NIT + > + RC maskedOuter( + Matrix< OutputType, hyperdags, RIT, CIT, NIT > &A, + const Matrix< MaskType, hyperdags, RIT, CIT, NIT > &mask, + const Vector< InputType1, hyperdags, Coords > &u, + const Vector< InputType2, hyperdags, Coords > &v, + const Operator &mul = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + !grb::is_object< MaskType >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + const RC ret = maskedOuter< descr >( + internal::getMatrix( A ), + internal::getMatrix( mask ), + internal::getVector( u ), + internal::getVector( v ), + mul, + phase + ); + if( ret != SUCCESS ) { return ret; } + if( phase != EXECUTE ) { return ret; } + if( nrows( A ) == 0 || ncols( A ) == 0 ) { return ret; } + std::array< const void *, 0 > sourcesP{}; + std::array< uintptr_t, 4 > sourcesC{ + getID( internal::getVector(u) ), + getID( internal::getVector(v) ), + getID( internal::getMatrix(mask) ), + getID( internal::getMatrix(A) ) + }; + std::array< uintptr_t, 1 > destinations{ + getID( internal::getMatrix(A) ) + }; + internal::hyperdags::generator.addOperation( + internal::hyperdags::MASKED_OUTER, + sourcesP.begin(), sourcesP.end(), + sourcesC.begin(), sourcesC.end(), + destinations.begin(), destinations.end() + ); + return ret; + } + template< Descriptor descr = descriptors::no_operation, typename OutputType, typename InputType1, typename InputType2, diff --git a/include/graphblas/hyperdags/hyperdags.hpp b/include/graphblas/hyperdags/hyperdags.hpp index 4ef0e0059..db5d33203 100644 --- a/include/graphblas/hyperdags/hyperdags.hpp +++ b/include/graphblas/hyperdags/hyperdags.hpp @@ -386,6 +386,8 @@ namespace grb { OUTER, + MASKED_OUTER, + UNZIP_VECTOR_VECTOR_VECTOR, ZIP_MATRIX_VECTOR_VECTOR_VECTOR, @@ -493,7 +495,7 @@ namespace grb { }; /** \internal How many operation vertex types exist. */ - const constexpr size_t numOperationVertexTypes = 106; + const constexpr size_t numOperationVertexTypes = 107; /** \internal An array of all operation vertex types. */ const constexpr enum OperationVertexType @@ -553,6 +555,7 @@ namespace grb { MXM_MATRIX_MATRIX_MATRIX_SEMIRING, MXM_MATRIX_MATRIX_MATRIX_MONOID, OUTER, + MASKED_OUTER, UNZIP_VECTOR_VECTOR_VECTOR, ZIP_MATRIX_VECTOR_VECTOR_VECTOR, ZIP_MATRIX_VECTOR_VECTOR, diff --git a/src/graphblas/hyperdags/hyperdags.cpp b/src/graphblas/hyperdags/hyperdags.cpp index 6000f3af7..da1c62b19 100644 --- a/src/graphblas/hyperdags/hyperdags.cpp +++ b/src/graphblas/hyperdags/hyperdags.cpp @@ -219,7 +219,10 @@ std::string grb::internal::hyperdags::toString( return "mxm( matrix, matrix, matrix, monoid, scalar, scalar )"; case OUTER: - return "outer( matrix, vector, vector, scalar, scalar )"; + return "outer( matrix, vector, vector, operation )"; + + case MASKED_OUTER: + return "maskedOuter( matrix, matrix, vector, vector, operation )"; case MXV_VECTOR_VECTOR_MATRIX_VECTOR_VECTOR_R: return "mxv( vector, vector, matrix, vector, vector, ring )"; From 6d3522beea90456833cc131f607a0a6b35a914e0 Mon Sep 17 00:00:00 2001 From: Benjamin Lozes Date: Fri, 16 Feb 2024 11:01:58 +0100 Subject: [PATCH 26/51] bsp1d+hybrid implementation --- include/graphblas/bsp1d/blas3.hpp | 44 +++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/include/graphblas/bsp1d/blas3.hpp b/include/graphblas/bsp1d/blas3.hpp index 386beb164..05f7d3c97 100644 --- a/include/graphblas/bsp1d/blas3.hpp +++ b/include/graphblas/bsp1d/blas3.hpp @@ -205,6 +205,50 @@ namespace grb { return internal::checkGlobalErrorStateOrClear( C, ret ); } + /** \internal Simply delegates to process-local backend */ + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType1, typename InputType2, + typename MaskType, typename OutputType, + typename Coords, + typename RIT, typename CIT, typename NIT + > + RC maskedOuter( + Matrix< OutputType, BSP1D, RIT, CIT, NIT > &A, + const Matrix< MaskType, BSP1D, RIT, CIT, NIT > &mask, + const Vector< InputType1, BSP1D, Coords > &u, + const Vector< InputType2, BSP1D, Coords > &v, + const Operator &mul = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + !grb::is_object< MaskType >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + assert( phase != TRY ); + RC ret = maskedOuter< descr, Operator >( + internal::getLocal( A ), + internal::getLocal( mask ), + internal::getLocal( u ), + internal::getLocal( v ), + mul, + phase + ); + if( phase == RESIZE ) { + if( collectives<>::allreduce( ret, operators::any_or< RC >() ) != SUCCESS ) { + return PANIC; + } else { + return SUCCESS; + } + } + assert( phase == EXECUTE ); + return internal::checkGlobalErrorStateOrClear( A, ret ); + } + } // namespace grb #endif From f15bc07ef190458eb5a99b9b183a72e51e4ade34 Mon Sep 17 00:00:00 2001 From: Benjamin Lozes Date: Fri, 16 Feb 2024 11:02:13 +0100 Subject: [PATCH 27/51] Enable test for all backends --- tests/unit/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index a4b41fdb7..0d61e2521 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -260,7 +260,7 @@ add_grb_executables( outer outer.cpp ) add_grb_executables( maskedOuter maskedOuter.cpp - BACKENDS reference reference_omp + BACKENDS reference reference_omp bsp1d hybrid hyperdags nonblocking ) # must generate the golden output for other tests From 4a171bd0d81a2dc67fb74a97fe62eeaf8dd25edf Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Fri, 16 Feb 2024 00:33:24 +0100 Subject: [PATCH 28/51] rename maskedOuter to outer --- include/graphblas/base/blas3.hpp | 4 ++-- include/graphblas/bsp1d/blas3.hpp | 4 ++-- include/graphblas/hyperdags/blas3.hpp | 4 ++-- include/graphblas/nonblocking/blas3.hpp | 6 +++--- include/graphblas/reference/blas3.hpp | 12 ++++++------ src/graphblas/hyperdags/hyperdags.cpp | 2 +- tests/unit/maskedOuter.cpp | 12 ++++++------ 7 files changed, 22 insertions(+), 22 deletions(-) diff --git a/include/graphblas/base/blas3.hpp b/include/graphblas/base/blas3.hpp index 95b2b439b..611c8c159 100644 --- a/include/graphblas/base/blas3.hpp +++ b/include/graphblas/base/blas3.hpp @@ -499,7 +499,7 @@ namespace grb { typename RIT, typename CIT, typename NIT, Backend backend > - RC maskedOuter( + RC outer( Matrix< OutputType, backend, RIT, CIT, NIT > &A, const Matrix< MaskType, backend, RIT, CIT, NIT > &mask, const Vector< InputType1, backend, Coords > &u, @@ -521,7 +521,7 @@ namespace grb { (void) op; (void) phase; #ifdef _DEBUG - std::cerr << "Selected backend does not implement grb::maskedOuter\n"; + std::cerr << "Selected backend does not implement grb::outer\n"; #endif #ifndef NDEBUG const bool selected_backend_does_not_support_masked_outer = false; diff --git a/include/graphblas/bsp1d/blas3.hpp b/include/graphblas/bsp1d/blas3.hpp index 05f7d3c97..6531dbeea 100644 --- a/include/graphblas/bsp1d/blas3.hpp +++ b/include/graphblas/bsp1d/blas3.hpp @@ -214,7 +214,7 @@ namespace grb { typename Coords, typename RIT, typename CIT, typename NIT > - RC maskedOuter( + RC outer( Matrix< OutputType, BSP1D, RIT, CIT, NIT > &A, const Matrix< MaskType, BSP1D, RIT, CIT, NIT > &mask, const Vector< InputType1, BSP1D, Coords > &u, @@ -230,7 +230,7 @@ namespace grb { void >::type * const = nullptr ) { assert( phase != TRY ); - RC ret = maskedOuter< descr, Operator >( + RC ret = outer< descr, Operator >( internal::getLocal( A ), internal::getLocal( mask ), internal::getLocal( u ), diff --git a/include/graphblas/hyperdags/blas3.hpp b/include/graphblas/hyperdags/blas3.hpp index aa74dc6be..16a52b9b9 100644 --- a/include/graphblas/hyperdags/blas3.hpp +++ b/include/graphblas/hyperdags/blas3.hpp @@ -265,7 +265,7 @@ namespace grb { typename Coords, typename RIT, typename CIT, typename NIT > - RC maskedOuter( + RC outer( Matrix< OutputType, hyperdags, RIT, CIT, NIT > &A, const Matrix< MaskType, hyperdags, RIT, CIT, NIT > &mask, const Vector< InputType1, hyperdags, Coords > &u, @@ -280,7 +280,7 @@ namespace grb { !grb::is_object< OutputType >::value, void >::type * const = nullptr ) { - const RC ret = maskedOuter< descr >( + const RC ret = outer< descr >( internal::getMatrix( A ), internal::getMatrix( mask ), internal::getVector( u ), diff --git a/include/graphblas/nonblocking/blas3.hpp b/include/graphblas/nonblocking/blas3.hpp index 07e47d00f..9cb925efd 100644 --- a/include/graphblas/nonblocking/blas3.hpp +++ b/include/graphblas/nonblocking/blas3.hpp @@ -424,7 +424,7 @@ namespace grb { typename Coords, typename RIT, typename CIT, typename NIT > - RC maskedOuter( + RC outer( Matrix< OutputType, nonblocking, RIT, CIT, NIT > &A, const Matrix< MaskType, nonblocking, RIT, CIT, NIT > &mask, const Vector< InputType1, nonblocking, Coords > &u, @@ -442,7 +442,7 @@ namespace grb { if( internal::NONBLOCKING::warn_if_not_native && config::PIPELINE::warn_if_not_native ) { - std::cerr << "Warning: maskedOuter (nonblocking) currently delegates " + std::cerr << "Warning: outer (nonblocking) currently delegates " << "to a blocking implementation.\n" << " Further similar such warnings will be suppressed.\n"; internal::NONBLOCKING::warn_if_not_native = false; @@ -453,7 +453,7 @@ namespace grb { internal::le.execution(); // second, delegate to the reference backend - return maskedOuter< descr, Operator >( + return outer< descr, Operator >( internal::getRefMatrix( A ), internal::getRefMatrix( mask ), internal::getRefVector( u ), diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index a514f8212..4540b7532 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -929,7 +929,7 @@ namespace grb { typename Coords, typename RIT, typename CIT, typename NIT > - RC maskedOuter( + RC outer( Matrix< OutputType, reference, RIT, CIT, NIT > &A, const Matrix< MaskType, reference, RIT, CIT, NIT > &mask, const Vector< InputType1, reference, Coords > &u, @@ -947,23 +947,23 @@ namespace grb { // static checks NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || std::is_same< typename Operator::D1, InputType1 >::value - ), "grb::maskedOuter", + ), "grb::outer", "called with a prefactor vector that does not match the first domain " "of the given multiplication operator" ); NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || std::is_same< typename Operator::D2, InputType2 >::value - ), "grb::maskedOuter", + ), "grb::outer", "called with a postfactor vector that does not match the first domain " "of the given multiplication operator" ); NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || std::is_same< typename Operator::D3, OutputType >::value - ), "grb::maskedOuter", + ), "grb::outer", "called with an output matrix that does not match the output domain of " "the given multiplication operator" ); static_assert( !(descr & descriptors::invert_mask), - "grb::maskedOuter: invert_mask descriptor cannot be used "); + "grb::outer: invert_mask descriptor cannot be used "); #ifdef _DEBUG - std::cout << "In grb::maskedOuter (reference)\n"; + std::cout << "In grb::outer (reference)\n"; #endif const size_t nrows = size( u ); diff --git a/src/graphblas/hyperdags/hyperdags.cpp b/src/graphblas/hyperdags/hyperdags.cpp index da1c62b19..011857b1f 100644 --- a/src/graphblas/hyperdags/hyperdags.cpp +++ b/src/graphblas/hyperdags/hyperdags.cpp @@ -222,7 +222,7 @@ std::string grb::internal::hyperdags::toString( return "outer( matrix, vector, vector, operation )"; case MASKED_OUTER: - return "maskedOuter( matrix, matrix, vector, vector, operation )"; + return "outer( matrix, matrix, vector, vector, operation )"; case MXV_VECTOR_VECTOR_MATRIX_VECTOR_VECTOR_R: return "mxv( vector, vector, matrix, vector, vector, ring )"; diff --git a/tests/unit/maskedOuter.cpp b/tests/unit/maskedOuter.cpp index ba6d7f8bd..1b4753666 100644 --- a/tests/unit/maskedOuter.cpp +++ b/tests/unit/maskedOuter.cpp @@ -83,8 +83,8 @@ void grbProgram( const void *, const size_t in_size, int &error ) { if( !error ) { - rc = grb::maskedOuter( Result1, Mask1, u, v, ring.getMultiplicativeOperator(), RESIZE ); - rc = rc ? rc : grb::maskedOuter( Result1, Mask1, u, v, ring.getMultiplicativeOperator() ); + rc = grb::outer( Result1, Mask1, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer( Result1, Mask1, u, v, ring.getMultiplicativeOperator() ); if( rc != grb::SUCCESS ) { std::cerr << "Unexpected return code from grb::maskedOuter: " << toString( rc ) << ".\n"; @@ -100,8 +100,8 @@ void grbProgram( const void *, const size_t in_size, int &error ) { } if( !error ) { - rc = grb::maskedOuter< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); - rc = rc ? rc : grb::maskedOuter< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator() ); + rc = grb::outer< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator() ); if( rc != grb::SUCCESS ) { std::cerr << "Unexpected return code from grb::maskedOuter: " << toString( rc ) << ".\n"; @@ -253,8 +253,8 @@ void grb_program_custom_size( const size_t &n, int &error ) { } if( !error ) { - rc = grb::maskedOuter( Result, Mask, u, v, ring.getMultiplicativeOperator(), RESIZE ); - rc = rc ? rc : grb::maskedOuter( Result, Mask, u, v, ring.getMultiplicativeOperator() ); + rc = grb::outer( Result, Mask, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer( Result, Mask, u, v, ring.getMultiplicativeOperator() ); if( rc != grb::SUCCESS ) { std::cerr << "Unexpected return code from grb::maskedOuter: " << toString( rc ) << ".\n"; From 98d62b21e2c6ee9fd6cf1905dbc21b3777b2dc6f Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Fri, 16 Feb 2024 15:28:20 +0100 Subject: [PATCH 29/51] rename --- include/graphblas/base/blas3.hpp | 2 +- tests/unit/maskedOuter.cpp | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/include/graphblas/base/blas3.hpp b/include/graphblas/base/blas3.hpp index 611c8c159..32c8d25ac 100644 --- a/include/graphblas/base/blas3.hpp +++ b/include/graphblas/base/blas3.hpp @@ -445,7 +445,7 @@ namespace grb { /** * Masked outer product of two vectors. Assuming vectors \a u and \a v are * oriented column-wise, the result matrix \a A will contain \f$ uv^T \f$, - * masked to non-zero values from \a mask.. This is an out-of-place + * masked to non-zero values from \a mask. This is an out-of-place * function and will be updated soon to be in-place instead. * * \internal Implemented via mxm as a multiplication of a column vector with diff --git a/tests/unit/maskedOuter.cpp b/tests/unit/maskedOuter.cpp index 1b4753666..6376e465e 100644 --- a/tests/unit/maskedOuter.cpp +++ b/tests/unit/maskedOuter.cpp @@ -86,7 +86,7 @@ void grbProgram( const void *, const size_t in_size, int &error ) { rc = grb::outer( Result1, Mask1, u, v, ring.getMultiplicativeOperator(), RESIZE ); rc = rc ? rc : grb::outer( Result1, Mask1, u, v, ring.getMultiplicativeOperator() ); if( rc != grb::SUCCESS ) { - std::cerr << "Unexpected return code from grb::maskedOuter: " + std::cerr << "Unexpected return code from grb::outer: " << toString( rc ) << ".\n"; error = 15; } @@ -103,7 +103,7 @@ void grbProgram( const void *, const size_t in_size, int &error ) { rc = grb::outer< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); rc = rc ? rc : grb::outer< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator() ); if( rc != grb::SUCCESS ) { - std::cerr << "Unexpected return code from grb::maskedOuter: " + std::cerr << "Unexpected return code from grb::outer: " << toString( rc ) << ".\n"; error = 25; } @@ -256,7 +256,7 @@ void grb_program_custom_size( const size_t &n, int &error ) { rc = grb::outer( Result, Mask, u, v, ring.getMultiplicativeOperator(), RESIZE ); rc = rc ? rc : grb::outer( Result, Mask, u, v, ring.getMultiplicativeOperator() ); if( rc != grb::SUCCESS ) { - std::cerr << "Unexpected return code from grb::maskedOuter: " + std::cerr << "Unexpected return code from grb::outer: " << toString( rc ) << ".\n"; error = 30; } From 7addcf72a550b08b76a5090fea0ae17ff755df8c Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Mon, 19 Feb 2024 15:43:19 +0100 Subject: [PATCH 30/51] descriptors support --- include/graphblas/reference/blas3.hpp | 146 +++++++++++++++++++++----- tests/unit/maskedOuter.cpp | 22 ++-- 2 files changed, 131 insertions(+), 37 deletions(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 4540b7532..b653cb94b 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -960,8 +960,6 @@ namespace grb { ), "grb::outer", "called with an output matrix that does not match the output domain of " "the given multiplication operator" ); - static_assert( !(descr & descriptors::invert_mask), - "grb::outer: invert_mask descriptor cannot be used "); #ifdef _DEBUG std::cout << "In grb::outer (reference)\n"; #endif @@ -996,19 +994,57 @@ namespace grb { const auto &mask_raw = internal::getCRS( mask ); + char * mask_arr = nullptr; + char * mask_buf = nullptr; + OutputType * mask_valbuf = nullptr; + internal::getMatrixBuffers( mask_arr, mask_buf, mask_valbuf, 1, mask ); + + internal::Coordinates< reference > mask_coors; + mask_coors.set( mask_arr, false, mask_buf, ncols ); + size_t nzc = 0; + + if( ( descr & descriptors::structural_complement ) != descriptors::structural_complement ) { + #ifdef _H_GRB_REFERENCE_OMP_BLAS3 - #pragma omp parallel for reduction(+:nzc) + #pragma omp parallel for reduction(+:nzc) #endif - for( size_t i = 0; i < nrows; ++i ) { - if( internal::getCoordinates( u ).assigned( i ) ) { - for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = mask_raw.row_index[ k ]; - if( internal::getCoordinates( v ).assigned( k_col ) ) { - nzc++; + for( size_t i = 0; i < nrows; ++i ) { + if( internal::getCoordinates( u ).assigned( i ) ) { + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + if( + internal::getCoordinates( v ).assigned( k_col ) && + utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) + ) { + + nzc++; + } + } + } + } + + } + else { + + for( size_t i = 0; i < nrows; ++i ) { + if( internal::getCoordinates( u ).assigned( i ) ) { + mask_coors.clear(); + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + mask_coors.assign( k_col ); + } + for( size_t j = 0; j < ncols; ++j ) { + if( + !mask_coors.assign( j ) && + internal::getCoordinates( v ).assigned( j ) + ) { + nzc++; + } } } } + } if( phase == RESIZE ) { @@ -1065,19 +1101,51 @@ namespace grb { nzc = 0; - for( size_t i = 0; i < nrows; ++i ) { - if( internal::getCoordinates( u ).assigned( i ) ) { - for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = mask_raw.row_index[ k ]; - if( internal::getCoordinates( v ).assigned( k_col ) ) { - (void) ++nzc; - if( !crs_only ) { - (void) ++CCS_raw.col_start[ k_col + 1 ]; + + if( ( descr & descriptors::structural_complement ) != descriptors::structural_complement ) { + for( size_t i = 0; i < nrows; ++i ) { + if( internal::getCoordinates( u ).assigned( i ) ) { + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + if( + internal::getCoordinates( v ).assigned( k_col ) && + utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) + ) { + (void) ++nzc; + if( !crs_only ) { + (void) ++CCS_raw.col_start[ k_col + 1 ]; + } } } } + CRS_raw.col_start[ i + 1 ] = nzc; } - CRS_raw.col_start[ i + 1 ] = nzc; + } + else { + + for( size_t i = 0; i < nrows; ++i ) { + if( internal::getCoordinates( u ).assigned( i ) ) { + mask_coors.clear(); + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + mask_coors.assign( k_col ); + } + for( size_t j = 0; j < ncols; ++j ) { + + if( + !mask_coors.assign( j ) && + internal::getCoordinates( v ).assigned( j ) + ) { + + (void) ++nzc; + if( !crs_only ) { + (void) ++CCS_raw.col_start[ j + 1 ]; + } + } + } + } + } + } if( !crs_only ) { @@ -1103,14 +1171,40 @@ namespace grb { for( size_t i = 0; i < nrows; ++i ) { if( internal::getCoordinates( u ).assigned( i ) ) { coors.clear(); - for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = mask_raw.row_index[ k ]; - if( internal::getCoordinates( v ).assigned( k_col ) ) { - coors.assign( k_col ); - grb::apply( valbuf[ k_col ], - x[i], - y[k_col], - mul ); + if( ( descr & descriptors::structural_complement ) != descriptors::structural_complement ) { + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + if( + internal::getCoordinates( v ).assigned( k_col ) && + utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) + ) { + coors.assign( k_col ); + grb::apply( valbuf[ k_col ], + x[i], + y[k_col], + mul ); + } + } + } + else { + + mask_coors.clear(); + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + mask_coors.assign( k_col ); + } + for( size_t j = 0; j < ncols; ++j ) { + + if( + !mask_coors.assign( j ) && + internal::getCoordinates( v ).assigned( j ) + ) { + coors.assign( j ); + grb::apply( valbuf[ j ], + x[i], + y[j], + mul ); + } } } } diff --git a/tests/unit/maskedOuter.cpp b/tests/unit/maskedOuter.cpp index 6376e465e..7440e252f 100644 --- a/tests/unit/maskedOuter.cpp +++ b/tests/unit/maskedOuter.cpp @@ -25,15 +25,15 @@ using namespace grb; // sample data -static const double m1[ 3 ] = { 0.5, 3.4, 5 }; +static const double m1[ 4 ] = { 0.5, 3.4, 5, 0 }; -static const size_t I1[ 3 ] = { 0, 1, 2 }; -static const size_t J1[ 3 ] = { 0, 1, 2 }; +static const size_t I1[ 4 ] = { 0, 1, 2, 0 }; +static const size_t J1[ 4 ] = { 0, 1, 2, 2 }; -static const double m2[ 6 ] = { 0.5, -1, 23, 34, -4.4, 3.14 }; +static const double m2[ 3 ] = { 0.5, -1, 3.14 }; -static const size_t I2[ 6 ] = { 0, 0, 1, 1, 2, 2 }; -static const size_t J2[ 6 ] = { 1, 2, 0, 2, 0, 1 }; +static const size_t I2[ 3 ] = { 0, 1, 2 }; +static const size_t J2[ 3 ] = { 0, 1, 2 }; static const double mask_test1_expect[ 3 ] = { 4, 10, 18 }; static const double mask_test2_expect[ 3 ] = { 11, 20, 27 }; @@ -73,7 +73,7 @@ void grbProgram( const void *, const size_t in_size, int &error ) { } if( !error ) { - rc = grb::buildMatrixUnique( Mask2, &( I2[ 0 ] ), &( J2[ 0 ] ), m2, 6, SEQUENTIAL ); + rc = grb::buildMatrixUnique( Mask2, &( I2[ 0 ] ), &( J2[ 0 ] ), m2, 3, SEQUENTIAL ); if( rc != SUCCESS ) { std::cerr << "\t second mask buildMatrix FAILED\n"; error = 10; @@ -100,8 +100,8 @@ void grbProgram( const void *, const size_t in_size, int &error ) { } if( !error ) { - rc = grb::outer< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); - rc = rc ? rc : grb::outer< descriptors::force_row_major >( Result2, Mask2, u, v, ring.getMultiplicativeOperator() ); + rc = grb::outer< descriptors::force_row_major | descriptors::structural_complement >( Result2, Mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer< descriptors::force_row_major | descriptors::structural_complement >( Result2, Mask2, u, v, ring.getMultiplicativeOperator() ); if( rc != grb::SUCCESS ) { std::cerr << "Unexpected return code from grb::outer: " << toString( rc ) << ".\n"; @@ -253,8 +253,8 @@ void grb_program_custom_size( const size_t &n, int &error ) { } if( !error ) { - rc = grb::outer( Result, Mask, u, v, ring.getMultiplicativeOperator(), RESIZE ); - rc = rc ? rc : grb::outer( Result, Mask, u, v, ring.getMultiplicativeOperator() ); + rc = grb::outer< descriptors::structural >( Result, Mask, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer< descriptors::structural >( Result, Mask, u, v, ring.getMultiplicativeOperator() ); if( rc != grb::SUCCESS ) { std::cerr << "Unexpected return code from grb::outer: " << toString( rc ) << ".\n"; From 296ca8514254749f643930d3dc905677fa9a0486 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Mon, 19 Feb 2024 17:00:36 +0100 Subject: [PATCH 31/51] unmerged conflict --- include/graphblas/base/blas3.hpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/include/graphblas/base/blas3.hpp b/include/graphblas/base/blas3.hpp index 9c26d27b8..32c8d25ac 100644 --- a/include/graphblas/base/blas3.hpp +++ b/include/graphblas/base/blas3.hpp @@ -445,11 +445,7 @@ namespace grb { /** * Masked outer product of two vectors. Assuming vectors \a u and \a v are * oriented column-wise, the result matrix \a A will contain \f$ uv^T \f$, -<<<<<<< HEAD * masked to non-zero values from \a mask. This is an out-of-place -======= - * masked to non-zero values from \a mask.. This is an out-of-place ->>>>>>> 8b4813514f3b87c9a0fe12d0808fd7a6bef80618 * function and will be updated soon to be in-place instead. * * \internal Implemented via mxm as a multiplication of a column vector with From e9012ba5ba3ec957fb511038e619e002a119beee Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 20 Feb 2024 15:56:21 +0100 Subject: [PATCH 32/51] code cleanup --- include/graphblas/reference/blas3.hpp | 63 +++++++++++++-------------- 1 file changed, 31 insertions(+), 32 deletions(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 5ce91e8e5..7cae53f3b 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -1076,17 +1076,9 @@ namespace grb { const InputType1 * __restrict__ const x = internal::getRaw( u ); const InputType2 * __restrict__ const y = internal::getRaw( v ); - char * arr = nullptr; - char * buf = nullptr; - OutputType * valbuf = nullptr; - internal::getMatrixBuffers( arr, buf, valbuf, 1, A ); config::NonzeroIndexType * A_col_index = internal::template getReferenceBuffer< typename config::NonzeroIndexType >( ncols + 1 ); - - internal::Coordinates< reference > coors; - coors.set( arr, false, buf, ncols ); - CRS_raw.col_start[ 0 ] = 0; if( !crs_only ) { @@ -1170,7 +1162,6 @@ namespace grb { nzc = 0; for( size_t i = 0; i < nrows; ++i ) { if( internal::getCoordinates( u ).assigned( i ) ) { - coors.clear(); if( ( descr & descriptors::structural_complement ) != descriptors::structural_complement ) { for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { const size_t k_col = mask_raw.row_index[ k ]; @@ -1178,11 +1169,23 @@ namespace grb { internal::getCoordinates( v ).assigned( k_col ) && utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) ) { - coors.assign( k_col ); - grb::apply( valbuf[ k_col ], - x[i], - y[k_col], + OutputType val; + grb::apply( val, + x[ i ], + y[ k_col ], mul ); + + CRS_raw.row_index[ nzc ] = k_col; + CRS_raw.setValue( nzc, val ); + // update CCS + if( !crs_only ) { + const size_t CCS_index = A_col_index[ k_col ]++ + CCS_raw.col_start[ k_col ]; + CCS_raw.row_index[ CCS_index ] = i; + CCS_raw.setValue( CCS_index, val ); + } + // update count + (void) ++nzc; + } } } @@ -1199,30 +1202,26 @@ namespace grb { !mask_coors.assign( j ) && internal::getCoordinates( v ).assigned( j ) ) { - coors.assign( j ); - grb::apply( valbuf[ j ], - x[i], - y[j], + OutputType val; + grb::apply( val, + x[ i ], + y[ j ], mul ); + + CRS_raw.row_index[ nzc ] = j; + CRS_raw.setValue( nzc, val ); + // update CCS + if( !crs_only ) { + const size_t CCS_index = A_col_index[ j ]++ + CCS_raw.col_start[ j ]; + CCS_raw.row_index[ CCS_index ] = i; + CCS_raw.setValue( CCS_index, val ); + } + // update count + (void) ++nzc; } } } } - for( size_t k = 0; k < coors.nonzeroes(); ++k ) { - assert( nzc < old_nzc ); - const size_t j = coors.index( k ); - // update CRS - CRS_raw.row_index[ nzc ] = j; - CRS_raw.setValue( nzc, valbuf[ j ] ); - // update CCS - if( !crs_only ) { - const size_t CCS_index = A_col_index[ j ]++ + CCS_raw.col_start[ j ]; - CCS_raw.row_index[ CCS_index ] = i; - CCS_raw.setValue( CCS_index, valbuf[ j ] ); - } - // update count - (void) ++nzc; - } CRS_raw.col_start[ i + 1 ] = nzc; } if( !crs_only ) { From 48594d18e68e77b32f529aca1b1d4d8835ddb7d0 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Wed, 21 Feb 2024 14:52:16 +0100 Subject: [PATCH 33/51] masked set matrix to matrix --- include/graphblas/reference/io.hpp | 200 +++++++++++++++++++++++++++++ 1 file changed, 200 insertions(+) diff --git a/include/graphblas/reference/io.hpp b/include/graphblas/reference/io.hpp index 891b13488..acd740760 100644 --- a/include/graphblas/reference/io.hpp +++ b/include/graphblas/reference/io.hpp @@ -1156,6 +1156,206 @@ namespace grb { } } + template< + Descriptor descr = descriptors::no_operation, + typename OutputType, typename MaskType, typename InputType, + typename RIT1, typename CIT1, typename NIT1, + typename RIT2, typename CIT2, typename NIT2 + > + RC set( + Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C, + const Matrix< MaskType, reference, RIT2, CIT2, NIT2 > &Mask, + const Matrix< InputType, reference, RIT2, CIT2, NIT2 > &A, + const Phase &phase = EXECUTE + ) noexcept { + static_assert( !std::is_void< OutputType >::value, + "grb::set (masked set to matrix): cannot have a pattern " + "matrix as output" ); + static_assert( std::is_convertible< InputType, OutputType >::value, + "grb::set (masked set to matrix): input type cannot be " + "converted to output type" + ); + static_assert( + ! ( descr & descriptors::structural && descr & descriptors::invert_mask ), + "grb::set can not be called with both descriptors::structural " + "and descriptors::invert_mask in the masked variant" + ); +#ifdef _DEBUG + std::cout << "Called grb::set (matrix-to-matrix-masked, reference)\n"; +#endif + // static checks + NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) || + std::is_same< InputType, OutputType >::value + ), "grb::set", + "called with non-matching value types" + ); + + // dynamic checks + assert( phase != TRY ); + + const size_t nrows = grb::nrows( C ); + const size_t ncols = grb::ncols( C ); + + const size_t m = grb::nrows( Mask ); + const size_t n = grb::ncols( Mask ); + + /*grb::Monoid< + grb::operators::left_assign< OutputType >, + grb::identities::zero + > dummyMonoid;*/ + + if( m == 0 || n == 0 ) { + // If the mask has a null size, it will be ignored + return set< descr >( C, A, phase ); + } + + if( nrows != grb::nrows( A ) || nrows != m ) { + return MISMATCH; + } + + if( ncols != grb::ncols( A ) || ncols != n ) { + return MISMATCH; + } + + const auto &mask_raw = internal::getCRS( Mask ); + + const auto &A_raw = internal::getCRS( A ); + + size_t nzc = 0; + + char * mask_arr = nullptr; + char * mask_buf = nullptr; + OutputType * mask_valbuf = nullptr; + internal::getMatrixBuffers( mask_arr, mask_buf, mask_valbuf, 1, Mask ); + + internal::Coordinates< reference > mask_coors; + mask_coors.set( mask_arr, false, mask_buf, ncols ); + + for( size_t i = 0; i < nrows; ++i ) { + mask_coors.clear(); + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + if( utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) ) { + mask_coors.assign( k_col ); + } + } + for( auto k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = A_raw.row_index[ k ]; + if( mask_coors.assigned( k_col ) ) { + nzc++; + } + } + } + + if( phase == RESIZE ) { + return resize( C, nzc ); + } + + assert( phase == EXECUTE ); + + if( capacity( C ) < nzc ) { +#ifdef _DEBUG + std::cout << "\t insufficient capacity to complete " + "requested masked set matrix to matrix computation\n"; +#endif + const RC clear_rc = clear( C ); + if( clear_rc != SUCCESS ) { + return PANIC; + } else { + return FAILED; + } + } + + auto &CRS_raw = internal::getCRS( C ); + auto &CCS_raw = internal::getCCS( C ); + + config::NonzeroIndexType * C_col_index = internal::template + getReferenceBuffer< typename config::NonzeroIndexType >( ncols + 1 ); + + CRS_raw.col_start[ 0 ] = 0; + +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for simd +#endif + for( size_t j = 0; j <= ncols; ++j ) { + CCS_raw.col_start[ j ] = 0; + } + + + nzc = 0; + for( size_t i = 0; i < nrows; ++i ) { + mask_coors.clear(); + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + if( utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) ) { + mask_coors.assign( k_col ); + } + } + for( auto k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = A_raw.row_index[ k ]; + if( mask_coors.assigned( k_col ) ) { + nzc++; + CCS_raw.col_start[ k_col + 1 ]++; + } + } + CRS_raw.col_start[ i + 1 ] = nzc; + } + + for( size_t j = 1; j < ncols; ++j ) { + CCS_raw.col_start[ j + 1 ] += CCS_raw.col_start[ j ]; + } + +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for simd +#endif + for( size_t j = 0; j < ncols; ++j ) { + C_col_index[ j ] = 0; + } + + + const size_t old_nzc = nzc; + + // use previously computed CCS offset array to update CCS during the + // computational phase + nzc = 0; + for( size_t i = 0; i < nrows; ++i ) { + mask_coors.clear(); + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + if( utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) ) { + mask_coors.assign( k_col ); + } + } + for( auto k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = A_raw.row_index[ k ]; + if( mask_coors.assigned( k_col ) ) { + OutputType val = A_raw.getValue( k, *( (OutputType*) nullptr ) ); + CRS_raw.row_index[ nzc ] = k_col; + CRS_raw.setValue( nzc, val ); + const size_t CCS_index = C_col_index[ k_col ]++ + CCS_raw.col_start[ k_col ]; + CCS_raw.row_index[ CCS_index ] = i; + CCS_raw.setValue( CCS_index, val ); + nzc++; + } + } + } + + for( size_t j = 0; j < ncols; ++j ) { + assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] == + C_col_index[ j ] ); + } + assert( old_nzc == nzc ); + + internal::setCurrentNonzeroes( C, nzc ); + + return SUCCESS; + } + + + + + + /** * Ingests raw data into a GraphBLAS vector. * From 0ef5c3d1df13cf7708c8b8ed6d6d372f0323cbe0 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Wed, 21 Feb 2024 14:52:25 +0100 Subject: [PATCH 34/51] tests --- tests/unit/matrixSet.cpp | 120 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 120 insertions(+) diff --git a/tests/unit/matrixSet.cpp b/tests/unit/matrixSet.cpp index dfa9c0c8c..7b205ad98 100644 --- a/tests/unit/matrixSet.cpp +++ b/tests/unit/matrixSet.cpp @@ -201,6 +201,126 @@ void grb_program( const size_t &n, grb::RC &rc ) { if( rc != SUCCESS ) { return; } + + //check masked matrix set + grb::Matrix< int > Mask( n, n ); + grb::Matrix< int > Output( n, n ); + grb::Matrix< int > Input( n, n ); + size_t I[ 2 * n - 1 ], J[ 2 * n - 1 ]; + int mask_vals [ 2 * n - 1 ]; + int input_vals [ 2 * n - 1 ]; + + for( size_t k = 0; k < n; ++k ) { + I[ k ] = J[ k ] = k; + mask_vals[ k ] = 1; + input_vals[ k ] = k; + if( k < n - 1 ) { + I[ n + k ] = k; + J[ n + k ] = k + 1; + mask_vals [ n + k ] = 0; + input_vals[ n + k ] = k; + } + } + + rc = grb::buildMatrixUnique( Mask, I, J, mask_vals, 2 * n - 1, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildMatrixUnique of mask matrix FAILED\n"; + return; + } + + rc = grb::buildMatrixUnique( Input, I, J, input_vals, 2 * n - 1, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildMatrixUnique of input matrix FAILED\n"; + return; + } + + rc = grb::resize(Output, 2 * n - 1 ); + if( rc != SUCCESS ) { + std::cerr << "\t grb::resize of output matrix for the matrix to matrix masked set FAILED\n"; + return; + } + + rc = grb::set< descriptors::structural >( Output, Mask, Input ); + if( rc != SUCCESS ) { + std::cerr << "\t grb::set structural (matrix to matrix masked) FAILED\n"; + return; + } + + if( grb::nnz( Output ) != 2 * n - 1 ) { + std::cerr << "\t unexpected number of output elements ( " << grb::nnz( Output ) << " ), expected " << 2 * n - 1 <<".\n"; + rc = FAILED; + } + + for( const auto & triplet : Output ) { + if( triplet.first.first != triplet.first.second && triplet.first.first != triplet.first.second - 1 ) { + std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ), value " << triplet.second << ".\n"; + rc = FAILED; + } if( triplet.first.first != triplet.second ) { + std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ) with value " << triplet.second; + std::cerr << ", expected value "<< triplet.first.first <<".\n"; + rc = FAILED; + } + } + + if( rc != SUCCESS ) { + return; + } + + + + rc = grb::set( Output, Mask, Input ); + if( rc != SUCCESS ) { + std::cerr << "\t grb::set (matrix to matrix masked) FAILED\n"; + return; + } + + if( grb::nnz( Output ) != n ) { + std::cerr << "\t unexpected number of output elements ( " << grb::nnz( Output ) << " ), expected " << n <<".\n"; + rc = FAILED; + } + + for( const auto & triplet : Output ) { + if( triplet.first.first != triplet.first.second ) { + std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ), value " << triplet.second << ".\n"; + rc = FAILED; + } if( triplet.first.first != triplet.second ) { + std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ) with value " << triplet.second; + std::cerr << ", expected value "<< triplet.first.first <<".\n"; + rc = FAILED; + } + } + + if( rc != SUCCESS ) { + return; + } + + + + rc = grb::set< descriptors::invert_mask >( Output, Mask, Input ); + if( rc != SUCCESS ) { + std::cerr << "\t grb::set invert mask (matrix to matrix masked) FAILED\n"; + return; + } + + if( grb::nnz( Output ) != n - 1 ) { + std::cerr << "\t unexpected number of output elements ( " << grb::nnz( Output ) << " ), expected " << n - 1 <<".\n"; + rc = FAILED; + } + + for( const auto & triplet : Output ) { + if( triplet.first.first != triplet.first.second - 1 ) { + std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ), value " << triplet.second << ".\n"; + rc = FAILED; + } if( triplet.first.first != triplet.second ) { + std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ) with value " << triplet.second; + std::cerr << ", expected value "<< triplet.first.first <<".\n"; + rc = FAILED; + } + } + + if( rc != SUCCESS ) { + return; + } } int main( int argc, char ** argv ) { From cabc337432b4d7d66d39947fd2798b80325cc733 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Wed, 21 Feb 2024 15:14:13 +0100 Subject: [PATCH 35/51] tests refactor --- tests/unit/matrixSet.cpp | 80 ++++++++++++++++++++-------------------- 1 file changed, 41 insertions(+), 39 deletions(-) diff --git a/tests/unit/matrixSet.cpp b/tests/unit/matrixSet.cpp index 7b205ad98..40ee0a69e 100644 --- a/tests/unit/matrixSet.cpp +++ b/tests/unit/matrixSet.cpp @@ -43,6 +43,12 @@ void grb_program( const size_t &n, grb::RC &rc ) { grb::Matrix< void > C( n, n ); grb::Matrix< void > D( n, n ); grb::Matrix< unsigned int > E( n, n ); + grb::Matrix< int > Mask( n, n ); + grb::Matrix< int > Output( n, n ); + grb::Matrix< int > Input( n, n ); + + + rc = grb::resize( A, 15 ); if( rc == SUCCESS ) { rc = grb::buildMatrixUnique( A, I, J, data1, 15, SEQUENTIAL ); @@ -61,6 +67,35 @@ void grb_program( const size_t &n, grb::RC &rc ) { } } } + + size_t I_mask[ 2 * n - 1 ], J_mask[ 2 * n - 1 ]; + int mask_vals [ 2 * n - 1 ]; + int input_vals [ 2 * n - 1 ]; + + for( size_t k = 0; k < n; ++k ) { + I_mask[ k ] = J_mask[ k ] = k; + mask_vals[ k ] = 1; + input_vals[ k ] = k; + if( k < n - 1 ) { + I_mask[ n + k ] = k; + J_mask[ n + k ] = k + 1; + mask_vals [ n + k ] = 0; + input_vals[ n + k ] = k; + } + } + + rc = grb::buildMatrixUnique( Mask, I_mask, J_mask, mask_vals, 2 * n - 1, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildMatrixUnique of mask matrix FAILED\n"; + return; + } + + rc = grb::buildMatrixUnique( Input, I_mask, J_mask, input_vals, 2 * n - 1, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t buildMatrixUnique of input matrix FAILED\n"; + return; + } + if( rc == SUCCESS ) { rc = grb::resize( B, 15 ); } @@ -73,6 +108,9 @@ void grb_program( const size_t &n, grb::RC &rc ) { if( rc == SUCCESS ) { rc = grb::resize( E, 15 ); } + if( rc == SUCCESS ) { + rc = grb::resize(Output, 2 * n - 1 ); + } if( rc != SUCCESS || grb::nnz( A ) != 15 ) { std::cerr << "\tinitialisation FAILED\n"; return; @@ -203,42 +241,6 @@ void grb_program( const size_t &n, grb::RC &rc ) { } //check masked matrix set - grb::Matrix< int > Mask( n, n ); - grb::Matrix< int > Output( n, n ); - grb::Matrix< int > Input( n, n ); - size_t I[ 2 * n - 1 ], J[ 2 * n - 1 ]; - int mask_vals [ 2 * n - 1 ]; - int input_vals [ 2 * n - 1 ]; - - for( size_t k = 0; k < n; ++k ) { - I[ k ] = J[ k ] = k; - mask_vals[ k ] = 1; - input_vals[ k ] = k; - if( k < n - 1 ) { - I[ n + k ] = k; - J[ n + k ] = k + 1; - mask_vals [ n + k ] = 0; - input_vals[ n + k ] = k; - } - } - - rc = grb::buildMatrixUnique( Mask, I, J, mask_vals, 2 * n - 1, SEQUENTIAL ); - if( rc != SUCCESS ) { - std::cerr << "\t buildMatrixUnique of mask matrix FAILED\n"; - return; - } - - rc = grb::buildMatrixUnique( Input, I, J, input_vals, 2 * n - 1, SEQUENTIAL ); - if( rc != SUCCESS ) { - std::cerr << "\t buildMatrixUnique of input matrix FAILED\n"; - return; - } - - rc = grb::resize(Output, 2 * n - 1 ); - if( rc != SUCCESS ) { - std::cerr << "\t grb::resize of output matrix for the matrix to matrix masked set FAILED\n"; - return; - } rc = grb::set< descriptors::structural >( Output, Mask, Input ); if( rc != SUCCESS ) { @@ -255,7 +257,7 @@ void grb_program( const size_t &n, grb::RC &rc ) { if( triplet.first.first != triplet.first.second && triplet.first.first != triplet.first.second - 1 ) { std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ), value " << triplet.second << ".\n"; rc = FAILED; - } if( triplet.first.first != triplet.second ) { + } if( (int) triplet.first.first != triplet.second ) { std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ) with value " << triplet.second; std::cerr << ", expected value "<< triplet.first.first <<".\n"; rc = FAILED; @@ -283,7 +285,7 @@ void grb_program( const size_t &n, grb::RC &rc ) { if( triplet.first.first != triplet.first.second ) { std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ), value " << triplet.second << ".\n"; rc = FAILED; - } if( triplet.first.first != triplet.second ) { + } if( (int) triplet.first.first != triplet.second ) { std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ) with value " << triplet.second; std::cerr << ", expected value "<< triplet.first.first <<".\n"; rc = FAILED; @@ -311,7 +313,7 @@ void grb_program( const size_t &n, grb::RC &rc ) { if( triplet.first.first != triplet.first.second - 1 ) { std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ), value " << triplet.second << ".\n"; rc = FAILED; - } if( triplet.first.first != triplet.second ) { + } if( (int) triplet.first.first != triplet.second ) { std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ) with value " << triplet.second; std::cerr << ", expected value "<< triplet.first.first <<".\n"; rc = FAILED; From 4040f69cab445451fb32670395a537addf2fd63c Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Wed, 21 Feb 2024 15:14:42 +0100 Subject: [PATCH 36/51] nonblocking implementation --- include/graphblas/nonblocking/io.hpp | 35 ++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/include/graphblas/nonblocking/io.hpp b/include/graphblas/nonblocking/io.hpp index 64e53172d..cabc006c3 100644 --- a/include/graphblas/nonblocking/io.hpp +++ b/include/graphblas/nonblocking/io.hpp @@ -1217,6 +1217,41 @@ namespace grb { } } + + template< + Descriptor descr = descriptors::no_operation, + typename OutputType, typename MaskType, typename InputType, + typename RIT1, typename CIT1, typename NIT1, + typename RIT2, typename CIT2, typename NIT2 + > + RC set( + Matrix< OutputType, nonblocking, RIT1, CIT1, NIT1 > &C, + const Matrix< MaskType, nonblocking, RIT2, CIT2, NIT2 > &Mask, + const Matrix< InputType, nonblocking, RIT2, CIT2, NIT2 > &A, + const Phase &phase = EXECUTE + ) noexcept { + if( internal::NONBLOCKING::warn_if_not_native && + config::PIPELINE::warn_if_not_native + ) { + std::cerr << "Warning: masked set (nonblocking) currently delegates " + << "to a blocking implementation.\n" + << " Further similar such warnings will be suppressed.\n"; + internal::NONBLOCKING::warn_if_not_native = false; + } + + // nonblocking execution is not supported + // first, execute any computation that is not completed + internal::le.execution(); + + // second, delegate to the reference backend + return set< descr >( + internal::getRefMatrix( C ), + internal::getRefMatrix( Mask ), + internal::getRefMatrix( A ), + phase + ); + } + template< Descriptor descr = descriptors::no_operation, typename InputType, From a5706905e489708e72d24200a1e59b821587989a Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Wed, 21 Feb 2024 15:15:09 +0100 Subject: [PATCH 37/51] hyperdags implementation --- include/graphblas/hyperdags/hyperdags.hpp | 4 ++- include/graphblas/hyperdags/io.hpp | 35 +++++++++++++++++++++++ src/graphblas/hyperdags/hyperdags.cpp | 3 ++ 3 files changed, 41 insertions(+), 1 deletion(-) diff --git a/include/graphblas/hyperdags/hyperdags.hpp b/include/graphblas/hyperdags/hyperdags.hpp index 4ef0e0059..5bfb3812a 100644 --- a/include/graphblas/hyperdags/hyperdags.hpp +++ b/include/graphblas/hyperdags/hyperdags.hpp @@ -380,6 +380,8 @@ namespace grb { SET_MATRIX_MATRIX_INPUT2, + SET_MATRIX_MATRIX_MASKED, + MXM_MATRIX_MATRIX_MATRIX_SEMIRING, MXM_MATRIX_MATRIX_MATRIX_MONOID, @@ -493,7 +495,7 @@ namespace grb { }; /** \internal How many operation vertex types exist. */ - const constexpr size_t numOperationVertexTypes = 106; + const constexpr size_t numOperationVertexTypes = 107; /** \internal An array of all operation vertex types. */ const constexpr enum OperationVertexType diff --git a/include/graphblas/hyperdags/io.hpp b/include/graphblas/hyperdags/io.hpp index 09b16634a..7e420e444 100644 --- a/include/graphblas/hyperdags/io.hpp +++ b/include/graphblas/hyperdags/io.hpp @@ -403,6 +403,41 @@ namespace grb { return ret; } + template< + Descriptor descr = descriptors::no_operation, + typename OutputType, typename MaskType, typename InputType, + typename RIT1, typename CIT1, typename NIT1, + typename RIT2, typename CIT2, typename NIT2 + > + RC set( + Matrix< OutputType, hyperdags, RIT1, CIT1, NIT1 > &C, + const Matrix< MaskType, hyperdags, RIT2, CIT2, NIT2 > &Mask, + const Matrix< InputType, hyperdags, RIT2, CIT2, NIT2 > &A, + const Phase &phase = EXECUTE + ) { + const RC ret = set< descr >( + internal::getMatrix( C ), internal::getMatrix( Mask ), + internal::getMatrix( A ), phase + ); + if( ret != SUCCESS ) { return ret; } + if( phase != EXECUTE ) { return ret; } + if( nrows( A ) == 0 || ncols( A ) == 0 ) { return ret; } + std::array< const void *, 0 > sourcesP{}; + std::array< uintptr_t, 3 > sourcesC{ + getID( internal::getMatrix(A) ), + getID( internal::getMatrix(Mask) ), + getID( internal::getMatrix(C) ) + }; + std::array< uintptr_t, 1 > destinations{ getID( internal::getMatrix(C) ) }; + internal::hyperdags::generator.addOperation( + internal::hyperdags::SET_MATRIX_MATRIX_MASKED, + sourcesP.begin(), sourcesP.end(), + sourcesC.begin(), sourcesC.end(), + destinations.begin(), destinations.end() + ); + return ret; + } + template< typename DataType, typename Coords > RC clear( Vector< DataType, hyperdags, Coords > &x ) { const RC ret = clear( internal::getVector( x ) ); diff --git a/src/graphblas/hyperdags/hyperdags.cpp b/src/graphblas/hyperdags/hyperdags.cpp index 6000f3af7..ff4c29235 100644 --- a/src/graphblas/hyperdags/hyperdags.cpp +++ b/src/graphblas/hyperdags/hyperdags.cpp @@ -215,6 +215,9 @@ std::string grb::internal::hyperdags::toString( case SET_MATRIX_MATRIX_INPUT2: return "set( matrix, matrix, scalar )"; + case SET_MATRIX_MATRIX_MASKED: + return "set( matrix, matrix, matrix )"; + case MXM_MATRIX_MATRIX_MATRIX_MONOID: return "mxm( matrix, matrix, matrix, monoid, scalar, scalar )"; From da9de374ec2913818cb5e6ac77de661f733ae178 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Sat, 24 Feb 2024 09:53:43 +0100 Subject: [PATCH 38/51] unsupport structural complement --- include/graphblas/reference/blas3.hpp | 178 +++++++------------------- tests/unit/maskedOuter.cpp | 12 +- 2 files changed, 53 insertions(+), 137 deletions(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 7cae53f3b..3f059f97d 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -960,6 +960,11 @@ namespace grb { ), "grb::outer", "called with an output matrix that does not match the output domain of " "the given multiplication operator" ); + static_assert( + ! ( descr & descriptors::structural && descr & descriptors::invert_mask ), + "grb::outer can not be called with both descriptors::structural " + "and descriptors::invert_mask in the masked variant" + ); #ifdef _DEBUG std::cout << "In grb::outer (reference)\n"; #endif @@ -1004,47 +1009,22 @@ namespace grb { size_t nzc = 0; - if( ( descr & descriptors::structural_complement ) != descriptors::structural_complement ) { - #ifdef _H_GRB_REFERENCE_OMP_BLAS3 - #pragma omp parallel for reduction(+:nzc) + #pragma omp parallel for reduction(+:nzc) #endif - for( size_t i = 0; i < nrows; ++i ) { - if( internal::getCoordinates( u ).assigned( i ) ) { - for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = mask_raw.row_index[ k ]; - if( - internal::getCoordinates( v ).assigned( k_col ) && - utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) - ) { - - nzc++; - } - } - } - } - - } - else { + for( size_t i = 0; i < nrows; ++i ) { + if( internal::getCoordinates( u ).assigned( i ) ) { + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + if( + internal::getCoordinates( v ).assigned( k_col ) && + utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) + ) { - for( size_t i = 0; i < nrows; ++i ) { - if( internal::getCoordinates( u ).assigned( i ) ) { - mask_coors.clear(); - for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = mask_raw.row_index[ k ]; - mask_coors.assign( k_col ); - } - for( size_t j = 0; j < ncols; ++j ) { - if( - !mask_coors.assign( j ) && - internal::getCoordinates( v ).assigned( j ) - ) { - nzc++; - } + nzc++; } } } - } if( phase == RESIZE ) { @@ -1094,50 +1074,22 @@ namespace grb { nzc = 0; - if( ( descr & descriptors::structural_complement ) != descriptors::structural_complement ) { - for( size_t i = 0; i < nrows; ++i ) { - if( internal::getCoordinates( u ).assigned( i ) ) { - for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = mask_raw.row_index[ k ]; - if( - internal::getCoordinates( v ).assigned( k_col ) && - utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) - ) { - (void) ++nzc; - if( !crs_only ) { - (void) ++CCS_raw.col_start[ k_col + 1 ]; - } - } - } - } - CRS_raw.col_start[ i + 1 ] = nzc; - } - } - else { - - for( size_t i = 0; i < nrows; ++i ) { - if( internal::getCoordinates( u ).assigned( i ) ) { - mask_coors.clear(); - for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = mask_raw.row_index[ k ]; - mask_coors.assign( k_col ); - } - for( size_t j = 0; j < ncols; ++j ) { - - if( - !mask_coors.assign( j ) && - internal::getCoordinates( v ).assigned( j ) - ) { - - (void) ++nzc; - if( !crs_only ) { - (void) ++CCS_raw.col_start[ j + 1 ]; - } + for( size_t i = 0; i < nrows; ++i ) { + if( internal::getCoordinates( u ).assigned( i ) ) { + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + if( + internal::getCoordinates( v ).assigned( k_col ) && + utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) + ) { + (void) ++nzc; + if( !crs_only ) { + (void) ++CCS_raw.col_start[ k_col + 1 ]; } } } - CRS_raw.col_start[ i + 1 ] = nzc; } + CRS_raw.col_start[ i + 1 ] = nzc; } if( !crs_only ) { @@ -1162,63 +1114,27 @@ namespace grb { nzc = 0; for( size_t i = 0; i < nrows; ++i ) { if( internal::getCoordinates( u ).assigned( i ) ) { - if( ( descr & descriptors::structural_complement ) != descriptors::structural_complement ) { - for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = mask_raw.row_index[ k ]; - if( - internal::getCoordinates( v ).assigned( k_col ) && - utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) - ) { - OutputType val; - grb::apply( val, - x[ i ], - y[ k_col ], - mul ); - - CRS_raw.row_index[ nzc ] = k_col; - CRS_raw.setValue( nzc, val ); - // update CCS - if( !crs_only ) { - const size_t CCS_index = A_col_index[ k_col ]++ + CCS_raw.col_start[ k_col ]; - CCS_raw.row_index[ CCS_index ] = i; - CCS_raw.setValue( CCS_index, val ); - } - // update count - (void) ++nzc; - - } - } - } - else { - - mask_coors.clear(); - for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = mask_raw.row_index[ k ]; - mask_coors.assign( k_col ); - } - for( size_t j = 0; j < ncols; ++j ) { - - if( - !mask_coors.assign( j ) && - internal::getCoordinates( v ).assigned( j ) - ) { - OutputType val; - grb::apply( val, - x[ i ], - y[ j ], - mul ); - - CRS_raw.row_index[ nzc ] = j; - CRS_raw.setValue( nzc, val ); - // update CCS - if( !crs_only ) { - const size_t CCS_index = A_col_index[ j ]++ + CCS_raw.col_start[ j ]; - CCS_raw.row_index[ CCS_index ] = i; - CCS_raw.setValue( CCS_index, val ); - } - // update count - (void) ++nzc; + for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = mask_raw.row_index[ k ]; + if( + internal::getCoordinates( v ).assigned( k_col ) && + utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) + ) { + OutputType val; + grb::apply( val, + x[ i ], + y[ k_col ], + mul ); + CRS_raw.row_index[ nzc ] = k_col; + CRS_raw.setValue( nzc, val ); + // update CCS + if( !crs_only ) { + const size_t CCS_index = A_col_index[ k_col ]++ + CCS_raw.col_start[ k_col ]; + CCS_raw.row_index[ CCS_index ] = i; + CCS_raw.setValue( CCS_index, val ); } + // update count + (void) ++nzc; } } } diff --git a/tests/unit/maskedOuter.cpp b/tests/unit/maskedOuter.cpp index 7440e252f..31ab6ab63 100644 --- a/tests/unit/maskedOuter.cpp +++ b/tests/unit/maskedOuter.cpp @@ -30,10 +30,10 @@ static const double m1[ 4 ] = { 0.5, 3.4, 5, 0 }; static const size_t I1[ 4 ] = { 0, 1, 2, 0 }; static const size_t J1[ 4 ] = { 0, 1, 2, 2 }; -static const double m2[ 3 ] = { 0.5, -1, 3.14 }; +static const double m2[ 8 ] = { 1, 1, 0, 0, 0, 0, 0, 0 }; -static const size_t I2[ 3 ] = { 0, 1, 2 }; -static const size_t J2[ 3 ] = { 0, 1, 2 }; +static const size_t I2[ 8 ] = { 0, 2, 0, 0, 1, 1, 2, 2 }; +static const size_t J2[ 8 ] = { 0, 2, 1, 2, 0, 2, 0, 1 }; static const double mask_test1_expect[ 3 ] = { 4, 10, 18 }; static const double mask_test2_expect[ 3 ] = { 11, 20, 27 }; @@ -73,7 +73,7 @@ void grbProgram( const void *, const size_t in_size, int &error ) { } if( !error ) { - rc = grb::buildMatrixUnique( Mask2, &( I2[ 0 ] ), &( J2[ 0 ] ), m2, 3, SEQUENTIAL ); + rc = grb::buildMatrixUnique( Mask2, &( I2[ 0 ] ), &( J2[ 0 ] ), m2, 8, SEQUENTIAL ); if( rc != SUCCESS ) { std::cerr << "\t second mask buildMatrix FAILED\n"; error = 10; @@ -100,8 +100,8 @@ void grbProgram( const void *, const size_t in_size, int &error ) { } if( !error ) { - rc = grb::outer< descriptors::force_row_major | descriptors::structural_complement >( Result2, Mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); - rc = rc ? rc : grb::outer< descriptors::force_row_major | descriptors::structural_complement >( Result2, Mask2, u, v, ring.getMultiplicativeOperator() ); + rc = grb::outer< descriptors::force_row_major | descriptors::invert_mask >( Result2, Mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer< descriptors::force_row_major | descriptors::invert_mask >( Result2, Mask2, u, v, ring.getMultiplicativeOperator() ); if( rc != grb::SUCCESS ) { std::cerr << "Unexpected return code from grb::outer: " << toString( rc ) << ".\n"; From 72892482e3b079301cb9a3f35af6f934944605a0 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Mon, 26 Feb 2024 11:31:12 +0100 Subject: [PATCH 39/51] fix print of the unit test --- tests/unit/maskedOuter.cpp | 31 +++++++++++++++---------------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/tests/unit/maskedOuter.cpp b/tests/unit/maskedOuter.cpp index 31ab6ab63..3bbc0cfe7 100644 --- a/tests/unit/maskedOuter.cpp +++ b/tests/unit/maskedOuter.cpp @@ -305,20 +305,6 @@ int main( int argc, char ** argv ) { (void)argc; std::cout << "Functional test executable: " << argv[ 0 ] << "\n"; - int error; - grb::Launcher< AUTOMATIC > launcher; - - if( launcher.exec( &grbProgram, nullptr, 0, error ) != SUCCESS ) { - std::cerr << "Test 1 failed to launch\n"; - error = 255; - } - if( error == 0 ) { - std::cout << "Test 1 OK\n" << std::endl; - } else { - std::cerr << std::flush; - std::cout << "Test 1 FAILED\n" << std::endl; - } - bool printUsage = false; size_t n = 100; @@ -347,14 +333,27 @@ int main( int argc, char ** argv ) { return 1; } - error = 0; + + + int error = 0; + grb::Launcher< AUTOMATIC > launcher; + + if( launcher.exec( &grbProgram, nullptr, 0, error ) != SUCCESS ) { + std::cerr << "Test 1 failed to launch\n"; + error = 255; + } + if( error != 0 ) { + std::cerr << std::flush; + std::cout << "Test 1 FAILED\n" << std::endl; + return 0; + } if( launcher.exec( &grb_program_custom_size, n, error ) != SUCCESS ) { std::cerr << "Launching test 2 FAILED\n"; error = 255; } if( error == 0 ) { - std::cout << "Test 2 OK\n" << std::endl; + std::cout << "Test OK\n" << std::endl; } else { std::cerr << std::flush; std::cout << "Test 2 FAILED\n" << std::endl; From c410d8f171c0af295cf6baaa606b6318a201f7e5 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Mon, 26 Feb 2024 11:31:31 +0100 Subject: [PATCH 40/51] code cleanup --- include/graphblas/reference/blas3.hpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 3f059f97d..3faee95b0 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -1106,9 +1106,6 @@ namespace grb { } } - - const size_t old_nzc = nzc; - // use previously computed CCS offset array to update CCS during the // computational phase nzc = 0; @@ -1146,7 +1143,6 @@ namespace grb { A_col_index[ j ] ); } } - assert( nzc == old_nzc ); internal::setCurrentNonzeroes( A, nzc ); From b06362484d0d2d75589a30e55e60ff2ece047d3b Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Mon, 26 Feb 2024 11:50:08 +0100 Subject: [PATCH 41/51] unit test fix --- tests/unit/maskedOuter.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit/maskedOuter.cpp b/tests/unit/maskedOuter.cpp index 3bbc0cfe7..7cfc8756c 100644 --- a/tests/unit/maskedOuter.cpp +++ b/tests/unit/maskedOuter.cpp @@ -345,7 +345,7 @@ int main( int argc, char ** argv ) { if( error != 0 ) { std::cerr << std::flush; std::cout << "Test 1 FAILED\n" << std::endl; - return 0; + return error; } if( launcher.exec( &grb_program_custom_size, n, error ) != SUCCESS ) { From 3de4e27a91e69229927730bd3bece263db1e76f1 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Mon, 26 Feb 2024 12:25:56 +0100 Subject: [PATCH 42/51] rename --- include/graphblas/hyperdags/io.hpp | 6 ++--- include/graphblas/nonblocking/io.hpp | 4 ++-- include/graphblas/reference/io.hpp | 13 ++++------ tests/unit/matrixSet.cpp | 36 ++++++++++++++-------------- 4 files changed, 28 insertions(+), 31 deletions(-) diff --git a/include/graphblas/hyperdags/io.hpp b/include/graphblas/hyperdags/io.hpp index 7e420e444..aff441bc3 100644 --- a/include/graphblas/hyperdags/io.hpp +++ b/include/graphblas/hyperdags/io.hpp @@ -411,12 +411,12 @@ namespace grb { > RC set( Matrix< OutputType, hyperdags, RIT1, CIT1, NIT1 > &C, - const Matrix< MaskType, hyperdags, RIT2, CIT2, NIT2 > &Mask, + const Matrix< MaskType, hyperdags, RIT2, CIT2, NIT2 > &M, const Matrix< InputType, hyperdags, RIT2, CIT2, NIT2 > &A, const Phase &phase = EXECUTE ) { const RC ret = set< descr >( - internal::getMatrix( C ), internal::getMatrix( Mask ), + internal::getMatrix( C ), internal::getMatrix( M ), internal::getMatrix( A ), phase ); if( ret != SUCCESS ) { return ret; } @@ -425,7 +425,7 @@ namespace grb { std::array< const void *, 0 > sourcesP{}; std::array< uintptr_t, 3 > sourcesC{ getID( internal::getMatrix(A) ), - getID( internal::getMatrix(Mask) ), + getID( internal::getMatrix(M) ), getID( internal::getMatrix(C) ) }; std::array< uintptr_t, 1 > destinations{ getID( internal::getMatrix(C) ) }; diff --git a/include/graphblas/nonblocking/io.hpp b/include/graphblas/nonblocking/io.hpp index cabc006c3..a429d8dd2 100644 --- a/include/graphblas/nonblocking/io.hpp +++ b/include/graphblas/nonblocking/io.hpp @@ -1226,7 +1226,7 @@ namespace grb { > RC set( Matrix< OutputType, nonblocking, RIT1, CIT1, NIT1 > &C, - const Matrix< MaskType, nonblocking, RIT2, CIT2, NIT2 > &Mask, + const Matrix< MaskType, nonblocking, RIT2, CIT2, NIT2 > &M, const Matrix< InputType, nonblocking, RIT2, CIT2, NIT2 > &A, const Phase &phase = EXECUTE ) noexcept { @@ -1246,7 +1246,7 @@ namespace grb { // second, delegate to the reference backend return set< descr >( internal::getRefMatrix( C ), - internal::getRefMatrix( Mask ), + internal::getRefMatrix( M ), internal::getRefMatrix( A ), phase ); diff --git a/include/graphblas/reference/io.hpp b/include/graphblas/reference/io.hpp index acd740760..c1c524725 100644 --- a/include/graphblas/reference/io.hpp +++ b/include/graphblas/reference/io.hpp @@ -1164,7 +1164,7 @@ namespace grb { > RC set( Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C, - const Matrix< MaskType, reference, RIT2, CIT2, NIT2 > &Mask, + const Matrix< MaskType, reference, RIT2, CIT2, NIT2 > &M, const Matrix< InputType, reference, RIT2, CIT2, NIT2 > &A, const Phase &phase = EXECUTE ) noexcept { @@ -1196,8 +1196,8 @@ namespace grb { const size_t nrows = grb::nrows( C ); const size_t ncols = grb::ncols( C ); - const size_t m = grb::nrows( Mask ); - const size_t n = grb::ncols( Mask ); + const size_t m = grb::nrows( M ); + const size_t n = grb::ncols( M ); /*grb::Monoid< grb::operators::left_assign< OutputType >, @@ -1217,7 +1217,7 @@ namespace grb { return MISMATCH; } - const auto &mask_raw = internal::getCRS( Mask ); + const auto &mask_raw = internal::getCRS( M ); const auto &A_raw = internal::getCRS( A ); @@ -1226,7 +1226,7 @@ namespace grb { char * mask_arr = nullptr; char * mask_buf = nullptr; OutputType * mask_valbuf = nullptr; - internal::getMatrixBuffers( mask_arr, mask_buf, mask_valbuf, 1, Mask ); + internal::getMatrixBuffers( mask_arr, mask_buf, mask_valbuf, 1, M ); internal::Coordinates< reference > mask_coors; mask_coors.set( mask_arr, false, mask_buf, ncols ); @@ -1313,8 +1313,6 @@ namespace grb { } - const size_t old_nzc = nzc; - // use previously computed CCS offset array to update CCS during the // computational phase nzc = 0; @@ -1344,7 +1342,6 @@ namespace grb { assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] == C_col_index[ j ] ); } - assert( old_nzc == nzc ); internal::setCurrentNonzeroes( C, nzc ); diff --git a/tests/unit/matrixSet.cpp b/tests/unit/matrixSet.cpp index 40ee0a69e..15eac3bd8 100644 --- a/tests/unit/matrixSet.cpp +++ b/tests/unit/matrixSet.cpp @@ -43,9 +43,9 @@ void grb_program( const size_t &n, grb::RC &rc ) { grb::Matrix< void > C( n, n ); grb::Matrix< void > D( n, n ); grb::Matrix< unsigned int > E( n, n ); - grb::Matrix< int > Mask( n, n ); - grb::Matrix< int > Output( n, n ); - grb::Matrix< int > Input( n, n ); + grb::Matrix< int > mask( n, n ); + grb::Matrix< int > output( n, n ); + grb::Matrix< int > input( n, n ); @@ -84,13 +84,13 @@ void grb_program( const size_t &n, grb::RC &rc ) { } } - rc = grb::buildMatrixUnique( Mask, I_mask, J_mask, mask_vals, 2 * n - 1, SEQUENTIAL ); + rc = grb::buildMatrixUnique( mask, I_mask, J_mask, mask_vals, 2 * n - 1, SEQUENTIAL ); if( rc != SUCCESS ) { std::cerr << "\t buildMatrixUnique of mask matrix FAILED\n"; return; } - rc = grb::buildMatrixUnique( Input, I_mask, J_mask, input_vals, 2 * n - 1, SEQUENTIAL ); + rc = grb::buildMatrixUnique( input, I_mask, J_mask, input_vals, 2 * n - 1, SEQUENTIAL ); if( rc != SUCCESS ) { std::cerr << "\t buildMatrixUnique of input matrix FAILED\n"; return; @@ -109,7 +109,7 @@ void grb_program( const size_t &n, grb::RC &rc ) { rc = grb::resize( E, 15 ); } if( rc == SUCCESS ) { - rc = grb::resize(Output, 2 * n - 1 ); + rc = grb::resize(output, 2 * n - 1 ); } if( rc != SUCCESS || grb::nnz( A ) != 15 ) { std::cerr << "\tinitialisation FAILED\n"; @@ -242,18 +242,18 @@ void grb_program( const size_t &n, grb::RC &rc ) { //check masked matrix set - rc = grb::set< descriptors::structural >( Output, Mask, Input ); + rc = grb::set< descriptors::structural >( output, mask, input ); if( rc != SUCCESS ) { std::cerr << "\t grb::set structural (matrix to matrix masked) FAILED\n"; return; } - if( grb::nnz( Output ) != 2 * n - 1 ) { - std::cerr << "\t unexpected number of output elements ( " << grb::nnz( Output ) << " ), expected " << 2 * n - 1 <<".\n"; + if( grb::nnz( output ) != 2 * n - 1 ) { + std::cerr << "\t unexpected number of output elements ( " << grb::nnz( output ) << " ), expected " << 2 * n - 1 <<".\n"; rc = FAILED; } - for( const auto & triplet : Output ) { + for( const auto & triplet : output ) { if( triplet.first.first != triplet.first.second && triplet.first.first != triplet.first.second - 1 ) { std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ), value " << triplet.second << ".\n"; rc = FAILED; @@ -270,18 +270,18 @@ void grb_program( const size_t &n, grb::RC &rc ) { - rc = grb::set( Output, Mask, Input ); + rc = grb::set( output, mask, input ); if( rc != SUCCESS ) { std::cerr << "\t grb::set (matrix to matrix masked) FAILED\n"; return; } - if( grb::nnz( Output ) != n ) { - std::cerr << "\t unexpected number of output elements ( " << grb::nnz( Output ) << " ), expected " << n <<".\n"; + if( grb::nnz( output ) != n ) { + std::cerr << "\t unexpected number of output elements ( " << grb::nnz( output ) << " ), expected " << n <<".\n"; rc = FAILED; } - for( const auto & triplet : Output ) { + for( const auto & triplet : output ) { if( triplet.first.first != triplet.first.second ) { std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ), value " << triplet.second << ".\n"; rc = FAILED; @@ -298,18 +298,18 @@ void grb_program( const size_t &n, grb::RC &rc ) { - rc = grb::set< descriptors::invert_mask >( Output, Mask, Input ); + rc = grb::set< descriptors::invert_mask >( output, mask, input ); if( rc != SUCCESS ) { std::cerr << "\t grb::set invert mask (matrix to matrix masked) FAILED\n"; return; } - if( grb::nnz( Output ) != n - 1 ) { - std::cerr << "\t unexpected number of output elements ( " << grb::nnz( Output ) << " ), expected " << n - 1 <<".\n"; + if( grb::nnz( output ) != n - 1 ) { + std::cerr << "\t unexpected number of output elements ( " << grb::nnz( output ) << " ), expected " << n - 1 <<".\n"; rc = FAILED; } - for( const auto & triplet : Output ) { + for( const auto & triplet : output ) { if( triplet.first.first != triplet.first.second - 1 ) { std::cerr << "\tunexpected entry at ( " << triplet.first.first << ", " << triplet.first.second << " ), value " << triplet.second << ".\n"; rc = FAILED; From 883479ef1134a33f515238f84eade2b8fb186114 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Mon, 26 Feb 2024 14:22:45 +0100 Subject: [PATCH 43/51] code cleanup --- include/graphblas/reference/blas3.hpp | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 3faee95b0..432ed53ca 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -1015,7 +1015,7 @@ namespace grb { for( size_t i = 0; i < nrows; ++i ) { if( internal::getCoordinates( u ).assigned( i ) ) { for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = mask_raw.row_index[ k ]; + const auto k_col = mask_raw.row_index[ k ]; if( internal::getCoordinates( v ).assigned( k_col ) && utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) @@ -1077,14 +1077,14 @@ namespace grb { for( size_t i = 0; i < nrows; ++i ) { if( internal::getCoordinates( u ).assigned( i ) ) { for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = mask_raw.row_index[ k ]; + const auto k_col = mask_raw.row_index[ k ]; if( internal::getCoordinates( v ).assigned( k_col ) && utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) ) { - (void) ++nzc; + nzc++; if( !crs_only ) { - (void) ++CCS_raw.col_start[ k_col + 1 ]; + CCS_raw.col_start[ k_col + 1 ]++; } } } @@ -1112,7 +1112,7 @@ namespace grb { for( size_t i = 0; i < nrows; ++i ) { if( internal::getCoordinates( u ).assigned( i ) ) { for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = mask_raw.row_index[ k ]; + const auto k_col = mask_raw.row_index[ k ]; if( internal::getCoordinates( v ).assigned( k_col ) && utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) @@ -1126,23 +1126,27 @@ namespace grb { CRS_raw.setValue( nzc, val ); // update CCS if( !crs_only ) { - const size_t CCS_index = A_col_index[ k_col ]++ + CCS_raw.col_start[ k_col ]; + const auto CCS_index = A_col_index[ k_col ] + CCS_raw.col_start[ k_col ]; + A_col_index[ k_col ]++; CCS_raw.row_index[ CCS_index ] = i; CCS_raw.setValue( CCS_index, val ); } // update count - (void) ++nzc; + nzc++; } } } CRS_raw.col_start[ i + 1 ] = nzc; } + +#ifndef NDEBUG if( !crs_only ) { for( size_t j = 0; j < ncols; ++j ) { assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] == A_col_index[ j ] ); } } +#endif internal::setCurrentNonzeroes( A, nzc ); From c36b4d9d4571339342e2042c9a779a83584cefcf Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Mon, 26 Feb 2024 15:33:13 +0100 Subject: [PATCH 44/51] handle void matrix mask --- include/graphblas/reference/blas3.hpp | 8 +- .../reference/compressed_storage.hpp | 16 ++ include/graphblas/utils.hpp | 33 ++++ tests/unit/maskedOuter.cpp | 162 +++++++++++++++++- 4 files changed, 208 insertions(+), 11 deletions(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 432ed53ca..236b41016 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -1001,7 +1001,7 @@ namespace grb { char * mask_arr = nullptr; char * mask_buf = nullptr; - OutputType * mask_valbuf = nullptr; + MaskType * mask_valbuf = nullptr; internal::getMatrixBuffers( mask_arr, mask_buf, mask_valbuf, 1, mask ); internal::Coordinates< reference > mask_coors; @@ -1018,7 +1018,7 @@ namespace grb { const auto k_col = mask_raw.row_index[ k ]; if( internal::getCoordinates( v ).assigned( k_col ) && - utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) + utils::interpretMatrixMask< descr, MaskType >( true, mask_raw.getValues(), k ) ) { nzc++; @@ -1080,7 +1080,7 @@ namespace grb { const auto k_col = mask_raw.row_index[ k ]; if( internal::getCoordinates( v ).assigned( k_col ) && - utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) + utils::interpretMatrixMask< descr, MaskType >( true, mask_raw.getValues(), k ) ) { nzc++; if( !crs_only ) { @@ -1115,7 +1115,7 @@ namespace grb { const auto k_col = mask_raw.row_index[ k ]; if( internal::getCoordinates( v ).assigned( k_col ) && - utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) + utils::interpretMatrixMask< descr, MaskType >( true, mask_raw.getValues(), k ) ) { OutputType val; grb::apply( val, diff --git a/include/graphblas/reference/compressed_storage.hpp b/include/graphblas/reference/compressed_storage.hpp index fa1aac97c..5d7939b87 100644 --- a/include/graphblas/reference/compressed_storage.hpp +++ b/include/graphblas/reference/compressed_storage.hpp @@ -425,6 +425,15 @@ namespace grb { return values; } + /** + * @returns The value array. + * + * \warning Does not check for NULL pointers. + */ + D * getValues() const noexcept { + return values; + } + /** * @returns The index array. * @@ -1105,6 +1114,13 @@ namespace grb { return nullptr; } + /** + * @returns A null pointer (since this is a pattern matrix). + */ + char * getValues() const noexcept { + return nullptr; + } + /** * @returns The index array. * diff --git a/include/graphblas/utils.hpp b/include/graphblas/utils.hpp index c5239afdc..02c0d3bd0 100644 --- a/include/graphblas/utils.hpp +++ b/include/graphblas/utils.hpp @@ -381,6 +381,39 @@ namespace grb { } } + /** Specialisation for void-valued matrice's masks */ + template< Descriptor descriptor, typename MatrixDataType, typename ValuesType > + static bool interpretMatrixMask( + const bool &assigned, + const ValuesType * const values, + const size_t k, + typename std::enable_if< + !std::is_void< MatrixDataType >::value + >::type * = nullptr + ) { + return interpretMask< descriptor, ValuesType >( assigned, values, k ); + } + + /** Specialisation for void-valued matrice's masks */ + template< Descriptor descriptor, typename MatrixDataType, typename ValuesType > + static bool interpretMatrixMask( + const bool &assigned, + const ValuesType * const, + const size_t, + typename std::enable_if< + std::is_void< MatrixDataType >::value + >::type * = nullptr + ) { + // set default mask to false + bool ret = assigned; + // check whether we should return the inverted value + if( descriptor & descriptors::invert_mask ) { + return !ret; + } else { + return ret; + } + } + } // namespace utils } // namespace grb diff --git a/tests/unit/maskedOuter.cpp b/tests/unit/maskedOuter.cpp index 7cfc8756c..b61ee4a0a 100644 --- a/tests/unit/maskedOuter.cpp +++ b/tests/unit/maskedOuter.cpp @@ -38,7 +38,7 @@ static const size_t J2[ 8 ] = { 0, 2, 1, 2, 0, 2, 0, 1 }; static const double mask_test1_expect[ 3 ] = { 4, 10, 18 }; static const double mask_test2_expect[ 3 ] = { 11, 20, 27 }; -void grbProgram( const void *, const size_t in_size, int &error ) { +void grb_program( const void *, const size_t in_size, int &error ) { error = 0; if( in_size != 0 ) { @@ -65,7 +65,7 @@ void grbProgram( const void *, const size_t in_size, int &error ) { grb::identities::zero, grb::identities::one > ring; - grb::RC rc = grb::buildMatrixUnique( Mask1, &( I1[ 0 ] ), &( J1[ 0 ] ), m1, 3, SEQUENTIAL ); + grb::RC rc = grb::buildMatrixUnique( Mask1, &( I1[ 0 ] ), &( J1[ 0 ] ), m1, 4, SEQUENTIAL ); if( rc != SUCCESS ) { std::cerr << "\t first mask buildMatrix FAILED\n"; @@ -301,6 +301,144 @@ void grb_program_custom_size( const size_t &n, int &error ) { } +static const double mask_void_test1_expect[ 3 ] = { 10, 10, 18 }; +static const double mask_void_test2_expect[ 3 ] = { 15, 20, 45 }; + +void grb_program_void( const void *, const size_t in_size, int &error ) { + error = 0; + + if( in_size != 0 ) { + (void)fprintf( stderr, "Unit tests called with unexpected input\n" ); + error = 1; + return; + } + + // allocate + grb::Vector< double > u = { 1, 2, 3 }; + grb::Vector< double > v = { 4, 5, 6 }; + grb::Matrix< double > Result1( 3, 3 ); + grb::Matrix< double > Result2( 3, 3 ); + grb::Matrix< void > Mask1( 3, 3 ); + grb::Matrix< void > Mask2( 3, 3 ); + grb::Vector< double > test1 = { 1, 1, 1 }; + grb::Vector< double > out1( 3 ); + grb::Vector< double > test2 = { 1, 1, 1 }; + grb::Vector< double > out2( 3 ); + + // semiring + grb::Semiring< + grb::operators::add< double >, grb::operators::mul< double >, + grb::identities::zero, grb::identities::one + > ring; + + grb::RC rc = grb::buildMatrixUnique( Mask1, &( I1[ 0 ] ), &( J1[ 0 ] ), 4, SEQUENTIAL ); + + if( rc != SUCCESS ) { + std::cerr << "\t first mask buildMatrix FAILED\n"; + error = 5; + } + + if( !error ) { + rc = grb::buildMatrixUnique( Mask2, &( I2[ 0 ] ), &( J2[ 0 ] ), 8, SEQUENTIAL ); + if( rc != SUCCESS ) { + std::cerr << "\t second mask buildMatrix FAILED\n"; + error = 10; + } + } + + + + if( !error ) { + rc = grb::outer( Result1, Mask1, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer( Result1, Mask1, u, v, ring.getMultiplicativeOperator() ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::outer: " + << toString( rc ) << ".\n"; + error = 15; + } + } + + + if( !error && grb::nnz( Result1 ) != 4 ) { + std::cerr << "\t Unexpected number of nonzeroes in matrix: " + << grb::nnz( Result1 ) << ", expected 4.\n"; + error = 20; + } + + if( !error ) { + rc = grb::outer< descriptors::force_row_major | descriptors::invert_mask >( Result2, Mask2, u, v, ring.getMultiplicativeOperator(), RESIZE ); + rc = rc ? rc : grb::outer< descriptors::force_row_major | descriptors::invert_mask >( Result2, Mask2, u, v, ring.getMultiplicativeOperator() ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from grb::outer: " + << toString( rc ) << ".\n"; + error = 25; + } + } + + + if( !error && grb::nnz( Result2 ) != 8 ) { + std::cerr << "\t Unexpected number of nonzeroes in matrix: " + << grb::nnz( Result2 ) << ", expected 8.\n"; + error = 30; + } + + if( !error ) { + rc = grb::vxm( out1, test1, Result1, ring ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from premultiplying Result1 by a vector (vxm): " + << toString( rc ) << ".\n"; + error = 35; + } + } + + if( !error ) { + if( grb::nnz( out1 ) != 3 ) { + std::cerr << "\t Unexpected number of nonzeroes (premultiply): " + << grb::nnz( out1 ) << ", expected 3\n"; + error = 40; + } + for( const auto &pair : out1 ) { + size_t i = pair.first; + if( pair.second != mask_test1_expect[ i ] ) { + std::cerr << "Premultiplying Result1 by a vector of all ones, " + << "unexpected value " << pair.second << " " + << "at coordinate " << i << ", expected " + << mask_test1_expect[ i ] << ".\n"; + error = 45; + break; + } + } + } + + if( !error ) { + rc = grb::vxm< grb::descriptors::transpose_matrix >( out2, test2, Result2, ring ); + if( rc != grb::SUCCESS ) { + std::cerr << "Unexpected return code from postmultiplying Result2 by a vector (vxm): " + << toString( rc ) << ".\n"; + error = 60; + } + } + + if( !error ) { + if( grb::nnz( out2 ) != 3 ) { + std::cerr << "\t Unexpected number of nonzeroes (postmultiply): " + << grb::nnz( out1 ) << ", expected 3\n"; + error = 65; + } + for( const auto &pair : out2 ) { + size_t i = pair.first; + if( pair.second != mask_test2_expect[ i ] ) { + std::cerr << "Postmultiplying Result2 by a vector of all ones, " + << "unexpected value " << pair.second << " " + << "at coordinate " << i << ", " + << "expected " << mask_test2_expect[ i ] << ".\n"; + error = 70; + break; + } + } + } +} + int main( int argc, char ** argv ) { (void)argc; std::cout << "Functional test executable: " << argv[ 0 ] << "\n"; @@ -338,7 +476,7 @@ int main( int argc, char ** argv ) { int error = 0; grb::Launcher< AUTOMATIC > launcher; - if( launcher.exec( &grbProgram, nullptr, 0, error ) != SUCCESS ) { + if( launcher.exec( &grb_program, nullptr, 0, error ) != SUCCESS ) { std::cerr << "Test 1 failed to launch\n"; error = 255; } @@ -352,14 +490,24 @@ int main( int argc, char ** argv ) { std::cerr << "Launching test 2 FAILED\n"; error = 255; } - if( error == 0 ) { - std::cout << "Test OK\n" << std::endl; - } else { + if( error != 0 ) { std::cerr << std::flush; std::cout << "Test 2 FAILED\n" << std::endl; + return error; + } + + if( launcher.exec( &grb_program_void, nullptr, 0, error ) != SUCCESS ) { + std::cerr << "Test 3 failed to launch\n"; + error = 255; + } + if( error != 0 ) { + std::cerr << std::flush; + std::cout << "Test 3 FAILED\n" << std::endl; + return error; } // done - return error; + std::cout << "Test OK\n" << std::endl; + return 0; } From e5223c1586a8c001e3ab11e172e43910788770be Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 27 Feb 2024 10:34:43 +0100 Subject: [PATCH 45/51] handle invert masks properly --- include/graphblas/utils.hpp | 14 ++++++++++---- tests/unit/maskedOuter.cpp | 10 +++++----- 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/include/graphblas/utils.hpp b/include/graphblas/utils.hpp index 02c0d3bd0..e62a4fb8b 100644 --- a/include/graphblas/utils.hpp +++ b/include/graphblas/utils.hpp @@ -323,7 +323,10 @@ namespace grb { ret = assigned; } else { // if based on value, if there is a value, cast it to bool - if( assigned ) { + if( !assigned ) { + return false; + } + else { ret = static_cast< bool >( val[ offset ] ); } // otherwise there is no value and false is assumed @@ -350,7 +353,10 @@ namespace grb { ret = assigned; } else { // if based on value, if there is a value, cast it to bool - if( assigned ) { + if( !assigned ) { + return false; + } + else { ret = static_cast< bool >( real( val [ offset ] ) ) || static_cast< bool >( imag( val [ offset ] ) ); } @@ -374,7 +380,7 @@ namespace grb { // set default mask to false bool ret = assigned; // check whether we should return the inverted value - if( descriptor & descriptors::invert_mask ) { + if( ( descriptor & descriptors::structural_complement ) == descriptors::structural_complement ) { return !ret; } else { return ret; @@ -407,7 +413,7 @@ namespace grb { // set default mask to false bool ret = assigned; // check whether we should return the inverted value - if( descriptor & descriptors::invert_mask ) { + if( ( descriptor & descriptors::structural_complement ) == descriptors::structural_complement ) { return !ret; } else { return ret; diff --git a/tests/unit/maskedOuter.cpp b/tests/unit/maskedOuter.cpp index b61ee4a0a..4256a9a3d 100644 --- a/tests/unit/maskedOuter.cpp +++ b/tests/unit/maskedOuter.cpp @@ -301,7 +301,7 @@ void grb_program_custom_size( const size_t &n, int &error ) { } -static const double mask_void_test1_expect[ 3 ] = { 10, 10, 18 }; +static const double mask_void_test1_expect[ 3 ] = { 4, 10, 24 }; static const double mask_void_test2_expect[ 3 ] = { 15, 20, 45 }; void grb_program_void( const void *, const size_t in_size, int &error ) { @@ -399,11 +399,11 @@ void grb_program_void( const void *, const size_t in_size, int &error ) { } for( const auto &pair : out1 ) { size_t i = pair.first; - if( pair.second != mask_test1_expect[ i ] ) { + if( pair.second != mask_void_test1_expect[ i ] ) { std::cerr << "Premultiplying Result1 by a vector of all ones, " << "unexpected value " << pair.second << " " << "at coordinate " << i << ", expected " - << mask_test1_expect[ i ] << ".\n"; + << mask_void_test1_expect[ i ] << ".\n"; error = 45; break; } @@ -427,11 +427,11 @@ void grb_program_void( const void *, const size_t in_size, int &error ) { } for( const auto &pair : out2 ) { size_t i = pair.first; - if( pair.second != mask_test2_expect[ i ] ) { + if( pair.second != mask_void_test2_expect[ i ] ) { std::cerr << "Postmultiplying Result2 by a vector of all ones, " << "unexpected value " << pair.second << " " << "at coordinate " << i << ", " - << "expected " << mask_test2_expect[ i ] << ".\n"; + << "expected " << mask_void_test2_expect[ i ] << ".\n"; error = 70; break; } From 46f58f37ab82457d62b12179fdf0aa5b0cf7d6b9 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 27 Feb 2024 10:44:51 +0100 Subject: [PATCH 46/51] code cleanup --- .../reference/compressed_storage.hpp | 17 ++++++++++ include/graphblas/reference/io.hpp | 28 ++++++++-------- include/graphblas/utils.hpp | 32 +++++++++++++++++++ 3 files changed, 64 insertions(+), 13 deletions(-) diff --git a/include/graphblas/reference/compressed_storage.hpp b/include/graphblas/reference/compressed_storage.hpp index fa1aac97c..8e5082b62 100644 --- a/include/graphblas/reference/compressed_storage.hpp +++ b/include/graphblas/reference/compressed_storage.hpp @@ -425,6 +425,16 @@ namespace grb { return values; } + /** + * @returns The value array. + * + * \warning Does not check for NULL pointers. + */ + D * getValues() const noexcept { + return values; + } + + /** * @returns The index array. * @@ -1105,6 +1115,13 @@ namespace grb { return nullptr; } + /** + * @returns A null pointer (since this is a pattern matrix). + */ + char * getValues() const noexcept { + return nullptr; + } + /** * @returns The index array. * diff --git a/include/graphblas/reference/io.hpp b/include/graphblas/reference/io.hpp index c1c524725..3dd529a17 100644 --- a/include/graphblas/reference/io.hpp +++ b/include/graphblas/reference/io.hpp @@ -1225,7 +1225,7 @@ namespace grb { char * mask_arr = nullptr; char * mask_buf = nullptr; - OutputType * mask_valbuf = nullptr; + MaskType * mask_valbuf = nullptr; internal::getMatrixBuffers( mask_arr, mask_buf, mask_valbuf, 1, M ); internal::Coordinates< reference > mask_coors; @@ -1234,13 +1234,13 @@ namespace grb { for( size_t i = 0; i < nrows; ++i ) { mask_coors.clear(); for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = mask_raw.row_index[ k ]; - if( utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) ) { + const auto k_col = mask_raw.row_index[ k ]; + if( utils::interpretMatrixMask< descr, MaskType >( true, mask_raw.getValues(), k ) ) { mask_coors.assign( k_col ); } } for( auto k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = A_raw.row_index[ k ]; + const auto k_col = A_raw.row_index[ k ]; if( mask_coors.assigned( k_col ) ) { nzc++; } @@ -1275,7 +1275,7 @@ namespace grb { CRS_raw.col_start[ 0 ] = 0; #ifdef _H_GRB_REFERENCE_OMP_BLAS3 - #pragma omp parallel for simd + #pragma omp parallel for simd #endif for( size_t j = 0; j <= ncols; ++j ) { CCS_raw.col_start[ j ] = 0; @@ -1286,13 +1286,13 @@ namespace grb { for( size_t i = 0; i < nrows; ++i ) { mask_coors.clear(); for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = mask_raw.row_index[ k ]; - if( utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) ) { + const auto k_col = mask_raw.row_index[ k ]; + if( utils::interpretMask< descr, MaskType >( true, mask_raw.getValues(), k ) ) { mask_coors.assign( k_col ); } } for( auto k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = A_raw.row_index[ k ]; + const auto k_col = A_raw.row_index[ k ]; if( mask_coors.assigned( k_col ) ) { nzc++; CCS_raw.col_start[ k_col + 1 ]++; @@ -1319,29 +1319,31 @@ namespace grb { for( size_t i = 0; i < nrows; ++i ) { mask_coors.clear(); for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = mask_raw.row_index[ k ]; - if( utils::interpretMask< descr, MaskType >( true, mask_raw.values, k ) ) { + const auto k_col = mask_raw.row_index[ k ]; + if( utils::interpretMatrixMask< descr, MaskType >( true, mask_raw.getValues(), k ) ) { mask_coors.assign( k_col ); } } for( auto k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = A_raw.row_index[ k ]; + const auto k_col = A_raw.row_index[ k ]; if( mask_coors.assigned( k_col ) ) { OutputType val = A_raw.getValue( k, *( (OutputType*) nullptr ) ); CRS_raw.row_index[ nzc ] = k_col; CRS_raw.setValue( nzc, val ); - const size_t CCS_index = C_col_index[ k_col ]++ + CCS_raw.col_start[ k_col ]; + const size_t CCS_index = C_col_index[ k_col ] + CCS_raw.col_start[ k_col ]; + C_col_index[ k_col ]++; CCS_raw.row_index[ CCS_index ] = i; CCS_raw.setValue( CCS_index, val ); nzc++; } } } - +#ifndef NDEBUG for( size_t j = 0; j < ncols; ++j ) { assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] == C_col_index[ j ] ); } +#endif internal::setCurrentNonzeroes( C, nzc ); diff --git a/include/graphblas/utils.hpp b/include/graphblas/utils.hpp index c5239afdc..e2590d882 100644 --- a/include/graphblas/utils.hpp +++ b/include/graphblas/utils.hpp @@ -380,6 +380,38 @@ namespace grb { return ret; } } + /** Specialisation for void-valued matrice's masks */ + template< Descriptor descriptor, typename MatrixDataType, typename ValuesType > + static bool interpretMatrixMask( + const bool &assigned, + const ValuesType * const values, + const size_t k, + typename std::enable_if< + !std::is_void< MatrixDataType >::value + >::type * = nullptr + ) { + return interpretMask< descriptor, ValuesType >( assigned, values, k ); + } + + /** Specialisation for void-valued matrice's masks */ + template< Descriptor descriptor, typename MatrixDataType, typename ValuesType > + static bool interpretMatrixMask( + const bool &assigned, + const ValuesType * const, + const size_t, + typename std::enable_if< + std::is_void< MatrixDataType >::value + >::type * = nullptr + ) { + // set default mask to false + bool ret = assigned; + // check whether we should return the inverted value + if( descriptor & descriptors::invert_mask ) { + return !ret; + } else { + return ret; + } + } } // namespace utils From cf2e317cc26ccbfc87c6407750767da7598c65f5 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 27 Feb 2024 10:58:39 +0100 Subject: [PATCH 47/51] omp reductions --- include/graphblas/reference/io.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/include/graphblas/reference/io.hpp b/include/graphblas/reference/io.hpp index 3dd529a17..127394e95 100644 --- a/include/graphblas/reference/io.hpp +++ b/include/graphblas/reference/io.hpp @@ -1239,6 +1239,9 @@ namespace grb { mask_coors.assign( k_col ); } } +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for reduction(+:nzc) +#endif for( auto k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { const auto k_col = A_raw.row_index[ k ]; if( mask_coors.assigned( k_col ) ) { From 5c639a8549624377c9e04917f9396e33a96c66a7 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 27 Feb 2024 16:05:24 +0100 Subject: [PATCH 48/51] draft of the submask selection functionality --- include/graphblas/reference/blas3.hpp | 69 +++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 868547183..d68eee6c3 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -799,6 +799,75 @@ namespace grb { return internal::matrix_zip_generic< descr, true >( A, x, y, x, phase ); } + /** + * Selecting a submatrix of matrix \a B based on the given vector of \a rows and \a cols + * Depends on the masked outer product and masked set. Currently, only the structural version is supported. + */ + + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename InputType, typename MaskType, typename OutputType, + typename Coords, + typename RIT, typename CIT, typename NIT + > + RC selectSubmatrix( + Matrix< OutputType, reference, RIT, CIT, NIT > &B, + const Matrix< InputType, reference, RIT, CIT, NIT > &A, + const Vector< MaskType, reference, Coords > &rows, + const Vector< MaskType, reference, Coords > &cols, + const typename std::enable_if< + grb::is_operator< Operator >::value && + !grb::is_object< InputType >::value && + !grb::is_object< MaskType >::value && + !grb::is_object< OutputType >::value, + void >::type * const = nullptr + ) { + //static asserts + static_assert( + ! ( descr & descriptors::transpose_matrix ), + "grb::selectSubmatrix can not be called with descriptors::transpose_matrix" + ); + static_assert( + ( descr & descriptors::structural ), + "Only the structural version is supported for grb::selectSubmatrix" + ); +#ifdef _DEBUG + std::cout << "In grb::selectSubmatrix (reference)\n"; +#endif + + const size_t nrows = size( rows ); + const size_t ncols = size( cols ); + if( nrows != grb::nrows( A ) || nrows != grb::nrows( B ) ) { + return MISMATCH; + } + + if( ncols != grb::ncols( A ) || ncols != grb::ncols( B ) ) { + return MISMATCH; + } + + //mask contains only those values that need to be selected from A + Matrix< MaskType, reference, RIT, CIT, NIT > mask( nrows, ncols ); + + ret = grb::maskedOuter( mask, A, rows, cols, grb::operators::zip< MaskType, MaskType, reference >(), Phase::RESIZE ); + if( ret != SUCCESS ) { + return ret; + } + ret = grb::maskedOuter( mask, A, rows, cols, grb::operators::zip< MaskType, MaskType, reference >() ); + if( ret != SUCCESS ) { + return ret; + } + + ret = grb::set( B, mask, A, Phase::RESIZE ); + if( ret != SUCCESS ) { + return ret; + } + ret = grb::set( B, mask, A ); + return ret; + } + + + /** * Outer product of two vectors. Assuming vectors \a u and \a v are oriented * column-wise, the result matrix \a A will contain \f$ uv^T \f$. This is an From 8ecab8ad3f98f54affd437a39ffb1668e00eb383 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 27 Feb 2024 16:08:27 +0100 Subject: [PATCH 49/51] rename --- include/graphblas/reference/blas3.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index d68eee6c3..52367c56d 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -849,11 +849,11 @@ namespace grb { //mask contains only those values that need to be selected from A Matrix< MaskType, reference, RIT, CIT, NIT > mask( nrows, ncols ); - ret = grb::maskedOuter( mask, A, rows, cols, grb::operators::zip< MaskType, MaskType, reference >(), Phase::RESIZE ); + ret = grb::outer( mask, A, rows, cols, grb::operators::zip< MaskType, MaskType, reference >(), Phase::RESIZE ); if( ret != SUCCESS ) { return ret; } - ret = grb::maskedOuter( mask, A, rows, cols, grb::operators::zip< MaskType, MaskType, reference >() ); + ret = grb::outer( mask, A, rows, cols, grb::operators::zip< MaskType, MaskType, reference >() ); if( ret != SUCCESS ) { return ret; } From eff8f6ccfe693d93079aad283f7736195820a6e4 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 27 Feb 2024 16:11:51 +0100 Subject: [PATCH 50/51] rename --- include/graphblas/reference/blas3.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 52367c56d..1dcc15b44 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -849,20 +849,20 @@ namespace grb { //mask contains only those values that need to be selected from A Matrix< MaskType, reference, RIT, CIT, NIT > mask( nrows, ncols ); - ret = grb::outer( mask, A, rows, cols, grb::operators::zip< MaskType, MaskType, reference >(), Phase::RESIZE ); + ret = outer( mask, A, rows, cols, grb::operators::zip< MaskType, MaskType, reference >(), Phase::RESIZE ); if( ret != SUCCESS ) { return ret; } - ret = grb::outer( mask, A, rows, cols, grb::operators::zip< MaskType, MaskType, reference >() ); + ret = outer( mask, A, rows, cols, grb::operators::zip< MaskType, MaskType, reference >() ); if( ret != SUCCESS ) { return ret; } - ret = grb::set( B, mask, A, Phase::RESIZE ); + ret = set( B, mask, A, Phase::RESIZE ); if( ret != SUCCESS ) { return ret; } - ret = grb::set( B, mask, A ); + ret = set( B, mask, A ); return ret; } From 0de1411b749323f72b33577cdad221b665982a87 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Tue, 27 Feb 2024 16:19:50 +0100 Subject: [PATCH 51/51] merge of required features --- include/graphblas/reference/blas3.hpp | 2 +- include/graphblas/utils.hpp | 32 --------------------------- 2 files changed, 1 insertion(+), 33 deletions(-) diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 68a4fc6cb..cb044fc4f 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -849,7 +849,7 @@ namespace grb { //mask contains only those values that need to be selected from A Matrix< MaskType, reference, RIT, CIT, NIT > mask( nrows, ncols ); - ret = outer( mask, A, rows, cols, grb::operators::zip< MaskType, MaskType, reference >(), Phase::RESIZE ); + RC ret = outer( mask, A, rows, cols, grb::operators::zip< MaskType, MaskType, reference >(), Phase::RESIZE ); if( ret != SUCCESS ) { return ret; } diff --git a/include/graphblas/utils.hpp b/include/graphblas/utils.hpp index 7052e6df8..e62a4fb8b 100644 --- a/include/graphblas/utils.hpp +++ b/include/graphblas/utils.hpp @@ -419,38 +419,6 @@ namespace grb { return ret; } } - /** Specialisation for void-valued matrice's masks */ - template< Descriptor descriptor, typename MatrixDataType, typename ValuesType > - static bool interpretMatrixMask( - const bool &assigned, - const ValuesType * const values, - const size_t k, - typename std::enable_if< - !std::is_void< MatrixDataType >::value - >::type * = nullptr - ) { - return interpretMask< descriptor, ValuesType >( assigned, values, k ); - } - - /** Specialisation for void-valued matrice's masks */ - template< Descriptor descriptor, typename MatrixDataType, typename ValuesType > - static bool interpretMatrixMask( - const bool &assigned, - const ValuesType * const, - const size_t, - typename std::enable_if< - std::is_void< MatrixDataType >::value - >::type * = nullptr - ) { - // set default mask to false - bool ret = assigned; - // check whether we should return the inverted value - if( descriptor & descriptors::invert_mask ) { - return !ret; - } else { - return ret; - } - } } // namespace utils