PMPP Chapter 3 Solutions

#include <cassert>
#include <cuda_runtime.h>

#include <cstdint>

namespace ch03::ex01 {

__global__ void square_matrix_multiply_by_row(const float *left,
                                              const float *right, float *output,
                                              std::uint32_t width) {
    // get output row index
    unsigned row_idx = blockDim.y * blockIdx.y + threadIdx.y;

    if (row_idx >= width)
        return;

    // for each element in the output row, calculate it

    for (unsigned dot_idx = 0; dot_idx < width; dot_idx++) {
        for (unsigned output_col = 0; output_col < width; output_col++) {

            unsigned left_idx = row_idx * width + dot_idx;
            unsigned right_idx = dot_idx * width + output_col;
            unsigned output_idx = row_idx * width + output_col;

            output[output_idx] += left[left_idx] * right[right_idx];
        }
    }
}

void launch_matrix_mul_by_row(const float *left, const float *right,
                              float *output, std::uint32_t width,
                              cudaStream_t stream) {
    constexpr unsigned THREADS_PER_BLOCK = 32;
    const unsigned num_blocks = (width - 1) / THREADS_PER_BLOCK + 1;

    auto grid_dim = dim3(1, num_blocks, 1);
    auto block_dim = dim3(1, THREADS_PER_BLOCK, 1);

    float *d_left, *d_right, *d_output;
    int n_bytes = width * width * sizeof(float);

    cudaMalloc(&d_left, n_bytes);
    cudaMalloc(&d_right, n_bytes);
    cudaMalloc(&d_output, n_bytes);

    cudaMemcpy(d_left, left, n_bytes, cudaMemcpyHostToDevice);
    cudaMemcpy(d_right, right, n_bytes, cudaMemcpyHostToDevice);
    cudaMemset(d_output, 0, n_bytes);

    square_matrix_multiply_by_row<<<grid_dim, block_dim, 0, stream>>>(
        d_left, d_right, d_output, width);

    cudaMemcpy(output, d_output, n_bytes, cudaMemcpyDeviceToHost);
    cudaFree(d_left);
    cudaFree(d_right);
    cudaFree(d_output);
}

__global__ void square_matrix_multiply_by_col(const float *left,
                                              const float *right, float *output,
                                              std::uint32_t width) {
    // get output row index
    unsigned col_idx = blockDim.x * blockIdx.x + threadIdx.x;

    if (col_idx >= width)
        return;

    // for each element in the output row, calculate it
    for (unsigned dot_idx = 0; dot_idx < width; dot_idx++) {
        for (unsigned output_row = 0; output_row < width; output_row++) {

            unsigned left_idx = output_row * width + dot_idx;
            unsigned right_idx = dot_idx * width + col_idx;
            unsigned output_idx = output_row * width + col_idx;

            output[output_idx] += left[left_idx] * right[right_idx];
        }
    }
}

void launch_matrix_mul_by_column(const float *left, const float *right,
                                 float *output, std::uint32_t width,
                                 cudaStream_t stream) {
    constexpr unsigned THREADS_PER_BLOCK = 32;
    const unsigned num_blocks = (width - 1) / THREADS_PER_BLOCK + 1;

    auto grid_dim = dim3(num_blocks, 1, 1);
    auto block_dim = dim3(THREADS_PER_BLOCK, 1, 1);

    float *d_left, *d_right, *d_output;
    int n_bytes = width * width * sizeof(float);

    cudaMalloc(&d_left, n_bytes);
    cudaMalloc(&d_right, n_bytes);
    cudaMalloc(&d_output, n_bytes);

    cudaMemcpy(d_left, left, n_bytes, cudaMemcpyHostToDevice);
    cudaMemcpy(d_right, right, n_bytes, cudaMemcpyHostToDevice);
    cudaMemset(d_output, 0, n_bytes);

    square_matrix_multiply_by_col<<<grid_dim, block_dim, 0, stream>>>(
        d_left, d_right, d_output, width);

    cudaMemcpy(output, d_output, n_bytes, cudaMemcpyDeviceToHost);
    cudaFree(d_left);
    cudaFree(d_right);
    cudaFree(d_output);
}
} // namespace ch03::ex01

#include <cuda_runtime.h>

#include <cstdint>

__global__ void matrix_vec_mult_kernel(const float *matrix, const float *vector,
                                       float *output, std::uint32_t width) {
    // you'd want to load the matrix row major first

    unsigned output_idx = blockDim.y * blockIdx.y + threadIdx.y;

    if (output_idx >= width)
        return;

    for (unsigned dot_idx = 0; dot_idx < width; dot_idx++) {
        output[output_idx] +=
            matrix[output_idx * width + dot_idx] * vector[dot_idx];
    }
}

void matrix_vector_multiply(float *output, const float *matrix,
                            const float *input, std::uint32_t width) {
    float *d_output, *d_matrix, *d_vector;

    unsigned n_matrix_bytes = width * width * sizeof(float);
    unsigned n_vector_bytes = width * sizeof(float);

    cudaMalloc(&d_matrix, n_matrix_bytes);
    cudaMalloc(&d_vector, n_vector_bytes);
    cudaMalloc(&d_output, n_vector_bytes);

    cudaMemcpy(d_matrix, matrix, n_matrix_bytes, cudaMemcpyHostToDevice);
    cudaMemcpy(d_vector, input, n_vector_bytes, cudaMemcpyHostToDevice);
    cudaMemset(d_output, 0, n_vector_bytes);

    matrix_vec_mult_kernel<<<dim3(1, unsigned((width - 1) / 32 + 1), 1),
                             dim3(1, 32, 1)>>>(d_matrix, d_vector, d_output,
                                               width);

    cudaMemcpy(output, d_output, n_vector_bytes, cudaMemcpyDeviceToHost);
    cudaFree(d_matrix);
    cudaFree(d_vector);
    cudaFree(d_output);
}