2

I am working on optimization of cuda program. So I first started with optimization of matrix multiplication program. Threading scheme which I have used for parallelization is Blocksize(1, 1),Gridsize(N ,N). I am using surface memory for memory optimization purpose(as use of shared memory is not possible for this threading scheme). When I compare the time after and before optimization, I found that execution takes double time after using surface memory(I have tried with different threading scheme but the problem remains same). From whatever I have read till now, global memory is slower than surface memory. So use of surface memory should take less time.Below I am giving matrix multiplication program with surface memory used. Can somebody tell me what is the problem?

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

//#define N 3

surface < void, 2 > a_surf;
surface < void, 2 > b_surf;
surface < void, 2 > c_surf;

void CUDA_SAFE_CALL(cudaError_t call, int line) {
    switch (call) {
    case cudaSuccess:
        break;
    default:
        printf("ERROR at line :%i.%d' ' %s\n",
            line, call, cudaGetErrorString(call));
        exit(-1);
        break;
    }

}

__global__ void mul(int N) {
    int a, b, c, temp;
    int i;

    unsigned int x = blockIdx.x * blockDim.x + (threadIdx.x);
    unsigned int y = blockIdx.y * blockDim.y + (threadIdx.y);
    if (x < N && y < N) {

        temp = 0;
        for (i = 0; i < N; i++) {
            surf2Dread( & a, a_surf, (x) * 4, i);
            surf2Dread( & b, b_surf, (i) * 4, y);
            temp += a * b;
        }
        c = temp;

        // Write to output surface
        surf2Dwrite(c, c_surf, x * 4, y);
    }
}

int main() {
    int N = 100;
    int a[N][N], b[N][N], c[N][N];
    int i, j;
    int temp;
    clock_t t1, t2;
    cudaArray * da, * db, * dc;
    cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc < int > ();

    dim3 dimBlock(1, 1);
    dim3 dimGrid(N, N);

    temp = 0;
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
            a[i][j] = ++temp;

    temp = 0;
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
            b[i][j] = ++temp;

    CUDA_SAFE_CALL(cudaMallocArray( & da, & channelDesc, N, N, cudaArraySurfaceLoadStore), __LINE__);
    CUDA_SAFE_CALL(cudaMallocArray( & db, & channelDesc, N, N, cudaArraySurfaceLoadStore), __LINE__);
    CUDA_SAFE_CALL(cudaMallocArray( & dc, & channelDesc, N, N, cudaArraySurfaceLoadStore), __LINE__);

    int s = N * N * sizeof(int);

    CUDA_SAFE_CALL(cudaMemcpyToArray(da, 0, 0, a, s, cudaMemcpyHostToDevice), __LINE__);
    CUDA_SAFE_CALL(cudaMemcpyToArray(db, 0, 0, b, s, cudaMemcpyHostToDevice), __LINE__);

    CUDA_SAFE_CALL(cudaBindSurfaceToArray(a_surf, da), __LINE__);
    CUDA_SAFE_CALL(cudaBindSurfaceToArray(b_surf, db), __LINE__);
    CUDA_SAFE_CALL(cudaBindSurfaceToArray(c_surf, dc), __LINE__);

    t1 = clock();
    mul <<<dimGrid, dimBlock>>> (N);
    t2 = clock();

    CUDA_SAFE_CALL(cudaMemcpyFromArray(c, dc, 0, 0, s, cudaMemcpyDeviceToHost), __LINE__);

    double t3 = (double) t2 - (double) t1;
    t3 = t3 / CLOCKS_PER_SEC;

    printf("\n CUDA time :%lf", t3);

    CUDA_SAFE_CALL(cudaFreeArray(da), __LINE__);
    CUDA_SAFE_CALL(cudaFreeArray(db), __LINE__);
    CUDA_SAFE_CALL(cudaFreeArray(dc), __LINE__);
}
8
  • 2
    You should note that the kernel calls are asynchronous. You are only timing the elapsed time of kernel launch.
    – pQB
    Apr 1, 2016 at 10:58
  • What operating system are you running this code on?
    – talonmies
    Apr 1, 2016 at 13:01
  • 1
    In addition to your timing issue, you'll get horribly slow performance this way: dim3 dimBlock(1, 1); That "threading scheme" you have chosen will leave approximately 97% of the GPU horsepower unused. Your optimization priorities are incorrect. Apr 1, 2016 at 14:57
  • @talonmies I am using ubuntu 14.04
    – mns
    Apr 1, 2016 at 15:10
  • OK so your time measurement is totally broken as well. POSIX clock() doesn't do what you think it does
    – talonmies
    Apr 1, 2016 at 15:13

1 Answer 1

6

Optimizing for caches is not a trivial matter. So such trivialized generalization as this:

From whatever I have read till now, global memory is slower than surface memory. So use of surface memory should take less time.

