source: CIVL/examples/cuda/matMult1.cvl@ e2877ba

#include <civlc.cvh>
/***********************************************************************
* FILENAME:  MM.cu
*            Matrix Multiplication
*            Matrix operands have row-major order.
*
* C = A * B
* Multiplies two square matrices (NxN * NxN).
* Matrix values have type double.
*
* A simple CUDA program has a basic workflow:
*     1)  Initialize matrix operands as double-precision arrays on host (CPU).
*     2)  Copy operands from host memory to GPU memory.
*     3)  Apply matrix operaton to operands on GPU
*     4)  Copy result from GPU memory to host memory.
*
*
* CUDA C Programming Guide Version 4.2 (3.2.3, p.22):
* http://developer.download.nvidia.com/compute/DevZone/docs/html/C/doc/CUDA_C_Programming_Guide.pdf
*
* MM with linearized matrix operands:
* http://www.hpcwire.com/hpcwire/2008-10-08/compilers_and_more_programming_gpus_today.html
*
*************************************************************************/
// online source: https://www.rcac.purdue.edu/userinfo/resources/carter/compile/MM.cu

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cuda.h>

$input int N;               /* size of square matrix                   */
$input int TILE_WIDTH;


/* MM kernel using global (not shared) memory.                          */
void _kernel_myMM_global (dim3 gridDim, dim3 blockDim, cudaStream_t s, const double * const A, const double * const B, double *C, int width) {

    void _kernel (_kernelInstance *this, cudaEvent_t e) {
 
        _waitInQueue(this, e);

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

            void _thread(uint3 threadIdx) {

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

                /* Get row and column from block and thread IDs                     */
                int row = (blockDim.y*blockIdx.y) + threadIdx.y;
                int col = (blockDim.x*blockIdx.x) + threadIdx.x;

                /* Initialize result of one element which one thread computes.      */
                double result=0.0;

                /* Compute one element of the matrix product.                       */
                for (int i = 0; i < width; ++i)
                    result += A[row*width + i] * B[i*width + col];

                /* Store the result of one matrix element in matrix C.              */
                C[row * width + col] = result;
                $barrier_destroy(_b);
            }
            _runProcs(blockDim, _thread);
            $gbarrier_destroy(_block_barrier);
        }
        _runProcs(gridDim, _block);
        _kernelFinish(this);
    }
    _enqueueKernel(s, _kernel);
}


/* MM kernel using shared memory.                                       */
void _kernel_myMM_shared (dim3 gridDim, dim3 blockDim, cudaStream_t s, const double * const A, const double * const B, double* C, int width) {

    void _kernel (_kernelInstance *this, cudaEvent_t e) {
 
        _waitInQueue(this, e);

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

            double A_shared[TILE_WIDTH][TILE_WIDTH];
            double B_shared[TILE_WIDTH][TILE_WIDTH];

            void _thread(uint3 threadIdx) {

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

                int bx = blockIdx.x;  int by = blockIdx.y;
                int tx = threadIdx.x; int ty = threadIdx.y;

                /* Identify the row and column of the C element to work on.         */
                int row = by * TILE_WIDTH + ty;
                int col = bx * TILE_WIDTH + tx;

                double result = 0.0;

                /* Loop over the A and B tiles required to compute the C element.   */
                for (int phase = 0; phase < width/TILE_WIDTH; ++phase) {
                    /* Shared effort: loading of A and B tiles into shared memory.  */

                    A_shared[ty][tx] = A[row*width + (phase*TILE_WIDTH + tx)];

                    B_shared[ty][tx] = B[col + (phase*TILE_WIDTH + ty)*width]; 

                    $barrier_call(_b); 

                    for (int k = 0; k < TILE_WIDTH; ++k) {
                        result += A_shared[ty][k] * B_shared[k][tx];
                    }

                    $barrier_call(_b); 

                }
                C[row*width+col] = result;
                $barrier_destroy(_b);
            }
            _runProcs(blockDim, _thread);
            $gbarrier_destroy(_block_barrier);
        }
        _runProcs(gridDim, _block);
        _kernelFinish(this);
    }
    _enqueueKernel(s, _kernel);
}


