БПФ в Matlab не позволяет выбирать, сколько потоков выполняет вычисления (http://stackoverflow.com/questions/9528833/matlabs-fftn-gets-slower-with-multithreading). По умолчанию он использует все ядра в автономном матлабе. Но в кластере каждый рабочий процесс по умолчанию запускается с одним процессором. Вы можете заставить его работать с большим количеством ядер (функция maxNumCompThreads). Это прекрасно работает с алгебраическими операциями, но функция БПФ остается (как ни странно?) одноядерной. Таким образом, я написал mex-файл, используя библиотеку fftw (как это делает Matlab), чтобы вычислить fft с желаемым количеством ядер. Но когда я пытаюсь сравнить коды с помощью планировщика FFTW_ESTIMATE (который используется по умолчанию в Matlab) и ясной мудрости, мой код остается в 3-4 раза медленнее, чем Matlab fft.
Вот код, который я использовал для мекса (применяется для 2D fft, назван FFT2mx):
#include <stdlib.h>
#include <stdio.h>
#include <mex.h>
#include <matrix.h>
#include <math.h>
#include </home/nicolas/Code/C/lib/include/fftw3.h>
void FFTNDSplit(int NumDims, const int N[], double *XReal, double *XImag, double *YReal, double *YImag, int Sign)
{
fftw_plan Plan;
fftw_iodim Dim[NumDims];
int k, NumEl;
for(k = 0, NumEl = 1; k < NumDims; k++)
{
Dim[NumDims - k - 1].n = N[k];
Dim[NumDims - k - 1].is = Dim[NumDims - k - 1].os = (k == 0) ? 1 : (N[k-1] * Dim[NumDims-k].is);
NumEl *= N[k];
}
//fftw_import_wisdom_from_filename("/home/nicolas/wisdom/wis");
if(!(Plan = fftw_plan_guru_split_dft(NumDims, Dim, 0, NULL, XReal,
XImag, YReal, YImag, FFTW_ESTIMATE)))
mexErrMsgTxt("FFTW3 failed to create plan.");
if(Sign == -1)
fftw_execute_split_dft(Plan, XReal, XImag, YReal, YImag);
else
{
fftw_execute_split_dft(Plan, XImag, XReal, YImag, YReal);
}
//if(!fftw_export_wisdom_to_filename("/home/nicolas/wisdom/wis"))
// mexErrMsgTxt("FFTW3 failed to save wisdom.");
fftw_destroy_plan(Plan);
return;
}
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[] )
{
int i, j,numCPU;
int NumDims;
const mwSize *N;
if (nrhs != 2) {
mexErrMsgIdAndTxt( "MATLAB:FFT2mx:invalidNumInputs",
"Two input argument required.");
}
if (!mxIsDouble(prhs[0])) {
mexErrMsgIdAndTxt( "MATLAB:FFT2mx:invalidNumInputs",
"Array must be double");
}
numCPU = (int) mxGetScalar(prhs[1]);
if (numCPU > 8) {
mexErrMsgIdAndTxt( "MATLAB:FFT2mx:invalidNumInputs",
"NumOfThreads < 8 requested");
}
/*if (!mxIsComplex(prhs[0])) {
mexErrMsgIdAndTxt( "MATLAB:FFT2mx:invalidNumInputs",
"Array must be complex");
}*/
NumDims = mxGetNumberOfDimensions(prhs[0]);
N = mxGetDimensions(prhs[0]);
plhs[0] = mxCreateDoubleMatrix(0, 0, mxCOMPLEX);
mxSetDimensions(plhs[0], N, NumDims);
mxSetData(plhs[0], mxMalloc( sizeof(double) * mxGetNumberOfElements(prhs[0]) ));
mxSetImagData(plhs[0], mxMalloc( sizeof(double) * mxGetNumberOfElements(prhs[0]) ));
fftw_init_threads();
fftw_plan_with_nthreads(numCPU);
FFTNDSplit(NumDims, N, (double *) mxGetPr(prhs[0]), (double *) mxGetPi(prhs[0]),
mxGetPr(plhs[0]), mxGetPi(plhs[0]), -1);
}
Связанный код Matlab:
function fft2mx(X,NumCPU)
FFT2mx(X,NumCPU)/sqrt(size(X,1)*size(X,2));
return;
Я компилирую код mex, используя статические библиотеки:
mex FFT2mx.cpp /home/nicolas/Code/C/lib/lib/libfftw3.a /home/nicolas/Code/C/lib/lib/libfftw3_threads.a
Все работает хорошо, только медленнее.
Библиотека FFTW была скомпилирована со следующими аргументами:
CC="gcc ${BUILD64} -fPIC" CXX="g++ ${BUILD64} -fPIC" \
./configure --prefix=/home/nicolas/Code/C/lib --enable-threads &&
make
make install
Я запускаю этот код на одном узле кластера с двумя четырехъядерными процессорами AMD Opteron(tm) и тестирую:
A = randn([2048 2048])+ i*randn([2048 2048]);
tic, fft2mx(A,8); toc;
tic, fftn(A); toc;
Ведьма возвращается:
Elapsed time is 0.482021 seconds.
Elapsed time is 0.151630 seconds.
Как можно настроить мой mex-код? Можно ли оптимизировать компиляцию библиотеки fftw? Есть ли способ ускорить алгоритм fftw, используя только планировщик ESTIMATE?
Я ищу любые идеи. Спасибо.
ИЗМЕНИТЬ:
Я принял во внимание то, что вы предложили (используя мудрость и статический план), и написал этот обновленный код:
# include <string.h>
# include <stdlib.h>
# include <stdio.h>
# include <mex.h>
# include <matrix.h>
# include <math.h>
# include </home/nicolas/Code/C/lib/include/fftw3.h>
char *Wisfile = NULL;
char *Wistemplate = "%s/.fftwis";
#define WISLEN 8
void set_wisfile(void)
{
char *home;
if (Wisfile) return;
home = getenv("HOME");
Wisfile = (char *)malloc(strlen(home) + WISLEN + 1);
sprintf(Wisfile, Wistemplate, home);
}
void cleanup(void) {
static fftw_plan PlanForward;
static int planlen;
static double *pr, *pi, *pr2, *pi2;
mexPrintf("MEX-file is terminating, destroying array\n");
fftw_destroy_plan(PlanForward);
fftw_free(pr2);
fftw_free(pi2);
fftw_free(pr);
fftw_free(pi);
}
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[] )
{
int i, j, numCPU, NumDims;
const mwSize *N;
fftw_complex *out, *in1;
static double *pr, *pi, *pr2, *pi2;
static int planlen = 0;
static fftw_plan PlanForward;
fftw_iodim Dim[NumDims];
int k, NumEl;
FILE *wisdom;
if (nrhs != 2) {
mexErrMsgIdAndTxt( "MATLAB:FFT2mx:invalidNumInputs",
"Two input argument required.");
}
if (!mxIsDouble(prhs[0])) {
mexErrMsgIdAndTxt( "MATLAB:FFT2mx:invalidNumInputs",
"Array must be double");
}
numCPU = (int) mxGetScalar(prhs[1]);
if (numCPU > 8) {
mexErrMsgIdAndTxt( "MATLAB:FFT2mx:invalidNumInputs",
"NumOfThreads < 8 requested");
}
if (!mxIsComplex(prhs[0])) {
mexErrMsgIdAndTxt( "MATLAB:FFT2mx:invalidNumInputs",
"Array must be complex");
}
NumDims = mxGetNumberOfDimensions(prhs[0]);
N = mxGetDimensions(prhs[0]);
for(k = 0, NumEl = 1; k < NumDims; k++)
{
Dim[NumDims - k - 1].n = N[k];
Dim[NumDims - k - 1].is = Dim[NumDims - k - 1].os = (k == 0) ? 1 : (N[k-1] * Dim[NumDims-k].is);
NumEl *= N[k];
}
/* If different size, free/destroy */
if(N[0] != planlen && planlen > 0) {
fftw_free(pr2);
fftw_free(pi2);
fftw_free(pr);
fftw_free(pi);
fftw_destroy_plan(PlanForward);
planlen = 0;
}
mexAtExit(cleanup);
/* Init */
fftw_init_threads();
// APPROACH 1
//pr = (double *) mxGetPr(prhs[0]);
//pi = (double *) mxGetPi(prhs[0]);
// APPROACH 2
pr = (double *) fftw_malloc( sizeof(double) * mxGetNumberOfElements(prhs[0]) );
pi = (double *) fftw_malloc( sizeof(double) * mxGetNumberOfElements(prhs[0]) );
tmp1 = (double *) mxGetPr(prhs[0]);
tmp2 = (double *) mxGetPi(prhs[0]);
for(k=0;k<mxGetNumberOfElements(prhs[0]);k++)
{
pr[k] = tmp1[k];
pi[k] = tmp2[k];
}
plhs[0] = mxCreateNumericMatrix(0, 0, mxDOUBLE_CLASS, mxCOMPLEX);
mxSetDimensions(plhs[0], N, NumDims);
mxSetData(plhs[0], (double* ) fftw_malloc( sizeof(double) * mxGetNumberOfElements(prhs[0]) ));
mxSetImagData(plhs[0], (double* ) fftw_malloc( sizeof(double) * mxGetNumberOfElements(prhs[0]) ));
pr2 = mxGetPr(plhs[0]);
pi2 = mxGetPi(plhs[0]);
fftw_init_threads();
fftw_plan_with_nthreads(numCPU);
/* Get any accumulated wisdom. */
set_wisfile();
wisdom = fopen(Wisfile, "r");
if (wisdom) {
fftw_import_wisdom_from_file(wisdom);
fclose(wisdom);
}
/* Compute plan */
//printf("%d",planlen);
if(planlen == 0 ) {
fftw_plan_with_nthreads(numCPU);
PlanForward = fftw_plan_guru_split_dft(NumDims, Dim, 0, NULL, pr, pi, pr2, pi2, FFTW_MEASURE);
planlen = N[0];
}
/* Save the wisdom. */
wisdom = fopen(Wisfile, "w");
if (wisdom) {
fftw_export_wisdom_to_file(wisdom);
fclose(wisdom);
}
/* execute */
fftw_execute_split_dft(PlanForward, pr, pi, pr2, pi2);
fftw_cleanup_threads();
}
Теперь я сталкиваюсь с некоторыми ошибками сегментации после нескольких вызовов (от 2 до 6) функции, и я не могу понять, почему. Я пробовал другой способ инициализации по указателю. Я также где-то читал, что указатель плана должен быть статическим, чтобы работать с соответствующим статическим планом. Что-нибудь вы видите, что я делаю неправильно?
Еще раз спасибо за ваши идеи.
