335 wrong fma perf test (#336)

The fma (fused multiply-add) performance test aims to replicate a simple triad benchmark, computing `z=alpha.*x+y` out-of-place.

The matching ALP/GraphBLAS primitive for that, once upon a time, was `grb::eWiseMulAdd`. Since that time, however, that primitive has undergone two changes: 1) it became in-place rather than out-of-place (i.e., computing `z+=alpha.*x+y`) while 2) since the release of the nonblocking backend this function has become deprecated (as the same could be computed using `grb::eWiseMul` followed by `grb::foldl` or `grb::foldr`, at negligible performance overheads compared to the one-shot `grb::eWiseMulAdd`).

Neither change was ever reflected in the `fma` test itself, however; and, at present, the original out-of-place `fma` computation `z=alpha.*x+y` similarly has no single matching (deprecated or not) ALP/GraphBLAS primitive. The nonblocking backend, however, could implement an fma by first calling `grb::set(z,y)` followed by the in-place `grb::eWiseMul`. 

Additionally, since the incorporation of the CMake build infra, the compiler-optimised performance tests were erroneously not compiled with performance flags, causing a further apple-to-oranges type comparison in this performance test.

The issues were detected by Alberto Scolari, who additionally indicated that adding performance tests to the CI and regularly reviewing those results would have caught this issue earlier (crossref #343 , visualisation only - the performance test now actually regularly trigger for `develop` on our internal CI). He also produced an initial set of fixes and we thank Alberto for both.

This MR fixes these issues regarding the fma performance test by:
 1. introducing the fma-blocking and fma-nonblocking shared-memory parallel performance test instead of a single fma-omp one. Only for fma-nonblocking is performance expected to match that with a compiler-optimised or lambda-generated fma;
 2. this distinction is now also printed to screen as part of the test run; and
 3. improved test initialisation checks and improved error detection.

This MR also introduced some changes to the `dot` and `reduce` performance tests, in tandem with changes to the `fma` performance test. Those "joint" changes are:
 4. the compiler-optimised shared-memory parallel fma (as well as the reduce and dot) performance tests now avoid critical sections (which is in line with standard HPC practices), and e.g. for reductions (reduce and dot) prefer the OpenMP `reduction`-clause;
 5. the compiler-optimised shared-memory parallel fma, reduce, and dot now use the same distribution of work across threads (essentially a static block-wise schedule);
 6. both sequential and shared-memory parallel compiler-optimised benchmarks are now given the proper performance flags during compilation;
 7. avoid mixing `fprintf` and `stdcerr/stdcout` in the fma, dot, and reduce performance tests as well as other more minor code style fixes;
 8. test output and test summary now distinguish between sequential and shared-memory parallel tests;
 9. use the `dense` descriptor where appropriate (potentially enhancing performance of ALP code);
 10. if any test auto-tunes the number of inner repetitions, the number is chosen so that a single experiment is expected to take no less than 100 milliseconds (down from 1 second);
 11. these tests now employ vectors of 100 million elements (up from 10 million); and, finally,
 12. it is double-checked that the numbers the new tests produce are on par with - or exceeding the performance of - the numbers published by Yzelman et al. (2020).

Author: Alberto Scolari

tests/performance/CMakeLists.txt

@@ -26,17 +26,24 @@ set( TEST_CATEGORY "performance" )
 add_library( bench_kernels OBJECT bench_kernels.c bench_kernels.h )
 add_library( bench_kernels_omp OBJECT bench_kernels.c bench_kernels.h )
 target_compile_definitions( bench_kernels_omp PRIVATE BENCH_KERNELS_OPENMP )
+target_link_libraries( bench_kernels PRIVATE test_performance_flags )
+target_link_libraries( bench_kernels_omp PRIVATE test_performance_flags OpenMP::OpenMP_C )
 
 add_grb_executables( fma fma.cpp $<TARGET_OBJECTS:bench_kernels>
 	BACKENDS reference NO_BACKEND_NAME
 	ADDITIONAL_LINK_LIBRARIES "rt"
 )
 
-add_grb_executables( fma-openmp fma.cpp $<TARGET_OBJECTS:bench_kernels_omp>
+add_grb_executables( fma-blocking fma.cpp $<TARGET_OBJECTS:bench_kernels_omp>
 	BACKENDS reference_omp NO_BACKEND_NAME
 	ADDITIONAL_LINK_LIBRARIES OpenMP::OpenMP_CXX "rt"
 )
 
+add_grb_executables( fma-nonblocking fma.cpp $<TARGET_OBJECTS:bench_kernels_omp>
+	BACKENDS nonblocking NO_BACKEND_NAME
+	ADDITIONAL_LINK_LIBRARIES OpenMP::OpenMP_CXX "rt"
+)
+
 add_grb_executables( reduce reduce.cpp $<TARGET_OBJECTS:bench_kernels>
 	BACKENDS reference NO_BACKEND_NAME
 )

tests/performance/bench_kernels.c

@@ -21,6 +21,8 @@
 
 #ifdef BENCH_KERNELS_OPENMP
 
+bool bench_kernels_parallel() { return true; }
+
 void bench_kernels_axpy(
 	double * restrict a,
 	const double alpha, const double * restrict x,
@@ -30,9 +32,25 @@ void bench_kernels_axpy(
 	assert( a != x );
 	assert( a != y );
 	assert( x != y );
-	#pragma omp parallel for schedule(static,8)
-	for( size_t i = 0; i < n; ++i ) {
-		a[ i ] = alpha * x[ i ] + y[ i ];
+	#pragma omp parallel
+	{
+		const size_t P = omp_get_num_threads();
+		const size_t s = omp_get_thread_num();
+		const size_t chunk = (n % P == 0) ? (n/P) : (n/P) + 1;
+		size_t start = chunk * s;
+		if( start > n - 1 ) {
+			start = n - 1;
+		}
+		size_t end = start + chunk;
+		if( end > n ) {
+			end = n;
+		}
+		assert( start <= end );
+		if( start != end ) {
+			for( size_t i = start; i < end; ++i ) {
+				a[ i ] = alpha * x[ i ] + y[ i ];
+			}
+		}
 	}
 }
 
@@ -45,7 +63,8 @@ void bench_kernels_dot(
 	assert( alpha != xr );
 	assert( alpha != yr );
 	*alpha = xr[ n - 1 ] * yr[ n - 1];
-	#pragma omp parallel
+	double global_alpha = 0;
+	#pragma omp parallel reduction(+:global_alpha)
 	{
 		const size_t P = omp_get_num_threads();
 		const size_t s = omp_get_thread_num();
@@ -64,20 +83,19 @@ void bench_kernels_dot(
 			for( size_t i = start; i < end - 1; ++i ) {
 				local_alpha += xr[ i ] * yr[ i ];
 			}
-			#pragma omp critical
-			{
-				*alpha += local_alpha;
-			}
+			global_alpha += local_alpha;
 		}
 	}
+	*alpha += global_alpha;
 }
 
 void bench_kernels_reduce(
 	double * restrict const alpha, const double * restrict xr, const size_t n
 ) {
 	assert( alpha != xr );
 	*alpha = xr[ n - 1 ];
-	#pragma omp parallel
+	double global_alpha = 0.0;
+	#pragma omp parallel reduction(+:global_alpha)
 	{
 		const size_t P = omp_get_num_threads();
 		const size_t s = omp_get_thread_num();
@@ -96,16 +114,16 @@ void bench_kernels_reduce(
 			for( size_t i = start; i < end - 1; ++i ) {
 				local_alpha += xr[ i ];
 			}
-			#pragma omp critical
-			{
-				*alpha += local_alpha;
-			}
+			global_alpha += local_alpha;
 		}
 	}
+	*alpha += global_alpha;
 }
 
 #else
 
+bool bench_kernels_parallel() { return false; }
+
 void bench_kernels_axpy(
 	double * restrict a,
 	const double alpha, const double * restrict x,

tests/performance/bench_kernels.h

@@ -18,7 +18,7 @@
 
 #include <omp.h>
 #include <assert.h>
-#include <stddef.h> //for size_t
+#include <stddef.h> // for size_t
 
 
 #ifdef __cplusplus
@@ -41,10 +41,14 @@ extern "C" {
 		double * __restrict__ const, const double * __restrict__, const size_t
 	);
 
+	bool bench_kernels_parallel();
+
 }
 
 #else
 
+#include <stdbool.h> // for bool
+
 /**
  * Executes \f$ a = \alpha x + y \f$ for \a a, \a x, and \a y vectors of
  * length \a n.
@@ -89,5 +93,8 @@ void bench_kernels_reduce(
 	double * restrict const alpha, const double * restrict x, const size_t n
 );
 
+/** @returns Whether the kernels defined here are (shared-memory) parallel. */
+bool bench_kernels_parallel();
+
 #endif