C++ Request of assessing a large scale matriks multiplication code

[wanderef]

Long story short I'm an amateur undergraduate progranmer who was tasked to create a program regrarding matriks.
While I succeded in creating the program. I have no means to test it's peformance or to assess which part is corect or wrong as this is only an extra project of additional point

So I humbly request for anyone to rate this program and tell me what i can improve on

Note: the scale of the program is corect it is a large scale program that require High-end GPU or supercomputer to use

#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <stdio.h>
#include <stdlib.h>
#include <chrono>
#include <fstream>

#define N 1048576 // Ukuran matriks besar
#define TILE_SIZE 256 // Ukuran tile optimal
#define STREAM_COUNT 16 // Jumlah stream untuk overlapping

void checkCublas(cublasStatus_t status) {
if (status != CUBLAS_STATUS_SUCCESS) {
printf("cuBLAS Error\n");
exit(EXIT_FAILURE);
}
}

global void fusedGemmKernel(half *A, half *B, float *C, int n) {
shared half tile_A[TILE_SIZE][TILE_SIZE];
shared half tile_B[TILE_SIZE][TILE_SIZE];

int tx = threadIdx.x, ty = threadIdx.y;
int row = blockIdx.y * TILE_SIZE + ty;
int col = blockIdx.x * TILE_SIZE + tx;
float value = 0.0f;

for (int k = 0; k < n / TILE_SIZE; ++k) {
tile_A[ty][tx] = A[row * n + (k * TILE_SIZE + tx)];
tile_B[ty][tx] = B[(k * TILE_SIZE + ty) * n + col];
__syncthreads();

for (int i = 0; i < TILE_SIZE; i++) {
value += __half2float(tile_A[ty]) * __half2float(tile_B[tx]);
}
__syncthreads();
}
C[row * n + col] = value;
}

int main() {
half *A, *B;
float *C;
cudaMallocHost(&A, N * N * sizeof(half));
cudaMallocHost(&B, N * N * sizeof(half));
cudaMallocHost(&C, N * N * sizeof(float));

half *d_A, *d_B;
float *d_C;
cudaMalloc(&d_A, N * N * sizeof(half));
cudaMalloc(&d_B, N * N * sizeof(half));
cudaMalloc(&d_C, N * N * sizeof(float));

cudaStream_t streams[STREAM_COUNT];
for (int i = 0; i < STREAM_COUNT; i++) {
cudaStreamCreate(&streams);
}

dim3 gridDim(N / TILE_SIZE, N / TILE_SIZE);
dim3 blockDim(TILE_SIZE, TILE_SIZE);

cudaMemcpyAsync(d_A, A, N * N * sizeof(half), cudaMemcpyHostToDevice, streams[0]);
cudaMemcpyAsync(d_B, B, N * N * sizeof(half), cudaMemcpyHostToDevice, streams[1]);

auto start = std::chrono::high_resolution_clock::now();
fusedGemmKernel<<<gridDim, blockDim, 0, streams[2]>>>(d_A, d_B, d_C, N);
cudaDeviceSynchronize();
auto end = std::chrono::high_resolution_clock::now();
std::chrono::duration<float> duration = end - start;

cudaMemcpyAsync(C, d_C, N * N * sizeof(float), cudaMemcpyDeviceToHost, streams[3]);
cudaDeviceSynchronize();

printf("Execution Time: %.5f s\n", duration.count());

// Menghitung TFLOPS
double num_operations = 2.0 * N * N * N;
double tflops = (num_operations / (duration.count() * 1.0e12));
printf("Performance: %.5f TFLOPS\n", tflops);

// Menyimpan hasil ke dalam file
std:😱fstream outfile("matrix_output.txt");
for (int i = 0; i < 10; i++) { // Simpan hanya sebagian kecil untuk efisiensi
for (int j = 0; j < 10; j++) {
outfile << C[i * N + j] << " ";
}
outfile << "\n";
}
outfile.close();
printf("Matrix output saved to matrix_output.txt\n");

for (int i = 0; i < STREAM_COUNT; i++) {
cudaStreamDestroy(streams);
}
cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
cudaFreeHost(A); cudaFreeHost(B); cudaFreeHost(C);
return 0;
}

This is the code in question