He publicado esto en matlab central, pero no recibo ninguna de las respuestas así que pensé en volver a publicar aquí.

Hace poco escribí una simple rutina en Matlab que utiliza una FFT en un bucle; la FFT domina los cálculos. Escribí la misma rutina en mex sólo para la experimentación y llama a la FFTW 3.3 de la biblioteca. Resulta que la rutina de matlab ejecuta más rápido que el mex rutina para matrices muy grandes (alrededor de dos veces más rápido). La mex usos de rutina de la sabiduría y realiza la misma FFT cálculos. También sé que matlab utiliza FFTW, pero es posible que su versión es ligeramente más optimizada? Que incluso se utiliza el FFTW_EXHAUSTIVE bandera y su todavía alrededor de dos veces más lento para grandes conjuntos de MATLAB contraparte. Además, me aseguró que el matlab que se utiliza un único hilo con el «-singleCompThread» la bandera y el archivo mex he utilizado no fue en modo de depuración. Solo por curiosidad, si este era el caso – o si hay algunas optimizaciones de matlab se usa bajo el capó que no sé acerca de. Gracias.

Aquí está el mex porción:

void class_cg_toeplitz::analysis() {
//This method computes CG iterations using FFTs
//Check for wisdom
if(fftw_import_wisdom_from_filename("cd.wis") == 0) {
mexPrintf("wisdom not loaded.\n");
} else {
mexPrintf("wisdom loaded.\n");
}
//Set FFTW Plan - use interleaved FFTW
fftw_plan plan_forward_d_buffer;    
fftw_plan plan_forward_A_vec;       
fftw_plan plan_backward_Ad_buffer;
fftw_complex *A_vec_fft;
fftw_complex *d_buffer_fft;
A_vec_fft = fftw_alloc_complex(n);
d_buffer_fft = fftw_alloc_complex(n);
//CREATE MASTER PLAN - Do this on an empty vector as creating a plane 
//with FFTW_MEASURE will erase the contents; 
//Use d_buffer
//This is somewhat dangerous because Ad_buffer is a vector; but it does not
//get resized so &Ad_buffer[0] should work
plan_forward_d_buffer = fftw_plan_dft_r2c_1d(d_buffer.size(),&d_buffer[0],d_buffer_fft,FFTW_EXHAUSTIVE);
plan_forward_A_vec = fftw_plan_dft_r2c_1d(A_vec.height,A_vec.value,A_vec_fft,FFTW_WISDOM_ONLY);
//A_vec_fft.*d_buffer_fft will overwrite d_buffer_fft
plan_backward_Ad_buffer = fftw_plan_dft_c2r_1d(Ad_buffer.size(),d_buffer_fft,&Ad_buffer[0],FFTW_EXHAUSTIVE);
//Get A_vec_fft
fftw_execute(plan_forward_A_vec);
//Find initial direction - this is the initial residual
for (int i=0;i<n;i++) {
d_buffer[i] = b.value[i];
r_buffer[i] = b.value[i];
}    
//Start CG iterations
norm_ro = norm(r_buffer);
double fft_reduction = (double)Ad_buffer.size(); //Must divide by size of vector because inverse FFT does not do this
while (norm(r_buffer)/norm_ro > relativeresidual_cutoff) {        
//Find Ad - use fft
fftw_execute(plan_forward_d_buffer);    
//Get A_vec_fft.*fft(d) - A_vec_fft is only real, but d_buffer_fft
//has complex elements; Overwrite d_buffer_fft        
for (int i=0;i<n;i++) {
d_buffer_fft[i][0] = d_buffer_fft[i][0]*A_vec_fft[i][0]/fft_reduction;
d_buffer_fft[i][1] = d_buffer_fft[i][1]*A_vec_fft[i][0]/fft_reduction;
}        
fftw_execute(plan_backward_Ad_buffer); 
//Calculate r'*r
rtr_buffer = 0;
for (int i=0;i<n;i++) {
rtr_buffer = rtr_buffer + r_buffer[i]*r_buffer[i];
}    
//Calculate alpha
alpha = 0;
for (int i=0;i<n;i++) {
alpha = alpha + d_buffer[i]*Ad_buffer[i];
}    
alpha = rtr_buffer/alpha;
//Calculate new x
for (int i=0;i<n;i++) {
x[i] = x[i] + alpha*d_buffer[i];
}   
//Calculate new residual
for (int i=0;i<n;i++) {
r_buffer[i] = r_buffer[i] - alpha*Ad_buffer[i];
}   
//Calculate beta
beta = 0;
for (int i=0;i<n;i++) {
beta = beta + r_buffer[i]*r_buffer[i];
}  
beta = beta/rtr_buffer;
//Calculate new direction vector
for (int i=0;i<n;i++) {
d_buffer[i] = r_buffer[i] + beta*d_buffer[i];
}  
*total_counter = *total_counter+1;
if(*total_counter >= iteration_cutoff) {
//Set total_counter to -1, this indicates failure
*total_counter = -1;
break;
}
}
//Store Wisdom
fftw_export_wisdom_to_filename("cd.wis");
//Free fft alloc'd memory and plans
fftw_destroy_plan(plan_forward_d_buffer);
fftw_destroy_plan(plan_forward_A_vec);
fftw_destroy_plan(plan_backward_Ad_buffer);
fftw_free(A_vec_fft);
fftw_free(d_buffer_fft);
};

