diff --git a/include/graphblas/base/blas3.hpp b/include/graphblas/base/blas3.hpp
index 2aab1be2a..425f7bc7a 100644
--- a/include/graphblas/base/blas3.hpp
+++ b/include/graphblas/base/blas3.hpp
@@ -87,14 +87,16 @@ namespace grb {
 	template<
 		Descriptor descr = descriptors::no_operation,
 		typename OutputType, typename InputType1, typename InputType2,
-		typename CIT, typename RIT, typename NIT,
+		typename CIT1, typename RIT1, typename NIT1,
+		typename CIT2, typename RIT2, typename NIT2,
+		typename CIT3, typename RIT3, typename NIT3,
 		class Semiring,
 		Backend backend
 	>
 	RC mxm(
-		Matrix< OutputType, backend, CIT, RIT, NIT > &C,
-		const Matrix< InputType1, backend, CIT, RIT, NIT > &A,
-		const Matrix< InputType2, backend, CIT, RIT, NIT > &B,
+		Matrix< OutputType, backend, CIT1, RIT1, NIT1 > &C,
+		const Matrix< InputType1, backend, CIT2, RIT2, NIT2 > &A,
+		const Matrix< InputType2, backend, CIT3, RIT3, NIT3 > &B,
 		const Semiring &ring = Semiring(),
 		const Phase &phase = EXECUTE
 	) {
@@ -242,6 +244,204 @@ namespace grb {
 		return ret == SUCCESS ? UNSUPPORTED : ret;
 	}
 
+	/**
+	 * Computes \f$ C = A \odot B \f$, out of place, monoid variant.
+	 *
+	 * Calculates the element-wise operation on one scalar to elements of one
+	 * matrix, \f$ C = A \odot B \f$, using the given monoid's operator. The input
+	 * and output matrices must be of same dimension.
+	 *
+	 * Any old entries of \a C will be removed after a successful call to this
+	 * primitive; that is, this is an out-of-place primitive.
+	 *
+	 * After a successful call to this primitive, the nonzero structure of \a C
+	 * will match that of the union of \a A and \a B. An implementing backend may
+	 * skip processing rows \a i and columns \a j that are not in the union of the
+	 * nonzero structure of \a A and \a B.
+	 *
+	 * \note When applying element-wise operators on sparse matrices using
+	 *       semirings, there is a difference between interpreting missing
+	 *       values as an annihilating identity or as a neutral identity--
+	 *       intuitively, such identities are known as `zero' or `one',
+	 *       respectively. As a consequence, this functionality is provided by
+	 *       #grb::eWiseApply depending on whether a monoid or operator is
+	 *       provided:
+	 *        - #grb::eWiseApply using monoids (neutral),
+	 *        - #grb::eWiseApply using operators (annihilating).
+	 *
+	 * @tparam descr      The descriptor to be used. Optional; the default is
+	 *                    #grb::descriptors::no_operation.
+	 * @tparam Monoid     The monoid to use.
+	 * @tparam InputType1 The value type of the left-hand matrix.
+	 * @tparam InputType2 The value type of the right-hand matrix.
+	 * @tparam OutputType The value type of the ouput matrix.
+	 *
+	 * @param[out]  C    The output matrix.
+	 * @param[in]   A    The left-hand input matrix.
+	 * @param[in]   B    The right-hand input matrix.
+	 * @param[in] monoid The monoid structure containing \f$ \odot \f$.
+	 * @param[in] phase  The #grb::Phase the call should execute. Optional; the
+	 *                   default parameter is #grb::EXECUTE.
+	 *
+	 * @return #grb::SUCCESS  On successful completion of this call.
+	 * @return #grb::MISMATCH Whenever the dimensions of \a x, \a y and \a z do
+	 *                        not match. All input data containers are left
+	 *                        untouched if this exit code is returned; it will be
+	 *                        as though this call was never made.
+	 * @return #grb::FAILED   If \a phase is #grb::EXECUTE, indicates that the
+	 *                        capacity of \a z was insufficient. The output matrix
+	 *                        \a z 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 z is retained.
+	 * @return #grb::PANIC    A general unmitigable error has been encountered. If
+	 *                        returned, ALP enters an undefined state and the user
+	 *                        program is encouraged to exit as quickly as possible.
+	 *
+	 * \par Performance semantics
+	 *
+	 * Each backend must define performance semantics for this primitive.
+	 *
+	 * @see perfSemantics
+	 */
+	template<
+		Descriptor descr = descriptors::no_operation,
+		class Monoid,
+		typename OutputType, typename InputType1, typename InputType2,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2,
+		typename RIT3, typename CIT3, typename NIT3,
+		Backend backend
+	>
+	RC eWiseApply(
+		Matrix< OutputType, backend, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType1, backend, RIT2, CIT2, NIT2 > &A,
+		const Matrix< InputType2, backend, RIT3, CIT3, NIT3 > &B,
+		const Monoid &monoid,
+		const Phase phase = EXECUTE,
+		const typename std::enable_if<
+			!grb::is_object< OutputType >::value &&
+			!grb::is_object< InputType1 >::value &&
+			!grb::is_object< InputType2 >::value &&
+			grb::is_monoid< Monoid >::value,
+		void >::type * const = nullptr
+	) {
+		(void) C;
+		(void) A;
+		(void) B;
+		(void) phase;
+#ifdef _DEBUG
+		std::cerr << "Selected backend does not implement grb::eWiseApply\n";
+#endif
+#ifndef NDEBUG
+		const bool selected_backend_does_not_support_ewiseapply = false;
+		assert( selected_backend_does_not_support_ewiseapply );
+#endif
+		const RC ret = grb::clear( A );
+		return ret == SUCCESS ? UNSUPPORTED : ret;
+	}
+
+	/**
+	 * Computes \f$ C = A \odot B \f$, out of place, operator variant.
+	 *
+	 * Calculates the element-wise operation on one scalar to elements of one
+	 * matrix, \f$ C = A \odot B \f$, using the given operator. The input and
+	 * output matrices must be of same dimension.
+	 *
+	 * Any old entries of \a C will be removed after a successful call to this
+	 * primitive; that is, this primitive is out-of-place.
+	 *
+	 * After a successful call to this primitive, the nonzero structure of \a C
+	 * will match that of the intersection of \a A and \a B. An implementing
+	 * backend may skip processing rows \a i and columns \a j that are not in the
+	 * intersection of the nonzero structure of \a A and \a B.
+	 *
+	 * \note When applying element-wise operators on sparse matrices using
+	 *       semirings, there is a difference between interpreting missing
+	 *       values as an annihilating identity or as a neutral identity--
+	 *       intuitively, such identities are known as `zero' or `one',
+	 *       respectively. As a consequence, this functionality is provided by
+	 *       #grb::eWiseApply depending on whether a monoid or operator is
+	 *       provided:
+	 *        - #grb::eWiseApply using monoids (neutral),
+	 *        - #grb::eWiseApply using operators (annihilating).
+	 *
+	 * @tparam descr      The descriptor to be used. Optional; the default is
+	 *                    #grb::descriptors::no_operation.
+	 * @tparam Operator   The operator to use.
+	 * @tparam InputType1 The value type of the left-hand matrix.
+	 * @tparam InputType2 The value type of the right-hand matrix.
+	 * @tparam OutputType The value type of the ouput matrix.
+	 *
+	 * @param[out]  C      The output matrix.
+	 * @param[in]   A      The left-hand input matrix.
+	 * @param[in]   B      The right-hand input matrix.
+	 * @param[in]   op     The operator.
+	 * @param[in]   phase  The #grb::Phase the call should execute. Optional; the
+	 *                     default parameter is #grb::EXECUTE.
+	 *
+	 * @return #grb::SUCCESS  On successful completion of this call.
+	 * @return #grb::MISMATCH Whenever the dimensions of \a x, \a y and \a z 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 z was insufficient. The output
+	 *                        matrix \a z 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 z is retained.
+	 * @return #grb::PANIC    A general unmitigable error has been encountered. If
+	 *                        returned, ALP enters an undefined state and the user
+	 *                        program is encouraged to exit as quickly as possible.
+	 *
+	 * \par Performance semantics
+	 *
+	 * Each backend must define performance semantics for this primitive.
+	 *
+	 * @see perfSemantics
+	 */
+	template<
+		Descriptor descr = grb::descriptors::no_operation,
+		class Operator,
+		typename OutputType, typename InputType1, typename InputType2,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2,
+		typename RIT3, typename CIT3, typename NIT3,
+		Backend backend
+	>
+	RC eWiseApply(
+		Matrix< OutputType, backend, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType1, backend, RIT2, CIT2, NIT2 > &A,
+		const Matrix< InputType2, backend, RIT3, CIT3, NIT3 > &B,
+		const Operator &op,
+		const Phase phase = EXECUTE,
+		const typename std::enable_if<
+			!grb::is_object< OutputType >::value &&
+			!grb::is_object< InputType1 >::value &&
+			!grb::is_object< InputType2 >::value &&
+			grb::is_operator< Operator >::value,
+		void >::type * const = nullptr
+	) {
+		(void) C;
+		(void) A;
+		(void) B;
+		(void) phase;
+#ifdef _DEBUG
+		std::cerr << "Selected backend does not implement grb::eWiseApply\n";
+#endif
+#ifndef NDEBUG
+		const bool selected_backend_does_not_support_ewiseapply	= false;
+		assert( selected_backend_does_not_support_ewiseapply );
+#endif
+		const RC ret = grb::clear( A );
+		return ret == SUCCESS ? UNSUPPORTED : ret;
+	}
+
 	/**
 	 * @}
 	 */
diff --git a/include/graphblas/base/io.hpp b/include/graphblas/base/io.hpp
index c0ed7e1cc..120721feb 100644
--- a/include/graphblas/base/io.hpp
+++ b/include/graphblas/base/io.hpp
@@ -784,7 +784,8 @@ namespace grb {
 		typename RIT, typename CIT, typename NIT
 	>
 	RC resize(
-		Matrix< InputType, backend, RIT, CIT, NIT > &A, const size_t new_nz
+		Matrix< InputType, backend, RIT, CIT, NIT > &A,
+		const size_t new_nz
 	) noexcept {
 #ifndef NDEBUG
 		const bool should_not_call_base_matrix_resize = false;
@@ -1327,14 +1328,14 @@ namespace grb {
 	 */
 	template<
 		Descriptor descr = descriptors::no_operation,
-		typename InputType,
+		typename InputType, typename RIT, typename CIT, typename NIT,
 		typename fwd_iterator1 = const size_t * __restrict__,
 		typename fwd_iterator2 = const size_t * __restrict__,
 		typename fwd_iterator3 = const InputType * __restrict__,
 		Backend implementation = config::default_backend
 	>
 	RC buildMatrixUnique(
-		Matrix< InputType, implementation > &A,
+		Matrix< InputType, implementation, RIT, CIT, NIT > &A,
 		fwd_iterator1 I, const fwd_iterator1 I_end,
 		fwd_iterator2 J, const fwd_iterator2 J_end,
 		fwd_iterator3 V, const fwd_iterator3 V_end,
@@ -1359,14 +1360,14 @@ namespace grb {
 	 */
 	template<
 		Descriptor descr = descriptors::no_operation,
-		typename InputType,
+		typename InputType, typename RIT, typename CIT, typename NIT,
 		typename fwd_iterator1 = const size_t * __restrict__,
 		typename fwd_iterator2 = const size_t * __restrict__,
 		typename fwd_iterator3 = const InputType * __restrict__,
 		Backend implementation = config::default_backend
 	>
 	RC buildMatrixUnique(
-		Matrix< InputType, implementation > &A,
+		Matrix< InputType, implementation, RIT, CIT, NIT > &A,
 		fwd_iterator1 I, fwd_iterator2 J, fwd_iterator3 V,
 		const size_t nz, const IOMode mode
 	) {
diff --git a/include/graphblas/bsp1d/blas2.hpp b/include/graphblas/bsp1d/blas2.hpp
index 42c5875d9..9cc1da305 100644
--- a/include/graphblas/bsp1d/blas2.hpp
+++ b/include/graphblas/bsp1d/blas2.hpp
@@ -449,12 +449,13 @@ namespace grb {
 		typename InputType2 = typename Ring::D2,
 		typename InputType3 = bool,
 		typename InputType4 = bool,
+		typename RIT, typename CIT, typename NIT,
 		typename Coords
 	>
 	RC mxv(
 		Vector< IOType, BSP1D, Coords > &u,
 		const Vector< InputType3, BSP1D, Coords > &mask,
-		const Matrix< InputType2, BSP1D > &A,
+		const Matrix< InputType2, BSP1D, RIT, CIT, NIT > &A,
 		const Vector< InputType1, BSP1D, Coords > &v,
 		const Ring &ring = Ring(),
 		const Phase &phase = EXECUTE,
diff --git a/include/graphblas/bsp1d/blas3.hpp b/include/graphblas/bsp1d/blas3.hpp
index 091128cf3..386beb164 100644
--- a/include/graphblas/bsp1d/blas3.hpp
+++ b/include/graphblas/bsp1d/blas3.hpp
@@ -40,9 +40,12 @@ namespace grb {
 		 * cleared.
 		 * \endinternal
 		 */
-		template< typename DataType, Backend backend >
+		template<
+			typename DataType, Backend backend,
+			typename RIT, typename CIT, typename NIT
+		>
 		RC checkGlobalErrorStateOrClear(
-			Matrix< DataType, backend > &A,
+			Matrix< DataType, backend, RIT, CIT, NIT > &A,
 			const RC local_rc
 		) noexcept {
 			RC global_rc = local_rc;
@@ -67,11 +70,15 @@ namespace grb {
 	// of the use of the above internal convenience function
 
 	/** \internal No implementation details; simply delegates */
-	template< Descriptor descr = descriptors::no_operation,
-		typename DataType1, typename DataType2 >
+	template<
+		Descriptor descr = descriptors::no_operation,
+		typename DataType1, typename DataType2,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2
+	>
 	RC set(
-		Matrix< DataType1, BSP1D > &out,
-		const Matrix< DataType2, BSP1D > &in,
+		Matrix< DataType1, BSP1D, RIT1, CIT1, NIT1 > &out,
+		const Matrix< DataType2, BSP1D, RIT2, CIT2, NIT2 > &in,
 		const Phase &phase = EXECUTE
 	) noexcept {
 		assert( phase != TRY );
@@ -88,11 +95,13 @@ namespace grb {
 	/** \internal Simply delegates to process-local backend. */
 	template<
 		Descriptor descr = descriptors::no_operation,
-		typename DataType1, typename DataType2, typename DataType3
+		typename DataType1, typename DataType2, typename DataType3,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2
 	>
 	RC set(
-		Matrix< DataType1, BSP1D > &out,
-		const Matrix< DataType2, BSP1D > &mask,
+		Matrix< DataType1, BSP1D, RIT1, CIT1, NIT1 > &out,
+		const Matrix< DataType2, BSP1D, RIT2, CIT2, NIT2 > &mask,
 		const DataType3 &val,
 		const Phase &phase = EXECUTE
 	) noexcept {
@@ -111,15 +120,20 @@ namespace grb {
 	/** \internal Simply delegates to process-local backend */
 	template<
 		Descriptor descr = descriptors::no_operation,
+		class MulMonoid,
 		typename OutputType, typename InputType1, typename InputType2,
-		class MulMonoid
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2,
+		typename RIT3, typename CIT3, typename NIT3
 	>
-	RC eWiseApply( Matrix< OutputType, BSP1D > &C,
-		const Matrix< InputType1, BSP1D > &A,
-		const Matrix< InputType2, BSP1D > &B,
+	RC eWiseApply(
+		Matrix< OutputType, BSP1D, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType1, BSP1D, RIT2, CIT2, NIT2 > &A,
+		const Matrix< InputType2, BSP1D, RIT3, CIT3, NIT3 > &B,
 		const MulMonoid &mul,
 		const Phase phase = EXECUTE,
-		const typename std::enable_if< !grb::is_object< OutputType >::value &&
+		const typename std::enable_if<
+			!grb::is_object< OutputType >::value &&
 			!grb::is_object< InputType1 >::value &&
 			!grb::is_object< InputType2 >::value &&
 			grb::is_monoid< MulMonoid >::value,
@@ -154,17 +168,22 @@ namespace grb {
 	/** \internal Simply delegates to process-local backend */
 	template<
 		Descriptor descr = descriptors::no_operation,
+		class Operator,
 		typename OutputType, typename InputType1, typename InputType2,
-		class Operator
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2,
+		typename RIT3, typename CIT3, typename NIT3
 	>
-	RC eWiseApply( Matrix< OutputType, BSP1D > &C,
-		const Matrix< InputType1, BSP1D > &A,
-		const Matrix< InputType2, BSP1D > &B,
+	RC eWiseApply(
+		Matrix< OutputType, BSP1D, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType1, BSP1D, RIT2, CIT2, NIT2 > &A,
+		const Matrix< InputType2, BSP1D, RIT3, CIT3, NIT3 > &B,
 		const Operator &op,
 		const Phase phase = EXECUTE,
-		const typename std::enable_if< !grb::is_object< OutputType >::value &&
-			!grb::is_object< InputType1 >::value && !
-			grb::is_object< InputType2 >::value &&
+		const typename std::enable_if<
+			!grb::is_object< OutputType >::value &&
+			!grb::is_object< InputType1 >::value &&
+			!grb::is_object< InputType2 >::value &&
 			grb::is_operator< Operator >::value,
 		void >::type * const = nullptr
 	) {
diff --git a/include/graphblas/hyperdags/blas2.hpp b/include/graphblas/hyperdags/blas2.hpp
index e2e781be2..421b0b646 100644
--- a/include/graphblas/hyperdags/blas2.hpp
+++ b/include/graphblas/hyperdags/blas2.hpp
@@ -39,13 +39,14 @@ namespace grb {
 	template<
 		Descriptor descr = descriptors::no_operation, class Ring,
 		typename IOType, typename InputType1, typename InputType2,
-		typename InputType3, typename Coords
+		typename InputType3, typename Coords,
+		typename RIT, typename CIT, typename NIT
 	>
 	RC vxm(
 		Vector< IOType, hyperdags, Coords > &u,
 		const Vector< InputType3, hyperdags, Coords > &mask,
 		const Vector< InputType1, hyperdags, Coords > &v,
-		const Matrix< InputType2, hyperdags > &A,
+		const Matrix< InputType2, hyperdags, RIT, CIT, NIT > &A,
 		const Ring &ring = Ring(),
 		const Phase &phase = EXECUTE,
 		const typename std::enable_if<
@@ -87,13 +88,14 @@ namespace grb {
 		Descriptor descr = descriptors::no_operation,
 		class AdditiveMonoid, class MultiplicativeOperator,
 		typename IOType, typename InputType1, typename InputType2,
-		typename InputType3, typename Coords
+		typename InputType3, typename Coords,
+		typename RIT, typename CIT, typename NIT
 	>
 	RC vxm(
 		Vector< IOType, hyperdags, Coords > &u,
 		const Vector< InputType3, hyperdags, Coords > &mask,
 		const Vector< InputType1, hyperdags, Coords > &v,
-		const Matrix< InputType2, hyperdags > &A,
+		const Matrix< InputType2, hyperdags, RIT, CIT, NIT > &A,
 		const AdditiveMonoid &add = AdditiveMonoid(),
 		const MultiplicativeOperator &mul = MultiplicativeOperator(),
 		const Phase &phase = EXECUTE,
@@ -140,12 +142,13 @@ namespace grb {
 		typename IOType = typename Ring::D4,
 		typename InputType1 = typename Ring::D1,
 		typename InputType2 = typename Ring::D2,
-		typename Coords
+		typename Coords,
+		typename RIT, typename CIT, typename NIT
 	>
 	RC vxm(
 		Vector< IOType, hyperdags, Coords > &u,
 		const Vector< InputType1, hyperdags, Coords > &v,
-		const Matrix< InputType2, hyperdags > &A,
+		const Matrix< InputType2, hyperdags, RIT, CIT, NIT > &A,
 		const Ring &ring = Ring(),
 		const Phase &phase = EXECUTE,
 		const typename std::enable_if<
@@ -186,12 +189,13 @@ namespace grb {
 		typename InputType1 = typename Ring::D1,
 		typename InputType2 = typename Ring::D2,
 		typename InputType3 = bool,
-		typename Coords
+		typename Coords,
+		typename RIT, typename CIT, typename NIT
 	>
 	RC mxv(
 		Vector< IOType, hyperdags, Coords > &u,
 		const Vector< InputType3, hyperdags, Coords > &mask,
-		const Matrix< InputType2, hyperdags > &A,
+		const Matrix< InputType2, hyperdags, RIT, CIT, NIT > &A,
 		const Vector< InputType1, hyperdags, Coords > &v,
 		const Ring &ring,
 		const Phase &phase = EXECUTE,
@@ -237,12 +241,13 @@ namespace grb {
 		bool input_may_be_masked = true,
 		class Ring,
 		typename IOType, typename InputType1, typename InputType2,
-		typename InputType3, typename InputType4, typename Coords
+		typename InputType3, typename InputType4, typename Coords,
+		typename RIT, typename CIT, typename NIT
 	>
 	RC mxv(
 		Vector< IOType, hyperdags, Coords > &u,
 		const Vector< InputType3, hyperdags, Coords > &mask,
-		const Matrix< InputType2, hyperdags > &A,
+		const Matrix< InputType2, hyperdags, RIT, CIT, NIT > &A,
 		const Vector< InputType1, hyperdags, Coords > &v,
 		const Vector< InputType4, hyperdags, Coords > &v_mask,
 		const Ring &ring,
@@ -292,12 +297,13 @@ namespace grb {
 		bool input_may_be_masked = true,
 		class AdditiveMonoid, class MultiplicativeOperator,
 		typename IOType, typename InputType1, typename InputType2,
-		typename InputType3, typename InputType4, typename Coords
+		typename InputType3, typename InputType4, typename Coords,
+		typename RIT, typename CIT, typename NIT
 	>
 	RC mxv(
 		Vector< IOType, hyperdags, Coords > &u,
 		const Vector< InputType3, hyperdags, Coords > &mask,
-		const Matrix< InputType2, hyperdags > &A,
+		const Matrix< InputType2, hyperdags, RIT, CIT, NIT > &A,
 		const Vector< InputType1, hyperdags, Coords > &v,
 		const Vector< InputType4, hyperdags, Coords > &v_mask,
 		const AdditiveMonoid &add = AdditiveMonoid(),
@@ -350,11 +356,12 @@ namespace grb {
 		typename IOType = typename Ring::D4,
 		typename InputType1 = typename Ring::D1,
 		typename InputType2 = typename Ring::D2,
-		typename Coords
+		typename Coords,
+		typename RIT, typename CIT, typename NIT
 	>
 	RC mxv(
 		Vector< IOType, hyperdags, Coords > &u,
-		const Matrix< InputType2, hyperdags > &A,
+		const Matrix< InputType2, hyperdags, RIT, CIT, NIT > &A,
 		const Vector< InputType1, hyperdags, Coords > &v,
 		const Ring &ring,
 		const Phase &phase = EXECUTE,
@@ -392,11 +399,12 @@ namespace grb {
 	template<
 		Descriptor descr = descriptors::no_operation,
 		class AdditiveMonoid, class MultiplicativeOperator,
-		typename IOType, typename InputType1, typename InputType2, typename Coords
+		typename IOType, typename InputType1, typename InputType2, typename Coords,
+		typename RIT, typename CIT, typename NIT
 	>
 	RC mxv(
 		Vector< IOType, hyperdags, Coords > &u,
-		const Matrix< InputType2, hyperdags > &A,
+		const Matrix< InputType2, hyperdags, RIT, CIT, NIT > &A,
 		const Vector< InputType1, hyperdags, Coords > &v,
 		const AdditiveMonoid &add = AdditiveMonoid(),
 		const MultiplicativeOperator &mul = MultiplicativeOperator(),
@@ -436,11 +444,12 @@ namespace grb {
 
 	/** \internal Uses a direct implementation. */
 	template<
-		typename Func, typename DataType
+		typename Func, typename DataType,
+		typename RIT, typename CIT, typename NIT
 	>
 	RC eWiseLambda(
 		const Func f,
-		const Matrix< DataType, hyperdags > &A
+		const Matrix< DataType, hyperdags, RIT, CIT, NIT > &A
 	) {
 		const RC ret = eWiseLambda( f, internal::getMatrix(A) );
 		if( ret != SUCCESS ) { return ret; }
@@ -461,11 +470,12 @@ namespace grb {
 
 		/** \internal This is the end recursion */
 		template<
-			typename Func, typename DataType
+			typename Func, typename DataType,
+			typename RIT, typename CIT, typename NIT
 		>
 		RC hyperdag_ewisematrix(
 			const Func f,
-			const Matrix< DataType, grb::hyperdags > &A,
+			const Matrix< DataType, grb::hyperdags, RIT, CIT, NIT > &A,
 			std::vector< uintptr_t > &sources,
 			std::vector< uintptr_t > &destinations
 		) {
@@ -486,11 +496,13 @@ namespace grb {
 		/** \internal This is the base recursion */
 		template<
 			typename Func, typename DataType1, typename DataType2,
-			typename Coords, typename... Args
+			typename Coords,
+			typename RIT, typename CIT, typename NIT,
+			typename... Args
 		>
 		RC hyperdag_ewisematrix(
 			const Func f,
-			const Matrix< DataType1, grb::hyperdags > &A,
+			const Matrix< DataType1, grb::hyperdags, RIT, CIT, NIT > &A,
 			std::vector< uintptr_t > &sources,
 			std::vector< uintptr_t > &destinations,
 			const Vector< DataType2, grb::hyperdags, Coords > &x,
@@ -507,11 +519,13 @@ namespace grb {
 	template<
 		typename Func,
 		typename DataType1, typename DataType2,
-		typename Coords, typename... Args
+		typename Coords,
+		typename RIT, typename CIT, typename NIT,
+		typename... Args
 	>
 	RC eWiseLambda(
 		const Func f,
-		const Matrix< DataType1, hyperdags > &A,
+		const Matrix< DataType1, hyperdags, RIT, CIT, NIT > &A,
 		const Vector< DataType2, hyperdags, Coords > &x,
 		Args... args
 	) {
@@ -527,14 +541,15 @@ namespace grb {
 		bool input_may_be_masked = true,
 		class Ring,
 		typename IOType, typename InputType1, typename InputType2,
-		typename InputType3, typename InputType4, typename Coords
+		typename InputType3, typename InputType4, typename Coords,
+		typename RIT, typename CIT, typename NIT
 	>
 	RC vxm(
 		Vector< IOType, hyperdags, Coords > &u,
 		const Vector< InputType3, hyperdags, Coords > &mask,
 		const Vector< InputType1, hyperdags, Coords > &v,
 		const Vector< InputType4, hyperdags, Coords > &v_mask,
-		const Matrix< InputType2, hyperdags > &A,
+		const Matrix< InputType2, hyperdags, RIT, CIT, NIT > &A,
 		const Ring &ring = Ring(),
 		const Phase &phase = EXECUTE,
 		const typename std::enable_if<
@@ -583,14 +598,15 @@ namespace grb {
 		bool input_may_be_masked = true,
 		class AdditiveMonoid, class MultiplicativeOperator,
 		typename IOType, typename InputType1, typename InputType2,
-		typename InputType3, typename InputType4, typename Coords
+		typename InputType3, typename InputType4, typename Coords,
+		typename RIT, typename CIT, typename NIT
 	>
 	RC vxm(
 		Vector< IOType, hyperdags, Coords > &u,
 		const Vector< InputType3, hyperdags, Coords > &mask,
 		const Vector< InputType1, hyperdags, Coords > &v,
 		const Vector< InputType4, hyperdags, Coords > &v_mask,
-		const Matrix< InputType2, hyperdags > &A,
+		const Matrix< InputType2, hyperdags, RIT, NIT, CIT > &A,
 		const AdditiveMonoid &add = AdditiveMonoid(),
 		const MultiplicativeOperator &mul = MultiplicativeOperator(),
 		const Phase &phase = EXECUTE,
@@ -639,12 +655,13 @@ namespace grb {
 	template<
 		Descriptor descr = descriptors::no_operation,
 		class AdditiveMonoid, class MultiplicativeOperator,
-		typename IOType, typename InputType1, typename InputType2, typename Coords
+		typename IOType, typename InputType1, typename InputType2, typename Coords,
+		typename RIT, typename CIT, typename NIT
 	>
 	RC vxm(
 		Vector< IOType, hyperdags, Coords > &u,
 		const Vector< InputType1, hyperdags, Coords > &v,
-		const Matrix< InputType2, hyperdags > &A,
+		const Matrix< InputType2, hyperdags, RIT, CIT, NIT > &A,
 		const AdditiveMonoid &add = AdditiveMonoid(),
 		const MultiplicativeOperator &mul = MultiplicativeOperator(),
 		const Phase &phase = EXECUTE,
diff --git a/include/graphblas/hyperdags/blas3.hpp b/include/graphblas/hyperdags/blas3.hpp
index 9448f5f57..ee0c10f36 100644
--- a/include/graphblas/hyperdags/blas3.hpp
+++ b/include/graphblas/hyperdags/blas3.hpp
@@ -39,14 +39,16 @@ namespace grb {
 
 	template<
 		Descriptor descr = descriptors::no_operation,
+		class MulMonoid,
 		typename OutputType, typename InputType1, typename InputType2,
 		typename RIT, typename CIT, typename NIT,
-		class MulMonoid
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2
 	>
 	RC eWiseApply(
 		Matrix< OutputType, hyperdags, RIT, CIT, NIT > &C,
-		const Matrix< InputType1, hyperdags > &A,
-		const Matrix< InputType2, hyperdags > &B,
+		const Matrix< InputType1, hyperdags, RIT1, CIT1, NIT1 > &A,
+		const Matrix< InputType2, hyperdags, RIT2, CIT2, NIT2 > &B,
 		const MulMonoid &mulmono,
 		const Phase phase = EXECUTE,
 		const typename std::enable_if<
@@ -82,14 +84,16 @@ namespace grb {
 
 	template<
 		Descriptor descr = grb::descriptors::no_operation,
+		class Operator,
 		typename OutputType, typename InputType1, typename InputType2,
 		typename RIT, typename CIT, typename NIT,
-		class Operator
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2
 	>
 	RC eWiseApply(
 		Matrix< OutputType, hyperdags, RIT, CIT, NIT > &C,
-		const Matrix< InputType1, hyperdags, RIT, CIT, NIT > &A,
-		const Matrix< InputType2, hyperdags, RIT, CIT, NIT > &B,
+		const Matrix< InputType1, hyperdags, RIT1, CIT1, NIT1 > &A,
+		const Matrix< InputType2, hyperdags, RIT2, CIT2, NIT2 > &B,
 		const Operator &mulOp,
 		const Phase phase = EXECUTE,
 		const typename std::enable_if<
diff --git a/include/graphblas/hyperdags/io.hpp b/include/graphblas/hyperdags/io.hpp
index e68af3eb7..db1f09e54 100644
--- a/include/graphblas/hyperdags/io.hpp
+++ b/include/graphblas/hyperdags/io.hpp
@@ -103,10 +103,11 @@ namespace grb {
 
 	template<
 		Descriptor descr = descriptors::no_operation,
-		typename InputType, typename fwd_iterator
+		typename InputType, typename RIT, typename CIT, typename NIT,
+		typename fwd_iterator
 	>
 	RC buildMatrixUnique(
-		Matrix< InputType, hyperdags > &A,
+		Matrix< InputType, hyperdags, RIT, CIT, NIT > &A,
 		fwd_iterator start,
 		const fwd_iterator end,
 		const IOMode mode
@@ -334,11 +335,12 @@ namespace grb {
 	template<
 		Descriptor descr = descriptors::no_operation,
 		typename OutputType, typename InputType,
-		typename RIT, typename CIT, typename NIT
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2
 	>
 	RC set(
-		Matrix< OutputType, hyperdags, RIT, CIT, NIT > &C,
-		const Matrix< InputType, hyperdags, RIT, CIT, NIT > &A,
+		Matrix< OutputType, hyperdags, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType, hyperdags, RIT2, CIT2, NIT2 > &A,
 		const Phase &phase = EXECUTE
 	) {
 		const RC ret = set< descr >(
@@ -365,11 +367,12 @@ namespace grb {
 	template<
 		Descriptor descr = descriptors::no_operation,
 		typename OutputType, typename InputType1, typename InputType2,
-		typename RIT, typename CIT, typename NIT
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2
 	>
 	RC set(
-		Matrix< OutputType, hyperdags, RIT, CIT, NIT > &C,
-		const Matrix< InputType1, hyperdags, RIT, CIT, NIT > &A,
+		Matrix< OutputType, hyperdags, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType1, hyperdags, RIT2, CIT2, NIT2 > &A,
 		const InputType2 &val,
 		const Phase &phase = EXECUTE
 	) {
@@ -440,13 +443,13 @@ namespace grb {
 		return size (internal::getVector(x));
 	}
 
-	template< typename InputType >
-	size_t nrows( const Matrix< InputType, hyperdags > &A ) noexcept {
+	template< typename InputType, typename RIT, typename CIT, typename NIT >
+	size_t nrows( const Matrix< InputType, hyperdags, RIT, CIT, NIT > &A ) noexcept {
 		return nrows(internal::getMatrix(A));
 	}
 
-	template< typename InputType >
-	size_t ncols( const Matrix< InputType, hyperdags > &A ) noexcept {
+	template< typename InputType, typename RIT, typename CIT, typename NIT >
+	size_t ncols( const Matrix< InputType, hyperdags, RIT, CIT, NIT > &A ) noexcept {
 		return ncols(internal::getMatrix(A));
 	}
 
@@ -455,8 +458,8 @@ namespace grb {
 		return capacity(internal::getVector( x ));
 	}
 
-	template< typename DataType >
-	size_t capacity( const Matrix< DataType, hyperdags > &A ) noexcept {
+	template< typename DataType, typename RIT, typename CIT, typename NIT >
+	size_t capacity( const Matrix< DataType, hyperdags, RIT, CIT, NIT > &A ) noexcept {
 		return capacity(internal::getMatrix( A ));
 	}
 
@@ -465,8 +468,8 @@ namespace grb {
 		return nnz( internal::getVector( x ) );
 	}
 
-	template< typename InputType >
-	size_t nnz( const Matrix< InputType, hyperdags > &A ) noexcept {
+	template< typename InputType, typename RIT, typename CIT, typename NIT >
+	size_t nnz( const Matrix< InputType, hyperdags, RIT, CIT, NIT > &A ) noexcept {
 		return nnz(internal::getMatrix(A));
 	}
 
@@ -475,8 +478,8 @@ namespace grb {
 		return getID(internal::getVector( x ));
 	}
 
-	template< typename InputType >
-	uintptr_t getID( const Matrix< InputType, hyperdags > &A ) {
+	template< typename InputType, typename RIT, typename CIT, typename NIT >
+	uintptr_t getID( const Matrix< InputType, hyperdags, RIT, CIT, NIT > &A ) {
 		return getID(internal::getMatrix( A ));
 	}
 
@@ -506,9 +509,9 @@ namespace grb {
 		return ret;
 	}
 
-	template< typename InputType >
+	template< typename InputType, typename RIT, typename CIT, typename NIT >
 	RC resize(
-		Matrix< InputType, hyperdags > &A,
+		Matrix< InputType, hyperdags, RIT, CIT, NIT > &A,
 		const size_t new_nz
 	) noexcept {
 		const RC ret = resize( internal::getMatrix(A), new_nz );
@@ -549,9 +552,12 @@ namespace grb {
 	}
 
 	/** \internal Dispatch to base wait implementation */
-	template< typename InputType, typename... Args >
+	template<
+		typename InputType, typename RIT, typename CIT, typename NIT,
+		typename... Args
+	>
 	RC wait(
-		const Matrix< InputType, hyperdags > &A,
+		const Matrix< InputType, hyperdags, RIT, CIT, NIT > &A,
 		const Args &... args
 	) {
 		(void) A;
diff --git a/include/graphblas/nonblocking/blas3.hpp b/include/graphblas/nonblocking/blas3.hpp
index 02afce1d6..5a222c7f2 100644
--- a/include/graphblas/nonblocking/blas3.hpp
+++ b/include/graphblas/nonblocking/blas3.hpp
@@ -78,28 +78,26 @@ namespace grb {
 		template<
 			bool allow_void,
 			Descriptor descr,
-			class MulMonoid,
-			typename OutputType,
-			typename InputType1,
-			typename InputType2,
-			typename RIT,
-			typename CIT,
-			typename NIT,
+			class Monoid,
 			class Operator,
-			class Monoid
+			class MulMonoid,
+			typename OutputType, typename InputType1, typename InputType2,
+			typename RIT1, typename CIT1, typename NIT1,
+			typename RIT2, typename CIT2, typename NIT2,
+			typename RIT3, typename CIT3, typename NIT3
 		>
 		RC mxm_generic(
-			Matrix< OutputType, nonblocking, RIT, CIT, NIT > &C,
-			const Matrix< InputType1, nonblocking, RIT, CIT, NIT > &A,
-			const Matrix< InputType2, nonblocking, RIT, CIT, NIT > &B,
+			Matrix< OutputType, nonblocking, RIT1, CIT1, NIT1 > &C,
+			const Matrix< InputType1, nonblocking, RIT2, CIT2, NIT2 > &A,
+			const Matrix< InputType2, nonblocking, RIT3, CIT3, NIT3 > &B,
 			const Operator &oper,
 			const Monoid &monoid,
 			const MulMonoid &mulMonoid,
 			const Phase &phase,
 			const typename std::enable_if<
 				!grb::is_object< OutputType >::value &&
-				!grb::is_object< InputType1 >::value && !
-				grb::is_object< InputType2 >::value &&
+				!grb::is_object< InputType1 >::value &&
+				!grb::is_object< InputType2 >::value &&
 				grb::is_operator< Operator >::value &&
 				grb::is_monoid< Monoid >::value,
 			void >::type * const = nullptr
@@ -109,10 +107,7 @@ namespace grb {
 			le.execution();
 
 			// second, delegate to the reference backend
-			return mxm_generic<
-					allow_void, descr, MulMonoid, OutputType,
-					InputType1, InputType2, RIT, CIT, NIT, Operator, Monoid
-				>(
+			return mxm_generic< allow_void, descr >(
 					getRefMatrix( C ), getRefMatrix( A ), getRefMatrix( B ),
 					oper, monoid, mulMonoid, phase
 				);
@@ -122,18 +117,16 @@ namespace grb {
 
 	template<
 		Descriptor descr = descriptors::no_operation,
-		typename OutputType,
-		typename InputType1,
-		typename InputType2,
-		typename RIT,
-		typename CIT,
-		typename NIT,
-		class Semiring
+		class Semiring,
+		typename OutputType, typename InputType1, typename InputType2,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2,
+		typename RIT3, typename CIT3, typename NIT3
 	>
 	RC mxm(
-		Matrix< OutputType, nonblocking, RIT, CIT, NIT > &C,
-		const Matrix< InputType1, nonblocking, RIT, CIT, NIT > &A,
-		const Matrix< InputType2, nonblocking, RIT, CIT, NIT > &B,
+		Matrix< OutputType, nonblocking, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType1, nonblocking, RIT2, CIT2, NIT2 > &A,
+		const Matrix< InputType2, nonblocking, RIT3, CIT3, NIT3 > &B,
 		const Semiring &ring = Semiring(),
 		const Phase &phase = EXECUTE,
 		const typename std::enable_if<
@@ -183,19 +176,17 @@ namespace grb {
 
 	template<
 		Descriptor descr = grb::descriptors::no_operation,
-		typename OutputType,
-		typename InputType1,
-		typename InputType2,
-		typename RIT,
-		typename CIT,
-		typename NIT,
 		class Operator,
-		class Monoid
+		class Monoid,
+		typename OutputType, typename InputType1, typename InputType2,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2,
+		typename RIT3, typename CIT3, typename NIT3
 	>
 	RC mxm(
-		Matrix< OutputType, nonblocking, RIT, CIT, NIT > &C,
-		const Matrix< InputType1, nonblocking, RIT, CIT, NIT > &A,
-		const Matrix< InputType2, nonblocking, RIT, CIT, NIT > &B,
+		Matrix< OutputType, nonblocking, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType1, nonblocking, RIT2, CIT2, NIT2 > &A,
+		const Matrix< InputType2, nonblocking, RIT3, CIT3, NIT3 > &B,
 		const Monoid &addM,
 		const Operator &mulOp,
 		const Phase &phase = EXECUTE,
@@ -265,13 +256,12 @@ namespace grb {
 			Descriptor descr = descriptors::no_operation,
 			bool matrix_is_void,
 			typename OutputType,
-			typename InputType1,
-			typename InputType2,
-			typename InputType3,
+			typename InputType1, typename InputType2, typename InputType3,
+			typename RIT, typename CIT, typename NIT,
 			typename Coords
 		>
 		RC matrix_zip_generic(
-			Matrix< OutputType, nonblocking > &A,
+			Matrix< OutputType, nonblocking, RIT, CIT, NIT > &A,
 			const Vector< InputType1, nonblocking, Coords > &x,
 			const Vector< InputType2, nonblocking, Coords > &y,
 			const Vector< InputType3, nonblocking, Coords > &z,
@@ -291,11 +281,7 @@ namespace grb {
 			le.execution();
 
 			// second, delegate to the reference backend
-			return matrix_zip_generic<
-					descr, matrix_is_void,
-					OutputType, InputType1, InputType2, InputType3,
-					Coords
-				>(
+			return matrix_zip_generic< descr, matrix_is_void >(
 					getRefMatrix( A ), getRefVector( x ), getRefVector( y ), getRefVector( z ),
 					phase
 				);
@@ -306,13 +292,12 @@ namespace grb {
 	template<
 		Descriptor descr = descriptors::no_operation,
 		typename OutputType,
-		typename InputType1,
-		typename InputType2,
-		typename InputType3,
+		typename InputType1, typename InputType2, typename InputType3,
+		typename RIT, typename CIT, typename NIT,
 		typename Coords
 	>
 	RC zip(
-		Matrix< OutputType, nonblocking > &A,
+		Matrix< OutputType, nonblocking, RIT, CIT, NIT > &A,
 		const Vector< InputType1, nonblocking, Coords > &x,
 		const Vector< InputType2, nonblocking, Coords > &y,
 		const Vector< InputType3, nonblocking, Coords > &z,
@@ -355,12 +340,12 @@ namespace grb {
 
 	template<
 		Descriptor descr = descriptors::no_operation,
-		typename InputType1,
-		typename InputType2,
+		typename InputType1, typename InputType2,
+		typename RIT, typename CIT, typename NIT,
 		typename Coords
 	>
 	RC zip(
-		Matrix< void, nonblocking > &A,
+		Matrix< void, nonblocking, RIT, CIT, NIT > &A,
 		const Vector< InputType1, nonblocking, Coords > &x,
 		const Vector< InputType2, nonblocking, Coords > &y,
 		const Phase &phase = EXECUTE
@@ -392,25 +377,23 @@ namespace grb {
 
 	template<
 		Descriptor descr = descriptors::no_operation,
-		typename InputType1,
-		typename InputType2,
-		typename OutputType,
+		class Operator,
+		typename InputType1, typename InputType2, typename OutputType,
 		typename Coords,
-		class Operator
+		typename RIT, typename CIT, typename NIT
 	>
 	RC outer(
-		Matrix< OutputType, nonblocking > &A,
+		Matrix< OutputType, nonblocking, RIT, CIT, NIT > &A,
 		const Vector< InputType1, nonblocking, Coords > &u,
 		const Vector< InputType2, nonblocking, Coords > &v,
 		const Operator &mul = Operator(),
 		const Phase &phase = EXECUTE,
 		const typename std::enable_if<
 			grb::is_operator< Operator >::value &&
-				!grb::is_object< InputType1 >::value &&
-				!grb::is_object< InputType2 >::value &&
-				!grb::is_object< OutputType >::value,
-			void
-		>::type * const = nullptr
+			!grb::is_object< InputType1 >::value &&
+			!grb::is_object< InputType2 >::value &&
+			!grb::is_object< OutputType >::value,
+		void >::type * const = nullptr
 	) {
 		if( internal::NONBLOCKING::warn_if_not_native &&
 			config::PIPELINE::warn_if_not_native
@@ -426,13 +409,11 @@ namespace grb {
 		internal::le.execution();
 
 		// second, delegate to the reference backend
-		return outer<
-				descr, InputType1, InputType2, OutputType, Coords, Operator
-			>(
-				internal::getRefMatrix( A ),
-				internal::getRefVector( u ), internal::getRefVector( v ),
-				mul, phase
-			);
+		return outer< descr, Operator >(
+			internal::getRefMatrix( A ),
+			internal::getRefVector( u ), internal::getRefVector( v ),
+			mul, phase
+		);
 	}
 
 	namespace internal {
@@ -440,16 +421,16 @@ namespace grb {
 		template<
 			bool allow_void,
 			Descriptor descr,
-			class MulMonoid,
-			typename OutputType,
-			typename InputType1,
-			typename InputType2,
-			class Operator
+			class MulMonoid, class Operator,
+			typename OutputType, typename InputType1, typename InputType2,
+			typename RIT1, typename CIT1, typename NIT1,
+			typename RIT2, typename CIT2, typename NIT2,
+			typename RIT3, typename CIT3, typename NIT3
 		>
 		RC eWiseApply_matrix_generic(
-			Matrix< OutputType, nonblocking > &C,
-			const Matrix< InputType1, nonblocking > &A,
-			const Matrix< InputType2, nonblocking > &B,
+			Matrix< OutputType, nonblocking, RIT1, CIT1, NIT1 > &C,
+			const Matrix< InputType1, nonblocking, RIT2, CIT2, NIT2 > &A,
+			const Matrix< InputType2, nonblocking, RIT3, CIT3, NIT3 > &B,
 			const Operator &oper,
 			const MulMonoid &mulMonoid,
 			const Phase &phase,
@@ -475,7 +456,8 @@ namespace grb {
 
 			// second, delegate to the reference backend
 			return eWiseApply_matrix_generic<
-					allow_void, descr, MulMonoid, OutputType, InputType1, InputType2, Operator
+					allow_void, descr,
+					MulMonoid, Operator
 				>(
 					getRefMatrix( C ), getRefMatrix( A ), getRefMatrix( B ),
 					oper, mulMonoid, phase
@@ -486,15 +468,16 @@ namespace grb {
 
 	template<
 		Descriptor descr = descriptors::no_operation,
-		typename OutputType,
-		typename InputType1,
-		typename InputType2,
-		class MulMonoid
+		class MulMonoid,
+		typename OutputType, typename InputType1, typename InputType2,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2,
+		typename RIT3, typename CIT3, typename NIT3
 	>
 	RC eWiseApply(
-		Matrix< OutputType, nonblocking > &C,
-		const Matrix< InputType1, nonblocking > &A,
-		const Matrix< InputType2, nonblocking > &B,
+		Matrix< OutputType, nonblocking, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType1, nonblocking, RIT2, CIT2, NIT2 > &A,
+		const Matrix< InputType2, nonblocking, RIT3, CIT3, NIT3 > &B,
 		const MulMonoid &mulmono,
 		const Phase phase = EXECUTE,
 		const typename std::enable_if< !grb::is_object< OutputType >::value &&
@@ -534,15 +517,16 @@ namespace grb {
 
 	template<
 		Descriptor descr = grb::descriptors::no_operation,
-		typename OutputType,
-		typename InputType1,
-		typename InputType2,
-		class Operator
+		class Operator,
+		typename OutputType, typename InputType1, typename InputType2,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2,
+		typename RIT3, typename CIT3, typename NIT3
 	>
 	RC eWiseApply(
-		Matrix< OutputType, nonblocking > &C,
-		const Matrix< InputType1, nonblocking > &A,
-		const Matrix< InputType2, nonblocking > &B,
+		Matrix< OutputType, nonblocking, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType1, nonblocking, RIT2, CIT2, NIT2 > &A,
+		const Matrix< InputType2, nonblocking, RIT3, CIT3, NIT3 > &B,
 		const Operator &mulOp,
 		const Phase phase = EXECUTE,
 		const typename std::enable_if< !grb::is_object< OutputType >::value &&
diff --git a/include/graphblas/nonblocking/io.hpp b/include/graphblas/nonblocking/io.hpp
index 44b7f3a4d..ff40be8dd 100644
--- a/include/graphblas/nonblocking/io.hpp
+++ b/include/graphblas/nonblocking/io.hpp
@@ -1107,12 +1107,13 @@ namespace grb {
 			bool A_is_mask,
 			Descriptor descr,
 			typename OutputType,
-			typename InputType1,
-			typename InputType2 = const OutputType
+			typename InputType1, typename InputType2 = const OutputType,
+			typename RIT1, typename CIT1, typename NIT1,
+			typename RIT2, typename CIT2, typename NIT2
 		>
 		RC set(
-			Matrix< OutputType, nonblocking > &C,
-			const Matrix< InputType1, nonblocking > &A,
+			Matrix< OutputType, nonblocking, RIT1, CIT1, NIT1 > &C,
+			const Matrix< InputType1, nonblocking, RIT2, CIT2, NIT2 > &A,
 			const InputType2 * __restrict__ id = nullptr
 		) noexcept {
 			if( internal::NONBLOCKING::warn_if_not_native &&
@@ -1137,12 +1138,13 @@ namespace grb {
 
 	template<
 		Descriptor descr = descriptors::no_operation,
-		typename OutputType,
-		typename InputType
+		typename OutputType, typename InputType,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2
 	>
 	RC set(
-		Matrix< OutputType, nonblocking > &C,
-		const Matrix< InputType, nonblocking > &A,
+		Matrix< OutputType, nonblocking, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType, nonblocking, RIT2, CIT2, NIT2 > &A,
 		const Phase &phase = EXECUTE
 	) noexcept {
 		static_assert( std::is_same< OutputType, void >::value ||
@@ -1175,13 +1177,13 @@ namespace grb {
 
 	template<
 		Descriptor descr = descriptors::no_operation,
-		typename OutputType,
-		typename InputType1,
-		typename InputType2
+		typename OutputType, typename InputType1, typename InputType2,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2
 	>
 	RC set(
-		Matrix< OutputType, nonblocking > &C,
-		const Matrix< InputType1, nonblocking > &A,
+		Matrix< OutputType, nonblocking, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType1, nonblocking, RIT2, CIT2, NIT2 > &A,
 		const InputType2 &val,
 		const Phase &phase = EXECUTE
 	) noexcept {
@@ -1317,10 +1319,11 @@ namespace grb {
 	/** \internal Dispatch to base wait implementation */
 	template<
 		typename InputType,
+		typename RIT, typename CIT, typename NIT,
 		typename... Args
 	>
 	RC wait(
-		const Matrix< InputType, nonblocking > &A,
+		const Matrix< InputType, nonblocking, RIT, CIT, NIT > &A,
 		const Args &... args
 	) {
 		(void) A;
@@ -1330,7 +1333,7 @@ namespace grb {
 		return wait( args... );
 	}
 
-	template< typename InputType >
+	template< typename InputType, typename RIT, typename CIT, typename NIT >
 	RC wait( const Matrix< InputType, nonblocking > &A ) {
 		(void) A;
 		//TODO: currently, matrices are read only and no action is required
diff --git a/include/graphblas/reference/blas3.hpp b/include/graphblas/reference/blas3.hpp
index 2e330f9ea..f3f918734 100644
--- a/include/graphblas/reference/blas3.hpp
+++ b/include/graphblas/reference/blas3.hpp
@@ -67,22 +67,26 @@ namespace grb {
 		template<
 			bool allow_void,
 			Descriptor descr,
+			class Monoid,
+			class Operator,
 			class MulMonoid,
 			typename OutputType, typename InputType1, typename InputType2,
-			typename RIT, typename CIT, typename NIT,
-			class Operator, class Monoid
+			typename RIT1, typename CIT1, typename NIT1,
+			typename RIT2, typename CIT2, typename NIT2,
+			typename RIT3, typename CIT3, typename NIT3
 		>
 		RC mxm_generic(
-			Matrix< OutputType, reference, RIT, CIT, NIT > &C,
-			const Matrix< InputType1, reference, RIT, CIT, NIT > &A,
-			const Matrix< InputType2, reference, RIT, CIT, NIT > &B,
+			Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C,
+			const Matrix< InputType1, reference, RIT2, CIT2, NIT2 > &A,
+			const Matrix< InputType2, reference, RIT3, CIT3, NIT3 > &B,
 			const Operator &oper,
 			const Monoid &monoid,
 			const MulMonoid &mulMonoid,
 			const Phase &phase,
-			const typename std::enable_if< !grb::is_object< OutputType >::value &&
-				!grb::is_object< InputType1 >::value && !
-				grb::is_object< InputType2 >::value &&
+			const typename std::enable_if<
+				!grb::is_object< OutputType >::value &&
+				!grb::is_object< InputType1 >::value &&
+				!grb::is_object< InputType2 >::value &&
 				grb::is_operator< Operator >::value &&
 				grb::is_monoid< Monoid >::value,
 			void >::type * const = nullptr
@@ -338,13 +342,15 @@ namespace grb {
 	template<
 		Descriptor descr = descriptors::no_operation,
 		typename OutputType, typename InputType1, typename InputType2,
-		typename RIT, typename CIT, typename NIT,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2,
+		typename RIT3, typename CIT3, typename NIT3,
 		class Semiring
 	>
 	RC mxm(
-		Matrix< OutputType, reference, RIT, CIT, NIT > &C,
-		const Matrix< InputType1, reference, RIT, CIT, NIT > &A,
-		const Matrix< InputType2, reference, RIT, CIT, NIT > &B,
+		Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType1, reference, RIT2, CIT2, NIT2 > &A,
+		const Matrix< InputType2, reference, RIT3, CIT3, NIT3 > &B,
 		const Semiring &ring = Semiring(),
 		const Phase &phase = EXECUTE,
 		const typename std::enable_if<
@@ -390,13 +396,15 @@ namespace grb {
 	template<
 		Descriptor descr = grb::descriptors::no_operation,
 		typename OutputType, typename InputType1, typename InputType2,
-		typename RIT, typename CIT, typename NIT,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2,
+		typename RIT3, typename CIT3, typename NIT3,
 		class Operator, class Monoid
 	>
 	RC mxm(
-		Matrix< OutputType, reference, RIT, CIT, NIT > &C,
-		const Matrix< InputType1, reference, RIT, CIT, NIT > &A,
-		const Matrix< InputType2, reference, RIT, CIT, NIT > &B,
+		Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType1, reference, RIT2, CIT2, NIT2 > &A,
+		const Matrix< InputType2, reference, RIT3, CIT3, NIT3 > &B,
 		const Monoid &addM,
 		const Operator &mulOp,
 		const Phase &phase = EXECUTE,
@@ -449,7 +457,6 @@ namespace grb {
 		return internal::mxm_generic< false, descr >(
 			C, A, B, mulOp, addM, Monoid(), phase
 		);
-
 	}
 
 	namespace internal {
@@ -459,10 +466,11 @@ namespace grb {
 			bool matrix_is_void,
 			typename OutputType, typename InputType1,
 			typename InputType2, typename InputType3,
+			typename RIT, typename CIT, typename NIT,
 			typename Coords
 		>
 		RC matrix_zip_generic(
-			Matrix< OutputType, reference > &A,
+			Matrix< OutputType, reference, RIT, CIT, NIT > &A,
 			const Vector< InputType1, reference, Coords > &x,
 			const Vector< InputType2, reference, Coords > &y,
 			const Vector< InputType3, reference, Coords > &z,
@@ -709,10 +717,11 @@ namespace grb {
 		Descriptor descr = descriptors::no_operation,
 		typename OutputType, typename InputType1,
 		typename InputType2, typename InputType3,
+		typename RIT, typename CIT, typename NIT,
 		typename Coords
 	>
 	RC zip(
-		Matrix< OutputType, reference > &A,
+		Matrix< OutputType, reference, RIT, CIT, NIT > &A,
 		const Vector< InputType1, reference, Coords > &x,
 		const Vector< InputType2, reference, Coords > &y,
 		const Vector< InputType3, reference, Coords > &z,
@@ -756,10 +765,11 @@ namespace grb {
 	template<
 		Descriptor descr = descriptors::no_operation,
 		typename InputType1, typename InputType2,
+		typename RIT, typename CIT, typename NIT,
 		typename Coords
 	>
 	RC zip(
-		Matrix< void, reference > &A,
+		Matrix< void, reference, RIT, CIT, NIT > &A,
 		const Vector< InputType1, reference, Coords > &x,
 		const Vector< InputType2, reference, Coords > &y,
 		const Phase &phase = EXECUTE
@@ -799,21 +809,23 @@ namespace grb {
 	 */
 	template<
 		Descriptor descr = descriptors::no_operation,
+		class Operator,
 		typename InputType1, typename InputType2, typename OutputType,
-		typename Coords, class Operator
+		typename Coords,
+		typename RIT, typename CIT, typename NIT
 	>
-	RC outer( Matrix< OutputType, reference > &A,
+	RC outer(
+		Matrix< OutputType, reference, RIT, CIT, NIT > &A,
 		const Vector< InputType1, reference, Coords > &u,
 		const Vector< InputType2, reference, Coords > &v,
 		const Operator &mul = Operator(),
 		const Phase &phase = EXECUTE,
 		const typename std::enable_if<
 			grb::is_operator< Operator >::value &&
-				!grb::is_object< InputType1 >::value &&
-				!grb::is_object< InputType2 >::value &&
-				!grb::is_object< OutputType >::value,
-			void
-		>::type * const = nullptr
+			!grb::is_object< InputType1 >::value &&
+			!grb::is_object< InputType2 >::value &&
+			!grb::is_object< OutputType >::value,
+			void >::type * const = nullptr
 	) {
 		// static checks
 		NO_CAST_ASSERT( ( !(descr & descriptors::no_casting) ||
@@ -920,13 +932,16 @@ namespace grb {
 		template<
 			bool allow_void,
 			Descriptor descr,
-			class MulMonoid,
+			class MulMonoid, class Operator,
 			typename OutputType, typename InputType1, typename InputType2,
-			class Operator
+			typename RIT1, typename CIT1, typename NIT1,
+			typename RIT2, typename CIT2, typename NIT2,
+			typename RIT3, typename CIT3, typename NIT3
 		>
-		RC eWiseApply_matrix_generic( Matrix< OutputType, reference > &C,
-			const Matrix< InputType1, reference > &A,
-			const Matrix< InputType2, reference > &B,
+		RC eWiseApply_matrix_generic(
+			Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C,
+			const Matrix< InputType1, reference, RIT2, CIT2, NIT2 > &A,
+			const Matrix< InputType2, reference, RIT3, CIT3, NIT3 > &B,
 			const Operator &oper,
 			const MulMonoid &mulMonoid,
 			const Phase &phase,
@@ -1200,12 +1215,16 @@ namespace grb {
 	 */
 	template<
 		Descriptor descr = descriptors::no_operation,
+		class MulMonoid,
 		typename OutputType, typename InputType1, typename InputType2,
-		class MulMonoid
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2,
+		typename RIT3, typename CIT3, typename NIT3
 	>
-	RC eWiseApply( Matrix< OutputType, reference > &C,
-		const Matrix< InputType1, reference > &A,
-		const Matrix< InputType2, reference > &B,
+	RC eWiseApply(
+		Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType1, reference, RIT2, CIT2, NIT2 > &A,
+		const Matrix< InputType2, reference, RIT3, CIT3, NIT3 > &B,
 		const MulMonoid &mulmono,
 		const Phase phase = EXECUTE,
 		const typename std::enable_if< !grb::is_object< OutputType >::value &&
@@ -1253,12 +1272,16 @@ namespace grb {
 
 	template<
 		Descriptor descr = grb::descriptors::no_operation,
+		class Operator,
 		typename OutputType, typename InputType1, typename InputType2,
-		class Operator
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2,
+		typename RIT3, typename CIT3, typename NIT3
 	>
-	RC eWiseApply( Matrix< OutputType, reference > &C,
-		const Matrix< InputType1, reference > &A,
-		const Matrix< InputType2, reference > &B,
+	RC eWiseApply(
+		Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType1, reference, RIT2, CIT2, NIT2 > &A,
+		const Matrix< InputType2, reference, RIT3, CIT3, NIT3 > &B,
 		const Operator &mulOp,
 		const Phase phase = EXECUTE,
 		const typename std::enable_if< !grb::is_object< OutputType >::value &&
diff --git a/include/graphblas/reference/io.hpp b/include/graphblas/reference/io.hpp
index d482384b5..10229d1c6 100644
--- a/include/graphblas/reference/io.hpp
+++ b/include/graphblas/reference/io.hpp
@@ -262,7 +262,8 @@ namespace grb {
 	 */
 	template< typename InputType, typename RIT, typename CIT, typename NIT >
 	RC resize(
-		Matrix< InputType, reference, RIT, CIT, NIT > &A, const size_t new_nz
+		Matrix< InputType, reference, RIT, CIT, NIT > &A,
+		const size_t new_nz
 	) noexcept {
 #ifdef _DEBUG
 		std::cerr << "In grb::resize (matrix, reference)\n"
@@ -969,11 +970,13 @@ namespace grb {
 			bool A_is_mask,
 			Descriptor descr,
 			typename OutputType, typename InputType1,
-			typename InputType2 = const OutputType
+			typename InputType2 = const OutputType,
+			typename RIT1, typename CIT1, typename NIT1,
+			typename RIT2, typename CIT2, typename NIT2
 		>
 		RC set(
-			Matrix< OutputType, reference > &C,
-			const Matrix< InputType1, reference > &A,
+			Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C,
+			const Matrix< InputType1, reference, RIT2, CIT2, NIT2 > &A,
 			const InputType2 * __restrict__ id = nullptr
 		) noexcept {
 #ifdef _DEBUG
@@ -1000,7 +1003,7 @@ namespace grb {
 				return MISMATCH;
 			}
 			if( A_is_mask ) {
-				assert( id != NULL );
+				assert( id != nullptr );
 			}
 
 			// catch trivial cases
@@ -1075,11 +1078,13 @@ namespace grb {
 
 	template<
 		Descriptor descr = descriptors::no_operation,
-		typename OutputType, typename InputType
+		typename OutputType, typename InputType,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2
 	>
 	RC set(
-		Matrix< OutputType, reference > &C,
-		const Matrix< InputType, reference > &A,
+		Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType, reference, RIT2, CIT2, NIT2 > &A,
 		const Phase &phase = EXECUTE
 	) noexcept {
 		static_assert( std::is_same< OutputType, void >::value ||
@@ -1112,11 +1117,13 @@ namespace grb {
 
 	template<
 		Descriptor descr = descriptors::no_operation,
-		typename OutputType, typename InputType1, typename InputType2
+		typename OutputType, typename InputType1, typename InputType2,
+		typename RIT1, typename CIT1, typename NIT1,
+		typename RIT2, typename CIT2, typename NIT2
 	>
 	RC set(
-		Matrix< OutputType, reference > &C,
-		const Matrix< InputType1, reference > &A,
+		Matrix< OutputType, reference, RIT1, CIT1, NIT1 > &C,
+		const Matrix< InputType1, reference, RIT2, CIT2, NIT2 > &A,
 		const InputType2 &val,
 		const Phase &phase = EXECUTE
 	) noexcept {
@@ -1260,11 +1267,13 @@ namespace grb {
 	 *      processes.
 	 * \endparblock
 	 */
-	template< Descriptor descr = descriptors::no_operation,
+	template<
+		Descriptor descr = descriptors::no_operation,
 		typename InputType, typename fwd_iterator, typename Coords,
 		class Dup = operators::right_assign< InputType >
 	>
-	RC buildVector( Vector< InputType, reference, Coords > &x,
+	RC buildVector(
+		Vector< InputType, reference, Coords > &x,
 		fwd_iterator start, const fwd_iterator end,
 		const IOMode mode, const Dup & dup = Dup()
 	) {
@@ -1280,7 +1289,7 @@ namespace grb {
 #ifndef NDEBUG
 		assert( mode == SEQUENTIAL || mode == PARALLEL );
 #else
-		(void)mode;
+		(void) mode;
 #endif
 
 		// declare temporary to meet delegate signature
@@ -1445,7 +1454,8 @@ namespace grb {
 		typename InputType, typename fwd_iterator1, typename fwd_iterator2,
 		typename Coords, class Dup = operators::right_assign< InputType >
 	>
-	RC buildVector( Vector< InputType, reference, Coords > &x,
+	RC buildVector(
+		Vector< InputType, reference, Coords > &x,
 		fwd_iterator1 ind_start, const fwd_iterator1 ind_end,
 		fwd_iterator2 val_start, const fwd_iterator2 val_end,
 		const IOMode mode,
@@ -1593,9 +1603,12 @@ namespace grb {
 	}
 
 	/** \internal Dispatch to base wait implementation */
-	template< typename InputType, typename... Args >
+	template<
+		typename InputType, typename RIT, typename CIT, typename NIT,
+		 typename... Args
+	>
 	RC wait(
-		const Matrix< InputType, reference > &A,
+		const Matrix< InputType, reference, RIT, CIT, NIT > &A,
 		const Args &... args
 	) {
 		(void) A;