c++ - La forma más rápida de calcular los valores abs() de una matriz compleja
arrays complex-numbers (6)
Quiero calcular los valores absolutos de los elementos de una matriz compleja en C o C ++. La forma más fácil sería
for(int i = 0; i < N; i++)
{
b[i] = cabs(a[i]);
}
Pero para grandes vectores que serán lentos. ¿Hay una manera de acelerar eso (utilizando la paralelización, por ejemplo)? El lenguaje puede ser C o C ++.
Además, puede usar std :: future y std :: async (son parte de C ++ 11), tal vez sea una forma más clara de lograr lo que quiere hacer:
#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;
}
Primero creamos procedimientos asíncronos y luego esperamos hasta que se calcule todo.
Aquí utilizo sqrt en lugar de taxis porque no sé qué son los taxis. Estoy seguro de que no importa.
Además, quizás encuentre útil este enlace: cplusplus.com
Dado que todas las iteraciones de bucle son independientes, puede usar el siguiente código para la paralelización:
#pragma omp parallel for
for(int i = 0; i < N; i++)
{
b[i] = cabs(a[i]);
}
Por supuesto, para usar esto, debe habilitar el soporte de OpenMP mientras compila su código (generalmente utilizando el indicador / openmp o configurando las opciones del proyecto).
Puedes encontrar varios ejemplos de uso de OpenMP en wiki .
O use Concurrency :: parallele_for así:
Concurrency::parallel_for(0, N, [&a, &b](int i)
{
b[i] = cabs(a[i]);
});
Si está utilizando un compilador moderno (GCC 5, por ejemplo), puede usar Cilk+ , que le dará una buena notación de matriz, uso automático de las instrucciones SIMD y paralelización.
Entonces, si quieres ejecutarlos en paralelo harías:
#include <cilk/cilk.h>
cilk_for(int i = 0; i < N; i++)
{
b[i] = cabs(a[i]);
}
o si quieres probar SIMD:
#pragma simd
for(int i = 0; i < N; i++)
{
b[i] = cabs(a[i]);
}
Pero, la mejor parte de Cilk es que puedes hacer:
b[:] = cabs(a[:])
En este caso, el compilador y el entorno de tiempo de ejecución decidirán a qué nivel se debe SIMDed y qué se debe paralizar (la forma óptima es aplicar SIMD en grandes trozos en paralelo). Dado que esto lo decide un programador de trabajo en tiempo de ejecución, Intel afirma que es capaz de proporcionar una programación casi óptima y que debería poder hacer un uso óptimo de la memoria caché.
Usar #pragma simd (incluso con -Ofast ) o confiar en la auto-vectorización de los compiladores es un ejemplo más de por qué es una mala idea esperar ciegamente que su compilador implemente SIMD de manera eficiente. Para poder utilizar SIMD de manera eficiente para esto, necesita usar una serie de estructuras de matrices. Por ejemplo, para un solo flotador con un ancho de SIMD de 4 podría usar
//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
};
Aquí está el código que muestra cómo podría hacer esto con SSE para el conjunto de instrucciones 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("");
}
Puede ayudar a desenrollar el bucle unas cuantas veces. En cualquier caso, todo esto es discutible para una N grande porque la operación está vinculada al ancho de banda de la memoria. Para N grande (lo que significa que cuando el uso de la memoria es mucho mayor que el caché de último nivel), aunque #pragma omp parallel puede ayudar, la mejor solución es no hacer esto para N grande. En lugar de eso, hágalo en trozos que encajen en la más baja. caché de nivel junto con otras operaciones de cómputo. Me refiero a algo como esto
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.
}
No implementé una matriz de estructura de matriz aquí, pero debería ser fácil ajustar el código para eso.
Usar operaciones vectoriales.
Si tiene glibc 2.22 (bastante reciente), puede usar las capacidades SIMD de OpenMP 4.0 para operar en vectores / arrays .
Libmvec es una biblioteca matemática vectorial agregada en Glibc 2.22.
Se agregó la biblioteca matemática vectorial para admitir las construcciones SIMD de OpenMP4.0 (# 2.8 en http://www.openmp.org/mp-documents/OpenMP4.0.0.pdf ) agregando implementaciones vectoriales de funciones matemáticas vectoriales.
Las funciones matemáticas vectoriales son variantes vectoriales de las operaciones matemáticas escalares correspondientes implementadas utilizando las extensiones SIMD ISA (por ejemplo, SSE o AVX para x86_64). Toman argumentos vectoriales empaquetados, realizan la operación en cada elemento del argumento vector empaquetado y devuelven un resultado vectorial empaquetado. Usar funciones matemáticas vectoriales es más rápido que llamar repetidamente las rutinas de matemáticas escalares.
También, vea Parallel for vs omp simd: ¿cuándo usar cada uno?
Si está ejecutando Solaris, puede usar explícitamente vhypot () de la biblioteca de vectores matemáticos libmvec.so para operar en un vector de números complejos para obtener el valor absoluto de cada uno:
Descripción
Estas funciones evalúan la función hypot (x, y) para un vector completo de valores a la vez. ...
El código fuente de libmvec se puede encontrar en http://src.illumos.org/source/xref/illumos-gate/usr/src/lib/libmvec/ y el código vhypot() específicamente en http://src.illumos.org/source/xref/illumos-gate/usr/src/lib/libmvec/common/__vhypot.c No recuerdo si Sun Microsystems alguna vez proporcionó una versión Linux de libmvec.so o no.