352 introduce a numa aware sptrsv

include/graphblas/base/blas2.hpp

@@ -337,6 +337,149 @@ namespace grb {
 		return UNSUPPORTED;
 	}
 
+	/**
+	 * Computes \f$ x \f$ from \f$ Tx=b \f$.
+	 *
+	 * \warning This is an experimental feature.
+	 *
+	 * Here, \f$ T, b \f$ are given while \f$ T \f$ additionally must be either
+	 * lower- or upper-triangular. The output \f$ x \f$ may furthermore be masked.
+	 *
+	 * @param[in,out] xb On input: the vector b. On output: the vector x.
+	 *                   When #grb::EXECUTE is given this vector must have
+	 *                   sufficient capacity to store both any pre-existing
+	 *                   nonzeroes in \f$ b \f$ plus the nonzeroes in the (masked)
+	 *                   result of \f$ x \f$.
+	 *                   When #grb::RESIZE is given, the vector contents are
+	 *                   unchanged, however its capacity may be enlarged in order
+	 *                   to ensure the above-described condition.
+	 *
+	 * \note On input, \a xb may contain both implicit and explicit zeroes.
+	 *
+	 * \note On output, \a xb may contain both implicit and explicit zeroes, even
+	 *       if on input it did not.
+	 *
+	 * @param[in] mask   The mask that acts on \f$ x \f$. This vector has either
+	 *                   size zero or size equal to that of \a x.
+	 *
+	 * \note If \a mask has size zero, it is interpreted as though the operation
+	 *       unmasked.
+	 *
+	 * @param[in] T       The upper- or lower-triangular input matrix. This must be
+	 *                    a square matrix with size equal to that of \a xb.
+	 * @param[in] forward Whether to perform forward substitution (i.e., \f$ T \f$
+	 *                    is lower-triangular), or to perform backward substitution
+	 *                    instead (i.e., \f$ T \f$ is upper-triangular).
+	 *
+	 * \note With the structural information ALP/Dense passes as template
+	 *       information, \a forward would not be required.
+	 *
+	 * @param[in] semiring    The semiring under which to perform the SpTrsv.
+	 * @param[in] subtraction The inverse of the additive operator of \a semiring.
+	 * @param[in] division    The inverse of the multiplicative operator of
+	 *                        \a semiring.
+	 *
+	 * \note That is, this operation in fact requires a field structure and not a
+	 *       semiring. A future extension to ALP may provide such a structure
+	 *       explicitly.
+	 *
+	 * @param[in] phase       The requested phase of the computation. Only
+	 *                        #grb::EXECUTE and #grb::RESIZE are supported.
+	 *
+	 * \todo Expand documentation.
+	 */
+	template<
+		Descriptor descr = descriptors::no_operation,
+		class Semiring,
+		class Subtraction,
+		class Division,
+		typename IOType, typename InputType1, typename InputType2,
+		typename Coords, typename RIT, typename CIT, typename NIT,
+		Backend backend
+	>
+	RC sptrsv(
+		Vector< IOType, backend, Coords > &xb,
+		const Vector< InputType2, backend, Coords > &mask,
+		const Matrix< InputType1, backend, RIT, CIT, NIT > &T,
+		const bool forward,
+		const Semiring &semiring = Semiring(),
+		const Subtraction &subtraction = Subtraction(),
+		const Division &division = Division(),
+		const Phase &phase = EXECUTE,
+		const typename std::enable_if<
+			grb::is_semiring< Semiring >::value &&
+			grb::is_operator< Subtraction >::value &&
+			grb::is_operator< Division >::value &&
+			!grb::is_object< IOType >::value &&
+			!grb::is_object< InputType1 >::value &&
+			!grb::is_object< InputType2 >::value,
+		void >::type * const = nullptr
+	) {
+#ifdef _DEBUG
+		std::cerr << "Selected backend does not implement masked grb::sptrsv\n";
+#endif
+#ifndef NDEBUG
+		const bool selected_backend_does_not_support_masked_sptrsv = false;
+		assert( selected_backend_does_not_support_masked_sptrsv );
+#endif
+		(void) xb;
+		(void) mask;
+		(void) T;
+		(void) semiring;
+		(void) subtraction;
+		(void) division;
+		(void) phase;
+		return UNSUPPORTED;
+	}
+
+	/**
+	 * Computes \f$ x \f$ from \f$ Tx = b \f$, unmasked variant.
+	 *
+	 * \warning This is an experimental feature.
+	 *
+	 * \todo Extend documentation
+	 */
+	template<
+		Descriptor descr = descriptors::no_operation,
+		class Semiring,
+		class Subtraction,
+		class Division,
+		typename IOType, typename InputType1,
+		typename Coords, typename RIT, typename CIT, typename NIT,
+		Backend backend
+	>
+	RC sptrsv(
+		Vector< IOType, backend, Coords > &xb,
+		const Matrix< InputType1, backend, RIT, CIT, NIT > &T,
+		const bool forward,
+		const Semiring &semiring = Semiring(),
+		const Subtraction &subtraction = Subtraction(),
+		const Division &division = Division(),
+		const Phase &phase = EXECUTE,
+		const typename std::enable_if<
+			grb::is_semiring< Semiring >::value &&
+			grb::is_operator< Subtraction >::value &&
+			grb::is_operator< Division >::value &&
+			!grb::is_object< IOType >::value &&
+			!grb::is_object< InputType1 >::value,
+		void >::type * const = nullptr
+	) {
+#ifdef _DEBUG
+		std::cerr << "Selected backend does not implement masked grb::sptrsv\n";
+#endif
+#ifndef NDEBUG
+		const bool selected_backend_does_not_support_unmasked_sptrsv = false;
+		assert( selected_backend_does_not_support_unmasked_sptrsv );
+#endif
+		(void) xb;
+		(void) T;
+		(void) semiring;
+		(void) subtraction;
+		(void) division;
+		(void) phase;
+		return UNSUPPORTED;
+	}
+
 	/**
 	 * Executes an arbitrary element-wise user-defined function \a f on all
 	 * nonzero elements of a given matrix \a A.

include/graphblas/reference/SptrsvSchedule.hpp

@@ -0,0 +1,236 @@
+
+/*
+ *   Copyright 2021 Huawei Technologies Co., Ltd.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/*
+ * @author A. N. Yzelman
+ * @date 4th of April, 2025
+ */
+
+#ifndef _H_GRB_REFERENCE_SPTRSVSCHEDULE
+#define _H_GRB_REFERENCE_SPTRSVSCHEDULE
+
+#include <limits>
+#include <vector>
+
+#include <graphblas/utils/autodeleter.hpp>
+
+
+namespace grb {
+
+	namespace internal {
+
+		template< typename NIT >
+		class SptrsvSchedule {
+
+			private:
+
+				/** One deleter for each chunk of thread-local data. */
+				std::vector< utils::AutoDeleter< char > > _deleters;
+
+				/**
+				 * Fixed-size buffer for realising the default schedule.
+				 *
+				 * The first two positions are reserved for the range of the default
+				 * schedule, while the last position is reserved for the end-position of
+				 * the default schedule.
+				 */
+				NIT default_schedule[ 3 ];
+
+				/** Computes and initialises a trivial schedule. */
+				void initTrivial( const NIT n ) {
+					// dynamic sanity check
+					assert( data[ 0 ] == nullptr );
+					// set trivial schedule
+					default_schedule[ 0 ] = 0;
+					default_schedule[ 1 ] = n;
+					default_schedule[ 2 ] = 1;
+					nRanges[ 0 ] = 1;
+					// set pointers to trivial schedule
+					data[ 0 ] = reinterpret_cast< char * >( &(default_schedule[0]) );
+					endPositions[ 0 ] = reinterpret_cast< char * >( &(default_schedule[2]) );
+					// verify trivial schedule
+#ifndef NDEBUG
+					{
+						NIT * const interpreted = reinterpret_cast< NIT * >(data[0]);
+						assert( interpreted[ 0 ] == 0 );
+						assert( interpreted[ 1 ] == n );
+					}
+					{
+						NIT * const interpreted = reinterpret_cast< NIT * >(endPositions[0]);
+						assert( *interpreted == 1 );
+					}
+#endif
+					// note that this is a trivial schedule
+					is_simple = true;
+
+					// the default schedule does not (cannot) sort input matrices
+					is_sorted = false;
+				}
+
+				/**
+				 * Implementation of move construction and assignment.
+				 */
+				void moveImpl( SptrsvSchedule &&toMove ) {
+					_deleters = std::move( toMove._deleters );
+					default_schedule = std::move( toMove.default_schedule );
+					is_sorted = toMove.is_sorted;
+					is_simple = toMove.is_simple;
+					supersteps = toMove.supersteps;
+					nThreads = toMove.nThreads;
+					data = std::move( toMove.data );
+					endPositions = std::move( toMove.endPositions );
+					nRanges = std::move( toMove.nRanges );
+					toMove.is_sorted = true;
+					toMove.is_simple = false;
+					toMove.supersteps = 0;
+					toMove.nThreads = 0;
+				}
+
+
+			public:
+
+				/** Whether the matrix is sorted as part of tuning. */
+				bool is_sorted;
+
+				/** Whether the schedule is simple.*/
+				bool is_simple;
+
+				/** Number of schedule steps. */
+				size_t supersteps;
+
+				/**
+				 * The number of threads the schedule is designed for.
+				 */
+				size_t nThreads;
+
+				/** One data pointer per thread. */
+				std::vector< char * > data;
+
+				/** One end-position array per thread. */
+				std::vector< char * > endPositions;
+
+				/** The total number of ranges per thread. */
+				std::vector< NIT > nRanges;
+
+				/** Move constructor. */
+				SptrsvSchedule( SptrsvSchedule &&toMove ) {
+					moveImpl( toMove );
+				}
+
+				/**
+				 * Base constructor.
+				 *
+				 * @param[in] n The size of the linear system.
+				 *
+				 * Initialises the schedule to the trivial always-valid schedule, which is
+				 * a single-thread single-superstep schedule that assigns all work to the
+				 * first (and only) thread.
+				 */
+				SptrsvSchedule( const NIT n ) :
+					_deleters( 1 ), is_sorted( true ), is_simple( false ),
+					supersteps( 1 ), nThreads( 1 ),
+					data( 1 ), endPositions( 1 ), nRanges( 1 )
+				{
+					data[ 0 ] = nullptr;
+					endPositions[ 0 ] = nullptr;
+					initTrivial( n );
+				}
+
+				/**
+				 * Specialised constructor for multiple threads.
+				 *
+				 * @param[in] n The size of the linear system.
+				 * @param[in] T The number of threads a schedule should be created for.
+				 *
+				 * Initialises the schedule to the trivial always-valid schedule, which is
+				 * a single-superstep schedule that assigns all work to the first thread.
+				 */
+				SptrsvSchedule( const NIT n, const size_t T ) :
+					_deleters( 2 * T ), is_sorted( true ), is_simple( false ),
+					supersteps( 1 ), nThreads( T ),
+					data( T, nullptr ), endPositions( T, nullptr ), nRanges( T, 0 )
+				{
+					if( T == 0 || T > std::numeric_limits< int >::max() ) {
+						throw std::runtime_error( "Invalid number of threads" );
+					}
+					initTrivial( n );
+				}
+
+				/**
+				 * Base destructor.
+				 *
+				 * Dynamic memory is freed through the #_deleters.
+				 */
+				~SptrsvSchedule() {
+					// sanity checks only
+					assert( nThreads != 0 );
+					if( nThreads == 1 ) {
+						assert( supersteps == 1 );
+					}
+				}
+
+				/** Move assignment. */
+				SptrsvSchedule& operator=( SptrsvSchedule &&toMove ) {
+					moveImpl( toMove );
+					return *this;
+				}
+
+				/**
+				 * Allocates a thread-local chunk of data.
+				 *
+				 * Must be called from within the thread that will use it(!)
+				 */
+				void alloc( const size_t s, const size_t nRanges_in ) {
+					assert( s < nThreads );
+					assert( supersteps > 0 );
+					assert( _deleters.size() >= s );
+					assert( data.size() >= s );
+					assert( data[ s ] == nullptr );
+					assert( endPositions[ s ] == nullptr );
+					if( nRanges[ s ] != nRanges_in ) {
+						throw std::runtime_error( "Given nRanges differs from the one given at "
+							"allocation time" );
+					}
+					grb::RC rc = grb::SUCCESS;
+					if( nRanges_in > 0 ) {
+						rc = utils::alloc(
+							"grb::internal::SptrsvSchedule (default constructor)",
+							"default thread-local data allocation, variant I",
+							data[ s ], 2 * nRanges_in * sizeof(NIT), false, _deleters[ s ],
+							endPositions[ s ], supersteps * sizeof( NIT ), false, _deleters[ s + nThreads ]
+						);
+					} else {
+						data[ s ] = nullptr;
+						rc = utils::alloc(
+							"grb::internal::SptrsvSchedule (default constructor)",
+							"default thread-local data allocation, variant II",
+							endPositions[ s ], supersteps * sizeof( NIT ), false, _deleters[ s + nThreads ]
+						);
+					}
+					if( rc != grb::SUCCESS ) {
+						throw std::bad_alloc();
+					}
+				}
+
+		};
+
+	} // end namespace `grb::internal'
+
+} // end namespace `grb'
+
+#endif // _H_GRB_REFERENCE_SPTRSVSCHEDULE
+