Katyusha's blog

1
wget https://developer.download.nvidia.com/compute/cuda/repos/wsl-ubuntu/x86_64/cuda-wsl-ubuntu.pin
2
sudo mv cuda-wsl-ubuntu.pin /etc/apt/preferences.d/cuda-repository-pin-600
3
wget https://developer.download.nvidia.com/compute/cuda/13.1.1/local_installers/cuda-repo-wsl-ubuntu-13-1-local_13.1.1-1_amd64.deb
4
sudo dpkg -i cuda-repo-wsl-ubuntu-13-1-local_13.1.1-1_amd64.deb
5
sudo cp /var/cuda-repo-wsl-ubuntu-13-1-local/cuda-*-keyring.gpg /usr/share/keyrings/
6
sudo apt-get update
7
sudo apt-get -y install cuda-toolkit-13-1

1
export CUDA_HOME=/usr/local/cuda-13.1
2
export PATH=$CUDA_HOME/bin:$PATH
3
export LD_LIBRARY_PATH=/usr/local/cuda/lib64:$LD_LIBRARY_PATH

1
void saxpyCuda(int N, float alpha, float* xarray, float* yarray, float* resultarray) {
2

3
    // must read both input arrays (xarray and yarray) and write to
4
    // output array (resultarray)
5
    int totalBytes = sizeof(float) * 3 * N;
6

7
    // compute number of blocks and threads per block.  In this
8
    // application we've hardcoded thread blocks to contain 512 CUDA
9
    // threads.
10
    const int threadsPerBlock = 512;
11

12
    // Notice the round up here.  The code needs to compute the number
13
    // of threads blocks needed such that there is one thread per
14
    // element of the arrays.  This code is written to work for values
15
    // of N that are not multiples of threadPerBlock.
16
    const int blocks = (N + threadsPerBlock - 1) / threadsPerBlock;
17

18
    // These are pointers that will be pointers to memory allocated
19
    // *one the GPU*.  You should allocate these pointers via
20
    // cudaMalloc.  You can access the resulting buffers from CUDA
21
    // device kernel code (see the kernel function saxpy_kernel()
22
    // above) but you cannot access the contents these buffers from
23
    // this thread. CPU threads cannot issue loads and stores from GPU
24
    // memory!
25
    float* device_x = nullptr;
26
    float* device_y = nullptr;
27
    float* device_result = nullptr;
28

29
    //
30
    // CS149 TODO: allocate device memory buffers on the GPU using cudaMalloc.
31
    //
32
    // We highly recommend taking a look at NVIDIA's
33
    // tutorial, which clearly walks you through the few lines of code
34
    // you need to write for this part of the assignment:
35
    //
36
    // https://devblogs.nvidia.com/easy-introduction-cuda-c-and-c/
37
    //
38
    cudaMalloc(&device_x, totalBytes / 3);
39
    cudaMalloc(&device_y, totalBytes / 3);
40
    cudaMalloc(&device_result, totalBytes / 3);
41
    // start timing after allocation of device memory
42
    double startTime = CycleTimer::currentSeconds();
43

44
    //
45
    // CS149 TODO: copy input arrays to the GPU using cudaMemcpy
46
    //
47

48
    cudaMemcpy(device_x, xarray, totalBytes / 3, cudaMemcpyHostToDevice);
49
    cudaMemcpy(device_y, yarray, totalBytes / 3, cudaMemcpyHostToDevice);
50

51
    // run CUDA kernel. (notice the <<< >>> brackets indicating a CUDA
52
    // kernel launch) Execution on the GPU occurs here.
53
    double kernel_time = CycleTimer::currentSeconds();
54
    saxpy_kernel<<<blocks, threadsPerBlock>>>(N, alpha, device_x, device_y, device_result);
55
    cudaDeviceSynchronize();
56
    kernel_time = CycleTimer::currentSeconds() - kernel_time;
57
    //
58
    // CS149 TODO: copy result from GPU back to CPU using cudaMemcpy
59
    //
60

61
    cudaMemcpy(resultarray, device_result, totalBytes / 3, cudaMemcpyDeviceToHost);
62

63
    // end timing after result has been copied back into host memory
64
    double endTime = CycleTimer::currentSeconds();
65

66
    cudaError_t errCode = cudaPeekAtLastError();
67
    if (errCode != cudaSuccess) {
68
        fprintf(stderr, "WARNING: A CUDA error occured: code=%d, %s\n",
69
    errCode, cudaGetErrorString(errCode));
70
    }
71

72
    double overallDuration = endTime - startTime;
73
    printf("Effective BW of kernel: %.3f ms\t\t[%.3f GB/s]\n", 1000.f * kernel_time, GBPerSec(totalBytes, kernel_time));
74
    printf("Effective BW by CUDA saxpy: %.3f ms\t\t[%.3f GB/s]\n", 1000.f * overallDuration, GBPerSec(totalBytes, overallDuration));
75

76
    //
77
    // CS149 TODO: free memory buffers on the GPU using cudaFree
78
    //
79
    cudaFree(device_x); cudaFree(device_y); cudaFree(device_result);
80
}

