Является ли моя реализация Strassen медленнее, чем наивное умножение матриц? - PullRequest
Новая
       26

Является ли моя реализация Strassen медленнее, чем наивное умножение матриц?

0 голосов
/ 28 апреля 2020

Итак, я реализовал алгоритм Штрассена в C, и даже для больших матриц он не может превзойти алгоритм classi c ikj. Я использую только Матрицы степеней 2, заполненные двойными числами. Я читал, что для «больших» n алгоритм Штрассена превосходит класс c, например, для n> = 1000. Для меня это не так. Я не использую точку отсечения, поскольку она не является частью исходного алгоритма. Если я это сделаю, то на самом деле он быстрее классического c, но без этого он работает намного хуже. Является ли это ожидаемым поведением, и источники, которые я нашел в Интернете об алгоритме Штрассена, на самом деле используют точку отсечения?

Я пытался оптимизировать Штрассен, используя только две матрицы, хранящиеся в виде одномерного массива, передавая их в качестве аргументов и на самом деле не создание подматриц, если это возможно, а просто передача нового размера и указателя на две ячейки, чтобы я не выделял столько памяти. Но, возможно, для больших Матриц это на самом деле плохо, так как я обращаюсь ко всей матрице c и мне нужно хранить ее где-то в основной памяти.

Вот ссылка на Github: https://github.com/AlexLoitzl/matmul Код еще не задокументирован должным образом ... А вот реализация Штрассена: 

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]);
}
...