0003-Do-optimizaton.patch

From e6932c3be3dbe1bb791fab4aa143019a23e73c30 Mon Sep 17 00:00:00 2001
From: wangyon1 <yong4.wang@intel.com>
Date: Thu, 4 Feb 2021 17:00:16 +0800
Subject: [PATCH 3/4] Do optimizaton

---
 .../Beamformer/HilbertFirEnvelope.dp.cpp      |  53 +++++-
 src/SupraLib/Beamformer/RxBeamformerCommon.h  |   6 +
 .../Beamformer/RxBeamformerCuda.dp.cpp        | 104 +++++++++++-
 .../RxSampleBeamformerDelayAndSum.h           | 105 ++++++++++++
 src/SupraLib/Beamformer/WindowFunction.h      |   8 +-
 src/SupraLib/Beamformer/helper.h              | 159 ++++++++++++++++++
 6 files changed, 425 insertions(+), 10 deletions(-)
 create mode 100644 src/SupraLib/Beamformer/helper.h

diff --git a/src/SupraLib/Beamformer/HilbertFirEnvelope.dp.cpp b/src/SupraLib/Beamformer/HilbertFirEnvelope.dp.cpp
index af16284..99cf7be 100644
--- a/src/SupraLib/Beamformer/HilbertFirEnvelope.dp.cpp
+++ b/src/SupraLib/Beamformer/HilbertFirEnvelope.dp.cpp
@@ -14,6 +14,7 @@
 #include "HilbertFirEnvelope.h"
 #include <utilities/utility.h>
 #include <utilities/FirFilterFactory.h>
+#include "helper.h"
 
 #include <dpct/dpl_utils.hpp>
 #include <oneapi/dpl/execution>
@@ -62,6 +63,54 @@ namespace supra
 
 	}
 
+	const int H_VEC_SIZE = 4;
+	template <typename InputType, typename OutputType>
+	void vec_kernelFilterDemodulation(
+		const InputType* __restrict__ signal,
+		const HilbertFirEnvelope::WorkType * __restrict__ filter,
+		OutputType * __restrict__ out,
+		const int numSamples,
+		const int numScanlines,
+		const int filterLength,
+		sycl::nd_item<3> item_ct1) {
+
+		int scanlineIdx = item_ct1.get_local_range().get(2) * item_ct1.get_group(2) + item_ct1.get_local_id(2);
+		int sampleIdx = item_ct1.get_local_range().get(1) * item_ct1.get_group(1) + item_ct1.get_local_id(1);
+
+		scanlineIdx *= H_VEC_SIZE;
+		if (scanlineIdx < numScanlines && sampleIdx < numSamples)
+		{
+			sycl::vec<HilbertFirEnvelope::WorkType, H_VEC_SIZE> accumulator(0.0);
+				
+			int startPoint = sampleIdx - filterLength / 2;
+			int endPoint = sampleIdx + filterLength / 2;
+			int currentFilterElement = 0;
+
+			for (int currentSample = startPoint;
+				currentSample <= endPoint;
+				currentSample ++, currentFilterElement++)
+			{
+				if (currentSample >= 0 && currentSample < numSamples)
+				{
+					sycl::vec<HilbertFirEnvelope::WorkType, H_VEC_SIZE> vec_sample(0.0);
+					#pragma unroll
+					for (int c = 0; c < H_VEC_SIZE; c++) {
+						vec_sample[c] = static_cast<HilbertFirEnvelope::WorkType>(signal[scanlineIdx + c +  currentSample * numScanlines]) 
+								* filter[currentFilterElement];
+					}
+					accumulator += vec_sample;
+				}
+			}
+			#pragma unroll
+			for (int c = 0; c < H_VEC_SIZE; c++) {
+				HilbertFirEnvelope::WorkType signalValue = static_cast<HilbertFirEnvelope::WorkType>(signal[scanlineIdx + c + sampleIdx*numScanlines]);
+				out[ scanlineIdx + c + sampleIdx * numScanlines ] = sycl::sqrt(squ(signalValue) + squ(accumulator[c]));
+			}
+		}
+
+	}
+
+
 	HilbertFirEnvelope::HilbertFirEnvelope(size_t filterLength)
 		: m_filterLength(filterLength)
 		, m_hilbertFilter(nullptr)
