Bugfix of BFS_levels

Author: byjtew

include/graphblas/algorithms/bfs.hpp

@@ -36,7 +36,7 @@
 
 #include <graphblas.hpp>
 
-// #define BFS_DEBUG
+// #define _DEBUG
 
 namespace grb {
 
@@ -52,7 +52,7 @@ namespace grb {
 				(void)name;
 				(void)os;
 
-#ifdef BFS_DEBUG
+#ifdef _DEBUG
 				if( rows > 64 || cols > 64 ) {
 					return;
 				}
@@ -85,7 +85,7 @@ namespace grb {
 				(void)name;
 				(void)os;
 
-#ifdef BFS_DEBUG
+#ifdef _DEBUG
 				if( rows > 64 || cols > 64 ) {
 					return;
 				}
@@ -113,7 +113,7 @@ namespace grb {
 			void printSparseMatrix( const Matrix< D > & mat, const std::string & name ) {
 				(void)mat;
 				(void)name;
-#ifdef BFS_DEBUG
+#ifdef _DEBUG
 				wait( mat );
 				printSparseMatrixIterator( nrows( mat ), ncols( mat ), mat.cbegin(), mat.cend(), name, std::cout );
 #endif
@@ -123,7 +123,7 @@ namespace grb {
 			void printSparseMatrix< void >( const Matrix< void > & mat, const std::string & name ) {
 				(void)mat;
 				(void)name;
-#ifdef BFS_DEBUG
+#ifdef _DEBUG
 				wait( mat );
 				printSparsePatternMatrixIterator( nrows( mat ), ncols( mat ), mat.cbegin(), mat.cend(), name, std::cout );
 #endif
@@ -133,7 +133,7 @@ namespace grb {
 			void printSparseVector( const Vector< D > & v, const std::string & name ) {
 				(void)v;
 				(void)name;
-#ifdef BFS_DEBUG
+#ifdef _DEBUG
 				if( size( v ) > 64 ) {
 					return;
 				}
@@ -165,7 +165,7 @@ namespace grb {
 			void printStdVector( const std::vector< T > & vector, const std::string & name ) {
 				(void)vector;
 				(void)name;
-#ifdef BFS_DEBUG
+#ifdef _DEBUG
 				if( vector.size() > 64 ) {
 					return;
 				}
@@ -229,7 +229,7 @@ namespace grb {
 				identities::infinity
 			>,
 			class SetFalseMonoid = Monoid<
-				operators::logical_and<bool>, 
+				operators::logical_nand< bool >, 
 				identities::logical_false
 			>
 		>
@@ -253,74 +253,74 @@ namespace grb {
 				void
 			> * const = nullptr
 		) {
+			max_level = 0;
+			explored_all = false;
+
 			RC rc = SUCCESS;
 			const size_t nvertices = nrows( A );
+			
+			// Frontier vectors
+			rc = rc ? rc : setElement( x, true, root );
 
-			{
-				// Frontier vectors
-				rc = rc ? rc : setElement( x, true, root );
-
-				utils::printSparseMatrix( A, "A" );
-				utils::printSparseVector( x, "x" );
+			utils::printSparseMatrix( A, "A" );
+			utils::printSparseVector( x, "x" );
 
-				// Output vector containing the minimum level at which each vertex is reached
-				rc = rc ? rc : setElement( levels, static_cast< T >( 0 ), root );
-				utils::printSparseVector( levels, "levels" );
+			// Output vector containing the minimum level at which each vertex is reached
+			rc = rc ? rc : setElement( levels, static_cast< T >( 0 ), root );
+			utils::printSparseVector( levels, "levels" );
 
-				// Vector of unvisited vertices
-				rc = rc ? rc : set( not_visited, true );
-				rc = rc ? rc : setElement( not_visited, false, root );
+			// Vector of unvisited vertices
+			rc = rc ? rc : set( not_visited, true );
+			rc = rc ? rc : setElement( not_visited, false, root );
 
-				size_t max_iter = max_iterations < 0 ? nvertices : max_iterations;
-				for( size_t level = 1; level <= max_iter; level++ ) {
+			size_t max_iter = max_iterations < 0 ? nvertices : max_iterations;
+			for( size_t level = 1; level <= max_iter; level++ ) {
 #ifdef _DEBUG
-					std::cout << "** Level " << level << ":" << std::endl << std::flush;
+				std::cout << "** Level " << level << ":" << std::endl << std::flush;
 #endif
-					max_level = level;
-
-					// Multiply the current frontier by the adjacency matrix
-					utils::printSparseVector( x, "x" );
-					utils::printSparseVector( not_visited, "not_visited" );
-					rc = rc ? rc : resize( y, 0UL );
-
-					rc = rc ? rc : vxm( y, not_visited, x, A, bool_semiring, Phase::RESIZE );
-					rc = rc ? rc : vxm( y, not_visited, x, A, bool_semiring, Phase::EXECUTE );
-					utils::printSparseVector( y, "y" );
+				// Multiply the current frontier by the adjacency matrix
+				utils::printSparseVector( x, "x" );
+				utils::printSparseVector( not_visited, "not_visited" );
+				rc = rc ? rc : clear( y );
+				rc = rc ? rc : vxm( y, not_visited, x, A, bool_semiring, Phase::RESIZE );
+				rc = rc ? rc : vxm( y, not_visited, x, A, bool_semiring, Phase::EXECUTE );
+				utils::printSparseVector( y, "y" );
 
-					// Update not_visited vector
-					rc = rc ? rc : foldl( not_visited, y, y, not_visited_monoid, Phase::RESIZE );
-					rc = rc ? rc : foldl( not_visited, y, y, not_visited_monoid, Phase::EXECUTE );
-					// Assign the current level to the newly discovered vertices only
+				// Update not_visited vector
+				rc = rc ? rc : foldl( not_visited, y, not_visited_monoid, Phase::RESIZE );
+				rc = rc ? rc : foldl( not_visited, y, not_visited_monoid, Phase::EXECUTE );
 
-					rc = rc ? rc : foldl( levels, y, level, min_monoid, Phase::RESIZE );
-					rc = rc ? rc : foldl( levels, y, level, min_monoid, Phase::EXECUTE );
-					utils::printSparseVector( levels, "levels" );
+				// Assign the current level to the newly discovered vertices only
+				rc = rc ? rc : foldl( levels, y, level, min_monoid, Phase::RESIZE );
+				rc = rc ? rc : foldl( levels, y, level, min_monoid, Phase::EXECUTE );
+				utils::printSparseVector( levels, "levels" );
 
-					// Check if all vertices have been discovered, equivalent of an std::all on the frontier
-					explored_all = nnz( levels ) == nvertices;
-					if( explored_all ) {
-						// If all vertices are discovered, stop
+				// Check if all vertices have been discovered, equivalent of an std::all on the frontier
+				explored_all = nnz( levels ) == nvertices;
+				if( explored_all ) {
+					max_level = level;
+					// If all vertices are discovered, stop
 #ifdef _DEBUG
-						std::cout << "Explored " << level << " levels to discover all of the "
-									<< nvertices << " vertices.\n" << std::flush;
+					std::cout << "Explored " << level << " levels to discover all of the "
+								<< nvertices << " vertices.\n" << std::flush;
 #endif
-						return rc;
-					}
-					bool can_continue = nnz( y ) > 0;
-					if( ! can_continue ) {
-						max_level = level - 1;
-						// If no new vertices are discovered, stop
+					return rc;
+				}
+				bool can_continue = nnz( y ) > 0;
+				if( !can_continue ) {
+					max_level = level - 1;
+					// If no new vertices are discovered, stop
 #ifdef _DEBUG
-						std::cout << "Explored " << level << " levels to discover "
-									<< nnz( levels ) << " vertices.\n" << std::flush;
+					std::cout << "Explored " << level << " levels to discover "
+								<< nnz( levels ) << " vertices.\n" << std::flush;
 #endif
-						break;
-					}
-
-					// Swap the frontier, avoid a copy
-					std::swap( x, y );
+					break;
 				}
+
+				// Swap the frontier, avoid a copy
+				std::swap( x, y );
 			}
+		
 
 			// Maximum number of iteration passed, not every vertex has been discovered
 #ifdef _DEBUG

include/graphblas/base/internalops.hpp

@@ -42,6 +42,114 @@ namespace grb {
 		/** Core implementations of the standard operators in #grb::operators. */
 		namespace internal {
 
+			/**
+			 * Standard negation operator.
+			 *
+			 * Assumes native availability of ! on the given data types or assumes that
+			 * the relevant operators are properly overloaded.
+			 *
+			 * @tparam Op The Operator class to negate. 
+			 * 			Requires the following typedefs:
+			 * 			- \b D1: The left-hand input domain.
+			 * 			- \b D2: The right-hand input domain.
+			 * 			- \b D3: The output domain.	
+			 * 			- \b operator_type: The internal::operator type to negate.	
+			 */
+			template< 
+				class Op,
+				enum Backend implementation = config::default_backend
+			>
+			class negate {
+				public:
+
+					/** Alias to the left-hand input data type. */
+					typedef typename Op::D1 left_type;
+
+					/** Alias to the right-hand input data type. */
+					typedef typename Op::D2 right_type;
+
+					/** Alias to the output data type. */
+					typedef typename Op::D3 result_type;
+
+					/** Whether this operator has an inplace foldl. */
+					static constexpr bool has_foldl = Op::operator_type::has_foldl;
+
+					/** Whether this operator has an inplace foldr. */
+					static constexpr bool has_foldr = Op::operator_type::has_foldr;
+
+					/**
+					 * Whether this operator is \em mathematically associative; that is,
+					 * associative when assuming equivalent data types for \a IN1, \a IN2,
+					 * and \a OUT, as well as assuming exact arithmetic, no overflows, etc.
+					 */
+					static constexpr bool is_associative = Op::operator_type::is_associative;
+
+					/**
+					 * Whether this operator is \em mathematically commutative; that is,
+					 * commutative when assuming equivalent data types for \a IN1, \a IN2,
+					 * and \a OUT, as well as assuming exact arithmetic, no overflows, etc.
+					 */
+					static constexpr bool is_commutative = Op::operator_type::is_commutative;
+
+					/**
+					 * Out-of-place application of the operator.
+					 *
+					 * @param[in]  a The left-hand side input. Must be pre-allocated and
+					 *               initialised.
+					 * @param[in]  b The right-hand side input. Must be pre-allocated and
+					 *               initialised.
+					 * @param[out] c The output. Must be pre-allocated.
+					 */
+					static void apply(
+						const left_type * __restrict__ const a,
+						const right_type * __restrict__ const b,
+						result_type * __restrict__ const c
+					) {
+						
+						Op::operator_type::apply( a, b, c );
+						*c = !*c;
+
+					}
+
+					/**
+					 * In-place left-to-right folding.
+					 *
+					 * @param[in]     a Pointer to the left-hand side input data.
+					 * @param[in,out] c Pointer to the right-hand side input data. This also
+					 *                  dubs as the output memory area.
+					 */
+					static void foldr(
+						const left_type * __restrict__ const a,
+						result_type * __restrict__ const c
+					) {
+						
+						Op::operator_type::foldr( a, c );
+						*c = !*c;
+
+					}
+
+					/**
+					 * In-place right-to-left folding.
+					 *
+					 * @param[in,out] c Pointer to the left-hand side input data. This also
+					 *                  dubs as the output memory area.
+					 * @param[in]     b Pointer to the right-hand side input data.
+					 */
+					static void foldl(
+						result_type * __restrict__ const c,
+						const right_type * __restrict__ const b
+					) {
+						
+						Op::operator_type::foldl( c, b );
+						*c = !*c;
+
+					}
+			};
+
+			template< class Op >
+			class not_op : public negate< Op > {};
+
+
 			/**
 			 * Standard argmin operator.
 			 *
@@ -4179,6 +4287,9 @@ namespace grb {
 					/** The output domain of this operator. */
 					typedef typename OperatorBase< OP >::D3 D3;
 
+					/** The type of the operator OP. */
+					typedef OP operator_type;
+
 					/**
 					 * Reduces a vector of type \a InputType into a value in \a IOType
 					 * by repeated application of this operator. The \a IOType is cast