are simply so broad as to be incorrect, in my opinion. It will be frequently true, but not always true. The specifics matter, and proper programming practice matters, too.

Surface memory is nothing more than global memory with an intervening cache. But global memory (on all GPUs supported by current CUDA versions) already has support from L2 (and in some cases L1) cache(s).

The code you have proposed for test/comparison has a number of issues that I would point out:

  1. Your timing methodology is incorrect. This:

    t1 = clock();
    mul <<<dimGrid, dimBlock>>> (N);
    t2 = clock();
    

    will time the duration of the kernel launch not the duration of the kernel execution. So this is almost never the correct way to time things. We can fix this by putting a cudaDeviceSynchronize(); call in the timing region, to force completion of the kernel before timing closure.

  2. This is a particularly bad construct if you are interested in performance:

    dim3 dimBlock(1, 1);
    

    because 31 out of every 32 threads in every GPU warp will be inactive, you are leaving 31/32 of the performance of the GPU unused. This has wide-ranging implications. I have no interest in studying the performance of such a scenario, and you shouldn't either (as it is not reflective of real-world performance on well-written codes), unless you are interested in microbenchmarking (not comparative benchmarking). So your code should be fixed to handle at least 32, and ideally 256 or more threads per block.

  3. You've provided no "global memory" comparison case. So I shall provide one.

  4. You've not stated many other factors important for comparative benchmarking, or perf analysis, such as the GPU and platform you are running on, as well as the compile command.

  5. In my opinion, the problem size is too small. A matrix multiply of 100x100 matrices is on the edge of a code that could reasonably occupy the GPU, or test it's performance limits. So I shall make the problem size larger.

With respect to the problem size argument, this is important for the cache discussion. First of all, the surface cache tends to be a spatially-optimized cache, whereas the ordinary L1 and L2 caches are linearly (cache-line) optimized. For very large 2D problems, the surface cache might give better behavior than the L2. But for very small problems, the difference will be less pronounced. Secondly, the surface cache is in addition to the L1 and L2 caches, so a good optimization strategy is to funnel some data through L1 and L2, and other data through surface, to maximize the available cache lines. In fact, since your input matrices are read-only, a further optimization might be to use textures rather than surface for those. But from a contrary point of view, if my problem is so small as to completely fit in the L2 cache, then the surface cache is not likely to give a significant improvement. Your original problem size included 3 matrices of 100x100 int quantities, so about 40Kbytes each, or 120K bytes total. This problem size will fit in the L2 cache of most GPUs. By increasing the problem size (as we shall see - to about 12MB total) we can severely handicap the global-memory-only case.

Here's a code and fully worked example, that has been modified to address most of the above issues. When I run this code on my Quadro5000 GPU on CUDA 7.5/Fedora 20, I observe the surface case to be about 8x faster than the global memory case:

$ cat t1129.cu
#include <stdio.h>
#include <iostream>

typedef int mytype;
const int blk_dim=16;

#define my_N 1000
#define A_VAL 1
#define B_VAL 2

surface < void, 2 > a_surf;
surface < void, 2 > b_surf;
surface < void, 2 > c_surf;

void CUDA_SAFE_CALL(cudaError_t call, int line) {
    switch (call) {
    case cudaSuccess:
        break;
    default:
        printf("ERROR at line :%i.%d' ' %s\n",
            line, call, cudaGetErrorString(call));
        exit(-1);
        break;
    }

}

#ifdef USE_GLOBAL
__global__ void mul(const mytype * __restrict__ d_a, const mytype * __restrict__ d_b, mytype * __restrict__ d_c, const int N)
#else
__global__ void mul(const int N)
#endif
{
    mytype a, b, c, temp;
    int i;

    unsigned int x = blockIdx.x * blockDim.x + (threadIdx.x);
    unsigned int y = blockIdx.y * blockDim.y + (threadIdx.y);
    if (x < N && y < N) {

        temp = 0;
        for (i = 0; i < N; i++) {
#ifdef USE_GLOBAL
            a = d_a[x*N+i];
            b = d_b[i*N+y];
#else
            surf2Dread( & a, a_surf, (x) * sizeof(mytype), i);
            surf2Dread( & b, b_surf, (i) * sizeof(mytype), y);
#endif
            temp += a * b;
        }
        c = temp;
#ifdef USE_GLOBAL
        d_c[x*N+y] = c;
#else
        // Write to output surface
        surf2Dwrite(c, c_surf, x * sizeof(mytype), y);
#endif
    }
}

