Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
6536 Discussions

## FFT and MKL Problems Valued Contributor II
235 Views

The program does not return from the call from the main program to the subroutine FFT().

MKLVARS.BAT does not seem to set the path or environment variables.

This is the latest preview and VS 2015.

1 Solution Employee
235 Views

Hi John,

You may fix the issue  with one  tiny place  :)

Double precision :: Y(n+2), U,  SUM, DF, W, SUMW

Because Real to Complex Conjugate-even format :  the storage of a one-dimensional (1D) size-N conjugate-even sequence in a
real array can be  the CCS, PACK, and PERM packed formats. The CCS format requires an array of size N+2.

Best Regards,
Ying

10 Replies Employee
235 Views

Could you please let us know which arguments you pass to FFT()?

It is declared as FFT(File1, FileOut, N, time, X) in FFT1.f90 and called like FFT() in Source1.f90. Valued Contributor II
235 Views

FFT is called from the main WULF program --

call FFT(sm, SOUT, NFA, time, X)

SM, sout are simply integers used to indicate files ie write(sm,*)

NFA is an integer to tell you the FFT array size

X is a

```    complex(kind=dp), ALLOCATABLE :: X(:)

INTEGER, PARAMETER :: dp = selected_real_kind(15, 307)```
```   !****************************************************************************
!
!  PROGRAM: Wulf
!
!  PURPOSE:  Entry point for the console application.
!
!****************************************************************************

program Wulf
use BASE
use FastFourierTransform
use S
use Scotia

implicit none
integer i,ISTAT,count, dummy, NG
REAL (KIND=dp) time
complex(kind=dp), ALLOCATABLE :: X(:)
REAL (KIND=dp), ALLOCATABLE :: aXYZ(:,:)

call Files(1,NG)

ALLOCATE (aXYZ(count1,count),X(NFA),STAT=istat)
IF (istat.NE.0) THEN
WRITE (*, *) '*** Could not allocate some arrays in WULF'
STOP
END IF

time = 0.1d0

do i =1,NFA
if(i .eq. 1) then
x(i) = (0.70154d0,0.0d0)
elseif(i .eq. 2) then
x(i) =(-2.5018d0,0.0d0)
elseif(i .eq. 3) then
x(i) = (-0.35385d0,0.0d0)
elseif(i .eq. 4) then
x(i) = (-0.82359d0,0.0d0)
elseif(i .eq. 5) then
x(i) = (-1.5771d0,0.0d0)
elseif(i .eq. 6) then
x(i) = (0.50797d0,0.0d0)
elseif(i .eq. 7) then
x(i) = (0.28198d0,0.0d0)
elseif(i .eq. 8) then
x(i) = (0.03348d0, 0.0d0)
end if
end do

call pardiso_unsym(aXYZ,count,1)
call FFT(sm, SOUT, NFA, time, X)

stop 'here'

end program Wulf
``` Employee
235 Views

The stack is being corrupted by the call to DftiComputeForward in these lines in routine FFT:

```    Status = DftiCreateDescriptor(My_Desc2_Handle, DFTI_DOUBLE, DFTI_REAL, 1, N)
Status = DftiCommitDescriptor(My_Desc2_Handle)
Status = DftiComputeForward(My_Desc2_Handle, Y)
```

I am not familiar enough with MKL to know why this goes wrong. Aren't there supposed to be calls to DftiSetValue somewhere? Valued Contributor II
235 Views

Steve:

Thanks - it is the 3rd line that causes the problem.

This code came from INTEL samples, I worked this out and turned off this code whilst I awaited the response.

When I ran the original INTEL code unchanged the next line after 3 crashed all the time.

I can of without this code - I worked around - but it is interesting.

Perhaps the MKL people will work it out - beyond my skill set.

But it does work well, sorry for your needing your help - but I could not post to the MKL forum.

The CHDIR command set does not work with MKL as far as I can see -- interesting - it works if MKL turned off. Something between IFPORT AND MKL use

John Employee
236 Views

Hi John,

You may fix the issue  with one  tiny place  :)

Double precision :: Y(n+2), U,  SUM, DF, W, SUMW

Because Real to Complex Conjugate-even format :  the storage of a one-dimensional (1D) size-N conjugate-even sequence in a
real array can be  the CCS, PACK, and PERM packed formats. The CCS format requires an array of size N+2.

