I have a matrix class with a static method to multiply two matrices with 10k rows and 10k columns dimensions. I have used a library called ILGPU to try and execute the operation through my GPU and to be honest, for a trillion operations algorithm, the GPU took 20 seconds to complete the process(credits to the thousands of cuda cores on it for parallelism), the CPU on the other hand takes forever since it has not completed since it started, I have tried to optimize the code to run faster on cpu with no sucess, can somebody help me bring the performance of the cpu based algorithm on par with the gpu based algorithm.

C#
//this method uses the gpu to do the multiplication with nuget package(ILGPU)
 public static Matrix MultiplyGpu(Matrix A, Matrix B)
 {
     if (A.Cols != B.Rows)
         throw new Exception("Matrices cannot be multiplied, they are not congruent");

     int n = A.Rows;                    // assuming square matrices for now
     var result = new Matrix(n, B.Cols);

     var sw = Stopwatch.StartNew();

     using var context = Context.CreateDefault();
     using var accelerator = context.GetPreferredDevice(preferCPU: false).CreateAccelerator(context);

     using var bufferA = accelerator.Allocate2DDenseX(new Index2D(n, n));
     using var bufferB = accelerator.Allocate2DDenseX(new Index2D(n, n));
     using var bufferC = accelerator.Allocate2DDenseX(new Index2D(n, n));

     // Convert double -> float and copy to GPU
     bufferA.CopyFromCPU(A.ToFloatArray());
     bufferB.CopyFromCPU(B.ToFloatArray());

     // Load kernel
     var kernel = accelerator.LoadAutoGroupedStreamKernel<
         Index2D,
         ArrayView2D,
         ArrayView2D,
         ArrayView2D,
         int>(MatrixMulKernel);

     kernel(new Index2D(n, n), bufferA.View, bufferB.View, bufferC.View, n);
     accelerator.Synchronize();

     // Copy result back
     var resultFloat = new float[n, n];
     bufferC.CopyToCPU(resultFloat);
     result.LoadFromFloatArray(resultFloat);

     Console.WriteLine($"GPU Multiply on RTX 4060 completed in: {sw.Elapsed}");

     return result;
 }

 static void MatrixMulKernel(
     Index2D index,
     ArrayView2D a,
     ArrayView2D b,
     ArrayView2D c,
     int n)
 {
     int row = index.X;
     int col = index.Y;
     float sum = 0f;
     for (int k = 0; k < n; k++)
         sum += a[row, k] * b[k, col];
     c[row, col] = sum;
 }


//THIS IS THE CPU BASED METHOD FOR CARRYING OUT THE MULTIPLICATION
 public static Matrix MultiplyParallel(Matrix A, Matrix B)
 {
     if (A.Cols != B.Rows)
     {
         throw new Exception("Matrices cannot be multiplied, they are not congruent");
     }

     int n = A.Rows;
     int m = B.Cols;
     var result = new Matrix(n, m);

     const int BLOCK_SIZE = 64;   // ← Tune this (try 32, 64, 128)

     var sw = Stopwatch.StartNew();

     Parallel.For(0, n, i =>
     {
         for (int kk = 0; kk < A.Cols; kk += BLOCK_SIZE)
         {
             for (int jj = 0; jj < m; jj += BLOCK_SIZE)
             {
                 for (int k = kk; k < Math.Min(kk + BLOCK_SIZE, A.Cols); k++)
                 {
                     double a_ik = A[i, k];        // hoist value

                     for (int j = jj; j < Math.Min(jj + BLOCK_SIZE, m); j++)
                     {
                         result[i, j] += a_ik * B[k, j];
                     }
                 }
             }
         }
     });

     Console.WriteLine($"Optimized MultiplyParallel took: {sw.Elapsed}");
     return result;
 }

I read from two large csv containing numbers and load them into the matrices, please help.