Я хочу рассчитать абсолютные значения элементов сложного массива в C или С++. Самый простой способ -
for(int i = 0; i < N; i++)
{
b[i] = cabs(a[i]);
}
Но для больших векторов, которые будут медленными. Есть ли способ ускорить это (например, используя распараллеливание)? Язык может быть C или С++.
Учитывая, что все итерации цикла независимы, вы можете использовать следующий код для распараллеливания:
#pragma omp parallel for
for(int i = 0; i < N; i++)
{
b[i] = cabs(a[i]);
}
Конечно, для этого вы должны включить поддержку OpenMP при компиляции кода (обычно с помощью флага /openmp или установки параметров проекта).
Вы можете найти несколько примеров использования OpenMP в wiki.
Используйте векторные операции.
Если у вас есть glibc 2.22 (довольно недавно), вы можете использовать SIMD-возможности OpenMP 4.0 для работать с векторами/массивами.
Libmvec - это векторная математическая библиотека, добавленная в Glibc 2.22.
Была добавлена векторная математическая библиотека для поддержки SIMD-конструкций OpenMP4.0 (# 2.8 в http://www.openmp.org/mp-documents/OpenMP4.0.0.pdf), добавив векторных реализаций векторных математических функций.
Векторные математические функции являются векторными вариантами соответствующей скалярной математики операции, реализованные с использованием SIMD ISA-расширений (например, SSE или AVX для x86_64). Они принимают упакованные векторные аргументы, выполняют операцию по каждый элемент упакованного векторного аргумента и вернуть упакованный вектор результат. Использование векторных математических функций выполняется быстрее, чем повторное вызов скалярные математические процедуры.
Также см. Параллель для vs omp simd: когда использовать каждый?
Если вы работаете в Solaris, вы можете явно использовать vhypot() из библиотеки векторов math libmvec.so для работы на вектор комплексных чисел, чтобы получить абсолютное значение каждого из них:
Описание
Эти функции оценивают функцию hypot (x, y) для целого вектора значений сразу....
Исходный код для libmvec можно найти в http://src.illumos.org/source/xref/illumos-gate/usr/src/lib/libmvec/ и код vhypot(), в частности, в http://src.illumos.org/source/xref/illumos-gate/usr/src/lib/libmvec/common/__vhypot.c Я не помню, предоставила ли Sun Microsystems версию Linux libmvec.so или нет.
Или используйте Concurrency:: parallele_for так:
Concurrency::parallel_for(0, N, [&a, &b](int i)
{
b[i] = cabs(a[i]);
});
Если вы используете современный компилятор (например, GCC 5), вы можете использовать Cilk +, который даст вам отличную нотацию массива, автоматически используя инструкции SIMD и распараллеливание.
Итак, если вы хотите запустить их параллельно, выполните следующие действия:
#include <cilk/cilk.h>
cilk_for(int i = 0; i < N; i++)
{
b[i] = cabs(a[i]);
}
или если вы хотите протестировать SIMD:
#pragma simd
for(int i = 0; i < N; i++)
{
b[i] = cabs(a[i]);
}
Но самая приятная часть Cilk - это то, что вы можете просто сделать:
b[:] = cabs(a[:])
В этом случае компилятор и среда выполнения будут решать, на каком уровне он должен быть SIMDED и что должно быть парализованным (оптимальным способом является одновременное применение SIMD на крупномасштабных кусках). Поскольку это определяется планировщиком работы во время выполнения, Intel заявляет, что способна обеспечить почти оптимальное планирование и что он должен иметь возможность оптимально использовать кеш.
Использование #pragma simd (даже с -Ofast) или использование автоинъекции компиляторов - пример того, почему это плохая идея, чтобы слепо ожидать, что ваш компилятор будет эффективно реализовывать SIMD. Чтобы эффективно использовать SIMD для этого, вам нужно использовать массив структур массивов. Например, для одиночного поплавка с шириной SIMD 4 вы можете использовать
//struct of arrays of four complex numbers
struct c4 {
float x[4]; // real values of four complex numbers
float y[4]; // imaginary values of four complex numbers
};
Вот код, показывающий, как вы могли бы сделать это с помощью SSE для набора инструкций x86.
#include <stdio.h>
#include <x86intrin.h>
#define N 10
struct c4{
float x[4];
float y[4];
};
static inline void cabs_soa4(struct c4 *a, float *b) {
__m128 x4 = _mm_loadu_ps(a->x);
__m128 y4 = _mm_loadu_ps(a->y);
__m128 b4 = _mm_sqrt_ps(_mm_add_ps(_mm_mul_ps(x4,x4), _mm_mul_ps(y4,y4)));
_mm_storeu_ps(b, b4);
}
int main(void)
{
int n4 = ((N+3)&-4)/4; //choose next multiple of 4 and divide by 4
printf("%d\n", n4);
struct c4 a[n4]; //array of struct of arrays
for(int i=0; i<n4; i++) {
for(int j=0; j<4; j++) { a[i].x[j] = 1, a[i].y[j] = -1;}
}
float b[4*n4];
for(int i=0; i<n4; i++) {
cabs_soa4(&a[i], &b[4*i]);
}
for(int i = 0; i<N; i++) printf("%.2f ", b[i]); puts("");
}
Это может помочь развернуть цикл несколько раз. В любом случае, все это спорно для большого N, потому что операция связана пропускная способность памяти. Для больших N (что означает, что использование памяти намного больше, чем кеш последнего уровня), хотя #pragma omp parallel может помочь некоторым, лучшим решением является не выполнение этого для больших N. Вместо этого сделайте это в кусках, которые подходят на самом низком уровне кеш вместе с другими вычислительными операциями. Я имею в виду что-то вроде этого
for(int i = 0; i < nchunks; i++) {
for(int j = 0; j < chunk_size; j++) {
b[i*chunk_size+j] = cabs(a[i*chunk_size+j]);
}
foo(&b[i*chunck_size]); // foo is computationally intensive.
}
Я не реализовал массив структуры массива здесь, но для этого нужно легко настроить код.
Кроме того, вы можете использовать std:: future и std:: async (они являются частью С++ 11), возможно, это более четкий способ достижения того, что вы хотите сделать:
#include <future>
...
int main()
{
...
// Create async calculations
std::future<void> *futures = new std::future<void>[N];
for (int i = 0; i < N; ++i)
{
futures[i] = std::async([&a, &b, i]
{
b[i] = std::sqrt(a[i]);
});
}
// Wait for calculation of all async procedures
for (int i = 0; i < N; ++i)
{
futures[i].get();
}
...
return 0;
}
Сначала создадим асинхронные процедуры, а затем подождем, пока все не будет подсчитано.
Здесь я использую sqrt вместо кабин, потому что я просто не знаю, что такое такси. Я уверен, что это неважно.
Кроме того, возможно, вы найдете эту ссылку полезной: cplusplus.com