FFT and MKL Problems

[JohnNichols]

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.

[Evgueni_P_Intel]

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.

[JohnNichols]

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
    

[Steven_L_Intel1]

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?

[JohnNichols]

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

[Ying_H_Intel · Accepted solution]

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

[JohnNichols]

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

[Evgueni_P_Intel]

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.

[JohnNichols]

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

[JohnNichols]

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

[JohnNichols]

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

I like this manual it is good.