1
Found 1 CUDA devices
2
Device 0: NVIDIA GeForce RTX 4060 Laptop GPU
3
   SMs:        24
4
   Global mem: 8188 MB
5
   CUDA Cap:   8.9
6
---------------------------------------------------------
7
Running 3 timing tests:
8
Effective BW of kernel: 13.045 ms               [85.673 GB/s]
9
Effective BW by CUDA saxpy: 153.652 ms          [7.273 GB/s]
10
Effective BW of kernel: 5.262 ms                [212.398 GB/s]
11
Effective BW by CUDA saxpy: 135.543 ms          [8.245 GB/s]
12
Effective BW of kernel: 5.248 ms                [212.972 GB/s]
13
Effective BW by CUDA saxpy: 136.620 ms          [8.180 GB/s]

1
[saxpy serial]:         [52.145] ms     [28.577] GB/s   [3.835] GFLOPS
2
[saxpy ispc]:           [51.940] ms     [28.689] GB/s   [3.851] GFLOPS
3
[saxpy task ispc]:      [31.122] ms     [47.879] GB/s   [6.426] GFLOPS

1
wget https://developer.download.nvidia.com/compute/cuda/repos/wsl-ubuntu/x86_64/cuda-wsl-ubuntu.pin
2
sudo mv cuda-wsl-ubuntu.pin /etc/apt/preferences.d/cuda-repository-pin-600
3
wget https://developer.download.nvidia.com/compute/cuda/12.8.1/local_installers/cuda-repo-wsl-ubuntu-12-8-local_12.8.1-1_amd64.deb
4
sudo dpkg -i cuda-repo-wsl-ubuntu-12-8-local_12.8.1-1_amd64.deb
5
sudo cp /var/cuda-repo-wsl-ubuntu-12-8-local/cuda-*-keyring.gpg /usr/share/keyrings/
6
sudo apt update

1
__global__ void Kernel_up_sweep(int *arrary, int threads, int stride) {
2
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
3
    if (idx < threads) {
4
        int id1 = idx * stride + (stride >> 1), id2 = (idx + 1) * stride;
5
        id1--; id2--;
6
        arrary[id2] = arrary[id1] + arrary[id2];
7
    }
8
}
9

10
__global__ void Kernel_down_sweep(int *arrary, int threads, int stride) {
11
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
12
    if (idx < threads) {
13
        int id1 = idx * stride + (stride >> 1), id2 = (idx + 1) * stride;
14
        id1--; id2--;
15
        int x = arrary[id1], y = arrary[id2];
16
        arrary[id1] = y; arrary[id2] = x + y;
17
    }
18
}
19