@@ -90,7 +139,7 @@ namespace supra
 		auto pEnv = make_shared<Container<OutputType> >(LocationGpu, inImageData->getStream(), numScanlines*numSamples);
 		sycl::range<3> blockSizeFilter(1, 8, 16);
 		sycl::range<3> gridSizeFilter(1, static_cast<unsigned int>((numSamples + blockSizeFilter[ 1 ] - 1) / blockSizeFilter[ 1 ]),
-									  static_cast<unsigned int>((numScanlines + blockSizeFilter[ 2 ] - 1) / blockSizeFilter[ 2 ]));
+									  static_cast<unsigned int>((numScanlines + blockSizeFilter[ 2 ] - 1) / blockSizeFilter[ 2 ] / H_VEC_SIZE));
 
 				static long hilbert_call_count = 0;
 
@@ -101,7 +150,7 @@ namespace supra
 						auto m_filterLength_ct5 = ( int )m_filterLength;
 
 						cgh.parallel_for(sycl::nd_range<3>(gridSizeFilter * blockSizeFilter, blockSizeFilter), [ = ](sycl::nd_item<3> item_ct1) {
-								kernelFilterDemodulation(inImageData_get_ct0, m_hilbertFilter_get_ct1, pEnv_get_ct2, numSamples, numScanlines, m_filterLength_ct5, item_ct1);
+								vec_kernelFilterDemodulation(inImageData_get_ct0, m_hilbertFilter_get_ct1, pEnv_get_ct2, numSamples, numScanlines, m_filterLength_ct5, item_ct1);
 						});
 				});
 
diff --git a/src/SupraLib/Beamformer/RxBeamformerCommon.h b/src/SupraLib/Beamformer/RxBeamformerCommon.h
index f86b5b9..35b05e0 100644
--- a/src/SupraLib/Beamformer/RxBeamformerCommon.h
+++ b/src/SupraLib/Beamformer/RxBeamformerCommon.h
@@ -34,6 +34,12 @@ namespace supra
 			(dirY*z)*(dirY*z)) + z;
 	}
 
+	const int Vec_SIZE = 8;
+	sycl::vec<float, Vec_SIZE> vec_computeDelayDTSPACE_D(float dirX, float dirY, float dirZ, sycl::vec<float, Vec_SIZE> x_element, float x, float z)
+	{
+		return sycl::sqrt(((x_element - (x + dirX * z)) * (x_element - (x + dirX * z)) + (dirY * z) * (dirY * z))) + z;
+	}
+
 	template <typename T>
 	inline T computeDelayDTSPACE3D_D(T dirX, T dirY, T dirZ, T x_element, T z_element, T x, T z, T d)
 	{
diff --git a/src/SupraLib/Beamformer/RxBeamformerCuda.dp.cpp b/src/SupraLib/Beamformer/RxBeamformerCuda.dp.cpp
index d97c382..7c7a4b7 100644
--- a/src/SupraLib/Beamformer/RxBeamformerCuda.dp.cpp
+++ b/src/SupraLib/Beamformer/RxBeamformerCuda.dp.cpp
@@ -305,6 +305,87 @@ namespace supra
 		}
 	}
 
