<commit_before>/*
 * %%%%%%%%%%%%%%%%%%%%%%%Original development%%%%%%%%%%%%%%%%%%%%%%%%%
 *  Copyright (c) 2016 Inria and University Pierre and <NAME>
 *  All rights reserved.
 * %%%%%%%%%%%%%%%%%%%%%%%Modifications and further additions%%%%%%%%%%
 *  <NAME>, 2017, within FELTOR and EXBLAS licenses
 */
/**
 *  @file exdot_cuda.cuh
 *  @brief CUDA version of exdot
 *
 *  @authors
 *    Developers : \n
 *        <NAME>  -- <EMAIL> \n
 *        <NAME> -- <EMAIL> \n
 *        <NAME> -- <EMAIL>
 */
#pragma once
#include "thrust/device_vector.h"
#include "accumulate.cuh"

namespace exblas{
///@cond
namespace gpu{



////////////////////////////////////////////////////////////////////////////////
// Auxiliary functions
////////////////////////////////////////////////////////////////////////////////
__device__
double TwoProductFMA(double a, double b, double *d) {
    double p = a * b;
    *d = fma(a, b, -p);
    return p;
}

__device__
double KnuthTwoSum(double a, double b, double *s) {
    double r = a + b;
    double z = r - a;
    *s = (a - (r - z)) + (b - z);
    return r;
}


template<uint NBFPE, uint WARP_COUNT>
__global__ void ExDOT(
    int64_t *d_PartialSuperaccs,
    const double *d_a,
    const double *d_b,
    const uint NbElements
) {
    __shared__ int64_t l_sa[WARP_COUNT * BIN_COUNT]; //shared variables live for a thread block (39 rows, 16 columns!)
    int64_t *l_workingBase = l_sa + (threadIdx.x & (WARP_COUNT - 1)); //the bitwise & with 15 is a modulo operation: threadIdx.x % 16
    //Initialize superaccs
    for (uint i = 0; i < BIN_COUNT; i++)
        l_workingBase[i * WARP_COUNT] = 0;
    __syncthreads();

    //Read data from global memory and scatter it to sub-superaccs
    double a[NBFPE] = {0.0};
    for(uint pos = blockIdx.x*blockDim.x+threadIdx.x; pos < NbElements; pos += gridDim.x*blockDim.x) {
        double r = 0.0;
        double x = TwoProductFMA(d_a[pos], d_b[pos], &r);
        //double x = d_a[pos]*d_b[pos];//ATTENTION: if we write it like this, cpu compiler might generate an fma from this while nvcc does not...

        #pragma unroll
        for(uint i = 0; i != NBFPE; ++i) {
            double s;
            a[i] = KnuthTwoSum(a[i], x, &s);
            x = s;
        }
        if (x != 0.0) {
            Accumulate(l_workingBase, x, WARP_COUNT);
            // Flush FPEs to superaccs
            #pragma unroll
            for(uint i = 0; i != NBFPE; ++i) {
                Accumulate(l_workingBase, a[i], WARP_COUNT);
                a[i] = 0.0;
            }
        }

        if (r != 0.0) {//add the rest r in the same manner
            #pragma unroll
            for(uint i = 0; i != NBFPE; ++i) {
                double s;
                a[i] = KnuthTwoSum(a[i], r, &s);
                r = s;
            }
            if (r != 0.0) {
                Accumulate(l_workingBase, r, WARP_COUNT);
                // Flush FPEs to superaccs
                #pragma unroll
                for(uint i = 0; i != NBFPE; ++i) {
                    Accumulate(l_workingBase, a[i], WARP_COUNT);
                    a[i] = 0.0;
                }
            }
        }
    }
    //Flush FPEs to superaccs
    #pragma unroll
    for(uint i = 0; i != NBFPE; ++i)
        Accumulate(l_workingBase, a[i], WARP_COUNT);
    __syncthreads();

    //Merge sub-superaccs into work-group partial-accumulator ( ATTENTION: PartialSuperacc is transposed!)
    uint pos = threadIdx.x;
    if(pos < WARP_COUNT) {
        int imin = IMIN, imax = IMAX;
        Normalize( l_workingBase, imin, imax, WARP_COUNT);
    }
    __syncthreads();

    if (pos < BIN_COUNT) {
        int64_t sum = 0;

        for(uint i = 0; i < WARP_COUNT; i++)
            sum += l_sa[pos * WARP_COUNT + i];

        d_PartialSuperaccs[blockIdx.x * BIN_COUNT + pos] = sum;
    }

    __syncthreads();
    if (pos == 0) {
        int imin = IMIN, imax = IMAX;
        Normalize(&d_PartialSuperaccs[blockIdx.x * BIN_COUNT], imin, imax);
    }
}

template<uint NBFPE, uint WARP_COUNT>
__global__ void ExDOT(
    int64_t *d_PartialSuperaccs,
    const double *d_a,
    const double *d_b,
    const double *d_c,
    const uint NbElements
) {
    __shared__ int64_t l_sa[WARP_COUNT * BIN_COUNT]; //shared variables live for a thread block (39 rows, 16 columns!)
    int64_t *l_workingBase = l_sa + (threadIdx.x & (WARP_COUNT - 1)); //the bitwise & with 15 is a modulo operation: threadIdx.x % 16
    //Initialize superaccs
    for (uint i = 0; i < BIN_COUNT; i++)
        l_workingBase[i * WARP_COUNT] = 0;
    __syncthreads();

    //Read data from global memory and scatter it to sub-superaccs
    double a[NBFPE] = {0.0};
    for(uint pos = blockIdx.x*blockDim.x+threadIdx.x; pos < NbElements; pos += gridDim.x*blockDim.x) {
        //double x2 = d_a[pos]*d_c[pos]*d_b[pos];
        //double r  = 0.0, r2 = 0.0;
        //double x  = TwoProductFMA(d_a[pos], d_b[pos], &r);
        //double x2 = TwoProductFMA(x , d_c[pos], &r2);
        double x1 = fma( d_a[pos], d_b[pos], 0);
        double x2 = fma( x1      , d_c[pos], 0);

        if( x2 != 0.0) {//accumulate x2
            #pragma unroll
            for(uint i = 0; i != NBFPE; ++i) {
                double s;
                a[i] = KnuthTwoSum(a[i], x2, &s);
                x2 = s;
            }
            if (x2 != 0.0) {
                Accumulate(l_workingBase, x2, WARP_COUNT);
                // Flush FPEs to superaccs
                #pragma unroll
                for(uint i = 0; i != NBFPE; ++i) {
                    Accumulate(l_workingBase, a[i], WARP_COUNT);
                    a[i] = 0.0;
                }
            }
        }
        //if (r2 != 0.0) {//add the rest r2
        //    #pragma unroll
        //    for(uint i = 0; i != NBFPE; ++i) {
        //        double s;
        //        a[i] = KnuthTwoSum(a[i], r2, &s);
        //        r2 = s; //error was here r = s
        //    }
        //    if (r2 != 0.0) { //error was here r != 0.0
        //        Accumulate(l_workingBase, r2, WARP_COUNT);
        //        // Flush FPEs to superaccs<commit_msg>HOTFIX: changed fma to intrinsic __fma_rn in exdot_cuda.cuh<commit_after>/*
 * %%%%%%%%%%%%%%%%%%%%%%%Original development%%%%%%%%%%%%%%%%%%%%%%%%%
 *  Copyright (c) 2016 Inria and University Pierre and <NAME>
 *  All rights reserved.
 * %%%%%%%%%%%%%%%%%%%%%%%Modifications and further additions%%%%%%%%%%
 *  <NAME>, 2017, within FELTOR and EXBLAS licenses
 */
/**
 *  @file exdot_cuda.cuh
 *  @brief CUDA version of exdot
 *
 *  @authors
 *    Developers : \n
 *        <NAME>  -- <EMAIL> \n
 *        <NAME> -- <EMAIL> \n
 *        <NAME> -- <EMAIL>
 */
#pragma once
#include "thrust/device_vector.h"
#include "accumulate.cuh"

namespace exblas{
///@cond
namespace gpu{



////////////////////////////////////////////////////////////////////////////////
// Auxiliary functions
////////////////////////////////////////////////////////////////////////////////
__device__
double TwoProductFMA(double a, double b, double *d) {
    double p = a * b;
    *d = __fma_rn(a, b, -p);
    return p;
}

__device__
double KnuthTwoSum(double a, double b, double *s) {
    double r = a + b;
    double z = r - a;
    *s = (a - (r - z)) + (b - z);
    return r;
}


template<uint NBFPE, uint WARP_COUNT>
__global__ void ExDOT(
    int64_t *d_PartialSuperaccs,
    const double *d_a,
    const double *d_b,
    const uint NbElements
) {
    __shared__ int64_t l_sa[WARP_COUNT * BIN_COUNT]; //shared variables live for a thread block (39 rows, 16 columns!)
    int64_t *l_workingBase = l_sa + (threadIdx.x & (WARP_COUNT - 1)); //the bitwise & with 15 is a modulo operation: threadIdx.x % 16
    //Initialize superaccs
    for (uint i = 0; i < BIN_COUNT; i++)
        l_workingBase[i * WARP_COUNT] = 0;
    __syncthreads();

    //Read data from global memory and scatter it to sub-superaccs
    double a[NBFPE] = {0.0};
    for(uint pos = blockIdx.x*blockDim.x+threadIdx.x; pos < NbElements; pos += gridDim.x*blockDim.x) {
        double r = 0.0;
        double x = TwoProductFMA(d_a[pos], d_b[pos], &r);
        //double x = d_a[pos]*d_b[pos];//ATTENTION: if we write it like this, cpu compiler might generate an fma from this while nvcc does not...

        #pragma unroll
        for(uint i = 0; i != NBFPE; ++i) {
            double s;
            a[i] = KnuthTwoSum(a[i], x, &s);
            x = s;
        }
        if (x != 0.0) {
            Accumulate(l_workingBase, x, WARP_COUNT);
            // Flush FPEs to superaccs
            #pragma unroll
            for(uint i = 0; i != NBFPE; ++i) {
                Accumulate(l_workingBase, a[i], WARP_COUNT);
                a[i] = 0.0;
            }
        }

        if (r != 0.0) {//add the rest r in the same manner
            #pragma unroll
            for(uint i = 0; i != NBFPE; ++i) {
                double s;
                a[i] = KnuthTwoSum(a[i], r, &s);
                r = s;
            }
            if (r != 0.0) {
                Accumulate(l_workingBase, r, WARP_COUNT);
                // Flush FPEs to superaccs
                #pragma unroll
                for(uint i = 0; i != NBFPE; ++i) {
                    Accumulate(l_workingBase, a[i], WARP_COUNT);
                    a[i] = 0.0;
                }
            }
        }
    }
    //Flush FPEs to superaccs
    #pragma unroll
    for(uint i = 0; i != NBFPE; ++i)
        Accumulate(l_workingBase, a[i], WARP_COUNT);
    __syncthreads();

    //Merge sub-superaccs into work-group partial-accumulator ( ATTENTION: PartialSuperacc is transposed!)
    uint pos = threadIdx.x;
    if(pos < WARP_COUNT) {
        int imin = IMIN, imax = IMAX;
        Normalize( l_workingBase, imin, imax, WARP_COUNT);
    }
    __syncthreads();

    if (pos < BIN_COUNT) {
        int64_t sum = 0;

        for(uint i = 0; i < WARP_COUNT; i++)
            sum += l_sa[pos * WARP_COUNT + i];

        d_PartialSuperaccs[blockIdx.x * BIN_COUNT + pos] = sum;
    }

    __syncthreads();
    if (pos == 0) {
        int imin = IMIN, imax = IMAX;
        Normalize(&d_PartialSuperaccs[blockIdx.x * BIN_COUNT], imin, imax);
    }
}

template<uint NBFPE, uint WARP_COUNT>
__global__ void ExDOT(
    int64_t *d_PartialSuperaccs,
    const double *d_a,
    const double *d_b,
    const double *d_c,
    const uint NbElements
) {
    __shared__ int64_t l_sa[WARP_COUNT * BIN_COUNT]; //shared variables live for a thread block (39 rows, 16 columns!)
    int64_t *l_workingBase = l_sa + (threadIdx.x & (WARP_COUNT - 1)); //the bitwise & with 15 is a modulo operation: threadIdx.x % 16
    //Initialize superaccs
    for (uint i = 0; i < BIN_COUNT; i++)
        l_workingBase[i * WARP_COUNT] = 0;
    __syncthreads();

    //Read data from global memory and scatter it to sub-superaccs
    double a[NBFPE] = {0.0};
    for(uint pos = blockIdx.x*blockDim.x+threadIdx.x; pos < NbElements; pos += gridDim.x*blockDim.x) {
        //double x2 = d_a[pos]*d_c[pos]*d_b[pos];
        //double r  = 0.0, r2 = 0.0;
        //double x  = TwoProductFMA(d_a[pos], d_b[pos], &r);
        //double x2 = TwoProductFMA(x , d_c[pos], &r2);
        double x1 = __fma_rn( d_a[pos], d_b[pos], 0);
        double x2 = __fma_rn( x1      , d_c[pos], 0);

        if( x2 != 0.0) {//accumulate x2
            #pragma unroll
            for(uint i = 0; i != NBFPE; ++i) {
                double s;
                a[i] = KnuthTwoSum(a[i], x2, &s);
                x2 = s;
            }
            if (x2 != 0.0) {
                Accumulate(l_workingBase, x2, WARP_COUNT);
                // Flush FPEs to superaccs
                #pragma unroll
                for(uint i = 0; i != NBFPE; ++i) {
                    Accumulate(l_workingBase, a[i], WARP_COUNT);
                    a[i] = 0.0;
                }
            }
        }
        //if (r2 != 0.0) {//add the rest r2
        //    #pragma unroll
        //    for(uint i = 0; i != NBFPE; ++i) {
        //        double s;
        //        a[i] = KnuthTwoSum(a[i], r2, &s);
        //        r2 = s; //error was here r = s
        //    }
        //    if (r2 != 0.0) { //error was here r != 0.0
        //        Accumulate(l_workingBase, r2, WARP_COUNT);
        //        // Flush FPEs to superaccs
{'fraction_non_alphanumeric': 0.16692603816856416, 'fraction_numerical': 0.01769186665574576, 'mean_word_length': 1.6876513317191284, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3294477', 'n_tokens_mistral': 4489, 'n_tokens_neox': 4265, 'n_words': 1156}