Resolving bug in grb::dot with nonblocking backend (#398)

The implementation of the dot-product for the nonblocking backend could iterate over different nonzero indices in its input vectors. This bug did not materialise always, and did not materialise if the dot-product was called with the dense descriptor. This MR extends the unit test of the dot product to test more challenging sparse nonzero patterns, especially patterns that differ amongst the inputs, and add tests with two vectors with partial overlap as well as with zero overlap. The thus-extended unit test was able to trigger the pre-existing bug, which this MR also addresses. In addressing the bug, furthermore, issue #408 was uncovered.

Thanks to @GiovaGa for finding the bug, and for proposing an initial fix and unit test extension!

Author: GiovaGa

include/graphblas/nonblocking/blas1.hpp

@@ -10398,6 +10398,9 @@ namespace grb {
 	// internal namespace for implementation of grb::dot
 	namespace internal {
 
+		/**
+		 * Invariant: x should have fewer nonzeroes than y
+		 */
 		template<
 			Descriptor descr,
 #ifdef GRB_BOOLEAN_DISPATCHER
@@ -10432,6 +10435,7 @@ namespace grb {
 #else
 			(void) upper_bound;
 #endif
+			assert( local_x.nonzeroes() <= local_y.nonzeroes() );
 
 			// get raw alias
 			const InputType1 * __restrict__ a = internal::getRaw( x );
@@ -10450,8 +10454,7 @@ namespace grb {
 
 					// prepare registers
 					for( size_t k = 0; k < AnyOp::blocksize; ++k, ++i ) {
-						mask[ k ] = already_dense_input_x ||
-							local_x.assigned( already_dense_input_y ? i : local_y.index( i ) );
+						mask[ k ] = already_dense_input_x || local_y.assigned(local_x.index( i ));
 					}
 
 					// rewind
@@ -10461,9 +10464,9 @@ namespace grb {
 					for( size_t k = 0; k < AnyOp::blocksize; ++k, ++i ) {
 						if( mask[ k ] ) {
 							xx[ k ] = static_cast< typename AnyOp::D1 >(
-								a[ ( already_dense_input_y ? i : local_y.index( i ) ) + lower_bound ] );
+								a[ ( already_dense_input_x ? i : local_x.index( i ) ) + lower_bound ] );
 							yy[ k ] = static_cast< typename AnyOp::D2 >(
-								b[ ( already_dense_input_y ? i : local_y.index( i ) ) + lower_bound ] );
+								b[ ( already_dense_input_x ? i : local_x.index( i ) ) + lower_bound ] );
 						}
 					}
 
@@ -10504,9 +10507,9 @@ namespace grb {
 				for( ; i < local_nz; ++i ) {
 					typename AddMonoid::D3 temp =
 						addMonoid.template getIdentity< typename AddMonoid::D3 >();
-					const size_t index = ( already_dense_input_y ? i : local_y.index( i ) ) +
+					const size_t index = ( already_dense_input_x ? i : local_x.index( i ) ) +
 						lower_bound;
-					if( already_dense_input_x || local_x.assigned( index - lower_bound ) ) {
+					if( already_dense_input_y || local_y.assigned( index - lower_bound ) ) {
 						apply( temp, a[ index ], b[ index ], anyOp );
 						foldr( temp, thread_local_output, addMonoid.getOperator() );
 					}
@@ -10658,7 +10661,9 @@ namespace grb {
 										already_dense_input_y, already_dense_input_x,
 										array_reduced[ thread_id ],
 										lower_bound, upper_bound,
-										local_y, local_x, x, y, local_y_nz,
+										local_y, local_x,
+										x, y,
+										local_y_nz,
 										addMonoid, anyOp
 									);
 							}

tests/unit/dot.cpp

@@ -27,6 +27,7 @@ void grb_program( const size_t &n, grb::RC &rc ) {
 
 	// repeatedly used containers
 	grb::Vector< bool > even_mask( n );
+	grb::Vector< bool > odd_mask( n );
 	grb::Vector< size_t > temp( n );
 	grb::Vector< double > left( n );
 	grb::Vector< double > right( n );
@@ -42,7 +43,14 @@ void grb_program( const size_t &n, grb::RC &rc ) {
 		}, temp );
 	rc = rc ? rc : grb::set( even_mask, temp, true );
 	if( rc != grb::SUCCESS ) {
-		std::cerr << "\t initialisation of mask FAILED\n";
+		std::cerr << "\t initialisation of even mask FAILED\n";
+		return;
+	}
+
+	// create odd mask
+	rc = grb::set< grb::descriptors::invert_mask >( odd_mask, temp, true );
+	if( rc != grb::SUCCESS ) {
+		std::cerr << "\t initialisation of odd mask FAILED\n";
 		return;
 	}
 
@@ -160,6 +168,75 @@ void grb_program( const size_t &n, grb::RC &rc ) {
 		return;
 	}
 
+	// test 5, init
+	rc = rc ? rc :  grb::clear( x );
+	rc = rc ? rc :  grb::clear( y );
+	assert( rc == grb::SUCCESS );
+	int true_dot = 0;
+	alpha = 0;
+
+	rc = rc ? rc : grb::set< grb::descriptors::use_index >( y, 0 );
+	for( size_t i : {2,3,5,7,13,17,19,23,29} ){
+		if( i >= n ) break;
+		rc = rc ? rc :  grb::setElement( x, 1, i );
+		true_dot += i;
+	}
+	assert( rc == grb::SUCCESS );
+
+	// test 5, exec
+	rc = grb::dot( alpha, x, y, intRing );
+	if( rc != SUCCESS ) {
+		std::cerr << "\t test 5 (sparse non constant-value vectors) dot FAILED\n";
+		return;
+	}
+
+	// test 5, check
+	if( alpha != true_dot ) {
+		std::cerr << "\t test 5 (sparse non constant-value vectors) unexpected value "
+			<< alpha << ", expected " << true_dot << ".\n";
+		rc = FAILED;
+		return;
+	}
+
+	// test 6, init
+	rc = grb::set( y, x );
+	rc = rc ? rc : grb::set< grb::descriptors::use_index >( x, 0 );
+	alpha = 0;
+
+	// test 6, exec
+	rc = rc ? rc : grb::dot( alpha, x, y, intRing );
+	if( rc != SUCCESS ) {
+		std::cerr << "\t test 6 (swapped version of test 5) dot FAILED\n";
+		return;
+	}
+
+	// test 6, check
+	if( alpha != true_dot ) {
+		std::cerr << "\t test 6 (swapped version of test 5) unexpected value "
+			<< alpha << ", expected " << true_dot << ".\n";
+		rc = FAILED;
+		return;
+	}
+
+	// test 7, init
+	rc = grb::set< grb::descriptors::use_index >( x, odd_mask, 0 );
+	rc = rc ? rc : grb::set< grb::descriptors::use_index >( y, even_mask, 0 );
+	alpha = 0;
+
+	// test 7, exec
+	rc = rc ? rc : grb::dot( alpha, x, y, intRing );
+	if( rc != SUCCESS ) {
+		std::cerr << "\t test 7 (sparse dot no overlap) FAILED\n";
+		return;
+	}
+
+	// test 7, check
+	if( alpha != 0 ) {
+		std::cerr << "\t test 7 (sparse dot no overlap) unexpected value "
+			<< alpha << ", expected 0.\n";
+		rc = FAILED;
+		return;
+	}
 }
 
 int main( int argc, char ** argv ) {