Skip to content
Snippets Groups Projects
Commit c5db99ed authored by Jan Thorbecke's avatar Jan Thorbecke
Browse files

MKL library calls for FFT; much faster

parent fc57d903
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 333 additions and 17 deletions
FFTlib/cc1fft.c
+ 30
− 0
View file @ c5db99ed
#include "genfft.h"
#ifdef MKL
#include "mkl_dfti.h"
void dfti_status_print(MKL_LONG status);
#endif
/**
* NAME: cc1fft
......@@ -56,6 +60,10 @@ void cc1fft(complex *data, int n, int sign)
static complex *work;
REAL scl;
complex *y;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
MKL_LONG Status;
#endif
#if defined(HAVE_LIBSCS)
......@@ -90,6 +98,28 @@ void cc1fft(complex *data, int n, int sign)
nprev = n;
}
acmlcc1fft(sign, scl, inpl, n, data, 1, y, 1, work, &isys);
#elif defined(MKL)
if (n != nprev) {
DftiFreeDescriptor(&handle);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_COMPLEX, 1, (MKL_LONG)n);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiCommitDescriptor(handle);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n;
}
if (sign < 0) {
Status = DftiComputeBackward(handle, data);
}
else {
Status = DftiComputeForward(handle, data);
}
#else
cc1_fft(data, n, sign);
#endif
......
FFTlib/cc2dfft.c
+ 32
− 0
View file @ c5db99ed
#include "genfft.h"
#ifdef MKL
#include "mkl_dfti.h"
void dfti_status_print(MKL_LONG status);
#endif
/**
* NAME: cc2dfft
......@@ -60,6 +64,10 @@ void cc2dfft(complex *data, int nx, int ny, int ldx, int sign)
REAL scl;
complex *y;
complex *tmp;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nxprev=0, nyprev=0;
MKL_LONG Status, N[2];
#else
int i, j;
complex *tmp;
......@@ -110,6 +118,30 @@ void cc2dfft(complex *data, int nx, int ny, int ldx, int sign)
else {
cc1fft(data, nx, sign);
}
#elif defined(MKL)
if (nx != nxprev || ny != nyprev) {
DftiFreeDescriptor(&handle);
N[0] = nx; N[1] = ny;
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_COMPLEX, 2, N);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiCommitDescriptor(handle);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nxprev = nx;
nyprev = ny;
}
if (sign < 0) {
Status = DftiComputeBackward(handle, data);
}
else {
Status = DftiComputeForward(handle, data);
}
#else
if (ny != 1) {
ccmfft(data, nx, ny, ldx, sign);
......
FFTlib/ccmfft.c
+ 34
− 0
View file @ c5db99ed
#include "genfft.h"
#ifdef MKL
#include "mkl_dfti.h"
void dfti_status_print(MKL_LONG status);
#endif
/**
* NAME: ccmfft
......@@ -59,6 +63,10 @@ void ccmfft(complex *data, int n1, int n2, int ld1, int sign)
static complex *work;
REAL scl;
complex *y;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
MKL_LONG Status;
#endif
#if defined(HAVE_LIBSCS)
......@@ -89,6 +97,32 @@ void ccmfft(complex *data, int n1, int n2, int ld1, int sign)
nprev = n1;
}
acmlccmfft(sign, scl, inpl, n2, n1, data, 1, ld1, y, 1, ld1, work, &isys);
#elif defined(MKL)
if (n1 != nprev) {
DftiFreeDescriptor(&handle);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_COMPLEX, 1, (MKL_LONG)n1);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiCommitDescriptor(handle);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n1;
}
if (sign < 0) {
for (int j=0; j<n2; j++) {
Status = DftiComputeBackward(handle, &data[j*ld1]);
}
}
else {
for (int j=0; j<n2; j++) {
Status = DftiComputeForward(handle, &data[j*ld1]);
}
}
#else
ccm_fft(data, n1, n2, ld1, sign);
#endif
......
FFTlib/cr1fft.c
+ 49
− 4
View file @ c5db99ed
#include "genfft.h"
#ifdef MKL
#include "mkl_dfti.h"
void dfti_status_print(MKL_LONG status);
#endif
/**
* NAME: cr1fft
......@@ -52,10 +56,11 @@ void cr1fft(complex *cdata, REAL *rdata, int n, int sign)
static int isys;
static REAL *work, *table, scale=1.0;
REAL scl;
#elif defined(FFTW3)
static int nprev=0;
int iopt, ier;
static float *work;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
REAL *tmp;
MKL_LONG Status;
#endif
#if defined(HAVE_LIBSCS)
......@@ -91,6 +96,46 @@ void cr1fft(complex *cdata, REAL *rdata, int n, int sign)
nprev = n;
}
acmlcrfft(one, n, rdata, work, &isys);
#elif defined(MKL)
if (n != nprev) {
DftiFreeDescriptor(&handle);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiSetValue(handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
//This options is what we would like, but is depreciated in the future
//Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_REAL);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n;
}
tmp = (float *)malloc(n*sizeof(float));
Status = DftiComputeBackward(handle, (MKL_Complex8 *)cdata, tmp);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiComputeBackward FAIL\n");
}
rdata[0] = tmp[0];
for (int i=1; i<n; i++) {
rdata[i] = -sign*tmp[n-i];
}
free(tmp);
#else
cr1_fft(cdata, rdata, n, sign);
#endif
......
FFTlib/crmfft.c
+ 51
− 0
View file @ c5db99ed
#include "genfft.h"
#include <string.h>
#ifdef MKL
#include "mkl_dfti.h"
void dfti_status_print(MKL_LONG status);
#endif
/**
* NAME: crmfft
......@@ -66,6 +70,11 @@ void crmfft(complex *cdata, REAL *rdata, int n1, int n2, int ldc, int ldr, int s
static int isys;
static REAL *work;
REAL scl, *data;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
REAL *tmp;
MKL_LONG Status;
#endif
#if defined(HAVE_LIBSCS)
......@@ -127,6 +136,48 @@ void crmfft(complex *cdata, REAL *rdata, int n1, int n2, int ldc, int ldr, int s
for (j=0; j<n2; j++) {
memcpy(&rdata[j*ldr],&data[j*n1],n1*sizeof(REAL));
}
#elif defined(MKL)
if (n1 != nprev) {
DftiFreeDescriptor(&handle);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n1);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiSetValue(handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
//This options is what we would like, but is depreciated in the future
//Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_REAL);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n1;
}
tmp = (float *)malloc(n1*sizeof(float));
for (int j=0; j<n2; j++) {
Status = DftiComputeBackward(handle, (MKL_Complex8 *)&cdata[j*ldc], tmp);
rdata[j*ldr] = tmp[0];
for (int i=1; i<n1; i++) {
rdata[j*ldr+i] = -sign*tmp[n1-i];
}
}
free(tmp);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiComputeBackward FAIL\n");
}
#else
crm_fft(cdata, rdata, n1, n2, ldc, ldr, sign);
#endif
......
FFTlib/rc1fft.c
+ 56
− 0
View file @ c5db99ed
#include "genfft.h"
#include <string.h>
#ifdef MKL
#include "mkl_dfti.h"
void dfti_status_print(MKL_LONG status);
#endif
/**
* NAME: rc1fft
......@@ -52,6 +56,10 @@ void rc1fft(REAL *rdata, complex *cdata, int n, int sign)
static int isys;
static REAL *work, *table, scale=1.0;
REAL scl, *data;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
MKL_LONG Status;
#endif
#if defined(HAVE_LIBSCS)
......@@ -92,6 +100,42 @@ void rc1fft(REAL *rdata, complex *cdata, int n, int sign)
}
cdata[n/2].i=0.0;
free(data);
#elif defined(MKL)
if (n != nprev) {
DftiFreeDescriptor(&handle);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiSetValue(handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n;
}
Status = DftiComputeForward(handle, rdata, (MKL_Complex8 *)cdata);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiComputeForward FAIL\n");
}
for (int i=1; i<((n-1)/2)+1; i++) {
cdata[i].i *= -sign;
}
#else
rc1_fft(rdata, cdata, n, sign);
#endif
......@@ -99,6 +143,18 @@ void rc1fft(REAL *rdata, complex *cdata, int n, int sign)
return;
}
#ifdef MKL
void dfti_status_print(MKL_LONG status)
{
MKL_LONG class_error;
char* error_message;
error_message = DftiErrorMessage(status);
printf("error_message = %s \n", error_message);
return;
}
#endif
/****************** NO COMPLEX DEFINED ******************/
......
FFTlib/rcmfft.c
+ 46
− 0
View file @ c5db99ed
#include "genfft.h"
#include <string.h>
#ifdef MKL
#include "mkl_dfti.h"
void dfti_status_print(MKL_LONG status);
#endif
/**
* NAME: rcmfft
......@@ -64,6 +68,10 @@ void rcmfft(REAL *rdata, complex *cdata, int n1, int n2, int ldr, int ldc, int s
static int isys;
static REAL *work;
REAL scl, *data;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
MKL_LONG Status;
#endif
#if defined(HAVE_LIBSCS)
......@@ -120,6 +128,44 @@ void rcmfft(REAL *rdata, complex *cdata, int n1, int n2, int ldr, int ldc, int s
cdata[j*ldc+n1/2].i=0.0;
}
free(data);
#elif defined(MKL)
if (n1 != nprev) {
DftiFreeDescriptor(&handle);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n1);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiSetValue(handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n1;
}
Status = DftiComputeForward(handle, rdata, (MKL_Complex8 *)cdata);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiComputeForward FAIL\n");
}
for (int j=0; j<n2; j++) {
Status = DftiComputeForward(handle, &rdata[j*ldr], (MKL_Complex8 *)&cdata[j*ldc]);
for (int i=1; i<((n1-1)/2)+1; i++) {
cdata[j*ldc+i].i *= -sign;
}
}
#else
rcm_fft(rdata, cdata, n1, n2, ldr, ldc, sign);
#endif
......
FFTlib/test/Makefile
+ 1
− 1
View file @ c5db99ed
......@@ -2,7 +2,7 @@
include ../../Make_include
LIBS += -L$L -lgenfft -lm
LIBS += -L$L -lgenfft $(FFT) -lm
ALL: cc1test cc2dtest ccmtest rc1test rc2dtest rcmtest rc1loop cc1loop
......
FFTlib/test/cc1loop.c
+ 2
− 2
View file @ c5db99ed
#include <genfft.h>
#include <time.h>
#include <kiss_fft.h>
main () {
int main () {
int j,i,n,sign, isign;
int N, Nmax=600, Nitcc;
......@@ -65,5 +64,6 @@ main () {
N += 1;
}
return 0;
}
FFTlib/test/cc1test.c
+ 2
− 1
View file @ c5db99ed
#include <genfft.h>
#include <time.h>
main () {
int main () {
int j,i,n,sign, isign;
int N, Nmax=8192, Nitcc;
......@@ -79,5 +79,6 @@ main () {
k += 1.0;
}
return 0;
}
FFTlib/test/cc2dtest.c
+ 2
− 1
View file @ c5db99ed
#include <genfft.h>
#include <time.h>
main () {
int main () {
int j,i,n,sign, isign;
int N, Nmax=1024, Nitcc, Nlot=1024, ld1;
......@@ -142,5 +142,6 @@ main () {
Nlot *= 2;
}
return 0;
}
FFTlib/test/ccmtest.c
+ 2
− 1
View file @ c5db99ed
#include <genfft.h>
#include <time.h>
main () {
int main () {
int j,i,n,sign, isign;
int N, Nmax=4097, Nitcc, Nlot=513, ld1;
......@@ -128,5 +128,6 @@ main () {
k += 1.0;
}
return 0;
}
FFTlib/test/rc1loop.c
+ 2
− 1
View file @ c5db99ed
......@@ -5,7 +5,7 @@ void crdft(complex *cdata, REAL *rdata, int n, int sign);
void rcdft(REAL *rdata, complex *cdata, int n, int sign);
main () {
int main () {
int j,i,n,k,sign, isign;
int N, Nmax=600, Nitcc;
......@@ -94,5 +94,6 @@ main () {
N += 1;
}
return 0;
}
FFTlib/test/rc1test.c
+ 20
− 3
View file @ c5db99ed
#include <genfft.h>
#include <time.h>
main () {
int main () {
int j,i,n,sign, isign;
int N, Nmax=8192, Nitcc;
......@@ -37,12 +37,28 @@ main () {
// c_data[0] = 1.0;
t = 0.0;
/*
N=16;
scl = 1.0/(float)N;
for (j=0;j<N;j++) data[j] = c_data[j];
rc1fft(data, cdata, N, sign);
cr1fft(cdata, data, N, isign);
for (i=0; i<N; i++)
fprintf(stderr,"%s: i = %d data = %f Ref-data = %f Complex=%f,%f\n", machine, i, data[i]*scl, c_data[i],cdata[i/2].r,cdata[i/2].i);
for (j=0;j<N;j++) data[j] = c_data[j];
rc1_fft(data, cdata, N, sign);
cr1_fft(cdata, data, N, isign);
for (i=0; i<N; i++)
fprintf(stderr,"%s: i = %d data = %f Ref-data = %f Complex=%f,%f\n", machine, i, data[i]*scl, c_data[i],cdata[i/2].r,cdata[i/2].i);
return;
*/
/* FFT */
for (i=0; i<2500; i++) {
for (j=0;j<N;j++) data[j] = c_data[j];
t0 = wallclock_time();
rc1fft(data, cdata, N, sign);
cr1fft(cdata, data, N, isign);
// cr1_fft(cdata, data, N, isign);
t1 = wallclock_time();
t += t1-t0;
}
......@@ -51,9 +67,9 @@ main () {
scl = 1.0/(float)N;
for (i=0; i<N; i++) {
/*
fprintf(stderr,"%s: i = %d data = %f C-data = %f\n", machine, i, data[i].r*scl, c_data[i].r);
fprintf(stderr,"%s: i = %d data = %f C-data = %f\n", machine, i, data[i].i*scl, c_data[i].i);
fprintf(stderr,"%s: i = %d data = %f Ref-data = %f Complex=%f,%f\n", machine, i, data[i]*scl, c_data[i],cdata[i/2].r,cdata[i/2].i);
*/
if (c_data[i] != 0.0) {
diff = fabs((data[i]*scl - c_data[i]) / c_data[i]);
}
......@@ -77,5 +93,6 @@ main () {
k += 1.0;
}
return 0;
}
FFTlib/test/rc2dtest.c
+ 2
− 2
View file @ c5db99ed
#include <genfft.h>
#include <time.h>
main () {
int main () {
int l,j,i,n,sign, isign, ldr, ldc;
int N, Nmax=2048, Nlot=150, Nitcc;
......@@ -141,6 +141,6 @@ main () {
Nlot *= 2;
}
return 0;
}
FFTlib/test/rcmtest.c
+ 2
− 1
View file @ c5db99ed
#include <genfft.h>
#include <time.h>
main () {
int main () {
int l,j,i,n,sign, isign, ldr, ldc;
int N, Nmax=8193, Nlot=150, Nitcc;
......@@ -82,5 +82,6 @@ main () {
k += 1.0;
}
return 0;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment