-
Notifications
You must be signed in to change notification settings - Fork 3
/
zfft3d.f
133 lines (133 loc) · 3.63 KB
/
zfft3d.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
C
C FFTE: A FAST FOURIER TRANSFORM PACKAGE
C
C (C) COPYRIGHT SOFTWARE, 2000-2004, 2008-2011, 2020
C ALL RIGHTS RESERVED
C BY
C DAISUKE TAKAHASHI
C CENTER FOR COMPUTATIONAL SCIENCES
C UNIVERSITY OF TSUKUBA
C 1-1-1 TENNODAI, TSUKUBA, IBARAKI 305-8577, JAPAN
C E-MAIL: [email protected]
C
C
C 3-D COMPLEX FFT ROUTINE
C
C FORTRAN77 SOURCE PROGRAM
C
C CALL ZFFT3D(A,NX,NY,NZ,IOPT)
C
C A(NX,NY,NZ) IS COMPLEX INPUT/OUTPUT VECTOR (COMPLEX*16)
C NX IS THE LENGTH OF THE TRANSFORMS IN THE X-DIRECTION (INTEGER*4)
C NY IS THE LENGTH OF THE TRANSFORMS IN THE Y-DIRECTION (INTEGER*4)
C NZ IS THE LENGTH OF THE TRANSFORMS IN THE Z-DIRECTION (INTEGER*4)
C ------------------------------------
C NX = (2**IP) * (3**IQ) * (5**IR)
C NY = (2**JP) * (3**JQ) * (5**JR)
C NZ = (2**KP) * (3**KQ) * (5**KR)
C ------------------------------------
C IOPT = 0 FOR INITIALIZING THE COEFFICIENTS (INTEGER*4)
C = -1 FOR FORWARD TRANSFORM
C = +1 FOR INVERSE TRANSFORM
C
C WRITTEN BY DAISUKE TAKAHASHI
C
SUBROUTINE ZFFT3D(A,NX,NY,NZ,IOPT)
IMPLICIT REAL*8 (A-H,O-Z)
INCLUDE 'param.h'
COMPLEX*16 A(*)
COMPLEX*16 B((NDA3+NP)*NBLK),C(NDA3)
COMPLEX*16 WX(NDA3),WY(NDA3),WZ(NDA3)
DIMENSION LNX(3),LNY(3),LNZ(3)
SAVE WX,WY,WZ
C
IF (IOPT .EQ. 0) THEN
CALL SETTBL(WX,NX)
CALL SETTBL(WY,NY)
CALL SETTBL(WZ,NZ)
RETURN
END IF
C
IF (IOPT .EQ. 1) THEN
!$OMP PARALLEL DO
!DIR$ VECTOR ALIGNED
DO 10 I=1,NX*NY*NZ
A(I)=DCONJG(A(I))
10 CONTINUE
END IF
C
CALL FACTOR(NX,LNX)
CALL FACTOR(NY,LNY)
CALL FACTOR(NZ,LNZ)
C
!$OMP PARALLEL PRIVATE(B,C)
CALL ZFFT3D0(A,B,B,C,WX,WY,WZ,NX,NY,NZ,LNX,LNY,LNZ)
!$OMP END PARALLEL
C
IF (IOPT .EQ. 1) THEN
DN=1.0D0/(DBLE(NX)*DBLE(NY)*DBLE(NZ))
!$OMP PARALLEL DO
!DIR$ VECTOR ALIGNED
DO 20 I=1,NX*NY*NZ
A(I)=DCONJG(A(I))*DN
20 CONTINUE
END IF
RETURN
END
SUBROUTINE ZFFT3D0(A,BY,BZ,C,WX,WY,WZ,NX,NY,NZ,LNX,LNY,LNZ)
IMPLICIT REAL*8 (A-H,O-Z)
INCLUDE 'param.h'
COMPLEX*16 A(NX,NY,*)
COMPLEX*16 BY(NY+NP,*),BZ(NZ+NP,*),C(*)
COMPLEX*16 WX(*),WY(*),WZ(*)
DIMENSION LNX(*),LNY(*),LNZ(*)
C
!$OMP DO PRIVATE(I,II,K,KK)
DO 80 J=1,NY
DO 70 II=1,NX,NBLK
DO 30 KK=1,NZ,NBLK
DO 20 I=II,MIN0(II+NBLK-1,NX)
!DIR$ VECTOR ALIGNED
DO 10 K=KK,MIN0(KK+NBLK-1,NZ)
BZ(K,I-II+1)=A(I,J,K)
10 CONTINUE
20 CONTINUE
30 CONTINUE
DO 40 I=II,MIN0(II+NBLK-1,NX)
CALL FFT235(BZ(1,I-II+1),C,WZ,NZ,LNZ)
40 CONTINUE
DO 60 K=1,NZ
!DIR$ VECTOR ALIGNED
DO 50 I=II,MIN0(II+NBLK-1,NX)
A(I,J,K)=BZ(K,I-II+1)
50 CONTINUE
60 CONTINUE
70 CONTINUE
80 CONTINUE
!$OMP DO PRIVATE(I,II,J,JJ)
DO 170 K=1,NZ
DO 150 II=1,NX,NBLK
DO 110 JJ=1,NY,NBLK
DO 100 I=II,MIN0(II+NBLK-1,NX)
!DIR$ VECTOR ALIGNED
DO 90 J=JJ,MIN0(JJ+NBLK-1,NY)
BY(J,I-II+1)=A(I,J,K)
90 CONTINUE
100 CONTINUE
110 CONTINUE
DO 120 I=II,MIN0(II+NBLK-1,NX)
CALL FFT235(BY(1,I-II+1),C,WY,NY,LNY)
120 CONTINUE
DO 140 J=1,NY
!DIR$ VECTOR ALIGNED
DO 130 I=II,MIN0(II+NBLK-1,NX)
A(I,J,K)=BY(J,I-II+1)
130 CONTINUE
140 CONTINUE
150 CONTINUE
DO 160 J=1,NY
CALL FFT235(A(1,J,K),C,WX,NX,LNX)
160 CONTINUE
170 CONTINUE
RETURN
END