Best Regards,
Ying Valued Contributor II
235 Views

Ying:

I copied that part of the program from your examples for MKL - perhaps you could get someone to fix the error there in the 1D FFT example, so this does not pop up again.

That line caused a problem from the first time I compiled.

John Employee
235 Views

Hi John,

We provide two examples showing 1D real-to-complex FFTs using CCE format: basic_dp_real_dft_1d.f90 and basic_sp_real_dft_1d.f90.

Both examples correctly allocate an array of N real numbers for input and an array of N/2+1 complex numbers for output.

Evgueni. Valued Contributor II
235 Views
```program basic_dp_complex_dft_1d

use MKL_DFTI, forget => DFTI_DOUBLE, DFTI_DOUBLE => DFTI_DOUBLE_R

! Size of 1D transform
integer, parameter :: N = 7

! Arbitrary harmonic used to verify FFT
integer, parameter :: H = 1

! Working precision is double precision
integer, parameter :: WP = selected_real_kind(15,307)

! Execution status
integer :: status = 0, ignored_status

! The data array
complex(WP), allocatable :: x (:)

! DFTI descriptor handle
type(DFTI_DESCRIPTOR), POINTER :: hand

hand => null()``` Valued Contributor II
235 Views

Evgueni.

There are another set of FFT examples for FFT - not just the ones in MKL Examples

I cannot remember where I found them on your website, but I think it was through a quicklink in the MKL 11.3? manual.

This code is different

John

```Subroutine FFT(File1, FileOut, N, time, X)
! Fortran example.
! 1D complex to complex, and real to conjugate-even
Use MKL_DFTI

implicit none
integer n
integer FileOut
Double Complex :: X(n), XC(N)
Double precision :: Y(n+2), U,  SUM, DF, W, SUMW, Res
Double precision :: YStart, time, FFTRatio, A,B,C, phi
integer i,j, File1,n1
type(DFTI_DESCRIPTOR), POINTER :: My_Desc1_Handle, My_Desc2_Handle
Integer :: Status

!...put input data into X(1),...,X(32); Y(1),...,Y(32)
! Perform a complex to complex transform
j = 1
YStart = 0.0d0
FFTRatio = time
DF = 1.0d0/(time*dfloat(n))
sum = 0.0d0
sumw = 0.0d0
phi = 0.0d0

do i =1,(N-1)
y(i) = real(x(i))
!write(*,*)i,x(i)
end do
Status = DftiCreateDescriptor( My_Desc1_Handle, DFTI_DOUBLE,DFTI_COMPLEX, 1, N )
Status = DftiCommitDescriptor( My_Desc1_Handle )
Status = DftiComputeForward( My_Desc1_Handle, X )
Status = DftiFreeDescriptor(My_Desc1_Handle)

do i =1,((N-1)/2)+1
x(i) = FFTRatio*x(i)
XC(i) = CONJG(x(i))
U = x(i)*xc(I)
sum = sum+U*df
sumw = sumw+(y(i)*y(i))*time
n1 = 1
a = 1.0d0
b = 2.0d0
call vdhypot( n1, real(x(i)), aimag(x(i)), c )
call vdAtan2(n1,aimag(x(i)),real(x(i)),phi)
write(FileOUT,10)i,(df*float(i)),(real(x(i))),aimag(x(i)),y(i),sum,sumw,c,phi,c*c
10       format(i5,2(",",F14.8),",",F14.8,",",F14.8,",",F14.8,",",F14.8,",",F14.8,",",F14.8,",",F14.8)
end do
! result is given by {X(1),X(2),...,X(32)}
! Perform a real to complex conjugate-even transform
Status = DftiCreateDescriptor(My_Desc2_Handle, DFTI_DOUBLE, DFTI_REAL, 1, N)
Status = DftiCommitDescriptor(My_Desc2_Handle)
Status = DftiComputeForward(My_Desc2_Handle, Y)
Status = DftiFreeDescriptor(My_Desc2_Handle)
! result is given in CCS format.
return
end subroutine FFT
END module``` Valued Contributor II
235 Views

Page 3166 in the manual is the code I used -- see manual attached

I like this manual it is good. 