Оптимизация вычисления матрицы в C

Question

Я пытаюсь оптимизировать некоторый написанный на С код, который умножает Матрицу на ее транспонирование и хотел бы знать, может ли кто-нибудь увидеть что-нибудь еще, что я могу сделать, чтобы сделать мой метод более эффективным с точки зрения времени / часов вычислений циклы.

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

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

Неизмененный пример кода:

void print_unmodified()
{
    unsigned int n = 64;
    unsigned int A[N_MAX];  // N_MAX = 64*64
    unsigned int B[N_MAX];

    // Initialising the A Matrix with values 1->64^2
    for (int i = 0; i < (n * n); i++)
    {
        A[i] = i + 1;
    }

    // Initialising the B Matrix with zeros
    for (int i = 0; i < (n * n); i++)
    {
        B[i] = 0;
    }

    // Matrix Multiplicaiton B = A*A'
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            for (int k = 0; k < n; k++)
            {
                B[i + n * j] = B[i + n * j] + A[i + n * k] * A[j + n * k];
            }
        }
    }
}

Мой модифицированный код:

void print_modified()
{
    unsigned int n = 64;
    unsigned int A[N_MAX];
    unsigned int B[N_MAX];
    unsigned int temp = 0;

    // Initialising the A Matrix with values 1->64^2
    for (int i = 0; i < (n * n); i++)
    {
        A[i] = i + 1;
    }

    // Matrix Multiplicaiton B = A*A'
    for (int i = 0; i < n; i++)
    {   

        for (int j = 0; j < n; j++)
        {
            temp = 0;

            for (int k = 0; k < (n*n); k+=n)
            {
                temp += A[j + k] * A[i + k];
            }

            B[j + n*i] = temp;
        }
    }
}

Оба метода заканчиваются полученной B-матрицей, поэтому я знаю, что мой текущий метод математически корректен.

РЕДАКТИРОВАТЬ: исправлен тип - я не хочу умножать А на его обратное, но умножить А на его транспонирование.

Answer 1

Размещенный код дает правильный ответ, матрица B , которая является произведением другой матрицы A путем ее транспонирования, но обе матрицы хранятся в паре массивов в порядок столбцов-мажоров (или, по крайней мере, это то, что ОП указывает в комментариях).

Учитывая этот фрагмент:

for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++) {
        for (int k = 0; k < n; k++) {
            B[i + n * j] = B[i + n * j] + A[i + n * k] * A[j + n * k];
        //                                  ^^^^^^^^^      ^^^^^^^^^
        }
    }
}

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

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

// Calculate the matrix B = A times A transposed. A is an n x n matrix.
// A and B are stored in plain arrays ('src' and 'dest') in column-major order 
void mult_AAT(size_t n, value_t const *src, value_t *dest)
{
    for (size_t i = 0, end = n * n; i < end; i++)
    {
        dest[i] = 0;
    }
    for (size_t k = 0; k < n; k++)         // <--
    {
        for (size_t j = 0; j < n; j++)     // <-
        {
            for (size_t i = 0; i < n; i++)
            {
                dest[j * n + i] += src[k * n + i] * src[k * n + j];
                //   ^^^^^^^^^         ^^^^^^^^^        ^^^^^^^^^     
            }
        }
    }
}

Live ЗДЕСЬ .

Answer 2

Вы можете улучшить свое первое решение, инициализировав матрицу A и B в том же цикле for.

n * k также рассчитывается два раза, вы можете сохранить его в переменнойчтобы сэкономить время.

Лучше использовать B[i + n * j] += ... вместо B[i + n * j] = B[i + n * j] + ..., потому что в первом случае B[i + n * j] читается один раз, а во втором - дважды.

void print_unmodified()
{
    unsigned int n = 64u;
    unsigned int A[N_MAX];
    unsigned int B[N_MAX];

    /* Initializing the A and B matrix with values 1->64^2 */
    for (int i = 0; i < (n * n); i++)
    {
        A[i] = i + 1;
        B[i] = 0u;
    }

    /* Matrix Multiplication B = A*A' */
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            int index = i + n*j;
            for (int k = 0; k < n; k++)
            {
                int p = n * k;
                B[index] += A[i + p] * A[j + p];
            }
        }
    }
}