Community
cancel
Showing results for
Did you mean: Beginner
140 Views

## Fortran FFT with MKL

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)

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

2 Replies Employee
140 Views

Hi

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

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. 