Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : MODULE fft_lib
8 :
9 : USE fft_kinds, ONLY: dp
10 : USE fft_plan, ONLY: fft_plan_type
11 : USE fftsg_lib, ONLY: fftsg1dm,&
12 : fftsg3d,&
13 : fftsg_do_cleanup,&
14 : fftsg_do_init,&
15 : fftsg_get_lengths
16 : USE fftw3_lib, ONLY: &
17 : fft_alloc => fftw_alloc, fft_dealloc => fftw_dealloc, fftw31dm, fftw33d, &
18 : fftw3_create_plan_1dm, fftw3_create_plan_3d, fftw3_destroy_plan, fftw3_do_cleanup, &
19 : fftw3_do_init, fftw3_get_lengths
20 : #include "../../base/base_uses.f90"
21 :
22 : IMPLICIT NONE
23 : PRIVATE
24 :
25 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fft_lib'
26 :
27 : PUBLIC :: fft_do_cleanup, fft_do_init, fft_get_lengths, fft_create_plan_3d
28 : PUBLIC :: fft_create_plan_1dm, fft_1dm, fft_library, fft_3d, fft_destroy_plan
29 : PUBLIC :: fft_alloc, fft_dealloc
30 :
31 : CONTAINS
32 : ! **************************************************************************************************
33 : !> \brief Interface to FFT libraries
34 : !> \param fftlib ...
35 : !> \return ...
36 : !> \par History
37 : !> IAB 09-Jan-2009 : Modified to use fft_plan_type
38 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2009 on behalf of the HECToR project
39 : !> \author JGH
40 : ! **************************************************************************************************
41 9909 : FUNCTION fft_library(fftlib) RESULT(flib)
42 :
43 : CHARACTER(len=*), INTENT(IN) :: fftlib
44 : INTEGER :: flib
45 :
46 : SELECT CASE (fftlib)
47 : CASE DEFAULT
48 12 : flib = -1
49 : CASE ("FFTSG")
50 12 : flib = 1
51 : CASE ("FFTW3")
52 9909 : flib = 3
53 : END SELECT
54 :
55 9909 : END FUNCTION fft_library
56 :
57 : ! **************************************************************************************************
58 : !> \brief ...
59 : !> \param fft_type ...
60 : !> \param wisdom_file ...
61 : ! **************************************************************************************************
62 9897 : SUBROUTINE fft_do_init(fft_type, wisdom_file)
63 : INTEGER, INTENT(IN) :: fft_type
64 : CHARACTER(LEN=*), INTENT(IN) :: wisdom_file
65 :
66 9897 : SELECT CASE (fft_type)
67 : CASE DEFAULT
68 0 : CPABORT("fft_do_init")
69 : CASE (1)
70 12 : CALL fftsg_do_init()
71 : CASE (3)
72 9897 : CALL fftw3_do_init(wisdom_file)
73 : END SELECT
74 :
75 9897 : END SUBROUTINE fft_do_init
76 :
77 : ! **************************************************************************************************
78 : !> \brief ...
79 : !> \param fft_type ...
80 : !> \param wisdom_file ...
81 : !> \param ionode ...
82 : ! **************************************************************************************************
83 9687 : SUBROUTINE fft_do_cleanup(fft_type, wisdom_file, ionode)
84 : INTEGER, INTENT(IN) :: fft_type
85 : CHARACTER(LEN=*), INTENT(IN) :: wisdom_file
86 : LOGICAL, INTENT(IN) :: ionode
87 :
88 9687 : SELECT CASE (fft_type)
89 : CASE DEFAULT
90 0 : CPABORT("fft_do_cleanup")
91 : CASE (1)
92 12 : CALL fftsg_do_cleanup()
93 : CASE (3)
94 9687 : CALL fftw3_do_cleanup(wisdom_file, ionode)
95 : END SELECT
96 :
97 9687 : END SUBROUTINE fft_do_cleanup
98 :
99 : ! **************************************************************************************************
100 : !> \brief ...
101 : !> \param fft_type ...
102 : !> \param DATA ...
103 : !> \param max_length ...
104 : ! **************************************************************************************************
105 115694 : SUBROUTINE fft_get_lengths(fft_type, DATA, max_length)
106 :
107 : INTEGER, INTENT(IN) :: fft_type
108 : INTEGER, DIMENSION(*) :: DATA
109 : INTEGER, INTENT(INOUT) :: max_length
110 :
111 115694 : SELECT CASE (fft_type)
112 : CASE DEFAULT
113 0 : CPABORT("fft_get_lengths")
114 : CASE (1)
115 115166 : CALL fftsg_get_lengths(DATA, max_length)
116 : CASE (3)
117 115694 : CALL fftw3_get_lengths(DATA, max_length)
118 : END SELECT
119 :
120 115694 : END SUBROUTINE fft_get_lengths
121 :
122 : ! **************************************************************************************************
123 :
124 : ! **************************************************************************************************
125 : !> \brief ...
126 : !> \param plan ...
127 : !> \param fft_type ...
128 : !> \param fft_in_place ...
129 : !> \param fsign ...
130 : !> \param n ...
131 : !> \param zin ...
132 : !> \param zout ...
133 : !> \param plan_style ...
134 : ! **************************************************************************************************
135 56436 : SUBROUTINE fft_create_plan_3d(plan, fft_type, fft_in_place, fsign, n, zin, zout, plan_style)
136 :
137 : TYPE(fft_plan_type), INTENT(INOUT) :: plan
138 : INTEGER, INTENT(IN) :: fft_type
139 : LOGICAL, INTENT(IN) :: fft_in_place
140 : INTEGER, INTENT(IN) :: fsign
141 : INTEGER, DIMENSION(3), INTENT(IN) :: n
142 : COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zin, zout
143 : INTEGER, INTENT(IN) :: plan_style
144 :
145 56436 : plan%fft_type = fft_type
146 56436 : plan%fsign = fsign
147 56436 : plan%fft_in_place = fft_in_place
148 225744 : plan%n_3d = n
149 56436 : !$ plan%need_alt_plan = .FALSE.
150 :
151 : ! Planning only needed for FFTW3
152 56436 : IF (fft_type == 3) THEN
153 56292 : CALL fftw3_create_plan_3d(plan, zin, zout, plan_style)
154 56292 : plan%valid = .TRUE.
155 : END IF
156 :
157 56436 : END SUBROUTINE fft_create_plan_3d
158 :
159 : !
160 : ! really ugly, plan is intent out, because plan%fsign is also a status flag
161 : ! if something goes wrong, plan%fsign is set to zero, and the plan becomes invalid
162 : !
163 : ! **************************************************************************************************
164 : !> \brief ...
165 : !> \param plan ...
166 : !> \param scale ...
167 : !> \param zin ...
168 : !> \param zout ...
169 : !> \param stat ...
170 : ! **************************************************************************************************
171 653108 : SUBROUTINE fft_3d(plan, scale, zin, zout, stat)
172 : TYPE(fft_plan_type), INTENT(IN) :: plan
173 : REAL(KIND=dp), INTENT(IN) :: scale
174 : COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zin, zout
175 : INTEGER, INTENT(OUT) :: stat
176 :
177 653108 : stat = plan%fsign
178 653108 : IF (plan%n_3d(1)*plan%n_3d(2)*plan%n_3d(3) > 0) THEN
179 653108 : SELECT CASE (plan%fft_type)
180 : CASE DEFAULT
181 0 : CPABORT("fft_3d")
182 : CASE (1)
183 2884 : CALL fftsg3d(plan%fft_in_place, stat, scale, plan%n_3d, zin, zout)
184 : CASE (3)
185 653108 : CALL fftw33d(plan, scale, zin, zout, stat)
186 : END SELECT
187 : END IF
188 : ! stat is set to zero on error, -1,+1 are OK
189 653108 : IF (stat == 0) THEN
190 0 : stat = 1
191 : ELSE
192 653108 : stat = 0
193 : END IF
194 :
195 653108 : END SUBROUTINE fft_3d
196 :
197 : ! **************************************************************************************************
198 :
199 : ! **************************************************************************************************
200 : !> \brief ...
201 : !> \param plan ...
202 : !> \param fft_type ...
203 : !> \param fsign ...
204 : !> \param trans ...
205 : !> \param n ...
206 : !> \param m ...
207 : !> \param zin ...
208 : !> \param zout ...
209 : !> \param plan_style ...
210 : ! **************************************************************************************************
211 163476 : SUBROUTINE fft_create_plan_1dm(plan, fft_type, fsign, trans, n, m, zin, zout, plan_style)
212 : TYPE(fft_plan_type), INTENT(INOUT) :: plan
213 : INTEGER, INTENT(IN) :: fft_type, fsign
214 : LOGICAL, INTENT(IN) :: trans
215 : INTEGER, INTENT(IN) :: n, m
216 : COMPLEX(KIND=dp), DIMENSION(*), INTENT(IN) :: zin, zout
217 : INTEGER, INTENT(IN) :: plan_style
218 :
219 163476 : plan%fft_type = fft_type
220 163476 : plan%fsign = fsign
221 163476 : plan%trans = trans
222 163476 : plan%n = n
223 163476 : plan%m = m
224 163476 : !$ plan%need_alt_plan = .FALSE.
225 :
226 : ! Planning only needed for FFTW3
227 163476 : IF ((fft_type == 3) .AND. (n*m /= 0)) THEN
228 163260 : CALL fftw3_create_plan_1dm(plan, zin, zout, plan_style)
229 163260 : plan%valid = .TRUE.
230 : ELSE
231 216 : plan%valid = .FALSE.
232 : END IF
233 :
234 163476 : END SUBROUTINE fft_create_plan_1dm
235 :
236 : ! **************************************************************************************************
237 : !> \brief ...
238 : !> \param plan ...
239 : ! **************************************************************************************************
240 365772 : SUBROUTINE fft_destroy_plan(plan)
241 : TYPE(fft_plan_type), INTENT(INOUT) :: plan
242 :
243 : ! Planning only needed for FFTW3
244 :
245 365772 : IF (plan%valid .AND. plan%fft_type == 3) THEN
246 219552 : CALL fftw3_destroy_plan(plan)
247 219552 : plan%valid = .FALSE.
248 : END IF
249 :
250 365772 : END SUBROUTINE fft_destroy_plan
251 :
252 : ! **************************************************************************************************
253 : !> \brief ...
254 : !> \param plan ...
255 : !> \param zin ...
256 : !> \param zout ...
257 : !> \param scale ...
258 : !> \param stat ...
259 : ! **************************************************************************************************
260 8458050 : SUBROUTINE fft_1dm(plan, zin, zout, scale, stat)
261 : TYPE(fft_plan_type), INTENT(IN) :: plan
262 : COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zin, zout
263 : REAL(KIND=dp), INTENT(IN) :: scale
264 : INTEGER, INTENT(OUT) :: stat
265 :
266 8458050 : stat = plan%fsign
267 8458050 : IF (plan%n*plan%m > 0) THEN
268 8458050 : SELECT CASE (plan%fft_type)
269 : CASE DEFAULT
270 0 : CPABORT("fft_1dm")
271 : CASE (1)
272 16458 : CALL fftsg1dm(stat, plan%trans, plan%n, plan%m, zin, zout, scale)
273 : CASE (3)
274 8458050 : CALL fftw31dm(plan, zin, zout, scale, stat)
275 : END SELECT
276 : END IF
277 : ! stat is set to zero on error, -1,+1 are OK
278 8458050 : IF (stat == 0) THEN
279 0 : stat = 1
280 : ELSE
281 8458050 : stat = 0
282 : END IF
283 :
284 8458050 : END SUBROUTINE fft_1dm
285 :
286 : END MODULE fft_lib
287 :
|
