LINUX.ORG.RU

Перемножение матриц на C


0

0

Есть две программы перемножения матриц. Первая умножает «наивным» алгоритмом ijk, вторая - по ikj (по строкам).

mainijk.c

#include <stdio.h>
#include <stdlib.h>

const int SIZE = 512;

float * generateMatrix(float* a, int n)
{
    int i, j;
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            *(a + i*n + j) = (i - j)*(rand()/10000.0f);
        }
    }
    return a;
}

float * multSquareMatrixIJK(float *a, float *b, float *c, int n)
{
    int i, j, k;
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            for (k = 0; k < n; k++)
            {
                *(c + i*n + j) += *(a + i*n + k) * *(b + k*n + j);
            }
        }
    }
    return c;
}

int main(int argc, char** argv)
{
    float *a = calloc(SIZE*SIZE, sizeof(float));
    float *b = calloc(SIZE*SIZE, sizeof(float));
    float *c = calloc(SIZE*SIZE, sizeof(float));

    generateMatrix(a, SIZE);
    generateMatrix(b, SIZE);
    multSquareMatrixIJK(a, b, c, SIZE);
    
    free(c);
    free(b);
    free(a);
    return 0;
}

mainikj.c

#include <stdio.h>
#include <stdlib.h>

const int SIZE = 512;

float * generateMatrix(float* a, int n)
{
    int i, j;
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            *(a + i*n + j) = (i - j) * (rand()/10000.0f);
        }
    }
    return a;
}

float * multSquareMatrixIKJ(float *a, float *b, float *c, int n)
{
    int i, j, k;
    for (i = 0; i < n; i++)
    {
        for (k = 0; k < n; k++)
        {
            for (j = 0; j < n; j++)
            {
                *(c + i*n + j) += *(a + i*n + k) * *(b + k*n + j);
            }
        }
    }
    return c;
}

int main(int argc, char** argv)
{
    float *a = calloc(SIZE*SIZE, sizeof(float));
    float *b = calloc(SIZE*SIZE, sizeof(float));
    float *c = calloc(SIZE*SIZE, sizeof(float));

    generateMatrix(a, SIZE);
    generateMatrix(b, SIZE);
    multSquareMatrixIKJ(a, b, c, SIZE);
    
    free(c);
    free(b);
    free(a);
    return 0;
}

Внимание, вопрос.

muha@arch:matrix$ time ./ijk
./ijk  9,18s user 0,24s system 93% cpu 10,131 total
muha@arch:matrix$ time ./ikj
./ikj  0,55s user 0,01s system 89% cpu 0,619 total

Почему так, если фактически обе программы отличаются лишь переставленными циклами (т.е., количество циклов одно и то же).

В первом случае есть

> *(b + k*n + j);

во внутреннем цикле по k. Указатель прыгает далеко и делает неэффективным процессорный кэш.

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

zzo
()
Ответ на: комментарий от zzo

То есть, это связано с тем, что в С массивы хранятся по строкам? А если бы они хранились по столбцам - то оба варианта были бы одинаково медленными, как я понимаю?

rexadecimal
() автор топика
Ответ на: комментарий от rexadecimal

> А если бы они хранились по столбцам - то оба варианта были бы одинаково медленными, как я понимаю?
Неправильно ты считаешь. Если бы по столбцам - просто обменялись бы проги своим временем.
Просто процессор работает не напрямую с оперативкой, а через свой кеш. Кеш отображает себе некоторый отрезок в RAM и работает намного быстрее. Если ты обращаешься к памяти, которая выходит за пределы кеша - он загружается заново из RAM. Так вот в первом случае ты "гуляешь" по сильно разным адресам в RAM, каждый раз заставляя загружать данные в кеш, а во втором случае идешь последовательно и кеш грузится только data_size/cache_size раз.

snizovtsev ★★★★★
()

Если компилятор достаточно умный, может свернуть \sum a_ik b_kj в a_ik \sum b_kj, уменьшив количество умножений c n^3 до n^2 (gcc умеет это делать по крайней мере в простых случаях).

grob ★★★★★
()
Ответ на: комментарий от grob

За разработку алгоритма умножения матриц, который делает n^2 умножений надо сразу звание академика давать. Самые лучше известные алгоритмы работают как n^{2.8}.

Reset ★★★★★
()
Ответ на: комментарий от Reset

> Самые лучше известные алгоритмы работают как n^{2.8}.

около 2.3

2.8 это Штрассен

dilmah ★★★★★
()
Ответ на: комментарий от Reset

да я уже понял, что х-ню сказал, левая часть меняется,

grob ★★★★★
()
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.