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 11685 : 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 14 : flib = -1
51 : CASE ("FFTSG")
52 14 : flib = 1
53 : CASE ("FFTW3")
54 11685 : flib = 3
55 : END SELECT
56 :
57 11685 : 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 11603 : 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 11603 : fft_type = fft_library(fftlib)
81 0 : SELECT CASE (fft_type)
82 : CASE DEFAULT
83 0 : CPABORT("fft_do_init")
84 : CASE (1)
85 14 : CALL fftsg_do_init()
86 : CASE (3)
87 11603 : CALL fftw3_do_init(wisdom_file, plan_style)
88 : END SELECT
89 :
90 11603 : END SUBROUTINE fft_do_init
91 :
92 : ! **************************************************************************************************
93 : !> \brief ...
94 : !> \param wisdom_file ...
95 : !> \param ionode ...
96 : ! **************************************************************************************************
97 11393 : SUBROUTINE fft_do_cleanup(wisdom_file, ionode)
98 : CHARACTER(LEN=*), INTENT(IN) :: wisdom_file
99 : LOGICAL, INTENT(IN) :: ionode
100 :
101 0 : SELECT CASE (fft_type)
102 : CASE DEFAULT
103 0 : CPABORT("fft_do_cleanup")
104 : CASE (1)
105 14 : CALL fftsg_do_cleanup()
106 : CASE (3)
107 11393 : CALL fftw3_do_cleanup(wisdom_file, ionode)
108 : END SELECT
109 :
110 11393 : END SUBROUTINE fft_do_cleanup
111 :
112 : ! **************************************************************************************************
113 : !> \brief ...
114 : !> \param DATA ...
115 : !> \param max_length ...
116 : ! **************************************************************************************************
117 0 : SUBROUTINE fft_get_lengths(DATA, max_length)
118 : INTEGER, DIMENSION(*) :: DATA
119 : INTEGER, INTENT(INOUT) :: max_length
120 :
121 0 : CALL fftw3_get_lengths(DATA, max_length)
122 :
123 0 : END SUBROUTINE fft_get_lengths
124 :
125 : ! **************************************************************************************************
126 :
127 : ! **************************************************************************************************
128 : !> \brief ...
129 : !> \param plan ...
130 : !> \param fft_in_place ...
131 : !> \param fsign ...
132 : !> \param n ...
133 : !> \param zin ...
134 : !> \param zout ...
135 : ! **************************************************************************************************
136 63484 : SUBROUTINE fft_create_plan_3d(plan, fft_in_place, fsign, n, zin, zout)
137 :
138 : TYPE(fft_plan_type), INTENT(INOUT) :: plan
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 :
144 63484 : plan%fsign = fsign
145 63484 : plan%fft_in_place = fft_in_place
146 253936 : plan%n_3d = n
147 63484 : !$ plan%need_alt_plan = .FALSE.
148 :
149 : ! Planning only needed for FFTW3
150 63484 : IF (fft_type == 3) THEN
151 63332 : CALL fftw3_create_plan_3d(plan, zin, zout)
152 63332 : plan%valid = .TRUE.
153 : END IF
154 :
155 63484 : END SUBROUTINE fft_create_plan_3d
156 :
157 : !
158 : ! really ugly, plan is intent out, because plan%fsign is also a status flag
159 : ! if something goes wrong, plan%fsign is set to zero, and the plan becomes invalid
160 : !
161 : ! **************************************************************************************************
162 : !> \brief ...
163 : !> \param plan ...
164 : !> \param scale ...
165 : !> \param zin ...
166 : !> \param zout ...
167 : !> \param stat ...
168 : ! **************************************************************************************************
169 678145 : SUBROUTINE fft_3d(plan, scale, zin, zout, stat)
170 : TYPE(fft_plan_type), INTENT(IN) :: plan
171 : REAL(KIND=dp), INTENT(IN) :: scale
172 : COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zin, zout
173 : INTEGER, INTENT(OUT) :: stat
174 :
175 678145 : stat = plan%fsign
176 678145 : IF (plan%n_3d(1)*plan%n_3d(2)*plan%n_3d(3) > 0) THEN
177 0 : SELECT CASE (fft_type)
178 : CASE DEFAULT
179 0 : CPABORT("fft_3d")
180 : CASE (1)
181 3008 : CALL fftsg3d(plan%fft_in_place, stat, scale, plan%n_3d, zin, zout)
182 : CASE (3)
183 678145 : CALL fftw33d(plan, scale, zin, zout, stat)
184 : END SELECT
185 : END IF
186 : ! stat is set to zero on error, -1,+1 are OK
187 678145 : IF (stat == 0) THEN
188 0 : stat = 1
189 : ELSE
190 678145 : stat = 0
191 : END IF
192 :
193 678145 : END SUBROUTINE fft_3d
194 :
195 : ! **************************************************************************************************
196 :
197 : ! **************************************************************************************************
198 : !> \brief ...
199 : !> \param plan ...
200 : !> \param fsign ...
201 : !> \param trans_in ...
202 : !> \param trans_out ...
203 : !> \param ldx_in ...
204 : !> \param ldx_out ...
205 : !> \param n ...
206 : !> \param m ...
207 : !> \param zin ...
208 : !> \param zout ...
209 : ! **************************************************************************************************
210 469166 : SUBROUTINE fft_create_plan_1d(plan, fsign, trans_in, trans_out, ldx_in, ldx_out, n, m, zin, zout)
211 : TYPE(fft_plan_type), INTENT(INOUT) :: plan
212 : INTEGER, INTENT(IN) :: fsign
213 : LOGICAL, INTENT(IN) :: trans_in, trans_out
214 : INTEGER, INTENT(IN) :: ldx_in, ldx_out, n, m
215 : COMPLEX(KIND=dp), DIMENSION(*), INTENT(IN) :: zin, zout
216 :
217 469166 : plan%fsign = fsign
218 469166 : plan%trans_in = trans_in
219 469166 : plan%trans_out = trans_out
220 469166 : IF (plan%trans_in) THEN
221 362999 : plan%ldx_in = ldx_in
222 362999 : plan%ldy_in = n
223 : ELSE
224 106167 : plan%ldx_in = ldx_in
225 106167 : plan%ldy_in = m
226 : END IF
227 469166 : IF (plan%trans_out) THEN
228 362999 : plan%ldx_out = ldx_out
229 362999 : plan%ldy_out = n
230 : ELSE
231 106167 : plan%ldx_out = ldx_out
232 106167 : plan%ldy_out = m
233 : END IF
234 :
235 469166 : plan%n = n
236 469166 : plan%m = m
237 469166 : !$ plan%need_alt_plan = .FALSE.
238 :
239 : ! Planning only needed for FFTW3
240 469166 : IF ((fft_type == 3) .AND. (n*m /= 0)) THEN
241 467762 : CALL fftw3_create_plan_1d(plan, zin, zout)
242 467762 : plan%valid = .TRUE.
243 : ELSE
244 1404 : plan%valid = .FALSE.
245 : END IF
246 :
247 469166 : END SUBROUTINE fft_create_plan_1d
248 :
249 : ! **************************************************************************************************
250 : !> \brief ...
251 : !> \param plan ...
252 : ! **************************************************************************************************
253 714926 : SUBROUTINE fft_destroy_plan(plan)
254 : TYPE(fft_plan_type), INTENT(INOUT) :: plan
255 :
256 : ! Planning only needed for FFTW3
257 :
258 714926 : IF (plan%valid .AND. fft_type == 3) THEN
259 531094 : CALL fftw3_destroy_plan(plan)
260 531094 : plan%valid = .FALSE.
261 : END IF
262 :
263 714926 : END SUBROUTINE fft_destroy_plan
264 :
265 : ! **************************************************************************************************
266 : !> \brief ...
267 : !> \param plan ...
268 : !> \param zin ...
269 : !> \param zout ...
270 : !> \param scale ...
271 : !> \param stat ...
272 : ! **************************************************************************************************
273 20388664 : SUBROUTINE fft_1d(plan, zin, zout, scale, stat)
274 : TYPE(fft_plan_type), INTENT(IN) :: plan
275 : COMPLEX(KIND=dp), DIMENSION(*), INTENT(INOUT) :: zin, zout
276 : REAL(KIND=dp), INTENT(IN) :: scale
277 : INTEGER, INTENT(OUT) :: stat
278 :
279 20388664 : stat = plan%fsign
280 20388664 : IF (plan%n*plan%m > 0) THEN
281 0 : SELECT CASE (fft_type)
282 : CASE DEFAULT
283 0 : CPABORT("fft_1d")
284 : CASE (1)
285 : CALL fftsg1dm(stat, plan%trans_in, plan%trans_out, plan%n, plan%m, &
286 30346 : plan%ldx_in, plan%ldy_in, plan%ldx_out, plan%ldy_out, zin, zout, scale)
287 : CASE (3)
288 20388664 : CALL fftw31dm(plan, zin, zout, scale, stat)
289 : END SELECT
290 : END IF
291 : ! stat is set to zero on error, -1,+1 are OK
292 20388664 : IF (stat == 0) THEN
293 0 : stat = 1
294 : ELSE
295 20388664 : stat = 0
296 : END IF
297 :
298 20388664 : END SUBROUTINE fft_1d
299 :
300 : END MODULE fft_lib
|
