What a selection operator could look like

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")

include/graphblas.hpp

@@ -513,6 +513,7 @@ namespace grb {
 #include <graphblas/ops.hpp>
 #include <graphblas/monoid.hpp>
 #include <graphblas/semiring.hpp>
+#include <graphblas/selection_ops.hpp>
 
 // Then include containers. If containers rely on ALP/GraphBLAS primitives that
 // are defined as free functions, then container implementations must forward-

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 <tt>true</tt> 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 <tt>unsigned int</tt> 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 <tt>size_t</tt>, but such use may result in a (slight) performance
+	 * penalty due to internal casting between possibly different index types.
+	 *
+	 * For <tt>void</tt> matrices, the select operator will assume a <tt>bool</tt>
+	 * for \a T. The operator will always receive <tt>true</tt> 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.

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

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
+