Как ускорить матричное умножение в С++?

Ответ 1

Передайте параметры по ссылке const, чтобы начать с:

matrix mult_std(matrix const& a, matrix const& b) {

Чтобы предоставить вам более подробную информацию, нам нужно знать подробности других используемых методов.
И чтобы ответить, почему исходный метод в 4 раза быстрее, нам нужно будет увидеть оригинальный метод.

Проблема, несомненно, ваша, поскольку эта проблема была решена миллион раз раньше.

Также при запросе этого типа вопроса ВСЕГДА предоставить компилируемый источник с соответствующими входами, чтобы мы могли построить и запустить код и посмотреть, что происходит.

Без кода, который мы просто угадываем.

Изменить

После исправления основной ошибки в исходном коде C (перегрузка буфера)

У меня есть обновленный код для запуска теста бок о бок в справедливом сравнении:

 // INCLUDES -------------------------------------------------------------------
 #include <stdlib.h>
 #include <stdio.h>
 #include <sys/time.h>
 #include <time.h>

 // DEFINES -------------------------------------------------------------------
 // The original problem was here. The MAXDIM was 500. But we were using arrays
 // that had a size of 512 in each dimension. This caused a buffer overrun that
 // the dim variable and caused it to be reset to 0. The result of this was causing
 // the multiplication loop to fall out before it had finished (as the loop was
 // controlled by this global variable.
 //
 // Everything now uses the MAXDIM variable directly.
 // This of course gives the C code an advantage as the compiler can optimize the
 // loop explicitly for the fixed size arrays and thus unroll loops more efficiently.
 #define MAXDIM 512
 #define RUNS 10

 // MATRIX FUNCTIONS ----------------------------------------------------------
 class matrix
 {
 public:
 matrix(int dim)
       : dim_(dim)
 {
         data_ = new int[dim_ * dim_];

 }

     inline int dim() const {
                         return dim_;
                 }
                 inline int& operator()(unsigned row, unsigned col) {
                         return data_[dim_*row + col];
                 }

                 inline int operator()(unsigned row, unsigned col) const {
                         return data_[dim_*row + col];
                 }

 private:
     int dim_;
     int* data_;
 };

// ---------------------------------------------------
 void random_matrix(int (&matrix)[MAXDIM][MAXDIM]) {
         for (int r = 0; r < MAXDIM; r++)
                 for (int c = 0; c < MAXDIM; c++)
                         matrix[r][c] = rand() % 100;
 }
 void random_matrix_class(matrix& matrix) {
         for (int r = 0; r < matrix.dim(); r++)
                 for (int c = 0; c < matrix.dim(); c++)
                         matrix(r, c) = rand() % 100;
 }

 template<typename T, typename M>
 float run(T f, M const& a, M const& b, M& c)
 {
         float time = 0;

         for (int i = 0; i < RUNS; i++) {
                 struct timeval start, end;
                 gettimeofday(&start, NULL);
                 f(a,b,c);
                 gettimeofday(&end, NULL);

                 long s = start.tv_sec * 1000 + start.tv_usec / 1000;
                 long e = end.tv_sec * 1000 + end.tv_usec / 1000;

                 time += e - s;
         }
         return time / RUNS;
 }
 // SEQ MULTIPLICATION ----------------------------------------------------------
  int* mult_seq(int const(&a)[MAXDIM][MAXDIM], int const(&b)[MAXDIM][MAXDIM], int (&z)[MAXDIM][MAXDIM]) {
          for (int r = 0; r < MAXDIM; r++) {
                  for (int c = 0; c < MAXDIM; c++) {
                          z[r][c] = 0;
                          for (int i = 0; i < MAXDIM; i++)
                                  z[r][c] += a[r][i] * b[i][c];
                  }
          }
  }
  void mult_std(matrix const& a, matrix const& b, matrix& z) {
          for (int r = 0; r < a.dim(); r++) {
                  for (int c = 0; c < a.dim(); c++) {
                          z(r,c) = 0;
                          for (int i = 0; i < a.dim(); i++)
                                  z(r,c) += a(r,i) * b(i,c);
                  }
          }
  }

  // MAIN ------------------------------------------------------------------------
  using namespace std;
  int main(int argc, char* argv[]) {
          srand(time(NULL));

          int matrix_a[MAXDIM][MAXDIM];
          int matrix_b[MAXDIM][MAXDIM];
          int matrix_c[MAXDIM][MAXDIM];
          random_matrix(matrix_a);
          random_matrix(matrix_b);
          printf("%d ", MAXDIM);
          printf("%f \n", run(mult_seq, matrix_a, matrix_b, matrix_c));

          matrix a(MAXDIM);
          matrix b(MAXDIM);
          matrix c(MAXDIM);
          random_matrix_class(a);
          random_matrix_class(b);
          printf("%d ", MAXDIM);
          printf("%f \n", run(mult_std, a, b, c));

          return 0;
  }

Теперь результаты:

$ g++ t1.cpp
$ ./a.exe
512 1270.900000
512 3308.800000

$ g++ -O3 t1.cpp
$ ./a.exe
512 284.900000
512 622.000000

Из этого мы видим, что код C примерно в два раза быстрее, чем код С++, когда он полностью оптимизирован. Я не вижу причины в коде.

Ответ 2

Говоря об ускорении, ваша функция будет более удобна для кэширования, если вы поменяете порядок итераций цикла k и j:

matrix mult_std(matrix a, matrix b) {
 matrix c(a.dim(), false, false);
 for (int i = 0; i < a.dim(); i++)
  for (int k = 0; k < a.dim(); k++)
   for (int j = 0; j < a.dim(); j++)  // swapped order
    c(i,j) += a(i,k) * b(k,j);

 return c;
}

Это потому, что индекс k в цикле самого внутреннего цикла приведет к провалу кеша в b на каждой итерации. С j как индексом внутреннего индекса, как c, так и b обращаются к нему, а a остается помещенным.

Ответ 3

Убедитесь, что члены dim() и operator()() объявлены встроенными и что оптимизация компилятора включена. Затем играйте с такими параметрами, как -funroll-loops (on gcc).

Насколько велик a.dim()? Если строка матрицы не вписывается только в пару строк кэша, вам будет лучше с шаблоном доступа к блоку, а не по полной строке в любое время.

Ответ 4

Вы говорите, что не хотите изменять алгоритм, но что это значит?

Развертывает ли счет цикла как "модификацию алгоритма"? Как насчет использования SSE/VMX в зависимости от того, какие инструкции SIMD доступны на вашем процессоре? Как насчет использования какой-либо формы блокировки для улучшения локальности кэша?

Если вы вообще не хотите перестраивать свой код, я сомневаюсь, что вы можете сделать больше, чем изменения, которые вы уже сделали. Все остальное становится компромиссом незначительных изменений в алгоритме для повышения производительности.

Конечно, вы все равно должны взглянуть на asm, сгенерированный компилятором. Это скажет вам гораздо больше о том, что можно сделать для ускорения кода.

Ответ 5

• Используйте SIMD, если сможете. Вы абсолютно должны использовать что-то вроде регистров VMX, если вы делаете обширную векторную математику, предполагая, что используете платформу, способную это сделать, иначе вы понесете огромный успех.
• Не передавать сложные типы, такие как matrix по значению - используйте ссылку const.
• Не вызывать функцию в каждой итерационной кеше dim() вне ваших циклов.
• Хотя компиляторы обычно оптимизируют это эффективно, часто бывает хорошей идеей, чтобы вызывающая сторона предоставляла ссылку на матрицу для вашей функции, а не возвращала матрицу по типу. В некоторых случаях это может привести к дорогостоящей операции копирования.

Ответ 6

Вот моя реализация алгоритма быстрого простого умножения для квадратных матриц с плавающей запятой (2D-массивы). Он должен быть немного быстрее, чем код chrisaycock, так как он поддерживает некоторые приращения.

static void fastMatrixMultiply(const int dim, float* dest, const float* srcA, const float* srcB)
{
    memset( dest, 0x0, dim * dim * sizeof(float) );

    for( int i = 0; i < dim; i++ ) {
        for( int k = 0; k < dim; k++ ) 
        {
            const float* a = srcA + i * dim + k;
            const float* b = srcB + k * dim;
            float* c = dest + i * dim;

            float* cMax = c + dim;
            while( c < cMax ) 
            {   
                *c++ += (*a) * (*b++);
            }
        }
    }
}

Ответ 7

Я беру здесь дикое предположение, но если вы динамически выделяете матрицы, это имеет огромное значение, возможно, проблема заключается в фрагментации. Опять же, я не знаю, как реализуется базовая матрица.

Почему вы не выделяете память для матриц вручную, гарантируя ее непрерывность и самостоятельно создавая структуру указателя?

Кроме того, имеет ли метод dim() дополнительную сложность? Я также объявляю его встроенным.