Наиболее эффективный способ вычисления экспоненты каждого элемента матрицы

Question

Я мигрирую из Matlab в C + GSL и хотел бы знать , каков наиболее эффективный способ вычисления матрицы B, для которого:

B[i][j] = exp(A[i][j])

где i в [0, Ny] и j в [0, Nx].

Обратите внимание, что это отличается от экспоненциальной матрицы:

B = exp(A)

, что может быть выполнено с помощью некоторого нестабильного / неподдерживаемого кода в GSL (linalg.h).

Я только что нашел решение для грубой силы (пара циклов for), но есть ли более разумный способ сделать это?

EDIT

Результаты поста решения Дрю Холла

Все результаты получены из цикла 1024x1024 for(for), в котором на каждой итерации назначаются два значения double (комплексное число). Время - это усредненное время за 100 казней .

• Результаты при учете режима {Row, Column} -Major для хранения матрицы:
  • 226,56 мс при циклическом цикле по строке во внутреннем цикле в режиме Row-Major (вариант 1).
  • 223,22 мс при циклическом перемещении по столбцу во внутреннем цикле в режиме Row-Major (вариант 2).
  • 224.60 мс при использовании функции gsl_matrix_complex_set, предоставляемой GSL (случай 3).

Исходный код для дела 1 :

for(i=0; i<Nx; i++)
{
    for(j=0; j<Ny; j++)
    {
        /* Operations to obtain c_value (including exponentiation) */
        matrix[2*(i*s_tda + j)] = GSL_REAL(c_value);
        matrix[2*(i*s_tda + j)+1] = GSL_IMAG(c_value);
    }
}

Исходный код для дела 2 :

for(i=0; i<Nx; i++)
{
    for(j=0; j<Ny; j++)
    {
        /* Operations to obtain c_value (including exponentiation) */
        matrix->data[2*(j*s_tda + i)] = GSL_REAL(c_value);
        matrix->data[2*(j*s_tda + i)+1] = GSL_IMAG(c_value);
    }
}

Исходный код для дела 3 :

for(i=0; i<Nx; i++)
{
    for(j=0; j<Ny; j++)
    {
        /* Operations to obtain c_value (including exponentiation) */
        gsl_matrix_complex_set(matrix, i, j, c_value);
    }
}

Answer 1

Нет способа избежать итерации по всем элементам и вызова exp() или эквивалентного для каждого. Но есть более быстрые и медленные способы итерации.

В частности, ваша цель должна заключаться в том, чтобы минимизировать количество кеш-пропусков. Выясните, хранятся ли ваши данные в порядке основной строки или основной столбца, и убедитесь, что ваши циклы расположены так, что внутренний цикл выполняет итерации по элементам, хранящимся в памяти, а внешний цикл переносит большой шаг к следующей строке если основной ряд) или столбец (если основной столбец). Хотя это кажется тривиальным, оно может ОГРОМНО повлиять на производительность (в зависимости от размера вашей матрицы).

Как только вы обработаете кеш, ваша следующая цель - удалить издержки цикла. Первый шаг (если ваш матричный API поддерживает это) - перейти от вложенных циклов (границы M & N) к одному циклу, проходящему по основным данным (граница M N). Для этого вам нужно получить необработанный указатель на базовый блок памяти (то есть двойной , а не двойной **).

Наконец, добавьте некоторое развертывание цикла (то есть сделайте 8 или 16 элементов для каждой итерации цикла), чтобы еще больше сократить издержки цикла, и это, вероятно, настолько быстро, насколько вы можете это сделать. Вероятно, вам понадобится последний оператор switch с откатом, чтобы очистить остальные элементы (когда ваш размер массива% block size! = 0).

Answer 2

Нет, если только у меня нет какой-то странной математической причуды, о которой я не слышал, вам просто нужно пройтись по элементам с двумя циклами for.

Answer 3

Если вы просто хотите применить exp к массиву чисел, ярлык на самом деле не существует. Ты должен это назвать (Nx * Ny) раз. Если некоторые из элементов матрицы простые, например 0, или есть повторяющиеся элементы, может помочь некоторая запоминание.

Однако, если вам действительно нужна экспоненциальная матрица (что очень полезно), алгоритм, на который мы опираемся, - DGPADM . Он в Фортране, но вы можете использовать f2c , чтобы преобразовать его в C. Вот бумага на нем.

Answer 4

Поскольку содержимое цикла не было показано, бит, который вычисляет значение c_value, мы не знаем, ограничена ли производительность кода пропускной способностью памяти или ЦП.Единственный способ узнать наверняка, это использовать профилировщик, и при этом сложный.Он должен быть в состоянии измерить задержку памяти, то есть количество времени, в течение которого процессор простаивал, ожидая поступления данных из оперативной памяти.

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

Если вы ограничены процессором, тогда вам доступно еще несколько вариантов.Использование SIMD является одним из вариантов, как и ручное кодирование кода с плавающей запятой (компилятор C / C ++ не очень хорош в коде FPU по многим причинам).Если бы это был я, и код во внутреннем цикле позволяет это сделать, я бы имел два указателя в массиве, один в начале и второй на 4/5 пути.Каждая итерация, операция SIMD будет выполняться с использованием первого указателя и скалярных операций FPU с использованием второго указателя, так что каждая итерация цикла дает пять значений.Затем я бы чередовал инструкции SIMD с инструкциями FPU, чтобы уменьшить затраты времени ожидания.Это не должно влиять на ваши кеши, так как (по крайней мере, на Pentium) MMU может одновременно передавать до четырех потоков данных (т.е. предварительно выбирать данные для вас без каких-либо подсказок или специальных инструкций).