topic Page 3166 in the manual is in Intel® oneAPI Math Kernel Library

FFT and MKL Problems

JohnNichols — Mon, 30 May 2016 21:57:32 GMT

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.

Could please let us know

Evgueni_P_Intel — Wed, 01 Jun 2016 19:41:00 GMT

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.

FFT is called from the main

JohnNichols — Wed, 01 Jun 2016 21:12:01 GMT

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)
    read(sRA,*)count

    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

The stack is being corrupted

Steven_L_Intel1 — Thu, 02 Jun 2016 20:21:33 GMT

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?

Steve:

JohnNichols — Thu, 02 Jun 2016 21:37:45 GMT

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

Hi John,

Ying_H_Intel — Fri, 03 Jun 2016 04:21:46 GMT

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

Ying:

JohnNichols — Fri, 03 Jun 2016 14:14:50 GMT

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

Hi John,

Evgueni_P_Intel — Fri, 03 Jun 2016 17:09:50 GMT

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.

program basic_dp_complex_dft

JohnNichols — Fri, 03 Jun 2016 18:08:10 GMT

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()

Evgueni.

JohnNichols — Fri, 03 Jun 2016 18:13:04 GMT

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

Page 3166 in the manual is

JohnNichols — Fri, 03 Jun 2016 18:18:09 GMT

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

I like this manual it is good.