On the surface, calculating a dot product is almost embarassingly parallel. But there are nuances to consider. Because the operation to be parallelized is only a simple add, the cost of transfering the data from CPU host to GPU device can quickly become prohibitive.

Blocked Multistream Version

On this system the host-device transfers go through a PCI-Express 3.0 x16 interface, which has a theoretical bandwidth of 15.5 GB/s. The maximum that has been achieved on this system in NVIDIA's publicly available bandwidthtest.cu test is 12.6 GB/s for host-to-device and 11.8 GB/s for device-to-host transfers.

Optimizing for memory transfers is therefore paramount for the GPU dot product to even have a chance of running competitively against the CPU dot products. As the AVX (multithreaded) (double) throughput is around 360 GB/s even at array sizes of 268M elements there is no hope for the GPU version compete. The single threaded AVX bandwidth is much lower however, already being at only 28 GB/s at 2M elements.

The approach will be the roughly following:

1. Split the input n into n=n_cpu+n_gpu. n_gpu is the largest power of two inside n.

2. Sum-reduce input arrays x,y of respective lengths n_gpu in parallel on the GPU.

3. Sum input arrays x,y of respective lengths n_cpu on the CPU via dotproduct_avx.

4. Add results from 2 and 3 together.

To parallelize memory transfers and computation as much as possible, async memory transfers and kernel launches via streams are helpful facilities. The computation via dotproduct_avx on the host therefore happens in parallel to the kernel computations on the GPU. This again introduces a lot of parameters whose optimal value depend on the computer:

1. MAX_STREAMS which determines the maximum amount of streams
2. MIN_N_PER_STREAM_SIZE which determines the minimum amount of elements needed in order to create a new streams
3. MAX_BLOCKS_PER_STREAM the maximum blocks a kernel within a stream is launched with
4. THREADS_PER_BLOCK the maximum threads a kernel within a stream is launched with

These values are set to the following:

1. MAX_STREAMS = 64
2. MIN_N_PER_STREAM_SIZE = 1<<12
3. MAX_BLOCKS_PER_STREAM = 32
4. THREADS_PER_BLOCK = 1024

Strictly speaking, the astute reader will immediately notice this probably won't help much on this computer but we will throw nevertheless the kitchen sink at this problem.

