// Define the tile sizes and thread configurations
#define TILE_SIZE_M 128
#define TILE_SIZE_N 128
#define TILE_SIZE_K 8
#define NUM_THREADS 256  // Number of threads per block
#define WARP_SIZE 32     // Warp size

// Shared memory leading dimensions (adjust as necessary)
#define SMEM_LDA (TILE_SIZE_K + 4)  // Adjusted to ensure proper alignment
#define SMEM_LDB TILE_SIZE_N
#define SMEM_LDC (TILE_SIZE_M / 2)

// Helper functions
__device__ __forceinline__ void ldg32_nc_0(float& reg, const void* ptr) {
    asm volatile("{.reg .pred p;\n"
        "mov.b32 %0, 0;\n"
#if __CUDACC_VER_MAJOR__ >= 11 && __CUDACC_VER_MINOR__ >= 4 &&                 \
    __CUDA_ARCH__ >= 750
        "ld.global.nc.L2::128B.f32 %0, [%1];}\n"
#else
        "ld.global.nc.f32 %0, [%1];}\n"
#endif
        : "=f"(reg)
        : "l"(ptr));
}

__device__ __forceinline__ uint32_t smem_u32addr(const void* smem_ptr) {
    uint32_t addr;
    asm("{.reg .u64 u64addr;\n"
        " cvta.to.shared.u64 u64addr, %1;\n"
        " cvt.u32.u64 %0, u64addr;}\n"
        : "=r"(addr)
        : "l"(smem_ptr));

    return addr;
}

__device__ __forceinline__ void lds128(float& reg0, float& reg1, float& reg2,
    float& reg3, const uint32_t& addr) {
    asm volatile("ld.shared.v4.f32 {%0, %1, %2, %3}, [%4];\n"
        : "=f"(reg0), "=f"(reg1), "=f"(reg2), "=f"(reg3)
        : "r"(addr));
}

__device__ __forceinline__ void sts128(const float& reg0, const float& reg1,
    const float& reg2, const float& reg3,
    const uint32_t& addr) {
    asm volatile("st.shared.v4.f32 [%0], {%1, %2, %3, %4};\n"
        :
    : "r"(addr), "f"(reg0), "f"(reg1), "f"(reg2), "f"(reg3));
}

__device__ __forceinline__ void stg128(const float& reg0, const float& reg1,
    const float& reg2, const float& reg3,
    const float* addr) {
    asm volatile("st.global.v4.f32 [%0], {%1, %2, %3, %4};\n"
        :
    : "l"(addr), "f"(reg0), "f"(reg1), "f"(reg2), "f"(reg3));
}

__device__ __forceinline__ void stg32(const float& reg, const void* ptr) {
    asm volatile("{.reg .pred p;\n"
        " st.global.f32 [%0], %1;}\n"
        :
    : "l"(ptr), "f"(reg));
}

__device__ __forceinline__ void sts32(const float& reg, const uint32_t& addr) {
    asm volatile("st.shared.f32 [%0], %1;\n" : : "r"(addr), "f"(reg));
}