+
+	const int ROW_SIZE = 2;
+	template <class SampleBeamformer, bool interpolateRFlines, bool interpolateBetweenTransmits, typename RFType, typename ResultType, typename LocationType>
+	void mixed_rxBeamformingDTSPACEKernel(size_t numTransducerElements, size_t numReceivedChannels, size_t numTimesteps, const RFType* __restrict__ RF, size_t numTxScanlines, size_t numRxScanlines,
+		const ScanlineRxParameters3D* __restrict__ scanlinesDT, size_t numDs, const LocationType* __restrict__ dsDT, const LocationType* __restrict__ x_elemsDT, LocationType speedOfSound,
+		LocationType dt, uint32_t additionalOffset, LocationType F, const WindowFunctionGpu* __restrict__ windowFunction, ResultType* __restrict__ s,   sycl::nd_item<3> item_ct1, float* __restrict__ mdataGpu)
+	{
+		int r = item_ct1.get_local_range().get(1) * item_ct1.get_group(1) + item_ct1.get_local_id(1);			 //@suppress("Symbol is not resolved")
+																												 //@suppress("Field cannot be resolved")
+		int scanlineIdx = item_ct1.get_local_range().get(2) * item_ct1.get_group(2) + item_ct1.get_local_id(2);	 //@suppress("Symbol is not resolved")
+		
+																												 //@suppress("Field cannot be resolved")
+		const int row_size = ROW_SIZE;
+		r = r * row_size;
+		if (r < numDs && scanlineIdx < numRxScanlines)
+		{
+			sycl::vec<LocationType, row_size> d(0.0);
+			#pragma unroll
+			for (int i = 0; i < row_size; i++) {
+				d[i] = dsDT[r + i];
+			}
+			// TODO should this also depend on the angle?
+			sycl::vec<float, row_size> aDT = d / (2*F);
+			ScanlineRxParameters3D scanline = scanlinesDT[ scanlineIdx ];
+			LocationType scanline_x = scanline.position.x;
+			LocationType  dirX = scanline.direction.x;
+			LocationType  dirY = scanline.direction.y;
+			LocationType  dirZ = scanline.direction.z;
+			LocationType  maxElementDistance = static_cast<LocationType>(scanline.maxElementDistance.x);
+
+			sycl::vec<float, row_size> sInterp(0.0);
+			
+			int highestWeightIndex;
+			if (!interpolateBetweenTransmits)
+			{
+				highestWeightIndex = 0;
+				float highestWeight = scanline.txWeights[ 0 ];
+				for (int k = 1; k < std::extent<decltype(scanline.txWeights)>::value; k++)
+				{
+					if (scanline.txWeights[ k ] > highestWeight)
+					{
+						highestWeight = scanline.txWeights[ k ];
+						highestWeightIndex = k;
+					}
+				}
+			}
+
+			// now iterate over all four txScanlines to interpolate beamformed scanlines from those transmits
+			for (int k = (interpolateBetweenTransmits ? 0 : highestWeightIndex);
+				 (interpolateBetweenTransmits && k < std::extent<decltype(scanline.txWeights)>::value) || (!interpolateBetweenTransmits && k == highestWeightIndex); k++)
+			{
+				if (scanline.txWeights[ k ] > 0.0)
+				{
+					ScanlineRxParameters3D::TransmitParameters txParams = scanline.txParameters[ k ];
+					uint32_t txScanlineIdx = txParams.txScanlineIdx;
+					if (txScanlineIdx >= numTxScanlines)
+					{
+						// ERROR!
+						return;
+					}
+
+					#pragma unroll
+					for (int i = 0; i < row_size; i++) {
+						LocationType invMaxElementDistance = 1.f / sycl::min(aDT[i], maxElementDistance);
+						sInterp[i] = SampleBeamformer::template vec_sampleBeamform2D<interpolateRFlines, RFType, float, LocationType>(txParams, RF, numTransducerElements, 
+						numReceivedChannels, numTimesteps, x_elemsDT, scanline_x, dirX, dirY, dirZ, aDT[i], d[i], invMaxElementDistance , speedOfSound, dt, additionalOffset, 
+						windowFunction, mdataGpu);
+					}
+					
+				}
+			}
+			
+			#pragma unroll
+			for (int i = 0; i < row_size; i++) {
+				s[scanlineIdx + (r + i) * numRxScanlines] = clampCast<ResultType>(sInterp[i]);
+			}
+		}
+	}
+
+
+
 	template <class SampleBeamformer, unsigned int maxWindowFunctionNumel, typename RFType, typename ResultType, typename LocationType>
 	void rxBeamformingDTspaceCuda3D(bool interpolateRFlines, bool interpolateBetweenTransmits, size_t numTransducerElements, vec2s elementLayout, size_t numReceivedChannels, size_t numTimesteps,
 									const RFType* RF, size_t numTxScanlines, size_t numRxScanlines, const ScanlineRxParameters3D* scanlines, size_t numZs, const LocationType* zs,
@@ -395,10 +476,11 @@ namespace supra
 	template <class SampleBeamformer, typename RFType, typename ResultType, typename LocationType>
 	void rxBeamformingDTspaceCuda(bool interpolateRFlines, bool interpolateBetweenTransmits, size_t numTransducerElements, size_t numReceivedChannels, size_t numTimesteps, const RFType* RF,
 								  size_t numTxScanlines, size_t numRxScanlines, const ScanlineRxParameters3D* scanlines, size_t numZs, const LocationType* zs, const LocationType* x_elems,
-								  LocationType speedOfSound, LocationType dt, uint32_t additionalOffset, LocationType F, const WindowFunctionGpu windowFunction, sycl::queue* stream, ResultType* s)
+								  LocationType speedOfSound, LocationType dt, uint32_t additionalOffset, LocationType F, const WindowFunctionGpu windowFunction, sycl::queue* stream, ResultType* s, LocationType *mdataGpu)
 	{
 		sycl::range<3> blockSize(1, 256, 1);
-		sycl::range<3> gridSize(1, static_cast<unsigned int>((numZs + blockSize[ 1 ] - 1) / blockSize[ 1 ]), static_cast<unsigned int>((numRxScanlines + blockSize[ 2 ] - 1) / blockSize[ 2 ]));
+		sycl::range<3> gridSize(1, static_cast<unsigned int>((numZs + blockSize[ 1 ] - 1) / blockSize[ 1 ] / ROW_SIZE), 
+								static_cast<unsigned int>((numRxScanlines + blockSize[ 2 ] - 1) / blockSize[ 2 ]));
 		if (interpolateRFlines)
 		{
 			if (interpolateBetweenTransmits)
@@ -417,8 +499,13 @@ namespace supra
 				
 						sycl::event beam_event = stream->submit([ & ](sycl::handler& cgh) {
 							  cgh.parallel_for(sycl::nd_range<3>(gridSize * blockSize, blockSize), [ = ](sycl::nd_item<3> item_ct1) {
-									rxBeamformingDTSPACEKernel<SampleBeamformer, true, false>(numTransducerElements, numReceivedChannels, numTimesteps, RF, numTxScanlines, numRxScanlines, scanlines,
-																							  numZs, zs, x_elems, speedOfSound, dt, additionalOffset, F, windowFunction, s, item_ct1);
+									/*rxBeamformingDTSPACEKernel<SampleBeamformer, true, false>(numTransducerElements, numReceivedChannels, numTimesteps, RF, numTxScanlines, numRxScanlines, scanlines,
+																							  numZs, zs, x_elems, speedOfSound, dt, additionalOffset, F, windowFunction, s, item_ct1);*/
+							  
+							  		mixed_rxBeamformingDTSPACEKernel<SampleBeamformer, true, false>(numTransducerElements, numReceivedChannels, numTimesteps, RF, numTxScanlines, numRxScanlines,
+																								   scanlines, numZs, zs, x_elems, speedOfSound, dt, additionalOffset, F, &windowFunction, s, item_ct1, mdataGpu);
+									 
+							  
 							  });
 						});
 
@@ -480,6 +567,12 @@ namespace supra
 			m_windowFunction = std::unique_ptr<WindowFunction>(new WindowFunction(windowType, windowParameter, m_windowFunctionNumEntries));
 		}
 
+		auto mdataGpu = (float*) sycl::malloc_device(m_windowFunctionNumEntries * sizeof(windowType), gRawData->getStream()->get_device(), gRawData->getStream()->get_context());
+		gRawData->getStream()->submit([&] (sycl::handler &h) {
+			h.memcpy(mdataGpu, m_windowFunction->m_data.data(), m_windowFunctionNumEntries * sizeof(WindowType));
+		});
+		gRawData->getStream()->wait();
+
 		auto beamformingFunction3D = &rxBeamformingDTspaceCuda3D<RxSampleBeamformerDelayAndSum, m_windowFunctionNumEntries, ChannelDataType, ImageDataType, LocationType>;
 		auto beamformingFunction2D = &rxBeamformingDTspaceCuda<RxSampleBeamformerDelayAndSum, ChannelDataType, ImageDataType, LocationType>;
 		
@@ -550,7 +643,8 @@ namespace supra
 				static_cast<LocationType>(fNumber),
 				*(m_windowFunction->getGpu()),
 				gRawData->getStream(),
-				pData->get()
+				pData->get(),
+				mdataGpu
 				);
 		}
 
