Как ускорить умножение матриц в C ++? - PullRequest
17 голосов
/ 29 ноября 2010

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

Сравнивая это решение с моим первым со статическими массивами, оно работает в 4 раза медленнее.Что я могу сделать, чтобы ускорить доступ к данным?Я не хочу менять алгоритм.

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

 return c;
}


РЕДАКТИРОВАТЬ
Я исправил свой вопрос avove! Я добавил полный исходный кодниже и попробовал некоторые из ваших советов:
  • поменялся местами k и j итерации цикла -> улучшение производительности
  • объявлено dim() и operator()() как inline -> улучшение производительности
  • передача аргументов по константной ссылке -> потеря производительности! почему?поэтому я им не пользуюсь.

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

Но у меня есть другая проблема: я получаю ошибку памяти в функции mult_strassen(...).Зачем?
terminate called after throwing an instance of 'std::bad_alloc' <br> what(): std::bad_alloc


СТАРАЯ ПРОГРАММА
main.c http://pastebin.com/qPgDWGpW

c99 main.c -o matrix -O3


НОВАЯ ПРОГРАММА
матрица.h http://pastebin.com/TYFYCTY7
matrix.cpp http://pastebin.com/wYADLJ8Y
main.cpp http://pastebin.com/48BSqGJr

g++ main.cpp matrix.cpp -o matrix -O3.


РЕДАКТИРОВАТЬ
Вот некоторые результаты.Сравнение стандартного алгоритма (std), порядка обмена j и k цикла (swap) и заблокированного алгоритма с размером блока 13 (блок).alt text

Ответы [ 7 ]

26 голосов
/ 29 ноября 2010

Говоря о ускорении, ваша функция будет более дружественной к кэшу, если вы поменяете местами порядок циклов 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 остается на месте.

4 голосов
/ 29 ноября 2010

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

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

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

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

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

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

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

1 голос
/ 18 июля 2014

Вот моя реализация алгоритма быстрого простого умножения для матриц с плавающей запятой (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++);
            }
        }
    }
}
1 голос
/ 29 ноября 2010

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

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 примерно в два раза быстрее, чем код C ++ при полной оптимизации.Я не вижу причины в коде.

0 голосов
/ 29 ноября 2010

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

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

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

...