Shared Code

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
#include <vector>
#include <cmath>

#define BORDE_SUPERIOR 100.0f
#define BORDES 0.0f

#define BLOCK_W 16
#define BLOCK_H 16
#define INNER_W 14
#define INNER_H 14

int MatrixSide = 1024;
int N_Steps = 5;

// ------------------- Kernel original -------------------
__global__ void UpdateMatrixOriginal(const int n, float* input, float* buffer)
{
    int tx = threadIdx.x;
    int ty = threadIdx.y;
    int idx_col = blockIdx.x * blockDim.x + tx;
    int idx_row = blockIdx.y * blockDim.y + ty;

    __shared__ float shMemHalo[BLOCK_H + 2][BLOCK_W + 2];

    int sharedX = tx + 1;
    int sharedY = ty + 1;

    float value = 0.0f;
    if (idx_col < n && idx_row < n)
        value = input[idx_row * n + idx_col];
    shMemHalo[sharedY][sharedX] = value;

    // Halos
    if (tx == 0)
        shMemHalo[sharedY][0] = (idx_col > 0) ? input[idx_row * n + idx_col - 1] : BORDES;
    if (tx == blockDim.x - 1 || idx_col == n - 1)
        shMemHalo[sharedY][sharedX + 1] = (idx_col < n - 1) ? input[idx_row * n + idx_col + 1] : BORDES;
    if (ty == 0)
        shMemHalo[0][sharedX] = (idx_row > 0) ? input[(idx_row - 1) * n + idx_col] : BORDE_SUPERIOR;
    if (ty == blockDim.y - 1 || idx_row == n - 1)
        shMemHalo[sharedY + 1][sharedX] = (idx_row < n - 1) ? input[(idx_row + 1) * n + idx_col] : BORDES;

    __syncthreads();

    if (idx_col > 0 && idx_col < n - 1 && idx_row > 0 && idx_row < n - 1)
    {
        buffer[idx_row * n + idx_col] = 0.25f * (
            shMemHalo[sharedY - 1][sharedX] +
            shMemHalo[sharedY + 1][sharedX] +
            shMemHalo[sharedY][sharedX - 1] +
            shMemHalo[sharedY][sharedX + 1]
        );
    }
}

// ------------------- Kernel optimizado 14x14 (usando 16x16 tile con halo) -------------------
__global__ void UpdateMatrixTile14x14(const int n, float* input, float* buffer)
{
    __shared__ float tile[BLOCK_H + 2][BLOCK_W + 2];

    int tx = threadIdx.x;
    int ty = threadIdx.y;

    int global_x = blockIdx.x * INNER_W + tx - 1;
    int global_y = blockIdx.y * INNER_H + ty - 1;

    // Cargar valor dentro de tile (con halo)
    float value = 0.0f;
    if (global_x >= 0 && global_x < n && global_y >= 0 && global_y < n)
        value = input[global_y * n + global_x];

    tile[ty][tx] = value;

    __syncthreads();

    // Región interna del bloque (14x14)
    if (tx >= 1 && tx <= INNER_W && ty >= 1 && ty <= INNER_H &&
        global_x > 0 && global_x < n - 1 && global_y > 0 && global_y < n - 1)
    {
        buffer[global_y * n + global_x] = 0.25f * (
            tile[ty - 1][tx] +
            tile[ty + 1][tx] +
            tile[ty][tx - 1] +
            tile[ty][tx + 1]
        );
    }
}

// ------------------- Checksum -------------------
float checksum(const std::vector<float>& mat)
{
    double sum = 0.0;
    for (float v : mat)
        sum += v;
    return static_cast<float>(sum);
}

// ------------------- Main -------------------
int main()
{
    size_t bytes = MatrixSide * MatrixSide * sizeof(float);
    std::vector<float> h_matrix(MatrixSide * MatrixSide, 0.0f);

    // Borde superior
    for (int col = 0; col < MatrixSide; col++)
        h_matrix[col] = BORDE_SUPERIOR;

    float* d_in;
    float* d_out;
    cudaMalloc(&d_in, bytes);
    cudaMalloc(&d_out, bytes);
    cudaMemcpy(d_in, h_matrix.data(), bytes, cudaMemcpyHostToDevice);

    dim3 dimBlock(BLOCK_W, BLOCK_H);
    dim3 dimGrid((MatrixSide + INNER_W - 1) / INNER_W, (MatrixSide + INNER_H - 1) / INNER_H);

    // Warmup
    UpdateMatrixOriginal<<<dimGrid, dimBlock>>>(MatrixSide, d_in, d_out);
    cudaDeviceSynchronize();
    UpdateMatrixTile14x14<<<dimGrid, dimBlock>>>(MatrixSide, d_in, d_out);
    cudaDeviceSynchronize();

    // --- Kernel original ---
    cudaMemcpy(d_in, h_matrix.data(), bytes, cudaMemcpyHostToDevice);
    cudaEvent_t start1, stop1;
    cudaEventCreate(&start1);
    cudaEventCreate(&stop1);
    cudaEventRecord(start1);
    for (int step = 0; step < N_Steps; step++)
    {
        UpdateMatrixOriginal<<<dimGrid, dimBlock>>>(MatrixSide, d_in, d_out);
        cudaDeviceSynchronize();
        std::swap(d_in, d_out);
    }
    cudaEventRecord(stop1);
    cudaEventSynchronize(stop1);
    float timeOriginal = 0;
    cudaEventElapsedTime(&timeOriginal, start1, stop1);
    std::cout << "Tiempo kernel original: " << timeOriginal << " ms\n";

    std::vector<float> h_result_original(MatrixSide * MatrixSide);
    cudaMemcpy(h_result_original.data(), d_in, bytes, cudaMemcpyDeviceToHost);
    float checksum_original = checksum(h_result_original);
    std::cout << "Checksum kernel original: " << checksum_original << "\n";

    // --- Kernel optimizado ---
    cudaMemcpy(d_in, h_matrix.data(), bytes, cudaMemcpyHostToDevice);
    cudaEvent_t start2, stop2;
    cudaEventCreate(&start2);
    cudaEventCreate(&stop2);
    cudaEventRecord(start2);
    for (int step = 0; step < N_Steps; step++)
    {
        UpdateMatrixTile14x14<<<dimGrid, dimBlock>>>(MatrixSide, d_in, d_out);
        cudaDeviceSynchronize();
        std::swap(d_in, d_out);
    }
    cudaEventRecord(stop2);
    cudaEventSynchronize(stop2);
    float timeTile = 0;
    cudaEventElapsedTime(&timeTile, start2, stop2);
    std::cout << "Tiempo kernel optimizado (14x14): " << timeTile << " ms\n";

    std::vector<float> h_result_tile(MatrixSide * MatrixSide);
    cudaMemcpy(h_result_tile.data(), d_in, bytes, cudaMemcpyDeviceToHost);
    float checksum_tile = checksum(h_result_tile);
    std::cout << "Checksum kernel optimizado: " << checksum_tile << "\n";

    // --- Comparación ---
    float epsilon = 1e-4f;
    if (std::fabs(checksum_original - checksum_tile) < epsilon)
        std::cout << "✅ Los resultados coinciden!\n";
    else
        std::cout << "❌ Diferencia en resultados!\n";

    cudaFree(d_in);
    cudaFree(d_out);
    return 0;
}