Fortran FFT with MKL

Matt_S_ · ‎01-21-2015

I'm currently messing with some code that uses an old F77 looking FFT subroutine full of goto's and other nasty stuff. I have been looking around at using MKL's FFT routines, but I am having trouble finding some examples of how they are used. My current FFT routine uses only double precision data, so I don't need a FFT dealing with complex. Can anyone point me to any example code using a 2D MKL FFT routine? (possibly within /opt/intel/mkl/examples)

Thanks in advance!

(For any interested, here is the current FFT routine being used)

	SUBROUTINE FFT(X,Y,N,ITYPE)
	implicit none
	INTEGER N,M,ITYPE
	real*8 X(*), Y(*)
	INTEGER   I, J, K, N1, N2, N4, IS, ID, I0, I1, I2, I3
	real*8   TWOPI, E, A, A3, CC1, SS1, CC3, SS3
	real*8   R1, R2, S1, S2, S3, XT
	INTRINSIC   SIN, COS
	PARAMETER   ( TWOPI = 6.2831853071795864769 )

       IF ( N .EQ. 1 ) RETURN

	   M=nint(log(real(N))/log(2.0))

       IF ( ITYPE .EQ. -1 ) THEN
	   DO 1, I = 1, N
       Y(I) = - Y(I)
1      CONTINUE
       ENDIF

       N2 = 2 * N
       DO 10, K = 1, M-1
		N2 = N2 / 2
		N4 = N2 / 4
		E = TWOPI / N2
		A = 0.0
		DO 20, J = 1, N4
		A3 = 3 * A
		CC1 = COS( A )
		SS1 = SIN( A )
		CC3 = COS( A3 )
		SS3 = SIN( A3 )
		A = J * E
		IS = J
		ID = 2 * N2
40       DO 30, I0 = IS, N-1, ID
         I1 = I0 + N4
         I2 = I1 + N4
         I3 = I2 + N4
         R1 = X(I0) - X(I2)
         X(I0) = X(I0) + X(I2)
         R2 = X(I1) - X(I3)
         X(I1) = X(I1) + X(I3)
         S1 = Y(I0) - Y(I2)
         Y(I0) = Y(I0) + Y(I2)
         S2 = Y(I1) - Y(I3)
         Y(I1) = Y(I1) + Y(I3)
         S3 = R1 - S2
         R1 = R1 + S2
         S2 = R2 - S1
         R2 = R2 + S1
         X(I2) = R1 * CC1 - S2 * SS1
         Y(I2) = - S2 * CC1 - R1 * SS1
         X(I3) = S3 * CC3 + R2 * SS3
         Y(I3) = R2 * CC3 - S3 * SS3
30       CONTINUE
		IS = 2 * ID - N2 + J
		ID = 4 * ID
		IF ( IS .LT. N ) GOTO 40
20      CONTINUE
10    CONTINUE

!C--------LAST STAGE, LENGTH-2 BUTTERFLY ----------------------C
       IS = 1
       ID = 4
50     DO 60, I0 = IS, N, ID
		I1 = I0 + 1
		R1 = X(I0)
		X(I0) = R1 + X(I1)
		X(I1) = R1 - X(I1)
		R1 = Y(I0)
		Y(I0) = R1 + Y(I1)
		Y(I1) = R1 - Y(I1)
60    CONTINUE
       IS = 2 * ID - 1
       ID = 4 * ID
       IF ( IS .LT. N ) GOTO 50

!C-------BIT REVERSE COUNTER-----------------------------------C
100    J = 1
       N1 = N - 1
       DO 104, I = 1, N1
		IF ( I .GE. J ) GOTO 101
		XT = X(J)
		X(J) = X(I)
		X(I) = XT
		XT = Y(J)
		Y(J) = Y(I)
		Y(I) = XT
101     K = N / 2
102     IF ( K .GE. J ) GOTO 103
		J = J - K
		K = K / 2
		GOTO 102
103     J = J + K
104   CONTINUE

	IF ( ITYPE .EQ. -1 ) THEN
	do I = 1, N
      Y(I) = - Y(I) 
	end do
	else
	do i=1,n
	x(i)=x(i)/n
	y(i)=y(i)/n
	end do
    ENDIF

    RETURN
    END

I am unfamiliar with "length 2 butterfly" and "bit reverse counter" in FFT. If there is any specific MKL FFT routine that performs these same calculations?

Zhang_Z_Intel · ‎01-21-2015

Hi

Looks like this page is exactly what you need. 2D FFT examples are found in the middle of the page.

https://software.intel.com/en-us/node/522290#AC30C86F-9E01-4A5F-B086-CE7FBCE2126A ;

Matt_S_ · ‎01-23-2015

Zhang, thank you for the link! I have discovered where the FFT currently used came from and it is no standard FFT. I could not replicate the results from any of the MKL FFT routine examples. The FFT used can be found here: https://github.com/syntheticpp/benchfft/blob/master/benchees/burrus/sffteu.f and it is described as, "a slight modification of a complex split radix FFT routine presented by C.S. Burrus."

Does anyone know if a variation of MKL's FFT performs anything of the sort?

Thanks again.