Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library
- Fast Discrete Fourier Transform with MKL

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

Harry_G_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-30-2014
05:25 PM

181 Views

Fast Discrete Fourier Transform with MKL

Hi all,

In R programming, there is the "fft" function: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/fft.html

> x <- c(102, 55, 89, 12, 3, 45, 9) > fft(x) [1] 315.00000+ 0.00000i 98.57101-82.76603i -23.61882-18.71932i [4] 124.54781+ 5.66758i 124.54781- 5.66758i -23.61882+18.71932i [7] 98.57101+82.76603i > Re(fft(x)) [1] 315.00000 98.57101 -23.61882 124.54781 124.54781 -23.61882 98.57101

What I have to is to replicate the *Re(fft(x)) *output in Fortran. I wonder if Intel MKL could do the job.

! testfft.f90 program myfft implicit none integer, dimension(7) :: x double precision, dimension(7) :: Refftx x = (/ 102, 55, 89, 12, 3, 45, 9 /) ! ============================================================ ! ! Here, I need to use a fft function to return the result in R ! ! ============================================================ ! ! For example, Refftx = theFFTfunctionIdonotKnow(x) print*, Refftx end program myfft

How do I have to build/compile the "testfft.f90" file?

The machine I use is equipped with **Intel Parallel Studio XE 2013**. I generally compile .f90 files by typing *ifort myfile.f90 *on the **Intel 64 Visual Studio 2008 mode **command prompt. The machine runs Window 7 64 bit.

Thank you very much!

Link Copied

5 Replies

Ying_H_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-30-2014
06:47 PM

181 Views

Hi Harry,

As i understand, your question may include several subquestions

1) first , replicate the *Re(fft(x)) *output in Fortran using MKL function

The step is not complex. MKL provide example code (both c and fortran code) in MKL install directory, for example, dftf\source\basic_dp_real_dft_1d.f90.

And the sample also include compile and link command, you can follow the setting to build the code.

simply, it is >ifort mkl_dfti.f90, basic_dp_real_dft_1d.f90 -mkl

I change the code of basic_dp_real_dft_1d.f90.

print *,"Initialize data for real-to-complex FFT"

! call init_r(x_real, N, H)

x_real = (/ 102, 55, 89, 12, 3, 45, 9 /)

print *, x_real

print *,"Compute forward transform"

status = DftiComputeForward(hand, x_real, x_cmplx)

if (0 /= status) goto 999

print *,"Verify the result"

print *, x_cmplx

print *, Real(x_cmplx)

Then you will see the result ( CCE storage). and how to call MKL function to implement the FFT computation.

2) second questions, integrate MKL to R environment.

we have two usage model:

1. build R with MKL, it will use MKL blas and Lapack funtion

2 If you want to use MKL FFT in R, there may some ways, the below is one of them we had tried.

Writing R extensions with MKL.

–KB article: Extending R with Intel® MKL

http://software.intel.com/en-us/articles/extending-r-with-intel-mkl

http://software.intel.com/en-us/articles/extending-r-with-intel-mkl

It is c language, but it should be same with Fortran language.

Hope it helps and please share if you have any progress

Best Regards,

Ying

Harry_G_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-30-2014
10:23 PM

181 Views

**<Question section 1:>**

I can modify the **basic_dp_real_dft_1d.f90** file here in **C:\Program Files (x86)\Intel\Composer XE 2013\mkl\examples\dftf\source**, as you suggested. Also, I noticed that the **mkl_dfti.f90 **file is located in **C:\Program Files (x86)\Intel\Composer XE 2013\mkl\include**. And the Command Prompt looks like:

Now, (1) which folder do I have to move to? (2) Then, is it the command > **ifort mkl_dfti.f90, basic_dp_real_dft_1d.f90 -mkl** , including a comma between the two .f90 files ? Or, do I have to copy the two files to my working folder and **ifort **there**?**

**<Question section 2:>**

I have a bunch of tasks before calculating the *Re(fft(x))*. Isn't there a way to retain my existing *main.f90 *code and call the FFT operation from the *main.f90*.

**EDIT:**

I have an update for the <Question section 1>. I happened to find https://software.intel.com/en-us/forums/topic/336695. The author used /Qmkl. I don't know what the heck is /Qmkl. The **/mkl** alone without the **Q** did not work. Anyway, the **ifort **command goes like:

**>ifort "C:\Program Files (x86)\Intel\Composer XE 2013\mkl\include\mkl_dfti.f90" dft1d.f90 -o dft1d /Qmkl**

The above is to replicate Re(fft(c(102,55,89,12,3,45,9,9,55))) in R. Something seemingly goes wrong. I attach the **dft1d.f90 **file, hoping someone is out there who is interested in this problem.

Thank you.

Ying_H_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

11-30-2014
11:46 PM

181 Views

Hi Harry,

Right, you can copy the example code to your work folder and compile it with /Qmkl option. ( it is compiler option)

And the FFT result you are seeing is stored by CCE format. You may see the details in MKL manual =>

DFTI_CONJUGATE_EVEN_STORAGE: storage scheme for a conjugate-even domain

The Intel MKL FFT interface supports two configuration values for this parameter:

Then just complete the result.

Best Regards,

Ying

Dmitry_B_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

12-01-2014
01:37 AM

181 Views

Harry,

That 'Verification of result' failed after you modified the example is no worry. Verifier checks the output of the transform assuming the input was specially crafted. Since you give your own input, the verifier complaints.

Unlike R, in Fortran the result contains only a half of full data, because the other half can be easily restored by conjugate-even symmetry of the result (follow Ying's note for details).

Thanks

Dima

Ying_H_Intel

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

01-03-2015
05:34 PM

181 Views

Hi Harry,

Just complete the thread,

after get FFT result which is CCE format, then use the below code to get the whole output or Re(z). please see MKL manual => chapter 12 => DFTI_CONJUGATE_EVEN_STORAGE: storage scheme for a conjugate-even domain.

The following examples illustrate usage of the DFTI_COMPLEX_COMPLEXstorage for a conjugate-even

domain:

**! Fortran
real :: AR(N1,N2,M) ! Array containing values of R
complex :: AZ(N1/2+1,N2,M) ! Array containing values of Z
...
! on input: R{k1,k2,m} = AR(k1,k2,m)
status = DftiComputeForward( desc, AR(:,1,1), AZ(:,1,1) )
! on output:
! for k1=1 … N1/2+1: Z{k1,k2,m} = AZ(k1,k2,m)
! for k1=N1/2+2 … N1: Z{k1,k2,m} = conj(AZ(mod(N1-k1+1,N1)+1,mod(N2-k2+1,N2)+1,m))**

// C/C++

float *AR = malloc( N1*N2*M * sizeof(AR[0]) );

complex *AZ = malloc( N1*(N2/2+1)*M * sizeof(AZ[0]) );

MKL_LONG is[3], os[3], idist, odist; // input and output strides and distance

...

// on input: R{k1,k2,m}

// = AR[is[0] + k1*is[1] + k2*is[2] + m*idist]

status = DftiComputeForward( desc, R, C );

// on output:

// for k2=0…N2/2: Z{k1,k2,m} = AZ[os[0]+k1*os[1]+k2*os[2]+m*odist]

// for k2=N2/2+1…N2-1: Z{k1,k2,m} = conj(AZ[os[0]+(N1-k1)%N1*os[1]

// +(N2-k2)%N2*os[2]+m*odist])

Dima give one example to convert the general input to CCE formart as below., you may do the inverse operation.

Best Regards,

Ying

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

For more complete information about compiler optimizations, see our Optimization Notice.