/**
 // Generalized SGEMM kernel with customizable parameters
*/
__global__ __launch_bounds__(NUM_THREADS, 2) void sgemm_kernel(int m, int n, int k,
    const float* a,
    const float* b,
    float* c) {

    // Adjust shared memory size based on tile sizes
    extern __shared__ char smem[];
    float* ashare = reinterpret_cast<float*>(smem);
    float* bshare = reinterpret_cast<float*>(smem + TILE_SIZE_M * SMEM_LDA * sizeof(float));

    float sum[8][8] = { 0 };
    float panelA[2][8] = { 0 }, panelB[2][8] = { 0 };

    int tid = threadIdx.x;

    int block_row = blockIdx.y * TILE_SIZE_M;
    int block_col = blockIdx.x * TILE_SIZE_N;

    // Calculate global memory offsets for loading A and B
    int from_a = (block_row + (tid / 8) * 4) * k + (tid % 8);
    int from_b = (tid / WARP_SIZE) * n + block_col + (tid % WARP_SIZE);

    float a_ldg_reg[4], b_ldg_reg[4];

    uint32_t a_sts_addr = smem_u32addr(ashare + (tid % 8) * SMEM_LDA + (tid / 8) * 4);
    uint32_t b_sts_addr = smem_u32addr(bshare + (tid / WARP_SIZE) * SMEM_LDB + (tid % WARP_SIZE));

    uint32_t aptr_base = smem_u32addr(ashare + (tid / 16) * 4);
    uint32_t bptr_base = smem_u32addr(bshare + (tid % 16) * 4);

    {
        // Load first panels of A and B from global memory to shared memory
#pragma unroll
        for (int i = 0; i < 4; ++i)
        {
            ldg32_nc_0(a_ldg_reg[i], (const char*)(a + from_a) + i * k * sizeof(float));
        }
        sts128(a_ldg_reg[0], a_ldg_reg[1], a_ldg_reg[2], a_ldg_reg[3], a_sts_addr);

#pragma unroll
        for (int i = 0; i < 4; ++i)
        {
            ldg32_nc_0(b_ldg_reg[i], (const char*)(b + from_b) + i * WARP_SIZE * sizeof(float));
        }
#pragma unroll
        for (int i = 0; i < 4; ++i)
        {
            sts32(b_ldg_reg[i], b_sts_addr + i * WARP_SIZE * sizeof(float));
        }
        __syncthreads();
        // Update offsets and flip flags for double-buffering
        from_a += TILE_SIZE_K;
        from_b += TILE_SIZE_K * n;

        a_sts_addr ^= 0x2000;
        b_sts_addr ^= 0x1000;
    }

    // Load first panel
    lds128(panelA[0][0], panelA[0][1], panelA[0][2], panelA[0][3], aptr_base);
    lds128(panelA[0][4], panelA[0][5], panelA[0][6], panelA[0][7], aptr_base + 64 * sizeof(float));

    lds128(panelB[0][0], panelB[0][1], panelB[0][2], panelB[0][3], bptr_base);
    lds128(panelB[0][4], panelB[0][5], panelB[0][6], panelB[0][7], bptr_base + 64 * sizeof(float));

    for (int loop = 0; loop < k; loop += TILE_SIZE_K)
    {
        // Compute
#pragma unroll
        for (int subk = 0; subk < TILE_SIZE_K; ++subk)
        {
            if (subk == TILE_SIZE_K - 1 && loop < k - TILE_SIZE_K)
            {
                // If there is more data, load next panels
                sts128(a_ldg_reg[0], a_ldg_reg[1], a_ldg_reg[2], a_ldg_reg[3], a_sts_addr);

#pragma unroll
                for (int i = 0; i < 4; ++i)
                {
                    sts32(b_ldg_reg[i], b_sts_addr + i * WARP_SIZE * sizeof(float));
                }
                __syncthreads();
                from_a += TILE_SIZE_K;
                from_b += TILE_SIZE_K * n;

                aptr_base ^= 0x2000;
                bptr_base ^= 0x1000;
                a_sts_addr ^= 0x2000;
                b_sts_addr ^= 0x1000;
            }

            const int pp = (subk + 1) % 2; // Ping-pong index
            lds128(panelA[pp][0], panelA[pp][1], panelA[pp][2], panelA[pp][3], aptr_base + ((subk + 1) % TILE_SIZE_K) * SMEM_LDA * sizeof(float));
            lds128(panelA[pp][4], panelA[pp][5], panelA[pp][6], panelA[pp][7], aptr_base + (((subk + 1) % TILE_SIZE_K) * SMEM_LDA + 64) * sizeof(float));

            lds128(panelB[pp][0], panelB[pp][1], panelB[pp][2], panelB[pp][3], bptr_base + ((subk + 1) % TILE_SIZE_K) * SMEM_LDB * sizeof(float));
            lds128(panelB[pp][4], panelB[pp][5], panelB[pp][6], panelB[pp][7], bptr_base + (((subk + 1) % TILE_SIZE_K) * SMEM_LDB + 64) * sizeof(float));

            if (subk == 0 && loop < k - TILE_SIZE_K)
            {
#pragma unroll
                for (int i = 0; i < 4; ++i)
                {
                    ldg32_nc_0(a_ldg_reg[i],
                        (const char*)(a + from_a) + i * k * sizeof(float));
                }
#pragma unroll
                for (int i = 0; i < 4; ++i)
                {
                    ldg32_nc_0(b_ldg_reg[i],
                        (const char*)(b + from_b) + i * WARP_SIZE * sizeof(float));
                }
            }

#pragma unroll
            for (int i = 0; i < 8; ++i)
            {
#pragma unroll
                for (int j = 0; j < 8; ++j)
                {
                    sum[i][j] += panelA[subk % 2][i] * panelB[subk % 2][j];
                }
            }
        }
    }

    int write_offset = (block_row + (tid / 16) * 4) * n + block_col + (tid % 16) * 4;
#pragma unroll
    for (int i = 0; i < 4; ++i)
    {
        stg128(sum[i][0], sum[i][1], sum[i][2], sum[i][3], c + write_offset + i * n);
        stg128(sum[i][4], sum[i][5], sum[i][6], sum[i][7], c + write_offset + i * n + 64);
        stg128(sum[i + 4][0], sum[i + 4][1], sum[i + 4][2], sum[i + 4][3], c + write_offset + (i + 64) * n);
        stg128(sum[i + 4][4], sum[i + 4][5], sum[i + 4][6], sum[i + 4][7], c + write_offset + (i + 64) * n + 64);
    }
}

// Helper function to initialize a matrix with random values
void initializeMatrix(float* matrix, int rows, int cols) {
    for (int i = 0; i < rows * cols; ++i) {
        matrix[i] = static_cast<float>(rand()) / RAND_MAX;
    }
}

// Helper function to verify the result
void verifyResult(const float* a, const float* b, const float* c, int m, int n, int k) {
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
            float sum = 0.0f;
            for (int l = 0; l < k; ++l) {
                sum += a[i * k + l] * b[l * n + j];
            }
            float diff = std::abs(c[i * n + j] - sum);
            if (diff > 1e-4) {
                std::cout << "Mismatch at (" << i << ", " << j << "): "
                    << c[i * n + j] << " vs " << sum << std::endl;
                return;
            }
        }
    }
    std::cout << "Result verified successfully!" << std::endl;
}

void MY_MMult_2(int m, int n, int k, float* d_A,
    float* d_B, float* d_C) {

    dim3 grid((n + TILE_SIZE_N - 1) / TILE_SIZE_N, (m + TILE_SIZE_M - 1) / TILE_SIZE_M);
    dim3 block(NUM_THREADS);

    // Adjust shared memory size based on tile sizes
    size_t sharedMemSize = (TILE_SIZE_M * SMEM_LDA + TILE_SIZE_K * SMEM_LDB) * sizeof(float);

    sgemm_kernel<<<grid, block, sharedMemSize>>>(m, n, k, d_A, d_B, d_C); hmm, what about the tiling and stuff like panel size and sum size, can you please examine the code carefully and see if anything else hardcoded can be extracted. And also write a simple test to test the correctness of the code and also a method to print the data output and input matrices in a beautiful way, keep the data size to a minimum so that we can print them in a pretty way.
}