diff --git a/include/graphblas/bsp1d/blas3.hpp b/include/graphblas/bsp1d/blas3.hpp index 386beb164..25f5902ab 100644 --- a/include/graphblas/bsp1d/blas3.hpp +++ b/include/graphblas/bsp1d/blas3.hpp @@ -131,20 +131,21 @@ namespace grb { const Matrix< InputType1, BSP1D, RIT2, CIT2, NIT2 > &A, const Matrix< InputType2, BSP1D, RIT3, CIT3, NIT3 > &B, const MulMonoid &mul, - const Phase phase = EXECUTE, + const Phase &phase = EXECUTE, const typename std::enable_if< - !grb::is_object< OutputType >::value && - !grb::is_object< InputType1 >::value && - !grb::is_object< InputType2 >::value && - grb::is_monoid< MulMonoid >::value, - void >::type * const = nullptr + !grb::is_object< OutputType >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + grb::is_monoid< MulMonoid >::value, void + >::type * const = nullptr ) { assert( phase != TRY ); RC local_rc = SUCCESS; if( phase == RESIZE ) { RC ret = eWiseApply< descr >( internal::getLocal( C ), - internal::getLocal( A ), internal::getLocal( B ), + internal::getLocal( A ), + internal::getLocal( B ), mul, RESIZE ); @@ -157,7 +158,8 @@ namespace grb { assert( phase == EXECUTE ); local_rc = eWiseApply< descr >( internal::getLocal( C ), - internal::getLocal( A ), internal::getLocal( B ), + internal::getLocal( A ), + internal::getLocal( B ), mul, EXECUTE ); @@ -179,31 +181,40 @@ namespace grb { const Matrix< InputType1, BSP1D, RIT2, CIT2, NIT2 > &A, const Matrix< InputType2, BSP1D, RIT3, CIT3, NIT3 > &B, const Operator &op, - const Phase phase = EXECUTE, + const Phase &phase = EXECUTE, const typename std::enable_if< - !grb::is_object< OutputType >::value && - !grb::is_object< InputType1 >::value && - !grb::is_object< InputType2 >::value && - grb::is_operator< Operator >::value, - void >::type * const = nullptr + !grb::is_object< OutputType >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + grb::is_operator< Operator >::value, void + >::type * const = nullptr ) { assert( phase != TRY ); - RC ret = eWiseApply< descr >( - internal::getLocal( C ), - internal::getLocal( A ), internal::getLocal( B ), - op, - phase - ); + RC local_rc = SUCCESS; if( phase == RESIZE ) { + RC ret = eWiseApply< descr >( + internal::getLocal( C ), + internal::getLocal( A ), + internal::getLocal( B ), + op, + RESIZE + ); if( collectives<>::allreduce( ret, operators::any_or< RC >() ) != SUCCESS ) { return PANIC; } else { - return SUCCESS; + return ret; } + } else { + assert( phase == EXECUTE ); + local_rc = eWiseApply< descr >( + internal::getLocal( C ), + internal::getLocal( A ), + internal::getLocal( B ), + op, + EXECUTE + ); } - assert( phase == EXECUTE ); - return internal::checkGlobalErrorStateOrClear( C, ret ); - } + return internal::checkGlobalErrorStateOrClear( C, local_rc ); } // namespace grb diff --git a/include/graphblas/nonblocking/blas3.hpp b/include/graphblas/nonblocking/blas3.hpp index c52c9aaff..50c640ae0 100644 --- a/include/graphblas/nonblocking/blas3.hpp +++ b/include/graphblas/nonblocking/blas3.hpp @@ -416,56 +416,6 @@ namespace grb { ); } - namespace internal { - - template< - bool allow_void, - Descriptor descr, - class MulMonoid, class Operator, - typename OutputType, typename InputType1, typename InputType2, - typename RIT1, typename CIT1, typename NIT1, - typename RIT2, typename CIT2, typename NIT2, - typename RIT3, typename CIT3, typename NIT3 - > - RC eWiseApply_matrix_generic( - Matrix< OutputType, nonblocking, RIT1, CIT1, NIT1 > &C, - const Matrix< InputType1, nonblocking, RIT2, CIT2, NIT2 > &A, - const Matrix< InputType2, nonblocking, RIT3, CIT3, NIT3 > &B, - const Operator &oper, - const MulMonoid &mulMonoid, - const Phase &phase, - const typename std::enable_if< - !grb::is_object< OutputType >::value && - !grb::is_object< InputType1 >::value && - !grb::is_object< InputType2 >::value && - grb::is_operator< Operator >::value, - void >::type * const = nullptr - ) { - if( internal::NONBLOCKING::warn_if_not_native && - config::PIPELINE::warn_if_not_native - ) { - std::cerr << "Warning: eWiseApply (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 - le.execution(); - - // second, delegate to the reference backend - return eWiseApply_matrix_generic< - allow_void, descr, - MulMonoid, Operator - >( - getRefMatrix( C ), getRefMatrix( A ), getRefMatrix( B ), - oper, mulMonoid, phase - ); - } - - } // namespace internal - template< Descriptor descr = descriptors::no_operation, class MulMonoid, @@ -507,11 +457,26 @@ namespace grb { ); #ifdef _DEBUG - std::cout << "In grb::eWiseApply_matrix_generic (nonblocking, monoid)\n"; + std::cout << "In grb::eWiseApply (nonblocking, monoid)\n"; #endif + if( internal::NONBLOCKING::warn_if_not_native && config::PIPELINE::warn_if_not_native ) { + std::cerr << "Warning: eWiseApply (nonblocking) currently delegates to a " + << "blocking implementation.\n" + << " Further similar such warnings will be suppressed.\n"; + internal::NONBLOCKING::warn_if_not_native = false; + } - return internal::eWiseApply_matrix_generic< true, descr >( - C, A, B, mulmono.getOperator(), mulmono, phase + // nonblocking execution is not supported + // first, execute any computation that is not completed + internal::le.execution(); + + // second, delegate to the reference backend + return eWiseApply< descr >( + internal::getRefMatrix( C ), + internal::getRefMatrix( A ), + internal::getRefMatrix( B ), + mulmono, + phase ); } @@ -535,6 +500,17 @@ namespace grb { grb::is_operator< Operator >::value, void >::type * const = nullptr ) { + if( internal::NONBLOCKING::warn_if_not_native && config::PIPELINE::warn_if_not_native ) { + std::cerr << "Warning: eWiseApply (nonblocking) currently delegates to a " + << "blocking implementation.\n" + << " Further similar such warnings will be suppressed.\n"; + internal::NONBLOCKING::warn_if_not_native = false; + } + +#ifdef _DEBUG + std::cout << "In grb::eWiseApply (nonblocking, op)\n"; +#endif + // static checks NO_CAST_ASSERT( ( !( descr & descriptors::no_casting ) || std::is_same< typename Operator::D1, InputType1 >::value ), @@ -554,20 +530,18 @@ namespace grb { "called with an output matrix C that does not match the output domain " "of the given multiplication operator" ); - static_assert( ( !( - std::is_same< InputType1, void >::value || - std::is_same< InputType2, void >::value ) - ), "grb::eWiseApply (nonblocking, matrix <- matrix x matrix, operator): " - "the operator version of eWiseApply cannot be used if either of the " - "input matrices is a pattern matrix (of type void)" - ); - typename grb::Monoid< - grb::operators::mul< double >, - grb::identities::one - > dummyMonoid; - return internal::eWiseApply_matrix_generic< false, descr >( - C, A, B, mulOp, dummyMonoid, phase + // nonblocking execution is not supported + // first, execute any computation that is not completed + internal::le.execution(); + + // second, delegate to the reference backend + return eWiseApply< descr >( + internal::getRefMatrix( C ), + internal::getRefMatrix( A ), + internal::getRefMatrix( B ), + mulOp, + phase ); } diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 868547183..563241ca3 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -57,6 +57,8 @@ "********************************************************************" \ "******************************\n" ); + + namespace grb { namespace internal { @@ -137,7 +139,7 @@ namespace grb { const auto &B_raw = !trans_right ? internal::getCRS( B ) : internal::getCCS( B ); - auto &C_raw = internal::getCRS( C ); + auto &CRS_raw = internal::getCRS( C ); auto &CCS_raw = internal::getCCS( C ); char * arr = nullptr; @@ -175,7 +177,7 @@ namespace grb { if( crs_only && phase == RESIZE ) { // we are using an auxialiary CRS that we cannot resize ourselves // instead, we update the offset array only - C_raw.col_start[ 0 ] = 0; + CRS_raw.col_start[ 0 ] = 0; } // if crs_only, then the below implements its resize phase // if not crs_only, then the below is both crucial for the resize phase, @@ -202,7 +204,7 @@ namespace grb { if( crs_only && phase == RESIZE ) { // we are using an auxialiary CRS that we cannot resize ourselves // instead, we update the offset array only - C_raw.col_start[ i + 1 ] = nzc; + CRS_raw.col_start[ i + 1 ] = nzc; } } } @@ -259,7 +261,7 @@ namespace grb { // use previously computed CCS offset array to update CCS during the // computational phase nzc = 0; - C_raw.col_start[ 0 ] = 0; + CRS_raw.col_start[ 0 ] = 0; for( size_t i = 0; i < m; ++i ) { coors.clear(); for( auto k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { @@ -302,8 +304,8 @@ namespace grb { assert( nzc < old_nzc ); const size_t j = coors.index( k ); // update CRS - C_raw.row_index[ nzc ] = j; - C_raw.setValue( nzc, valbuf[ j ] ); + CRS_raw.row_index[ nzc ] = j; + CRS_raw.setValue( nzc, valbuf[ j ] ); // update CCS if( !crs_only ) { const size_t CCS_index = C_col_index[ j ]++ + CCS_raw.col_start[ j ]; @@ -313,7 +315,7 @@ namespace grb { // update count (void) ++nzc; } - C_raw.col_start[ i + 1 ] = nzc; + CRS_raw.col_start[ i + 1 ] = nzc; } #ifndef NDEBUG @@ -918,135 +920,99 @@ namespace grb { namespace internal { - /** - * \internal general elementwise matrix application that all eWiseApply - * variants refer to. - * @param[in] oper The operator corresponding to \a mulMonoid if - * \a allow_void is true; otherwise, an arbitrary operator - * under which to perform the eWiseApply. - * @param[in] mulMonoid The monoid under which to perform the eWiseApply if - * \a allow_void is true; otherwise, will be ignored. - * \endinternal - */ - template< - bool allow_void, Descriptor descr, - class MulMonoid, class Operator, + class Operator, typename OutputType, typename InputType1, typename InputType2, typename RIT1, typename CIT1, typename NIT1, typename RIT2, typename CIT2, typename NIT2, typename RIT3, typename CIT3, typename NIT3 > - RC eWiseApply_matrix_generic( + RC eWiseApply_matrix_generic_intersection( Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C, const Matrix< InputType1, reference, RIT2, CIT2, NIT2 > &A, const Matrix< InputType2, reference, RIT3, CIT3, NIT3 > &B, const Operator &oper, - const MulMonoid &mulMonoid, const Phase &phase, const typename std::enable_if< - !grb::is_object< OutputType >::value && - !grb::is_object< InputType1 >::value && - !grb::is_object< InputType2 >::value && - grb::is_operator< Operator >::value, - void >::type * const = nullptr + !is_object< OutputType >::value && + !is_object< InputType1 >::value && + !is_object< InputType2 >::value && + is_operator< Operator >::value + >::type * const = nullptr ) { - assert( !(descr & descriptors::force_row_major ) ); - static_assert( allow_void || - ( !( - std::is_same< InputType1, void >::value || - std::is_same< InputType2, void >::value - ) ), - "grb::internal::eWiseApply_matrix_generic: the non-monoid version of " - "elementwise mxm can only be used if neither of the input matrices " - "is a pattern matrix (of type void)" ); - assert( phase != TRY ); - #ifdef _DEBUG - std::cout << "In grb::internal::eWiseApply_matrix_generic\n"; + std::cout << "In grb::internal::eWiseApply_matrix_generic_intersection\n"; #endif + assert( phase != TRY ); + constexpr bool crs_only = descr & descriptors::force_row_major; // get whether the matrices should be transposed prior to execution constexpr bool trans_left = descr & descriptors::transpose_left; constexpr bool trans_right = descr & descriptors::transpose_right; // run-time checks - const size_t m = grb::nrows( C ); - const size_t n = grb::ncols( C ); - const size_t m_A = !trans_left ? grb::nrows( A ) : grb::ncols( A ); - const size_t n_A = !trans_left ? grb::ncols( A ) : grb::nrows( A ); - const size_t m_B = !trans_right ? grb::nrows( B ) : grb::ncols( B ); - const size_t n_B = !trans_right ? grb::ncols( B ) : grb::nrows( B ); + const size_t m = nrows( C ); + const size_t n = ncols( C ); + const size_t m_A = !trans_left ? nrows( A ) : ncols( A ); + const size_t n_A = !trans_left ? ncols( A ) : nrows( A ); + const size_t m_B = !trans_right ? nrows( B ) : ncols( B ); + const size_t n_B = !trans_right ? ncols( B ) : nrows( B ); + + if( crs_only && (trans_left || trans_right) ) { +#ifdef _DEBUG + std::cerr << "grb::descriptors::force_row_major and " + << "grb::descriptors::transpose_left/right are mutually " + << "exclusive\n"; +#endif + return ILLEGAL; + } if( m != m_A || m != m_B || n != n_A || n != n_B ) { +#ifdef _DEBUG + std::cerr << "grb::eWiseApply: dimensions of input matrices do not match\n"; +#endif return MISMATCH; } + if( getID(A) == getID(C) || getID(B) == getID(C) ) { +#ifdef _DEBUG + std::cerr << "grb::eWiseApply: The output matrix can not simultaneously be " + << "one of the input matrices\n"; +#endif + } + const auto &A_raw = !trans_left ? internal::getCRS( A ) : internal::getCCS( A ); const auto &B_raw = !trans_right ? internal::getCRS( B ) : internal::getCCS( B ); - auto &C_raw = internal::getCRS( C ); + auto &CRS_raw = internal::getCRS( C ); auto &CCS_raw = internal::getCCS( C ); - -#ifdef _DEBUG - std::cout << "\t\t A offset array = { "; - for( size_t i = 0; i <= m_A; ++i ) { - std::cout << A_raw.col_start[ i ] << " "; - } - std::cout << "}\n"; - for( size_t i = 0; i < m_A; ++i ) { - for( size_t k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { - std::cout << "\t\t ( " << i << ", " << A_raw.row_index[ k ] << " ) = " - << A_raw.getPrintValue( k ) << "\n"; - } - } - std::cout << "\t\t B offset array = { "; - for( size_t j = 0; j <= m_B; ++j ) { - std::cout << B_raw.col_start[ j ] << " "; - } - std::cout << "}\n"; - for( size_t j = 0; j < m_B; ++j ) { - for( size_t k = B_raw.col_start[ j ]; k < B_raw.col_start[ j + 1 ]; ++k ) { - std::cout << "\t\t ( " << B_raw.row_index[ k ] << ", " << j << " ) = " - << B_raw.getPrintValue( k ) << "\n"; - } - } -#endif + const auto dummy_identity = identities::zero< OutputType >::value(); // retrieve buffers - char * arr1, * arr2, * arr3, * buf1, * buf2, * buf3; - arr1 = arr2 = buf1 = buf2 = nullptr; + char * arr1 = nullptr, * arr3 = nullptr, * buf1 = nullptr, * buf3 = nullptr; InputType1 * vbuf1 = nullptr; - InputType2 * vbuf2 = nullptr; OutputType * valbuf = nullptr; internal::getMatrixBuffers( arr1, buf1, vbuf1, 1, A ); - internal::getMatrixBuffers( arr2, buf2, vbuf2, 1, B ); internal::getMatrixBuffers( arr3, buf3, valbuf, 1, C ); // end buffer retrieval - // initialisations - internal::Coordinates< reference > coors1, coors2; + // initialisations of the coordinates + Coordinates coors1; coors1.set( arr1, false, buf1, n ); - coors2.set( arr2, false, buf2, n ); + // end initialisations of the coordinates + + if( !crs_only ) { #ifdef _H_GRB_REFERENCE_OMP_BLAS3 - #pragma omp parallel - { - size_t start, end; - config::OMP::localRange( start, end, 0, n + 1 ); -#else - const size_t start = 0; - const size_t end = n + 1; + #pragma omp parallel for simd #endif - for( size_t j = start; j < end; ++j ) { + for( size_t j = 0; j <= n; ++j ) { CCS_raw.col_start[ j ] = 0; } -#ifdef _H_GRB_REFERENCE_OMP_BLAS3 } -#endif // end initialisations // nonzero count @@ -1054,6 +1020,7 @@ namespace grb { // symbolic phase if( phase == RESIZE ) { + nzc = 0; for( size_t i = 0; i < m; ++i ) { coors1.clear(); for( size_t k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { @@ -1063,22 +1030,33 @@ namespace grb { for( size_t l = B_raw.col_start[ i ]; l < B_raw.col_start[ i + 1 ]; ++l ) { const size_t l_col = B_raw.row_index[ l ]; if( coors1.assigned( l_col ) ) { - (void)++nzc; + (void) ++nzc; } } } const RC ret = grb::resize( C, nzc ); - if( ret != SUCCESS ) { - return ret; - } +#ifdef _DEBUG + std::cout << "grb::resize( C, " << nzc << " ) = " << ret << "\n"; +#endif + return ret; } // computational phase if( phase == EXECUTE ) { + nzc = 0; // retrieve additional buffer - config::NonzeroIndexType * const C_col_index = internal::template - getReferenceBuffer< typename config::NonzeroIndexType >( n + 1 ); + config::NonzeroIndexType * const C_col_index = + getReferenceBuffer< config::NonzeroIndexType >( n + 1 ); + + if( !crs_only ) { +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for simd +#endif + for( size_t j = 0; j < n+1; ++j ) { + C_col_index[ j ] = 0; + } + } // perform column-wise nonzero count for( size_t i = 0; i < m; ++i ) { @@ -1091,7 +1069,9 @@ namespace grb { const size_t l_col = B_raw.row_index[ l ]; if( coors1.assigned( l_col ) ) { (void) ++nzc; - (void) ++CCS_raw.col_start[ l_col + 1 ]; + if( !crs_only ) { + (void) ++CCS_raw.col_start[ l_col + 1 ]; + } } } } @@ -1111,92 +1091,373 @@ namespace grb { } // prefix sum for CCS_raw.col_start - assert( CCS_raw.col_start[ 0 ] == 0 ); - for( size_t j = 1; j < n; ++j ) { - CCS_raw.col_start[ j + 1 ] += CCS_raw.col_start[ j ]; - } - assert( CCS_raw.col_start[ n ] == nzc ); - - // set C_col_index to all zero -#ifdef _H_GRB_REFERENCE_OMP_BLAS3 - #pragma omp parallel - { - size_t start, end; - config::OMP::localRange( start, end, 0, n ); -#else - const size_t start = 0; - const size_t end = n; -#endif - for( size_t j = start; j < end; ++j ) { - C_col_index[ j ] = 0; + if( !crs_only ) { + assert( CCS_raw.col_start[ 0 ] == 0 ); + for( size_t j = 1; j < n; ++j ) { + CCS_raw.col_start[ j + 1 ] += CCS_raw.col_start[ j ]; } -#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + assert( CCS_raw.col_start[ n ] == nzc ); } -#endif // do computations - size_t nzc = 0; - C_raw.col_start[ 0 ] = 0; + nzc = 0; + CRS_raw.col_start[ 0 ] = 0; for( size_t i = 0; i < m; ++i ) { coors1.clear(); - coors2.clear(); + for( size_t k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { + const auto k_col = A_raw.row_index[ k ]; + if( !coors1.assign( k_col ) ) { + utils::assignValue( valbuf, k_col, A_raw.getValue( k, dummy_identity ) ); + } + } + + for( size_t l = B_raw.col_start[ i ]; l < B_raw.col_start[ i + 1 ]; ++l ) { + const auto j = B_raw.row_index[ l ]; + if( !coors1.assigned( j ) ) { // Union case: ignored + continue; + } + + const auto valbuf_value_before = valbuf[ j ]; + OutputType result_value; + (void) grb::apply( result_value, valbuf_value_before, B_raw.getValue( l, dummy_identity ), oper ); + + // update CRS + CRS_raw.row_index[ nzc ] = j; + CRS_raw.setValue( nzc, result_value ); + + // update CCS + if( !crs_only ) { + ++C_col_index[ j ]; + const size_t CCS_index = CCS_raw.col_start[ j+1 ] - C_col_index[ j ]; + CCS_raw.row_index[ CCS_index ] = i; + CCS_raw.setValue( CCS_index, result_value ); + } + + // update count + (void) ++nzc; + } + + CRS_raw.col_start[ i + 1 ] = nzc; + } + +#ifndef NDEBUG + if( !crs_only ) { + for( size_t j = 0; j < n; ++j ) { + assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] == C_col_index[ j ] ); + } + } +#endif + + // set final number of nonzeroes in output matrix #ifdef _DEBUG - std::cout << "\t The elements "; + std::cout << "internal::setCurrentNonzeroes( C, " << nzc << " )\n"; #endif - for( size_t k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { - const size_t k_col = A_raw.row_index[ k ]; - coors1.assign( k_col ); - valbuf[ k_col ] = A_raw.getValue( k, - mulMonoid.template getIdentity< typename Operator::D1 >() ); + internal::setCurrentNonzeroes( C, nzc ); + } + + // done + return SUCCESS; + } + + template< + Descriptor descr, + class Monoid, + typename OutputType, typename InputType1, typename InputType2, + typename RIT1, typename CIT1, typename NIT1, + typename RIT2, typename CIT2, typename NIT2, + typename RIT3, typename CIT3, typename NIT3 + > + RC eWiseApply_matrix_generic_union( + Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C, + const Matrix< InputType1, reference, RIT2, CIT2, NIT2 > &A, + const Matrix< InputType2, reference, RIT3, CIT3, NIT3 > &B, + const Monoid &monoid, + const Phase &phase, + const typename std::enable_if< + !is_object< OutputType >::value && + !is_object< InputType1 >::value && + !is_object< InputType2 >::value && + is_monoid< Monoid >::value + >::type * const = nullptr + ) { + #ifdef _DEBUG - std::cout << "A( " << i << ", " << k_col << " ) = " << A_raw.getValue( k, - mulMonoid.template getIdentity< typename Operator::D1 >() ) << ", "; + std::cout << "In grb::internal::eWiseApply_matrix_generic_union\n"; #endif - } + assert( phase != TRY ); + constexpr bool crs_only = descr & descriptors::force_row_major; + // get whether the matrices should be transposed prior to execution + constexpr bool trans_left = descr & descriptors::transpose_left; + constexpr bool trans_right = descr & descriptors::transpose_right; + + if( crs_only && (trans_left || trans_right) ) { #ifdef _DEBUG - std::cout << "are multiplied pairwise with "; + std::cerr << "grb::descriptors::force_row_major and " + << "grb::descriptors::transpose_left/right are mutually " + << "exclusive\n"; #endif + return ILLEGAL; + } + + if( getID(A) == getID(C) || getID(B) == getID(C) ) { +#ifdef _DEBUG + std::cerr << "grb::eWiseApply: The output matrix can not simultaneously be " + << "one of the input matrices\n"; +#endif + } + + // run-time checks + const size_t m = nrows( C ); + const size_t n = ncols( C ); + const size_t m_A = !trans_left ? nrows( A ) : ncols( A ); + const size_t n_A = !trans_left ? ncols( A ) : nrows( A ); + const size_t m_B = !trans_right ? nrows( B ) : ncols( B ); + const size_t n_B = !trans_right ? ncols( B ) : nrows( B ); + + // Identities + const auto identity_A = monoid.template getIdentity< OutputType >(); + const auto identity_B = monoid.template getIdentity< OutputType >(); + + if( m != m_A || m != m_B || n != n_A || n != n_B ) { +#ifdef _DEBUG + std::cerr << "grb::eWiseApply: dimensions of input matrices do not match\n"; +#endif + return MISMATCH; + } + + const auto oper = monoid.getOperator(); + const auto &A_raw = !trans_left ? + internal::getCRS( A ) : + internal::getCCS( A ); + const auto &B_raw = !trans_right ? + internal::getCRS( B ) : + internal::getCCS( B ); + auto &CRS_raw = internal::getCRS( C ); + auto &CCS_raw = internal::getCCS( C ); + + // retrieve buffers + char *arr1 = nullptr, *arr3 = nullptr; + char *buf1 = nullptr, *buf3 = nullptr; + InputType1 * vbuf1 = nullptr; + OutputType * vbuf3 = nullptr; + internal::getMatrixBuffers( arr1, buf1, vbuf1, 1, A ); + internal::getMatrixBuffers( arr3, buf3, vbuf3, 1, C ); + // end buffer retrieval + + // initialisations of the coordinates + // Note: By using the buffer of the output matrix C, we can + // allow A and B to be the same matrix (with the same buffer) + Coordinates< reference > coors1, coors2; + coors1.set( arr1, false, buf1, n ); + coors2.set( arr3, false, buf3, n ); + // end initialisations of the coordinates + + if( !crs_only ) { +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for simd +#endif + for( size_t j = 0; j < n + 1; ++j ) { + CCS_raw.col_start[ j ] = 0; + } + } + // end initialisations + + // nonzero count + size_t nzc = 0; + + // symbolic phase + if( phase == RESIZE ) { + nzc = 0; + for( size_t i = 0; i < m; ++i ) { + coors1.clear(); + for( size_t k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { + const auto k_col = A_raw.row_index[ k ]; + if( !coors1.assign( k_col ) ) { + (void) ++nzc; + } + } for( size_t l = B_raw.col_start[ i ]; l < B_raw.col_start[ i + 1 ]; ++l ) { const size_t l_col = B_raw.row_index[ l ]; - if( coors1.assigned( l_col ) ) { - coors2.assign( l_col ); - (void)grb::apply( valbuf[ l_col ], valbuf[ l_col ], B_raw.getValue( l, - mulMonoid.template getIdentity< typename Operator::D2 >() ), oper ); + if( !coors1.assigned( l_col ) ) { + (void) ++nzc; + } + } + } + + const RC ret = grb::resize( C, nzc ); #ifdef _DEBUG - std::cout << "B( " << i << ", " << l_col << " ) = " << B_raw.getValue( l, - mulMonoid.template getIdentity< typename Operator::D2 >() ) - << " to yield C( " << i << ", " << l_col << " ), "; + std::cout << "grb::resize( C, " << nzc << " ) = " << ret << "\n"; #endif + return ret; + } + + // computational phase + if( phase == EXECUTE ) { + // retrieve additional buffer + config::NonzeroIndexType * const C_col_index = + getReferenceBuffer< config::NonzeroIndexType >( n + 1 ); + + // perform column-wise nonzero count + nzc = 0; + for( size_t i = 0; i < m; ++i ) { + coors1.clear(); + for( size_t k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { + const size_t k_col = A_raw.row_index[ k ]; + (void) coors1.assign( k_col ); + (void) ++nzc; + + if( !crs_only ) { + (void) ++CCS_raw.col_start[ k_col + 1 ]; + } + } + for( size_t l = B_raw.col_start[ i ]; l < B_raw.col_start[ i + 1 ]; ++l ) { + const size_t l_col = B_raw.row_index[ l ]; + if( !coors1.assigned( l_col ) ) { + (void) ++nzc; + if( !crs_only ) { + (void) ++CCS_raw.col_start[ l_col + 1 ]; + } } } + } + + // check capacity + if( nzc > capacity( C ) ) { #ifdef _DEBUG - std::cout << "\n"; + std::cout << "\t detected insufficient capacity " + << "for requested operation\n"; #endif - for( size_t k = 0; k < coors2.nonzeroes(); ++k ) { - const size_t j = coors2.index( k ); + const RC clear_rc = clear( C ); + if( clear_rc != SUCCESS ) { + return PANIC; + } + return FAILED; + } + + // prefix sum for CCS_raw.col_start + if( !crs_only ) { + for( size_t j = 1; j < n; ++j ) { + CCS_raw.col_start[ j + 1 ] += CCS_raw.col_start[ j ]; + } +#ifndef NDEBUG + assert( CCS_raw.col_start[ 0 ] == 0 ); + assert( CCS_raw.col_start[ n ] == nzc ); +#endif + } + + // set C_col_index to all zero + if( !crs_only ) { +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel for simd +#endif + for( size_t j = 0; j < n; j++ ) { + C_col_index[ j ] = 0; + } + } + + // do computations + nzc = 0; + CRS_raw.col_start[ 0 ] = 0; + for( size_t i = 0; i < m; ++i ) { +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp parallel +#endif + { +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp single nowait +#endif + { + coors1.clear(); + for( size_t k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) { + const auto k_col = A_raw.row_index[ k ]; + (void) coors1.assign( k_col ); + utils::assignValue( vbuf1, k_col, A_raw.getValue( k, identity_A ) ); + } + } +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + #pragma omp single nowait +#endif + { + coors2.clear(); + for( size_t k = B_raw.col_start[ i ]; k < B_raw.col_start[ i + 1 ]; ++k ) { + const auto k_col = B_raw.row_index[ k ]; + (void) coors2.assign( k_col ); + utils::assignValue( vbuf3, k_col, B_raw.getValue( k, identity_B ) ); + } + } + } + + for( size_t k = 0; k < coors1.nonzeroes(); ++k ) { + const auto j = coors1.index( k ); + const auto A_val = utils::getValue(vbuf1, j, identity_A); + const auto B_val = coors2.assigned(j) ? utils::getValue(vbuf3, j, identity_B) : identity_B; + + OutputType result_value; + (void) grb::apply( result_value, A_val, B_val, oper ); + // update CRS - C_raw.row_index[ nzc ] = j; - C_raw.setValue( nzc, valbuf[ j ] ); + CRS_raw.row_index[ nzc ] = j; + CRS_raw.setValue( nzc, result_value ); + // update CCS - const size_t CCS_index = C_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 = CCS_raw.col_start[ j+1 ] - ++C_col_index[ j ]; +#ifdef NDEBUG + assert( CCS_index < capacity( C ) ); + assert( CCS_index < CCS_raw.col_start[ j+1 ] ); + assert( CCS_index >= CCS_raw.col_start[ j ] ); +#endif + CCS_raw.row_index[ CCS_index ] = i; + CCS_raw.setValue( CCS_index, result_value ); + } // update count - (void)++nzc; + (void) ++nzc; } - C_raw.col_start[ i + 1 ] = nzc; -#ifdef _DEBUG - std::cout << "\n"; + for( size_t k = 0; k < coors2.nonzeroes(); ++k ) { + const auto j = coors2.index( k ); + if( coors1.assigned(j) ) { // Intersection case: already handled + continue; + } + const auto A_val = coors1.assigned(j) ? utils::getValue(vbuf1, j, identity_A) : identity_A; + const auto B_val = utils::getValue(vbuf3, j, identity_B); + + OutputType result_value; + (void) grb::apply( result_value, A_val, B_val, oper ); + + // update CRS + CRS_raw.row_index[ nzc ] = j; + CRS_raw.setValue( nzc, result_value ); + + // update CCS + if( !crs_only ) { + const size_t CCS_index = CCS_raw.col_start[ j+1 ] - ++C_col_index[ j ]; +#ifdef NDEBUG + assert( CCS_index < capacity( C ) ); + assert( CCS_index < CCS_raw.col_start[ j+1 ] ); + assert( CCS_index >= CCS_raw.col_start[ j ] ); #endif + CCS_raw.row_index[ CCS_index ] = i; + CCS_raw.setValue( CCS_index, result_value ); + } + // update count + (void) ++nzc; + } + + CRS_raw.col_start[ i + 1 ] = nzc; } #ifndef NDEBUG - for( size_t j = 0; j < n; ++j ) { - assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] == C_col_index[ j ] ); + if( !crs_only ) { + for( size_t j = 0; j < n; ++j ) { + assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] == C_col_index[ j ] ); + } } #endif // set final number of nonzeroes in output matrix +#ifdef _DEBUG + std::cout << "internal::setCurrentNonzeroes( C, " << nzc << " )\n"; +#endif internal::setCurrentNonzeroes( C, nzc ); } @@ -1207,11 +1468,11 @@ namespace grb { } // namespace internal /** - * Computes \f$ C = A . B \f$ for a given monoid. + * Computes \f$ C = A . B \f$ for a given monoid (union pattern). * * \internal Allows pattern matrix inputs. * - * \internal Dispatches to internal::eWiseApply_matrix_generic + * \internal Dispatches to internal::eWiseApply_matrix_generic_union */ template< Descriptor descr = descriptors::no_operation, @@ -1226,12 +1487,14 @@ namespace grb { const Matrix< InputType1, reference, RIT2, CIT2, NIT2 > &A, const Matrix< InputType2, reference, RIT3, CIT3, NIT3 > &B, const MulMonoid &mulmono, - const Phase phase = EXECUTE, - const typename std::enable_if< !grb::is_object< OutputType >::value && - !grb::is_object< InputType1 >::value && - !grb::is_object< InputType2 >::value && - grb::is_monoid< MulMonoid >::value, - void >::type * const = nullptr + const Phase &phase = EXECUTE, + const typename std::enable_if< + !grb::is_object< OutputType >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + grb::is_monoid< MulMonoid >::value, + void + >::type * const = nullptr ) { // static checks NO_CAST_ASSERT( ( !( descr & descriptors::no_casting ) || @@ -1254,22 +1517,21 @@ namespace grb { ); #ifdef _DEBUG - std::cout << "In grb::eWiseApply_matrix_generic (reference, monoid)\n"; + std::cout << "In grb::eWiseApply( reference, monoid )\n"; #endif - return internal::eWiseApply_matrix_generic< true, descr >( - C, A, B, mulmono.getOperator(), mulmono, phase + return internal::eWiseApply_matrix_generic_union< descr >( + C, A, B, mulmono, phase ); } /** - * Computes \f$ C = A . B \f$ for a given binary operator. + * Computes \f$ C = A . B \f$ for a given operator (intersection pattern). * - * \internal Pattern matrices not allowed + * \internal Allows pattern matrix inputs. * - * \internal Dispatches to internal::eWiseApply_matrix_generic + * \internal Dispatches to internal::eWiseApply_matrix_generic_intersection */ - template< Descriptor descr = grb::descriptors::no_operation, class Operator, @@ -1282,23 +1544,33 @@ namespace grb { Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C, const Matrix< InputType1, reference, RIT2, CIT2, NIT2 > &A, const Matrix< InputType2, reference, RIT3, CIT3, NIT3 > &B, - const Operator &mulOp, - const Phase phase = EXECUTE, - const typename std::enable_if< !grb::is_object< OutputType >::value && - !grb::is_object< InputType1 >::value && - !grb::is_object< InputType2 >::value && - grb::is_operator< Operator >::value, - void >::type * const = nullptr + const Operator &op, + const Phase &phase = EXECUTE, + const typename std::enable_if< + !grb::is_object< OutputType >::value && + !grb::is_object< InputType1 >::value && + !grb::is_object< InputType2 >::value && + grb::is_operator< Operator >::value, + void + >::type * const = nullptr ) { + typedef typename std::conditional< + std::is_void::value, + typename Operator::D1, + InputType1>::type ActualInputType1; + typedef typename std::conditional< + std::is_void::value, + typename Operator::D2, + InputType1>::type ActualInputType2; // static checks NO_CAST_ASSERT( ( !( descr & descriptors::no_casting ) || - std::is_same< typename Operator::D1, InputType1 >::value ), + std::is_same< typename Operator::D1, ActualInputType1 >::value ), "grb::eWiseApply (reference, matrix <- matrix x matrix, operator)", "called with a prefactor input matrix A 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 ), + std::is_same< typename Operator::D2, ActualInputType2 >::value ), "grb::eWiseApply (reference, matrix <- matrix x matrix, operator)", "called with a postfactor input matrix B that does not match the first " "domain of the given multiplication operator" @@ -1309,20 +1581,19 @@ namespace grb { "called with an output matrix C that does not match the output domain " "of the given multiplication operator" ); - static_assert( ( !( - std::is_same< InputType1, void >::value || - std::is_same< InputType2, void >::value ) - ), "grb::eWiseApply (reference, matrix <- matrix x matrix, operator): " - "the operator version of eWiseApply cannot be used if either of the " - "input matrices is a pattern matrix (of type void)" + static_assert( + !std::is_void< OutputType >::value || + ( std::is_void< InputType1 >::value && std::is_void< InputType2 >::value ), + "grb::eWiseApply: the elementwise mxm only support" + " output pattern-matrix (of type void) if both" + " input matrices are also pattern matrices" ); +#ifdef _DEBUG + std::cout << "In grb::eWiseApply( reference, operator )\n"; +#endif - typename grb::Monoid< - grb::operators::mul< double >, - grb::identities::one - > dummyMonoid; - return internal::eWiseApply_matrix_generic< false, descr >( - C, A, B, mulOp, dummyMonoid, phase + return internal::eWiseApply_matrix_generic_intersection< descr >( + C, A, B, op, phase ); } diff --git a/include/graphblas/utils.hpp b/include/graphblas/utils.hpp index c5239afdc..a82aa3125 100644 --- a/include/graphblas/utils.hpp +++ b/include/graphblas/utils.hpp @@ -54,6 +54,24 @@ namespace grb { */ namespace utils { + template< typename D, typename T > + static void assignValue( + D *array, size_t i, const T& value, + typename std::enable_if< !std::is_void< D >::value >::type * const = nullptr + ) { array[i] = value; } + + template< typename T > + static void assignValue( void *, size_t, const T& ) { /* do nothing */ } + + template< typename D, typename T > + static T getValue( + const D *array, size_t i, const T&, + typename std::enable_if< !std::is_void< D >::value >::type * const = nullptr + ) { return array[i]; } + + template< typename T > + static T getValue( const void *, size_t, const T& identity ) { return identity; } + /** * Checks whether two values are equal. * diff --git a/tests/performance/bench_kernels.h b/tests/performance/bench_kernels.h index 8fea223b3..29901a1dc 100644 --- a/tests/performance/bench_kernels.h +++ b/tests/performance/bench_kernels.h @@ -62,7 +62,7 @@ void bench_kernels_axpy( const size_t n ); -/** +/** * Executes the inner-product computation \f$ alpha = (x,y) \f$ with \a x and * \a y vectors of length \a n. * diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index ab0dac498..bb568a427 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -89,6 +89,11 @@ add_grb_executables( ewiseapply ewiseapply.cpp BACKENDS reference reference_omp bsp1d hybrid hyperdags nonblocking ) +add_grb_executables( eWiseApplyMatrix_variants eWiseApplyMatrix_variants.cpp + BACKENDS reference reference_omp bsp1d hybrid hyperdags nonblocking + ADDITIONAL_LINK_LIBRARIES test_utils_headers +) + add_grb_executables( eWiseMatrix eWiseMatrix.cpp BACKENDS reference reference_omp bsp1d hybrid hyperdags nonblocking ) diff --git a/tests/unit/eWiseApplyMatrixReference.cpp b/tests/unit/eWiseApplyMatrixReference.cpp index aac355450..1a1982b5d 100644 --- a/tests/unit/eWiseApplyMatrixReference.cpp +++ b/tests/unit/eWiseApplyMatrixReference.cpp @@ -15,255 +15,350 @@ * limitations under the License. */ -#include -#include +#include #include +#include +#include #include -#ifdef _DEBUG - #include -#endif +using namespace grb; // static data corresponding to small matrices -static const size_t I_A[] = { 0, 0, 1, 1, 2, 2, 3, 3 }; -static const size_t J_A[] = { 0, 2, 1, 2, 2, 3, 0, 2 }; -static const double V_A[] = { 1, 3, 4, 2, 6, 7, 5, 8 }; - -static const size_t I_B[] = { 0, 0, 1, 2, 3, 3 }; -static const size_t J_B[] = { 0, 3, 1, 1, 2, 3 }; -static const double V_B[] = { 9, 10, 11, 12, 14, 13 }; - -static const size_t I_C[] = { 0, 1, 3 }; -static const size_t J_C[] = { 0, 1, 2 }; -static const double V_C[] = { 9, 44, 112 }; - -static const size_t rowlens[] = { 1, 1, 0, 1 }; -static const size_t collens[] = { 1, 1, 1, 0 }; - -static const double expect1_CRS[] = { 9, 44, 112 }; -static const double expect1_CCS[] = { 9, 44, 112 }; - -static const double expect2_CRS[] = { 1, 4, 8 }; -static const double expect2_CCS[] = { 1, 4, 8 }; +/** + * A: + * 1 _ 3 _ + * _ 4 2 _ + * _ _ 6 7 + * 5 _ _ 8 + */ +static const std::vector< size_t > I_A { 0, 0, 1, 1, 2, 2, 3, 3 }; +static const std::vector< size_t > J_A { 0, 2, 1, 2, 2, 3, 0, 2 }; +static const std::vector< int > V_A { 1, 3, 4, 2, 6, 7, 5, 8 }; + +/** + * B: + * 9 __ __ __ + * __ 11 12 __ + * __ 14 __ __ + * __ __ __ 13 + */ +static const std::vector< size_t > I_B { 0, 0, 1, 2, 3, 3 }; +static const std::vector< size_t > J_B { 0, 3, 1, 1, 2, 3 }; +static const std::vector< int > V_B { 9, 10, 11, 12, 14, 13 }; + +/** + * C_intersection: + * 9 ___ ___ ___ + * ___ 44 ___ ___ + * ___ ___ ___ ___ + * ___ ___ 112 ___ + */ +static const std::vector< size_t > I_C_intersection { 0, 1, 3 }; +static const std::vector< size_t > J_C_intersection { 0, 1, 2 }; +static const std::vector< int > V_C_intersection { 9, 44, 112 }; + +/** + * C_union_A_B: + * 9 ___ 3 10 + * ___ 44 2 ___ + * ___ 12 6 7 + * 5 ___ 112 13 + */ -static const double expect3_CRS[] = { 9, 11, 14 }; -static const double expect3_CCS[] = { 9, 11, 14 }; +static const std::vector< size_t > I_C_union { 0, 0, 0, 1, 1, 2, 2, 2, 3, 3, 3 }; +static const std::vector< size_t > J_C_union { 0, 2, 3, 1, 2, 1, 2, 3, 0, 2, 3 }; +static const std::vector< int > V_C_union_A_B { 9, 3, 10, 44, 2, 12, 6, 7, 5, 112, 13 }; -static const double expect4_CRS[] = { 1, 1, 1 }; -static const double expect4_CCS[] = { 1, 1, 1 }; +/** + * C_union_A_B_pattern: + * 1 _ 3 1 + * _ 4 2 _ + * _ 1 6 7 + * 5 _ 8 1 + */ +static const std::vector< size_t > I_C_union_A_B_pattern { 0, 0, 0, 1, 1, 2, 2, 2, 3, 3, 3 }; +static const std::vector< size_t > J_C_union_A_B_pattern { 0, 2, 3, 1, 2, 1, 2, 3, 0, 2, 3 }; +static const std::vector< int > V_C_union_A_B_pattern { 1, 3, 1, 4, 2, 1, 6, 7, 5, 8, 1 }; + +/** + * C_union_A_pattern_B: + * 9 __ 1 10 + * __ 11 1 __ + * __ 12 1 1 + * 1 __ 14 13 + */ +static const std::vector< size_t > I_C_union_A_pattern_B { 0, 0, 0, 1, 1, 2, 2, 2, 3, 3, 3 }; +static const std::vector< size_t > J_C_union_A_pattern_B { 0, 2, 3, 1, 2, 1, 2, 3, 0, 2, 3 }; +static const std::vector< int > V_C_union_A_pattern_B { 9, 1, 10, 11, 1, 12, 1, 1, 1, 14, 13 }; + +/** + * C_union_A_pattern_B_pattern: + * 1 _ 1 1 + * _ 1 1 _ + * _ 1 1 1 + * 1 _ 1 1 + */ +static const std::vector< size_t > I_C_union_A_pattern_B_pattern { 0, 0, 0, 1, 1, 2, 2, 2, 3, 3, 3 }; +static const std::vector< size_t > J_C_union_A_pattern_B_pattern { 0, 2, 3, 1, 2, 1, 2, 3, 0, 2, 3 }; +static const std::vector< int > V_C_union_A_pattern_B_pattern { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }; // helper function to check internal data structures // of the reference backend template< typename T > -void checkCRSandCCS( const grb::Matrix< T > &C, - const size_t n, - const size_t * rlens, - const size_t * clens, - const size_t * I, - const size_t * J, - const double * expect_CRS, - const double * expect_CCS, - grb::RC &rc +void checkCRSandCCS( + const Matrix< T > & obtained, + const Matrix< T > & expected, + RC & rc ) { - // check CRS output - const auto &crs1 = grb::internal::getCRS( C ); - for( size_t i = 0; i < n; ++i ) { - const size_t entries = crs1.col_start[ i + 1 ] - crs1.col_start[ i ]; - if( entries != rlens[ i ] ) { - std::cerr << "Error: unexpected number of entries " << entries << ", " - << " expected " << rlens[ i ] << " (CRS).\n"; - rc = grb::FAILED; - } - for( size_t k = crs1.col_start[ i ]; k < crs1.col_start[ i + 1 ]; ++k ) { - if( crs1.row_index[ k ] != J[ k ] ) { - std::cerr << "Error: unexpected entry at ( " << i << ", " - << crs1.row_index[ k ] << " ), " - << "expected one at ( " << i << ", " << J[ k ] << " ) " - << "instead (CRS).\n"; - rc = grb::FAILED; - } - if( crs1.values[ k ] != expect_CRS[ k ] ) { - std::cerr << "Error: unexpected value " << crs1.values[ k ] << "; " - << "expected " << expect_CRS[ k ] << " (CRS).\n"; - rc = grb::FAILED; + if( nnz( obtained ) != nnz( expected ) ) { + std::cerr << "Error: unexpected number of non-zero entries; " + << "expected " << nnz( expected ) << ", " + << "obtained " << nnz( obtained ) << ".\n"; + rc = FAILED; + } + + { // check CRS output + const auto & crsObtained = internal::getCRS( obtained ); + const auto & crsExpected = internal::getCRS( expected ); + for( size_t i = 0; i < nrows( obtained ); ++i ) { + for( size_t k = crsObtained.col_start[ i ]; k < crsObtained.col_start[ i + 1 ]; ++k ) { + const auto nValuesInRow = crsObtained.col_start[ i + 1 ] - crsObtained.col_start[ i ]; + const auto expectedValuesInRow = crsExpected.col_start[ i + 1 ] - crsExpected.col_start[ i ]; + if( nValuesInRow != expectedValuesInRow ) { + std::cerr << "Error: unexpected number of non-zero entries in row " << i << "; " + << "expected " << expectedValuesInRow << ", " + << "obtained " << nValuesInRow << " (CRS).\n"; + rc = FAILED; + } + const auto searchedJ = crsObtained.row_index[ k ]; + const auto searchedV = crsObtained.values[ k ]; + bool found = false; + for( size_t l = crsExpected.col_start[ i ]; l < crsExpected.col_start[ i + 1 ]; ++l ) { + found |= ( crsExpected.row_index[ l ] == searchedJ ) && ( crsExpected.values[ l ] == searchedV ); + } + if( !found ) { + std::cerr << "Error: Can not found entry " + << "( " << i << ", " << searchedJ << " ) = " + << searchedV << " (CRS).\n"; + rc = FAILED; + } } } } - // check CCS output - const auto &ccs1 = grb::internal::getCCS( C ); - for( size_t j = 0; j < n; ++j ) { - const size_t entries = ccs1.col_start[ j + 1 ] - ccs1.col_start[ j ]; - if( entries != clens[ j ] ) { - std::cerr << "Error: unexpected number of entries " << entries << ", " - << "expected " << clens[ j ] << " (CCS).\n"; - rc = grb::FAILED; - } - for( size_t k = ccs1.col_start[ j ]; k < ccs1.col_start[ j + 1 ]; ++k ) { - if( ccs1.row_index[ k ] != I[ k ] ) { - std::cerr << "Error: unexpected entry at " - << "( " << ccs1.row_index[ k ] << ", " << j << " ), " - << "expected one at ( " << I[ k ] << ", " << j << " ) " - << "instead (CCS).\n"; - rc = grb::FAILED; - } - if( ccs1.values[ k ] != expect_CCS[ k ] ) { - std::cerr << "Error: unexpected value " << ccs1.values[ k ] << "; " - << "expected " << expect_CCS[ k ] << " (CCS).\n"; - rc = grb::FAILED; + { // check CCS output + const auto & ccsObtained = internal::getCCS( obtained ); + const auto & ccsExpected = internal::getCCS( expected ); + for( size_t i = 0; i < ncols( obtained ); ++i ) { + for( size_t k = ccsExpected.col_start[ i ]; k < ccsExpected.col_start[ i + 1 ]; ++k ) { + const auto nValuesInRow = ccsObtained.col_start[ i + 1 ] - ccsObtained.col_start[ i ]; + const auto expectedValuesInRow = ccsExpected.col_start[ i + 1 ] - ccsExpected.col_start[ i ]; + if( nValuesInRow != expectedValuesInRow ) { + std::cerr << "Error: unexpected number of non-zero entries in row " << i << "; " + << "expected " << expectedValuesInRow << ", " + << "obtained " << nValuesInRow << " (CCS).\n"; + rc = FAILED; + } + const auto searchedJ = ccsObtained.row_index[ k ]; + const auto searchedV = ccsObtained.values[ k ]; + bool found = false; + for( size_t l = ccsExpected.col_start[ i ]; l < ccsExpected.col_start[ i + 1 ]; ++l ) { + found |= ( ccsExpected.row_index[ l ] == searchedJ ) && ( ccsExpected.values[ l ] == searchedV ); + } + if( !found ) { + std::cerr << "Error: Can not found entry " + << "( " << i << ", " << searchedJ << " ) = " + << searchedV << " (CCS).\n"; + rc = FAILED; + } } } } } -void grbProgram( const void *, const size_t, grb::RC &rc ) { -#ifdef _DEBUG - constexpr const size_t smax = std::numeric_limits< size_t >::max(); -#endif +void grbProgram( const void *, const size_t, RC & rc ) { - // initialise test - grb::Monoid< grb::operators::mul< double >, grb::identities::one > mulmono; + // initialize test + const grb::Monoid< grb::operators::mul< int >, + grb::identities::one > mulmono; const size_t n = 4; const size_t nelts_A = 8; const size_t nelts_B = 6; - grb::Matrix< double > A( n, n ); - grb::Matrix< double > B( n, n ); - grb::Matrix< void > A_pattern( n, n ); - grb::Matrix< void > B_pattern( n, n ); - grb::Matrix< double > C( n, n ); - - rc = grb::resize( A, nelts_A ); - if( rc == grb::SUCCESS ) { - rc = grb::buildMatrixUnique( A, I_A, J_A, V_A, nelts_A, grb::SEQUENTIAL ); -#ifdef _DEBUG - print_matrix( A, smax, "A" ); -#endif - } - if( rc == grb::SUCCESS ) { - rc = grb::resize( B, nelts_B ); - } - if( rc == grb::SUCCESS ) { - rc = grb::buildMatrixUnique( B, I_B, J_B, V_B, nelts_B, grb::SEQUENTIAL ); -#ifdef _DEBUG - print_matrix( B, smax, "B" ); -#endif - } - if( rc == grb::SUCCESS ) { - rc = grb::resize( A_pattern, nelts_A ); - } - if( rc == grb::SUCCESS ) { - rc = grb::buildMatrixUnique( A_pattern, I_A, J_A, nelts_A, grb::SEQUENTIAL ); -#ifdef _DEBUG - print_matrix( A_pattern, smax, "A_pattern" ); -#endif - } - if( rc == grb::SUCCESS ) { - rc = grb::resize( B_pattern, nelts_B ); - } - if( rc == grb::SUCCESS ) { - rc = grb::buildMatrixUnique( B_pattern, I_B, J_B, nelts_B, grb::SEQUENTIAL ); -#ifdef _DEBUG - print_matrix( B_pattern, smax, "B_pattern" ); -#endif - } - if( rc != grb::SUCCESS ) { - std::cerr << "\tinitialisation FAILED\n"; - return; - } - - // test 1: compute with the monoid mxm_elementwise - std::cout << "\t Verifying the monoid version of mxm_elementwise, " - << "A and B value matrices\n"; - rc = grb::eWiseApply( C, A, B, mulmono, grb::RESIZE ); - rc = rc ? rc : grb::eWiseApply( C, A, B, mulmono ); - if( rc != grb::SUCCESS ) { - std::cerr << "Call to grb::eWiseApply FAILED\n"; - return; - } - - checkCRSandCCS( C, n, rowlens, collens, I_C, J_C, expect1_CRS, expect1_CCS, rc ); - - if( rc != grb::SUCCESS ) { - return; - } - - // test 2: compute with the monoid mxm_elementwise, A value matrix, B pattern matrix \n"; - std::cout << "\t Verifying the monoid version of mxm_elementwise, " - << "A value matrix, B pattern matrix\n"; - rc = grb::eWiseApply( C, A, B_pattern, mulmono, grb::RESIZE ); - rc = rc ? rc : grb::eWiseApply( C, A, B_pattern, mulmono ); - if( rc != grb::SUCCESS ) { - std::cerr << "Call to grb::eWiseApply FAILED\n"; - return; - } - - checkCRSandCCS( C, n, rowlens, collens, I_C, J_C, expect2_CRS, expect2_CCS, rc ); - - if( rc != grb::SUCCESS ) { - return; - } - - // test 3: compute with the monoid mxm_elementwise, A pattern matrix, B value matrix \n"; - std::cout << "\t Verifying the monoid version of mxm_elementwise, " - << "A pattern matrix, B value matrix\n"; - rc = grb::eWiseApply( C, A_pattern, B, mulmono, grb::RESIZE ); - rc = rc ? rc : grb::eWiseApply( C, A_pattern, B, mulmono ); - if( rc != grb::SUCCESS ) { - std::cerr << "Call to grb::eWiseApply FAILED\n"; - return; + Matrix< int > A( n, n ); + Matrix< int > B( n, n ); + Matrix< void > A_pattern( n, n ); + Matrix< void > B_pattern( n, n ); + Matrix< int > C( n, n ); + + assert( SUCCESS == resize( A, nelts_A ) ); + assert( SUCCESS == + buildMatrixUnique( A, I_A.data(), J_A.data(), V_A.data(), nelts_A, SEQUENTIAL ) + ); + assert( SUCCESS == resize( B, nelts_B ) ); + assert( SUCCESS == + buildMatrixUnique( B, I_B.data(), J_B.data(), V_B.data(), nelts_B, SEQUENTIAL ) + ); + assert( SUCCESS == resize( A_pattern, nelts_A ) ); + assert( SUCCESS == + buildMatrixUnique( A_pattern, I_A.data(), J_A.data(), nelts_A, SEQUENTIAL ) + ); + assert( SUCCESS == resize( B_pattern, nelts_B ) ); + assert( SUCCESS == + buildMatrixUnique( B_pattern, I_B.data(), J_B.data(), nelts_B, SEQUENTIAL ) + ); + + { // test 1: compute with the monoid mxm_elementwise + std::cout << "\t Verifying the monoid version of mxm_elementwise, " + << "A and B value matrices\n"; + clear( C ); + rc = grb::eWiseApply( C, A, B, mulmono, RESIZE ); + rc = rc ? rc : grb::eWiseApply( C, A, B, mulmono ); + if( rc != SUCCESS ) { + std::cerr << "Call to grb::eWiseApply FAILED\n"; + return; + } + Matrix< int > union_A_B( n, n ); + assert( SUCCESS == + buildMatrixUnique( + union_A_B, + I_C_union.data(), + J_C_union.data(), + V_C_union_A_B.data(), + I_C_union.size(), + SEQUENTIAL ) + ); + checkCRSandCCS( C, union_A_B, rc ); + + if( rc != SUCCESS ) { + return; + } } - checkCRSandCCS( C, n, rowlens, collens, I_C, J_C, expect3_CRS, expect3_CCS, rc ); - - if( rc != grb::SUCCESS ) { - return; + { // test 2: compute with the monoid mxm_elementwise, A value matrix, B pattern matrix \n"; + std::cout << "\t Verifying the monoid version of mxm_elementwise, " + << "A value matrix, B pattern matrix\n"; + clear( C ); + rc = grb::eWiseApply( C, A, B_pattern, mulmono, RESIZE ); + rc = rc ? rc : grb::eWiseApply( C, A, B_pattern, mulmono ); + if( rc != SUCCESS ) { + std::cerr << "Call to grb::eWiseApply FAILED\n"; + return; + } + Matrix< int > union_A_B_pattern( n, n ); + assert( SUCCESS == + buildMatrixUnique( + union_A_B_pattern, + I_C_union_A_B_pattern.data(), + J_C_union_A_B_pattern.data(), + V_C_union_A_B_pattern.data(), + I_C_union_A_B_pattern.size(), + SEQUENTIAL ) + ); + checkCRSandCCS( C, union_A_B_pattern, rc ); + + if( rc != SUCCESS ) { + return; + } } - // test 4: compute with the monoid mxm_elementwise, A pattern matrix, B pattern matrix \n"; - std::cout << "\t Verifying the monoid version of mxm_elementwise, " - << "A pattern matrix, B pattern matrix\n"; - rc = grb::eWiseApply( C, A_pattern, B_pattern, mulmono, grb::RESIZE ); - rc = rc ? rc : grb::eWiseApply( C, A_pattern, B_pattern, mulmono ); - if( rc != grb::SUCCESS ) { - std::cerr << "Call to grb::eWiseApply FAILED\n"; - return; + { // test 3: compute with the monoid mxm_elementwise, A pattern matrix, B value matrix \n"; + std::cout << "\t Verifying the monoid version of mxm_elementwise, " + << "A pattern matrix, B value matrix\n"; + clear( C ); + rc = grb::eWiseApply( C, A_pattern, B, mulmono, RESIZE ); + rc = rc ? rc : grb::eWiseApply( C, A_pattern, B, mulmono ); + if( rc != SUCCESS ) { + std::cerr << "Call to grb::eWiseApply FAILED\n"; + return; + } + Matrix< int > union_A_pattern_B( n, n ); + assert( SUCCESS == + buildMatrixUnique( + union_A_pattern_B, + I_C_union_A_pattern_B.data(), + J_C_union_A_pattern_B.data(), + V_C_union_A_pattern_B.data(), + I_C_union_A_pattern_B.size(), + SEQUENTIAL ) + ); + checkCRSandCCS( C, union_A_pattern_B, rc ); + + if( rc != SUCCESS ) { + return; + } } - checkCRSandCCS( C, n, rowlens, collens, I_C, J_C, expect4_CRS, expect4_CCS, rc ); - - if( rc != grb::SUCCESS ) { - return; + { // test 4: compute with the monoid mxm_elementwise, A pattern matrix, B pattern matrix \n"; + std::cout << "\t Verifying the monoid version of mxm_elementwise, " + << "A pattern matrix, B pattern matrix\n"; + clear( C ); + rc = grb::eWiseApply( C, A_pattern, B_pattern, mulmono, RESIZE ); + rc = rc ? rc : grb::eWiseApply( C, A_pattern, B_pattern, mulmono ); + if( rc != SUCCESS ) { + std::cerr << "Call to grb::eWiseApply FAILED\n"; + return; + } + Matrix< int > union_A_pattern_B_pattern( n, n ); + assert( SUCCESS == + buildMatrixUnique( + union_A_pattern_B_pattern, + I_C_union_A_pattern_B_pattern.data(), + J_C_union_A_pattern_B_pattern.data(), + V_C_union_A_pattern_B_pattern.data(), + I_C_union_A_pattern_B_pattern.size(), + SEQUENTIAL ) + ); + checkCRSandCCS( C, union_A_pattern_B_pattern, rc ); + + if( rc != SUCCESS ) { + return; + } } - // test 5: compute with the operator mxm_elementwise (pattern matrices not allowed) \n"; - std::cout << "\t Verifying the operator version of mxm_elementwise " - << "(only value matrices)\n"; - rc = grb::eWiseApply( C, A, B, mulmono.getOperator(), grb::RESIZE ); - rc = rc ? rc : grb::eWiseApply( C, A, B, mulmono.getOperator() ); - if( rc != grb::SUCCESS ) { - std::cerr << "Call to grb::eWiseApply FAILED\n"; - return; + { // test 5: compute with the operator mxm_elementwise (pattern matrices not allowed) \n"; + std::cout << "\t Verifying the operator version of mxm_elementwise " + << "(only value matrices)\n"; + clear( C ); + rc = grb::eWiseApply( C, A, B, mulmono.getOperator(), RESIZE ); + rc = rc ? rc : grb::eWiseApply( C, A, B, mulmono.getOperator() ); + if( rc != SUCCESS ) { + std::cerr << "Call to grb::eWiseApply FAILED\n"; + return; + } + Matrix< int > intersection_A_B( n, n ); + assert( SUCCESS == + buildMatrixUnique( + intersection_A_B, + I_C_intersection.data(), + J_C_intersection.data(), + V_C_intersection.data(), + I_C_intersection.size(), + SEQUENTIAL ) + ); + checkCRSandCCS( C, intersection_A_B, rc ); + + if( rc != SUCCESS ) { + return; + } } - - checkCRSandCCS( C, n, rowlens, collens, I_C, J_C, expect1_CRS, expect1_CCS, rc ); } int main( int argc, char ** argv ) { (void)argc; std::cout << "Functional test executable: " << argv[ 0 ] << "\n"; - grb::RC rc; + RC rc; grb::Launcher< grb::AUTOMATIC > launcher; - if( launcher.exec( &grbProgram, NULL, 0, rc ) != grb::SUCCESS ) { + if( launcher.exec( &grbProgram, NULL, 0, rc ) != SUCCESS ) { std::cerr << "Test failed to launch\n"; - rc = grb::FAILED; + rc = FAILED; } - if( rc == grb::SUCCESS ) { + if( rc == SUCCESS ) { std::cout << "Test OK\n" << std::endl; } else { std::cerr << std::flush; @@ -273,4 +368,3 @@ int main( int argc, char ** argv ) { // done return 0; } - diff --git a/tests/unit/eWiseApplyMatrix_variants.cpp b/tests/unit/eWiseApplyMatrix_variants.cpp new file mode 100644 index 000000000..8b8cc7437 --- /dev/null +++ b/tests/unit/eWiseApplyMatrix_variants.cpp @@ -0,0 +1,528 @@ + +/* + * 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. + */ + +/* + * @author Benjamin Lozes + * @date 24th of May, 2023 + * + * @brief Test for eWiseApply(Matrix, Monoid) + * and eWiseApply(Matrix, Operator) variants + * + * This test is meant to ensure the behaviour of the eWiseApply(Matrix, Monoid) + * and eWiseApply(Matrix, Operator) variants is correct. Precisely, we expect + * the following behaviour: + * - eWiseApply(Matrix, Monoid) should apply the monoid to all elements of + * the two matrices, INCLUDING the couples (non_zero, zero), using the + * provided identity value for the zero elements. + * - eWiseApply(Matrix, Operator) should apply the operator to all elements + * of the two matrices, EXCLUDING the couples (non_zero, zero) + * + */ + +#include +#include +#include +#include + +#include +#include + +using namespace grb; + +using nz_type = int; + +constexpr nz_type A_INITIAL_VALUE = 1; +constexpr nz_type B_INITIAL_VALUE = 3; + +// #define _DEBUG + + +template< typename D > +bool equals_matrix( + const Matrix< D > &A, + const Matrix< D > &B +) { + if( nrows( A ) != nrows( B ) || + ncols( A ) != ncols( B ) || + nnz( A ) != nnz( B ) + ) { + return false; + } + + wait( A ); + wait( B ); + + std::vector< + std::pair< std::pair< size_t, size_t >, D > + > A_vec( A.cbegin(), A.cend() ); + std::vector< + std::pair< std::pair< size_t, size_t >, D > + > B_vec( B.cbegin(), B.cend() ); + return std::is_permutation( A_vec.cbegin(), A_vec.cend(), B_vec.cbegin() ); +} + + +template< + class Monoid, + typename ValueTypeA, + typename ValueTypeB, + typename ValueTypeC, + Descriptor descr = descriptors::no_operation +> +struct input_t { + const Matrix< ValueTypeA > &A; + const Matrix< ValueTypeB > &B; + const Matrix< ValueTypeC > &C_monoid; + const Matrix< ValueTypeC > &C_operator; + const Monoid &monoid; + + input_t( + const Matrix< ValueTypeA > &A = {0,0}, + const Matrix< ValueTypeB > &B = {0,0}, + const Matrix< ValueTypeC > &C_monoid = {0,0}, + const Matrix< ValueTypeC > &C_operator = {0,0}, + const Monoid &monoid = Monoid() + ) : A( A ), + B( B ), + C_monoid( C_monoid ), + C_operator( C_operator ), + monoid( monoid ) {} +}; + + +struct output_t { + RC rc; +}; + +template< + class Monoid, + typename ValueTypeA, + typename ValueTypeB, + typename ValueTypeC, + Descriptor descr +> +void grb_program( + const input_t< Monoid, ValueTypeA, ValueTypeB, ValueTypeC, descr > &input, + output_t &output +) { + static_assert( is_monoid< Monoid >::value, "Monoid required" ); + const auto &op = input.monoid.getOperator(); + + RC &rc = output.rc; + + { // Operator variant + std::cout << " -- eWiseApply using Operator, supposed to be" + << " annihilating non-zeroes -> INTERSECTION\n"; + Matrix< ValueTypeC > C( nrows( input.A ), ncols( input.A ) ); + + rc = eWiseApply( C, input.A, input.B, op, RESIZE ); + if( rc != SUCCESS ) { + std::cerr << "Error: Phase::RESIZE\n"; + return; + } + if( capacity( C ) < nnz( input.C_operator ) ) { + std::cerr << "Error: Capacity should be at least " << nnz( input.C_operator ) << "\n"; + rc = FAILED; + return; + } + + rc = eWiseApply( C, input.A, input.B, op, EXECUTE ); + if( rc != SUCCESS ) { + std::cerr << "Error: Phase::EXECUTE\n"; + return; + } + if( !equals_matrix( C, input.C_operator ) ) { + std::cerr << "Error: Wrong result\n"; + rc = FAILED; + return; + } + + std::cout << "Result (operator) is correct\n"; + } + + { // Monoid variant + std::cout << " -- eWiseApply using Monoid, supposed to consider" + << " non-zeroes as the identity -> UNION\n"; + Matrix< ValueTypeC > C( nrows( input.A ), ncols( input.A ) ); + + rc = eWiseApply( C, input.A, input.B, input.monoid, RESIZE ); + if( rc != SUCCESS ) { + std::cerr << "Error: Phase::RESIZE\n"; + return; + } + if( capacity( C ) < nnz( input.C_operator ) ) { + std::cerr << "Error: Capacity should be at least " << nnz( input.C_monoid ) << "\n"; + rc = FAILED; + return; + } + + rc = eWiseApply( C, input.A, input.B, input.monoid, EXECUTE ); + if( rc != SUCCESS ) { + std::cerr << "Error: Phase::EXECUTE\n"; + return; + } + if( !equals_matrix( C, input.C_monoid ) ) { + std::cerr << "Error: Wrong result\n"; + rc = FAILED; + return; + } + + std::cout << "Result (monoid) is correct\n"; + } + + rc = SUCCESS; +} + +void test_program( const size_t& N, size_t& ) { + /** Matrix A: Matrix filled with A_INITIAL_VALUE + * X X X X X + * _ _ _ _ _ + * _ _ _ _ _ (...) + * _ _ _ _ _ + * _ _ _ _ _ + * (...) + */ + Matrix< nz_type > A( N, N, N ); + Matrix< void > A_void( N, N, N ); + { + std::vector< size_t > A_rows( N, 0 ), A_cols( N, 0 ); + std::vector< nz_type > A_values( N, A_INITIAL_VALUE ); + std::iota( A_cols.begin(), A_cols.end(), 0 ); + if( + SUCCESS != + buildMatrixUnique( A, A_rows.data(), A_cols.data(), A_values.data(), A_values.size(), SEQUENTIAL ) + ) { + throw std::runtime_error("(LINE " + std::to_string(__LINE__) + + ": Test FAILED: buildMatrixUnique" ); + } + if( + SUCCESS != + buildMatrixUnique( A_void, A_rows.data(), A_cols.data(), A_rows.size(), SEQUENTIAL ) + ) { + throw std::runtime_error("(LINE " + std::to_string(__LINE__) + + ": Test FAILED: buildMatrixUnique" ); + } + } + + + /** Matrix B: Column matrix filled with B_INITIAL_VALUE + * Y _ _ _ _ + * Y _ _ _ _ + * Y _ _ _ _ (...) + * Y _ _ _ _ + * Y _ _ _ _ + * (...) + */ + Matrix< nz_type > B( N, N, N ); + Matrix< void > B_void( N, N, N ); + { + std::vector< size_t > B_rows( N, 0 ), B_cols( N, 0 ); + std::vector< nz_type > B_values( N, B_INITIAL_VALUE ); + std::iota( B_rows.begin(), B_rows.end(), 0 ); + if( SUCCESS != + buildMatrixUnique( B, B_rows.data(), B_cols.data(), B_values.data(), B_values.size(), SEQUENTIAL) + ) { + throw std::runtime_error("(LINE " + std::to_string(__LINE__) + + ": Test FAILED: buildMatrixUnique" ); + } + if( + SUCCESS != + buildMatrixUnique( B_void, B_rows.data(), B_cols.data(), B_rows.size(), SEQUENTIAL ) + ) { + throw std::runtime_error("(LINE " + std::to_string(__LINE__) + + ": Test FAILED: buildMatrixUnique" ); + } + } + + { // C = A .+ B + std::cout << "-- Test C = A .+ B\n"; + /** Matrix C_monoid_truth: Union of A and B + * X+Y X X X X + * Y ___ ___ ___ ___ + * Y ___ ___ ___ ___ (...) + * Y ___ ___ ___ ___ + * Y ___ ___ ___ ___ + * (...) + */ + Matrix< nz_type > C_monoid_truth( N, N ); + size_t nvalues = nrows( A ) + ncols( B ) - 1; + std::vector< size_t > C_monoid_truth_rows( nvalues, 0 ), C_monoid_truth_cols( nvalues, 0 ); + std::vector< nz_type > C_monoid_truth_values( nvalues, 0 ); + C_monoid_truth_values[ 0 ] = A_INITIAL_VALUE + B_INITIAL_VALUE; + std::iota( C_monoid_truth_rows.begin() + nrows( A ), C_monoid_truth_rows.end(), 1 ); + std::iota( C_monoid_truth_cols.begin() + 1, C_monoid_truth_cols.begin() + nrows( A ), 1 ); + std::fill( C_monoid_truth_values.begin() + 1, C_monoid_truth_values.begin() + nrows( A ), A_INITIAL_VALUE ); + std::fill( C_monoid_truth_values.begin() + nrows( A ), C_monoid_truth_values.end(), B_INITIAL_VALUE ); + if( SUCCESS != + buildMatrixUnique( + C_monoid_truth, + C_monoid_truth_rows.data(), + C_monoid_truth_cols.data(), + C_monoid_truth_values.data(), + C_monoid_truth_values.size(), + SEQUENTIAL) + ) { + throw std::runtime_error("(LINE " + std::to_string(__LINE__) + + ": Test FAILED: buildMatrixUnique" ); + } + + /** Matrix C_op_truth: Intersection of A and B + * X+Y ___ ___ ___ ___ + * ___ ___ ___ ___ ___ + * ___ ___ ___ ___ ___ (...) + * ___ ___ ___ ___ ___ + * ___ ___ ___ ___ ___ + * (...) + */ + Matrix< nz_type > C_op_truth( N, N ); + std::vector< size_t > C_op_truth_rows( 1, 0 ), C_op_truth_cols( 1, 0 ); + std::vector< nz_type > C_op_truth_values( 1, A_INITIAL_VALUE + B_INITIAL_VALUE ); + if( SUCCESS != + buildMatrixUnique( + C_op_truth, + C_op_truth_rows.data(), + C_op_truth_cols.data(), + C_op_truth_values.data(), + C_op_truth_values.size(), + SEQUENTIAL) + ) { + throw std::runtime_error("(LINE " + std::to_string(__LINE__) + + ": Test FAILED: buildMatrixUnique" ); + } + + input_t< + Monoid< operators::add< nz_type >, identities::zero >, + nz_type, nz_type, nz_type + > input { A, B, C_monoid_truth, C_op_truth }; + output_t output { SUCCESS }; + // Run the test + grb_program(input, output ); + if( output.rc != SUCCESS ) { + throw std::runtime_error("(LINE " + std::to_string(__LINE__) + + ": Test FAILED (" + toString( output.rc ) + ")" ); + } + } + + { // C = A .+ A + std::cout << "-- Test C = A .+ A\n"; + /** Matrix C_truth: Union/intersection of A and A + * X+X X+X X+X X+X X+X + * ___ ___ ___ ___ ___ + * ___ ___ ___ ___ ___(...) + * ___ ___ ___ ___ ___ + * ___ ___ ___ ___ ___ + * (...) + */ + Matrix< nz_type > C_truth( N, N ); + size_t nvalues = ncols( A ); + std::vector< size_t > C_truth_rows( nvalues, 0 ), C_truth_cols( nvalues, 0 ); + std::vector< nz_type > C_truth_values( nvalues, A_INITIAL_VALUE+A_INITIAL_VALUE ); + std::iota( C_truth_cols.begin(), C_truth_cols.end(), 0 ); + if( SUCCESS != + buildMatrixUnique( + C_truth, + C_truth_rows.data(), + C_truth_cols.data(), + C_truth_values.data(), + C_truth_values.size(), + SEQUENTIAL + ) + ) { + throw std::runtime_error("(LINE " + std::to_string(__LINE__) + + ": Test FAILED: buildMatrixUnique" ); + } + + input_t< + Monoid< operators::add< nz_type >, identities::zero >, + nz_type, nz_type, nz_type + > input { A, A, C_truth, C_truth }; + output_t output { SUCCESS }; + // Run the test + grb_program(input, output ); + if( output.rc != SUCCESS ) { + throw std::runtime_error("(LINE " + std::to_string(__LINE__) + + ": Test FAILED (" + toString( output.rc ) + ")" ); + } + } + + { // C = A .+ A(void) + std::cout << "-- Test C = A .+ A(void)\n"; + /** Matrix C_truth: Union/intersection of A and A + * X+0 X+0 X+0 X+0 X+0 + * ___ ___ ___ ___ ___ + * ___ ___ ___ ___ ___(...) + * ___ ___ ___ ___ ___ + * ___ ___ ___ ___ ___ + * (...) + */ + const Matrix< nz_type >& C_truth = A; + + input_t< + Monoid< operators::add< nz_type >, identities::zero >, + nz_type, void, nz_type + > input { A, A_void, C_truth, C_truth }; + output_t output { SUCCESS }; + // Run the test + grb_program(input, output ); + if( output.rc != SUCCESS ) { + throw std::runtime_error("(LINE " + std::to_string(__LINE__) + + ": Test FAILED (" + toString( output.rc ) + ")" ); + } + } + + { // C = A(void) .+ A + std::cout << "-- Test C = A(void) .+ A\n"; + /** Matrix C_truth: Union/intersection of A and A + * 0+X 0+X 0+X 0+X 0+X + * ___ ___ ___ ___ ___ + * ___ ___ ___ ___ ___(...) + * ___ ___ ___ ___ ___ + * ___ ___ ___ ___ ___ + * (...) + */ + const Matrix< nz_type >& C_truth = A; + + input_t< + Monoid< operators::add< nz_type >, identities::zero >, + void, nz_type, nz_type + > input { A_void, A, C_truth, C_truth }; + output_t output { SUCCESS }; + // Run the test + grb_program(input, output ); + if( output.rc != SUCCESS ) { + throw std::runtime_error("(LINE " + std::to_string(__LINE__) + + ": Test FAILED (" + toString( output.rc ) + ")" ); + } + } + + { // C = A(void) .+ A + std::cout << "-- Test C = A(void) .+ A(void)\n"; + /** Matrix C_truth: Union/intersection of A and A + * 0+0 0+0 0+0 0+0 0+0 + * ___ ___ ___ ___ ___ + * ___ ___ ___ ___ ___(...) + * ___ ___ ___ ___ ___ + * ___ ___ ___ ___ ___ + * (...) + */ + Matrix< nz_type > C_truth( N, N ); + size_t nvalues = ncols( A ); + std::vector< size_t > C_truth_rows( nvalues, 0 ), C_truth_cols( nvalues, 0 ); + std::vector< nz_type > C_truth_values( nvalues, 0 ); + std::iota( C_truth_cols.begin(), C_truth_cols.end(), 0 ); + if( SUCCESS != + buildMatrixUnique( + C_truth, + C_truth_rows.data(), + C_truth_cols.data(), + C_truth_values.data(), + C_truth_values.size(), + SEQUENTIAL + ) + ) { + throw std::runtime_error("(LINE " + std::to_string(__LINE__) + + ": Test FAILED: buildMatrixUnique" ); + } + + input_t< + Monoid< operators::add< nz_type >, identities::zero >, + void, void, nz_type + > input { A_void, A_void, C_truth, C_truth }; + output_t output { SUCCESS }; + // Run the test + grb_program(input, output ); + if( output.rc != SUCCESS ) { + throw std::runtime_error("(LINE " + std::to_string(__LINE__) + + ": Test FAILED (" + toString( output.rc ) + ")" ); + } + } + + + { // C = A .+ Bt + std::cout << "-- Test C = A .+ Bt\n"; + /** Matrix C_truth: Union/intersection of A and Bt + * X+Y X+Y X+Y X+Y X+Y + * ___ ___ ___ ___ ___ + * ___ ___ ___ ___ ___(...) + * ___ ___ ___ ___ ___ + * ___ ___ ___ ___ ___ + * (...) + */ + Matrix< nz_type > C_truth( N, N ); + size_t nvalues = ncols( A ); + std::vector< size_t > C_truth_rows( nvalues, 0 ), C_truth_cols( nvalues, 0 ); + std::vector< nz_type > C_truth_values( nvalues, A_INITIAL_VALUE+B_INITIAL_VALUE ); + std::iota( C_truth_cols.begin(), C_truth_cols.end(), 0 ); + if( SUCCESS != + buildMatrixUnique( + C_truth, + C_truth_rows.data(), + C_truth_cols.data(), + C_truth_values.data(), + C_truth_values.size(), + SEQUENTIAL + ) + ) { + throw std::runtime_error("(LINE " + std::to_string(__LINE__) + + ": Test FAILED: buildMatrixUnique" ); + } + + input_t< + Monoid< operators::add< nz_type >, identities::zero >, + nz_type, nz_type, nz_type, + descriptors::transpose_right + > input { A, B, C_truth, C_truth }; + output_t output { SUCCESS }; + // Run the test + grb_program(input, output ); + if( output.rc != SUCCESS ) { + throw std::runtime_error("(LINE " + std::to_string(__LINE__) + + ": Test FAILED (" + toString( output.rc ) + ")" ); + } + } + + +} + +int main( int argc, char ** argv ) { + (void) argc; + (void) argv; + + size_t N = 1000; + + if( argc > 2 ) { + std::cout << "Usage: " << argv[ 0 ] << " [n=" << N << "]" << std::endl; + return 1; + } + if( argc == 2 ) { + N = std::stoul( argv[ 1 ] ); + } + + std::cout << "This is functional test " << argv[ 0 ] << std::endl << std::flush; + + // Launch the test + Launcher< AUTOMATIC > launcher; + RC rc = launcher.exec( &test_program, N, N, true ); + if( rc != SUCCESS ) { + std::cout << "Test FAILED (" << grb::toString( rc ) << ")" << std::endl; + return static_cast( rc ); + } + + std::cerr << std::flush; + std::cout << std::flush << "Test OK" << std::endl; + return 0; +} diff --git a/tests/unit/eWiseApply_matrix.cpp b/tests/unit/eWiseApply_matrix.cpp index 5f24b0365..01f92e9bf 100644 --- a/tests/unit/eWiseApply_matrix.cpp +++ b/tests/unit/eWiseApply_matrix.cpp @@ -94,4 +94,3 @@ int main( int argc, char ** argv ) { return 0; } } - diff --git a/tests/unit/spy.cpp b/tests/unit/spy.cpp index 03dac8d17..f773bae82 100644 --- a/tests/unit/spy.cpp +++ b/tests/unit/spy.cpp @@ -82,7 +82,8 @@ void grb_program( const void * const fn_p, const size_t fn_length, grb::RC &rc ) if( rc == grb::SUCCESS ) { grb::Matrix< double > chk( p, q ); rc = rc ? rc : grb::resize( chk, grb::nnz( spy ) ); - rc = rc ? rc : grb::eWiseApply( chk, spy, spy2, ring.getMultiplicativeOperator() ); + rc = rc ? rc : grb::eWiseApply( chk, spy, spy2, ring.getMultiplicativeMonoid(), grb::Phase::RESIZE ); + rc = rc ? rc : grb::eWiseApply( chk, spy, spy2, ring.getMultiplicativeMonoid() ); if( rc == grb::SUCCESS && grb::nnz( chk ) != grb::nnz( spy ) ) { std::cerr << "Unexpected number of nonzeroes for chk: " << grb::nnz(chk) << ", expected " << grb::nnz(spy) << "\n"; rc = grb::FAILED; @@ -114,7 +115,8 @@ void grb_program( const void * const fn_p, const size_t fn_length, grb::RC &rc ) if( rc == grb::SUCCESS ) { grb::Matrix< double > chk( p, q ); rc = rc ? rc : grb::resize( chk, nnz( spy ) ); - rc = rc ? rc : grb::eWiseApply( chk, spy, spy2, ring.getMultiplicativeOperator() ); + rc = rc ? rc : grb::eWiseApply( chk, spy, spy2, ring.getMultiplicativeMonoid(), grb::Phase::RESIZE ); + rc = rc ? rc : grb::eWiseApply( chk, spy, spy2, ring.getMultiplicativeMonoid() ); if( rc == grb::SUCCESS && grb::nnz( chk ) != grb::nnz( spy ) ) { std::cerr << "Unexpected number of nonzeroes for chk (pattern): " << grb::nnz(chk) << ", expected " << grb::nnz(spy) << "\n"; rc = grb::FAILED; @@ -146,7 +148,8 @@ void grb_program( const void * const fn_p, const size_t fn_length, grb::RC &rc ) if( rc == grb::SUCCESS ) { grb::Matrix< double > chk( p, q ); rc = rc ? rc : grb::resize( chk, nnz( spy ) ); - rc = rc ? rc : grb::eWiseApply( chk, spy, spy2, ring.getMultiplicativeOperator() ); + rc = rc ? rc : grb::eWiseApply( chk, spy, spy2, ring.getMultiplicativeMonoid(), grb::Phase::RESIZE ); + rc = rc ? rc : grb::eWiseApply( chk, spy, spy2, ring.getMultiplicativeMonoid() ); if( rc == grb::SUCCESS && grb::nnz( chk ) != grb::nnz( spy ) ) { std::cerr << "Unexpected number of nonzeroes for chk (boolean): " << grb::nnz(chk) << ", expected " << grb::nnz(spy) << "\n"; rc = grb::FAILED; diff --git a/tests/unit/unittests.sh b/tests/unit/unittests.sh index 5b7dfc16a..fdb4dc6c2 100755 --- a/tests/unit/unittests.sh +++ b/tests/unit/unittests.sh @@ -597,12 +597,26 @@ for MODE in ${MODES}; do grep 'Test OK' ${TEST_OUT_DIR}/eWiseApply_matrix_${MODE}_${BACKEND}_${P}_${T} || echo "Test FAILED" echo " " - echo ">>> [x] [ ] Testing grb::eWiseLambda (matrices)" - $runner ${TEST_BIN_DIR}/eWiseMatrix_${MODE}_${BACKEND} &> ${TEST_OUT_DIR}/eWiseMatrix_${MODE}_${BACKEND}_${P}_${T}.log - head -1 ${TEST_OUT_DIR}/eWiseMatrix_${MODE}_${BACKEND}_${P}_${T}.log - grep 'Test OK' ${TEST_OUT_DIR}/eWiseMatrix_${MODE}_${BACKEND}_${P}_${T}.log || echo "Test FAILED" + echo ">>> [x] [ ] Testing grb::id on vectors and matrices" + $runner ${TEST_BIN_DIR}/id_${MODE}_${BACKEND} &> ${TEST_OUT_DIR}/id_${MODE}_${BACKEND}_${P}_${T}.log + head -1 ${TEST_OUT_DIR}/id_${MODE}_${BACKEND}_${P}_${T}.log + grep 'Test OK' ${TEST_OUT_DIR}/id_${MODE}_${BACKEND}_${P}_${T}.log || echo "Test FAILED" echo " " + echo ">>> [x] [ ] Testing grb::eWiseApply (matrices, Monoid / Operator)" + echo " using small matrices (100x100)" + $runner ${TEST_BIN_DIR}/eWiseApplyMatrix_variants_${MODE}_${BACKEND} 100 &> ${TEST_OUT_DIR}/eWiseApplyMatrix_variants_small_${MODE}_${BACKEND}_${P}_${T}.log + head -1 ${TEST_OUT_DIR}/eWiseApplyMatrix_variants_small_${MODE}_${BACKEND}_${P}_${T}.log + grep 'Test OK' ${TEST_OUT_DIR}/eWiseApplyMatrix_variants_small_${MODE}_${BACKEND}_${P}_${T}.log || echo "Test FAILED" + echo " " + + echo ">>> [x] [ ] Testing grb::eWiseApply (matrices, Monoid / Operator)" + echo " using large matrices (10'000x10'000)" + $runner ${TEST_BIN_DIR}/eWiseApplyMatrix_variants_${MODE}_${BACKEND} 10000 &> ${TEST_OUT_DIR}/eWiseApplyMatrix_variants_large_${MODE}_${BACKEND}_${P}_${T}.log + head -1 ${TEST_OUT_DIR}/eWiseApplyMatrix_variants_large_${MODE}_${BACKEND}_${P}_${T}.log + grep 'Test OK' ${TEST_OUT_DIR}/eWiseApplyMatrix_variants_large_${MODE}_${BACKEND}_${P}_${T}.log || echo "Test FAILED" + echo " " + echo ">>> [x] [ ] Testing grb::zip on two vectors of doubles and" echo " ints of size 10 000 000." $runner ${TEST_BIN_DIR}/zip_${MODE}_${BACKEND} 10000000 &> ${TEST_OUT_DIR}/zip_large_${MODE}_${BACKEND}_${P}_${T}