Быстрое умножение матрицы 4x4 на C

Я пытаюсь найти оптимизированную реализацию C или Assembler для функции, которая умножает две матрицы 4x4 друг на друга. Платформа представляет собой iPhone или iPod на базе ARM6 или ARM7.

В настоящее время я использую довольно стандартный подход - только немного развернутый цикл.

#define O(y,x) (y + (x<<2))

static inline void Matrix4x4MultiplyBy4x4 (float *src1, float *src2, float *dest)
{
    *(dest+O(0,0)) = (*(src1+O(0,0)) * *(src2+O(0,0))) + (*(src1+O(0,1)) * *(src2+O(1,0))) + (*(src1+O(0,2)) * *(src2+O(2,0))) + (*(src1+O(0,3)) * *(src2+O(3,0))); 
    *(dest+O(0,1)) = (*(src1+O(0,0)) * *(src2+O(0,1))) + (*(src1+O(0,1)) * *(src2+O(1,1))) + (*(src1+O(0,2)) * *(src2+O(2,1))) + (*(src1+O(0,3)) * *(src2+O(3,1))); 
    *(dest+O(0,2)) = (*(src1+O(0,0)) * *(src2+O(0,2))) + (*(src1+O(0,1)) * *(src2+O(1,2))) + (*(src1+O(0,2)) * *(src2+O(2,2))) + (*(src1+O(0,3)) * *(src2+O(3,2))); 
    *(dest+O(0,3)) = (*(src1+O(0,0)) * *(src2+O(0,3))) + (*(src1+O(0,1)) * *(src2+O(1,3))) + (*(src1+O(0,2)) * *(src2+O(2,3))) + (*(src1+O(0,3)) * *(src2+O(3,3))); 
    *(dest+O(1,0)) = (*(src1+O(1,0)) * *(src2+O(0,0))) + (*(src1+O(1,1)) * *(src2+O(1,0))) + (*(src1+O(1,2)) * *(src2+O(2,0))) + (*(src1+O(1,3)) * *(src2+O(3,0))); 
    *(dest+O(1,1)) = (*(src1+O(1,0)) * *(src2+O(0,1))) + (*(src1+O(1,1)) * *(src2+O(1,1))) + (*(src1+O(1,2)) * *(src2+O(2,1))) + (*(src1+O(1,3)) * *(src2+O(3,1))); 
    *(dest+O(1,2)) = (*(src1+O(1,0)) * *(src2+O(0,2))) + (*(src1+O(1,1)) * *(src2+O(1,2))) + (*(src1+O(1,2)) * *(src2+O(2,2))) + (*(src1+O(1,3)) * *(src2+O(3,2))); 
    *(dest+O(1,3)) = (*(src1+O(1,0)) * *(src2+O(0,3))) + (*(src1+O(1,1)) * *(src2+O(1,3))) + (*(src1+O(1,2)) * *(src2+O(2,3))) + (*(src1+O(1,3)) * *(src2+O(3,3))); 
    *(dest+O(2,0)) = (*(src1+O(2,0)) * *(src2+O(0,0))) + (*(src1+O(2,1)) * *(src2+O(1,0))) + (*(src1+O(2,2)) * *(src2+O(2,0))) + (*(src1+O(2,3)) * *(src2+O(3,0))); 
    *(dest+O(2,1)) = (*(src1+O(2,0)) * *(src2+O(0,1))) + (*(src1+O(2,1)) * *(src2+O(1,1))) + (*(src1+O(2,2)) * *(src2+O(2,1))) + (*(src1+O(2,3)) * *(src2+O(3,1))); 
    *(dest+O(2,2)) = (*(src1+O(2,0)) * *(src2+O(0,2))) + (*(src1+O(2,1)) * *(src2+O(1,2))) + (*(src1+O(2,2)) * *(src2+O(2,2))) + (*(src1+O(2,3)) * *(src2+O(3,2))); 
    *(dest+O(2,3)) = (*(src1+O(2,0)) * *(src2+O(0,3))) + (*(src1+O(2,1)) * *(src2+O(1,3))) + (*(src1+O(2,2)) * *(src2+O(2,3))) + (*(src1+O(2,3)) * *(src2+O(3,3))); 
    *(dest+O(3,0)) = (*(src1+O(3,0)) * *(src2+O(0,0))) + (*(src1+O(3,1)) * *(src2+O(1,0))) + (*(src1+O(3,2)) * *(src2+O(2,0))) + (*(src1+O(3,3)) * *(src2+O(3,0))); 
    *(dest+O(3,1)) = (*(src1+O(3,0)) * *(src2+O(0,1))) + (*(src1+O(3,1)) * *(src2+O(1,1))) + (*(src1+O(3,2)) * *(src2+O(2,1))) + (*(src1+O(3,3)) * *(src2+O(3,1))); 
    *(dest+O(3,2)) = (*(src1+O(3,0)) * *(src2+O(0,2))) + (*(src1+O(3,1)) * *(src2+O(1,2))) + (*(src1+O(3,2)) * *(src2+O(2,2))) + (*(src1+O(3,3)) * *(src2+O(3,2))); 
    *(dest+O(3,3)) = (*(src1+O(3,0)) * *(src2+O(0,3))) + (*(src1+O(3,1)) * *(src2+O(1,3))) + (*(src1+O(3,2)) * *(src2+O(2,3))) + (*(src1+O(3,3)) * *(src2+O(3,3))); 
};

Я бы выиграл от использования алгоритма Strassen или Coppersmith-Winograd?

Ответ 1

