From eff687b7bbb8a68239730401d3e9dccd14ce6d61 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Fri, 26 Jan 2024 13:47:44 +0100 Subject: [PATCH 01/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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/41] 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 da9de374ec2913818cb5e6ac77de661f733ae178 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Sat, 24 Feb 2024 09:53:43 +0100 Subject: [PATCH 33/41] 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 34/41] 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 35/41] 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 36/41] 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 883479ef1134a33f515238f84eade2b8fb186114 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Mon, 26 Feb 2024 14:22:45 +0100 Subject: [PATCH 37/41] 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 38/41] 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 39/41] 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 8e71a735a3510fcafab9aa90ab62164242c622c4 Mon Sep 17 00:00:00 2001 From: Aleksa Milisavljevic Date: Wed, 28 Feb 2024 11:56:17 +0100 Subject: [PATCH 40/41] add name to notice --- NOTICE | 2 ++ 1 file changed, 2 insertions(+) diff --git a/NOTICE b/NOTICE index 54090d018..6646d92b1 100644 --- a/NOTICE +++ b/NOTICE @@ -34,6 +34,8 @@ to Huawei Technologies Co., Ltd. or one of its subsidiaries: - Denis Jelovina, Huawei Technologies Switzerland AG; 2022-current. - Benjamin Lozes, Huawei Technologies Switzerland AG; 2023. + + - Aleksa Milisavljevic, Huawei Technologies Switzerland AG; 2024. The experimental banshee backend has been developed in collaboration with Prof. Luca Benini at ETH Zuerich and his group. In particular this backend From ff8466274ef274a291fcbe1379569717afac7c05 Mon Sep 17 00:00:00 2001 From: Benjamin Lozes <62449444+byjtew@users.noreply.github.com> Date: Wed, 28 Feb 2024 17:28:00 +0100 Subject: [PATCH 41/41] Update NOTICE --- NOTICE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NOTICE b/NOTICE index 6646d92b1..342839658 100644 --- a/NOTICE +++ b/NOTICE @@ -33,7 +33,7 @@ to Huawei Technologies Co., Ltd. or one of its subsidiaries: - Denis Jelovina, Huawei Technologies Switzerland AG; 2022-current. - - Benjamin Lozes, Huawei Technologies Switzerland AG; 2023. + - Benjamin Lozes, Huawei Technologies Switzerland AG; 2023-2024. - Aleksa Milisavljevic, Huawei Technologies Switzerland AG; 2024.