diff --git a/include/CMakeLists.txt b/include/CMakeLists.txt index cb6d5ba97..2511528ee 100644 --- a/include/CMakeLists.txt +++ b/include/CMakeLists.txt @@ -52,7 +52,7 @@ set( root_files "graphblas/rc.hpp" "graphblas/semiring.hpp" "graphblas/spmd.hpp" "graphblas/tags.hpp" "graphblas/type_traits.hpp" "graphblas/utils.hpp" "graphblas/vector.hpp" "graphblas/synchronizedNonzeroIterator.hpp" - "graphblas/nonzeroStorage.hpp" + "graphblas/nonzeroStorage.hpp" "graphblas/selection_ops.hpp" ) set( GRB_INCLUDE_INSTALL_DIR "${INCLUDE_INSTALL_DIR}/graphblas") diff --git a/include/graphblas.hpp b/include/graphblas.hpp index c1374327e..029f77bf7 100644 --- a/include/graphblas.hpp +++ b/include/graphblas.hpp @@ -513,6 +513,7 @@ namespace grb { #include #include #include +#include // Then include containers. If containers rely on ALP/GraphBLAS primitives that // are defined as free functions, then container implementations must forward- diff --git a/include/graphblas/base/blas3.hpp b/include/graphblas/base/blas3.hpp index 7c91b871f..4eb40e334 100644 --- a/include/graphblas/base/blas3.hpp +++ b/include/graphblas/base/blas3.hpp @@ -442,6 +442,133 @@ namespace grb { return ret == SUCCESS ? UNSUPPORTED : ret; } + /** + * Selects elements from a matrix based on a given selection boolean operator. + * + * This function template is used to select elements from a given input matrix + * based on a provided selection operator. The selected elements are then stored + * in a different output matrix, making this an out-of-place operation. The + * output matrix is cleared before the selection. + * + * After a successful call to this primitive, the nonzero structure of \a B + * will match the one of \a A without the elements that were not matched by + * the selection operator. Any values at those positions are copied from \a A + * to \a B. All the elements of \a B will normally return true when + * applied to the selection operator. + * + * \note An exception to the last point may occur if the value types of \a A + * and \a B do not match, while the selection operator depends on those + * values in a way that makes it behave differently. + * + * @tparam descr The descriptor to be used. Optional; the default + * is #grb::descriptors::no_operation. + * @tparam SelectionOperator The selection operator type, a function with the + * following signature: + * `bool( const RIT &, const CIT &, const T & )`. + * Here, + * - RIT: The row index type of the input matrix, + * or a type that is convertible from it. + * - CIT: The column index type of the input matrix, + * or a type that is convertible from it. + * - T: The value type of the input matrix, or a + * type that is convertible to it. + * + * The types for \a RIT and \a CIT are given by + * -# grb::config::RowIndexType and + * -# grb::config::ColIndexType, + * respectively. For most use cases, the default is unsigned int for + * both types. The safest and most performant choice is therefore to supply an + * operator with the aforementioned two configuration types for \a RIT and + * \a CIT. The most generic safe choice that does not depend on configured + * types is size_t, but such use may result in a (slight) performance + * penalty due to internal casting between possibly different index types. + * + * For void matrices, the select operator will assume a bool + * for \a T. The operator will always receive true as the value + * corresponding to the sparse pattern. + * + * @tparam Tin The value type of the input matrix. + * @tparam RITin The row index type of the input matrix. + * @tparam CITin The column index type of the input matrix. + * @tparam NITin The nonzero index type of the input matrix. + * @tparam Tout The value type of the output matrix. + * @tparam RITout The row index type of the output matrix. + * @tparam CITout The column index type of the output matrix. + * @tparam NITout The nonzero index type of the output matrix. + * @tparam backend The backend to use for the operation. + * + * @param[out] B The output matrix. Will be cleared before the selection. + * @param[in] A The input matrix. + * @param[in] op The selection boolean operator. + * @param[in] phase The #grb::Phase the call should execute. Optional; the + * default parameter is #grb::EXECUTE. + * + * \note Pre-defined selection operators can be found in the namespace + * #grb::operators::select. + * + * @return #grb::SUCCESS On successful completion of this call. + * @return #grb::MISMATCH Whenever the dimensions of \a A and \a B 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 B was insufficient. The output + * matrix \a B 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 B 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. + * + * \parblock + * \par Descriptors + * + * Only #grb::descriptors::no_casting is accepted. + * \endparblock + * + * \par Performance semantics + * + * Each backend must define performance semantics for this primitive. + * + * @see perfSemantics + */ + template< + Descriptor descr = descriptors::no_operation, + class SelectionOperator, + typename Tin, + typename RITin, typename CITin, typename NITin, + typename Tout, + typename RITout, typename CITout, typename NITout, + Backend backend + > + RC select( + Matrix< Tout, backend, RITout, CITout, NITout > &B, + const Matrix< Tin, backend, RITin, CITin, NITin > &A, + const SelectionOperator &op, + const Phase &phase = EXECUTE, + const typename std::enable_if< + !is_object< Tin >::value && + !is_object< Tout >::value + >::type * const = nullptr + ) { + (void) descr; + (void) B; + (void) A; + (void) op; + (void) phase; +#ifdef _DEBUG + std::cerr << "Selected backend does not implement grb::select\n"; +#endif +#ifndef NDEBUG + const bool selected_backend_does_not_support_select = false; + assert( selected_backend_does_not_support_select ); +#endif + return UNSUPPORTED; + } + /** * Reduces, or \em folds, a matrix into a scalar, according to a commutative * monoid. diff --git a/include/graphblas/base/internalops.hpp b/include/graphblas/base/internalops.hpp index 9a5e8ca7d..fb4bf30b6 100644 --- a/include/graphblas/base/internalops.hpp +++ b/include/graphblas/base/internalops.hpp @@ -4276,6 +4276,163 @@ namespace grb { } // namespace operators + namespace operators::select::internal { + + /** + * This class takes a generic operator implementation and exposes a more + * convenient operator() function based on it. This function allows arbitrary + * data types being passed as parameters, and automatically handles any + * casting required for the raw operator. + * + * @tparam OP The generic operator implementation. + * + * @see Operator for full details. + */ + template< typename OP, typename > + struct MatrixSelectionOperatorBase { + typedef typename OP::value_type D; + typedef typename OP::row_type RIT; + typedef typename OP::column_type CIT; + + template< typename RIT1, typename CIT1, typename D1 > + bool operator()( + const RIT1 &x, const CIT1 &y, const D1 &v + ) const noexcept { + const RIT a = static_cast< RIT >( x ); + const CIT b = static_cast< CIT >( y ); + const D c = static_cast< D >( v ); + return OP::apply( &a, &b, &c ); + } + + /** + * This is the high-performance version of apply() in the sense that no + * casting is required. This version will be automatically called whenever + * possible (non-void variant). + */ + bool operator()( + const RIT &x, const CIT &y, const D &v + ) const noexcept { + return OP::apply( &x, &y, &v ); + } + + }; + + /** This is the void value type variant. */ + template< typename OP > + struct MatrixSelectionOperatorBase< OP, void > { + typedef typename OP::row_type RIT; + typedef typename OP::column_type CIT; + + template< typename RIT1, typename CIT1, typename D1 > + bool operator()( + const RIT1 &x, const CIT1 &y, const D1 &v + ) const noexcept { + (void) v; + const RIT a = static_cast< RIT >( x ); + const CIT b = static_cast< CIT >( y ); + return OP::apply( &a, &b, nullptr ); + } + + }; + + /** + * Implements the is_diagonal matrix selector. + */ + template< + typename D, typename RIT, typename CIT + > + struct is_diagonal { + typedef D value_type; + typedef RIT row_type; + typedef CIT column_type; + + static bool apply( + const row_type * __restrict__ const x, + const column_type * __restrict__ const y, + const value_type * __restrict__ const + ) { + return *x == *y; + } + }; + + /** + * Implements the strictly lower triangular matrix selector. + */ + template< + typename D, typename RIT, typename CIT + > + struct is_strictly_lower { + typedef D value_type; + typedef RIT row_type; + typedef CIT column_type; + + static bool apply( + const row_type * __restrict__ const x, + const column_type * __restrict__ const y, + const value_type * __restrict__ const + ) { + return *x > *y; + } + }; + + /** + * Implements the lower triangular matrix selector. + */ + template< + typename D, typename RIT, typename CIT + > + struct is_lower_or_diagonal { + typedef D value_type; + typedef RIT row_type; + typedef CIT column_type; + + static bool apply( + const row_type * __restrict__ const x, + const column_type * __restrict__ const y, + const value_type * __restrict__ const + ) { + return *x >= *y; + } + }; + + /** + * Implements the strictly upper triangular matrix selector. + */ + template< + typename D, typename RIT, typename CIT + > + struct is_strictly_upper { + typedef D value_type; + typedef RIT row_type; + typedef CIT column_type; + + static bool apply( + const row_type * __restrict__ const x, + const column_type * __restrict__ const y, + const value_type * __restrict__ const v + ) { return !is_lower_or_diagonal< D, RIT, CIT >::apply( x, y, v ); } + }; + + /** + * Implements the upper triangular matrix selector. + */ + template< + typename D, typename RIT, typename CIT + > + struct is_upper_or_diagonal { + typedef D value_type; + typedef RIT row_type; + typedef CIT column_type; + + static bool apply( + const row_type * __restrict__ const x, + const column_type * __restrict__ const y, + const value_type * __restrict__ const v + ) { return !is_strictly_lower< D, RIT, CIT >::apply( x, y, v ); } + }; + + } // namespace operators::select::internal + } // namespace grb #endif // _H_GRB_INTERNAL_OPERATORS_BASE diff --git a/include/graphblas/bsp1d/blas3.hpp b/include/graphblas/bsp1d/blas3.hpp index f3526480c..6db61258c 100644 --- a/include/graphblas/bsp1d/blas3.hpp +++ b/include/graphblas/bsp1d/blas3.hpp @@ -200,13 +200,57 @@ namespace grb { if( collectives<>::allreduce( ret, operators::any_or< RC >() ) != SUCCESS ) { return PANIC; } else { - return SUCCESS; + return ret; } } assert( phase == EXECUTE ); return internal::checkGlobalErrorStateOrClear( C, ret ); } + /** \internal Simply delegates to process-local backend */ + template< + Descriptor descr = descriptors::no_operation, + class SelectionOperator, + typename Tin, + typename RITin, typename CITin, typename NITin, + typename Tout, + typename RITout, typename CITout, typename NITout + > + RC select( + Matrix< Tout, BSP1D, RITout, CITout, NITout > &out, + const Matrix< Tin, BSP1D, RITin, CITin, NITin > &in, + const SelectionOperator &op, + const Phase &phase = EXECUTE, + const typename std::enable_if< + !is_object< Tin >::value && + !is_object< Tout >::value + >::type * const = nullptr + ) { + assert( phase != TRY ); + + const auto coordinatesTranslationFunctions = + in.getLocalToGlobalCoordinatesTranslationFunctions(); + + RC ret = internal::select_generic< descr >( + internal::getLocal( out ), + internal::getLocal( in ), + op, + std::get<0>(coordinatesTranslationFunctions), + std::get<1>(coordinatesTranslationFunctions), + phase + ); + + if( phase == RESIZE ) { + if( collectives<>::allreduce( ret, operators::any_or< RC >() ) != SUCCESS ) { + return PANIC; + } else { + return ret; + } + } + assert( phase == EXECUTE ); + return internal::checkGlobalErrorStateOrClear( out, ret ); + } + template< Descriptor descr = descriptors::no_operation, class Monoid, @@ -923,3 +967,4 @@ namespace grb { } // namespace grb #endif + diff --git a/include/graphblas/bsp1d/matrix.hpp b/include/graphblas/bsp1d/matrix.hpp index f773939f0..feccd7061 100644 --- a/include/graphblas/bsp1d/matrix.hpp +++ b/include/graphblas/bsp1d/matrix.hpp @@ -150,6 +150,29 @@ namespace grb { const Ring &, const Phase & ); + /* ********************* + BLAS-2 friends + ********************* */ + + template< + Descriptor, + class SelectionOperator, + typename Tin, + typename RITin, typename CITin, typename NITin, + typename Tout, + typename RITout, typename CITout, typename NITout + > + friend RC select( + Matrix< Tout, BSP1D, RITout, CITout, NITout > &, + const Matrix< Tin, BSP1D, RITin, CITin, NITin > &, + const SelectionOperator &, + const Phase &, + const typename std::enable_if< + !is_object< Tin >::value && + !is_object< Tout >::value + >::type * const + ); + /* ********************* Internal friends ********************* */ @@ -334,6 +357,42 @@ namespace grb { } + protected: + + /** + * Helper functions to get the global coordinates of this matrix + * from local coordinates. + * + * \return [ + * 0: local row index to global row index, + * 1: local column index to global column index, + * ] + */ + std::tuple< + std::function< size_t( size_t ) >, + std::function< size_t( size_t ) > + > getLocalToGlobalCoordinatesTranslationFunctions() const noexcept { + const auto &lpf_data = internal::grb_BSP1D.cload(); + const size_t rows = nrows( *this ); + const size_t columns = ncols( *this ); + + return std::make_tuple( + [ lpf_data, rows ]( const size_t i ) -> size_t { + return internal::Distribution< BSP1D >::local_index_to_global( + i, rows, lpf_data.s, lpf_data.P ); + }, + [ lpf_data, columns ]( const size_t j ) -> size_t { + const size_t col_pid = internal::Distribution<>::offset_to_pid( + j, columns, lpf_data.P ); + const size_t col_off = internal::Distribution<>::local_offset( + columns, col_pid, lpf_data.P ); + return internal::Distribution< BSP1D >::local_index_to_global( + j - col_off, columns, col_pid, lpf_data.P ); + } + ); + } + + public: /** @see Matrix::value_type */ diff --git a/include/graphblas/hyperdags/blas3.hpp b/include/graphblas/hyperdags/blas3.hpp index e97995204..44f9daa86 100644 --- a/include/graphblas/hyperdags/blas3.hpp +++ b/include/graphblas/hyperdags/blas3.hpp @@ -334,6 +334,50 @@ namespace grb { return ret; } + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename Tin, + typename RITin, typename CITin, typename NITin, + typename Tout, + typename RITout, typename CITout, typename NITout + > + RC select( + Matrix< Tout, hyperdags, RITout, CITout, NITout > &out, + const Matrix< Tin, hyperdags, RITin, CITin, NITin > &in, + const Operator &op = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + !is_object< Tin >::value && + !is_object< Tout >::value + >::type * const = nullptr + ) { + const RC ret = select< descr >( + internal::getMatrix( out ), + internal::getMatrix( in ), + op, + phase + ); + if( ret != SUCCESS ) { return ret; } + if( phase != EXECUTE ) { return ret; } + if( nrows( out ) == 0 || ncols( out ) == 0 ) { return ret; } + std::array< const void *, 0 > sourcesP{}; + std::array< uintptr_t, 2 > sourcesC{ + getID( internal::getMatrix( in ) ), + getID( internal::getMatrix( out ) ) + }; + std::array< uintptr_t, 1 > destinations{ + getID( internal::getMatrix( out ) ) + }; + internal::hyperdags::generator.addOperation( + internal::hyperdags::SELECT_MATRIX_MATRIX, + sourcesP.begin(), sourcesP.end(), + sourcesC.begin(), sourcesC.end(), + destinations.begin(), destinations.end() + ); + return ret; + } + template< Descriptor descr = descriptors::no_operation, class MonoidOrSemiring, diff --git a/include/graphblas/hyperdags/hyperdags.hpp b/include/graphblas/hyperdags/hyperdags.hpp index d032c8e38..5ecc4dbef 100644 --- a/include/graphblas/hyperdags/hyperdags.hpp +++ b/include/graphblas/hyperdags/hyperdags.hpp @@ -490,17 +490,19 @@ namespace grb { EWISELAMBDA_FUNC_VECTOR, + SELECT_MATRIX_MATRIX, + FOLDL_SCALAR_MATRIX_MASK_MONOID, FOLDL_SCALAR_MATRIX_MONOID, FOLDR_SCALAR_MATRIX_MASK_MONOID, - FOLDR_SCALAR_MATRIX_MONOID, + FOLDR_SCALAR_MATRIX_MONOID }; /** \internal How many operation vertex types exist. */ - const constexpr size_t numOperationVertexTypes = 110; + const constexpr size_t numOperationVertexTypes = 111; /** \internal An array of all operation vertex types. */ const constexpr enum OperationVertexType @@ -612,10 +614,11 @@ namespace grb { EWISEMUL_VECTOR_VECTOR_VECTOR_BETA_RING, EWISEMUL_VECTOR_VECTOR_ALPHA_BETA_RING, EWISELAMBDA_FUNC_VECTOR, + SELECT_MATRIX_MATRIX, FOLDL_SCALAR_MATRIX_MASK_MONOID, FOLDL_SCALAR_MATRIX_MONOID, FOLDR_SCALAR_MATRIX_MASK_MONOID, - FOLDR_SCALAR_MATRIX_MONOID, + FOLDR_SCALAR_MATRIX_MONOID }; /** \internal @returns The operation vertex type as a string. */ diff --git a/include/graphblas/nonblocking/blas3.hpp b/include/graphblas/nonblocking/blas3.hpp index dacbbdc77..5f7e2c63f 100644 --- a/include/graphblas/nonblocking/blas3.hpp +++ b/include/graphblas/nonblocking/blas3.hpp @@ -571,6 +571,49 @@ namespace grb { ); } + template< + Descriptor descr = descriptors::no_operation, + class Operator, + typename Tin, + typename RITin, typename CITin, typename NITin, + typename Tout, + typename RITout, typename CITout, typename NITout + > + RC select( + Matrix< Tout, nonblocking, RITout, CITout, NITout > &out, + const Matrix< Tin, nonblocking, RITin, CITin, NITin > &in, + const Operator &op = Operator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + !is_object< Tin >::value && + !is_object< Tout >::value + >::type * const = nullptr + ) { +#ifdef _DEBUG + std::cout << "In grb::select( nonblocking )\n"; +#endif + if( internal::NONBLOCKING::warn_if_not_native && + config::PIPELINE::warn_if_not_native + ) { + std::cerr << "Warning: select (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 select< descr, Operator >( + internal::getRefMatrix( out ), + internal::getRefMatrix( in ), + op, + phase + ); + } + template< Descriptor descr = descriptors::no_operation, class Monoid, diff --git a/include/graphblas/ops.hpp b/include/graphblas/ops.hpp index e5ff732d5..81b480929 100644 --- a/include/graphblas/ops.hpp +++ b/include/graphblas/ops.hpp @@ -53,7 +53,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class left_assign : public internal::Operator< @@ -62,7 +62,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = left_assign< A, B, C, D >; left_assign() {} @@ -80,7 +80,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class left_assign_if : public internal::Operator< @@ -89,7 +89,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = left_assign_if< A, B, C, D >; left_assign_if() {} @@ -110,7 +110,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class right_assign : public internal::Operator< internal::right_assign< D1, D2, D3, implementation > @@ -118,7 +118,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = right_assign< A, B, C, D >; right_assign() {} @@ -136,7 +136,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class right_assign_if : public internal::Operator< internal::right_assign_if< D1, D2, D3, implementation > @@ -144,7 +144,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = right_assign_if< A, B, C, D >; right_assign_if() {} @@ -170,7 +170,7 @@ namespace grb { // [Operator Wrapping] template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class add : public internal::Operator< internal::add< D1, D2, D3, implementation > @@ -178,7 +178,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = add< A, B, C, D >; add() {} @@ -203,7 +203,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class mul : public internal::Operator< internal::mul< D1, D2, D3, implementation > @@ -211,7 +211,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = mul< A, B, C, D >; mul() {} @@ -236,7 +236,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class max : public internal::Operator< internal::max< D1, D2, D3, implementation > @@ -244,7 +244,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = max< A, B, C, D >; max() {} @@ -269,7 +269,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class min : public internal::Operator< internal::min< D1, D2, D3, implementation > @@ -277,7 +277,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = min< A, B, C, D >; min() {} @@ -296,7 +296,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class subtract : public internal::Operator< internal::substract< D1, D2, D3, implementation > @@ -304,7 +304,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = subtract< A, B, C, D >; subtract() {} @@ -323,7 +323,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class divide : public internal::Operator< internal::divide< D1, D2, D3, implementation > @@ -331,7 +331,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = divide< A, B, C, D >; divide() {} @@ -348,7 +348,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class divide_reverse : public internal::Operator< internal::divide_reverse< D1, D2, D3, implementation > @@ -356,7 +356,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = divide_reverse< A, B, C, D >; divide_reverse() {} @@ -374,7 +374,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class equal : public internal::Operator< internal::equal< D1, D2, D3, implementation > @@ -382,7 +382,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = equal< A, B, C, D >; equal() {} @@ -400,7 +400,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class not_equal : public internal::Operator< internal::not_equal< D1, D2, D3, implementation > @@ -408,7 +408,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = not_equal< A, B, C, D >; not_equal() {} @@ -429,7 +429,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class any_or : public internal::Operator< internal::any_or< D1, D2, D3, implementation > @@ -437,7 +437,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = any_or< A, B, C, D >; any_or() {} @@ -457,7 +457,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class logical_or : public internal::Operator< internal::logical_or< D1, D2, D3, implementation > @@ -465,7 +465,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = logical_or< A, B, C, D >; logical_or() {} @@ -485,7 +485,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class logical_and : public internal::Operator< internal::logical_and< D1, D2, D3, implementation > @@ -493,7 +493,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = logical_and< A, B, C, D >; logical_and() {} @@ -509,7 +509,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class relu : public internal::Operator< internal::relu< D1, D2, D3, implementation > @@ -517,7 +517,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = relu< A, B, C, D >; relu() {} @@ -536,7 +536,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class abs_diff : public internal::Operator< internal::abs_diff< D1, D2, D3, implementation > @@ -544,7 +544,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = abs_diff< A, B, C, D >; abs_diff() {} @@ -618,7 +618,7 @@ namespace grb { */ template< typename D1, typename D2, typename D3, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class square_diff : public internal::Operator< internal::square_diff< D1, D2, D3, implementation > @@ -626,7 +626,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = square_diff< A, B, C, D >; square_diff() {} @@ -644,7 +644,7 @@ namespace grb { */ template< typename IN1, typename IN2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class zip : public internal::Operator< internal::zip< IN1, IN2, implementation > @@ -652,7 +652,7 @@ namespace grb { public: - template< typename A, typename B, enum Backend D > + template< typename A, typename B, Backend D > using GenericOperator = zip< A, B, D >; zip() {} @@ -673,7 +673,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class equal_first : public internal::Operator< internal::equal_first< D1, D2, D3, implementation > @@ -681,7 +681,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = equal_first< A, B, C, D >; equal_first() {} @@ -702,7 +702,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class less_than : public internal::Operator< internal::lt< D1, D2, D3, implementation > @@ -710,7 +710,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = less_than< A, B, C, D >; less_than() {} @@ -731,7 +731,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class leq : public internal::Operator< internal::leq< D1, D2, D3, implementation > @@ -739,7 +739,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = leq< A, B, C, D >; leq() {} @@ -760,7 +760,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class greater_than: public internal::Operator< internal::gt< D1, D2, D3, implementation > @@ -768,7 +768,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = greater_than< A, B, C, D >; greater_than() {} @@ -789,7 +789,7 @@ namespace grb { */ template< typename D1, typename D2 = D1, typename D3 = D2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class geq : public internal::Operator< internal::geq< D1, D2, D3, implementation > @@ -797,7 +797,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = geq< A, B, C, D >; geq() {} @@ -862,7 +862,7 @@ namespace grb { */ template< typename IN1, typename IN2, typename OUT, bool conj_left, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class conjugate_mul : public operators::internal::Operator< internal::conjugate_mul< IN1, IN2, OUT, conj_left, implementation > @@ -870,7 +870,7 @@ namespace grb { public: - template< typename A, typename B, typename C, bool D, enum Backend E > + template< typename A, typename B, typename C, bool D, Backend E > using GenericOperator = conjugate_mul< A, B, C, D, E >; conjugate_mul() {} @@ -913,7 +913,7 @@ namespace grb { */ template< typename IN1, typename IN2 = IN1, typename OUT = IN2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class conjugate_right_mul : public operators::internal::Operator< internal::conjugate_mul< IN1, IN2, OUT, false, implementation > @@ -921,7 +921,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = conjugate_right_mul< A, B, C, D >; conjugate_right_mul() {} @@ -964,7 +964,7 @@ namespace grb { */ template< typename IN1, typename IN2 = IN1, typename OUT = IN2, - enum Backend implementation = config::default_backend + Backend implementation = config::default_backend > class conjugate_left_mul : public operators::internal::Operator< internal::conjugate_mul< IN1, IN2, OUT, true, implementation > @@ -972,7 +972,7 @@ namespace grb { public: - template< typename A, typename B, typename C, enum Backend D > + template< typename A, typename B, typename C, Backend D > using GenericOperator = conjugate_left_mul< A, B, C, D >; conjugate_left_mul() {} @@ -981,215 +981,215 @@ namespace grb { } // namespace operators - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::left_assign_if< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::right_assign_if< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::left_assign< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::right_assign< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; // [Operator Type Traits] - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::add< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; // [Operator Type Traits] - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::mul< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::max< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::min< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::subtract< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::divide< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::divide_reverse< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::equal< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::not_equal< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::any_or< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::logical_or< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::logical_and< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::abs_diff< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::relu< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; template< typename IType, typename VType > struct is_operator< operators::argmin< IType, VType > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; template< typename IType, typename VType > struct is_operator< operators::argmax< IType, VType > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::square_diff< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename IN1, typename IN2, enum Backend implementation > + template< typename IN1, typename IN2, Backend implementation > struct is_operator< operators::zip< IN1, IN2, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::equal_first< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::less_than< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::leq< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::greater_than< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::geq< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; template< typename D1, typename D2, typename D3, - bool cl, enum Backend implementation + bool cl, Backend implementation > struct is_operator< operators::conjugate_mul< D1, D2, D3, cl, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::conjugate_left_mul< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct is_operator< operators::conjugate_right_mul< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; template< typename D1, typename D2, typename D3 > struct is_idempotent< operators::min< D1, D2, D3 >, void > { - static const constexpr bool value = true; + static constexpr bool value = true; }; template< typename D1, typename D2, typename D3 > struct is_idempotent< operators::max< D1, D2, D3 >, void > { - static const constexpr bool value = true; + static constexpr bool value = true; }; template< typename D1, typename D2, typename D3 > struct is_idempotent< operators::any_or< D1, D2, D3 >, void > { - static const constexpr bool value = true; + static constexpr bool value = true; }; template< typename D1, typename D2, typename D3 > struct is_idempotent< operators::logical_or< D1, D2, D3 >, void > { - static const constexpr bool value = true; + static constexpr bool value = true; }; template< typename D1, typename D2, typename D3 > struct is_idempotent< operators::logical_and< D1, D2, D3 >, void > { - static const constexpr bool value = true; + static constexpr bool value = true; }; template< typename D1, typename D2, typename D3 > struct is_idempotent< operators::relu< D1, D2, D3 >, void > { - static const constexpr bool value = true; + static constexpr bool value = true; }; template< typename D1, typename D2, typename D3 > struct is_idempotent< operators::left_assign_if< D1, D2, D3 >, void > { - static const constexpr bool value = true; + static constexpr bool value = true; }; template< typename D1, typename D2, typename D3 > struct is_idempotent< operators::right_assign_if< D1, D2, D3 >, void > { - static const constexpr bool value = true; + static constexpr bool value = true; }; template< typename IType, typename VType > struct is_idempotent< operators::argmin< IType, VType >, void > { - static const constexpr bool value = true; + static constexpr bool value = true; }; template< typename IType, typename VType > struct is_idempotent< operators::argmax< IType, VType >, void > { - static const constexpr bool value = true; + static constexpr bool value = true; }; template< typename OP > @@ -1197,7 +1197,7 @@ namespace grb { OP, typename std::enable_if< is_operator< OP >::value, void >::type > { - static constexpr const bool value = OP::is_associative(); + static constexpr bool value = OP::is_associative(); }; template< typename OP > @@ -1205,21 +1205,21 @@ namespace grb { OP, typename std::enable_if< is_operator< OP >::value, void >::type > { - static constexpr const bool value = OP::is_commutative(); + static constexpr bool value = OP::is_commutative(); }; // internal type traits follow namespace internal { - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct maybe_noop< operators::left_assign_if< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; - template< typename D1, typename D2, typename D3, enum Backend implementation > + template< typename D1, typename D2, typename D3, Backend implementation > struct maybe_noop< operators::right_assign_if< D1, D2, D3, implementation > > { - static const constexpr bool value = true; + static constexpr bool value = true; }; } // end namespace grb::internal diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp index 171bd32ee..744ab24e0 100644 --- a/include/graphblas/reference/blas3.hpp +++ b/include/graphblas/reference/blas3.hpp @@ -62,6 +62,203 @@ namespace grb { namespace internal { + /** + * \internal general select implementation that + * all select variants refer to + */ + template< + Descriptor descr, + class SelectionOperator, + typename Tin, + typename RITin, typename CITin, typename NITin, + typename Tout, + typename RITout, typename CITout, typename NITout + > + RC select_generic( + Matrix< Tout, reference, RITout, CITout, NITout > &out, + const Matrix< Tin, reference, RITin, CITin, NITin > &in, + const SelectionOperator &op, + const std::function< size_t( size_t ) > &row_l2g, + const std::function< size_t( size_t ) > &col_l2g, + const Phase &phase + ) { + typedef typename std::conditional< + std::is_void< Tin >::value, + bool, + Tin + >::type InValuesType; + // the identity will only be used for void input matrices + // if unused, we initialise it via default-construction, which is always ok + // to do since value types must be POD. + const InValuesType identity = std::is_void< Tin >::value + ? true + : InValuesType(); + + constexpr bool crs_only = descr & descriptors::force_row_major; + constexpr bool transpose_input = descr & descriptors::transpose_matrix; + static_assert( !(crs_only && transpose_input), + "The descriptors::force_row_major and descriptors::transpose_matrix " + "flags cannot be used simultaneously" ); + + const auto &in_raw = transpose_input + ? internal::getCCS( in ) + : internal::getCRS( in ); + auto &out_crs = internal::getCRS( out ); + auto &out_ccs = internal::getCCS( out ); + const size_t + m = transpose_input ? ncols( in ) : nrows( in ), + n = transpose_input ? nrows( in ) : ncols( in ); + const size_t + m_out = nrows( out ), + n_out = ncols( out ); + + typedef typename std::conditional< transpose_input, CITin, RITin >::type + EffectiveRowType; + typedef typename std::conditional< transpose_input, RITin, CITin >::type + EffectiveColumnType; + + // Check if the dimensions fit + if( m != m_out || n != n_out ) { + return MISMATCH; + } + + if( m == 0 || n == 0 || nnz(in ) == 0 ) { + return SUCCESS; + } + + if( phase == RESIZE ) { + size_t nzc = 0; + +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + const size_t nthreads = m < config::OMP::minLoopSize() + ? 1 + : config::OMP::threads(); + #pragma omp parallel reduction(+: nzc) num_threads( nthreads ) +#endif + { + size_t start_row = 0; + size_t end_row = m; +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + config::OMP::localRange( start_row, end_row, 0, m ); +#endif + size_t local_nzc = 0; + for( size_t i = start_row; i < end_row; ++i ) { + for( + size_t k = in_raw.col_start[ i ]; + k < in_raw.col_start[ i + 1 ]; + ++k + ) { + const auto j = in_raw.row_index[ k ]; + const auto global_row = + static_cast< EffectiveRowType >( row_l2g( i ) ); + const auto global_col = + static_cast< EffectiveColumnType >( col_l2g( j ) ); + const auto value = in_raw.getValue( k, identity ); + if( op( global_row, global_col, value ) ) { + (void) ++local_nzc; + } + } + } + nzc += local_nzc; + } + return grb::resize( out, nzc ); + } + + // Execute phase only from here on + assert( phase == EXECUTE ); + + // Declare the column counter array + config::NonzeroIndexType * col_counter = nullptr; + + if( !crs_only ) { + // Allocate the column counter array + char *arr = nullptr, *buf = nullptr; + Tin *valbuf = nullptr; + internal::getMatrixBuffers( arr, buf, valbuf, 1, out ); + col_counter = + internal::getReferenceBuffer< config::NonzeroIndexType >( n + 1 ); + +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + const size_t nthreads = ( n + 1 ) < config::OMP::minLoopSize() + ? 1 + : config::OMP::threads(); + #pragma omp parallel for simd num_threads( nthreads ) +#endif + for( size_t j = 0; j < n + 1; ++j ) { + out_ccs.col_start[ j ] = 0; + } + + // Fill CCS.col_start with the number of nonzeros in each column + for( size_t i = 0; i < m; ++i ) { + for( size_t k = in_raw.col_start[ i ]; k < in_raw.col_start[ i + 1 ]; ++k ) { + const auto j = in_raw.row_index[ k ]; + const auto global_row = static_cast< EffectiveRowType >( row_l2g( i ) ); + const auto global_col = + static_cast< EffectiveColumnType >( col_l2g( j ) ); + const auto value = in_raw.getValue( k, identity ); + if( op( global_row, global_col, value ) ) { + (void) ++out_ccs.col_start[ j + 1 ]; + } + } + } + + // Prefix sum of CCS.col_start + for( size_t j = 1; j < n + 1; ++j ) { + out_ccs.col_start[ j ] += out_ccs.col_start[ j - 1 ]; + } + + { // Initialise the column counter array with zeros +#ifdef _H_GRB_REFERENCE_OMP_BLAS3 + const size_t nthreads = ( n + 1 ) < config::OMP::minLoopSize() + ? 1 + : config::OMP::threads(); + #pragma omp parallel for simd num_threads( nthreads ) +#endif + for( size_t j = 0; j < n + 1; ++j ) { + col_counter[ j ] = 0; + } + } + } + + out_crs.col_start[ 0 ] = 0; + size_t nzc = 0; + for( size_t i = 0; i < m; ++i ) { + for( size_t k = in_raw.col_start[ i ]; k < in_raw.col_start[ i + 1 ]; + ++k + ) { + const auto j = in_raw.row_index[ k ]; + const auto global_row = static_cast< EffectiveRowType >( row_l2g( i ) ); + const auto global_col = static_cast< EffectiveColumnType >( col_l2g( j ) ); + const auto value = in_raw.getValue( k, identity ); + if( !op( global_row, global_col, value ) ) { + continue; + } +#ifdef _DEBUG + std::cout << "\tKeeping value at: " << row_l2g( i ) << ", " + << col_l2g( j ) << " -> idx=" << nzc << "\n"; +#endif + + // Update CCS + if( !crs_only ) { + const auto idx = out_ccs.col_start[ j ] + col_counter[ j ]; + (void) ++col_counter[ j ]; + out_ccs.row_index[ idx ] = i; + out_ccs.setValue( idx, value ); + } + + // Update CRS + out_crs.row_index[ nzc ] = j; + out_crs.setValue( nzc, value ); + (void) ++nzc; + } + out_crs.col_start[ i + 1 ] = nzc; + } + + internal::setCurrentNonzeroes( out, nzc ); + + return SUCCESS; + } + /** * \internal general mxm implementation that all mxm variants refer to */ @@ -1663,6 +1860,48 @@ namespace grb { ); } + template< + Descriptor descr = descriptors::no_operation, + class SelectionOperator, + typename Tin, + typename RITin, typename CITin, typename NITin, + typename Tout, + typename RITout, typename CITout, typename NITout + > + RC select( + Matrix< Tout, reference, RITout, CITout, NITout > &out, + const Matrix< Tin, reference, RITin, CITin, NITin > &in, + const SelectionOperator &op = SelectionOperator(), + const Phase &phase = EXECUTE, + const typename std::enable_if< + !is_object< Tin >::value && + !is_object< Tout >::value + >::type * const = nullptr + ) { +#ifdef _DEBUG + std::cout << "In grb::select( reference )\n"; +#endif + // static sanity checks + static_assert( + !(descr & descriptors::no_casting) && ( + (std::is_void< Tin >::value && std::is_void< Tout >::value) || + std::is_same< Tin, Tout >::value + ), + "grb::select (reference): " + "input and output matrix types must match" + ); + + // dispatch + return internal::select_generic< descr >( + out, + in, + op, + []( size_t i ) { return i; }, + []( size_t j ) { return j; }, + phase + ); + } + template< Descriptor descr = descriptors::no_operation, class Monoid, diff --git a/include/graphblas/selection_ops.hpp b/include/graphblas/selection_ops.hpp new file mode 100644 index 000000000..efa188a5d --- /dev/null +++ b/include/graphblas/selection_ops.hpp @@ -0,0 +1,126 @@ + +/* + * Copyright 2024 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. + */ + +/** + * @file + * + * Provides a set of standard binary selection operators. + * + * @author Benjamin D. Lozes + * @date 27th of February 2024 + */ + +#ifndef _H_GRB_SELECTION_OPERATORS +#define _H_GRB_SELECTION_OPERATORS + +#include "internalops.hpp" +#include "type_traits.hpp" + +namespace grb { + + /** + * This namespace holds various standard selection-operators such + * as #grb::operators::is_diagonal or #grb::operators::is_strictly_lower. + */ + namespace operators::select { + + template< + typename D, + typename RIT = config::RowIndexType, + typename CIT = config::ColIndexType + > + struct is_diagonal : public internal::MatrixSelectionOperatorBase< + internal::is_diagonal< D, RIT, CIT >, D + > {}; + + template< + typename D, + typename RIT = config::RowIndexType, + typename CIT = config::ColIndexType + > + struct is_strictly_lower : public internal::MatrixSelectionOperatorBase< + internal::is_strictly_lower< D, RIT, CIT >, D + > {}; + + template< + typename D, + typename RIT = config::RowIndexType, + typename CIT = config::ColIndexType + > + struct is_lower_or_diagonal : public internal::MatrixSelectionOperatorBase< + internal::is_lower_or_diagonal< D, RIT, CIT >, D + > {}; + + template< + typename D, + typename RIT = config::RowIndexType, + typename CIT = config::ColIndexType + > + struct is_strictly_upper : public internal::MatrixSelectionOperatorBase< + internal::is_strictly_upper< D, RIT, CIT >, D + > {}; + + template< + typename D, + typename RIT = config::RowIndexType, + typename CIT = config::ColIndexType + > + struct is_upper_or_diagonal : public internal::MatrixSelectionOperatorBase< + internal::is_upper_or_diagonal< D, RIT, CIT >, D + > {}; + + } // namespace operators::select + + template< typename D1, typename D2, typename D3 > + struct is_matrix_selection_operator< + operators::select::is_diagonal< D1, D2, D3 > + > { + static constexpr bool value = true; + }; + + template< typename D1, typename D2, typename D3 > + struct is_matrix_selection_operator< + operators::select::is_strictly_lower< D1, D2, D3 > + > { + static constexpr bool value = true; + }; + + template< typename D1, typename D2, typename D3 > + struct is_matrix_selection_operator< + operators::select::is_lower_or_diagonal< D1, D2, D3 > + > { + static constexpr bool value = true; + }; + + template< typename D1, typename D2, typename D3 > + struct is_matrix_selection_operator< + operators::select::is_strictly_upper< D1, D2, D3 > + > { + static constexpr bool value = true; + }; + + template< typename D1, typename D2, typename D3 > + struct is_matrix_selection_operator< + operators::select::is_upper_or_diagonal< D1, D2, D3 > + > { + static constexpr bool value = true; + }; + +} // end namespace grb + +#endif // end ``_H_GRB_SELECTION_OPERATORS'' + diff --git a/include/graphblas/type_traits.hpp b/include/graphblas/type_traits.hpp index 1c7fdbd7e..44a94e402 100644 --- a/include/graphblas/type_traits.hpp +++ b/include/graphblas/type_traits.hpp @@ -111,6 +111,24 @@ namespace grb { static const constexpr bool value = false; }; + /** + * Used to inspect whether a given type is an ALP matrix selection operator. + * + * @tparam T The type to inspect. + * + * \ingroup typeTraits + */ + template< typename T > + struct is_matrix_selection_operator { + + /** + * Whether \a T is an ALP operator. + * + * \internal Base case: an arbitrary type is not an ALP operator. + */ + static constexpr bool value = false; + }; + /** * Used to inspect whether a given type is an ALP/GraphBLAS object. * diff --git a/src/graphblas/hyperdags/hyperdags.cpp b/src/graphblas/hyperdags/hyperdags.cpp index 90774d0e2..256b31fba 100644 --- a/src/graphblas/hyperdags/hyperdags.cpp +++ b/src/graphblas/hyperdags/hyperdags.cpp @@ -380,6 +380,9 @@ std::string grb::internal::hyperdags::toString( case GETID_MATRIX: return "getID( matrix )"; + case SELECT_MATRIX_MATRIX: + return "select( matrix, matrix, selection_operator )"; + case FOLDL_SCALAR_MATRIX_MASK_MONOID: return "foldl( scalar, matrix, matrix, monoid )"; diff --git a/src/graphblas/reference/config.cpp b/src/graphblas/reference/config.cpp index d957f899f..cdda0a5af 100644 --- a/src/graphblas/reference/config.cpp +++ b/src/graphblas/reference/config.cpp @@ -17,6 +17,7 @@ #include "graphblas/reference/config.hpp" + using namespace grb; using namespace config; @@ -29,7 +30,7 @@ std::string grb::config::toString( const ALLOC_MODE mode ) { } assert( false ); std::cerr << "Warning: unknown memory allocation mode passed to " - "grb::config::toString." - << std::endl; + << "grb::config::toString." << std::endl; return "unknown"; } + diff --git a/src/graphblas/reference/init.cpp b/src/graphblas/reference/init.cpp index c5e035a98..03352356e 100644 --- a/src/graphblas/reference/init.cpp +++ b/src/graphblas/reference/init.cpp @@ -43,7 +43,10 @@ size_t grb::internal::reference_bufsize = 0; static grb::utils::AutoDeleter< size_t > privateSizetOMP_deleter; template<> -grb::RC grb::init< grb::reference >( const size_t s, const size_t P, void * const data ) { +grb::RC grb::init< grb::reference >( + const size_t s, const size_t P, + void * const data +) { // we don't use any implementation-specific init data (void)data; // print output @@ -79,15 +82,25 @@ grb::RC grb::init< grb::reference_omp >( const size_t s, const size_t P, void * RC rc = grb::SUCCESS; // print output const auto T = config::OMP::threads(); - std::cerr << "Info: grb::init (reference_omp) called. OpenMP is set to utilise " << T << " threads.\n"; - rc = grb::utils::alloc( "", "", grb::internal::privateSizetOMP, T * sizeof( grb::config::CACHE_LINE_SIZE::value() ) * sizeof( size_t ), true, privateSizetOMP_deleter ); + std::cerr << "Info: grb::init (reference_omp) called. OpenMP is set to " + << "utilise " << T << " threads.\n"; + rc = grb::utils::alloc( + "", + "", + grb::internal::privateSizetOMP, + T * sizeof( grb::config::CACHE_LINE_SIZE::value() ) * sizeof( size_t ), + true, + privateSizetOMP_deleter + ); // use same initialisation procedure as sequential implementation if( rc == grb::SUCCESS ) { rc = grb::init< grb::reference >( s, P, data ); } // pre-reserve a minimum buffer equal to the combined L1 cache size if( rc == grb::SUCCESS ) { - if( ! internal::template ensureReferenceBufsize< char >( T * config::MEMORY::l1_cache_size() ) ) { + if( !internal::template ensureReferenceBufsize< char >( + T * config::MEMORY::l1_cache_size() ) + ) { rc = grb::OUTOFMEM; } } 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 548f1c17c..fe6ddee06 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -65,6 +65,11 @@ add_grb_executables( distribution_bsp1d distribution_bsp1d.cpp BACKENDS reference NO_BACKEND_NAME ) +add_grb_executables( distribution_matrix_bsp1d distribution_matrix_bsp1d.cpp + BACKENDS bsp1d NO_BACKEND_NAME + ADDITIONAL_LINK_LIBRARIES test_utils +) + add_grb_executables( distribution distribution.cpp BACKENDS bsp1d NO_BACKEND_NAME ) @@ -179,6 +184,11 @@ add_grb_executables( RBGaussSeidel RBGaussSeidel.cpp BACKENDS reference reference_omp bsp1d hybrid hyperdags nonblocking ) +add_grb_executables( selectMatrix selectMatrix.cpp + BACKENDS reference reference_omp bsp1d hybrid hyperdags nonblocking + ADDITIONAL_LINK_LIBRARIES test_utils_headers +) + add_grb_executables( set set.cpp BACKENDS reference reference_omp bsp1d hybrid hyperdags nonblocking ) diff --git a/tests/unit/distribution_matrix_bsp1d.cpp b/tests/unit/distribution_matrix_bsp1d.cpp new file mode 100644 index 000000000..e9feb4c0b --- /dev/null +++ b/tests/unit/distribution_matrix_bsp1d.cpp @@ -0,0 +1,153 @@ +/* + * Copyright 2024 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 + +#include + + +using namespace grb; + +template< typename LpfData > +void print_local( const LpfData &lpf_data, std::stringstream &local_ss ) { + for( size_t p = 0; p < lpf_data.P; ++p ) { + spmd< BSP1D >::barrier(); + if( p == lpf_data.s ) { + if( local_ss.str().empty() ) { + std::cerr << "Process " << lpf_data.s << ": [nothing to print]\n"; + } else { + std::cerr << "Process " << lpf_data.s << ":\n" << local_ss.str() << "\n"; + } + } + usleep( 1000 ); + } + local_ss.clear(); +} + +template< typename D > +void grb_program( const size_t &n, RC &rc ) { + rc = SUCCESS; + Matrix< D, BSP1D > I_distributed( n, n, n ); + { + // Build the identity matrix + std::vector< D > values( n, 1 ); + std::vector< size_t > iota_indices( n, 0 ); + std::iota( iota_indices.begin(), iota_indices.end(), 0 ); + rc = buildMatrixUnique( I_distributed, iota_indices.data(), + iota_indices.data(), values.data(), n, SEQUENTIAL ); + } + + // Each process now have to check if the global + // coordinates match, since it should always be: i == j + const Matrix< D, reference > &local_matrix = + internal::getLocal( I_distributed ); + const auto &lpf_data = internal::grb_BSP1D.cload(); + + std::stringstream local_ss; + { + // Check for the CRS + const auto &crs = internal::getCRS( local_matrix ); + for(size_t i = 0; i < nrows( local_matrix ); ++i ) { + for( size_t k = crs.col_start[i]; k < crs.col_start[ i + 1 ]; ++k ) { + const auto j = crs.row_index[k]; // Local AND global since the distribution is 1D + + const size_t col_pid = internal::Distribution<>::offset_to_pid( j, + ncols(I_distributed), lpf_data.P ); + const size_t col_off = internal::Distribution<>::local_offset( + ncols(I_distributed), col_pid, lpf_data.P ); + const auto global_i = internal::Distribution<>::local_index_to_global( i, + nrows(I_distributed), lpf_data.s, lpf_data.P ); + const auto global_j = internal::Distribution<>::local_index_to_global( + j - col_off, ncols(I_distributed), col_pid, lpf_data.P ); + + if (global_i != global_j) { + local_ss << " Wrong coordinate in CRS found at: ( " + << std::setw(3) << i << ", " << std::setw(3) << global_j << " ) " + << "--(mapped to global)--> ( " + << std::setw(3) << global_i << ", "<< std::setw(3) << global_j + << " )\n"; + rc = FAILED; + } + } + } + } + + print_local( lpf_data, local_ss ); + + { // Check for the CCS + const auto &ccs = internal::getCCS( local_matrix ); + for( size_t j = 0; j < ncols(local_matrix); ++j ) { + for( size_t k = ccs.col_start[j]; k < ccs.col_start[j+1]; ++k ) { + const auto i = ccs.row_index[k]; + + const size_t col_pid = internal::Distribution<>::offset_to_pid( j, + ncols( I_distributed ), lpf_data.P ); + const size_t col_off = internal::Distribution<>::local_offset( + ncols(I_distributed), col_pid, lpf_data.P ); + const auto global_i = internal::Distribution<>::local_index_to_global( i, + nrows(I_distributed), lpf_data.s, lpf_data.P ); + const auto global_j = internal::Distribution<>::local_index_to_global( + j - col_off, ncols(I_distributed), col_pid, lpf_data.P ); + + if( global_i != global_j ) { + local_ss << " Wrong coordinate in CCS found at: ( " + << std::setw(3) << i << ", " << std::setw(3) << global_j << " ) " + << "--(mapped to global)--> ( " + << std::setw(3) << global_i << ", " << std::setw(3) << global_j << " )\n"; + rc = RC::FAILED; + } + } + } + } + + print_local( lpf_data, local_ss ); + + if( collectives<>::allreduce( rc, operators::any_or() ) != SUCCESS ) { + rc = PANIC; + } +} + +int main( int argc, char * * argv ) { + // defaults + if(argc != 2 ) { + std::cerr << "Usage: " << argv[0] << " \n"; + return 0; + } + std::cout << "This is functional test " << argv[0] << "\n"; + + Launcher< AUTOMATIC > launcher; + RC out = SUCCESS; + size_t n = std::strtoul( argv[1], nullptr, 10 ); + + const RC launch_rc = launcher.exec( &grb_program< int >, n, out, true ); + if( launch_rc != SUCCESS ) { + std::cerr << "Launch test failed" << std::endl; + std::cout << "Test FAILED\n" << std::endl; + return 10; + } + + if( out != SUCCESS ) { + std::cerr << std::flush; + std::cout << "Test FAILED (" << toString(out) << ")\n" << std::endl; + return 10 + static_cast< int >(launch_rc); + } + + std::cout << "Test OK\n" << std::endl; + return 0; +} + diff --git a/tests/unit/fold_matrix_to_scalar.cpp b/tests/unit/fold_matrix_to_scalar.cpp index e97f67ff7..00a1f524e 100644 --- a/tests/unit/fold_matrix_to_scalar.cpp +++ b/tests/unit/fold_matrix_to_scalar.cpp @@ -27,9 +27,9 @@ * Tests whether the foldl and foldl API calls produce the expected results. * * The test cases are focused on the following aspects: - * * The types of the result, the matrix values and the operator - * * The initial value of the reduction result - * * The order of the operands (foldr, foldl) + * - The types of the result, the matrix values and the operator + * - The initial value of the reduction result + * - The order of the operands (foldr, foldl) */ #include diff --git a/tests/unit/selectMatrix.cpp b/tests/unit/selectMatrix.cpp new file mode 100644 index 000000000..c6d80e93d --- /dev/null +++ b/tests/unit/selectMatrix.cpp @@ -0,0 +1,653 @@ +/* + * Copyright 2021 Huawei Technologies Co., Ltd. + * + * Licensed under the Apache License, Version 2.0 (the "License" ); + * you may !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 +#include +#include + +#include +#include + +using namespace grb; + +#define STDERR_WITH_LINE std::cerr << "[Line " << __LINE__ << "] " + +constexpr bool Debug = false; + +namespace { + + template< bool enabled, typename D > + void printSparseMatrix( + const Matrix< D > &mat, + const std::string &name = "", + std::ostream &os = std::cerr + ) { + if( !enabled ) return; + wait( mat ); + print_matrix( mat, 256, name, os ); + } + +} // namespace + +template< typename D > +std::tuple< + const typename config::RowIndexType, + const typename config::ColIndexType, + D +> getMatrixEntry( + const std::pair< + const std::pair< config::RowIndexType, config::ColIndexType >, + D + > &entry, + const typename std::enable_if< + !std::is_void< D >::value + >::type * = nullptr +) { + return std::make_tuple( + entry.first.first, + entry.first.second, + entry.second + ); +} + +template< typename D > +std::tuple< const size_t, const size_t, bool > getMatrixEntry( + const std::pair< const size_t, const size_t > &entry, + const typename std::enable_if< std::is_void< D >::value >::type * = nullptr +) { + return std::make_tuple( + entry.first, + entry.second, + true + ); +} + +template< + typename T, typename RIT, typename CIT, Backend implementation, + typename Predicate +> +size_t count_nnz_if( + const Matrix< T, implementation, RIT, CIT > &mat, + Predicate predicate, + const typename std::enable_if< !std::is_void< T >::value >::type * = nullptr +) { + size_t count = 0; + for( const auto &each : mat ) { + const auto entry = getMatrixEntry< T >( each ); + const auto r = std::get<0>( entry ); + const auto c = std::get<1>( entry ); + const auto v = std::get<2>( entry ); + if( predicate( r, c, v ) ) { + (void) ++count; + } + } + const RC rc = collectives<>::allreduce( count, operators::add< size_t >() ); + if( rc != SUCCESS ) { + throw std::runtime_error( "count_nnz_if: could not all-reduce final count" ); + } + return count; +} + +template< + typename T, typename RIT, typename CIT, Backend implementation, + typename Predicate +> +size_t count_nnz_if( + const Matrix< T, implementation, RIT, CIT > &mat, + Predicate predicate, + const typename std::enable_if< std::is_void< T >::value >::type * = nullptr +) { + size_t count = 0; + for( const auto &each : mat ) { + const auto entry = getMatrixEntry< T >( each ); + const auto r = std::get< 0 >( entry ); + const auto c = std::get< 1 >( entry ); + const auto v = std::get< 2 >( entry ); + if( predicate( r, c, v ) ) { + (void) ++count; + } + } + const RC rc = collectives<>::allreduce( count, operators::add< size_t >() ); + if( rc != SUCCESS ) { + throw std::runtime_error( "count_nnz_if: could not all-reduce final count" ); + } + return count; +} + +template< + typename T, + typename D, + typename Func, + typename RIT, typename CIT, + Backend implementation +> +bool matrix_validate_predicate( + const Matrix< T, implementation, RIT, CIT > &src, + const Matrix< D, implementation, RIT, CIT > &obtained, + Func predicate, + typename std::enable_if< !std::is_void< D >::value >::type * = nullptr +) { + const size_t expected_nvals = count_nnz_if( src, predicate ); + if( expected_nvals != nnz( obtained ) ) { + std::cerr << " /!\\ Expected " << expected_nvals + << " non-zero entries, but obtained " << nnz( obtained ) + << std::endl; + return false; + } + + bool valid = true; + for( const auto &each : obtained ) { + const auto entry = getMatrixEntry< D >(each); + const auto r = std::get<0>( entry ); + const auto c = std::get<1>( entry ); + const auto v = std::get<2>( entry ); + const bool match = predicate( r, c, v ); + if( !match ) { + std::cerr << " /!\\ Predicate failed for (" + << std::get<0>( entry ) << ", " << std::get<1>( entry ) + << ", " << std::get<2>( entry ) << ")" + << std::endl; + valid = false; + break; + } + } + if( + collectives<>::allreduce( valid, operators::logical_and< bool >() ) + != SUCCESS + ) { + return false; + } + + return valid; +} + +template< + typename T, + typename D, + typename Func, + typename RIT, typename CIT, + Backend implementation +> +bool matrix_validate_predicate( + const Matrix< T, implementation, RIT, CIT > &src, + const Matrix< D, implementation, RIT, CIT > &obtained, + Func predicate, + typename std::enable_if< std::is_void< D >::value >::type * = nullptr +) { + const size_t expected_nvals = count_nnz_if( src, predicate ); + if( expected_nvals != nnz( obtained ) ) { + std::cerr << " /!\\ Expected " << expected_nvals + << " non-zero entries, but obtained " << nnz( obtained ) + << std::endl; + return false; + } + + bool valid = true; + for( const auto &each : obtained ) { + const auto entry = getMatrixEntry< D >( each ); + const auto r = std::get<0>( entry ); + const auto c = std::get<1>( entry ); + const auto v = std::get<2>( entry ); + const bool match = predicate( r, c, v ); + if( !match ) { + std::cerr << " /!\\ Predicate failed for (" + << std::get<0>( entry ) << ", " << std::get<1>( entry ) + << ", " << std::get<2>( entry ) << ")" + << std::endl; + valid = false; + break; + } + } + if( + collectives<>::allreduce( valid, operators::logical_and< bool >() ) + != SUCCESS + ) { + return false; + } + + return valid; +} + +template< + typename D, + typename SelectionOperator, + typename RIT, typename CIT, + Backend implementation +> +RC test_case( + const Matrix< D, implementation, RIT, CIT > &input, + const SelectionOperator op, + const std::string &test_name +) { + std::cout << test_name << std::endl; + + Matrix< D, implementation, RIT, CIT > output( + nrows( input ), ncols( input ), 0 + ); + + RC rc = select( output, input, op, RESIZE ); + if( rc != SUCCESS ) { + std::cerr << "(non-lambda variant): RESIZE phase of test <" + << test_name << "> failed, rc is \"" + << toString(rc) << "\"" << std::endl; + return rc; + } + + rc = select( output, input, op, EXECUTE ); + if( rc != SUCCESS ) { + std::cerr << "(non-lambda variant): EXECUTE phase of test <" + << test_name << "> failed, rc is \"" + << toString(rc) << "\"" << std::endl; + return rc; + } + + printSparseMatrix< Debug >( output ); + + const bool valid = matrix_validate_predicate( input, output, op ); + if( !valid ) { + std::cerr << "(non-lambda variant): Test <" + << test_name << "> failed, output matrix is invalid" + << std::endl; + return FAILED; + } + + return SUCCESS; +} + +template< + typename D, + typename RIT, typename CIT, typename NIT, + Backend implementation +> +RC buildMatrixUniqueWrapper( + Matrix< D, implementation, RIT, CIT, NIT > &mat, + const size_t * row_indices, + const size_t * col_indices, + const size_t nvals, + const IOMode io_mode, + const typename std::enable_if< + !std::is_void< D >::value + >::type * const = nullptr +) { + std::vector< D > values( nvals, 1 ); + return buildMatrixUnique( mat, row_indices, + col_indices, values.data(), nvals, io_mode ); +} + +template< + typename D, + typename RIT, typename CIT, typename NIT, + Backend implementation +> +RC buildMatrixUniqueWrapper( + Matrix< D, implementation, RIT, CIT, NIT > &mat, + const size_t * row_indices, + const size_t * col_indices, + const size_t nvals, + const IOMode io_mode, + const typename std::enable_if< + std::is_void< D >::value + >::type * const = nullptr +) { + return buildMatrixUnique( mat, row_indices, + col_indices, nvals, io_mode ); +} + +template< typename D > +RC buildMatrices( + grb::Matrix< D > &I, + grb::Matrix< D > &I_transposed, + grb::Matrix< D > &One_row, + grb::Matrix< D > &One_col, + const size_t n +) { + std::vector< size_t > const_indices_zero( n, 0 ); + std::vector< size_t > iota_indices( n, 0 ); + std::iota( iota_indices.begin(), iota_indices.end(), 0 ); + std::vector< size_t > reverse_iota_indices( n, 0 ); + for ( size_t i = 0; i < n; ++i ) { + reverse_iota_indices[i] = n - i - 1; + } + + RC rc = buildMatrixUniqueWrapper( + I, iota_indices.data(), + iota_indices.data(), n, SEQUENTIAL + ); + printSparseMatrix< Debug >( I, "identity" ); + + rc = rc + ? rc + : buildMatrixUniqueWrapper( + I_transposed, iota_indices.data(), + reverse_iota_indices.data(), n, SEQUENTIAL + ); + printSparseMatrix< Debug >( I_transposed, "transposed-identity" ); + + rc = rc + ? rc + : buildMatrixUniqueWrapper( + One_row, const_indices_zero.data(), + iota_indices.data(), n, SEQUENTIAL + ); + printSparseMatrix< Debug >( One_row, "one-row" ); + + rc = rc + ? rc + : buildMatrixUniqueWrapper( + One_col, iota_indices.data(), + const_indices_zero.data(), n, SEQUENTIAL + ); + printSparseMatrix< Debug >( One_col, "one-column" ); + + return rc; +} + +template< typename D > +void grb_program_operators( const size_t &n, RC &rc ) { + std::cerr << "Building matrices" << std::endl; + const std::string D_name = std::is_void< D >::value + ? "void" + : "non-void"; + + Matrix< D > + I( n, n, n ), + I_transposed( n, n, n ), + One_row( n, n, n ), + One_col( n, n, n ); + + // Build matrices + rc = buildMatrices( I, I_transposed, One_row, One_col, n ); + std::cerr << "Matrices built" << std::endl; + + // Test 01: Select + rc = rc ? rc : test_case( I, operators::select::is_diagonal< D >(), + "Test 01: Select > out of " ); + rc = rc ? rc : test_case( I_transposed, operators::select::is_diagonal< D >(), + "Test 01: Select > out of " ); + rc = rc ? rc : test_case( One_row, operators::select::is_diagonal< D >(), + "Test 01: Select > out of " ); + rc = rc ? rc : test_case( One_col, operators::select::is_diagonal< D >(), + "Test 01: Select > out of " ); + + // Test 02: Select + rc = rc ? rc : test_case( I, operators::select::is_strictly_lower< D >(), + "Test 02: Select > out of " ); + rc = rc ? rc : test_case( I_transposed, + operators::select::is_strictly_lower< D >(), + "Test 02: Select > out of " ); + rc = rc ? rc : test_case( One_row, operators::select::is_strictly_lower< D >(), + "Test 02: Select > out of " ); + rc = rc ? rc : test_case( One_col, operators::select::is_strictly_lower< D >(), + "Test 02: Select > out of " ); + + // Test 03: Select + rc = rc ? rc : test_case( I, operators::select::is_strictly_upper< D >(), + "Test 03: Select > out of " ); + rc = rc ? rc : test_case( I_transposed, + operators::select::is_strictly_upper< D >(), + "Test 03: Select > out of " ); + rc = rc ? rc : test_case( One_row, operators::select::is_strictly_upper< D >(), + "Test 03: Select > out of " ); + rc = rc ? rc : test_case( One_col, operators::select::is_strictly_upper< D >(), + "Test 03: Select > out of " ); + + // Test 04: Select + rc = rc ? rc : test_case( I, operators::select::is_lower_or_diagonal< D >(), + "Test 04: Select > out of " ); + rc = rc ? rc : test_case( I_transposed, + operators::select::is_lower_or_diagonal< D >(), + "Test 04: Select > out of " ); + rc = rc ? rc : test_case( One_row, + operators::select::is_lower_or_diagonal< D >(), + "Test 04: Select > out of " ); + rc = rc ? rc : test_case( One_col, + operators::select::is_lower_or_diagonal< D >(), + "Test 04: Select > out of " ); + + // Test 05: Select + rc = rc ? rc : test_case( I, operators::select::is_upper_or_diagonal< D >(), + "Test 05: Select > out of " ); + rc = rc ? rc : test_case( I_transposed, + operators::select::is_upper_or_diagonal< D >(), + "Test 05: Select > out of " ); + rc = rc ? rc : test_case( One_row, + operators::select::is_upper_or_diagonal< D >(), + "Test 05: Select > out of " ); + rc = rc ? rc : test_case( One_col, + operators::select::is_upper_or_diagonal< D >(), + "Test 05: Select > out of " ); + + // done + if( collectives<>::allreduce( rc, operators::any_or< RC >() ) != SUCCESS ) { + rc = PANIC; + } +} + +template< typename D, typename RIT, typename CIT > +void grb_program_lambdas( const size_t &n, RC &rc ) { + typedef typename std::conditional< + std::is_void< D >::value, + bool, + D + >::type D_safe; + const std::string D_name = "non-void"; + + Matrix< D > + I( n, n, n ), + I_transposed( n, n, n ), + One_row( n, n, n ), + One_col( n, n, n ); + + // Build matrices + rc = buildMatrices( I, I_transposed, One_row, One_col, n ); + + // Test 06: Select using lambda function + rc = rc ? rc : test_case( I, + [](const RIT &i, const CIT &j, const D_safe &val) -> bool { + assert( !std::is_void< D >::value || val ); + (void) val; + return i <= j; + }, "Test 06: Select > out of " + ); + rc = rc ? rc : test_case( I_transposed, + [](const RIT &i, const CIT &j, const D_safe &val) -> bool { + assert( !std::is_void< D >::value || val ); + (void) val; + return i >= j; + }, "Test 06: Select > out of " + ); + rc = rc ? rc : test_case( One_row, + [](const RIT &i, const CIT &j, const D_safe &val) -> bool { + assert( !std::is_void< D >::value || val ); + (void) val; + return i >= j; + }, "Test 06: Select > out of " + ); + rc = rc ? rc : test_case( One_col, + [](const RIT &i, const CIT &j, const D_safe &val) -> bool { + assert( !std::is_void< D >::value || val ); + (void) val; + return i >= j; + }, "Test 06: Select > out of " + ); + + // Test 07: Select using lambda function + rc = rc ? rc : test_case( I, + [](const RIT &i, const CIT &j, const D_safe &val) -> bool { + assert( !std::is_void< D >::value || val ); + (void) val; + return i > j; + }, "Test 07: Select > out of " + ); + rc = rc ? rc : test_case( I_transposed, + [](const RIT &i, const CIT &j, const D_safe &val) -> bool { + assert( !std::is_void< D >::value || val ); + (void) val; + return i > j; + }, "Test 07: Select > out of " + ); + rc = rc ? rc : test_case( One_row, + [](const RIT &i, const CIT &j, const D_safe &val) -> bool { + assert( !std::is_void< D >::value || val ); + (void) val; + return i > j; + }, "Test 07: Select > out of " + ); + rc = rc ? rc : test_case( One_col, + [](const RIT &i, const CIT &j, const D_safe &val) -> bool { + assert( !std::is_void< D >::value || val ); + (void) val; + return i > j; + }, "Test 07: Select > out of " + ); + + // done + if( collectives<>::allreduce( rc, operators::any_or< RC >() ) != SUCCESS ) { + rc = PANIC; + } +} + +template< typename RIT, typename CIT, typename LauncherT > +RC test_lambda_versions( const size_t &n, LauncherT &launcher ) { + RC out = SUCCESS; + { + std::cout << "-- -- Running lambda test with using matrix-type: double" + << std::endl; + if( launcher.exec( &grb_program_lambdas< + double, RIT, CIT + >, n, out, true ) != SUCCESS + ) { + STDERR_WITH_LINE << "Launching test FAILED\n"; + return PANIC; + } + if( out != SUCCESS ) { + STDERR_WITH_LINE << "Test FAILED (" << toString(out) << ")" << std::endl; + return out; + } + } + + { + std::cout << "-- -- Running lambda test with using matrix-type: void" + << std::endl; + if( launcher.exec( &grb_program_lambdas< + void, RIT, CIT + >, n, out, true ) != SUCCESS + ) { + STDERR_WITH_LINE << "Launching test FAILED\n"; + return PANIC; + } + if( out != SUCCESS ) { + STDERR_WITH_LINE << "Test FAILED (" << toString(out) << ")" << std::endl; + return out; + } + } + return out; +} + +int main( int argc, char** argv ) { + (void) argc; + (void) argv; + + RC out = SUCCESS; + + std::cout << "This is functional test " << argv[0] << "\n"; + Launcher< AUTOMATIC > launcher; + + const size_t n = argc > 1 + ? std::strtoul( argv[1], nullptr, 10 ) + : 1000; + + { + std::cout << "-- -- Running test with using matrix-type: int" << std::endl; + if( launcher.exec( &grb_program_operators< int >, n, out, true ) + != SUCCESS + ) { + STDERR_WITH_LINE << "Launching test FAILED\n"; + return 10; + } + if( out != SUCCESS ) { + STDERR_WITH_LINE << "Test FAILED (" << toString(out) << ")" << std::endl; + return 20 + static_cast< int >(out); + } + } + + { + std::cout << "-- -- Running test with using matrix-type: void" << std::endl; + if( launcher.exec( &grb_program_operators< void >, n, out, true ) + != SUCCESS + ) { + STDERR_WITH_LINE << "Launching test FAILED\n"; + return 30; + } + if( out != SUCCESS ) { + STDERR_WITH_LINE << "Test FAILED (" << toString(out) << ")" << std::endl; + return 40 + static_cast< int >(out); + } + } + + { + std::cout << "-- -- Running test using lambda functions (matching index types)" + << std::endl; + assert( out == SUCCESS ); + out = test_lambda_versions< + grb::config::RowIndexType, + grb::config::ColIndexType + >( n, launcher ); + if( out != SUCCESS ) { + STDERR_WITH_LINE << "Test FAILED (" << toString(out) << ")" << std::endl; + return 50 + static_cast< int >(out); + } + } + + { + std::cout << "-- -- Running test using lambda functions (mis-matching index " + << "types)" << std::endl; + assert( out == SUCCESS ); + out = test_lambda_versions< size_t, size_t >( n, launcher ); + if( out != SUCCESS ) { + STDERR_WITH_LINE << "Test FAILED (" << toString(out) << ")" << std::endl; + return 60 + static_cast< int >(out); + } + } + + std::cerr << std::flush; + std::cout << "Test OK" << std::endl; + return 0; +} + diff --git a/tests/unit/unittests.sh b/tests/unit/unittests.sh index f140b27b3..f70d486a2 100755 --- a/tests/unit/unittests.sh +++ b/tests/unit/unittests.sh @@ -616,6 +616,18 @@ for MODE in ${MODES}; do grep 'Test OK' ${TEST_OUT_DIR}/zip_large_${MODE}_${BACKEND}_${P}_${T} || echo "Test FAILED" echo " " + echo ">>> [x] [ ] Testing grb::select on matrices of integers and of size 3" + $runner ${TEST_BIN_DIR}/selectMatrix_${MODE}_${BACKEND} 3 &> ${TEST_OUT_DIR}/selectMatrix_${MODE}_${BACKEND}_${P}_${T}_3 + head -1 ${TEST_OUT_DIR}/selectMatrix_${MODE}_${BACKEND}_${P}_${T}_3 + grep 'Test OK' ${TEST_OUT_DIR}/selectMatrix_${MODE}_${BACKEND}_${P}_${T}_3 || echo "Test FAILED" + echo " " + + echo ">>> [x] [ ] Testing grb::select on matrices of integers and of size 5'000" + $runner ${TEST_BIN_DIR}/selectMatrix_${MODE}_${BACKEND} 5000 &> ${TEST_OUT_DIR}/selectMatrix_${MODE}_${BACKEND}_${P}_${T} + head -1 ${TEST_OUT_DIR}/selectMatrix_${MODE}_${BACKEND}_${P}_${T} + grep 'Test OK' ${TEST_OUT_DIR}/selectMatrix_${MODE}_${BACKEND}_${P}_${T} || echo "Test FAILED" + echo " " + echo ">>> [x] [ ] Testing copy-constructor of square pattern matrices" echo " of size 1003." $runner ${TEST_BIN_DIR}/copyVoidMatrices_${MODE}_${BACKEND} 1003 &> ${TEST_OUT_DIR}/copyVoidMatrices_${MODE}_${BACKEND}_${P}_${T} @@ -761,10 +773,16 @@ for MODE in ${MODES}; do if [ "$BACKEND" = "bsp1d" ]; then echo "Additional unit tests for the BSP1D backend:" echo " " - echo ">>> [x] [ ] Testing BSP1D distribution for a vector of size 100000" + echo ">>> [x] [ ] Testing BSP1D distribution for a vector of size 100 000" echo " " ${TEST_BIN_DIR}/distribution_bsp1d_${MODE} + echo ">>> [x] [ ] Testing BSP1D distribution for an identity matrix of size" + echo " 7777 x 7777. The test evaluates whether the internal data" + echo " structures match the BSP1D distribution" + echo " " + ${TEST_BIN_DIR}/distribution_matrix_bsp1d_${MODE} 7777 + echo ">>> [x] [ ] Testing dense vector times matrix using the double (+,*)" echo " semiring where matrix elements are doubles and vector" echo " elements ints. The input matrix is taken from west0497."