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

thread-safe MKL FFT's

parent b8cbc2f4
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
.gitignore
+ 1
− 1
View file @ f8db7205
......@@ -4,7 +4,7 @@ bin/*
*.su
*.bin
*/.DS*
*/*/.DS*
**/.DS*
.DS*
Make_include
fdelmodc_cuda.tgz
......
FFTlib/cc1fft.c
+ 16
− 9
View file @ f8db7205
......@@ -61,9 +61,16 @@ void cc1fft(complex *data, int n, int sign)
REAL scl;
complex *y;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
static DFTI_DESCRIPTOR_HANDLE handle[MAX_NUMTHREADS];
static int nprev[MAX_NUMTHREADS];
MKL_LONG Status;
int id;
#endif
#ifdef _OPENMP
id = omp_get_thread_num();
#else
id = 0;
#endif
#if defined(HAVE_LIBSCS)
......@@ -99,26 +106,26 @@ void cc1fft(complex *data, int n, int sign)
}
acmlcc1fft(sign, scl, inpl, n, data, 1, y, 1, work, &isys);
#elif defined(MKL)
if (n != nprev) {
DftiFreeDescriptor(&handle);
if (n != nprev[id]) {
DftiFreeDescriptor(&handle[id]);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_COMPLEX, 1, (MKL_LONG)n);
Status = DftiCreateDescriptor(&handle[id], 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);
Status = DftiCommitDescriptor(handle[id]);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n;
nprev[id] = n;
}
if (sign < 0) {
Status = DftiComputeBackward(handle, data);
Status = DftiComputeBackward(handle[id], data);
}
else {
Status = DftiComputeForward(handle, data);
Status = DftiComputeForward(handle[id], data);
}
#else
cc1_fft(data, n, sign);
......
FFTlib/ccmfft.c
+ 16
− 10
View file @ f8db7205
......@@ -64,10 +64,16 @@ void ccmfft(complex *data, int n1, int n2, int ld1, int sign)
REAL scl;
complex *y;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
static DFTI_DESCRIPTOR_HANDLE handle[MAX_NUMTHREADS];
static int nprev[MAX_NUMTHREADS];
MKL_LONG Status;
int j;
int j, id;
#endif
#ifdef _OPENMP
id = omp_get_thread_num();
#else
id = 0;
#endif
#if defined(HAVE_LIBSCS)
......@@ -99,29 +105,29 @@ void ccmfft(complex *data, int n1, int n2, int ld1, int sign)
}
acmlccmfft(sign, scl, inpl, n2, n1, data, 1, ld1, y, 1, ld1, work, &isys);
#elif defined(MKL)
if (n1 != nprev) {
DftiFreeDescriptor(&handle);
if (n1 != nprev[id]) {
DftiFreeDescriptor(&handle[id]);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_COMPLEX, 1, (MKL_LONG)n1);
Status = DftiCreateDescriptor(&handle[id], 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);
Status = DftiCommitDescriptor(handle[id]);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n1;
nprev[id] = n1;
}
if (sign < 0) {
for (j=0; j<n2; j++) {
Status = DftiComputeBackward(handle, &data[j*ld1]);
Status = DftiComputeBackward(handle[id], &data[j*ld1]);
}
}
else {
for (j=0; j<n2; j++) {
Status = DftiComputeForward(handle, &data[j*ld1]);
Status = DftiComputeForward(handle[id], &data[j*ld1]);
}
}
#else
......
FFTlib/cr1fft.c
+ 18
− 11
View file @ f8db7205
......@@ -57,13 +57,20 @@ void cr1fft(complex *cdata, REAL *rdata, int n, int sign)
static REAL *work, *table, scale=1.0;
REAL scl;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
static DFTI_DESCRIPTOR_HANDLE handle[MAX_NUMTHREADS];
static int nprev[MAX_NUMTHREADS];
REAL *tmp;
MKL_LONG Status;
int i;
int i, id;
#endif
#ifdef _OPENMP
id = omp_get_thread_num();
#else
id = 0;
#endif
#if defined(HAVE_LIBSCS)
if (n != nprev) {
isys = 0;
......@@ -98,34 +105,34 @@ void cr1fft(complex *cdata, REAL *rdata, int n, int sign)
}
acmlcrfft(one, n, rdata, work, &isys);
#elif defined(MKL)
if (n != nprev) {
DftiFreeDescriptor(&handle);
if (n != nprev[id]) {
DftiFreeDescriptor(&handle[id]);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n);
Status = DftiCreateDescriptor(&handle[id], 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);
Status = DftiSetValue(handle[id], 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);
Status = DftiSetValue(handle[id], DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
Status = DftiCommitDescriptor(handle[id]);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n;
nprev[id] = n;
}
tmp = (float *)malloc(n*sizeof(float));
Status = DftiComputeBackward(handle, (MKL_Complex8 *)cdata, tmp);
Status = DftiComputeBackward(handle[id], (MKL_Complex8 *)cdata, tmp);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiComputeBackward FAIL\n");
......
FFTlib/crmfft.c
+ 19
− 12
View file @ f8db7205
......@@ -71,13 +71,20 @@ void crmfft(complex *cdata, REAL *rdata, int n1, int n2, int ldc, int ldr, int s
static REAL *work;
REAL scl, *data;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
static DFTI_DESCRIPTOR_HANDLE handle[MAX_NUMTHREADS];
static int nprev[MAX_NUMTHREADS];
REAL *tmp;
MKL_LONG Status;
int i, j;
int i, j, id;
#endif
#ifdef _OPENMP
id = omp_get_thread_num();
#else
id = 0;
#endif
#if defined(HAVE_LIBSCS)
nmp = mp_my_threadnum();
if(nmp>=MAX_NUMTHREADS) {
......@@ -138,37 +145,37 @@ void crmfft(complex *cdata, REAL *rdata, int n1, int n2, int ldc, int ldr, int s
memcpy(&rdata[j*ldr],&data[j*n1],n1*sizeof(REAL));
}
#elif defined(MKL)
if (n1 != nprev) {
DftiFreeDescriptor(&handle);
if (n1 != nprev[id]) {
DftiFreeDescriptor(&handle[id]);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n1);
Status = DftiCreateDescriptor(&handle[id], 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);
Status = DftiSetValue(handle[id], 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);
Status = DftiSetValue(handle[id], 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);
//Status = DftiSetValue(handle[id], DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_REAL);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
Status = DftiCommitDescriptor(handle[id]);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n1;
nprev[id] = n1;
}
tmp = (float *)malloc(n1*sizeof(float));
for (j=0; j<n2; j++) {
Status = DftiComputeBackward(handle, (MKL_Complex8 *)&cdata[j*ldc], tmp);
Status = DftiComputeBackward(handle[id], (MKL_Complex8 *)&cdata[j*ldc], tmp);
rdata[j*ldr] = tmp[0];
if (sign < 0) {
for (i=1; i<n1; i++) rdata[j*ldr+i] = -sign*tmp[n1-i];
......
FFTlib/rc1fft.c
+ 18
− 19
View file @ f8db7205
......@@ -57,10 +57,16 @@ void rc1fft(REAL *rdata, complex *cdata, int n, int sign)
static REAL *work, *table, scale=1.0;
REAL scl, *data;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
static DFTI_DESCRIPTOR_HANDLE handle[MAX_NUMTHREADS];
static int nprev[MAX_NUMTHREADS];
MKL_LONG Status;
int i;
int i, id;
#endif
#ifdef _OPENMP
id = omp_get_thread_num();
#else
id = 0;
#endif
#if defined(HAVE_LIBSCS)
......@@ -102,40 +108,33 @@ 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);
if (n != nprev[id]) {
DftiFreeDescriptor(&handle[id]);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n);
Status = DftiCreateDescriptor(&handle[id], 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);
Status = DftiSetValue(handle[id], DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
/*
Status = DftiSetValue(handle, DFTI_FORWARD_DOMAIN, DFTI_REAL);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
*/
Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
Status = DftiSetValue(handle[id], DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
Status = DftiCommitDescriptor(handle[id]);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n;
nprev[id] = n;
}
Status = DftiComputeForward(handle, rdata, (MKL_Complex8 *)cdata);
Status = DftiComputeForward(handle[id], rdata, (MKL_Complex8 *)cdata);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiComputeForward FAIL\n");
......
FFTlib/rcmfft.c
+ 19
− 12
View file @ f8db7205
......@@ -69,12 +69,19 @@ void rcmfft(REAL *rdata, complex *cdata, int n1, int n2, int ldr, int ldc, int s
static REAL *work;
REAL scl, *data;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
static DFTI_DESCRIPTOR_HANDLE handle[MAX_NUMTHREADS];
static int nprev[MAX_NUMTHREADS];
MKL_LONG Status;
int i,j;
int i, j, id;
#endif
#ifdef _OPENMP
id = omp_get_thread_num();
#else
id = 0;
#endif
#if defined(HAVE_LIBSCS)
nmp = mp_my_threadnum();
if(nmp>=MAX_NUMTHREADS) {
......@@ -130,39 +137,39 @@ void rcmfft(REAL *rdata, complex *cdata, int n1, int n2, int ldr, int ldc, int s
}
free(data);
#elif defined(MKL)
if (n1 != nprev) {
DftiFreeDescriptor(&handle);
if (n1 != nprev[id]) {
DftiFreeDescriptor(&handle[id]);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n1);
Status = DftiCreateDescriptor(&handle[id], 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);
Status = DftiSetValue(handle[id], 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);
Status = DftiSetValue(handle[id], DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
Status = DftiCommitDescriptor(handle[id]);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n1;
nprev[id] = n1;
}
Status = DftiComputeForward(handle, rdata, (MKL_Complex8 *)cdata);
Status = DftiComputeForward(handle[id], rdata, (MKL_Complex8 *)cdata);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiComputeForward FAIL\n");
}
for (j=0; j<n2; j++) {
Status = DftiComputeForward(handle, &rdata[j*ldr], (MKL_Complex8 *)&cdata[j*ldc]);
Status = DftiComputeForward(handle[id], &rdata[j*ldr], (MKL_Complex8 *)&cdata[j*ldc]);
for (i=1; i<((n1-1)/2)+1; i++) {
cdata[j*ldc+i].i *= -sign;
}
......
marchenko/Makefile
+ 24
− 4
View file @ f8db7205
......@@ -7,7 +7,7 @@ include ../Make_include
#ALL: fmute marchenko marchenko2
ALL: fmute marchenko
ALL: fmute marchenko marchenko_primaries
SRCJ = fmute.c \
getFileInfo.c \
......@@ -35,6 +35,20 @@ SRCH = marchenko.c \
docpkge.c \
getpars.c
SRCP = marchenko_primaries.c \
getFileInfo.c \
readData.c \
readShotData.c \
readTinvData.c \
writeData.c \
writeDataIter.c \
wallclock_time.c \
name_ext.c \
verbosepkg.c \
atopkge.c \
docpkge.c \
getpars.c
OBJJ = $(SRCJ:%.c=%.o)
fmute: $(OBJJ)
......@@ -45,22 +59,28 @@ OBJH = $(SRCH:%.c=%.o)
marchenko: $(OBJH)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko $(OBJH) $(LIBS)
OBJP = $(SRCP:%.c=%.o)
marchenko_primaries: $(OBJP)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko_primaries $(OBJP) $(LIBS)
OBJH2 = $(SRCH2:%.c=%.o)
marchenko2: $(OBJH2)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko2 $(OBJH2) $(LIBS)
install: fmute marchenko
install: fmute marchenko marchenko_primaries
cp fmute $B
cp marchenko $B
cp marchenko_primaries $B
# cp marchenko2 $B
clean:
rm -f core fmute $(OBJJ) marchenko $(OBJH) marchenko2 $(OBJH2)
rm -f core fmute $(OBJJ) marchenko $(OBJH) marchenko2 $(OBJH2) marchenko_primaries $(OBJP)
realclean: clean
rm -f $B/fmute $B/marchenko $B/marchenko2
rm -f $B/fmute $B/marchenko $B/marchenko2 $B/marchenko_primaries
......
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