diff --git a/src/SupraLib/Beamformer/RxSampleBeamformerDelayAndSum.h b/src/SupraLib/Beamformer/RxSampleBeamformerDelayAndSum.h
index f5330d5..b9c574d 100644
--- a/src/SupraLib/Beamformer/RxSampleBeamformerDelayAndSum.h
+++ b/src/SupraLib/Beamformer/RxSampleBeamformerDelayAndSum.h
@@ -17,6 +17,7 @@
 #include "USImageProperties.h"
 #include "WindowFunction.h"
 #include "RxBeamformerCommon.h"
+#include "helper.h"
 
 //TODO ALL ELEMENT/SCANLINE Y positons are actually Z! Change all variable names accordingly
 namespace supra
@@ -180,6 +181,110 @@ namespace supra
 				return 0;
 			}
 		}
+
+		template <bool interpolateRFlines, typename RFType, typename ResultType, typename LocationType>
+		static ResultType vec_sampleBeamform2D(
+			ScanlineRxParameters3D::TransmitParameters txParams, 
+			const RFType* RF, 
+			uint32_t numTransducerElements, 
+			uint32_t numReceivedChannels,
+			uint32_t numTimesteps, 
+			const LocationType* x_elemsDT, 
+			LocationType scanline_x, 
+			LocationType dirX, 
+			LocationType dirY, 
+			LocationType dirZ, 
+			LocationType aDT,
+			LocationType depth, 
+			LocationType invMaxElementDistance, 
+			LocationType speedOfSound, 
+			LocationType dt, 
+			int32_t additionalOffset,
+			const WindowFunctionGpu* __restrict__ windowFunction,
+			const float* mdataGpu)
+		{
+			const int VEC_SIZE = Vec_SIZE;
+			float sampleAcum = 0.0f;	
+			float weightAcum = 0.0f;
+			int numAdds = 0;
+			LocationType initialDelay = txParams.initialDelay;
+			uint32_t txScanlineIdx = txParams.txScanlineIdx;
+
+			for (int32_t elemIdxX = txParams.firstActiveElementIndex.x; elemIdxX < txParams.lastActiveElementIndex.x; elemIdxX += VEC_SIZE)
+			{
+				sycl::vec<int, VEC_SIZE> channelIdx;
+				sycl::vec<LocationType, VEC_SIZE> x_elem;
+
+				#pragma unroll
+				for (int i = 0; i < VEC_SIZE; i +=2) {
+					channelIdx[i] = (elemIdxX + i) % numReceivedChannels;
+					channelIdx[i+1] = (elemIdxX + i + 1) % numReceivedChannels;
+					x_elem[i] = x_elemsDT[elemIdxX + i];
+					x_elem[i + 1] = x_elemsDT[elemIdxX + i + 1];
+				}
+				sycl::vec<float, VEC_SIZE> sample;
+				sycl::vec<int, VEC_SIZE> mask = (sycl::fabs(x_elem - scanline_x) <= aDT);
+				/*sycl spec1.2.1 mentioned: true return  -1, false return 0*/
+				mask *= -1;
+				numAdds += utils<int, VEC_SIZE>::add_vec(mask);
+
+				sycl::vec<float, VEC_SIZE> relativeIndex = (x_elem - scanline_x) * invMaxElementDistance;
+				sycl::vec<float, VEC_SIZE> relativeIndexClamped = sycl::min(sycl::max(relativeIndex, -1.0f), 1.0f);	
+				sycl::vec<float, VEC_SIZE> absoluteIndex = windowFunction->m_scale * (relativeIndexClamped + 1.0f);	
+				sycl::vec<int, VEC_SIZE> absoluteIndex_int = absoluteIndex.convert<int, sycl::rounding_mode::automatic>();
+				sycl::vec<float, VEC_SIZE> weight;
+
+				#pragma unroll
+				for (int i = 0; i < VEC_SIZE; i += 2 ) {
+					weight[i] = mdataGpu[absoluteIndex_int[i]];
+					weight[i + 1] = mdataGpu[absoluteIndex_int[i + 1]];
+				}
+
+				weight *= mask.convert<float, sycl::rounding_mode::automatic>();
+				weightAcum += utils<float, VEC_SIZE>::add_vec(weight);
+
+				sycl::vec<LocationType, VEC_SIZE> delayf = initialDelay +
+					vec_computeDelayDTSPACE_D(dirX, dirY, dirZ, x_elem, scanline_x, depth) + additionalOffset;
+				sycl::vec<float, VEC_SIZE> delay = sycl::floor(delayf);
+				sycl::vec<int, VEC_SIZE> delay_index = delay.convert<int, sycl::rounding_mode::automatic>() + channelIdx*numTimesteps + 
+					txScanlineIdx*numReceivedChannels*numTimesteps;
+				
+				delayf -= delay;
+
+				sycl::vec<float, VEC_SIZE> RF_data;
+				sycl::vec<float, VEC_SIZE> RF_data_one;
+				
+				#pragma unroll
+				for (int j = 0; j < VEC_SIZE; j += 2) {
+					RF_data[j] = (float)RF[delay_index[j]];
+					RF_data[j + 1] = (float)RF[delay_index[j + 1]];
+					RF_data_one[j] = (float)RF[delay_index[j]+1];
+					RF_data_one[j + 1] = (float)RF[delay_index[j + 1]+1];
+				}
+
+				sycl::vec<int, VEC_SIZE> mask1 = (delay < (numTimesteps - 1));
+				mask1 *= -1;
+				sample = weight * ((1.0f - delayf) * RF_data + 
+					delayf  * RF_data_one) * mask1.convert<float, sycl::rounding_mode::automatic>();
+				sampleAcum += utils<float, VEC_SIZE>::add_vec(sample);
+				
+				sycl::vec<int, VEC_SIZE> mask2 = (delay < numTimesteps && delayf == 0.0);
+				mask2 *= -1;
+				sample = weight * RF_data * mask2.convert<float, sycl::rounding_mode::automatic>();
+				sampleAcum += utils<float, VEC_SIZE>::add_vec(sample);
+				
+			}
+
+			if (numAdds > 0)
+			{
+				return sampleAcum / weightAcum * numAdds;
+			}
+			else
+			{
+				return 0;
+			}
+		}
+
 	};
 }
 