/************************************************************************/
/************************************************************************/
/************************************************************************/ 
 

int main (int argc, char** argv) {

        int _main ( void ) {

                /* Set device based on input from command line            */
        /*
                if (argc > 1) {
                        if (cudaSetDevice(atoi(argv[1])) != cudaSuccess) {
                                int num_devices;
                                cudaGetDeviceCount(&num_devices);
                                fprintf(stderr, "Error initializing device %s,\
         device value must be 0-%d\n", argv[1], (num_devices-1));
                                return 0;
                        }
                } else {
                        fprintf(stderr, "Usage: %s gpu_device\n", argv[0]);
                        return 0;
                }
        */

                /* Declare CPU arrays.                                              */
                double A[N*N],B[N*N],C[N*N];       /* linearized CPU double arrays  */ 
        double cpuResult[N*N], gpuGlobalResult[N*N], gpuSharedResult[N*N];
                int r,c;

                /* Declare GPU arrays.                                              */
                double *G_A,*G_B,*G_C;             /* linearized GPU double arrays  */
                size_t size_a,size_b,size_c;       /* size of linearized array in bytes */
                size_a = size_b = size_c = N*N;

                /* Setup a clock.                                                   */
                cudaEvent_t start, stop;
                float CPU_elapsedtime, GPU_global_elapsedtime, GPU_shared_elapsedtime;
                cudaEventCreate(&start);
                cudaEventCreate(&stop);




                /* 1)  Initialize matrix operands as double-precision arrays on host (CPU). */
                for (r=0;r<N;++r)
                for (c=0;c<N;++c) {
                        A[r*N+c] = 1.0;
                        B[r*N+c] = 1.0;
                        C[r*N+c] = 0.0;
                }


        /*-----------------------------------------------------------------------*/

                /* MM on a CPU.                                                      */
                cudaEventRecord(start,0);
                for (int r = 0; r < N; ++r )
                for (int c = 0; c < N; ++c )
                for (int k = 0; k < N; ++k )
                        C[r*N+c] += A[r*N+c] * B[k*N+c];
                cudaEventRecord(stop,0);
                cudaEventSynchronize(stop);
                cudaEventElapsedTime(&CPU_elapsedtime,start,stop);
        memcpy(cpuResult, C, size_c*sizeof(double));
                printf("                                                            speedup\n");
                printf("                                                            -------\n");
                printf("Elapsed time in CPU:                   %7.1f milliseconds\n", CPU_elapsedtime); 
        /*-----------------------------------------------------------------------*/

                /* MM on Global Memory of GPGPU.                                     */
                cudaEventRecord(start,0);

                /* 2)  Copy operands from CPU memory to GPGPU memory.                */
                //cudaMalloc((void**)&G_A,size_a*sizeof(double));  /* alloc A in GPGPU */
        G_A = (double*)$malloc($root, size_a*sizeof(double));
                //cudaMalloc((void**)&G_B,size_b*sizeof(double));  /* alloc B in GPGPU */
        G_B = (double*)$malloc($root, size_b*sizeof(double));
                //cudaMalloc((void**)&G_C,size_c*sizeof(double));  /* alloc C in GPGPU */
        G_C = (double*)$malloc($root, size_c*sizeof(double));
                cudaMemcpy(G_A,A,size_a*sizeof(double),cudaMemcpyHostToDevice);
                cudaMemcpy(G_B,B,size_b*sizeof(double),cudaMemcpyHostToDevice);

                /* 3)  Apply matrix operation to operands on GPGPU                   */
                /*     There is no partial final block in this example.              */
                dim3 block = { TILE_WIDTH, TILE_WIDTH, 1 };      /* using a 2D block: 16,16,1 */
                dim3 grid = { N/TILE_WIDTH, N/TILE_WIDTH, 1 };   /* as many 16x16-thread blocks as needed: */
                _kernel_myMM_global(grid, block, 0, G_A, G_B, G_C, N);  /* grid(16,16,1)  */ 

                /* 4)  Copy result from GPGPU memory to CPU memory.                  */
                cudaMemcpy(C,G_C,size_c*sizeof(double),cudaMemcpyDeviceToHost);

                /* Deallocate memory on GPGPU.                                       */
                cudaFree(G_A);
                cudaFree(G_B);
                cudaFree(G_C);

                cudaEventRecord(stop,0);
                cudaEventSynchronize(stop);
                cudaEventElapsedTime(&GPU_global_elapsedtime,start,stop);
        memcpy(gpuGlobalResult, C, size_c*sizeof(double));
                printf("Elapsed time in GPU (global memory):   %7.1f milliseconds  %5.1f\n",
                           GPU_global_elapsedtime,CPU_elapsedtime/GPU_global_elapsedtime);
        /*
                printf("\nGLOBAL MEMORY:\n");
                for (r=0;r<10;++r)
                for (c=0;c<10;++c) {
                        printf("%2d,%2d   %g\n", r,c,C[r*N+c]);
                }
        */
        /*-----------------------------------------------------------------------*/

                /* MM on Shared Memory of GPGPU.                                     */
                cudaEventRecord(start,0);

                /* 2)  Copy operands from CPU memory to GPGPU memory.                */
        G_A = (double*)$malloc($root, size_a*sizeof(double));
        G_B = (double*)$malloc($root, size_b*sizeof(double));
        G_C = (double*)$malloc($root, size_c*sizeof(double));
                cudaMemcpy(G_A,A,size_a*sizeof(double),cudaMemcpyHostToDevice);
                cudaMemcpy(G_B,B,size_b*sizeof(double),cudaMemcpyHostToDevice);

                /* 3)  Apply matrix operation to operands on GPGPU                   */
                /*     There is not partial final block in this example.             */
                /*     Use the same grid and block from the previous case.           */
        printf("a\n");
                _kernel_myMM_shared(grid, block, 0, G_A, G_B, G_C, N);
        printf("b\n");

                /* 4)  Copy result from GPGPU memory to CPU memory.                  */
                cudaMemcpy(C,G_C,size_c*sizeof(double),cudaMemcpyDeviceToHost);

                /* Deallocate memory on GPGPU.                                       */
                cudaFree(G_A);
                cudaFree(G_B);
                cudaFree(G_C);

                cudaEventRecord(stop,0);
                cudaEventSynchronize(stop);
                cudaEventElapsedTime(&GPU_shared_elapsedtime,start,stop);
        memcpy(gpuSharedResult, C, size_c*sizeof(double));
                printf("Elapsed time in GPU (shared memory):   %7.1f milliseconds  %5.1f\n",
                           GPU_shared_elapsedtime,CPU_elapsedtime/GPU_shared_elapsedtime);
        /*
                printf("\nSHARED MEMORY:\n");
                for (r=0;r<10;++r)
                for (c=0;c<10;++c) {
                        printf("%2d,%2d   %g\n", r,c,C[r*N+c]);
                }
        */
        /*-----------------------------------------------------------------------*/ 

                /* Deallocate the clock.                                             */
                cudaEventDestroy(start);
                cudaEventDestroy(stop);
        //$assert(($forall {int i | i >= 0 && i < N*N} cpuResult[i] == gpuGlobalResult[i] && cpuResult[i] == gpuSharedResult[i]));
        for (int i = 0; i < N*N; i ++) {
            printf("i = %d\n", i);
            printf("cpu = %f, gpuGlobal = %f, gpuShared = %f\n", cpuResult[i], gpuGlobalResult[i], gpuSharedResult[i]);
            $assert((cpuResult[i] == gpuGlobalResult[i]));
            $assert((cpuResult[i] == gpuSharedResult[i]));
        }

                return 0; 
        }

    printf("1\n");
    _cudaInit();
    printf("2\n");
    _main();
    printf("3\n");
    _cudaFinalize();
    printf("4\n");
    return 0;
}