- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi,
I am working on software for acoustics simulation, which requires 3D Discrete Cosine Transforms to be performed. FFTW3 provides the fftwf_plan_r2r_3d function to create plans for just this purpose:
[cpp]
fftwf_plan plan = fftwf_plan_r2r_3d(Nx, Ny, Nz, gInData, gOutData, FFTW_REDFT10, FFTW_REDFT10, FFTW_REDFT10, FFTW_PRESERVE_INPUT);
fftwf_execute(plan);
[/cpp]
However, MKL does not support 3D DCTs. So I tried to implement the 3D DCT in terms of multiple 1D DCTs. For some of these 1D DCTs, I need to specify strides between consecutive elements. Again, FFTW3 allows this via the function fftwf_plan_many_r2r, which also doesn't seem to be supported in MKL (the returned plans from MKL's FFTW3 wrappers are always NULL).
[cpp]
fftwf_r2r_kind kind = FFTW_REDFT10;
fftwf_plan planX = fftwf_plan_many_r2r(1, &Nx, 1, gInData, NULL, Ny*Nz, 1, gOutData, NULL, Ny*Nz, 1, &kind, FFTW_PRESERVE_INPUT);
fftwf_plan planY = fftwf_plan_many_r2r(1, &Ny, 1, gInData, NULL, Nz, 1, gOutData, NULL, Nz, 1, &kind, FFTW_PRESERVE_INPUT);
fftwf_plan planZ = fftwf_plan_many_r2r(1, &Nz, 1, gInData, NULL, 1, 1, gOutData, NULL, 1, 1, &kind, FFTW_PRESERVE_INPUT);
for (int x = 0; x < Nx; x++)
for (int y = 0; y < Ny; y++)
fftwf_execute_r2r(planZ, &gInData[x*Ny*Nz + y*Nz], &gOutData[x*Ny*Nz + y*Nz]);
for (int x = 0; x < Nx; x++)
for (int z = 0; z < Nz; z++)
fftwf_execute_r2r(planY, &gOutData[x*Ny*Nz + z], &gInData[x*Ny*Nz + z]);
for (int y = 0; y < Ny; y++)
for (int z = 0; z < Nz; z++)
fftwf_execute_r2r(planX, &gInData[y*Nz + z], &gOutData[y*Nz + z]);
[/cpp]
I next tried the raw Trigonometric Transforms functions, s_init_trig_transform and the like. These don't seem to offer any way of setting a stride between elements. Since s_commit_trig_transform returns a DFTI_DESCRIPTOR_HANDLE, I tried using DftiSetValue to specify a stride on the descriptor, followed by a call to DftiCommitDescriptor.
[cpp]
int tt_type = MKL_STAGGERED_COSINE_TRANSFORM;
int ir;
DFTI_DESCRIPTOR_HANDLE dftiX, dftiY, dftiZ;
int iparX[128];
float* sparX = new float[3*Nx/2 + 2];
s_init_trig_transform(&Nx, &tt_type, iparX, sparX, &ir);
s_commit_trig_transform(gInData, &dftiX, iparX, sparX, &ir);
int xStride[2] = {0, Ny*Nz};
DftiSetValue(dftiX, DFTI_INPUT_STRIDES, xStride);
DftiSetValue(dftiX, DFTI_OUTPUT_STRIDES, xStride);
DftiCommitDescriptor(dftiX);
int iparY[128];
float* sparY = new float[3*Ny/2 + 2];
s_init_trig_transform(&Ny, &tt_type, iparY, sparY, &ir);
s_commit_trig_transform(gInData, &dftiY, iparY, sparY, &ir);
int yStride[2] = {0, Nz};
DftiSetValue(dftiY, DFTI_INPUT_STRIDES, yStride);
DftiSetValue(dftiY, DFTI_OUTPUT_STRIDES, yStride);
DftiCommitDescriptor(dftiY);
int iparZ[128];
float* sparZ = new float[3*Nz/2 + 2];
s_init_trig_transform(&Nz, &tt_type, iparZ, sparZ, &ir);
s_commit_trig_transform(gInData, &dftiZ, iparZ, sparZ, &ir);
iparX[6] = 0;
iparY[6] = 0;
iparZ[6] = 0;
iparX[10] = 1;
iparY[10] = 1;
iparZ[10] = 1;
iparX[14] = Nx;
iparY[14] = Ny;
iparZ[14] = Nz;
memcpy(gOutData, gInData, Nx*Ny*Nz*sizeof(float));
for (int x = 0; x < Nx; x++)
for (int y = 0; y < Ny; y++)
s_backward_trig_transform(&gOutData[x*Ny*Nz + y*Nz], &dftiZ, iparZ, sparZ, &ir);
for (int x = 0; x < Nx; x++)
for (int z = 0; z < Nz; z++)
s_backward_trig_transform(&gOutData[x*Ny*Nz + z], &dftiY, iparY, sparY, &ir);
for (int y = 0; y < Ny; y++)
for (int z = 0; z < Nz; z++)
s_backward_trig_transform(&gOutData[y*Nz + z], &dftiX, iparX, sparX, &ir);
[/cpp]
However, this does not produce correct results (i.e., the results don't match at all with the FFTW 3D DCT or the FFTW multiple-1D-DCT results). I am using MKL 10.3 Update 9 on Windows (32-bit).
Have I missed something here? Is there a way to either a) specify strides for 1D DCTs, or b) perform 3D DCTs using MKL? In the absence of a way to use MKL to compute 3D DCTs, I will be forced to switch away from using MKL.
Thanks,
Lakulish Antani.
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page