Как переписать функцию из последовательной реализации в многопоточную в c ++?

Question

Вопрос о многопоточности. У меня есть последовательно реализованный алгоритм c ++ (основная функция прикреплена внизу), который умножает 2 матрицы, но делает это быстрее, чем обычный алгоритм умножения. Как я могу реализовать многопоточный алгоритм с учетом моего аппаратного обеспечения (2 ядра и 4 потока)?
Я думал, что рекурсивный вызов multiStrassen (p1-p7) может быть выполнен в 7 потоках, но он не может быть очень простым, потому что у меня есть максимум 4 потока, которые могут работать одновременно.

// The following functions are simple, so I did not attach them
// multiply - multiplication of two matrices in the usual way
// summation and subtraction - sum and subtract two matrices
// splitMatrix - divide the matrix into 4 submatrices of the same size each (square matrices)
// collectMatrix - collect a matrix of 4 submatrices (back splitMatrix)

vector<vector<int>> multiStrassen(const vector<vector<int>> & a, const vector<vector<int>> & b, int n) {
  if (n <= 64) {
      return multiply(a, b, n);
  }

  n = n >> 1;

  vector<vector<int>> a11, a12, a21, a22, b11, b12, b21, b22;

  a11.resize(n,vector<int>(n));
  a12.resize(n,vector<int>(n));
  a21.resize(n,vector<int>(n));
  a22.resize(n,vector<int>(n));

  b11.resize(n,vector<int>(n));
  b12.resize(n,vector<int>(n));
  b21.resize(n,vector<int>(n));
  b22.resize(n,vector<int>(n));

  splitMatrix(a, a11, a12, a21, a22);
  splitMatrix(b, b11, b12, b21, b22);

  vector<vector<int>> p1 = multiStrassen(summation(a11, a22), summation(b11, b22), n);
  vector<vector<int>> p2 = multiStrassen(summation(a21, a22), b11, n);
  vector<vector<int>> p3 = multiStrassen(a11, subtraction(b12, b22), n);
  vector<vector<int>> p4 = multiStrassen(a22, subtraction(b21, b11), n);
  vector<vector<int>> p5 = multiStrassen(summation(a11, a12), b22, n);
  vector<vector<int>> p6 = multiStrassen(subtraction(a21, a11), summation(b11, b12), n);
  vector<vector<int>> p7 = multiStrassen(subtraction(a12, a22), summation(b21, b22), n);

  vector<vector<int>> c11 = summation(summation(p1, p4), subtraction(p7, p5));
  vector<vector<int>> c12 = summation(p3, p5);
  vector<vector<int>> c21 = summation(p2, p4);
  vector<vector<int>> c22 = summation(subtraction(p1, p2), summation(p3, p6));

  return collectMatrix(c11, c12, c21, c22);
}

Если на 7 потоках, то это изменится:

vector<vector<int>> p1, p2, p3, p4, p5, p6, p7;
thread th1( [&]() { p1 = multiStrassen(summation(a11, a22), summation(b11, b22), n); } );
thread th2( [&]() { p2 = multiStrassen(summation(a21, a22), b11, n); } );
thread th3( [&]() { p3 = multiStrassen(a11, subtraction(b12, b22), n); } );
thread th4( [&]() { p4 = multiStrassen(a22, subtraction(b21, b11), n); } );
thread th5( [&]() { p5 = multiStrassen(summation(a11, a12), b22, n); } );
thread th6( [&]() { p6 = multiStrassen(subtraction(a21, a11), summation(b11, b12), n); } );
thread th7( [&]() { p7 = multiStrassen(subtraction(a12, a22), summation(b21, b22), n); } );

th1.join();
th2.join();
th3.join();
th4.join();
th5.join();
th6.join();
th7.join();

Это работает быстрее, но я думаю, что это должно быть лучше. Так что вопрос, как сделать лучше?