def atomic_add(API):
    """
    Return the atomic_add for the given API. 
    Overcome the missing atomic_add_float for OpenCL-1.2. 
    Note: will be checked if OpenCL 2.0 provided by all GPU vendors. 
    """
    ocl_add = """
    // #pragma OPENCL EXTENSION cl_khr_int64_base_atomics : enable
    KERNEL void atomic_add_float( 
            GLOBAL_MEM float *ptr, 
            const float temp) 
    {
    // The work-around of AtomicAdd for float
    // lockless add *source += operand 
    // Caution!!!!!!! Use with care! You have been warned!
    // http://simpleopencl.blogspot.com/2013/05/atomic-operations-and-floats-in-opencl.html
    // Source: https://github.com/clMathLibraries/clSPARSE/blob/master/src/library/kernels/csrmv_adaptive.cl
        union {
            unsigned int intVal;
            float floatVal;
        } newVal;
        
        union {
            unsigned int intVal;
            float floatVal;
        } prevVal;
        
        do {
            prevVal.floatVal = *ptr;
            newVal.floatVal = prevVal.floatVal + temp;
        } while (atomic_cmpxchg((__global volatile unsigned int *)ptr, prevVal.intVal, newVal.intVal) != prevVal.intVal);
    };  
    
    KERNEL void atomic_add_float2_old( 
            GLOBAL_MEM float2 *ptr, 
            const float2 temp) 
    {
    atomic_add_float( 
            (GLOBAL_MEM float*)ptr, 
            temp.x);
    atomic_add_float( 
    (GLOBAL_MEM float*)ptr+1, 
    temp.y);
    };  
    
    KERNEL void atomic_add_float2( 
            GLOBAL_MEM float2 *ptr, 
            const float2 temp, 
            GLOBAL_MEM float2 *res) 
        { // Add atomic two sum for real/imaginary parts
          // *res array stores the error
        union {
            unsigned int intVal;
            float floatVal;
        } newVal;
        
        union {
            unsigned int intVal;
            float floatVal;
        } prevVal;
        float x;
        float val;
        float newacc;
        float r;
            
        __global const float * realptr = ( __global float *)ptr;
        __global const float * imagptr = ( __global float *)ptr + 1;
        __global const float * realres = ( __global float *)res;
        __global const float * imagres = ( __global float *)res + 1;  
          
        do { // atomic add of the value
            prevVal.floatVal = *realptr ;
            newVal.floatVal = prevVal.floatVal + temp.x;
        } while (atomic_cmpxchg((__global volatile unsigned int *)realptr, prevVal.intVal, newVal.intVal) != prevVal.intVal);
        
        x = prevVal.floatVal;
        val = temp.x;
        
        if (fabs(x) < fabs(val))
        {
            const float swap = x;
            x = val;
            val = swap;
        }
        newacc = x + val;
        r = val - (newacc - x);    // Recover rounding error using Fast2Sum, assumes |oldacc|>=|val|
                
        do { // atomic_add(realres, r)
            prevVal.floatVal = *realres ;
            newVal.floatVal = prevVal.floatVal + r;
        } while (atomic_cmpxchg((__global volatile unsigned int *)realres, prevVal.intVal, newVal.intVal) != prevVal.intVal);
        
        // End of real part
        
        // Now do the imaginary part
        
        do {
            prevVal.floatVal = *imagptr;
            newVal.floatVal = prevVal.floatVal + temp.y;
        } while (atomic_cmpxchg((__global volatile unsigned int *)imagptr, prevVal.intVal, newVal.intVal) != prevVal.intVal);
        
                x = prevVal.floatVal; // this copy the old value to x
                val = temp.y; // this copy the  real part to val
                
                if (fabs(x) < fabs(val)) 
                {
                    const float swap = x;
                    x = val;
                    val = swap;
                }
                newacc = x + val;
                r = val - (newacc - x);    // Recover rounding error using Fast2Sum, assumes |oldacc|>=|val|
                
                do {
                    prevVal.floatVal = *imagres ;
                    newVal.floatVal = prevVal.floatVal + r;
                } while (atomic_cmpxchg((__global volatile unsigned int *)imagres, prevVal.intVal, newVal.intVal) != prevVal.intVal);
       
    };  
    
    """
    cuda_add = """
    
    __device__ void atomic_add_float_twosum(
            GLOBAL_MEM float * address,
            float val,
            GLOBAL_MEM float * residue)
    {
        
        // https://devtalk.nvidia.com/default/topic/817899/atomicadd-kahan-summation/
        // By Sylvain Collange
        // Also add the branch clsparse
        
        float x = atomicAdd(address, val);
        
        if (fabs(x) < fabs(val))
        {
            const float swap = x;
            x = val;
            val = swap;
        }
        float newacc = x + val;
        float r = val - (newacc - x);    // Recover rounding error using Fast2Sum, assumes |oldacc|>=|val|
        atomicAdd(residue, r);    // Accumulate error in residue
        
    }    
    
    __device__ void atomic_add_float( 
            GLOBAL_MEM float *ptr, 
            const float temp) 
    { // Wrapper around CUDA atomicAdd();
    atomicAdd(ptr, temp); 
    };   
    __device__ void atomic_add_float2( 
            GLOBAL_MEM float2 * ptr, 
            const float2 temp,
            GLOBAL_MEM float2 *res) 
    { // Wrapper around CUDA atomicAdd();
    // atomicAdd((float*)ptr, temp.x);
    // atomicAdd((float*)ptr+1, temp.y);
    atomic_add_float_twosum((float*)ptr, temp.x, (float*)res);
    atomic_add_float_twosum((float*)ptr+1, temp.y, (float*)res+1);
    };       
    
    """
    if API == 'cuda':
        code_text = cuda_add
    elif API == 'ocl':
        code_text = ocl_add
    return code_text
