268 Provide grb::set(Matrix out, Matrix mask, Matrix in)

include/graphblas/hyperdags/hyperdags.hpp

@@ -380,6 +380,8 @@ namespace grb {
 
 				SET_MATRIX_MATRIX_INPUT2,
 
+				SET_MATRIX_MATRIX_MASKED,
+
 				MXM_MATRIX_MATRIX_MATRIX_SEMIRING,
 
 				MXM_MATRIX_MATRIX_MATRIX_MONOID,
@@ -493,7 +495,7 @@ namespace grb {
 			};
 
 			/** \internal How many operation vertex types exist. */
-			const constexpr size_t numOperationVertexTypes = 106;
+			const constexpr size_t numOperationVertexTypes = 107;
 
 			/** \internal An array of all operation vertex types. */
 			const constexpr enum OperationVertexType

include/graphblas/hyperdags/io.hpp

@@ -403,6 +403,41 @@ namespace grb {
 		return ret;
 	}
 
+	template<
+		Descriptor descr = descriptors::no_operation,
+		typename OutputType, typename MaskType, typename InputType,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2
+	>
+	RC set(
+		Matrix< OutputType, hyperdags, RIT1, CIT1, NIT1 > &C,
+		const Matrix< MaskType, hyperdags, RIT2, CIT2, NIT2 > &M,
+		const Matrix< InputType, hyperdags, RIT2, CIT2, NIT2 > &A,
+		const Phase &phase = EXECUTE
+	) {
+		const RC ret = set< descr >(
+			internal::getMatrix( C ), internal::getMatrix( M ),
+			internal::getMatrix( A ), phase
+		);
+		if( ret != SUCCESS ) { return ret; }
+		if( phase != EXECUTE ) { return ret; }
+		if( nrows( A ) == 0 || ncols( A ) == 0 ) { return ret; }
+		std::array< const void *, 0 > sourcesP{};
+		std::array< uintptr_t, 3 > sourcesC{
+			getID( internal::getMatrix(A) ),
+			getID( internal::getMatrix(M) ),
+			getID( internal::getMatrix(C) )
+		};
+		std::array< uintptr_t, 1 > destinations{ getID( internal::getMatrix(C) ) };
+		internal::hyperdags::generator.addOperation(
+			internal::hyperdags::SET_MATRIX_MATRIX_MASKED,
+			sourcesP.begin(), sourcesP.end(),
+			sourcesC.begin(), sourcesC.end(),
+			destinations.begin(), destinations.end()
+		);
+		return ret;
+	}
+
 	template< typename DataType, typename Coords >
 	RC clear( Vector< DataType, hyperdags, Coords > &x ) {
 		const RC ret = clear( internal::getVector( x ) );

include/graphblas/nonblocking/io.hpp

@@ -1217,6 +1217,41 @@ namespace grb {
 		}
 	}
 
+
+	template<
+		Descriptor descr = descriptors::no_operation,
+		typename OutputType, typename MaskType, typename InputType,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2
+	>
+	RC set(
+		Matrix< OutputType, nonblocking, RIT1, CIT1, NIT1 > &C,
+		const Matrix< MaskType, nonblocking, RIT2, CIT2, NIT2 > &M,
+		const Matrix< InputType, nonblocking, RIT2, CIT2, NIT2 > &A,
+		const Phase &phase = EXECUTE
+	) noexcept {
+		if( internal::NONBLOCKING::warn_if_not_native &&
+			config::PIPELINE::warn_if_not_native
+		) {
+			std::cerr << "Warning: masked set (nonblocking) currently delegates "
+				<< "to a blocking implementation.\n"
+				<< "         Further similar such warnings will be suppressed.\n";
+			internal::NONBLOCKING::warn_if_not_native = false;
+		}
+
+		// nonblocking execution is not supported
+		// first, execute any computation that is not completed
+		internal::le.execution();
+
+		// second, delegate to the reference backend
+		return set< descr >(
+			internal::getRefMatrix( C ),
+			internal::getRefMatrix( M ),
+			internal::getRefMatrix( A ),
+			phase
+		);
+	}
+
 	template<
 		Descriptor descr = descriptors::no_operation,
 		typename InputType,

include/graphblas/reference/compressed_storage.hpp

@@ -425,6 +425,16 @@ namespace grb {
 					return values;
 				}
 
+				/**
+				 * @returns The value array.
+				 *
+				 * \warning Does not check for <tt>NULL</tt> pointers.
+				 */
+				D * getValues() const noexcept {
+					return values;
+				}
+
+
 				/**
 				 * @returns The index array.
 				 *
@@ -1105,6 +1115,13 @@ namespace grb {
 					return nullptr;
 				}
 
+				/**
+				 * @returns A null pointer (since this is a pattern matrix).
+				 */
+				char * getValues() const noexcept {
+					return nullptr;
+				}
+
 				/**
 				 * @returns The index array.
 				 *

include/graphblas/reference/io.hpp

@@ -1156,6 +1156,208 @@ namespace grb {
 		}
 	}
 
+	template<
+		Descriptor descr = descriptors::no_operation,
+		typename OutputType, typename MaskType, typename InputType,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2
+	>
+	RC set(
+		Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C,
+		const Matrix< MaskType, reference, RIT2, CIT2, NIT2 > &M,
+		const Matrix< InputType, reference, RIT2, CIT2, NIT2 > &A,
+		const Phase &phase = EXECUTE
+	) noexcept {
+		static_assert( !std::is_void< OutputType >::value,
+			"grb::set (masked set to matrix): cannot have a pattern "
+			"matrix as output" );
+		static_assert( std::is_convertible< InputType, OutputType >::value,
+			"grb::set (masked set to matrix): input type cannot be "
+			"converted to output type"
+		);
+		static_assert(
+			! ( descr & descriptors::structural && descr & descriptors::invert_mask ),
+			"grb::set can not be called with both descriptors::structural "
+			"and descriptors::invert_mask in the masked variant"
+		);
+#ifdef _DEBUG
+		std::cout << "Called grb::set (matrix-to-matrix-masked, reference)\n";
+#endif
+		// static checks
+		NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) ||
+				std::is_same< InputType, OutputType >::value
+			), "grb::set",
+			"called with non-matching value types"
+		);
+
+		// dynamic checks
+		assert( phase != TRY );
+
+		const size_t nrows = grb::nrows( C );
+		const size_t ncols = grb::ncols( C );
+
+		const size_t m = grb::nrows( M );
+		const size_t n = grb::ncols( M );
+
+		/*grb::Monoid<
+			grb::operators::left_assign< OutputType >,
+			grb::identities::zero
+		> dummyMonoid;*/
+
+		if( m == 0 || n == 0 ) {
+			// If the mask has a null size, it will be ignored
+			return set< descr >( C, A, phase );
+		}
+
+		if( nrows != grb::nrows( A ) || nrows != m ) {
+			return MISMATCH;
+		}
+
+		if( ncols != grb::ncols( A ) || ncols != n ) {
+			return MISMATCH;
+		}
+
+		const auto &mask_raw = internal::getCRS( M );
+
+		const auto &A_raw = internal::getCRS( A );
+
+		size_t nzc = 0;
+
+		char * mask_arr = nullptr;
+		char * mask_buf = nullptr;
+		MaskType * mask_valbuf = nullptr;
+		internal::getMatrixBuffers( mask_arr, mask_buf, mask_valbuf, 1, M );
+
+		internal::Coordinates< reference > mask_coors;
+		mask_coors.set( mask_arr, false, mask_buf, ncols );
+
+		for( size_t i = 0; i < nrows; ++i ) {
+			mask_coors.clear();
+			for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) {
+				const auto k_col = mask_raw.row_index[ k ];
+				if( utils::interpretMatrixMask< descr, MaskType >( true, mask_raw.getValues(), k ) ) {
+					mask_coors.assign( k_col );
+				}
+			}
+#ifdef _H_GRB_REFERENCE_OMP_BLAS3
+			#pragma omp parallel for reduction(+:nzc)
+#endif
+			for( auto k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) {
+				const auto k_col = A_raw.row_index[ k ];
+				if( mask_coors.assigned( k_col ) ) {
+					nzc++;
+				}
+			}
+		}
+
+		if( phase == RESIZE ) {
+			return resize( C, nzc );
+		}
+
+		assert( phase == EXECUTE );
+
+		if( capacity( C ) < nzc ) {
+#ifdef _DEBUG
+			std::cout << "\t insufficient capacity to complete "
+				"requested masked set matrix to matrix computation\n";
+#endif
+			const RC clear_rc = clear( C );
+			if( clear_rc != SUCCESS ) {
+				return PANIC;
+			} else {
+				return FAILED;
+			}
+		}
+
+		auto &CRS_raw = internal::getCRS( C );
+		auto &CCS_raw = internal::getCCS( C );
+
+		config::NonzeroIndexType * C_col_index = internal::template
+			getReferenceBuffer< typename config::NonzeroIndexType >( ncols + 1 );
+
+		CRS_raw.col_start[ 0 ] = 0;
+
+#ifdef _H_GRB_REFERENCE_OMP_BLAS3
+		#pragma omp parallel for simd
+#endif
+		for( size_t j = 0; j <= ncols; ++j ) {
+			CCS_raw.col_start[ j ] = 0;
+		}
+
+
+		nzc = 0;
+		for( size_t i = 0; i < nrows; ++i ) {
+			mask_coors.clear();
+			for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) {
+				const auto k_col = mask_raw.row_index[ k ];
+				if( utils::interpretMask< descr, MaskType >( true, mask_raw.getValues(), k ) ) {
+					mask_coors.assign( k_col );
+				}
+			}
+			for( auto k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) {
+				const auto k_col = A_raw.row_index[ k ];
+				if( mask_coors.assigned( k_col ) ) {
+					nzc++;
+					CCS_raw.col_start[ k_col + 1 ]++;
+				}
+			}
+			CRS_raw.col_start[ i + 1 ] = nzc;
+		}
+
+		for( size_t j = 1; j < ncols; ++j ) {
+			CCS_raw.col_start[ j + 1 ] += CCS_raw.col_start[ j ];
+		}
+
+#ifdef _H_GRB_REFERENCE_OMP_BLAS3
+		#pragma omp parallel for simd
+#endif
+		for( size_t j = 0; j < ncols; ++j ) {
+			C_col_index[ j ] = 0;
+		}
+
+
+		// use previously computed CCS offset array to update CCS during the
+		// computational phase
+		nzc = 0;
+		for( size_t i = 0; i < nrows; ++i ) {
+			mask_coors.clear();
+			for( auto k = mask_raw.col_start[ i ]; k < mask_raw.col_start[ i + 1 ]; ++k ) {
+				const auto k_col = mask_raw.row_index[ k ];
+				if( utils::interpretMatrixMask< descr, MaskType >( true, mask_raw.getValues(), k ) ) {
+					mask_coors.assign( k_col );
+				}
+			}
+			for( auto k = A_raw.col_start[ i ]; k < A_raw.col_start[ i + 1 ]; ++k ) {
+				const auto k_col = A_raw.row_index[ k ];
+				if( mask_coors.assigned( k_col ) ) {
+					OutputType val = A_raw.getValue( k, *( (OutputType*) nullptr ) );
+					CRS_raw.row_index[ nzc ] = k_col;
+					CRS_raw.setValue( nzc, val );
+					const size_t CCS_index = C_col_index[ k_col ] + CCS_raw.col_start[ k_col ];
+					C_col_index[ k_col ]++;
+					CCS_raw.row_index[ CCS_index ] = i;
+					CCS_raw.setValue( CCS_index, val );
+					nzc++;
+				}
+			}
+		}
+#ifndef NDEBUG
+		for( size_t j = 0; j < ncols; ++j ) {
+			assert( CCS_raw.col_start[ j + 1 ] - CCS_raw.col_start[ j ] ==
+				C_col_index[ j ] );
+		}
+#endif
+
+		internal::setCurrentNonzeroes( C, nzc );
+
+		return SUCCESS;
+	}
+
+
+
+
+
+
 	/**
 	 * Ingests raw data into a GraphBLAS vector.
 	 *

include/graphblas/utils.hpp

@@ -380,6 +380,38 @@ namespace grb {
 				return ret;
 			}
 		}
+		/** Specialisation for void-valued matrice's masks */
+		template< Descriptor descriptor, typename MatrixDataType, typename ValuesType >
+		static bool interpretMatrixMask(
+			const bool &assigned,
+			const ValuesType * const values,
+			const size_t k,
+			typename std::enable_if<
+				!std::is_void< MatrixDataType >::value
+			>::type * = nullptr
+		) {
+			return interpretMask< descriptor, ValuesType >( assigned, values, k );
+		}
+
+		/** Specialisation for void-valued matrice's masks */
+		template< Descriptor descriptor, typename MatrixDataType, typename ValuesType >
+		static bool interpretMatrixMask(
+			const bool &assigned,
+			const ValuesType * const,
+			const size_t,
+			typename std::enable_if<
+				std::is_void< MatrixDataType >::value
+			>::type * = nullptr
+		) {
+			// set default mask to false
+			bool ret = assigned;
+			// check whether we should return the inverted value
+			if( descriptor & descriptors::invert_mask ) {
+				return !ret;
+			} else {
+				return ret;
+			}
+		}
 
 	} // namespace utils
 

src/graphblas/hyperdags/hyperdags.cpp

@@ -215,6 +215,9 @@ std::string grb::internal::hyperdags::toString(
 		case SET_MATRIX_MATRIX_INPUT2:
 			return "set( matrix, matrix, scalar )";
 
+		case SET_MATRIX_MATRIX_MASKED:
+			return "set( matrix, matrix, matrix )";
+
 		case MXM_MATRIX_MATRIX_MATRIX_MONOID:
 			return "mxm( matrix, matrix, matrix, monoid, scalar, scalar )";