1
Matlab의 FFT는 계산을 수행하는 스레드 수를 선택할 수 없습니다 (http://stackoverflow.com/questions/9528833/matlabs-fftn-gets-slower-with-multithreading).). 기본적으로 독립 실행 형 matlab에있는 모든 코어를 사용합니다. 그러나 클러스터에서는 각 작업자가 기본적으로 단일 CPU로 시작됩니다. 당신은 더 많은 코어 (maxNumCompThreads 기능)와 함께 작동하도록 할 수 있습니다. 이것은 대수 연산과 완벽하게 작동하지만 FFT 함수는 (이상하게도?) 단일 코어로 남아 있습니다. 필자는 fftw 라이브러리 (matlab처럼)를 사용하여 원하는 코어 수와 함께 fft를 계산하는 mex 파일을 작성했습니다. 그러나 FFTW_ESTIMATE 플래너 (Matlab의 기본 설정)와 명확한 지혜를 사용하여 코드를 비교하려고하면 내 코드가 Matlab fft보다 3-4 배 느려집니다. 여기 FFTW를 Matlab FFT로 최적화
내가 MEX 사용 코드 (2 차원 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 코드를 컴파일 라이브러리 :
모든 것이 잘 작동하며 느리게 진행됩니다.
FFTW 라이브러리는 다음 인수로 컴파일하고있다 : 나는 2 쿼드 코어 AMD 옵테론 하나 개의 클러스터 노드에서이 코드를 실행하고 내가 함께 테스트
CC="gcc ${BUILD64} -fPIC" CXX="g++ ${BUILD64} -fPIC" \
./configure --prefix=/home/nicolas/Code/C/lib --enable-threads &&
make
make install
:
A = randn([2048 2048])+ i*randn([2048 2048]);
tic, fft2mx(A,8); toc;
tic, fftn(A); toc;
마녀 반환 : 내 MEX 코드가
Elapsed time is 0.482021 seconds.
Elapsed time is 0.151630 seconds.
을 튜닝 할 수 있습니까? fftw 라이브러리의 컴파일을 최적화 할 수 있습니까? ESTIMATE 계획 수립 자만 사용하여 fftw 알고리즘의 속도를 높일 수 있습니까?
나는 통찰력을 찾고있다. 고맙습니다.
편집 : 나는 지금 한 후 일부 세그멘테이션 오류가 발생하고
# 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 사이)하고 이유를 알 수 없습니다. 나는 포인터로 초기화하는 다른 방법을 시도했다. 필자는 어딘가에서 계획의 포인터가 해당 정적 계획을 사용하기 위해 정적이어야한다고 읽었습니다. 내가 뭘 잘못하고있는 걸 본거야?
귀하의 통찰력을 다시 한번 감사드립니다.
의견을 보내 주셔서 감사합니다. 함수에 대한 단일 호출로 코드를 테스트하기 때문에 계획 생성이 패널티가되어서는 안됩니다. 맞습니까? – Nicolas
계획 생성 *은 페널티입니다. 계획을 작성 (및 파기)하는 데는 시간이 걸립니다. MEX 내부의 루프에서 1000 개의 FFT를 수행하고 이것을 1000 개의 MATLAB FFT와 비교하십시오. 그러면 큰 차이가 있습니다. –
여러 FFT를 적용 할 때 지혜를 저장하고 다시로드 한 후 새 계획을 만드는 것은 암기에 의한 것입니다. 나는 계획 자체를 다시 matlab에 전달한 다음 mex 함수로 다시 전달하는 방법을 모르겠다. 하지만 여전히 단일 _FFT를 적용 할 때 코드가 maltab보다 훨씬 느린 이유를 이해하지 못합니다 (이는 계획도 계산해야 함). 그런데 FFT (arround fftw_execute)를 100 번 반복하면 Matlab은 11 초, FFTW 코드는 13 초, MATLAB은 56 초, FFTW는 76 초, 500 회 반복됩니다 ... 아직 뭔가가 없습니다. – Nicolas