Оптимизация умножения матриц с помощью транспонирования матриц
Я работаю над заданием, в котором я переставляю матрицу, чтобы уменьшить потери в кеше для операции умножения матриц. Из того, что я понимаю из одноклассников, я должен получить улучшение в 8 раз. Тем не менее, я получаю только 2 раза … что я могу делать не так?
void transpose(int size, matrix m) {
int i, j;
for (i = 0; i < size; i++)
for (j = 0; j < size; j++)
std::swap(m.element[i][j], m.element[j][i]);
}
void mm(matrix a, matrix b, matrix result) {
int i, j, k;
int size = a.size;
long long before, after;
before = wall_clock_time();
// Do the multiplication
transpose(size, b); // transpose the matrix to reduce cache miss
for (i = 0; i < size; i++)
for (j = 0; j < size; j++) {
int tmp = 0; // save memory writes
for(k = 0; k < size; k++)
tmp += a.element[i][k] * b.element[j][k];
result.element[i][j] = tmp;
}
after = wall_clock_time();
fprintf(stderr, "Matrix multiplication took %1.2f seconds\n", ((float)(after - before))/1000000000);
}
До сих пор я все делаю правильно?
К вашему сведению: следующая оптимизация, которую мне нужно сделать, это использовать SIMD / Intel SSE3
Решение
До сих пор я все делаю правильно?
У тебя проблемы с транспонированием. Вы должны были увидеть эту проблему, прежде чем начать беспокоиться о производительности. Когда вы делаете какие-либо взломы для оптимизации это всегда хорошая идея использовать наивную, но неоптимальную реализацию в качестве теста. Оптимизация, достигающая 100-кратного ускорения, бесполезна, если она не дает правильного ответа.
Еще одна оптимизация, которая поможет, это перейти по ссылке. Вы передаете копии. На самом деле, ваш matrix result может никогда не выйти, потому что вы передаете копии. Еще раз, вы должны были проверить.
Еще одна оптимизация, которая поможет ускорению, заключается в кэшировании некоторых указателей. Это все еще довольно медленно:
for(k = 0; k < size; k++)
tmp += a.element[i][k] * b.element[j][k];
result.element[i][j] = tmp;
Оптимизатор может обойти проблемы с указателями, но, вероятно, нет. По крайней мере, если вы не используете нестандартный __restrict__ ключевое слово, чтобы сообщить компилятору, что ваши матрицы не перекрываются. Указатели кэша, поэтому вам не нужно делать a.element[i], b.element[j], а также result.element[i], И все же это может помочь сказать компилятору, что эти массивы не пересекаются с __restrict__ ключевое слово.
добавление
После просмотра кода ему нужна помощь. Небольшой комментарий первым. Вы не пишете C ++. Ваш код на C с небольшим намеком на C ++. Вы используете struct скорее, чем class, malloc скорее, чем new, typedef struct а не просто structЗаголовки C, а не C ++.
Из-за вашей реализации вашего struct matrixМой комментарий о медлительности из-за конструкторов копирования был неверным. То, что это было неправильно, еще хуже! Использование неявно определенного конструктора копирования в сочетании с классами или структурами, которые содержат голые указатели, играет с огнем. Вы очень сильно обожжетесь, если кто-то позвонит m(a, a, a_squared) получить квадрат матрицы a, Вы обожгетесь еще хуже, если кто-то ожидает m(a, a, a) сделать вычисление на месте a2.
Математически ваш код покрывает лишь небольшую часть проблемы умножения матриц. Что если кто-то захочет умножить матрицу 100×1000 на матрицу 1000×200? Это совершенно верно, но ваш код не справляется с этим, потому что ваш код работает только с квадратными матрицами. С другой стороны, ваш код позволит кому-то умножить матрицу 100×100 на матрицу 200×200, что не имеет смысла.
Конструктивно ваш код имеет почти 100% гарантию того, что он будет медленным из-за того, что вы используете рваные массивы. malloc может разбрызгивать строки ваших матриц по всей памяти. Вы получите намного лучшую производительность, если матрица внутренне представлена в виде непрерывного массива, но доступ к ней осуществляется так, как если бы это была матрица NxM. C ++ предоставляет несколько хороших механизмов для этого.
Другие решения
Если ваше назначение подразумевает, что вы ДОЛЖНЫ транспонировать, тогда, конечно, вы должны исправить свою процедуру транспонирования. В существующем состоянии он выполняет транспонирование ДВА раза, в результате чего транспонирование вообще не происходит. Цикл j = не должен читать
j=0; j<size; j++
но
j=0; j<i; j++
Транспонирование не является необходимым, чтобы избежать обработки элементов одной из матриц-факторов в «неправильном» порядке. Просто поменяйте местами J-петлю и K-петлю. Оставляя в стороне любую (другую) настройку производительности, базовая структура цикла должна быть:
for (int i=0; i<size; i++)
{
for (int k=0; k<size; k++)
{
double tmp = a[i][k];
for (int j=0; j<size; j++)
{
result[i][j] += tmp * b[k][j];
}
}
}