Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 : USE fftw3_lib, ONLY: &
16 : fft_alloc => fftw_alloc, fft_dealloc => fftw_dealloc, fftw31dm, fftw33d, &
17 : fftw3_create_plan_1d, fftw3_create_plan_3d, fftw3_destroy_plan, fftw3_do_cleanup, &
18 : fftw3_do_init, fftw3_get_lengths
19 : #include "../../base/base_uses.f90"
20 :
21 : IMPLICIT NONE
22 : PRIVATE
23 :
24 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fft_lib'
25 :
26 : INTEGER, SAVE :: fft_type = 0
27 :
28 : PUBLIC :: fft_do_cleanup, fft_do_init, fft_get_lengths, fft_create_plan_3d
29 : PUBLIC :: fft_create_plan_1d, fft_1d, fft_library, fft_3d, fft_destroy_plan
30 : PUBLIC :: fft_supports_arbitrary_lengths
31 : PUBLIC :: fft_alloc, fft_dealloc
32 :
33 : CONTAINS
34 : ! **************************************************************************************************
35 : !> \brief Interface to FFT libraries
36 : !> \param fftlib ...
37 : !> \return ...
38 : !> \par History
39 : !> IAB 09-Jan-2009 : Modified to use fft_plan_type
40 : !> (c) The Numerical Algorithms Group (NAG) Ltd, 2009 on behalf of the HECToR project
41 : !> \author JGH
42 : ! **************************************************************************************************
43 11814 : FUNCTION fft_library(fftlib) RESULT(flib)
44 :
45 : CHARACTER(len=*), INTENT(IN) :: fftlib
46 : INTEGER :: flib
47 :
48 : SELECT CASE (fftlib)
49 : CASE DEFAULT
50 16 : flib = -1
51 : CASE ("FFTSG")
52 16 : flib = 1
53 : CASE ("FFTW3")
54 11814 : flib = 3
55 : END SELECT
56 :
57 11814 : END FUNCTION fft_library
58 :
59 : ! **************************************************************************************************
60 : !> \brief Query whether the active FFT backend supports arbitrary transform lengths.
61 : !> \return whether arbitrary positive transform lengths are supported
62 : ! **************************************************************************************************
63 82 : LOGICAL FUNCTION fft_supports_arbitrary_lengths() RESULT(supported)
64 :
65 82 : supported = fft_type == fft_library("FFTW3")
66 :
67 82 : END FUNCTION fft_supports_arbitrary_lengths
68 :
69 : ! **************************************************************************************************
70 : !> \brief ...
71 : !> \param fftlib ...
72 : !> \param plan_style ...
73 : !> \param wisdom_file ...
74 : ! **************************************************************************************************
75 11732 : SUBROUTINE fft_do_init(fftlib, plan_style, wisdom_file)
76 : CHARACTER(LEN=*), INTENT(IN) :: fftlib
77 : INTEGER, INTENT(IN) :: plan_style
78 : CHARACTER(LEN=*), INTENT(IN) :: wisdom_file
79 :
80 11732 : fft_type = fft_library(fftlib)
81 0 : SELECT CASE (fft_type)
82 : CASE DEFAULT
83 0 : CPABORT("fft_do_init")
84 : CASE (1)
85 16 : CALL fftsg_do_init()
86 : CASE (3)
87 11732 : CALL fftw3_do_init(wisdom_file, plan_style)
88 : END SELECT
89 :
90 11732 : END SUBROUTINE fft_do_init
91 :
92 : ! **************************************************************************************************
93 : !> \brief ...
94 : !> \param wisdom_file ...
95 : !> \param ionode ...
96 : !> \param export_wisdom Export the FFTW wisdom to wisdom_file (default: no export)
97 : ! **************************************************************************************************
98 11521 : SUBROUTINE fft_do_cleanup(wisdom_file, ionode, export_wisdom)
99 : CHARACTER(LEN=*), INTENT(IN) :: wisdom_file
100 : LOGICAL, INTENT(IN) :: ionode
101 : LOGICAL, INTENT(IN), OPTIONAL :: export_wisdom
102 :
103 0 : SELECT CASE (fft_type)
104 : CASE DEFAULT
105 0 : CPABORT("fft_do_cleanup")
106 : CASE (1)
107 16 : CALL fftsg_do_cleanup()
108 : CASE (3)
109 11521 : CALL fftw3_do_cleanup(wisdom_file, ionode, export_wisdom)
110 : END SELECT
111 :
112 11521 : END SUBROUTINE fft_do_cleanup
113 :
114 : ! **************************************************************************************************
115 : !> \brief ...
116 : !> \param DATA ...
117 : !> \param max_length ...
118 : ! **************************************************************************************************
119 0 : SUBROUTINE fft_get_lengths(DATA, max_length)
120 : INTEGER, DIMENSION(*) :: DATA
121 : INTEGER, INTENT(INOUT) :: max_length
122 :
123 0 : CALL fftw3_get_lengths(DATA, max_length)
124 :
125 0 : END SUBROUTINE fft_get_lengths
126 :
127 : ! **************************************************************************************************
128 :
129 : ! **************************************************************************************************
130 : !> \brief ...
131 : !> \param plan ...
132 : !> \param fft_in_place ...
133 : !> \param fsign ...
134 : !> \param n ...
135 : !> \param zin ...
136 : !> \param zout ...
137 : ! **************************************************************************************************
138 63976 : SUBROUTINE fft_create_plan_3d(plan, fft_in_place, fsign, n, zin, zout)
139 :
140 : TYPE(fft_plan_type), INTENT(INOUT) :: plan
141 : LOGICAL, INTENT(IN) :: fft_in_place
142 : INTEGER, INTENT(IN) :: fsign
143 : INTEGER, DIMENSION(3), INTENT(IN) :: n
144 : COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zin, zout
145 :
146 63976 : plan%fsign = fsign
147 63976 : plan%fft_in_place = fft_in_place
148 255904 : plan%n_3d = n
149 63976 : !$ plan%need_alt_plan = .FALSE.
150 :
151 : ! Planning only needed for FFTW3
152 63976 : IF (fft_type == 3) THEN
153 63792 : CALL fftw3_create_plan_3d(plan, zin, zout)
154 63792 : plan%valid = .TRUE.
155 : END IF
156 :
157 63976 : 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 698574 : 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 698574 : stat = plan%fsign
178 698574 : IF (plan%n_3d(1)*plan%n_3d(2)*plan%n_3d(3) > 0) THEN
179 0 : SELECT CASE (fft_type)
180 : CASE DEFAULT
181 0 : CPABORT("fft_3d")
182 : CASE (1)
183 6344 : CALL fftsg3d(plan%fft_in_place, stat, scale, plan%n_3d, zin, zout)
184 : CASE (3)
185 698574 : 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 698574 : IF (stat == 0) THEN
190 0 : stat = 1
191 : ELSE
192 698574 : stat = 0
193 : END IF
194 :
195 698574 : END SUBROUTINE fft_3d
196 :
197 : ! **************************************************************************************************
198 :
199 : ! **************************************************************************************************
200 : !> \brief ...
201 : !> \param plan ...
202 : !> \param fsign ...
203 : !> \param trans_in ...
204 : !> \param trans_out ...
205 : !> \param ldx_in ...
206 : !> \param ldx_out ...
207 : !> \param n ...
208 : !> \param m ...
209 : !> \param zin ...
210 : !> \param zout ...
211 : ! **************************************************************************************************
212 473796 : SUBROUTINE fft_create_plan_1d(plan, fsign, trans_in, trans_out, ldx_in, ldx_out, n, m, zin, zout)
213 : TYPE(fft_plan_type), INTENT(INOUT) :: plan
214 : INTEGER, INTENT(IN) :: fsign
215 : LOGICAL, INTENT(IN) :: trans_in, trans_out
216 : INTEGER, INTENT(IN) :: ldx_in, ldx_out, n, m
217 : COMPLEX(KIND=dp), DIMENSION(*), INTENT(IN) :: zin, zout
218 :
219 473796 : plan%fsign = fsign
220 473796 : plan%trans_in = trans_in
221 473796 : plan%trans_out = trans_out
222 473796 : IF (plan%trans_in) THEN
223 366900 : plan%ldx_in = ldx_in
224 366900 : plan%ldy_in = n
225 : ELSE
226 106896 : plan%ldx_in = ldx_in
227 106896 : plan%ldy_in = m
228 : END IF
229 473796 : IF (plan%trans_out) THEN
230 366900 : plan%ldx_out = ldx_out
231 366900 : plan%ldy_out = n
232 : ELSE
233 106896 : plan%ldx_out = ldx_out
234 106896 : plan%ldy_out = m
235 : END IF
236 :
237 473796 : plan%n = n
238 473796 : plan%m = m
239 473796 : !$ plan%need_alt_plan = .FALSE.
240 :
241 : ! Planning only needed for FFTW3
242 473796 : IF ((fft_type == 3) .AND. (n*m /= 0)) THEN
243 472392 : CALL fftw3_create_plan_1d(plan, zin, zout)
244 472392 : plan%valid = .TRUE.
245 : ELSE
246 1404 : plan%valid = .FALSE.
247 : END IF
248 :
249 473796 : END SUBROUTINE fft_create_plan_1d
250 :
251 : ! **************************************************************************************************
252 : !> \brief ...
253 : !> \param plan ...
254 : ! **************************************************************************************************
255 721818 : SUBROUTINE fft_destroy_plan(plan)
256 : TYPE(fft_plan_type), INTENT(INOUT) :: plan
257 :
258 : ! Planning only needed for FFTW3
259 :
260 721818 : IF (plan%valid .AND. fft_type == 3) THEN
261 536154 : CALL fftw3_destroy_plan(plan)
262 536154 : plan%valid = .FALSE.
263 : END IF
264 :
265 721818 : END SUBROUTINE fft_destroy_plan
266 :
267 : ! **************************************************************************************************
268 : !> \brief ...
269 : !> \param plan ...
270 : !> \param zin ...
271 : !> \param zout ...
272 : !> \param scale ...
273 : !> \param stat ...
274 : ! **************************************************************************************************
275 20876585 : SUBROUTINE fft_1d(plan, zin, zout, scale, stat)
276 : TYPE(fft_plan_type), INTENT(IN) :: plan
277 : COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zin, zout
278 : REAL(KIND=dp), INTENT(IN) :: scale
279 : INTEGER, INTENT(OUT) :: stat
280 :
281 20876585 : stat = plan%fsign
282 20876585 : IF (plan%n*plan%m > 0) THEN
283 0 : SELECT CASE (fft_type)
284 : CASE DEFAULT
285 0 : CPABORT("fft_1d")
286 : CASE (1)
287 : CALL fftsg1dm(stat, plan%trans_in, plan%trans_out, plan%n, plan%m, &
288 30946 : plan%ldx_in, plan%ldy_in, plan%ldx_out, plan%ldy_out, zin, zout, scale)
289 : CASE (3)
290 20876585 : CALL fftw31dm(plan, zin, zout, scale, stat)
291 : END SELECT
292 : END IF
293 : ! stat is set to zero on error, -1,+1 are OK
294 20876585 : IF (stat == 0) THEN
295 0 : stat = 1
296 : ELSE
297 20876585 : stat = 0
298 : END IF
299 :
300 20876585 : END SUBROUTINE fft_1d
301 :
302 : END MODULE fft_lib
|