Нет, алгоритм Strassen или Coppersmith-Winograd здесь не имеет большого значения. Они начинают окупаться только для больших матриц.

Если ваше матричное умножение действительно является узким местом, вы можете переписать алгоритм с помощью инструкций NEON SIMD. Это помогло бы только ARMv7, поскольку ARMv6 не имеет этого расширения, хотя.

Я ожидаю, что коэффициент 3 ускорится по скомпилированному C-коду для вашего случая.

EDIT: здесь вы можете найти хорошую реализацию в ARM-NEON: http://code.google.com/p/math-neon/

Для вашего C-кода есть две возможности, которые вы можете сделать для ускорения кода:

Не встраивайте функцию. Ваше умножение матрицы генерирует совсем немного кода при его развертке, и ARM имеет только крошечный кеш команд. Чрезмерная инкрустация может сделать ваш код медленнее, потому что CPU будет занят загрузкой кода в кеш вместо его выполнения.
Используйте ключевое слово ограничения, чтобы сообщить компилятору, что указатели источника и получателя не перекрываются в памяти. В настоящее время компилятор вынужден перезагружать каждое исходное значение из памяти всякий раз, когда результат записывается, поскольку он должен предположить, что источник и получатель могут перекрываться или даже указывать на одну и ту же память.

Ответ 2

Просто задирать. Интересно, почему люди по-прежнему сфакулируют свой код добровольно? C уже трудно читать, нет необходимости добавлять к нему.

static inline void Matrix4x4MultiplyBy4x4 (float src1[4][4], float src2[4][4], float dest[4][4])
{
dest[0][0] = src1[0][0] * src2[0][0] + src1[0][1] * src2[1][0] + src1[0][2] * src2[2][0] + src1[0][3] * src2[3][0]; 
dest[0][1] = src1[0][0] * src2[0][1] + src1[0][1] * src2[1][1] + src1[0][2] * src2[2][1] + src1[0][3] * src2[3][1]; 
dest[0][2] = src1[0][0] * src2[0][2] + src1[0][1] * src2[1][2] + src1[0][2] * src2[2][2] + src1[0][3] * src2[3][2]; 
dest[0][3] = src1[0][0] * src2[0][3] + src1[0][1] * src2[1][3] + src1[0][2] * src2[2][3] + src1[0][3] * src2[3][3]; 
dest[1][0] = src1[1][0] * src2[0][0] + src1[1][1] * src2[1][0] + src1[1][2] * src2[2][0] + src1[1][3] * src2[3][0]; 
dest[1][1] = src1[1][0] * src2[0][1] + src1[1][1] * src2[1][1] + src1[1][2] * src2[2][1] + src1[1][3] * src2[3][1]; 
dest[1][2] = src1[1][0] * src2[0][2] + src1[1][1] * src2[1][2] + src1[1][2] * src2[2][2] + src1[1][3] * src2[3][2]; 
dest[1][3] = src1[1][0] * src2[0][3] + src1[1][1] * src2[1][3] + src1[1][2] * src2[2][3] + src1[1][3] * src2[3][3]; 
dest[2][0] = src1[2][0] * src2[0][0] + src1[2][1] * src2[1][0] + src1[2][2] * src2[2][0] + src1[2][3] * src2[3][0]; 
dest[2][1] = src1[2][0] * src2[0][1] + src1[2][1] * src2[1][1] + src1[2][2] * src2[2][1] + src1[2][3] * src2[3][1]; 
dest[2][2] = src1[2][0] * src2[0][2] + src1[2][1] * src2[1][2] + src1[2][2] * src2[2][2] + src1[2][3] * src2[3][2]; 
dest[2][3] = src1[2][0] * src2[0][3] + src1[2][1] * src2[1][3] + src1[2][2] * src2[2][3] + src1[2][3] * src2[3][3]; 
dest[3][0] = src1[3][0] * src2[0][0] + src1[3][1] * src2[1][0] + src1[3][2] * src2[2][0] + src1[3][3] * src2[3][0]; 
dest[3][1] = src1[3][0] * src2[0][1] + src1[3][1] * src2[1][1] + src1[3][2] * src2[2][1] + src1[3][3] * src2[3][1]; 
dest[3][2] = src1[3][0] * src2[0][2] + src1[3][1] * src2[1][2] + src1[3][2] * src2[2][2] + src1[3][3] * src2[3][2]; 
dest[3][3] = src1[3][0] * src2[0][3] + src1[3][1] * src2[1][3] + src1[3][2] * src2[2][3] + src1[3][3] * src2[3][3]; 
};

Ответ 3

Вы уверены, что ваш развернутый код быстрее, чем явный подход на основе цикла? Имейте в виду, что компиляторы обычно лучше, чем люди, выполняющие оптимизации!

На самом деле, я бы поспорил, что есть больше шансов, что компилятор испустит автоматически инструкции SIMD из хорошо написанного цикла, чем из серии "несвязанных" заявлений...

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

Ответ 4

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

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

Но если вы хотите быстро, я бы использовал инструкции SIMD, если они доступны. Я был бы удивлен, если чипы ARM в эти дни их не будут. Если это так, вы можете управлять всеми продуктами в строке /colum в одной команде; если SIMD имеет ширину 8, вы можете управлять умножением 2 строки/столбца в одной команде. Установка операндов для выполнения этой команды может потребовать некоторых танцев; Инструкции SIMD легко подберут ваши строки (смежные значения), но не будут подбирать столбцы (несмежные). И может потребоваться некоторое количество усилий, чтобы вычислить сумму продуктов в строке/столбце.

Ответ 5

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

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

Пол