Для тех, кто обучался/обучается на математических или программистских факультетах вузов, я думаю, эта статья будет не в новость, но стало самому интересно протестировать скорость работы разных алгоритмов. Также её можно рассматривать, как некое пособие по оптимизации, но такую оптимизацию стоит проводить только, когда это действительно необходимо, т.к. читаемость кода рушится на глазах, да и отлаживать такое намного сложнее.
Наверняка большинству будет лень читать всю статью, но советую промотать вниз и почитать выводы — помойму интересные цифры там.
Итак задача: перемножить две большие матрицы double-ов (размерами третьего порядка). Для простоты будем рассматривать квадратные матрицы, хотя все алгоритмы подойдут и для прямоугольных. Алгоритм писался на C++, но классов нигде не использовал, так что можно считать код C-совместимым (возможно только cout использовал).
Не буду объяснять тут, что такое матрица и как их перемножать — тем, кто этого не знает, вряд ли будет интересно как ускорить перемножение…
Для начала определимся, что матрицы будем хранить в одномерном массиве размера n*n. К элементу A(i, j) будем обращаться как a[i*n+j]. Хранить в одномерном быстрее чем в двумерном массиве, т.к. одно целочисленное умножение и одно обращение к памяти работает быстрее, чем два обращения к памяти.
Первый самый логичный способ перемножения
Заметим, что во первых мы в цикле слишком много лишний раз умножаем i*n. Каждый раз обращаться к памяти прибавляя к c[i*n+j] необязательно — можно все это складывать в переменную, а уже потом приравнять её c[i*n+j]. Заметим еще, что сейчас компиляторы очень хорошо оптимизируют наш код, но проблемы возникают в точках ветвления (if, for и т.п.), для улучшения производительности нужно стремиться увеличить глубину конвейера инструкций, т.е. линейной части алгоритма (когда вычисления идут подряд, без условий или циклов), поэтому развернем наш основной цикл так, чтобы внутри третьего цикла проходило не одно, а 4 умножения. Ну и постараемся как можно реже обращаться к памяти и умножать, а сохранять все в локальных переменных и указателях, что также ускоряет работу. Делаем дубляж кода для четных и нечетных размеров матриц, чтобы не проверять на четность внутри цикла, а проверить только один раз. Получаем вот какую функцию:
Ну и последний штрих. Когда мы идем по матрице в цикле — нам так или иначе приходится брать информацию из оперативной памяти (особенно когда считываем по столбцу). Из оперативной памяти в кеш информация сохраняется строками (т.е. когда мы читаем какой-то элемент — в кэш кладется несколько следующих за ним элементов в памяти), а когда мы идем по столбцу — нам приходится каждый раз лазить в память. Чтобы такого не случалось, а работать чаще всего приходилось с кэшем, а не с оперативной памятью сделаем следующее: разобьем нашу матрицу на подматрицы размера примерно 100 на 100 (размер задаётся в аргументе функции). Если поделилось не нацело — не страшно, оставляем прямоугольные куски снизу и справа. Назовем наши подматрицы A(1,1),..., A(n,n) соответственно. И будем перемножать матрицу из подматриц по обычному правилу. Не сложно математически доказать, что ответ получится такой же как и при обычном перемножении, но теперь когда мы будем перемножать две маленькие матрицы — они обе поместятся в кэш-память и перемножатся намного быстрее. Еще для перемножения нам понадобятся функции «взять подматрицу» и «положить подматрицу добавив к тому что было». Функции обнуления матрицы, а также функция перемножения прямоугольных матриц. Т.к. прямоугольных умножений у нас не очень много, то не стал её оптимизировать так же, как и для квадратных. Сейчас кажется, что правильнее было конечно написать одну хорошую функцию для прямоугольных, а для квадратных сделать просто алиас, не помню почему, но я этого не сделал. Итак пишем окончательную функцию:
Исходной матрицей возьмем матрицу Гильберта, т.к. все элементы меньше единицы и очень удобно умножать (не будет переполнения).
Полный текст программы тут или тут
Файл matrix_test.cpp компилировался так:
Вот результаты тестирования на моих компьютерах (размер подматриц — 200 для iMac и Toshiba вроде бы оптимальный, для Athlon оптимальнее 115, но даже без подматриц он показал себя не с лучшей стороны):
Тестовые компьютеры:
Как видим прирост производительности колоссальный (ускорение более чем в 10 раз на самом современном процессоре). Стоит учесть, что это все считалось в одном процессе, а значит двухядерность Core2Duo не сыграла роли. В дальнейшем конечно правильно будет с помощью multithread или mpi распараллелить программу и получить прирост еще в 1.8-1.9 раза. Также стоит отметить, что прогресс не стоит на месте, и более современные процессоры на голову обгоняют модели 5-летней давности. Вобщем как видите код получился ужасен, но это стоит того. И важен не сам метод перемножения матриц, а возможные техники оптимизации программы.
Наверняка большинству будет лень читать всю статью, но советую промотать вниз и почитать выводы — помойму интересные цифры там.
Итак задача: перемножить две большие матрицы double-ов (размерами третьего порядка). Для простоты будем рассматривать квадратные матрицы, хотя все алгоритмы подойдут и для прямоугольных. Алгоритм писался на C++, но классов нигде не использовал, так что можно считать код C-совместимым (возможно только cout использовал).
Не буду объяснять тут, что такое матрица и как их перемножать — тем, кто этого не знает, вряд ли будет интересно как ускорить перемножение…
Для начала определимся, что матрицы будем хранить в одномерном массиве размера n*n. К элементу A(i, j) будем обращаться как a[i*n+j]. Хранить в одномерном быстрее чем в двумерном массиве, т.к. одно целочисленное умножение и одно обращение к памяти работает быстрее, чем два обращения к памяти.
Первый самый логичный способ перемножения
int multMatrixSqBad(double *a, const int n, double *b, double *c){
int i, j, k;
for (i=0; i<n; i++){
for (j=0; j<n; j++){
c[i*n+j]=0;
for (k=0; k<n; k++){
c[i*n+j]+=a[i*n+k]*b[k*n+j];
}
}
}
return 0;
}
Заметим, что во первых мы в цикле слишком много лишний раз умножаем i*n. Каждый раз обращаться к памяти прибавляя к c[i*n+j] необязательно — можно все это складывать в переменную, а уже потом приравнять её c[i*n+j]. Заметим еще, что сейчас компиляторы очень хорошо оптимизируют наш код, но проблемы возникают в точках ветвления (if, for и т.п.), для улучшения производительности нужно стремиться увеличить глубину конвейера инструкций, т.е. линейной части алгоритма (когда вычисления идут подряд, без условий или циклов), поэтому развернем наш основной цикл так, чтобы внутри третьего цикла проходило не одно, а 4 умножения. Ну и постараемся как можно реже обращаться к памяти и умножать, а сохранять все в локальных переменных и указателях, что также ускоряет работу. Делаем дубляж кода для четных и нечетных размеров матриц, чтобы не проверять на четность внутри цикла, а проверить только один раз. Получаем вот какую функцию:
int multMatrixSq(double *a, const int n, double *b, double *c){
double s1=0, s2=0, s3=0, s4=0;
double *f, *g, *h;
int i, j, k;
if (n%2==0){
for (i = 0; i < n-1; i+=2)
{
f=a+i*n;
for (j = 0; j < n-1; j+=2)
{
g=b+j;
for (k=0;k<n;k++){
s1+=f[k]*g[n*k];
s2+=f[n+k]*g[n*k];
s3+=f[k]*g[n*k+1];
s4+=f[n+k]*g[n*k+1];
}
h=c+i*n+j;
h[0]=s1; s1=0;
h[n]=s2; s2=0;
h[1]=s3; s3=0;
h[n+1]=s4; s4=0;
}
}
}
else {
for (i = 0; i < n-1; i+=2)
{
f=a+i*n;
for (j = 0; j < n-1; j+=2)
{
g=b+j;
for (k=0;k<n;k++){
s1+=f[k]*g[n*k];
s2+=f[n+k]*g[n*k];
s3+=f[k]*g[n*k+1];
s4+=f[n+k]*g[n*k+1];
}
h=c+i*n+j;
h[0]=s1; s1=0;
h[n]=s2; s2=0;
h[1]=s3; s3=0;
h[n+1]=s4; s4=0;
}
if (j==n-1){
g=b+j;
for (k=0;k<n;k++){
s1+=f[k]*g[n*k];
s2+=f[n+k]*g[n*k];
}
h=c+i*n+j;
h[0]=s1; s1=0;
h[n]=s2; s2=0;
}
}
if (i==n-1){
f=a+i*n;
for (j = 0; j < n-1; j+=2)
{
g=b+j;
for (k=0;k<n;k++){
s1+=f[k]*g[n*k];
s3+=f[k]*g[n*k+1];
}
h=c+i*n+j;
h[0]=s1; s1=0;
h[1]=s3; s3=0;
}
if (j==n-1){
g=b+j;
for (k=0;k<n;k++){
s1+=f[k]*g[n*k];
}
h=c+i*n+j;
h[0]=s1; s1=0;
}
}
}
return 0;
}
Ну и последний штрих. Когда мы идем по матрице в цикле — нам так или иначе приходится брать информацию из оперативной памяти (особенно когда считываем по столбцу). Из оперативной памяти в кеш информация сохраняется строками (т.е. когда мы читаем какой-то элемент — в кэш кладется несколько следующих за ним элементов в памяти), а когда мы идем по столбцу — нам приходится каждый раз лазить в память. Чтобы такого не случалось, а работать чаще всего приходилось с кэшем, а не с оперативной памятью сделаем следующее: разобьем нашу матрицу на подматрицы размера примерно 100 на 100 (размер задаётся в аргументе функции). Если поделилось не нацело — не страшно, оставляем прямоугольные куски снизу и справа. Назовем наши подматрицы A(1,1),..., A(n,n) соответственно. И будем перемножать матрицу из подматриц по обычному правилу. Не сложно математически доказать, что ответ получится такой же как и при обычном перемножении, но теперь когда мы будем перемножать две маленькие матрицы — они обе поместятся в кэш-память и перемножатся намного быстрее. Еще для перемножения нам понадобятся функции «взять подматрицу» и «положить подматрицу добавив к тому что было». Функции обнуления матрицы, а также функция перемножения прямоугольных матриц. Т.к. прямоугольных умножений у нас не очень много, то не стал её оптимизировать так же, как и для квадратных. Сейчас кажется, что правильнее было конечно написать одну хорошую функцию для прямоугольных, а для квадратных сделать просто алиас, не помню почему, но я этого не сделал. Итак пишем окончательную функцию:
int multMatrixBlock(double *a, const int n, double *b, const int m, double *res, double *c, double *d, double*e){
int n_block=n/m;
int m_small;
int i,j,k,o;
m_small=n-n_block*m;
clearMatrix(res, n*n);
for ( i = 0; i <= n_block; i++)
{
for ( j = 0; j <= n_block; j++)
{
for ( k = 0; k <=n_block; k++)
{
if (i<n_block && j<n_block && k<n_block) {
getSubMatrix(a, n, c, m, i, k);
getSubMatrix(b, n, d, m, k, j);
multMatrixSq(c, m, d, e);
setSubMatrixAdded(res, n, e, m, i, j);
}
else if (i<n_block && j<n_block && k==n_block) {
getSubMatrix(a, n, c, m, i, k);
getSubMatrix(b, n, d, m, k, j);
multMatrixRect(c, m, d,m_small, m, e);
setSubMatrixAdded(res, n, e, m, i, j);
}
else if (i==n_block && j<n_block && k<n_block) {
getSubMatrix(a, n, c, m, i, k);
getSubMatrix(b, n, d, m, k, j);
multMatrixRect(c, m_small, d,m, m, e);
setSubMatrixAdded(res, n, e, m, i, j);
}
else if (i==n_block && j<n_block && k==n_block) {
getSubMatrix(a, n, c, m, i, k);
getSubMatrix(b, n, d, m, k, j);
multMatrixRect(c, m_small, d,m_small, m, e);
setSubMatrixAdded(res, n, e, m, i, j);
}
else if (i<n_block && j==n_block && k<n_block) {
getSubMatrix(a, n, c, m, i, k);
getSubMatrix(b, n, d, m, k, j);
multMatrixRect(c, m, d,m, m_small, e);
setSubMatrixAdded(res, n, e, m, i, j);
}
else if (i<n_block && j==n_block && k==n_block) {
getSubMatrix(a, n, c, m, i, k);
getSubMatrix(b, n, d, m, k, j);
multMatrixRect(c, m, d,m_small, m_small, e);
setSubMatrixAdded(res, n, e, m, i, j);
}
else if (i==n_block && j==n_block && k<n_block) {
getSubMatrix(a, n, c, m, i, k);
getSubMatrix(b, n, d, m, k, j);
multMatrixRect(c, m_small, d,m, m_small, e);
setSubMatrixAdded(res, n, e, m, i, j);
}
else if (i==n_block && j==n_block && k==n_block) {
getSubMatrix(a, n, c, m, i, k);
getSubMatrix(b, n, d, m, k, j);
multMatrixSq(c, m_small, d, e);
setSubMatrixAdded(res, n, e, m, i, j);
}
}
}
}
return 0;
}
Исходной матрицей возьмем матрицу Гильберта, т.к. все элементы меньше единицы и очень удобно умножать (не будет переполнения).
double formula(int i, int j, int n){
return 1/(double(i+j+1));
}
Полный текст программы тут или тут
Файл matrix_test.cpp компилировался так:
g++ matrix_test.cpp -O3
Вот результаты тестирования на моих компьютерах (размер подматриц — 200 для iMac и Toshiba вроде бы оптимальный, для Athlon оптимальнее 115, но даже без подматриц он показал себя не с лучшей стороны):
Тестовые компьютеры:
Как видим прирост производительности колоссальный (ускорение более чем в 10 раз на самом современном процессоре). Стоит учесть, что это все считалось в одном процессе, а значит двухядерность Core2Duo не сыграла роли. В дальнейшем конечно правильно будет с помощью multithread или mpi распараллелить программу и получить прирост еще в 1.8-1.9 раза. Также стоит отметить, что прогресс не стоит на месте, и более современные процессоры на голову обгоняют модели 5-летней давности. Вобщем как видите код получился ужасен, но это стоит того. И важен не сам метод перемножения матриц, а возможные техники оптимизации программы.
