source: CIVL/examples/cuda/dot2.cvl@ 77ad37b

/*
 * Copyright 1993-2010 NVIDIA Corporation.  All rights reserved.
 *
 * NVIDIA Corporation and its licensors retain all intellectual property and 
 * proprietary rights in and to this software and related documentation. 
 * Any use, reproduction, disclosure, or distribution of this software 
 * and related documentation without an express license agreement from
 * NVIDIA Corporation is strictly prohibited.
 *
 * Please refer to the applicable NVIDIA end user license agreement (EULA) 
 * associated with this source code for terms and conditions that govern 
 * your use of this NVIDIA software.
 *
 * This is a version translated to CIVL.
 * 
 */

#include <civlc.cvh>
#include <stdio.h>
#include <cuda.h>

#define imin(a,b) (a<b?a:b)

_Bool isPowerOfTwo(int x) {
        if (x == 1) {
                return $true;
        } else {
                return x % 2 == 0 && isPowerOfTwo(x / 2);
        }
}


// the length of the vectors to dot product
$input int N;
// upper bound on N
$input int B;
$assume((0 <= N && N <= B));
$input int THREADS_PER_BLOCK; // thread number per block: must be a power of 2, due to the while loop at the end of gpuThread();
$input int THREADS_B; 
$assume((1 <= THREADS_PER_BLOCK && THREADS_PER_BLOCK <= THREADS_B));

const int threadsPerBlock = THREADS_PER_BLOCK;
// the number of blocks to create
const int blocksPerGrid =
      imin( 32, (N+threadsPerBlock-1) / threadsPerBlock );

void _kernel_dot(dim3 gridDim, dim3 blockDim, cudaStream_t s, float *a, float *b, float *c ) {

  void _kernel (_kernelInstance *this, cudaEvent_t e) {

    void _block (uint3 blockIdx) {
      // calculate total number of threads
      int _numThreads = blockDim.x * blockDim.y * blockDim.z;
      $gbarrier _block_barrier = $gbarrier_create($here, _numThreads);

      // shared memory for a block
      float cache[threadsPerBlock];

      void _thread (uint3 threadIdx) {

        // get the one dimensional index of this threads
        int _tid = _index(blockDim, threadIdx);
        $barrier _b = $barrier_create($here, _block_barrier, _tid);

        // get this thread's id
        int tid = threadIdx.x + blockIdx.x * blockDim.x;
        int cacheIndex = threadIdx.x;

        float temp = 0;
        int i = blockDim.x/2;

        while (tid < N) {
          temp += a[tid] * b[tid];
          tid += blockDim.x * gridDim.x;
        }
        
        // set the cache values
        cache[cacheIndex] = temp;
        
        // synchronize threads in this block
        
        $barrier_call(_b);

        // for reductions, threadsPerBlock must be a power of 2
        // because of the following code
        while (i != 0) {
          if (cacheIndex < i) {
            cache[cacheIndex] += cache[cacheIndex + i];
          }
          $barrier_call(_b);
          i /= 2;
        }

        if (cacheIndex == 0) {
          c[blockIdx.x] = cache[0];
        }
        

        $barrier_destroy(_b);
      }

      _runProcs(blockDim, _thread);
      $gbarrier_destroy(_block_barrier);
    }

    _waitInQueue(this, e);
    _runProcs(gridDim, _block);
    _kernelFinish(this);
  }
  _enqueueKernel(s, _kernel);
}

int main ( void ) {

  int _main( void ) {
    float   *a, *b, c, *partial_c;
    float   *dev_a, *dev_b, *dev_partial_c;

    // allocate memory on the cpu side
    a = (float*)$malloc($root,  N*sizeof(float) );
    b = (float*)$malloc($root,  N*sizeof(float) );
    partial_c = (float*)$malloc($root, blocksPerGrid*sizeof(float) );

    // allocate the memory on the GPU
    dev_a = (float *)$malloc($root, N*sizeof(float) );
    dev_b = (float *)$malloc($root, N*sizeof(float) );
    dev_partial_c = (float *)$malloc($root, blocksPerGrid*sizeof(float) );

    // fill in the host memory with data
    for (int i=0; i<N; i++) {
      a[i] = i;
      b[i] = i*2;
    }

    // copy the arrays 'a' and 'b' to the GPU
    cudaMemcpy( dev_a, a, N*sizeof(float), cudaMemcpyHostToDevice );
    cudaMemcpy( dev_b, b, N*sizeof(float), cudaMemcpyHostToDevice ); 

    dim3 _t1, _t2;
    _t1 = _toDim3(blocksPerGrid);
    _t2 = _toDim3(threadsPerBlock);
    _kernel_dot(_t1, _t2, 0, dev_a, dev_b, dev_partial_c );
// copy the array 'c' back from the GPU to the CPU
    cudaMemcpy( partial_c, dev_partial_c, blocksPerGrid*sizeof(float), cudaMemcpyDeviceToHost );

    // finish up on the CPU side
    c = 0;
    for (int i=0; i<blocksPerGrid; i++) {
      c += partial_c[i];
    }

    #define sum_squares(x)  (x*(x+1)*(2*x+1)/6)
    printf( "Does GPU value %.6g = %.6g?\n", c,
         2 * sum_squares( (float)(N - 1) ) );
    $assert((c == 2 * sum_squares( (float)(N - 1) ) ));

    // free memory on the gpu side
    cudaFree( dev_a ); // civl versionof cudaFree will check memory correctness, then call regular free
    cudaFree( dev_b );
    cudaFree( dev_partial_c );

    // free memory on the cpu side
    $free( a );
    $free( b );
    $free( partial_c );
    return 0;
  }

        _Bool valid = isPowerOfTwo(THREADS_PER_BLOCK);
        $assume((valid));

  _cudaInit();
  _main();
  _cudaFinalize();
}