source: CIVL/examples/translation/cuda/dot2.cvl@ 764d472

/*
 * 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.h"
#include "civlc-cuda.cvl"
#include "stdio.h"

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

$input int N;
$input int B;
$assume(0 <= N && N <= B);
const int threadsPerBlock = 2;
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 ($proc prevProc) {

                _tryWait(prevProc);

                void _block (uint3 blockIdx) {
                        int _numThreads = blockDim.x * blockDim.y * blockDim.z;
                        $gbarrier _block_barrier = $gbarrier_create($here, _numThreads);

                        float cache[threadsPerBlock];

                        void _thread (uint3 threadIdx) {

                                int _tid = _index(blockDim, threadIdx);
                                $barrier _b = $barrier_create($here, _block_barrier, _tid);

                                int tid = threadIdx.x + blockIdx.x * blockDim.x;
                                int cacheIndex = threadIdx.x;

                                float   temp = 0;
                                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
                                int i = blockDim.x/2;
                                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);
                        }
                        _createProcs(blockDim, _thread);
                        $gbarrier_destroy(_block_barrier);
                }
                _createProcs(gridDim, _block);
        }
        _initKernel(s, _kernel);
}

int main ( void ) {
        _host_scope = $here;

        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(_host_scope,  N*sizeof(float) );
                b = (float*)$malloc(_host_scope,  N*sizeof(float) );
                partial_c = (float*)$malloc(_host_scope, blocksPerGrid*sizeof(float) );

                // allocate the memory on the GPU
                dev_a = (float *)$malloc( _global_scope, N*sizeof(float) );
                dev_b = (float *)$malloc( _global_scope, N*sizeof(float) );
                dev_partial_c = (float *)$malloc( _global_scope, 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 ); // civl version of cudaMemcpy will check memory correctness, then call regular memcpy
                cudaMemcpy( dev_b, b, N*sizeof(float),
                                                                  cudaMemcpyHostToDevice ); 

                dim3 _t1, _t2;
                _t1 = _toDim3(blocksPerGrid);
                _t2 = _toDim3(threadsPerBlock);
                _kernel_dot(_t1, _t2, 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;
        }

        _main();
}