diff --git a/src/SupraLib/Beamformer/WindowFunction.h b/src/SupraLib/Beamformer/WindowFunction.h
index 4f69c89..71c2267 100644
--- a/src/SupraLib/Beamformer/WindowFunction.h
+++ b/src/SupraLib/Beamformer/WindowFunction.h
@@ -88,7 +88,7 @@ namespace supra
 		{
 			return m_numEntriesPerFunction;
 		}
-
+		float m_scale;
 	private:
 		friend WindowFunction;
 		WindowFunctionGpu(size_t numEntriesPerFunction, const ElementType* data)
@@ -96,7 +96,7 @@ namespace supra
 			, m_data(data)
 			, m_scale(static_cast<float>(numEntriesPerFunction - 1)*0.5f) {};
 
-		float m_scale;
+		
 		uint32_t m_numEntriesPerFunction;
 		const ElementType* m_data;
 	};
@@ -136,9 +136,11 @@ namespace supra
 				return 0;
 			}
 		}
+
+		std::vector<ElementType> m_data;
 	private:
 		size_t m_numEntriesPerFunction;
-		std::vector<ElementType> m_data;
+		
 		std::unique_ptr<Container<ElementType> > m_dataGpu;
 		ElementType m_scale;
 		WindowType m_type;
diff --git a/src/SupraLib/Beamformer/helper.h b/src/SupraLib/Beamformer/helper.h
new file mode 100644
index 0000000..15cfe21
--- /dev/null
+++ b/src/SupraLib/Beamformer/helper.h
@@ -0,0 +1,159 @@
+//==---------- helper.hpp - SYCL sub_group helper functions ----------------==//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+#include <CL/sycl.hpp>
+#include <cmath>
+#include <iostream>
+
+using namespace cl::sycl;
+
+template <typename T1, int N> struct utils {
+  static T1 add_vec(const vec<T1, N> &v);
+  static bool cmp_vec(const vec<T1, N> &v, const vec<T1, N> &r);
+  static std::string stringify_vec(const vec<T1, N> &v);
+};
+template <typename T2> struct utils<T2, 1> {
+  static T2 add_vec(const vec<T2, 1> &v) { return v.s0(); }
+  static bool cmp_vec(const vec<T2, 1> &v, const vec<T2, 1> &r) {
+    return v.s0() == r.s0();
+  }
+  static std::string stringify_vec(const vec<T2, 1> &v) {
+    return std::to_string((T2)v.s0());
+  }
+};
+template <typename T2> struct utils<T2, 2> {
+  static T2 add_vec(const vec<T2, 2> &v) { return v.s0() + v.s1(); }
+  static bool cmp_vec(const vec<T2, 2> &v, const vec<T2, 2> &r) {
+    return v.s0() == r.s0() && v.s1() == r.s1();
+  }
+  static std::string stringify_vec(const vec<T2, 2> &v) {
+    return std::string("(") + std::to_string((T2)v.s0()) + ", " +
+           std::to_string((T2)v.s1()) + " )";
+  }
+};
+template <typename T2> struct utils<T2, 4> {
+  static T2 add_vec(const vec<T2, 4> &v) {
+    return v.s0() + v.s1() + v.s2() + v.s3();
+  }
+  static bool cmp_vec(const vec<T2, 4> &v, const vec<T2, 4> &r) {
+    return v.s0() == r.s0() && v.s1() == r.s1() && v.s2() == r.s2() &&
+           v.s3() == r.s3();
+  }
+  static std::string stringify_vec(const vec<T2, 4> &v) {
+    return std::string("(") + std::to_string((T2)v.s0()) + ", " +
+           std::to_string((T2)v.s1()) + std::to_string((T2)v.s2()) + ", " +
+           std::to_string((T2)v.s3()) + " )";
+  }
+};
+template <typename T2> struct utils<T2, 8> {
+  static T2 add_vec(const vec<T2, 8> &v) {
+    return v.s0() + v.s1() + v.s2() + v.s3() + v.s4() + v.s5() + v.s6() +
+           v.s7();
+  }
+  static bool cmp_vec(const vec<T2, 8> &v, const vec<T2, 8> &r) {
+    return v.s0() == r.s0() && v.s1() == r.s1() && v.s2() == r.s2() &&
+           v.s3() == r.s3() && v.s4() == r.s4() && v.s5() == r.s5() &&
+           v.s6() == r.s6() && v.s7() == r.s7();
+  }
+  static std::string stringify_vec(const vec<T2, 8> &v) {
+    return std::string("(") + std::to_string((T2)v.s0()) + ", " +
+           std::to_string((T2)v.s1()) + std::to_string((T2)v.s2()) + ", " +
+           std::to_string((T2)v.s3()) + std::to_string((T2)v.s4()) + ", " +
+           std::to_string((T2)v.s5()) + std::to_string((T2)v.s6()) + ", " +
+           std::to_string((T2)v.s7()) + " )";
+  }
+};
+
+template <typename T2> struct utils<T2, 16> {
+  static T2 add_vec(const vec<T2, 16> &v) {
+    return v.s0() + v.s1() + v.s2() + v.s3() + v.s4() + v.s5() + v.s6() +
+           v.s7() + v.s8() + v.s9() + v.sA() + v.sB() + v.sC() + v.sD() +
+           v.sE() + v.sF();
+  }
+  static bool cmp_vec(const vec<T2, 16> &v, const vec<T2, 16> &r) {
+    return v.s0() == r.s0() && v.s1() == r.s1() && v.s2() == r.s2() &&
+           v.s3() == r.s3() && v.s4() == r.s4() && v.s5() == r.s5() &&
+           v.s6() == r.s6() && v.s7() == r.s7() && v.s8() == r.s8() &&
+           v.s9() == r.s9() && v.sA() == r.sA() && v.sB() == r.sB() &&
+           v.sC() == r.sC() && v.sD() == r.sD() && v.sE() == r.sE() &&
+           v.sF() == r.sF();
+  }
+  static std::string stringify_vec(const vec<T2, 16> &v) {
+    return std::string("(") + std::to_string((T2)v.s0()) + ", " +
+           std::to_string((T2)v.s1()) + std::to_string((T2)v.s2()) + ", " +
+           std::to_string((T2)v.s3()) + std::to_string((T2)v.s4()) + ", " +
+           std::to_string((T2)v.s5()) + std::to_string((T2)v.s6()) + ", " +
+           std::to_string((T2)v.s7()) + std::to_string((T2)v.s8()) + ", " +
+           std::to_string((T2)v.s9()) + std::to_string((T2)v.sA()) + ", " +
+           std::to_string((T2)v.sB()) + std::to_string((T2)v.sC()) + ", " +
+           std::to_string((T2)v.sE()) + std::to_string((T2)v.sD()) + ", " +
+           std::to_string((T2)v.sF()) + " )";
+  }
+};
+
+template <typename T> 
+inline void exit_if_not_equal(T val, T ref, const char *name) {
+  if (std::is_floating_point<T>::value) {
+    if (std::fabs(val - ref) > 0.01) {
+      std::cout << "Unexpected result for " << name << ": " << (double)val
+                << " expected value: " << (double)ref << std::endl;
+      exit(1);
+    }
+  } else {
+    if ((val - ref) != 0) {
+      std::cout << "Unexpected result for " << name << ": " << (long)val
+                << " expected value: " << (long)ref << std::endl;
+      exit(1);
+    }
+  }
+}
+
+template <> 
+inline void exit_if_not_equal(half val, half ref, const char *name) {
+  int16_t cmp_val = reinterpret_cast<int16_t&>(val);
+  int16_t cmp_ref = reinterpret_cast<int16_t&>(ref);
+  if (std::abs(cmp_val - cmp_ref) > 1) {
+    std::cout << "Unexpected result for " << name << ": " << (float)val
+              << " expected value: " << (float)ref << std::endl;
+    exit(1);
+  }
+}
+
+template <typename T, int N>
+inline void exit_if_not_equal_vec(vec<T, N> val, vec<T, N> ref, const char *name) {
+  if (!utils<T, N>::cmp_vec(ref, val)) {
+    std::cout << "Unexpected result for " << name << ": "
+              << utils<T, N>::stringify_vec(val)
+              << " expected value: " << utils<T, N>::stringify_vec(ref)
+              << std::endl;
+
+    exit(1);
+  }
+}
+
+/* CPU returns max number of SG, GPU returns max SG size for
+ * CL_DEVICE_MAX_NUM_SUB_GROUPS device parameter. This function aligns the
+ * value.
+ * */
+inline size_t get_sg_size(const device &Device) {
+  size_t max_num_sg = Device.get_info<info::device::max_num_sub_groups>();
+  if (Device.get_info<info::device::device_type>() == info::device_type::cpu) {
+    size_t max_wg_size = Device.get_info<info::device::max_work_group_size>();
+    return max_wg_size / max_num_sg;
+  }
+  if (Device.get_info<info::device::device_type>() == info::device_type::gpu) {
+    return max_num_sg;
+  }
+  std::cout << "Unexpected deive type" << std::endl;
+  exit(1);
+}
+
+inline bool core_sg_supported(const device &Device) {
+  return (Device.has_extension("cl_khr_subgroups") ||
+          Device.get_info<info::device::version>().find(" 2.1") !=
+              string_class::npos);
+}
-- 
2.17.1