This leads to the following attempt:

    template<class T>
    struct VECTORIZED_CVT;

    template<>
    struct VECTORIZED_CVT<double>{
        using type=double4;
    };

    template<>
    struct VECTORIZED_CVT<float>{
        using type=float4;
    };

    namespace helper{



        template<class T>
        __device__ T get_zero_vector();

        template<>
        __device__ double4 get_zero_vector() { 
            return make_double4(0.0, 0.0, 0.0, 0.0); 
        }

        template<>
        __device__ float4 get_zero_vector() { 
            return make_float4(0.0, 0.0, 0.0, 0.0); 
        }

        __device__
        __forceinline__ void fma(double4* x, double4* y, double4& value){
            value.x=::fma((*x).x,(*y).x,value.x); //vectorized fma would be better
            value.y=::fma((*x).y,(*y).y,value.y); //vectorized fma would be better
            value.z=::fma((*x).z,(*y).z,value.z); //vectorized fma would be better
            value.w=::fma((*x).w,(*y).w,value.w); //vectorized fma would be better
        }

        __device__
        __forceinline__ void fma(float4* x, float4* y, float4& value){
            value.x=::fma((*x).x,(*y).x,value.x); //vectorized fma would be better
            value.y=::fma((*x).y,(*y).y,value.y); //vectorized fma would be better
            value.z=::fma((*x).z,(*y).z,value.z); //vectorized fma would be better
            value.w=::fma((*x).w,(*y).w,value.w); //vectorized fma would be better
        }

        template<int OFFSET>
        __device__
        __forceinline__ void shared_memory_add_assignment_of_vector(double4* sdata){
                sdata[threadIdx.x].x+=sdata[threadIdx.x+OFFSET].x;
                sdata[threadIdx.x].y+=sdata[threadIdx.x+OFFSET].y;
                sdata[threadIdx.x].z+=sdata[threadIdx.x+OFFSET].z;
                sdata[threadIdx.x].w+=sdata[threadIdx.x+OFFSET].w;
        }

        template<int OFFSET>
        __device__
        __forceinline__ void shared_memory_add_assignment_of_vector(float4* sdata){
                sdata[threadIdx.x].x+=sdata[threadIdx.x+OFFSET].x;
                sdata[threadIdx.x].y+=sdata[threadIdx.x+OFFSET].y;
                sdata[threadIdx.x].z+=sdata[threadIdx.x+OFFSET].z;
                sdata[threadIdx.x].w+=sdata[threadIdx.x+OFFSET].w;
        }


        //no syncthreads needed because I am in a warp
        template<size_t THREADS_PER_BLOCK, typename VECTORIZED_T>
        __device__
        __forceinline__ void reduce_warp(VECTORIZED_T* sdata){
            if(THREADS_PER_BLOCK>=64){      
                    shared_memory_add_assignment_of_vector<32>(sdata);
                __syncthreads();
            }
            if(THREADS_PER_BLOCK>=32){
                    shared_memory_add_assignment_of_vector<16>(sdata);
                __syncthreads();
            }
            if(THREADS_PER_BLOCK>=16){
                    shared_memory_add_assignment_of_vector<8>(sdata);
                __syncthreads();
            }
            if(THREADS_PER_BLOCK>=8){
                shared_memory_add_assignment_of_vector<4>(sdata);
                __syncthreads();
            }
            if(THREADS_PER_BLOCK>=4){
                    shared_memory_add_assignment_of_vector<2>(sdata);
                __syncthreads();
            }
            if(THREADS_PER_BLOCK>=2){
                    shared_memory_add_assignment_of_vector<1>(sdata);
                __syncthreads();
            }           
        }

        double vectorized_type_to_scalar_sum(double4& x){
            return x.x+x.y+x.z+x.w;
        }

        float vectorized_type_to_scalar_sum(float4& x){
            return x.x+x.y+x.z+x.w;
        }
    }

    template<size_t THREADS_PER_BLOCK, typename VECTORIZED_T>
    __global__
    void reduce_k(size_t n, VECTORIZED_T* x, VECTORIZED_T* y, VECTORIZED_T* r){
        constexpr size_t memsize=(THREADS_PER_BLOCK<=64)?64:THREADS_PER_BLOCK;
        using T=decltype(x[0].x);
        static __shared__ VECTORIZED_T sdata[memsize];
        VECTORIZED_T value = helper::get_zero_vector<VECTORIZED_T>();
        sdata[threadIdx.x] = helper::get_zero_vector<VECTORIZED_T>();
        const int stride = gridDim.x*blockDim.x;

        for (int i=blockDim.x*blockIdx.x+threadIdx.x; i<n; i+=stride){
            helper::fma(x+i,y+i,value); 
        }
        sdata[threadIdx.x]=value;
        __syncthreads();

        if constexpr (THREADS_PER_BLOCK>=1024){
            if (threadIdx.x < 512) { 
                helper::shared_memory_add_assignment_of_vector<512>(sdata);
            }
            __syncthreads(); 
        }
        if constexpr (THREADS_PER_BLOCK>=512){
            if (threadIdx.x < 256) { 
                helper::shared_memory_add_assignment_of_vector<256>(sdata);
            }
            __syncthreads(); 
        }

        if constexpr (THREADS_PER_BLOCK>=256){
            if (threadIdx.x < 128) { 
                helper::shared_memory_add_assignment_of_vector<128>(sdata);
            }
            __syncthreads(); 
        }

        if constexpr (THREADS_PER_BLOCK>=128){
            if (threadIdx.x < 64) { 
                helper::shared_memory_add_assignment_of_vector<64>(sdata);
            }
            __syncthreads(); 
        }       
        if (threadIdx.x<32){
            helper::reduce_warp<THREADS_PER_BLOCK,VECTORIZED_T>(sdata);
            if (threadIdx.x == 0){ //If statement not needed because all threads because all threads write the same value
                r[blockIdx.x]=sdata[0];
            }
        }
    }

    template<size_t MAX_STREAMS, size_t MIN_N_PER_STREAM_SIZE, size_t MAX_BLOCKS_PER_STREAM, size_t THREADS_PER_BLOCK, typename T> requires dotprod::IsPositivePowerOfTwo<MAX_BLOCKS_PER_STREAM> && dotprod::IsPositivePowerOfTwo<THREADS_PER_BLOCK>
    __host__
    T reduce(size_t n, T* x_pinned, T* y_pinned){
        if (n == 0){
            return T(0);
        }
        using VECTORIZED_T=typename VECTORIZED_CVT<T>::type;
        int VECTORIZED_ELEMENTS=sizeof(VECTORIZED_T)/sizeof(T);

        T res=T(0);
        int n_power_of_two=1u<<dotprod::log2_int(n); 
        uint64_t streams_needed_no_clamp=(n_power_of_two+static_cast<long long>(MIN_N_PER_STREAM_SIZE-1LL))/MIN_N_PER_STREAM_SIZE;
        size_t streams_needed=(streams_needed_no_clamp<MAX_STREAMS)?streams_needed_no_clamp:MAX_STREAMS; //AND only works if MAX_STREAMS is power of two
        int n_per_stream=1u<<dotprod::log2_int(n_power_of_two/streams_needed);
        int n_vectorized_per_stream=n_per_stream/VECTORIZED_ELEMENTS;
        int n_gpu=n_per_stream*streams_needed;
        size_t n_cpu=n-n_gpu;
        int blocks_used=(((n_vectorized_per_stream+THREADS_PER_BLOCK-1)/THREADS_PER_BLOCK)<MAX_BLOCKS_PER_STREAM)?(n_vectorized_per_stream+THREADS_PER_BLOCK-1)/THREADS_PER_BLOCK:MAX_BLOCKS_PER_STREAM;        
        VECTORIZED_T* r_d;

        VECTORIZED_T* r_h = new VECTORIZED_T[MAX_STREAMS*MAX_BLOCKS_PER_STREAM]; 
        CHECK_CUDA_ERROR(cudaMalloc((void**)&r_d, sizeof(VECTORIZED_T)*streams_needed*blocks_used));

        cudaStream_t streams[MAX_STREAMS];
        for (size_t i=0;i<streams_needed;i++){
            size_t bytes_to_copy=sizeof(T)*n_per_stream;
            VECTORIZED_T* _x_d;
            VECTORIZED_T* _y_d;
            VECTORIZED_T* _r_d=r_d+i*blocks_used;
            cudaStreamCreate (&(streams[i]));
            CHECK_CUDA_ERROR(cudaMallocAsync((void**)&_x_d, bytes_to_copy, streams[i]));
            CHECK_CUDA_ERROR(cudaMallocAsync((void**)&_y_d, bytes_to_copy, streams[i]));
            CHECK_CUDA_ERROR(cudaMemcpyAsync (_x_d, reinterpret_cast<VECTORIZED_T*>(&(x_pinned[i*n_per_stream])), bytes_to_copy, cudaMemcpyHostToDevice, streams[i] ));
            CHECK_CUDA_ERROR(cudaMemcpyAsync (_y_d, reinterpret_cast<VECTORIZED_T*>(&(y_pinned[i*n_per_stream])), bytes_to_copy, cudaMemcpyHostToDevice, streams[i] ));
            reduce_k<THREADS_PER_BLOCK,VECTORIZED_T><<<dim3(blocks_used,1,1),dim3(THREADS_PER_BLOCK,1,1),0,streams[i]>>>(n_vectorized_per_stream,_x_d,_y_d,_r_d);
            cudaFreeAsync(_x_d,streams[i]);
            cudaFreeAsync(_y_d,streams[i]);

        }
        res+=avx::dp<0,0>(n_cpu,x_pinned+n_gpu,y_pinned+n_gpu);
        for (size_t i=0;i<streams_needed;i++){
            CHECK_CUDA_ERROR(cudaStreamSynchronize(streams[i]));
            CHECK_CUDA_ERROR(cudaStreamDestroy(streams[i]));
        }

        CHECK_CUDA_ERROR(cudaMemcpy(r_h, r_d, sizeof(VECTORIZED_T)*streams_needed*blocks_used, cudaMemcpyDeviceToHost));            
        for (int i=0;i<streams_needed*blocks_used;i++){
            res+=helper::vectorized_type_to_scalar_sum(r_h[i]);
        }
        delete[] r_h;
        cudaFree(r_d);
        return res;
    }