Aquí está el matlab porción:

% Take FFT of A_vec.
A_vec_fft = fft(A_vec); % Take fft once
% Find initial direction - this is the initial residual 
x = zeros(n,1); % search direction
r = zeros(n,1); % residual
d = zeros(n+(n-2),1); % search direction; pad to allow FFT
for i = 1:n
d(i) = b(i); 
r(i) = b(i); 
end
% Enter CG iterations
total_counter = 0;
rtr_buffer = 0;
alpha = 0;
beta = 0;
Ad_buffer = zeros(n+(n-2),1); % This holds the product of A*d - calculate this once per iteration and using FFT; only 1:n is used
norm_ro = norm(r);
while(norm(r)/norm_ro > 10^-6)
% Find Ad - use fft
Ad_buffer = ifft(A_vec_fft.*fft(d)); 
% Calculate rtr_buffer
rtr_buffer = r'*r;
% Calculate alpha    
alpha = rtr_buffer/(d(1:n)'*Ad_buffer(1:n));
% Calculate new x
x = x + alpha*d(1:n);
% Calculate new residual
r = r - alpha*Ad_buffer(1:n);
% Calculate beta
beta = r'*r/(rtr_buffer);
% Calculate new direction vector
d(1:n) = r + beta*d(1:n);      
% Update counter
total_counter = total_counter+1; 
end

En términos de tiempo, para N = 50000 y b = 1:n, se tarda alrededor de 10,5 segundos con mex y 4.4 segundos con matlab. Estoy usando R2011b. Gracias

  • ¿Cuáles son las dimensiones de sus datos, y cuáles son los absolutos de los tiempos?
  • Están ambos en el lugar fft?
  • usted podría ejecutar su código de Matlab con el analizador de encendido para obtener información más detallada sobre el tiempo dedicado a cada función (en por ciento), esto podría dar una pista fueron Matlab está optimizado
  • Me encontré con el profiler en el matlab parte; casi todos se gasta en la FFT. También corrí valgrind en la mex y básicamente todo lo que se gasta en la FFT así. De la 4.4 segundos para que el matlab parte, el analizador dice 4 segundos se gastan en la FFT en matlab. Para el mex parte, valgrind dice 84.99% se gasta en fftw_execute.
  • interesante. Acabo de comprobar: matlab también tiene un fftw comando que permite controlar los parámetros de optimización utilizado internamente para la fftw lib(->ayuda fftw). con este comando también se puede obtener la sabiduría de la base de datos de matlab se ha estado utilizando para los cálculos. sería interesante ver lo que los resultados que usted consigue cuando usted alimenta a matlabs la sabiduría de la base de datos para su programa de c++ y viceversa…
  • Voy a mirar en esto.
  • También he descubierto que matlab utiliza un no determinista número de iteraciones para converger. Por la n – 1000, se tarda entre 83-85 iteraciones… @daño también, mi matlab fftw versión y la versión que he instalado son diferentes (3.3.3 vs 3.2.2) por lo que la sabiduría no son compatibles por lo que he probado.
  • En Matlab bin/<PLATAFORMA> usted puede encontrar el archivo ‘fftw.spec», que especifica las diferentes bibliotecas de diferentes CPU – así que yo diría que las bibliotecas están especialmente optimizadas.
  • ¿Estás usando el mismo residual criterios para ambos? Veo 10e-6 en la secuencia de comandos de matlab, pero sólo relativeresidual_cutoff en la mex, cuya definición no se muestra.
  • Creo, que necesita para construir su fftw basado en multi-hilo LAPACK.
  • He comprobado bin/<plataforma> y no todos los archivos con nombre fftw.spec. @ CaptainMurphy Tanto de corte criterio son 10e-6. @ iampat yo quería una comparación directa con un único hilo de rendimiento. Matlab todavía supera a mex incluso con un único hilo de rendimiento.
  • Este puede ser trivial e irrelevante, pero yo aviso 2 llamadas a fftw_execute() en el .mex código; pero sólo 1 en el matlab. Supongo que hay algo obvio que me estoy perdiendo, pero he pensado que me gustaría comentar.
  • En matlab: Ad_buffer = ifft(A_vec_fft.*fft(d)); Esta línea es de dos llamadas. La primera es un avance de la fft y la segunda es la inversa. Thats whats que se realiza en la mex versión.
  • hacer los dos conjuntos de datos tienen una longitud igual a la potencia de 2? debido a que la velocidad de cálculo de la FFT.
  • usted puede estar interesado en saber que fft tiene más mejoras de rendimiento en la última MATLAB versión R2013a (en Cpu que soporte conjunto de instrucciones AVX): mathworks.com/help/matlab/release-notes.html#btsiwqu-1
  • Muy interesante. Estoy leyendo sobre el montaje ahora así que espero ser capaz de entender el AVX cosas en el futuro cercano. Tal vez esto también sugiere matlab implementa mejoras a sí mismos fuera de fftw.

InformationsquelleAutor Justin | 2013-03-08

4 Comentarios

  1. 13

    Un par de observaciones, más que una respuesta definitiva, pues yo no conozco a ninguno de los detalles de MATLAB FFT aplicación:

    • Basado en el código, puedo ver dos explicaciones para la diferencia de velocidad:
      • la diferencia de velocidad es explicado por las diferencias en los niveles de optimización de la FFT
      • el bucle while en MATLAB se ejecuta un número significativamente menor de los tiempos de

    Voy a suponer que usted ya se veía en el segundo problema, y que el número de iteraciones son comparables. (Si no, esto es más probable que para algunos la exactitud y la pena más investigaciones.)

    Ahora, con respecto a la FFT comparación de velocidad:

    • Sí, la teoría es que FFTW es más rápido que otros de alto nivel de la FFT de las implementaciones, pero solo es aplicable siempre que comparar manzanas con manzanas: aquí está la comparación de las implementaciones en un nivel más abajo, en el ámbito de la asamblea, donde no sólo la selección del algoritmo, pero su optimización real para un procesador específico y por los desarrolladores de software con diferentes habilidades que viene en el juego
    • He optimizado o revisado optimizado Fft en asamblea en muchos procesadores a lo largo del año (yo estaba en el benchmarking de la industria) y gran algoritmos son sólo una parte de la historia. Hay consideraciones que son muy específicos de la arquitectura de la codificación para (contabilidad para las latencias, la programación de las instrucciones, la optimización de registro de uso, disposición de los datos en la memoria, la contabilidad de la sucursal de la toma/no se toman las latencias, etc.) y que hacen que las diferencias tan importante como la selección del algoritmo.
    • Con N=500000, estamos hablando también de gran tamaño de los búferes de memoria: la otra puerta para más optimizaciones que rápidamente puede llegar a ser muy específicos a la plataforma de ejecutar su código en: cómo te las arreglas para evitar errores de caché de no ser dictada por el algoritmo, por lo tanto como la forma en que el flujo de datos y lo optimizaciones de un desarrollador de software puede tener utilizan para llevar los datos en la memoria de manera eficaz.
    • Aunque no sé los detalles de MATLAB implementación de la FFT, estoy bastante seguro de que un ejército de DSP ingenieros ha sido (y todavía es) rectificado en su optimización, ya que es la clave para muchos de los diseños. Esto podría muy bien significar que MATLAB tenía la combinación correcta de los desarrolladores para producir un mucho más rápido FFT.
    • Lolo, el quid de la cuestión es que MATLAB implementa FFTW. Curiosamente el número de iteraciones hasta que la convergencia en la rutina de MATLAB parece ser no determinista. para N = 1000, toma 83-85 mientras que el mex versión es constante (85 para N = 1000 si recuerdo correctamente). En este punto, he de ordenación de la que acaba de concluir matlab debe de estar haciendo algo «bajo el capó», que no estoy al tanto acerca de… o eso, o mi mex aplicación es más lento debido a que me perdí de una optimización en algún lugar. No estoy seguro.
    • Todo lo que usted dice de los puntos a la misma conclusión que yo: 83-85 vs 85 significa que la FFT, el rendimiento se explica la diferencia, y lo hace el 90% vs 84.99% de datos de perfiles. La implementación en MATLAB es simplemente mejor optimizado, que es plausible con un algoritmo como el que, con tantas oportunidades de optimización en cada etapa. Yo no las calificaría como «bajo el capó» trucos pero solo un tiempo bien invertido en la creación de una FFT aplicación que está optimizado en un mejor nivel que el de MEX contraparte uso. No creo que le falta algo en su mex aplicación.
    • Como Lolo es decir, nada de lo que está sucediendo «bajo el capó», Matlab tiene una mejor optimizado la aplicación de la MKL, revise mi respuesta, es la respuesta a su pregunta…
    • Está usted seguro de que MATLAB utiliza MKL para su fft (si este es el caso, voy a aceptar su respuesta)? Pensé que MATLAB utiliza fftw. Ellos no afirman explícitamente que utilizan en su documentación para fft pero que tienen citas para FFTW.org y a la FFTW de papel. También tienen una función de «fftw» que le permite meterse con sabiduría. Sin embargo, es posible que el uso de MKL para su fft tal vez para las personas con procesadores intel (estoy usando un i7). Sé FFTW utiliza codelets para diferentes tamaños de fft. Tal vez es posible que de matlab codelets son más optimizado que los proporcionados por fftw.org.
    • La cosa se complica: software.intel.com/en-us/articles/….
    • Sí, estoy bastante seguro de que, verificará el lunes para estar absolutamente seguro (sólo tengo Matlab en mi oficina)… Y el artículo que usted proporcione de Intel acaba de estados que Intel MKL apoya la FFTW interfaces, pero la implementación subyacente es específico de Intel, y saben cómo sus propios procesadores funcionan muy bien, por lo que optimizar de manera muy eficiente. Mejor que FFTW los desarrolladores pueden. Realmente creo que explica la diferencia de rendimiento.
    • Sólo es para su información, me han concedido la plena recompensa, pero esta respuesta es mera especulación…
    • He hecho más investigación y parece que mi respuesta es completamente equivocado… 🙂 Mis habilidades de investigación son por desgracia muy limitada para entender todo con 100% de certeza, pero parece que Matlab tiene su propia implementación de fftw en: libmwfftw.dll, que no parecen depender de libfftw3.dll (que está lleno de Matlab de todos modos…) Esto puede parecer extraño, pero como lo que yo puedo decir, no veo la dependencia a libfftw3.dll… Así que Lolo podría estar en lo correcto, o daños a la derecha, y es el Matlab sabiduría que acelera el código. Has probado a desactivar el Matlab sabiduría con fftw('wisdom', [])?

  2. 8

    Esto es clásico de la ganancia de rendimiento, gracias a su bajo nivel y la arquitectura específica de optimización.

    Matlab utiliza la FFT de la Intel MKL (Math Kernel Library) binario (mkl.dll). Estas son rutinas optimizadas (en el nivel de conjunto) por Intel para procesadores Intel. Incluso en los AMD parece dar buen rendimiento aumenta.

    FFTW parece normal biblioteca de c que no es tan optimizado. Por lo tanto la ganancia de rendimiento para el uso de la MKL.

    • MATLAB incluye su propia construcción de la fuente abierta FFTW biblioteca, compilado con soporte multihilo y ESS/AVX vectorizados instrucciones. Llamar version('-fftw') muestra FFTW-3.3.3-sse2-avx. Hay dos bibliotecas compartidas se encuentran en MATLAB bin que exportar el FFTW API de interfaz: libmwfftw3.dll y libmwfftw3f.dll (además de un tercer lib libmwmfl_fft.dll construido en la parte superior de las dos anteriores se pretende resumen el uso de FFTW planes). Así que a pesar de MATLAB utiliza Intel MKL como la optimización de BLAS/LAPACK aplicación, no está llamando a la FFTW interfaz de MKL como lo que puedo decir.
    • Gracias por la aclaración! BTW, a saber, ¿cómo te diste cuenta de que estos dos archivos binarios de exportación de la FFTW interfaz API? ¿Y sabes cuál es la diferencia entre ambos binarios? En mi R2010a sólo tengo uno libmwfftw.dll la biblioteca de…
    • Yo simplemente uso Dependency Walker para obtener una lista de funciones exportadas por cualquier DLL (verás funciones familiares como fftw_plan_dft_1d, fftw_execute, etc..). La primera DLL corresponde a la doblela precisión de versión de FFTW, la segunda es de la sola precisión de la versión (tengo la última MATLAB R2014a). Se me olvidó decir que también hay otros dos archivos DLL de la aplicación de las paralelas de memoria distribuida de la versión de FFTW que utiliza MPI (buscar libmwfftw3_mpi.dll y libmwfftw3f_mpi.dll)
    • también si usted tiene la PCT caja de herramientas, fft se pueden ejecutar en la GPU, que se implementa el uso de la cuFFT library (busque el cufft*.dll archivo)
    • Thx por la info
  3. 3

    He encontrado el siguiente comentario en el sitio web de the MathWorks [1]:

    Nota en las grandes potencias de 2: Para FFT dimensiones que son potencias de
    2, entre 2^14 y 2^22, el software MATLAB usos especiales precargado
    la información en su base de datos interna para optimizar el cálculo de la FFT.
    Ningún ajuste se realiza cuando la dimensión de la ITF es una potencia de 2,
    a menos que desactive la base de datos utilizando el comando fftw(‘sabiduría’, []).

    Aunque se refiera a potencias de 2, es posible sugerencia sobre que MATLAB emplea su propia «sabiduría» cuando el uso de FFTW para cierta (gran) variedad de tamaños. Considere la posibilidad de: 2^16 = 65536.

    [1] R2013b Documentación disponible a través de http://www.mathworks.de/de/help/matlab/ref/fftw.html (consultado el 29 de Octubre de 2013)

  4. 3

    EDICIÓN: @wakjah ‘s respuesta a esta respuesta es precisa: FFTW no admite división real e imaginaria de la memoria de almacenamiento a través de su Gurú de la interfaz. Mi afirmación acerca de la piratería es, por tanto, no precisa, pero puede muy bien aplicar si FFTW del Gurú de la interfaz no se utiliza – que es el caso por defecto, así que ten cuidado todavía!

    Primero, lo siento por ser un año de retraso. No estoy convencido de que el aumento de velocidad que vemos proviene de MKL o de otras optimizaciones. Hay algo fundamentalmente diferente entre FFTW y Matlab, y eso es lo complejo que los datos se almacenan en la memoria.

    En Matlab, las partes real e imaginaria de un complejo vector X son independientes de las matrices de Xre[i] y Xim[i] (lineal en la memoria, eficiente cuando se trabaja sobre cualquiera de ellos por separado).

    En FFTW, las partes real e imaginaria son entrelazados como double[2] por defecto, es decir, X[i][0] es la parte real, y X[i][1] es la parte imaginaria.

    Por lo tanto, el uso de la FFTW biblioteca en mex archivos uno no puede usar el Matlab array directamente, sino que debe asignar memoria nueva en primer lugar, a continuación, paquete de la entrada de Matlab en FFTW formato, y luego descomprimir la salida de FFTW en Matlab formato. es decir,

    X = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    Y = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

    luego

    for (size_t i=0; i<N; ++i) {
    X[i][0] = Xre[i];
    X[i][1] = Xim[i];
    }

    luego

    for (size_t i=0; i<N; ++i) {
    Yre[i] = Y[i][0];
    Yim[i] = Y[i][1];
    }

    Por lo tanto, esto requiere 2x asignaciones de memoria + 4x lee + 4x escribe — todos de tamaño N. Esto toma un peaje a la velocidad del sabio en grandes problemas.

    Tengo una corazonada de que Mathworks puede tener hackeado el FFTW3 código que le permita leer vectores de entrada directamente en el formato de Matlab, que evita todas las anteriores.

    En este escenario, sólo se puede asignar X y el uso de X de Y para ejecutar FFTW en el lugar (como fftw_plan_*(N, X, X, ...) en lugar de fftw_plan_*(N, X, Y, ...)), ya que va a ser copiado a la Yre y Yim Matlab vector, a menos que la aplicación requiere que el/los beneficios de mantener a X e y por separado.

    EDITAR: Busca en el consumo de memoria en tiempo real cuando se ejecuta Matlab del fft2() y mi código basado en la fftw3 biblioteca, muestra que Matlab sólo se asigna sólo una compleja matriz (la salida), mientras que la de mi código necesidades de estas dos matrices (el *fftw_complex búfer además de la salida de Matlab). En lugar de la conversión entre Matlab y fftw formatos no es posible porque el Matlab es real e imaginaria de las matrices no son consecutivos en la memoria. Esto sugiere que Mathworks hackeado la fftw3 biblioteca para leer/escribir los datos utilizando el formato de Matlab.

    Otro optimización de varias llamadas, es asignar de forma persistente (utilizando mexMakeMemoryPersistent()). No estoy seguro de si la implementación en Matlab hace tan bien.

    Saludos.

    p.s. Como una nota del lado, el Matlab complejo formato de almacenamiento de datos es más eficiente para operar en el real o imaginaria de los vectores por separado. En FFTW del formato que tendría que hacer ++2 lecturas de memoria.

    • Salvo que el FFTW Gurú de la Interfaz soporta dividir reales y complejos, matrices – es decir, el mismo que el formato de MATLAB – hacking necesario.
    • Acepto la corrección, +1 y gracias! He editado mi respuesta para reflejar su respuesta.

Dejar respuesta

Please enter your comment!
Please enter your name here