20
void exclusive_scan(int* input, int N, int* result)
21
{
22

23
    // CS149 TODO:
24
    //
25
    // Implement your exclusive scan implementation here.  Keep in
26
    // mind that although the arguments to this function are device
27
    // allocated arrays, this is a function that is running in a thread
28
    // on the CPU.  Your implementation will need to make multiple calls
29
    // to CUDA kernel functions (that you must write) to implement the
30
    // scan.
31

32
    int n = nextPow2(N);
33
    for (int stride = 2; stride <= n / 2; stride <<= 1) {
34
        int threads = n / stride, blocks = (threads - 1 + THREADS_PER_BLOCK) / THREADS_PER_BLOCK;
35
        Kernel_up_sweep<<<blocks, THREADS_PER_BLOCK>>>(result, threads, stride);
36
        cudaDeviceSynchronize();
37
    }
38
    cudaMemset(result + (n - 1), 0, sizeof(int));
39
    for (int stride = n; stride >= 2; stride >>= 1) {
40
        int threads = n / stride, blocks = (threads - 1 + THREADS_PER_BLOCK) / THREADS_PER_BLOCK;
41
        Kernel_down_sweep<<<blocks, THREADS_PER_BLOCK>>>(result, threads, stride);
42
        cudaDeviceSynchronize();
43
    }
44
}

1
Scan Score Table:
2
-------------------------
3
-------------------------------------------------------------------------
4
| Element Count   | Ref Time        | Student Time    | Score           |
5
-------------------------------------------------------------------------
6
| 1000000         | 0.968           | 0.999           | 1.25            |
7
| 10000000        | 10.075          | 8.529           | 1.25            |
8
| 20000000        | 20.351          | 13.911          | 1.25            |
9
| 40000000        | 38.031          | 26.655          | 1.25            |
10
-------------------------------------------------------------------------
11
|                                   | Total score:    | 5.0/5.0         |
12
-------------------------------------------------------------------------

1
__global__ void Kernel_if_repeat(int *arrary, int n, int *output) {
2
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
3
    if (idx + 1 < n) {
4
        output[idx] = (arrary[idx] == arrary[idx + 1]);
5
    }
6
}
7

8
__global__ void Kernel_find_repeats(int *flag, int *index, int n, int *output) {
9
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
10
    if (idx + 1 < n && flag[idx]) {
11
        output[index[idx]] = idx;
12
    }
13
}
14