While there is surely room for improvement, since memory transfers are the primary bottleneck hyper optimization of computational operations is not as urgent. In order to ensure certain parameters are a power of two, the IsPositivePowerOfTwo template constraint described in the appendix comes in handy.

The CUDA version is even slower than the naive loop. Through clever selection of the block and stream parameters the CUDA version might come closer to parity but due to the aforementioned bandwidth bottleneck there is little hope to become significantly faster. The GPU bandwidth is low and clearly limited by the PCI-Express 3.0 x16 bandwidth of 15.5 GB/s. The NVIDIA bandwidthtest.cu serves as a good lower limit because this is the minimum time needed just to copy the data to device.

Note: Since this is a CPU-GPU hybrid algorithm, for some array sizes the throughput may be slightly greater especially if n_cpu ≈ n_gpu.

As the double and float versions both max out at around the same throughput, the float version can on work a lot more elements in a given time.

An interesting side effect is that the chosen reduction approach has different numerical error bounds than the sequential additions in the single threaded AVX dot product before.

Special Case Optimization: x = y

As for the AVX version it is possible to optimize for sum of squares too.

template<size_t MAX_STREAMS, size_t MIN_N_PER_STREAM_SIZE, size_t MAX_BLOCKS_PER_STREAM, size_t THREADS_PER_BLOCK, typename T> requires dotprod::IsPositivePowerOfTwo<MAX_BLOCKS_PER_STREAM> && dotprod::IsPositivePowerOfTwo<THREADS_PER_BLOCK>
__host__
T reduce(size_t n, T* x_pinned){
    if (n == 0){
        return T(0);
    }
    using VECTORIZED_T=typename VECTORIZED_CVT<T>::type;
    int VECTORIZED_ELEMENTS=sizeof(VECTORIZED_T)/sizeof(T);

    T res=T(0);
    int n_power_of_two=1u<<dotprod::log2_int(n); //todo change
    uint64_t streams_needed_no_clamp=(n_power_of_two+static_cast<long long>(MIN_N_PER_STREAM_SIZE-1LL))/MIN_N_PER_STREAM_SIZE;
    size_t streams_needed=(streams_needed_no_clamp<MAX_STREAMS)?streams_needed_no_clamp:MAX_STREAMS; //AND only works if MAX_STREAMS is power of two
    int n_per_stream=1u<<dotprod::log2_int(n_power_of_two/streams_needed);
    int n_vectorized_per_stream=n_per_stream/VECTORIZED_ELEMENTS;
    int n_gpu=n_per_stream*streams_needed;
    size_t n_cpu=n-n_gpu;
    int blocks_used=(((n_vectorized_per_stream+THREADS_PER_BLOCK-1)/THREADS_PER_BLOCK)<MAX_BLOCKS_PER_STREAM)?(n_vectorized_per_stream+THREADS_PER_BLOCK-1)/THREADS_PER_BLOCK:MAX_BLOCKS_PER_STREAM;        
    VECTORIZED_T* r_d;

    VECTORIZED_T* r_h = new VECTORIZED_T[MAX_STREAMS*MAX_BLOCKS_PER_STREAM]; 
    CHECK_CUDA_ERROR(cudaMalloc((void**)&r_d, sizeof(VECTORIZED_T)*streams_needed*blocks_used));

    cudaStream_t streams[MAX_STREAMS];
    for (size_t i=0;i<streams_needed;i++){
        size_t bytes_to_copy=sizeof(T)*n_per_stream;
        VECTORIZED_T* _x_d;
        VECTORIZED_T* _r_d=r_d+i*blocks_used;
        cudaStreamCreate (&(streams[i]));
        CHECK_CUDA_ERROR(cudaMallocAsync((void**)&_x_d, bytes_to_copy, streams[i]));
        CHECK_CUDA_ERROR(cudaMemcpyAsync (_x_d, reinterpret_cast<VECTORIZED_T*>(&(x_pinned[i*n_per_stream])), bytes_to_copy, cudaMemcpyHostToDevice, streams[i] ));
        reduce_k<THREADS_PER_BLOCK,VECTORIZED_T><<<dim3(blocks_used,1,1),dim3(THREADS_PER_BLOCK,1,1),0,streams[i]>>>(n_vectorized_per_stream,_x_d,_x_d,_r_d);
        cudaFreeAsync(_x_d,streams[i]);

    }
    res+=avx::dp_sq<0,0>(n_cpu,x_pinned+n_gpu);
    for (size_t i=0;i<streams_needed;i++){
        CHECK_CUDA_ERROR(cudaStreamSynchronize(streams[i]));
        CHECK_CUDA_ERROR(cudaStreamDestroy(streams[i]));
    }

    CHECK_CUDA_ERROR(cudaMemcpy(r_h, r_d, sizeof(VECTORIZED_T)*streams_needed*blocks_used, cudaMemcpyDeviceToHost));            
    for (int i=0;i<streams_needed*blocks_used;i++){
        res+=helper::vectorized_type_to_scalar_sum(r_h[i]);
    }
    delete[] r_h;
    cudaFree(r_d);
    return res;
}

Again, the practical PCI-Express 3.0 x16 bandwidth cannot be exceeded for n_cpu ≪ n_gpu.