Is my implementation of Strassen Slower than naive Matrix Multiplication?

Question

So I've implemented Strassen'S Algorithm in C and even for large matrices it's not able to beat the classic ikj algorithm. I only use Matrices that er powers of 2, filled with doubles. I read that for "large" n, Strassen's Algorithm beats the classic one, e.g. for n >= 1000. For me this is not the case. I don't use a cut-off point since, it's not part of the original algorithm. If i do, it's actually faster than the classic one, but without it's performing severely worse. Is this the expected behaviour and the sources I found online about Strassen's Algorithm are actually using a cut off point?

I tried to optimize Strassen, by using only two matrices stored as one-dimensional array, passing them as arguments and actually not creating submatrices, if possible, but just passing the new size and a pointer two the cell so i don't allocate as much memory. But maybe for large Matrices this is actually bad since I'm accessing the whole matric and need to store it somewhere in main memory.

Here is a Github Link: https://github.com/AlexLoitzl/matmul The Code is not properly documented yet... And here is the Strassen Implementation:

void strassen_mult_base(int n, int x_row_length, const double X[], int y_row_length, const double Y[], int z_row_length, double Z[], int base, void (*inner_func)()) {

  if (n == 1){
    Z[0] = X[0] * Y[0];
    return;
  }

  const int k = n/2;
  //Split Matrices X and Y into 4 blocks without any memory allocation by using proper offset
  const double *A = X;
  const double *B = X + k;
  const double *C = X + k*x_row_length;
  const double *D = C + k;

  const double *E = Y;
  const double *F = Y + k;
  const double *G = Y + k*y_row_length;
  const double *H = G + k;

  //Allocate memory for temporary matrices P0 - P6 for our 7 Multiplications
  const int size = k*k*sizeof(double);
  double *P[7];
  for (int i = 0; i < 7; i++) {
    P[i] = (double *) malloc(size);
  }
  //Allocate memory for interim results of calculating P0 - P6
  double *S = (double *) malloc(size);
  double *R = (double *) malloc(size);

  // P0 = A*(F - H);
  sub(k, y_row_length, F, y_row_length, H, k, S);
  strassen_mult_base(k, x_row_length, A, k, S, k, P[0], base, inner_func);

  // P1 = (A + B)*H
  add(k, x_row_length, A, x_row_length, B, k, S);
  strassen_mult_base(k, k, S, y_row_length, H, k, P[1], base, inner_func);

  // P2 = (C + D)*E
  add(k, x_row_length, C, x_row_length, D, k, S);
  strassen_mult_base(k, k, S, y_row_length, E, k, P[2], base, inner_func);

  // P3 = D*(G - E);
  sub(k, y_row_length, G, y_row_length, E, k, S);
  strassen_mult_base(k, x_row_length, D, k, S, k, P[3], base, inner_func);

  // P4 = (A + D)*(E + H)
  add(k, x_row_length, A, x_row_length, D, k, S);
  add(k, y_row_length, E, y_row_length, H, k, R);
  strassen_mult_base(k, k, S, k, R, k, P[4], base, inner_func);

  // P5 = (B - D)*(G + H)
  sub(k, x_row_length, B, x_row_length, D, k, S);
  add(k, y_row_length, G, y_row_length, H, k, R);
  strassen_mult_base(k, k, S, k, R, k, P[5], base, inner_func);

  // P6 = (A - C)*(E + F)
  sub(k, x_row_length, A, x_row_length, C, k, S);
  add(k, y_row_length, E, y_row_length, F, k, R);
  strassen_mult_base(k, k, S, k, R, k, P[6], base, inner_func);

  // Z upper left = (P3 + P4) + (P5 - P1)
  add(k, k, P[4], k, P[3], k, S);
  sub(k, k, P[5], k, P[1], k, R);
  add(k, k, S, k, R, z_row_length, Z);

  // Z lower left = P2 + P3
  add(k, k, P[2], k, P[3], z_row_length, Z + k*z_row_length);

  // Z upper right = P0 + P1
  add(k, k, P[0], k, P[1], z_row_length, Z + k);

  // Z lower right = (P0 + P4) - (P2 + P6)
  add(k, k, P[0], k, P[4], k, S);
  add(k, k, P[2], k, P[6], k, R);
  sub(k, k, S, k, R, z_row_length, Z + k*(z_row_length + 1));

  free(R);  // deallocate temp matrices
  free(S);
  for (int i = 6; i >= 0; i--)
    free(P[i]);
}