int main() {
    const int N = my_N;
    mytype *a, *b, *c, *d_a, *d_b, *d_c;
    int i, j;
    clock_t t1, t2;
    cudaArray * da, * db, * dc;
    cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc < mytype > ();

    dim3 dimBlock(blk_dim, blk_dim);
    dim3 dimGrid((N+dimBlock.x-1)/dimBlock.x, (N+dimBlock.y-1)/dimBlock.y);
    int s = N * N * sizeof(mytype);

    a = (mytype *)malloc(s);
    b = (mytype *)malloc(s);
    c = (mytype *)malloc(s);

    CUDA_SAFE_CALL(cudaMalloc(&d_a, s), __LINE__);
    CUDA_SAFE_CALL(cudaMalloc(&d_b, s), __LINE__);
    CUDA_SAFE_CALL(cudaMalloc(&d_c, s), __LINE__);

    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
            a[i*N+j] = A_VAL;

    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
            b[i*N+j] = B_VAL;

    CUDA_SAFE_CALL(cudaMallocArray( & da, & channelDesc, N, N, cudaArraySurfaceLoadStore), __LINE__);
    CUDA_SAFE_CALL(cudaMallocArray( & db, & channelDesc, N, N, cudaArraySurfaceLoadStore), __LINE__);
    CUDA_SAFE_CALL(cudaMallocArray( & dc, & channelDesc, N, N, cudaArraySurfaceLoadStore), __LINE__);


    CUDA_SAFE_CALL(cudaMemcpyToArray(da, 0, 0, a, s, cudaMemcpyHostToDevice), __LINE__);
    CUDA_SAFE_CALL(cudaMemcpyToArray(db, 0, 0, b, s, cudaMemcpyHostToDevice), __LINE__);

    CUDA_SAFE_CALL(cudaBindSurfaceToArray(a_surf, da), __LINE__);
    CUDA_SAFE_CALL(cudaBindSurfaceToArray(b_surf, db), __LINE__);
    CUDA_SAFE_CALL(cudaBindSurfaceToArray(c_surf, dc), __LINE__);

#ifdef USE_GLOBAL
    CUDA_SAFE_CALL(cudaMemcpy(d_a, a, s, cudaMemcpyHostToDevice), __LINE__);
    CUDA_SAFE_CALL(cudaMemcpy(d_b, b, s, cudaMemcpyHostToDevice), __LINE__);
#endif
    t1 = clock();
#ifdef USE_GLOBAL
    mul <<<dimGrid, dimBlock>>> (d_a, d_b, d_c, N);
#else
    mul <<<dimGrid, dimBlock>>> (N);
#endif
    cudaDeviceSynchronize();
    t2 = clock();

    CUDA_SAFE_CALL(cudaMemcpyFromArray(c, dc, 0, 0, s, cudaMemcpyDeviceToHost), __LINE__);
#ifdef USE_GLOBAL
    CUDA_SAFE_CALL(cudaMemcpy(c, d_c, s, cudaMemcpyDeviceToHost), __LINE__);
#endif

    double t3 = (double) t2 - (double) t1;
    t3 = t3 / CLOCKS_PER_SEC;

    printf("\n CUDA time :%lf\n", t3);
    for (i=0; i < N*N; i++)
      if(c[i] != A_VAL*B_VAL*N) {std::cout << "mismatch at: " << i << ", was: " << c[i] << " should be: " << A_VAL*B_VAL*N << std::endl;  return 1;}

    CUDA_SAFE_CALL(cudaFreeArray(da), __LINE__);
    CUDA_SAFE_CALL(cudaFreeArray(db), __LINE__);
    CUDA_SAFE_CALL(cudaFreeArray(dc), __LINE__);
    std::cout << "Success!"  << std::endl;
    return 0;
}
[bob@cluster1 misc]$ nvcc -O3 -o t1129 t1129.cu
[bob@cluster1 misc]$ ./t1129

 CUDA time :0.028771
Success!
$ nvcc -O3 -DUSE_GLOBAL -o t1129 t1129.cu
$ ./t1129

 CUDA time :0.243635
Success!
$

As a final note, there are many other optimizations we could talk about, which would probably shift the comparison one way or the other. But if you actually want to do fast matrix multiply operations, you should use CUBLAS. You should not write your own matrix multiply routines.

2
  • Thank you Robert. Use of cudaDeviceSynchronize() worked, its use shows less time for surface memory than global memory. But it increases total time required for execution. I have GeForce GT 740M GPU, CUDA 7.0 ,Ubuntu 14.04.
    – mns
    Apr 3, 2016 at 10:42
  • When I run this test sample today on my GTX 1650, surfaces are only 2x faster than global memory. However, more interestingly, if I switch up the indexing of the linear memory (to a = d_a[y*N+i]; b = d_b[i*N+x]; and similar with d_c), the linear-memory approach is 2x the speed of surfaces. This change makes better use of coalescing, but is otherwise not an optimisation per se (it is how I would have personally indexed it anyway), so I find it interesting that surfaces didn't help at all there.
    – lxop
    Feb 28, 2023 at 2:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.