15
int find_repeats(int* device_input, int length, int* device_output) {
16

17
    // CS149 TODO:
18
    //
19
    // Implement this function. You will probably want to
20
    // make use of one or more calls to exclusive_scan(), as well as
21
    // additional CUDA kernel launches.
22
    //
23
    // Note: As in the scan code, the calling code ensures that
24
    // allocated arrays are a power of 2 in size, so you can use your
25
    // exclusive_scan function with them. However, your implementation
26
    // must ensure that the results of find_repeats are correct given
27
    // the actual array length.
28
    int n = nextPow2(length);
29
    int *flag = nullptr, *sum = nullptr;
30
    cudaMalloc(&flag, sizeof(int) * n);
31
    cudaMalloc(&sum, sizeof(int) * n);
32
    int blocks = (n + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK;
33
    Kernel_if_repeat <<<blocks, THREADS_PER_BLOCK>>> (device_input, n, flag);
34
    cudaDeviceSynchronize();
35
    cudaMemcpy(sum, flag, sizeof(int) * n, cudaMemcpyDeviceToDevice);
36
    exclusive_scan(sum, n, sum);
37
    Kernel_find_repeats <<<blocks, THREADS_PER_BLOCK>>> (flag, sum, n, device_output);
38
    cudaDeviceSynchronize();
39
    int cnt;
40
    cudaMemcpy(&cnt, sum + length - 1, sizeof(int), cudaMemcpyDeviceToHost);
41
    cudaFree(flag); cudaFree(sum);
42
    return cnt;
43
}

1
Find_repeats Score Table:
2
-------------------------
3
-------------------------------------------------------------------------
4
| Element Count   | Ref Time        | Student Time    | Score           |
5
-------------------------------------------------------------------------
6
| 1000000         | 1.887           | 2.02            | 1.25            |
7
| 10000000        | 16.067          | 14.997          | 1.25            |
8
| 20000000        | 28.853          | 26.017          | 1.25            |
9
| 40000000        | 56.414          | 50.418          | 1.25            |
10
-------------------------------------------------------------------------
11
|                                   | Total score:    | 5.0/5.0         |
12
-------------------------------------------------------------------------

1
sudo apt update
2
sudo apt install freeglut3-dev

1
对于每一帧：
2
  清空图像
3
  对于每一个圆：
4
    更新其位置和速度
5
  对于每一个圆：
6
    计算出能包含该圆的正方形
7
    对于这个正方形内的每一个坐标：
8
      计算这个坐标的中心
9
      如果这个中心在圆内：
10
        计算这个圆在此处的颜色
11
        使用这个颜色的贡献更新坐标像素上

1
   result_r = C_alpha * C_r + (1.0 - C_alpha) * P_r
2
   result_g = C_alpha * C_g + (1.0 - C_alpha) * P_g
3
   result_b = C_alpha * C_b + (1.0 - C_alpha) * P_b

1
__global__ void kernelRenderCircles() {
2

3
    int index = blockIdx.x * blockDim.x + threadIdx.x;
4

5
    if (index >= cuConstRendererParams.numCircles)
6
        return;
7

8
    int index3 = 3 * index;
9

10
    // read position and radius
11
    float3 p = *(float3*)(&cuConstRendererParams.position[index3]);
12
    float  rad = cuConstRendererParams.radius[index];
13

14
    // compute the bounding box of the circle. The bound is in integer
15
    // screen coordinates, so it's clamped to the edges of the screen.
16
    short imageWidth = cuConstRendererParams.imageWidth;
17
    short imageHeight = cuConstRendererParams.imageHeight;
18
    short minX = static_cast<short>(imageWidth * (p.x - rad));
19
    short maxX = static_cast<short>(imageWidth * (p.x + rad)) + 1;
20
    short minY = static_cast<short>(imageHeight * (p.y - rad));
21
    short maxY = static_cast<short>(imageHeight * (p.y + rad)) + 1;
22

23
    // a bunch of clamps.  Is there a CUDA built-in for this?
24
    short screenMinX = (minX > 0) ? ((minX < imageWidth) ? minX : imageWidth) : 0;
25
    short screenMaxX = (maxX > 0) ? ((maxX < imageWidth) ? maxX : imageWidth) : 0;
26
    short screenMinY = (minY > 0) ? ((minY < imageHeight) ? minY : imageHeight) : 0;
27
    short screenMaxY = (maxY > 0) ? ((maxY < imageHeight) ? maxY : imageHeight) : 0;
28

29
    float invWidth = 1.f / imageWidth;
30
    float invHeight = 1.f / imageHeight;
31

32
    // for all pixels in the bonding box
33
    for (int pixelY=screenMinY; pixelY<screenMaxY; pixelY++) {
34
        float4* imgPtr = (float4*)(&cuConstRendererParams.imageData[4 * (pixelY * imageWidth + screenMinX)]);
35
        for (int pixelX=screenMinX; pixelX<screenMaxX; pixelX++) {
36
            float2 pixelCenterNorm = make_float2(invWidth * (static_cast<float>(pixelX) + 0.5f),
37
                                                 invHeight * (static_cast<float>(pixelY) + 0.5f));
38
            shadePixel(index, pixelCenterNorm, p, imgPtr);
39
            imgPtr++;
40
        }
41
    }
42
}

1
__global__ void kernelRenderPixels() {
2
    int x = blockIdx.x * blockDim.x + threadIdx.x;
3
    int y = blockIdx.y * blockDim.y + threadIdx.y;
4
    int width = cuConstRendererParams.imageWidth, height = cuConstRendererParams.imageHeight;
5
    if (x >= width || y >= height) return;
6
    float invwidth = 1.f / width, invheight = 1.f / height;
7
    float2 pixelCenterNorm = make_float2(invwidth * (static_cast<float>(x) + 0.5f), invheight * (static_cast<float>(y) + 0.5f));
8
    float4 *imgPtr = (float4*)(&cuConstRendererParams.imageData[4 * (y * width + x)]);
9
    for (int i = 0; i < cuConstRendererParams.numCircles; i++) {
10
        float3 p = *(float3*)(&cuConstRendererParams.position[3 * i]);
11
        shadePixel(i, pixelCenterNorm, p, imgPtr);
12
    }
13
}

1
Score table:
2
------------
3
--------------------------------------------------------------------------
4
| Scene Name      | Ref Time (T_ref) | Your Time (T)   | Score           |
5
--------------------------------------------------------------------------
6
| rgb             | 0.3232           | 0.2669          | 9               |
7
| rand10k         | 2.8221           | 39.2085         | 2               |
8
| rand100k        | 23.3191          | 382.3064        | 2               |
9
| pattern         | 0.6382           | 5.5301          | 3               |
10
| snowsingle      | 15.0153          | 343.2084        | 2               |
11
| biglittle       | 13.4723          | 43.2001         | 5               |
12
| rand1M          | 154.3743         | 3913.552        | 2               |
13
| micro2M         | 294.7338         | 7892.6395       | 2               |
14
--------------------------------------------------------------------------
15
|                                    | Total score:    | 27/72           |
16
--------------------------------------------------------------------------

1
#define BLOCKX 16
2
#define BLOCKY 16
3
#define THREADX 16
4
#define THREADY 16
5

6
__global__ void kernelRenderPixels() {
7
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
8
    int idy = blockIdx.y * blockDim.y + threadIdx.y;
9
    int width = cuConstRendererParams.imageWidth, height = cuConstRendererParams.imageHeight;
10
    float invwidth = 1.f / width, invheight = 1.f / height;
11
    int X0 = idx * THREADX, X1 = X0 + THREADX;
12
    int Y0 = idy * THREADY, Y1 = Y0 + THREADY;
13
    if (X0 >= width || Y0 >= height) return;
14
    float normX0 = X0 * invwidth, normX1 = X1 * invwidth;
15
    float normY0 = Y0 * invheight, normY1 = Y1 * invheight;
16
    int num = cuConstRendererParams.numCircles;
17
    for (int i = 0; i < num; i++) {
18
        float cirx = cuConstRendererParams.position[3 * i], ciry = cuConstRendererParams.position[3 * i + 1];
19
        float cirr = cuConstRendererParams.radius[i];
20
        if (circleInBox(cirx, ciry, cirr, normX0, normX1, normY1, normY0)) {
21
            float3 p = *(float3*)(&cuConstRendererParams.position[3 * i]);
22
            for (int j = Y0; j < Y1; j++) {
23
                float _y = static_cast <float> (j) + 0.5f;
24
                float4* imgPtr = (float4*)(&cuConstRendererParams.imageData[4 * (j * width + X0)]);
25
                for (int k = X0; k < X1; k++) {
26
                    float _x = static_cast <float> (k) + 0.5f;
27
                    shadePixel(i, make_float2(invwidth * _x, invheight * _y), p, imgPtr);
28
                    ++imgPtr;
29
                }
30
            }
31
        }
32
    }
33
}
34

35
void
36
CudaRenderer::render() {
37

38
    // 256 threads per block is a healthy number
39
    dim3 blockDim(THREADX, THREADY);
40
    dim3 gridDim((image->width + blockDim.x - 1) / blockDim.x, (image->height + blockDim.y - 1) / blockDim.y);
41

42
    kernelRenderPixels<<<gridDim, blockDim>>>();
43
    cudaDeviceSynchronize();
44
}

1
Score table:
2
------------
3
--------------------------------------------------------------------------
4
| Scene Name      | Ref Time (T_ref) | Your Time (T)   | Score           |
5
--------------------------------------------------------------------------
6
| rgb             | 0.2807           | 0.8516          | 5               |
7
| rand10k         | 2.729            | 29.2422         | 2               |
8
| rand100k        | 23.1492          | 234.5919        | 2               |
9
| pattern         | 0.6532           | 1.1787          | 6               |
10
| snowsingle      | 14.8701          | 30.8657         | 6               |
11
| biglittle       | 13.9142          | 416.6634        | 2               |
12
| rand1M          | 154.5485         | 523.2159        | 5               |
13
| micro2M         | 296.6468         | 496.183         | 7               |
14
--------------------------------------------------------------------------
15
|                                    | Total score:    | 35/72           |
16
--------------------------------------------------------------------------

1
#define BLOCKX 16
2
#define BLOCKY 16
3
#define BLOCKSIZE 256
4
#define SCAN_BLOCK_DIM   BLOCKSIZE
5
#include "exclusiveScan.cu_inl"
6

7
__global__ void kernelRenderPixels() {
8
    __shared__ uint isinside[BLOCKSIZE], insidesum[BLOCKSIZE], tmp[BLOCKSIZE << 1], cir[BLOCKSIZE];
9
    int width = cuConstRendererParams.imageWidth, height = cuConstRendererParams.imageHeight;
10
    float invwidth = 1.0f / width, invheight = 1.0f / height;
11
    int X0 = blockIdx.x * BLOCKX, X1 = min(X0 + BLOCKX, width);
12
    int Y0 = blockIdx.y * BLOCKY, Y1 = min(Y0 + BLOCKY, height);
13
    float normX0 = X0 * invwidth, normX1 = X1 * invwidth;
14
    float normY0 = Y0 * invheight, normY1 = Y1 * invheight;
15
    int threadid = threadIdx.y * BLOCKX + threadIdx.x;
16
    int num = cuConstRendererParams.numCircles;
17
    int pixelx = X0 + threadIdx.x;
18
    int pixely = Y0 + threadIdx.y;
19
    float2 pixel = make_float2(invwidth * (static_cast<float>(pixelx) + 0.5f), invheight * (static_cast <float>(pixely) + 0.5f));
20
    float4 *imgPtr = (float4*)(&cuConstRendererParams.imageData[4 * (pixely * width + pixelx)]);
21
    for (int _i = 0; _i < num; _i += BLOCKSIZE) {
22
        int i = _i + threadid;
23
        if (i < num) {
24
            float cirx = cuConstRendererParams.position[3 * i], ciry = cuConstRendererParams.position[3 * i + 1];
25
            float cirr = cuConstRendererParams.radius[i];
26
            isinside[threadid] = circleInBox(cirx, ciry, cirr, normX0, normX1, normY1, normY0);
27
        }
28
        else isinside[threadid] = 0;
29

30
        __syncthreads();
31
        sharedMemExclusiveScan(threadid, isinside, insidesum, tmp, BLOCKSIZE);
32
        if (isinside[threadid]) {
33
            cir[insidesum[threadid]] = i;
34
        }
35
        __syncthreads();
36
        if (pixelx < width && pixely < height) {
37
            int cirnum = insidesum[BLOCKSIZE - 1] + isinside[BLOCKSIZE - 1];
38
            for (int j = 0; j < cirnum; j++) {
39
                int ciridx = cir[j];
40
                float3 p = *(float3*)(&cuConstRendererParams.position[ciridx * 3]);
41
                shadePixel(ciridx, pixel, p, imgPtr);
42
            }
43
        }
44
    }
45
}
46

47
void
48
CudaRenderer::render() {
49

50
    // 256 threads per block is a healthy number
51
    dim3 blockDim(BLOCKX, BLOCKY);
52
    dim3 gridDim((image->width + blockDim.x - 1) / blockDim.x, (image->height + blockDim.y - 1) / blockDim.y);
53

54
    kernelRenderPixels<<<gridDim, blockDim>>>();
55
    cudaDeviceSynchronize();
56
}

1
Score table:
2
------------
3
--------------------------------------------------------------------------
4
| Scene Name      | Ref Time (T_ref) | Your Time (T)   | Score           |
5
--------------------------------------------------------------------------
6
| rgb             | 0.4447           | 0.4018          | 9               |
7
| rand10k         | 3.3763           | 3.7078          | 9               |
8
| rand100k        | 28.7109          | 32.373          | 9               |
9
| pattern         | 0.5684           | 0.5763          | 9               |
10
| snowsingle      | 20.157           | 18.984          | 9               |
11
| biglittle       | 13.8375          | 32.426          | 5               |
12
| rand1M          | 153.9384         | 153.1018        | 9               |
13
| micro2M         | 295.5285         | 288.3962        | 9               |
14
--------------------------------------------------------------------------
15
|                                    | Total score:    | 68/72           |