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 :
8 : ! **************************************************************************************************
9 : !> \brief Creates the wavelet kernel for the wavelet based poisson solver.
10 : !> \author Florian Schiffmann (09.2007,fschiff)
11 : ! **************************************************************************************************
12 : MODULE ps_wavelet_base
13 : USE fft_lib, ONLY: fft_1d,&
14 : fft_create_plan_1d,&
15 : fft_destroy_plan
16 : USE fft_plan, ONLY: fft_plan_type
17 : USE fft_tools, ONLY: BWFFT,&
18 : FWFFT,&
19 : fft_alloc,&
20 : fft_dealloc
21 : USE kinds, ONLY: dp
22 : USE mathconstants, ONLY: pi
23 : USE message_passing, ONLY: mp_comm_type
24 : #include "../base/base_uses.f90"
25 :
26 : IMPLICIT NONE
27 :
28 : PRIVATE
29 :
30 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ps_wavelet_base'
31 :
32 : PUBLIC :: scramble_unpack, p_poissonsolver, s_poissonsolver, f_poissonsolver
33 :
34 : CONTAINS
35 :
36 : ! **************************************************************************************************
37 : !> \brief ...
38 : !> \param n1 ...
39 : !> \param n2 ...
40 : !> \param n3 ...
41 : !> \param nd1 ...
42 : !> \param nd2 ...
43 : !> \param nd3 ...
44 : !> \param md1 ...
45 : !> \param md2 ...
46 : !> \param md3 ...
47 : !> \param nproc ...
48 : !> \param iproc ...
49 : !> \param zf ...
50 : !> \param scal ...
51 : !> \param hx ...
52 : !> \param hy ...
53 : !> \param hz ...
54 : !> \param mpi_group ...
55 : ! **************************************************************************************************
56 17323 : SUBROUTINE P_PoissonSolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, zf &
57 : , scal, hx, hy, hz, mpi_group)
58 : INTEGER, INTENT(in) :: n1, n2, n3, nd1, nd2, nd3, md1, md2, &
59 : md3, nproc, iproc
60 : REAL(KIND=dp), DIMENSION(md1, md3, md2/nproc), &
61 : INTENT(inout) :: zf
62 : REAL(KIND=dp), INTENT(in) :: scal, hx, hy, hz
63 :
64 : CLASS(mp_comm_type), INTENT(in) :: mpi_group
65 :
66 : INTEGER, PARAMETER :: ncache_optimal = 8*1024
67 :
68 : INTEGER :: i1, i3, j, j2, &
69 : J2stb, J2stf, j3, Jp2stb, Jp2stf, lot1, lot2, lot3, &
70 : lzt, ma, mb, ncache, nfft, stat, &
71 : final_chunk_size3, final_chunk_size1, final_chunk_size2
72 17323 : COMPLEX(KIND=dp), POINTER, CONTIGUOUS, DIMENSION(:, :) :: zt
73 17323 : COMPLEX(KIND=dp), POINTER, CONTIGUOUS, DIMENSION(:) :: zw1, zw2
74 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: zmpi2
75 : COMPLEX(KIND=dp), ALLOCATABLE, &
76 17323 : DIMENSION(:, :, :, :) :: zmpi1
77 : TYPE(fft_plan_type) :: fft_plan_bw3, fft_plan_bw3_last, fft_plan_fw3, fft_plan_fw3_last, &
78 : fft_plan_bw1, fft_plan_bw1_last, fft_plan_fw1, fft_plan_fw1_last, &
79 : fft_plan_bw2, fft_plan_bw2_last, fft_plan_fw2, fft_plan_fw2_last
80 :
81 0 : IF (nd1 < n1/2 + 1) CPABORT("Parallel convolution:ERROR:nd1")
82 17323 : IF (nd2 < n2/2 + 1) CPABORT("Parallel convolution:ERROR:nd2")
83 17323 : IF (nd3 < n3/2 + 1) CPABORT("Parallel convolution:ERROR:nd3")
84 17323 : IF (md1 < n1) CPABORT("Parallel convolution:ERROR:md1")
85 17323 : IF (md2 < n2) CPABORT("Parallel convolution:ERROR:md2")
86 17323 : IF (md3 < n3) CPABORT("Parallel convolution:ERROR:md3")
87 17323 : IF (MOD(nd3, nproc) /= 0) CPABORT("Parallel convolution:ERROR:nd3")
88 17323 : IF (MOD(md2, nproc) /= 0) CPABORT("Parallel convolution:ERROR:md2")
89 :
90 : !defining work arrays dimensions
91 17323 : ncache = ncache_optimal
92 17323 : IF (ncache <= MAX(n1, n2, n3)*4) ncache = MAX(n1, n2, n3)*4
93 :
94 17323 : lzt = n2
95 17323 : IF (MOD(n2, 2) == 0) lzt = lzt + 1
96 :
97 : !Allocations
98 34646 : CALL fft_alloc(zw1, [ncache/4])
99 17323 : zw1 = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
100 34646 : CALL fft_alloc(zw2, [ncache/4])
101 17323 : zw2 = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
102 51969 : CALL fft_alloc(zt, [lzt, n1])
103 17323 : zt = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
104 86615 : ALLOCATE (zmpi2(n1, md2/nproc, nd3), SOURCE=CMPLX(0.0_dp, 0.0_dp, KIND=dp))
105 41353 : IF (nproc > 1) ALLOCATE (zmpi1(n1, md2/nproc, nd3/nproc, nproc), SOURCE=CMPLX(0.0_dp, 0.0_dp, KIND=dp))
106 :
107 : ! transform along z axis
108 : ! Leading dimensions in all three directions
109 17323 : lot1 = ncache/(4*n1)
110 17323 : lot2 = ncache/(4*n2)
111 17323 : lot3 = ncache/(4*n3)
112 :
113 : ! The size of the last chunk
114 17323 : final_chunk_size1 = MOD(n2, lot1)
115 17323 : final_chunk_size2 = MOD(n1, lot2)
116 17323 : final_chunk_size3 = MOD(n1, lot3)
117 :
118 : ! Prevent OOB-access if n2 < lot1
119 17323 : IF (n2 >= lot1) THEN
120 251 : CALL fft_create_plan_1d(fft_plan_bw1, BWFFT, .TRUE., .TRUE., lot1, lzt, n1, lot1, zw1, zt)
121 251 : CALL fft_create_plan_1d(fft_plan_fw1, FWFFT, .TRUE., .TRUE., lzt, lot1, n1, lot1, zt, zw1)
122 : END IF
123 17323 : IF (final_chunk_size1 > 0) THEN
124 : CALL fft_create_plan_1d(fft_plan_bw1_last, BWFFT, .TRUE., .TRUE., lot1, lzt, n1, &
125 17323 : final_chunk_size1, zw1, zt)
126 : CALL fft_create_plan_1d(fft_plan_fw1_last, FWFFT, .TRUE., .TRUE., lzt, lot1, n1, &
127 17323 : final_chunk_size1, zt, zw1)
128 : END IF
129 :
130 : ! Prevent OOB-access if n2 < lot1
131 17323 : IF (n1 >= lot2) THEN
132 251 : CALL fft_create_plan_1d(fft_plan_bw2, BWFFT, .TRUE., .TRUE., lot2, lot2, n2, lot2, zw1, zw2)
133 251 : CALL fft_create_plan_1d(fft_plan_fw2, FWFFT, .TRUE., .TRUE., lot2, lot2, n2, lot2, zw2, zw1)
134 : END IF
135 17323 : IF (final_chunk_size2 > 0) THEN
136 : CALL fft_create_plan_1d(fft_plan_bw2_last, BWFFT, .TRUE., .TRUE., lot2, lot2, n2, &
137 17323 : final_chunk_size2, zw1, zw2)
138 : CALL fft_create_plan_1d(fft_plan_fw2_last, FWFFT, .TRUE., .TRUE., lot2, lot2, n2, &
139 17323 : final_chunk_size2, zw2, zw1)
140 : END IF
141 :
142 : ! Prevent OOB-access if n1 < lot3
143 17323 : IF (n1 >= lot3) THEN
144 251 : CALL fft_create_plan_1d(fft_plan_fw3, FWFFT, .TRUE., .TRUE., lot3, lot3, n3, lot3, zw1, zw2)
145 251 : CALL fft_create_plan_1d(fft_plan_bw3, BWFFT, .TRUE., .TRUE., lot3, lot3, n3, lot3, zw1, zw2)
146 : END IF
147 17323 : IF (final_chunk_size3 > 0) THEN
148 : CALL fft_create_plan_1d(fft_plan_fw3_last, FWFFT, .TRUE., .TRUE., &
149 17323 : lot3, lot3, n3, final_chunk_size3, zw1, zw2)
150 : CALL fft_create_plan_1d(fft_plan_bw3_last, BWFFT, .TRUE., .TRUE., &
151 17323 : lot3, lot3, n3, final_chunk_size3, zw1, zw2)
152 : END IF
153 :
154 320326 : DO j2 = 1, md2/nproc
155 : !this condition ensures that we manage only the interesting part for the FFT
156 320326 : IF (iproc*(md2/nproc) + j2 <= n2) THEN
157 613678 : DO i1 = 1, n1, lot3
158 311348 : ma = i1
159 311348 : mb = MIN(i1 + (lot3 - 1), n1)
160 311348 : nfft = mb - ma + 1
161 : !inserting real data into complex array of half length
162 311348 : CALL P_fill_upcorn(md1, md3, lot3, nfft, n3, zf(i1, 1, j2), zw1)
163 :
164 : !performing FFT
165 : !input: I1,I3,J2,(Jp2)
166 : ! This is equivalent to
167 311348 : IF (nfft == lot3) THEN
168 9018 : CALL fft_1d(fft_plan_bw3, zw1, zw2, 1.0_dp, stat)
169 : ELSE
170 302330 : CALL fft_1d(fft_plan_bw3_last, zw1, zw2, 1.0_dp, stat)
171 : END IF
172 :
173 : !output: I1,i3,J2,(Jp2)
174 : !exchanging components
175 : !input: I1,i3,J2,(Jp2)
176 613678 : CALL scramble_P(i1, j2, lot3, nfft, n1, n3, md2, nproc, nd3, zw2, zmpi2)
177 : !output: I1,J2,i3,(Jp2)
178 : END DO
179 : END IF
180 : END DO
181 :
182 : !Interprocessor data transposition
183 : !input: I1,J2,j3,jp3,(Jp2)
184 17323 : IF (nproc > 1) THEN
185 : !communication scheduling
186 4806 : CALL mpi_group%alltoall(zmpi2, zmpi1, n1*(md2/nproc)*(nd3/nproc))
187 : END IF
188 : !output: I1,J2,j3,Jp2,(jp3)
189 :
190 : !now each process perform complete convolution of its planes
191 182731 : DO j3 = 1, nd3/nproc
192 : !this condition ensures that we manage only the interesting part for the FFT
193 182731 : IF (iproc*(nd3/nproc) + j3 <= n3/2 + 1) THEN
194 164470 : Jp2stb = 1
195 164470 : J2stb = 1
196 164470 : Jp2stf = 1
197 164470 : J2stf = 1
198 :
199 : ! transform along x axis
200 :
201 333616 : DO j = 1, n2, lot1
202 169146 : ma = j
203 169146 : mb = MIN(j + (lot1 - 1), n2)
204 169146 : nfft = mb - ma + 1
205 :
206 : !reverse index ordering, leaving the planes to be transformed at the end
207 : !input: I1,J2,j3,Jp2,(jp3)
208 169146 : IF (nproc == 1) THEN
209 136958 : CALL P_mpiswitch_upcorn(j3, nfft, Jp2stb, J2stb, lot1, n1, md2, nd3, nproc, zmpi2, zw1)
210 : ELSE
211 32188 : CALL P_mpiswitch_upcorn(j3, nfft, Jp2stb, J2stb, lot1, n1, md2, nd3, nproc, zmpi1, zw1)
212 : END IF
213 : !output: J2,Jp2,I1,j3,(jp3)
214 :
215 : !performing FFT
216 : !input: I2,I1,j3,(jp3)
217 : ! This is equivalent to
218 333616 : IF (nfft == lot1) THEN
219 4676 : CALL fft_1d(fft_plan_bw1, zw1, zt(j:, 1), 1.0_dp, stat)
220 : ELSE
221 164470 : CALL fft_1d(fft_plan_bw1_last, zw1, zt(j:, 1), 1.0_dp, stat)
222 : END IF
223 : !output: I2,i1,j3,(jp3)
224 : END DO
225 :
226 : !transform along y axis
227 :
228 333616 : DO j = 1, n1, lot2
229 169146 : ma = j
230 169146 : mb = MIN(j + (lot2 - 1), n1)
231 169146 : nfft = mb - ma + 1
232 :
233 : !reverse ordering
234 : !input: I2,i1,j3,(jp3)
235 169146 : CALL P_switch_upcorn(nfft, n2, lot2, n1, lzt, zt(:, j), zw1)
236 : !output: i1,I2,j3,(jp3)
237 :
238 : !performing FFT
239 : !input: i1,I2,j3,(jp3)
240 : ! This is equivalent to
241 169146 : IF (nfft == lot2) THEN
242 4676 : CALL fft_1d(fft_plan_bw2, zw1, zw2, 1.0_dp, stat)
243 : ELSE
244 164470 : CALL fft_1d(fft_plan_bw2_last, zw1, zw2, 1.0_dp, stat)
245 : END IF
246 : !output: i1,i2,j3,(jp3)
247 :
248 : !Multiply with kernel in fourier space
249 169146 : i3 = iproc*(nd3/nproc) + j3
250 169146 : CALL P_multkernel(n1, n2, n3, lot2, nfft, j, i3, zw2, hx, hy, hz)
251 :
252 : !TRANSFORM BACK IN REAL SPACE
253 :
254 : !transform along y axis
255 : !input: i1,i2,j3,(jp3)
256 : ! This is equivalent to
257 169146 : IF (nfft == lot2) THEN
258 4676 : CALL fft_1d(fft_plan_fw2, zw2, zw1, 1.0_dp, stat)
259 : ELSE
260 164470 : CALL fft_1d(fft_plan_fw2_last, zw2, zw1, 1.0_dp, stat)
261 : END IF
262 :
263 : !reverse ordering
264 : !input: i1,I2,j3,(jp3)
265 333616 : CALL P_unswitch_downcorn(nfft, n2, lot2, n1, lzt, zw1, zt(:, j))
266 : !output: I2,i1,j3,(jp3)
267 : END DO
268 :
269 : !transform along x axis
270 : !input: I2,i1,j3,(jp3)
271 333616 : DO j = 1, n2, lot1
272 169146 : ma = j
273 169146 : mb = MIN(j + (lot1 - 1), n2)
274 169146 : nfft = mb - ma + 1
275 :
276 : !performing FFT
277 : ! This is equivalent to
278 169146 : IF (nfft == lot1) THEN
279 4676 : CALL fft_1d(fft_plan_fw1, zt(j:, 1), zw2, 1.0_dp, stat)
280 : ELSE
281 164470 : CALL fft_1d(fft_plan_fw1_last, zt(j:, 1), zw2, 1.0_dp, stat)
282 : END IF
283 : !output: I2,I1,j3,(jp3)
284 :
285 : !reverse ordering
286 : !input: J2,Jp2,I1,j3,(jp3)
287 333616 : IF (nproc == 1) THEN
288 136958 : CALL P_unmpiswitch_downcorn(j3, nfft, Jp2stf, J2stf, lot1, n1, md2, nd3, nproc, zw2, zmpi2)
289 : ELSE
290 32188 : CALL P_unmpiswitch_downcorn(j3, nfft, Jp2stf, J2stf, lot1, n1, md2, nd3, nproc, zw2, zmpi1)
291 : END IF
292 : ! output: I1,J2,j3,Jp2,(jp3)
293 : END DO
294 : END IF
295 : END DO
296 :
297 : !Interprocessor data transposition
298 : !input: I1,J2,j3,Jp2,(jp3)
299 17323 : IF (nproc > 1) THEN
300 : !communication scheduling
301 4806 : CALL mpi_group%alltoall(zmpi1, zmpi2, n1*(md2/nproc)*(nd3/nproc))
302 : END IF
303 : !output: I1,J2,j3,jp3,(Jp2)
304 : !transform along z axis
305 : !input: I1,J2,i3,(Jp2)
306 320326 : DO j2 = 1, md2/nproc
307 : !this condition ensures that we manage only the interesting part for the FFT
308 320326 : IF (iproc*(md2/nproc) + j2 <= n2) THEN
309 613678 : DO i1 = 1, n1, lot3
310 311348 : ma = i1
311 311348 : mb = MIN(i1 + (lot3 - 1), n1)
312 311348 : nfft = mb - ma + 1
313 :
314 : !reverse ordering
315 : !input: I1,J2,i3,(Jp2)
316 311348 : CALL unscramble_P(i1, j2, lot3, nfft, n1, n3, md2, nproc, nd3, zmpi2, zw1)
317 : !output: I1,i3,J2,(Jp2)
318 :
319 : !performing FFT
320 : !input: I1,i3,J2,(Jp2)
321 : ! This is equivalent to
322 311348 : IF (nfft == lot3) THEN
323 9018 : CALL fft_1d(fft_plan_fw3, zw1, zw2, 1.0_dp, stat)
324 : ELSE
325 302330 : CALL fft_1d(fft_plan_fw3_last, zw1, zw2, 1.0_dp, stat)
326 : END IF
327 : !output: I1,I3,J2,(Jp2)
328 :
329 : !rebuild the output array
330 613678 : CALL P_unfill_downcorn(md1, md3, lot3, nfft, n3, zw2, zf(i1, 1, j2), scal)
331 :
332 : END DO
333 : END IF
334 : END DO
335 :
336 17323 : IF (n2 >= lot1) CALL fft_destroy_plan(fft_plan_bw1)
337 17323 : IF (final_chunk_size1 > 0) CALL fft_destroy_plan(fft_plan_bw1_last)
338 17323 : IF (n2 >= lot1) CALL fft_destroy_plan(fft_plan_fw1)
339 17323 : IF (final_chunk_size1 > 0) CALL fft_destroy_plan(fft_plan_fw1_last)
340 :
341 17323 : IF (n1 >= lot2) CALL fft_destroy_plan(fft_plan_bw2)
342 17323 : IF (final_chunk_size2 > 0) CALL fft_destroy_plan(fft_plan_bw2_last)
343 17323 : IF (n1 >= lot2) CALL fft_destroy_plan(fft_plan_fw2)
344 17323 : IF (final_chunk_size2 > 0) CALL fft_destroy_plan(fft_plan_fw2_last)
345 :
346 17323 : IF (n1 >= lot3) CALL fft_destroy_plan(fft_plan_fw3)
347 17323 : IF (final_chunk_size3 > 0) CALL fft_destroy_plan(fft_plan_fw3_last)
348 17323 : IF (n1 >= lot3) CALL fft_destroy_plan(fft_plan_bw3)
349 17323 : IF (final_chunk_size3 > 0) CALL fft_destroy_plan(fft_plan_bw3_last)
350 :
351 : !De-allocations
352 17323 : DEALLOCATE (zmpi2)
353 17323 : CALL fft_dealloc(zw1)
354 17323 : CALL fft_dealloc(zw2)
355 17323 : CALL fft_dealloc(zt)
356 17323 : IF (nproc > 1) DEALLOCATE (zmpi1)
357 658274 : END SUBROUTINE P_PoissonSolver
358 :
359 : ! **************************************************************************************************
360 : !> \brief ...
361 : !> \param j3 ...
362 : !> \param nfft ...
363 : !> \param Jp2stb ...
364 : !> \param J2stb ...
365 : !> \param lot ...
366 : !> \param n1 ...
367 : !> \param md2 ...
368 : !> \param nd3 ...
369 : !> \param nproc ...
370 : !> \param zmpi1 ...
371 : !> \param zw ...
372 : ! **************************************************************************************************
373 169146 : SUBROUTINE P_mpiswitch_upcorn(j3, nfft, Jp2stb, J2stb, lot, n1, md2, nd3, nproc, zmpi1, zw)
374 : INTEGER, INTENT(in) :: j3, nfft
375 : INTEGER, INTENT(inout) :: Jp2stb, J2stb
376 : INTEGER, INTENT(in) :: lot, n1, md2, nd3, nproc
377 : COMPLEX(KIND=dp), &
378 : DIMENSION(n1, md2/nproc, nd3/nproc, nproc), &
379 : INTENT(in) :: zmpi1
380 : COMPLEX(KIND=dp), DIMENSION(lot, n1), &
381 : INTENT(inout) :: zw
382 :
383 : INTEGER :: I1, J2, Jp2, mfft
384 :
385 169146 : mfft = 0
386 354650 : DO Jp2 = Jp2stb, nproc
387 3840810 : DO J2 = J2stb, md2/nproc
388 3655306 : mfft = mfft + 1
389 3655306 : IF (mfft > nfft) THEN
390 13478 : Jp2stb = Jp2
391 13478 : J2stb = J2
392 13478 : RETURN
393 : END IF
394 95925808 : DO I1 = 1, n1
395 95740304 : zw(mfft, I1) = zmpi1(I1, J2, j3, Jp2)
396 : END DO
397 : END DO
398 341172 : J2stb = 1
399 : END DO
400 : END SUBROUTINE P_mpiswitch_upcorn
401 :
402 : ! **************************************************************************************************
403 : !> \brief ...
404 : !> \param nfft ...
405 : !> \param n2 ...
406 : !> \param lot ...
407 : !> \param n1 ...
408 : !> \param lzt ...
409 : !> \param zt ...
410 : !> \param zw ...
411 : ! **************************************************************************************************
412 169146 : SUBROUTINE P_switch_upcorn(nfft, n2, lot, n1, lzt, zt, zw)
413 : INTEGER, INTENT(in) :: nfft, n2, lot, n1, lzt
414 : COMPLEX(KIND=dp), DIMENSION(lzt, n1), INTENT(in) :: zt
415 : COMPLEX(KIND=dp), DIMENSION(lot, n2), &
416 : INTENT(inout) :: zw
417 :
418 : INTEGER :: i, j
419 :
420 3810974 : DO j = 1, nfft
421 95909450 : DO i = 1, n2
422 95740304 : zw(j, i) = zt(i, j)
423 : END DO
424 : END DO
425 :
426 169146 : END SUBROUTINE P_switch_upcorn
427 :
428 : ! **************************************************************************************************
429 : !> \brief ...
430 : !> \param nfft ...
431 : !> \param n2 ...
432 : !> \param lot ...
433 : !> \param n1 ...
434 : !> \param lzt ...
435 : !> \param zw ...
436 : !> \param zt ...
437 : ! **************************************************************************************************
438 169146 : SUBROUTINE P_unswitch_downcorn(nfft, n2, lot, n1, lzt, zw, zt)
439 : INTEGER, INTENT(in) :: nfft, n2, lot, n1, lzt
440 : COMPLEX(KIND=dp), DIMENSION(lot, n2), INTENT(in) :: zw
441 : COMPLEX(KIND=dp), DIMENSION(lzt, n1), &
442 : INTENT(inout) :: zt
443 :
444 : INTEGER :: i, j
445 :
446 3810974 : DO j = 1, nfft
447 95909450 : DO i = 1, n2
448 95740304 : zt(i, j) = zw(j, i)
449 : END DO
450 : END DO
451 :
452 169146 : END SUBROUTINE P_unswitch_downcorn
453 :
454 : ! **************************************************************************************************
455 : !> \brief ...
456 : !> \param j3 ...
457 : !> \param nfft ...
458 : !> \param Jp2stf ...
459 : !> \param J2stf ...
460 : !> \param lot ...
461 : !> \param n1 ...
462 : !> \param md2 ...
463 : !> \param nd3 ...
464 : !> \param nproc ...
465 : !> \param zw ...
466 : !> \param zmpi1 ...
467 : ! **************************************************************************************************
468 169146 : SUBROUTINE P_unmpiswitch_downcorn(j3, nfft, Jp2stf, J2stf, lot, n1, md2, nd3, nproc, zw, zmpi1)
469 : INTEGER, INTENT(in) :: j3, nfft
470 : INTEGER, INTENT(inout) :: Jp2stf, J2stf
471 : INTEGER, INTENT(in) :: lot, n1, md2, nd3, nproc
472 : COMPLEX(KIND=dp), DIMENSION(lot, n1), INTENT(in) :: zw
473 : COMPLEX(KIND=dp), &
474 : DIMENSION(n1, md2/nproc, nd3/nproc, nproc), &
475 : INTENT(inout) :: zmpi1
476 :
477 : INTEGER :: I1, J2, Jp2, mfft
478 :
479 169146 : mfft = 0
480 354650 : DO Jp2 = Jp2stf, nproc
481 3840810 : DO J2 = J2stf, md2/nproc
482 3655306 : mfft = mfft + 1
483 3655306 : IF (mfft > nfft) THEN
484 13478 : Jp2stf = Jp2
485 13478 : J2stf = J2
486 13478 : RETURN
487 : END IF
488 95925808 : DO I1 = 1, n1
489 95740304 : zmpi1(I1, J2, j3, Jp2) = zw(mfft, I1)
490 : END DO
491 : END DO
492 341172 : J2stf = 1
493 : END DO
494 : END SUBROUTINE P_unmpiswitch_downcorn
495 :
496 : ! **************************************************************************************************
497 : !> \brief (Based on suitable modifications of S.Goedecker routines)
498 : !> Restore data into output array
499 : !> \param md1 Dimensions of the undistributed part of the real grid
500 : !> \param md3 Dimensions of the undistributed part of the real grid
501 : !> \param lot ...
502 : !> \param nfft number of planes
503 : !> \param n3 (twice the) dimension of the last FFTtransform
504 : !> \param zw FFT work array
505 : !> \param zf Original distributed density as well as
506 : !> Distributed solution of the poisson equation (inout)
507 : !> \param scal Needed to achieve unitarity and correct dimensions
508 : !> \date February 2006
509 : !> \author S. Goedecker, L. Genovese
510 : !> \note Assuming that high frequencies are in the corners
511 : !> and that n3 is multiple of 4
512 : !>
513 : !> RESTRICTIONS on USAGE
514 : !> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
515 : !> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
516 : !> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
517 : !> This file is distributed under the terms of the
518 : !> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
519 : ! **************************************************************************************************
520 311348 : SUBROUTINE P_unfill_downcorn(md1, md3, lot, nfft, n3, zw, zf, scal)
521 : INTEGER, INTENT(in) :: md1, md3, lot, nfft, n3
522 : COMPLEX(KIND=dp), DIMENSION(lot, n3), INTENT(in) :: zw
523 : REAL(KIND=dp), DIMENSION(md1, md3), INTENT(inout) :: zf
524 : REAL(KIND=dp), INTENT(in) :: scal
525 :
526 : INTEGER :: i1, i3
527 : REAL(KIND=dp) :: pot1
528 :
529 7534542 : DO i3 = 1, n3
530 179338112 : DO i1 = 1, nfft
531 171803570 : pot1 = scal*REAL(zw(i1, i3), dp)
532 179026764 : zf(i1, i3) = pot1
533 : END DO
534 : END DO
535 :
536 311348 : END SUBROUTINE P_unfill_downcorn
537 :
538 : ! **************************************************************************************************
539 : !> \brief ...
540 : !> \param md1 ...
541 : !> \param md3 ...
542 : !> \param lot ...
543 : !> \param nfft ...
544 : !> \param n3 ...
545 : !> \param zf ...
546 : !> \param zw ...
547 : ! **************************************************************************************************
548 311348 : SUBROUTINE P_fill_upcorn(md1, md3, lot, nfft, n3, zf, zw)
549 : INTEGER, INTENT(in) :: md1, md3, lot, nfft, n3
550 : REAL(KIND=dp), DIMENSION(md1, md3), INTENT(in) :: zf
551 : COMPLEX(KIND=dp), DIMENSION(lot, n3), &
552 : INTENT(inout) :: zw
553 :
554 : INTEGER :: i1, i3
555 :
556 7534542 : DO i3 = 1, n3
557 179338112 : DO i1 = 1, nfft
558 179026764 : zw(i1, i3) = CMPLX(zf(i1, i3), 0.0_dp, dp)
559 : END DO
560 : END DO
561 :
562 311348 : END SUBROUTINE P_fill_upcorn
563 :
564 : ! **************************************************************************************************
565 : !> \brief (Based on suitable modifications of S.Goedecker routines)
566 : !> Assign the correct planes to the work array zmpi2
567 : !> in order to prepare for interprocessor data transposition.
568 : !> \param i1 Starting points of the plane and number of remaining lines
569 : !> \param j2 Starting points of the plane and number of remaining lines
570 : !> \param lot Starting points of the plane and number of remaining lines
571 : !> \param nfft Starting points of the plane and number of remaining lines
572 : !> \param n1 logical dimension of the FFT transform, reference for work arrays
573 : !> \param n3 logical dimension of the FFT transform, reference for work arrays
574 : !> \param md2 Dimensions of real grid
575 : !> \param nproc ...
576 : !> \param nd3 Dimensions of the kernel
577 : !> \param zw Work array (input)
578 : !> \param zmpi2 Work array for multiprocessor manipulation (output)
579 : !> \date February 2006
580 : !> \author S. Goedecker, L. Genovese
581 : !> \note
582 : !> RESTRICTIONS on USAGE
583 : !> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
584 : !> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
585 : !> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
586 : !> This file is distributed under the terms of the
587 : !> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
588 : ! **************************************************************************************************
589 311348 : SUBROUTINE scramble_P(i1, j2, lot, nfft, n1, n3, md2, nproc, nd3, zw, zmpi2)
590 : INTEGER, INTENT(in) :: i1, j2, lot, nfft, n1, n3, md2, nproc, &
591 : nd3
592 : COMPLEX(KIND=dp), DIMENSION(lot, n3), INTENT(in) :: zw
593 : COMPLEX(KIND=dp), DIMENSION(n1, md2/nproc, nd3), &
594 : INTENT(inout) :: zmpi2
595 :
596 : INTEGER :: i, i3
597 :
598 4205680 : DO i3 = 1, n3/2 + 1
599 96304156 : DO i = 0, nfft - 1
600 95992808 : zmpi2(i1 + i, j2, i3) = zw(i + 1, i3)
601 : END DO
602 : END DO
603 :
604 311348 : END SUBROUTINE scramble_P
605 :
606 : ! **************************************************************************************************
607 : !> \brief (Based on suitable modifications of S.Goedecker routines)
608 : !> Insert the correct planes of the work array zmpi2
609 : !> in order to prepare for backward FFT transform
610 : !> \param i1 Starting points of the plane and number of remaining lines
611 : !> \param j2 Starting points of the plane and number of remaining lines
612 : !> \param lot Starting points of the plane and number of remaining lines
613 : !> \param nfft Starting points of the plane and number of remaining lines
614 : !> \param n1 logical dimension of the FFT transform, reference for work arrays
615 : !> \param n3 logical dimension of the FFT transform, reference for work arrays
616 : !> \param md2 Dimensions of real grid
617 : !> \param nproc ...
618 : !> \param nd3 Dimensions of the kernel
619 : !> \param zmpi2 Work array for multiprocessor manipulation (output)
620 : !> \param zw Work array (input)
621 : !> \date February 2006
622 : !> \author S. Goedecker, L. Genovese
623 : !> \note
624 : !> RESTRICTIONS on USAGE
625 : !> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
626 : !> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
627 : !> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
628 : !> This file is distributed under the terms of the
629 : !> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
630 : ! **************************************************************************************************
631 311348 : SUBROUTINE unscramble_P(i1, j2, lot, nfft, n1, n3, md2, nproc, nd3, zmpi2, zw)
632 : INTEGER, INTENT(in) :: i1, j2, lot, nfft, n1, n3, md2, nproc, &
633 : nd3
634 : COMPLEX(KIND=dp), DIMENSION(n1, md2/nproc, nd3), &
635 : INTENT(in) :: zmpi2
636 : COMPLEX(KIND=dp), DIMENSION(lot, n3), &
637 : INTENT(inout) :: zw
638 :
639 : INTEGER :: i, i3, j3
640 :
641 311348 : i3 = 1
642 7047570 : DO i = 0, nfft - 1
643 7047570 : zw(i + 1, i3) = zmpi2(i1 + i, j2, i3)
644 : END DO
645 :
646 3894332 : DO i3 = 2, n3/2 + 1
647 3582984 : j3 = n3 + 2 - i3
648 89256586 : DO i = 0, nfft - 1
649 85362254 : zw(i + 1, j3) = CONJG(zmpi2(i1 + i, j2, i3))
650 88945238 : zw(i + 1, i3) = zmpi2(i1 + i, j2, i3)
651 : END DO
652 : END DO
653 :
654 311348 : END SUBROUTINE unscramble_P
655 :
656 : ! **************************************************************************************************
657 : !> \brief (Based on suitable modifications of S.Goedecker routines)
658 : !> Multiply with the kernel taking into account its symmetry
659 : !> Conceived to be used into convolution loops
660 : !> \param n1 ...
661 : !> \param n2 ...
662 : !> \param n3 ...
663 : !> \param lot ...
664 : !> \param nfft ...
665 : !> \param jS ...
666 : !> \param i3 ...
667 : !> \param zw Work array (input/output)
668 : !> n1,n2: logical dimension of the FFT transform, reference for zw
669 : !> nd1,nd2: Dimensions of POT
670 : !> jS,j3,nfft: starting point of the plane and number of remaining lines
671 : !>
672 : !> RESTRICTIONS on USAGE
673 : !> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
674 : !> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
675 : !> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
676 : !> This file is distributed under the terms of the
677 : !> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
678 : !> \param hx ...
679 : !> \param hy ...
680 : !> \param hz ...
681 : !> \date February 2006
682 : !> \author S. Goedecker, L. Genovese
683 : ! **************************************************************************************************
684 169146 : SUBROUTINE P_multkernel(n1, n2, n3, lot, nfft, jS, i3, zw, hx, hy, hz)
685 : INTEGER, INTENT(in) :: n1, n2, n3, lot, nfft, jS, i3
686 : COMPLEX(KIND=dp), DIMENSION(lot, n2), &
687 : INTENT(inout) :: zw
688 : REAL(KIND=dp), INTENT(in) :: hx, hy, hz
689 :
690 : INTEGER :: i1, i2, j1, j2, j3
691 : REAL(KIND=dp) :: fourpi2, ker, mu3, p1, p2
692 :
693 169146 : fourpi2 = 4._dp*pi**2
694 169146 : j3 = i3 !n3/2+1-abs(n3/2+2-i3)
695 169146 : mu3 = REAL(j3 - 1, KIND=dp)/REAL(n3, KIND=dp)
696 169146 : mu3 = (mu3/hy)**2 !beware of the exchanged dimension
697 : !Body
698 : !generic case
699 4063478 : DO i2 = 1, n2
700 96161954 : DO i1 = 1, nfft
701 92098476 : j1 = i1 + jS - 1
702 92098476 : j1 = j1 - (j1/(n1/2 + 2))*n1 !n1/2+1-abs(n1/2+2-jS-i1)
703 92098476 : j2 = i2 - (i2/(n2/2 + 2))*n2 !n2/2+1-abs(n2/2+1-i2)
704 92098476 : p1 = REAL(j1 - 1, KIND=dp)/REAL(n1, KIND=dp)
705 92098476 : p2 = REAL(j2 - 1, KIND=dp)/REAL(n2, KIND=dp)
706 92098476 : ker = -fourpi2*((p1/hx)**2 + (p2/hz)**2 + mu3) !beware of the exchanged dimension
707 92098476 : IF (ker /= 0._dp) ker = 1._dp/ker
708 95992808 : zw(i1, i2) = zw(i1, i2)*ker
709 : END DO
710 : END DO
711 :
712 169146 : END SUBROUTINE P_multkernel
713 :
714 : ! **************************************************************************************************
715 : !> \brief (Based on suitable modifications of S.Goedecker routines)
716 : !> Multiply with the kernel taking into account its symmetry
717 : !> Conceived to be used into convolution loops
718 : !> \param nd1 ...
719 : !> \param nd2 ...
720 : !> \param n1 ...
721 : !> \param n2 ...
722 : !> \param lot ...
723 : !> \param nfft ...
724 : !> \param jS ...
725 : !> \param pot Kernel, symmetric and real, half the length
726 : !> \param zw Work array (input/output)
727 : !> n1,n2: logical dimension of the FFT transform, reference for zw
728 : !> nd1,nd2: Dimensions of POT
729 : !> jS, nfft: starting point of the plane and number of remaining lines
730 : !>
731 : !> RESTRICTIONS on USAGE
732 : !> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
733 : !> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
734 : !> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
735 : !> This file is distributed under the terms of the
736 : !> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
737 : !> \date February 2006
738 : !> \author S. Goedecker, L. Genovese
739 : ! **************************************************************************************************
740 2246238 : SUBROUTINE multkernel(nd1, nd2, n1, n2, lot, nfft, jS, pot, zw)
741 : INTEGER, INTENT(in) :: nd1, nd2, n1, n2, lot, nfft, jS
742 : REAL(KIND=dp), DIMENSION(nd1, nd2), INTENT(in) :: pot
743 : COMPLEX(KIND=dp), DIMENSION(lot, n2), &
744 : INTENT(inout) :: zw
745 :
746 : INTEGER :: i2, j, j1, j2
747 :
748 42306128 : DO j = 1, nfft
749 40059890 : j1 = j + jS - 1
750 40059890 : j1 = j1 + (j1/(n1/2 + 2))*(n1 + 2 - 2*j1)
751 42306128 : zw(j, 1) = zw(j, 1)*pot(j1, 1)
752 : END DO
753 :
754 : !generic case
755 123417162 : DO i2 = 2, n2/2
756 2060343990 : DO j = 1, nfft
757 1936926828 : j1 = j + jS - 1
758 1936926828 : j1 = j1 + (j1/(n1/2 + 2))*(n1 + 2 - 2*j1)
759 1936926828 : j2 = n2 + 2 - i2
760 1936926828 : zw(j, i2) = zw(j, i2)*pot(j1, i2)
761 2058097752 : zw(j, j2) = zw(j, j2)*pot(j1, i2)
762 : END DO
763 : END DO
764 :
765 : !case i2=n2/2+1
766 42306128 : DO j = 1, nfft
767 40059890 : j1 = j + jS - 1
768 40059890 : j1 = j1 + (j1/(n1/2 + 2))*(n1 + 2 - 2*j1)
769 40059890 : j2 = n2/2 + 1
770 42306128 : zw(j, j2) = zw(j, j2)*pot(j1, j2)
771 : END DO
772 :
773 2246238 : END SUBROUTINE multkernel
774 :
775 : ! **************************************************************************************************
776 : !> \brief !HERE POT MUST BE THE KERNEL (BEWARE THE HALF DIMENSION)
777 : !> ****h* BigDFT/S_PoissonSolver
778 : !> (Based on suitable modifications of S.Goedecker routines)
779 : !> Applies the local FFT space Kernel to the density in Real space.
780 : !> Does NOT calculate the LDA exchange-correlation terms
781 : !> \param n1 logical dimension of the transform.
782 : !> \param n2 logical dimension of the transform.
783 : !> \param n3 logical dimension of the transform.
784 : !> \param nd1 Dimension of POT
785 : !> \param nd2 Dimension of POT
786 : !> \param nd3 Dimension of POT
787 : !> \param md1 Dimension of ZF
788 : !> \param md2 Dimension of ZF
789 : !> \param md3 Dimension of ZF
790 : !> \param nproc number of processors used as returned by MPI_COMM_SIZE
791 : !> \param iproc [0:nproc-1] number of processor as returned by MPI_COMM_RANK
792 : !> \param pot Kernel, only the distributed part (REAL)
793 : !> POT(i1,i2,i3)
794 : !> i1=1,nd1 , i2=1,nd2 , i3=1,nd3/nproc
795 : !> \param zf Density (input/output)
796 : !> ZF(i1,i3,i2)
797 : !> i1=1,md1 , i2=1,md2/nproc , i3=1,md3
798 : !> \param scal factor of renormalization of the FFT in order to acheve unitarity
799 : !> and the correct dimension
800 : !> \param mpi_group ...
801 : !> \date October 2006
802 : !> \author S. Goedecker, L. Genovese
803 : !> \note
804 : !> RESTRICTIONS on USAGE
805 : !> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
806 : !> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
807 : !> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
808 : !> This file is distributed under the terms of the
809 : !> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
810 : ! **************************************************************************************************
811 54 : SUBROUTINE S_PoissonSolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, pot, zf, &
812 : scal, mpi_group)
813 : INTEGER, INTENT(in) :: n1, n2, n3, nd1, nd2, nd3, md1, md2, &
814 : md3, nproc, iproc
815 : REAL(KIND=dp), DIMENSION(nd1, nd2, nd3/nproc), &
816 : INTENT(in) :: pot
817 : REAL(KIND=dp), DIMENSION(md1, md3, md2/nproc), &
818 : INTENT(inout) :: zf
819 : REAL(KIND=dp), INTENT(in) :: scal
820 :
821 : CLASS(mp_comm_type), INTENT(in) :: mpi_group
822 :
823 : CHARACTER(len=*), PARAMETER :: routineN = 'S_PoissonSolver'
824 : INTEGER, PARAMETER :: ncache_optimal = 8*1024
825 :
826 : INTEGER :: handle, i1, i3, j, j2, J2stb, J2stf, j3, Jp2stb, &
827 : Jp2stf, lot1, lot2, lot3, lzt, ma, mb, ncache, nfft, stat, &
828 : final_chunk_size1, final_chunk_size2, final_chunk_size3
829 : REAL(kind=dp) :: twopion
830 54 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cosinarr
831 54 : COMPLEX(KIND=dp), POINTER, CONTIGUOUS, DIMENSION(:, :) :: zt
832 54 : COMPLEX(KIND=dp), POINTER, CONTIGUOUS, DIMENSION(:) :: zw1, zw2
833 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: zmpi2
834 : COMPLEX(KIND=dp), ALLOCATABLE, &
835 54 : DIMENSION(:, :, :, :) :: zmpi1
836 : TYPE(fft_plan_type) :: fft_plan_bw3, fft_plan_bw3_last, fft_plan_fw3, fft_plan_fw3_last, &
837 : fft_plan_bw1, fft_plan_bw1_last, fft_plan_fw1, fft_plan_fw1_last, &
838 : fft_plan_bw2, fft_plan_bw2_last, fft_plan_fw2, fft_plan_fw2_last
839 :
840 54 : CALL timeset(routineN, handle)
841 : ! check input
842 54 : IF (MOD(n3, 2) /= 0) CPABORT("Parallel convolution:ERROR:n3")
843 54 : IF (nd1 < n1/2 + 1) CPABORT("Parallel convolution:ERROR:nd1")
844 54 : IF (nd2 < n2/2 + 1) CPABORT("Parallel convolution:ERROR:nd2")
845 54 : IF (nd3 < n3/2 + 1) CPABORT("Parallel convolution:ERROR:nd3")
846 54 : IF (md1 < n1) CPABORT("Parallel convolution:ERROR:md1")
847 54 : IF (md2 < n2) CPABORT("Parallel convolution:ERROR:md2")
848 54 : IF (md3 < n3/2) CPABORT("Parallel convolution:ERROR:md3")
849 54 : IF (MOD(nd3, nproc) /= 0) CPABORT("Parallel convolution:ERROR:nd3")
850 54 : IF (MOD(md2, nproc) /= 0) CPABORT("Parallel convolution:ERROR:md2")
851 :
852 : !defining work arrays dimensions
853 54 : ncache = ncache_optimal
854 54 : IF (ncache <= MAX(n1, n2, n3/2)*4) ncache = MAX(n1, n2, n3/2)*4
855 :
856 54 : lzt = n2
857 54 : IF (MOD(n2, 2) == 0) lzt = lzt + 1
858 54 : IF (MOD(n2, 4) == 0) lzt = lzt + 1 !maybe this is useless
859 :
860 108 : CALL fft_alloc(zw1, [ncache/4])
861 54 : zw1 = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
862 108 : CALL fft_alloc(zw2, [ncache/4])
863 54 : zw2 = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
864 162 : CALL fft_alloc(zt, [lzt, n1])
865 54 : zt = CMPLX(0.0_dp, 0.0_dp, dp)
866 270 : ALLOCATE (zmpi2(n1, md2/nproc, nd3), SOURCE=CMPLX(0.0_dp, 0.0_dp, dp))
867 216 : ALLOCATE (cosinarr(2, n3/2))
868 324 : IF (nproc > 1) ALLOCATE (zmpi1(n1, md2/nproc, nd3/nproc, nproc), SOURCE=CMPLX(0.0_dp, 0.0_dp, dp))
869 :
870 : !Calculating array of phases for HalFFT decoding
871 54 : twopion = 8._dp*ATAN(1._dp)/REAL(n3, KIND=dp)
872 2970 : DO i3 = 1, n3/2
873 2916 : cosinarr(1, i3) = COS(twopion*(i3 - 1))
874 2970 : cosinarr(2, i3) = -SIN(twopion*(i3 - 1))
875 : END DO
876 :
877 : !initializing integral
878 :
879 : ! transform along z axis
880 54 : lot1 = ncache/(4*n1)
881 54 : lot2 = ncache/(4*n2)
882 54 : lot3 = ncache/(2*n3)
883 :
884 : ! The size of the last chunk
885 54 : final_chunk_size1 = MOD(n2, lot1)
886 54 : final_chunk_size2 = MOD(n1, lot2)
887 54 : final_chunk_size3 = MOD(n1, lot3)
888 :
889 : ! Prevent OOB-access if n2 < lot1
890 54 : IF (n2 >= lot1) THEN
891 54 : CALL fft_create_plan_1d(fft_plan_bw1, BWFFT, .TRUE., .TRUE., lot1, lzt, n1, lot1, zw1, zt)
892 54 : CALL fft_create_plan_1d(fft_plan_fw1, FWFFT, .TRUE., .TRUE., lzt, lot1, n1, lot1, zt, zw1)
893 : END IF
894 54 : IF (final_chunk_size1 > 0) THEN
895 : CALL fft_create_plan_1d(fft_plan_bw1_last, BWFFT, .TRUE., .TRUE., lot1, lzt, n1, &
896 54 : final_chunk_size1, zw1, zt)
897 : CALL fft_create_plan_1d(fft_plan_fw1_last, FWFFT, .TRUE., .TRUE., lzt, lot1, n1, &
898 54 : final_chunk_size1, zt, zw1)
899 : END IF
900 :
901 : ! Prevent OOB-access if n2 < lot1
902 54 : IF (n1 >= lot2) THEN
903 54 : CALL fft_create_plan_1d(fft_plan_bw2, BWFFT, .TRUE., .TRUE., lot2, lot2, n2, lot2, zw1, zw2)
904 54 : CALL fft_create_plan_1d(fft_plan_fw2, FWFFT, .TRUE., .TRUE., lot2, lot2, n2, lot2, zw2, zw1)
905 : END IF
906 54 : IF (final_chunk_size2 > 0) THEN
907 : CALL fft_create_plan_1d(fft_plan_bw2_last, BWFFT, .TRUE., .TRUE., lot2, lot2, n2, &
908 54 : final_chunk_size2, zw1, zw2)
909 : CALL fft_create_plan_1d(fft_plan_fw2_last, FWFFT, .TRUE., .TRUE., lot2, lot2, n2, &
910 54 : final_chunk_size2, zw2, zw1)
911 : END IF
912 :
913 : ! Prevent OOB-access if n1 < lot3
914 54 : IF (n1 >= lot3) THEN
915 54 : CALL fft_create_plan_1d(fft_plan_fw3, FWFFT, .TRUE., .TRUE., lot3, lot3, n3/2, lot3, zw1, zw2)
916 54 : CALL fft_create_plan_1d(fft_plan_bw3, BWFFT, .TRUE., .TRUE., lot3, lot3, n3/2, lot3, zw1, zw2)
917 : END IF
918 54 : IF (final_chunk_size3 > 0) THEN
919 : CALL fft_create_plan_1d(fft_plan_fw3_last, FWFFT, .TRUE., .TRUE., &
920 54 : lot3, lot3, n3/2, final_chunk_size3, zw1, zw2)
921 : CALL fft_create_plan_1d(fft_plan_bw3_last, BWFFT, .TRUE., .TRUE., &
922 54 : lot3, lot3, n3/2, final_chunk_size3, zw1, zw2)
923 : END IF
924 :
925 1512 : DO j2 = 1, md2/nproc
926 : !this condition ensures that we manage only the interesting part for the FFT
927 1512 : IF (iproc*(md2/nproc) + j2 <= n2) THEN
928 4374 : DO i1 = 1, n1, lot3
929 2916 : ma = i1
930 2916 : mb = MIN(i1 + (lot3 - 1), n1)
931 2916 : nfft = mb - ma + 1
932 :
933 : !inserting real data into complex array of half length
934 2916 : CALL halfill_upcorn(md1, md3, lot3, nfft, n3, zf(i1, 1, j2), zw1)
935 :
936 : !performing FFT
937 : !input: I1,I3,J2,(Jp2)
938 : ! This is equivalent to
939 2916 : IF (nfft == lot3) THEN
940 1458 : CALL fft_1d(fft_plan_bw3, zw1, zw2, 1.0_dp, stat)
941 : ELSE
942 1458 : CALL fft_1d(fft_plan_bw3_last, zw1, zw2, 1.0_dp, stat)
943 : END IF
944 : !output: I1,i3,J2,(Jp2)
945 : !unpacking FFT in order to restore correct result,
946 : !while exchanging components
947 : !input: I1,i3,J2,(Jp2)
948 4374 : CALL scramble_unpack(i1, j2, lot3, nfft, n1, n3, md2, nproc, nd3, zw2, zmpi2, cosinarr)
949 : !output: I1,J2,i3,(Jp2)
950 : END DO
951 : END IF
952 : END DO
953 : !Interprocessor data transposition
954 : !input: I1,J2,j3,jp3,(Jp2)
955 54 : IF (nproc > 1) THEN
956 54 : CALL mpi_group%alltoall(zmpi2, zmpi1, n1*(md2/nproc)*(nd3/nproc))
957 : END IF
958 : !output: I1,J2,j3,Jp2,(jp3)
959 :
960 : !now each process perform complete convolution of its planes
961 1566 : DO j3 = 1, nd3/nproc
962 : !this condition ensures that we manage only the interesting part for the FFT
963 1566 : IF (iproc*(nd3/nproc) + j3 <= n3/2 + 1) THEN
964 1485 : Jp2stb = 1
965 1485 : J2stb = 1
966 1485 : Jp2stf = 1
967 1485 : J2stf = 1
968 :
969 : ! transform along x axis
970 :
971 4455 : DO j = 1, n2, lot1
972 2970 : ma = j
973 2970 : mb = MIN(j + (lot1 - 1), n2)
974 2970 : nfft = mb - ma + 1
975 :
976 : !reverse index ordering, leaving the planes to be transformed at the end
977 : !input: I1,J2,j3,Jp2,(jp3)
978 2970 : IF (nproc == 1) THEN
979 0 : CALL S_mpiswitch_upcorn(j3, nfft, Jp2stb, J2stb, lot1, n1, md2, nd3, nproc, zmpi2, zw1)
980 : ELSE
981 2970 : CALL S_mpiswitch_upcorn(j3, nfft, Jp2stb, J2stb, lot1, n1, md2, nd3, nproc, zmpi1, zw1)
982 : END IF
983 : !output: J2,Jp2,I1,j3,(jp3)
984 :
985 : !performing FFT
986 : !input: I2,I1,j3,(jp3)
987 : ! This is equivalent to
988 4455 : IF (nfft == lot1) THEN
989 1485 : CALL fft_1d(fft_plan_bw1, zw1, zt(j:, 1), 1.0_dp, stat)
990 : ELSE
991 1485 : CALL fft_1d(fft_plan_bw1_last, zw1, zt(j:, 1), 1.0_dp, stat)
992 : END IF
993 : !output: I2,i1,j3,(jp3)
994 : END DO
995 :
996 : !transform along y axis
997 :
998 4455 : DO j = 1, n1, lot2
999 2970 : ma = j
1000 2970 : mb = MIN(j + (lot2 - 1), n1)
1001 2970 : nfft = mb - ma + 1
1002 :
1003 : !reverse ordering
1004 : !input: I2,i1,j3,(jp3)
1005 2970 : CALL S_switch_upcorn(nfft, n2, lot2, n1, lzt, zt(:, j), zw1)
1006 : !output: i1,I2,j3,(jp3)
1007 :
1008 : !performing FFT
1009 : !input: i1,I2,j3,(jp3)
1010 : ! This is equivalent to
1011 2970 : IF (nfft == lot2) THEN
1012 1485 : CALL fft_1d(fft_plan_bw2, zw1, zw2, 1.0_dp, stat)
1013 : ELSE
1014 1485 : CALL fft_1d(fft_plan_bw2_last, zw1, zw2, 1.0_dp, stat)
1015 : END IF
1016 : !output: i1,i2,j3,(jp3)
1017 :
1018 : !Multiply with kernel in fourier space
1019 2970 : CALL multkernel(nd1, nd2, n1, n2, lot2, nfft, j, pot(1, 1, j3), zw2)
1020 :
1021 : !TRANSFORM BACK IN REAL SPACE
1022 :
1023 : !transform along y axis
1024 : !input: i1,i2,j3,(jp3)
1025 : ! This is equivalent to
1026 2970 : IF (nfft == lot2) THEN
1027 1485 : CALL fft_1d(fft_plan_fw2, zw2, zw1, 1.0_dp, stat)
1028 : ELSE
1029 1485 : CALL fft_1d(fft_plan_fw2_last, zw2, zw1, 1.0_dp, stat)
1030 : END IF
1031 :
1032 : !reverse ordering
1033 : !input: i1,I2,j3,(jp3)
1034 4455 : CALL S_unswitch_downcorn(nfft, n2, lot2, n1, lzt, zw1, zt(:, j))
1035 : !output: I2,i1,j3,(jp3)
1036 : END DO
1037 :
1038 : !transform along x axis
1039 : !input: I2,i1,j3,(jp3)
1040 4455 : DO j = 1, n2, lot1
1041 2970 : ma = j
1042 2970 : mb = MIN(j + (lot1 - 1), n2)
1043 2970 : nfft = mb - ma + 1
1044 :
1045 : !performing FFT
1046 : ! This is equivalent to
1047 2970 : IF (nfft == lot1) THEN
1048 1485 : CALL fft_1d(fft_plan_fw1, zt(j:, 1), zw2, 1.0_dp, stat)
1049 : ELSE
1050 1485 : CALL fft_1d(fft_plan_fw1_last, zt(j:, 1), zw2, 1.0_dp, stat)
1051 : END IF
1052 : !output: I2,I1,j3,(jp3)
1053 :
1054 : !reverse ordering
1055 : !input: J2,Jp2,I1,j3,(jp3)
1056 4455 : IF (nproc == 1) THEN
1057 0 : CALL S_unmpiswitch_downcorn(j3, nfft, Jp2stf, J2stf, lot1, n1, md2, nd3, nproc, zw2, zmpi2)
1058 : ELSE
1059 2970 : CALL S_unmpiswitch_downcorn(j3, nfft, Jp2stf, J2stf, lot1, n1, md2, nd3, nproc, zw2, zmpi1)
1060 : END IF
1061 : ! output: I1,J2,j3,Jp2,(jp3)
1062 : END DO
1063 : END IF
1064 : END DO
1065 :
1066 : !Interprocessor data transposition
1067 : !input: I1,J2,j3,Jp2,(jp3)
1068 54 : IF (nproc > 1) THEN
1069 : !communication scheduling
1070 54 : CALL mpi_group%alltoall(zmpi1, zmpi2, n1*(md2/nproc)*(nd3/nproc))
1071 : END IF
1072 :
1073 : !output: I1,J2,j3,jp3,(Jp2)
1074 :
1075 : !transform along z axis
1076 : !input: I1,J2,i3,(Jp2)
1077 1512 : DO j2 = 1, md2/nproc
1078 : !this condition ensures that we manage only the interesting part for the FFT
1079 1512 : IF (iproc*(md2/nproc) + j2 <= n2) THEN
1080 4374 : DO i1 = 1, n1, lot3
1081 2916 : ma = i1
1082 2916 : mb = MIN(i1 + (lot3 - 1), n1)
1083 2916 : nfft = mb - ma + 1
1084 :
1085 : !reverse ordering and repack the FFT data in order to be backward HalFFT transformed
1086 : !input: I1,J2,i3,(Jp2)
1087 2916 : CALL unscramble_pack(i1, j2, lot3, nfft, n1, n3, md2, nproc, nd3, zmpi2, zw1, cosinarr)
1088 : !output: I1,i3,J2,(Jp2)
1089 :
1090 : !performing FFT
1091 : !input: I1,i3,J2,(Jp2)
1092 : ! This is equivalent to
1093 2916 : IF (nfft == lot3) THEN
1094 1458 : CALL fft_1d(fft_plan_fw3, zw1, zw2, 1.0_dp, stat)
1095 : ELSE
1096 1458 : CALL fft_1d(fft_plan_fw3_last, zw1, zw2, 1.0_dp, stat)
1097 : END IF
1098 : !output: I1,I3,J2,(Jp2)
1099 :
1100 : !rebuild the output array
1101 4374 : CALL unfill_downcorn(md1, md3, lot3, nfft, n3, zw2, zf(i1, 1, j2), scal)
1102 :
1103 : !integrate local pieces together
1104 : !ehartree=ehartree+0.5_dp*ehartreetmp*hx*hy*hz
1105 : END DO
1106 : END IF
1107 : END DO
1108 :
1109 54 : IF (n2 >= lot1) CALL fft_destroy_plan(fft_plan_bw1)
1110 54 : IF (final_chunk_size1 > 0) CALL fft_destroy_plan(fft_plan_bw1_last)
1111 54 : IF (n2 >= lot1) CALL fft_destroy_plan(fft_plan_fw1)
1112 54 : IF (final_chunk_size1 > 0) CALL fft_destroy_plan(fft_plan_fw1_last)
1113 :
1114 54 : IF (n1 >= lot2) CALL fft_destroy_plan(fft_plan_bw2)
1115 54 : IF (final_chunk_size2 > 0) CALL fft_destroy_plan(fft_plan_bw2_last)
1116 54 : IF (n1 >= lot2) CALL fft_destroy_plan(fft_plan_fw2)
1117 54 : IF (final_chunk_size2 > 0) CALL fft_destroy_plan(fft_plan_fw2_last)
1118 :
1119 54 : IF (n1 >= lot3) CALL fft_destroy_plan(fft_plan_fw3)
1120 54 : IF (final_chunk_size3 > 0) CALL fft_destroy_plan(fft_plan_fw3_last)
1121 54 : IF (n1 >= lot3) CALL fft_destroy_plan(fft_plan_bw3)
1122 54 : IF (final_chunk_size3 > 0) CALL fft_destroy_plan(fft_plan_bw3_last)
1123 :
1124 : !De-allocations
1125 54 : DEALLOCATE (zmpi2)
1126 54 : CALL fft_dealloc(zw1)
1127 54 : CALL fft_dealloc(zw2)
1128 54 : CALL fft_dealloc(zt)
1129 54 : DEALLOCATE (cosinarr)
1130 54 : IF (nproc > 1) DEALLOCATE (zmpi1)
1131 :
1132 54 : CALL timestop(handle)
1133 2052 : END SUBROUTINE S_PoissonSolver
1134 :
1135 : ! **************************************************************************************************
1136 : !> \brief ...
1137 : !> \param j3 ...
1138 : !> \param nfft ...
1139 : !> \param Jp2stb ...
1140 : !> \param J2stb ...
1141 : !> \param lot ...
1142 : !> \param n1 ...
1143 : !> \param md2 ...
1144 : !> \param nd3 ...
1145 : !> \param nproc ...
1146 : !> \param zmpi1 ...
1147 : !> \param zw ...
1148 : ! **************************************************************************************************
1149 2970 : SUBROUTINE S_mpiswitch_upcorn(j3, nfft, Jp2stb, J2stb, lot, n1, md2, nd3, nproc, zmpi1, zw)
1150 : INTEGER, INTENT(in) :: j3, nfft
1151 : INTEGER, INTENT(inout) :: Jp2stb, J2stb
1152 : INTEGER, INTENT(in) :: lot, n1, md2, nd3, nproc
1153 : COMPLEX(KIND=dp), &
1154 : DIMENSION(n1, md2/nproc, nd3/nproc, nproc), &
1155 : INTENT(in) :: zmpi1
1156 : COMPLEX(KIND=dp), DIMENSION(lot, n1), &
1157 : INTENT(inout) :: zw
1158 :
1159 : INTEGER :: I1, J2, Jp2, mfft
1160 :
1161 2970 : mfft = 0
1162 5940 : DO Jp2 = Jp2stb, nproc
1163 84645 : DO J2 = J2stb, md2/nproc
1164 81675 : mfft = mfft + 1
1165 81675 : IF (mfft > nfft) THEN
1166 1485 : Jp2stb = Jp2
1167 1485 : J2stb = J2
1168 1485 : RETURN
1169 : END IF
1170 4413420 : DO I1 = 1, n1
1171 4410450 : zw(mfft, I1) = zmpi1(I1, J2, j3, Jp2)
1172 : END DO
1173 : END DO
1174 4455 : J2stb = 1
1175 : END DO
1176 : END SUBROUTINE S_mpiswitch_upcorn
1177 :
1178 : ! **************************************************************************************************
1179 : !> \brief ...
1180 : !> \param nfft ...
1181 : !> \param n2 ...
1182 : !> \param lot ...
1183 : !> \param n1 ...
1184 : !> \param lzt ...
1185 : !> \param zt ...
1186 : !> \param zw ...
1187 : ! **************************************************************************************************
1188 2970 : SUBROUTINE S_switch_upcorn(nfft, n2, lot, n1, lzt, zt, zw)
1189 : INTEGER, INTENT(in) :: nfft, n2, lot, n1, lzt
1190 : COMPLEX(KIND=dp), DIMENSION(lzt, n1), INTENT(in) :: zt
1191 : COMPLEX(KIND=dp), DIMENSION(lot, n2), &
1192 : INTENT(inout) :: zw
1193 :
1194 : INTEGER :: i, j
1195 :
1196 83160 : DO j = 1, nfft
1197 4413420 : DO i = 1, n2
1198 4410450 : zw(j, i) = zt(i, j)
1199 : END DO
1200 : END DO
1201 2970 : END SUBROUTINE S_switch_upcorn
1202 :
1203 : ! **************************************************************************************************
1204 : !> \brief ...
1205 : !> \param nfft ...
1206 : !> \param n2 ...
1207 : !> \param lot ...
1208 : !> \param n1 ...
1209 : !> \param lzt ...
1210 : !> \param zw ...
1211 : !> \param zt ...
1212 : ! **************************************************************************************************
1213 2970 : SUBROUTINE S_unswitch_downcorn(nfft, n2, lot, n1, lzt, zw, zt)
1214 : INTEGER, INTENT(in) :: nfft, n2, lot, n1, lzt
1215 : COMPLEX(KIND=dp), DIMENSION(lot, n2), INTENT(in) :: zw
1216 : COMPLEX(KIND=dp), DIMENSION(lzt, n1), &
1217 : INTENT(inout) :: zt
1218 :
1219 : INTEGER :: i, j
1220 :
1221 83160 : DO j = 1, nfft
1222 4413420 : DO i = 1, n2
1223 4410450 : zt(i, j) = zw(j, i)
1224 : END DO
1225 : END DO
1226 2970 : END SUBROUTINE S_unswitch_downcorn
1227 :
1228 : ! **************************************************************************************************
1229 : !> \brief ...
1230 : !> \param j3 ...
1231 : !> \param nfft ...
1232 : !> \param Jp2stf ...
1233 : !> \param J2stf ...
1234 : !> \param lot ...
1235 : !> \param n1 ...
1236 : !> \param md2 ...
1237 : !> \param nd3 ...
1238 : !> \param nproc ...
1239 : !> \param zw ...
1240 : !> \param zmpi1 ...
1241 : ! **************************************************************************************************
1242 2970 : SUBROUTINE S_unmpiswitch_downcorn(j3, nfft, Jp2stf, J2stf, lot, n1, md2, nd3, nproc, zw, zmpi1)
1243 : INTEGER, INTENT(in) :: j3, nfft
1244 : INTEGER, INTENT(inout) :: Jp2stf, J2stf
1245 : INTEGER, INTENT(in) :: lot, n1, md2, nd3, nproc
1246 : COMPLEX(KIND=dp), DIMENSION(lot, n1), INTENT(in) :: zw
1247 : COMPLEX(KIND=dp), &
1248 : DIMENSION(n1, md2/nproc, nd3/nproc, nproc), &
1249 : INTENT(inout) :: zmpi1
1250 :
1251 : INTEGER :: I1, J2, Jp2, mfft
1252 :
1253 2970 : mfft = 0
1254 5940 : DO Jp2 = Jp2stf, nproc
1255 84645 : DO J2 = J2stf, md2/nproc
1256 81675 : mfft = mfft + 1
1257 81675 : IF (mfft > nfft) THEN
1258 1485 : Jp2stf = Jp2
1259 1485 : J2stf = J2
1260 1485 : RETURN
1261 : END IF
1262 4413420 : DO I1 = 1, n1
1263 4410450 : zmpi1(I1, J2, j3, Jp2) = zw(mfft, I1)
1264 : END DO
1265 : END DO
1266 4455 : J2stf = 1
1267 : END DO
1268 : END SUBROUTINE S_unmpiswitch_downcorn
1269 :
1270 : ! **************************************************************************************************
1271 : !> \brief (Based on suitable modifications of S.Goedecker routines)
1272 : !> Restore data into output array, calculating in the meanwhile
1273 : !> Hartree energy of the potential
1274 : !> \param md1 Dimensions of the undistributed part of the real grid
1275 : !> \param md3 Dimensions of the undistributed part of the real grid
1276 : !> \param lot ...
1277 : !> \param nfft number of planes
1278 : !> \param n3 (twice the) dimension of the last FFTtransform.
1279 : !> \param zw FFT work array
1280 : !> \param zf Original distributed density as well as
1281 : !> Distributed solution of the poisson equation (inout)
1282 : !> \param scal Needed to achieve unitarity and correct dimensions
1283 : !> \date February 2006
1284 : !> \author S. Goedecker, L. Genovese
1285 : !> \note Assuming that high frequencies are in the corners
1286 : !> and that n3 is multiple of 4
1287 : !>
1288 : !> RESTRICTIONS on USAGE
1289 : !> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
1290 : !> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
1291 : !> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
1292 : !> This file is distributed under the terms of the
1293 : !> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
1294 : ! **************************************************************************************************
1295 676327 : SUBROUTINE unfill_downcorn(md1, md3, lot, nfft, n3, zw, zf, scal)
1296 : INTEGER, INTENT(in) :: md1, md3, lot, nfft, n3
1297 : COMPLEX(KIND=dp), DIMENSION(lot, n3/2), INTENT(in) :: zw
1298 : REAL(KIND=dp), DIMENSION(md1, md3), INTENT(inout) :: zf
1299 : REAL(KIND=dp), INTENT(in) :: scal
1300 :
1301 : INTEGER :: i1, i3
1302 : REAL(KIND=dp) :: pot1
1303 :
1304 18469469 : DO i3 = 1, n3/4
1305 504746046 : DO i1 = 1, nfft
1306 486276577 : pot1 = scal*REAL(zw(i1, i3), dp)
1307 : !ehartreetmp =ehartreetmp + pot1* zf(i1,2*i3-1)
1308 486276577 : zf(i1, 2*i3 - 1) = pot1
1309 486276577 : pot1 = scal*AIMAG(zw(i1, i3))
1310 : !ehartreetmp =ehartreetmp + pot1* zf(i1,2*i3)
1311 504069719 : zf(i1, 2*i3) = pot1
1312 : END DO
1313 : END DO
1314 676327 : END SUBROUTINE unfill_downcorn
1315 :
1316 : ! **************************************************************************************************
1317 : !> \brief ...
1318 : !> \param md1 ...
1319 : !> \param md3 ...
1320 : !> \param lot ...
1321 : !> \param nfft ...
1322 : !> \param n3 ...
1323 : !> \param zf ...
1324 : !> \param zw ...
1325 : ! **************************************************************************************************
1326 676327 : SUBROUTINE halfill_upcorn(md1, md3, lot, nfft, n3, zf, zw)
1327 : INTEGER :: md1, md3, lot, nfft, n3
1328 : REAL(KIND=dp) :: zf(md1, md3)
1329 : COMPLEX(KIND=dp) :: zw(lot, n3/2)
1330 :
1331 : INTEGER :: i1, i3
1332 :
1333 18469469 : DO i3 = 1, n3/4
1334 : ! WARNING: Assuming that high frequencies are in the corners
1335 : ! and that n3 is multiple of 4
1336 : !in principle we can relax this condition
1337 504746046 : DO i1 = 1, nfft
1338 504069719 : zw(i1, i3) = CMPLX(0.0_dp, 0.0_dp, dp)
1339 : END DO
1340 : END DO
1341 18469469 : DO i3 = n3/4 + 1, n3/2
1342 504746046 : DO i1 = 1, nfft
1343 504069719 : zw(i1, i3) = CMPLX(zf(i1, 2*i3 - 1 - n3/2), zf(i1, 2*i3 - n3/2), dp)
1344 : END DO
1345 : END DO
1346 :
1347 676327 : END SUBROUTINE halfill_upcorn
1348 :
1349 : ! **************************************************************************************************
1350 : !> \brief (Based on suitable modifications of S.Goedecker routines)
1351 : !> Assign the correct planes to the work array zmpi2
1352 : !> in order to prepare for interprocessor data transposition.
1353 : !> In the meanwhile, it unpacks the data of the HalFFT in order to prepare for
1354 : !> multiplication with the kernel
1355 : !> \param i1 Starting points of the plane and number of remaining lines
1356 : !> \param j2 Starting points of the plane and number of remaining lines
1357 : !> \param lot Starting points of the plane and number of remaining lines
1358 : !> \param nfft Starting points of the plane and number of remaining lines
1359 : !> \param n1 logical dimension of the FFT transform, reference for work arrays
1360 : !> \param n3 logical dimension of the FFT transform, reference for work arrays
1361 : !> \param md2 Dimensions of real grid
1362 : !> \param nproc ...
1363 : !> \param nd3 Dimensions of the kernel
1364 : !> \param zw Work array (input)
1365 : !> \param zmpi2 Work array for multiprocessor manipulation (output)
1366 : !> \param cosinarr Array of the phases needed for unpacking
1367 : !> \date February 2006
1368 : !> \author S. Goedecker, L. Genovese
1369 : !> \note
1370 : !> RESTRICTIONS on USAGE
1371 : !> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
1372 : !> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
1373 : !> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
1374 : !> This file is distributed under the terms of the
1375 : !> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
1376 : ! **************************************************************************************************
1377 781605 : SUBROUTINE scramble_unpack(i1, j2, lot, nfft, n1, n3, md2, nproc, nd3, zw, zmpi2, cosinarr)
1378 : INTEGER, INTENT(in) :: i1, j2, lot, nfft, n1, n3, md2, nproc, &
1379 : nd3
1380 : COMPLEX(KIND=dp), DIMENSION(lot, n3/2), INTENT(in) :: zw
1381 : COMPLEX(KIND=dp), DIMENSION(n1, md2/nproc, nd3), &
1382 : INTENT(inout) :: zmpi2
1383 : REAL(KIND=dp), DIMENSION(2, n3/2), INTENT(in) :: cosinarr
1384 :
1385 : INTEGER :: i, i3, ind1, ind2
1386 : REAL(KIND=dp) :: a, b, c, cp, d, feI, feR, fI, foI, foR, &
1387 : fR, sp
1388 :
1389 : !case i3=1 and i3=n3/2+1
1390 :
1391 23214805 : DO i = 0, nfft - 1
1392 22433200 : a = REAL(zw(i + 1, 1), dp)
1393 22433200 : b = AIMAG(zw(i + 1, 1))
1394 22433200 : zmpi2(i1 + i, j2, 1) = CMPLX(a + b, 0.0_dp, dp)
1395 23214805 : zmpi2(i1 + i, j2, n3/2 + 1) = CMPLX(a - b, 0.0_dp, dp)
1396 : END DO
1397 : !case 2<=i3<=n3/2
1398 42417300 : DO i3 = 2, n3/2
1399 41635695 : ind1 = i3
1400 41635695 : ind2 = n3/2 - i3 + 2
1401 41635695 : cp = cosinarr(1, i3)
1402 41635695 : sp = cosinarr(2, i3)
1403 1182671990 : DO i = 0, nfft - 1
1404 1140254690 : a = REAL(zw(i + 1, ind1), dp)
1405 1140254690 : b = AIMAG(zw(i + 1, ind1))
1406 1140254690 : c = REAL(zw(i + 1, ind2), dp)
1407 1140254690 : d = AIMAG(zw(i + 1, ind2))
1408 1140254690 : feR = .5_dp*(a + c)
1409 1140254690 : feI = .5_dp*(b - d)
1410 1140254690 : foR = .5_dp*(a - c)
1411 1140254690 : foI = .5_dp*(b + d)
1412 1140254690 : fR = feR + cp*foI - sp*foR
1413 1140254690 : fI = feI - cp*foR - sp*foI
1414 1181890385 : zmpi2(i1 + i, j2, ind1) = CMPLX(fR, fI, dp)
1415 : END DO
1416 : END DO
1417 :
1418 781605 : END SUBROUTINE scramble_unpack
1419 :
1420 : ! **************************************************************************************************
1421 : !> \brief (Based on suitable modifications of S.Goedecker routines)
1422 : !> Insert the correct planes of the work array zmpi2
1423 : !> in order to prepare for backward FFT transform
1424 : !> In the meanwhile, it packs the data in order to be transformed with the HalFFT
1425 : !> procedure
1426 : !> \param i1 Starting points of the plane and number of remaining lines
1427 : !> \param j2 Starting points of the plane and number of remaining lines
1428 : !> \param lot Starting points of the plane and number of remaining lines
1429 : !> \param nfft Starting points of the plane and number of remaining lines
1430 : !> \param n1 logical dimension of the FFT transform, reference for work arrays
1431 : !> \param n3 logical dimension of the FFT transform, reference for work arrays
1432 : !> \param md2 Dimensions of real grid
1433 : !> \param nproc ...
1434 : !> \param nd3 Dimensions of the kernel
1435 : !> \param zmpi2 Work array for multiprocessor manipulation (output)
1436 : !> \param zw Work array (inout)
1437 : !> \param cosinarr Array of the phases needed for packing
1438 : !> \date February 2006
1439 : !> \author S. Goedecker, L. Genovese
1440 : !> \note
1441 : !> RESTRICTIONS on USAGE
1442 : !> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
1443 : !> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
1444 : !> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
1445 : !> This file is distributed under the terms of the
1446 : !> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
1447 : ! **************************************************************************************************
1448 676327 : SUBROUTINE unscramble_pack(i1, j2, lot, nfft, n1, n3, md2, nproc, nd3, zmpi2, zw, cosinarr)
1449 : INTEGER, INTENT(in) :: i1, j2, lot, nfft, n1, n3, md2, nproc, &
1450 : nd3
1451 : COMPLEX(KIND=dp), DIMENSION(n1, md2/nproc, nd3), &
1452 : INTENT(in) :: zmpi2
1453 : COMPLEX(KIND=dp), DIMENSION(lot, n3/2), &
1454 : INTENT(inout) :: zw
1455 : REAL(KIND=dp), DIMENSION(2, n3/2), INTENT(in) :: cosinarr
1456 :
1457 : INTEGER :: i, i3, indA, indB
1458 : REAL(KIND=dp) :: a, b, c, cp, d, ie, ih, io, re, rh, ro, &
1459 : sp
1460 :
1461 36262611 : DO i3 = 1, n3/2
1462 35586284 : indA = i3
1463 35586284 : indB = n3/2 + 2 - i3
1464 35586284 : cp = cosinarr(1, i3)
1465 35586284 : sp = cosinarr(2, i3)
1466 1008815765 : DO i = 0, nfft - 1
1467 972553154 : a = REAL(zmpi2(i1 + i, j2, indA), dp)
1468 972553154 : b = AIMAG(zmpi2(i1 + i, j2, indA))
1469 972553154 : c = REAL(zmpi2(i1 + i, j2, indB), dp)
1470 972553154 : d = -AIMAG(zmpi2(i1 + i, j2, indB))
1471 972553154 : re = (a + c)
1472 972553154 : ie = (b + d)
1473 972553154 : ro = (a - c)*cp - (b - d)*sp
1474 972553154 : io = (a - c)*sp + (b - d)*cp
1475 972553154 : rh = re - io
1476 972553154 : ih = ie + ro
1477 1008139438 : zw(i + 1, indA) = CMPLX(rh, ih, dp)
1478 : END DO
1479 : END DO
1480 :
1481 676327 : END SUBROUTINE unscramble_pack
1482 :
1483 : ! **************************************************************************************************
1484 : !> \brief (Based on suitable modifications of S.Goedecker routines)
1485 : !> Applies the local FFT space Kernel to the density in Real space.
1486 : !> Calculates also the LDA exchange-correlation terms
1487 : !> \param n1 logical dimension of the transform.
1488 : !> \param n2 logical dimension of the transform.
1489 : !> \param n3 logical dimension of the transform.
1490 : !> \param nd1 Dimension of POT
1491 : !> \param nd2 Dimension of POT
1492 : !> \param nd3 Dimension of POT
1493 : !> \param md1 Dimension of ZF
1494 : !> \param md2 Dimension of ZF
1495 : !> \param md3 Dimension of ZF
1496 : !> \param nproc number of processors used as returned by MPI_COMM_SIZE
1497 : !> \param iproc [0:nproc-1] number of processor as returned by MPI_COMM_RANK
1498 : !> \param pot Kernel, only the distributed part (REAL)
1499 : !> POT(i1,i2,i3)
1500 : !> i1=1,nd1 , i2=1,nd2 , i3=1,nd3/nproc
1501 : !> \param zf Density (input/output)
1502 : !> ZF(i1,i3,i2)
1503 : !> i1=1,md1 , i2=1,md2/nproc , i3=1,md3
1504 : !> \param scal factor of renormalization of the FFT in order to acheve unitarity
1505 : !> and the correct dimension
1506 : !> \param mpi_group ...
1507 : !> \date February 2006
1508 : !> \author S. Goedecker, L. Genovese
1509 : !> \note
1510 : !> RESTRICTIONS on USAGE
1511 : !> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
1512 : !> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
1513 : !> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
1514 : !> This file is distributed under the terms of the
1515 : !> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
1516 : ! **************************************************************************************************
1517 16922 : SUBROUTINE F_PoissonSolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, pot, zf, &
1518 : scal, mpi_group)
1519 : INTEGER, INTENT(in) :: n1, n2, n3, nd1, nd2, nd3, md1, md2, &
1520 : md3, nproc, iproc
1521 : REAL(KIND=dp), DIMENSION(nd1, nd2, nd3/nproc), &
1522 : INTENT(in) :: pot
1523 : REAL(KIND=dp), DIMENSION(md1, md3, md2/nproc), &
1524 : INTENT(inout) :: zf
1525 : REAL(KIND=dp), INTENT(in) :: scal
1526 :
1527 : CLASS(mp_comm_type), INTENT(in) :: mpi_group
1528 :
1529 : INTEGER, PARAMETER :: ncache_optimal = 8*1024
1530 :
1531 : INTEGER :: i1, i3, j, j2, &
1532 : J2stb, J2stf, j3, Jp2stb, Jp2stf, lot1, lot2, lot3, &
1533 : lzt, ma, mb, ncache, nfft, stat, &
1534 : final_chunk_size1, final_chunk_size2, final_chunk_size3
1535 : REAL(kind=dp) :: twopion
1536 16922 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cosinarr
1537 16922 : COMPLEX(KIND=dp), POINTER, CONTIGUOUS, DIMENSION(:, :) :: zt
1538 16922 : COMPLEX(KIND=dp), POINTER, CONTIGUOUS, DIMENSION(:) :: zw1, zw2
1539 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: zmpi2
1540 : COMPLEX(KIND=dp), ALLOCATABLE, &
1541 16922 : DIMENSION(:, :, :, :) :: zmpi1
1542 : TYPE(fft_plan_type) :: fft_plan_bw3, fft_plan_bw3_last, fft_plan_fw3, fft_plan_fw3_last, &
1543 : fft_plan_bw1, fft_plan_bw1_last, fft_plan_fw1, fft_plan_fw1_last, &
1544 : fft_plan_bw2, fft_plan_bw2_last, fft_plan_fw2, fft_plan_fw2_last
1545 :
1546 0 : IF (MOD(n1, 2) /= 0) CPABORT("Parallel convolution:ERROR:n1")
1547 16922 : IF (MOD(n2, 2) /= 0) CPABORT("Parallel convolution:ERROR:n2")
1548 16922 : IF (MOD(n3, 2) /= 0) CPABORT("Parallel convolution:ERROR:n3")
1549 16922 : IF (nd1 < n1/2 + 1) CPABORT("Parallel convolution:ERROR:nd1")
1550 16922 : IF (nd2 < n2/2 + 1) CPABORT("Parallel convolution:ERROR:nd2")
1551 16922 : IF (nd3 < n3/2 + 1) CPABORT("Parallel convolution:ERROR:nd3")
1552 16922 : IF (md1 < n1/2) CPABORT("Parallel convolution:ERROR:md1")
1553 16922 : IF (md2 < n2/2) CPABORT("Parallel convolution:ERROR:md2")
1554 16922 : IF (md3 < n3/2) CPABORT("Parallel convolution:ERROR:md3")
1555 16922 : IF (MOD(nd3, nproc) /= 0) CPABORT("Parallel convolution:ERROR:nd3")
1556 16922 : IF (MOD(md2, nproc) /= 0) CPABORT("Parallel convolution:ERROR:md2")
1557 :
1558 : !defining work arrays dimensions
1559 :
1560 16922 : ncache = ncache_optimal
1561 16922 : IF (ncache <= MAX(n1, n2, n3/2)*4) ncache = MAX(n1, n2, n3/2)*4
1562 16922 : lzt = n2/2
1563 16922 : IF (MOD(n2/2, 2) == 0) lzt = lzt + 1
1564 16922 : IF (MOD(n2/2, 4) == 0) lzt = lzt + 1
1565 :
1566 : !Allocations
1567 33844 : CALL fft_alloc(zw1, [ncache/4])
1568 16922 : zw1 = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1569 33844 : CALL fft_alloc(zw2, [ncache/4])
1570 16922 : zw2 = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1571 50766 : CALL fft_alloc(zt, [lzt, n1])
1572 16922 : zt = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
1573 84610 : ALLOCATE (zmpi2(n1, md2/nproc, nd3), SOURCE=CMPLX(0.0_dp, 0.0_dp, dp))
1574 67688 : ALLOCATE (cosinarr(2, n3/2), SOURCE=0.0_dp)
1575 63702 : IF (nproc > 1) ALLOCATE (zmpi1(n1, md2/nproc, nd3/nproc, nproc), SOURCE=CMPLX(0.0_dp, 0.0_dp, dp))
1576 :
1577 : !Calculating array of phases for HalFFT decoding
1578 16922 : twopion = 8._dp*ATAN(1._dp)/REAL(n3, KIND=dp)
1579 698018 : DO i3 = 1, n3/2
1580 681096 : cosinarr(1, i3) = COS(twopion*(i3 - 1))
1581 698018 : cosinarr(2, i3) = -SIN(twopion*(i3 - 1))
1582 : END DO
1583 :
1584 : ! transform along z axis
1585 16922 : lot1 = ncache/(4*n1)
1586 16922 : lot2 = ncache/(4*n2)
1587 16922 : lot3 = ncache/(2*n3)
1588 :
1589 : ! The size of the last chunk
1590 16922 : final_chunk_size1 = MOD(n2/2, lot1)
1591 16922 : final_chunk_size2 = MOD(n1, lot2)
1592 16922 : final_chunk_size3 = MOD(n1/2, lot3)
1593 :
1594 : ! Prevent OOB-access if n2 < lot1
1595 16922 : IF (n2/2 >= lot1) THEN
1596 11623 : CALL fft_create_plan_1d(fft_plan_bw1, BWFFT, .TRUE., .TRUE., lot1, lzt, n1, lot1, zw1, zt)
1597 11623 : CALL fft_create_plan_1d(fft_plan_fw1, FWFFT, .TRUE., .TRUE., lzt, lot1, n1, lot1, zt, zw1)
1598 : END IF
1599 16922 : IF (final_chunk_size1 > 0) THEN
1600 : CALL fft_create_plan_1d(fft_plan_bw1_last, BWFFT, .TRUE., .TRUE., lot1, lzt, n1, &
1601 14071 : final_chunk_size1, zw1, zt)
1602 : CALL fft_create_plan_1d(fft_plan_fw1_last, FWFFT, .TRUE., .TRUE., lzt, lot1, n1, &
1603 14071 : final_chunk_size1, zt, zw1)
1604 : END IF
1605 :
1606 : ! Prevent OOB-access if n2 < lot1
1607 16922 : IF (n1 >= lot2) THEN
1608 14154 : CALL fft_create_plan_1d(fft_plan_bw2, BWFFT, .TRUE., .TRUE., lot2, lot2, n2, lot2, zw1, zw2)
1609 14154 : CALL fft_create_plan_1d(fft_plan_fw2, FWFFT, .TRUE., .TRUE., lot2, lot2, n2, lot2, zw2, zw1)
1610 : END IF
1611 16922 : IF (final_chunk_size2 > 0) THEN
1612 : CALL fft_create_plan_1d(fft_plan_bw2_last, BWFFT, .TRUE., .TRUE., lot2, lot2, n2, &
1613 14071 : final_chunk_size2, zw1, zw2)
1614 : CALL fft_create_plan_1d(fft_plan_fw2_last, FWFFT, .TRUE., .TRUE., lot2, lot2, n2, &
1615 14071 : final_chunk_size2, zw2, zw1)
1616 : END IF
1617 :
1618 : ! Prevent OOB-access if n1 < lot3
1619 16922 : IF (n1/2 >= lot3) THEN
1620 4976 : CALL fft_create_plan_1d(fft_plan_fw3, FWFFT, .TRUE., .TRUE., lot3, lot3, n3/2, lot3, zw1, zw2)
1621 4976 : CALL fft_create_plan_1d(fft_plan_bw3, BWFFT, .TRUE., .TRUE., lot3, lot3, n3/2, lot3, zw1, zw2)
1622 : END IF
1623 16922 : IF (final_chunk_size3 > 0) THEN
1624 : CALL fft_create_plan_1d(fft_plan_fw3_last, FWFFT, .TRUE., .TRUE., &
1625 15900 : lot3, lot3, n3/2, final_chunk_size3, zw1, zw2)
1626 : CALL fft_create_plan_1d(fft_plan_bw3_last, BWFFT, .TRUE., .TRUE., &
1627 15900 : lot3, lot3, n3/2, final_chunk_size3, zw1, zw2)
1628 : END IF
1629 :
1630 464202 : DO j2 = 1, md2/nproc
1631 : !this condition ensures that we manage only the interesting part for the FFT
1632 464202 : IF (iproc*(md2/nproc) + j2 <= n2/2) THEN
1633 1119135 : DO i1 = 1, (n1/2), lot3
1634 673411 : ma = i1
1635 673411 : mb = MIN(i1 + (lot3 - 1), (n1/2))
1636 673411 : nfft = mb - ma + 1
1637 :
1638 : !inserting real data into complex array of half length
1639 673411 : CALL halfill_upcorn(md1, md3, lot3, nfft, n3, zf(i1, 1, j2), zw1)
1640 :
1641 : !performing FFT
1642 : !input: I1,I3,J2,(Jp2)
1643 : ! This is equivalent to
1644 673411 : IF (nfft == lot3) THEN
1645 265846 : CALL fft_1d(fft_plan_bw3, zw1, zw2, 1.0_dp, stat)
1646 : ELSE
1647 407565 : CALL fft_1d(fft_plan_bw3_last, zw1, zw2, 1.0_dp, stat)
1648 : END IF
1649 : !output: I1,i3,J2,(Jp2)
1650 :
1651 : !unpacking FFT in order to restore correct result,
1652 : !while exchanging components
1653 : !input: I1,i3,J2,(Jp2)
1654 1119135 : CALL scramble_unpack(i1, j2, lot3, nfft, n1/2, n3, md2, nproc, nd3, zw2, zmpi2, cosinarr)
1655 : !output: I1,J2,i3,(Jp2)
1656 : END DO
1657 : END IF
1658 : END DO
1659 :
1660 : !Interprocessor data transposition
1661 : !input: I1,J2,j3,jp3,(Jp2)
1662 16922 : IF (nproc > 1) THEN
1663 : !communication scheduling
1664 9356 : CALL mpi_group%alltoall(zmpi2, zmpi1, n1/2*(md2/nproc)*(nd3/nproc))
1665 : END IF
1666 : !output: I1,J2,j3,Jp2,(jp3)
1667 :
1668 : !now each process perform complete convolution of its planes
1669 490112 : DO j3 = 1, nd3/nproc
1670 : !this condition ensures that we manage only the interesting part for the FFT
1671 490112 : IF (iproc*(nd3/nproc) + j3 <= n3/2 + 1) THEN
1672 468512 : Jp2stb = 1
1673 468512 : J2stb = 1
1674 468512 : Jp2stf = 1
1675 468512 : J2stf = 1
1676 :
1677 : ! transform along x axis
1678 :
1679 1656377 : DO j = 1, n2/2, lot1
1680 1187865 : ma = j
1681 1187865 : mb = MIN(j + (lot1 - 1), n2/2)
1682 1187865 : nfft = mb - ma + 1
1683 :
1684 : !reverse index ordering, leaving the planes to be transformed at the end
1685 : !input: I1,J2,j3,Jp2,(jp3)
1686 1187865 : IF (nproc == 1) THEN
1687 421221 : CALL mpiswitch_upcorn(j3, nfft, Jp2stb, J2stb, lot1, n1, md2, nd3, nproc, zmpi2, zw1)
1688 : ELSE
1689 766644 : CALL mpiswitch_upcorn(j3, nfft, Jp2stb, J2stb, lot1, n1, md2, nd3, nproc, zmpi1, zw1)
1690 : END IF
1691 : !output: J2,Jp2,I1,j3,(jp3)
1692 :
1693 : !performing FFT
1694 : !input: I2,I1,j3,(jp3)
1695 : ! This is equivalent to
1696 1656377 : IF (nfft == lot1) THEN
1697 805720 : CALL fft_1d(fft_plan_bw1, zw1, zt(j:, 1), 1.0_dp, stat)
1698 : ELSE
1699 382145 : CALL fft_1d(fft_plan_bw1_last, zw1, zt(j:, 1), 1.0_dp, stat)
1700 : END IF
1701 : !output: I2,i1,j3,(jp3)
1702 : END DO
1703 :
1704 : !transform along y axis
1705 :
1706 2711780 : DO j = 1, n1, lot2
1707 2243268 : ma = j
1708 2243268 : mb = MIN(j + (lot2 - 1), n1)
1709 2243268 : nfft = mb - ma + 1
1710 :
1711 : !reverse ordering
1712 : !input: I2,i1,j3,(jp3)
1713 2243268 : CALL switch_upcorn(nfft, n2, lot2, n1, lzt, zt(:, j), zw1)
1714 : !output: i1,I2,j3,(jp3)
1715 :
1716 : !performing FFT
1717 : !input: i1,I2,j3,(jp3)
1718 : ! This is equivalent to
1719 2243268 : IF (nfft == lot2) THEN
1720 1861123 : CALL fft_1d(fft_plan_bw2, zw1, zw2, 1.0_dp, stat)
1721 : ELSE
1722 382145 : CALL fft_1d(fft_plan_bw2_last, zw1, zw2, 1.0_dp, stat)
1723 : END IF
1724 : !output: i1,i2,j3,(jp3)
1725 :
1726 : !Multiply with kernel in fourier space
1727 2243268 : CALL multkernel(nd1, nd2, n1, n2, lot2, nfft, j, pot(1, 1, j3), zw2)
1728 :
1729 : !TRANSFORM BACK IN REAL SPACE
1730 :
1731 : !transform along y axis
1732 : !input: i1,i2,j3,(jp3)
1733 : ! This is equivalent to
1734 2243268 : IF (nfft == lot2) THEN
1735 1861123 : CALL fft_1d(fft_plan_fw2, zw2, zw1, 1.0_dp, stat)
1736 : ELSE
1737 382145 : CALL fft_1d(fft_plan_fw2_last, zw2, zw1, 1.0_dp, stat)
1738 : END IF
1739 :
1740 : !reverse ordering
1741 : !input: i1,I2,j3,(jp3)
1742 2711780 : CALL unswitch_downcorn(nfft, n2, lot2, n1, lzt, zw1, zt(:, j))
1743 : !output: I2,i1,j3,(jp3)
1744 : END DO
1745 :
1746 : !transform along x axis
1747 : !input: I2,i1,j3,(jp3)
1748 1656377 : DO j = 1, n2/2, lot1
1749 1187865 : ma = j
1750 1187865 : mb = MIN(j + (lot1 - 1), n2/2)
1751 1187865 : nfft = mb - ma + 1
1752 :
1753 : !performing FFT
1754 : ! This is equivalent to
1755 1187865 : IF (nfft == lot1) THEN
1756 805720 : CALL fft_1d(fft_plan_fw1, zt(j:, 1), zw2, 1.0_dp, stat)
1757 : ELSE
1758 382145 : CALL fft_1d(fft_plan_fw1_last, zt(j:, 1), zw2, 1.0_dp, stat)
1759 : END IF
1760 : !output: I2,I1,j3,(jp3)
1761 :
1762 : !reverse ordering
1763 : !input: J2,Jp2,I1,j3,(jp3)
1764 1656377 : IF (nproc == 1) THEN
1765 421221 : CALL unmpiswitch_downcorn(j3, nfft, Jp2stf, J2stf, lot1, n1, md2, nd3, nproc, zw2, zmpi2)
1766 : ELSE
1767 766644 : CALL unmpiswitch_downcorn(j3, nfft, Jp2stf, J2stf, lot1, n1, md2, nd3, nproc, zw2, zmpi1)
1768 : END IF
1769 : ! output: I1,J2,j3,Jp2,(jp3)
1770 : END DO
1771 : END IF
1772 : END DO
1773 :
1774 : !Interprocessor data transposition
1775 : !input: I1,J2,j3,Jp2,(jp3)
1776 16922 : IF (nproc > 1) THEN
1777 : !communication scheduling
1778 9356 : CALL mpi_group%alltoall(zmpi1, zmpi2, n1/2*(md2/nproc)*(nd3/nproc))
1779 : !output: I1,J2,j3,jp3,(Jp2)
1780 : END IF
1781 :
1782 : !transform along z axis
1783 : !input: I1,J2,i3,(Jp2)
1784 464202 : DO j2 = 1, md2/nproc
1785 : !this condition ensures that we manage only the interesting part for the FFT
1786 464202 : IF (iproc*(md2/nproc) + j2 <= n2/2) THEN
1787 1119135 : DO i1 = 1, (n1/2), lot3
1788 673411 : ma = i1
1789 673411 : mb = MIN(i1 + (lot3 - 1), (n1/2))
1790 673411 : nfft = mb - ma + 1
1791 :
1792 : !reverse ordering and repack the FFT data in order to be backward HalFFT transformed
1793 : !input: I1,J2,i3,(Jp2)
1794 673411 : CALL unscramble_pack(i1, j2, lot3, nfft, n1/2, n3, md2, nproc, nd3, zmpi2, zw1, cosinarr)
1795 : !output: I1,i3,J2,(Jp2)
1796 :
1797 : !performing FFT
1798 : !input: I1,i3,J2,(Jp2)
1799 : ! This is equivalent to
1800 673411 : IF (nfft == lot3) THEN
1801 265846 : CALL fft_1d(fft_plan_fw3, zw1, zw2, 1.0_dp, stat)
1802 : ELSE
1803 407565 : CALL fft_1d(fft_plan_fw3_last, zw1, zw2, 1.0_dp, stat)
1804 : END IF
1805 : !output: I1,I3,J2,(Jp2)
1806 :
1807 : !calculates the exchange correlation terms locally and rebuild the output array
1808 1119135 : CALL unfill_downcorn(md1, md3, lot3, nfft, n3, zw2, zf(i1, 1, j2), scal)
1809 : END DO
1810 : END IF
1811 : END DO
1812 :
1813 16922 : IF (n2/2 >= lot1) CALL fft_destroy_plan(fft_plan_bw1)
1814 16922 : IF (final_chunk_size1 > 0) CALL fft_destroy_plan(fft_plan_bw1_last)
1815 16922 : IF (n2/2 >= lot1) CALL fft_destroy_plan(fft_plan_fw1)
1816 16922 : IF (final_chunk_size1 > 0) CALL fft_destroy_plan(fft_plan_fw1_last)
1817 :
1818 16922 : IF (n1 >= lot2) CALL fft_destroy_plan(fft_plan_bw2)
1819 16922 : IF (final_chunk_size2 > 0) CALL fft_destroy_plan(fft_plan_bw2_last)
1820 16922 : IF (n1 >= lot2) CALL fft_destroy_plan(fft_plan_fw2)
1821 16922 : IF (final_chunk_size2 > 0) CALL fft_destroy_plan(fft_plan_fw2_last)
1822 :
1823 16922 : IF (n1/2 >= lot3) CALL fft_destroy_plan(fft_plan_fw3)
1824 16922 : IF (final_chunk_size3 > 0) CALL fft_destroy_plan(fft_plan_fw3_last)
1825 16922 : IF (n1/2 >= lot3) CALL fft_destroy_plan(fft_plan_bw3)
1826 16922 : IF (final_chunk_size3 > 0) CALL fft_destroy_plan(fft_plan_bw3_last)
1827 :
1828 : !De-allocations
1829 16922 : DEALLOCATE (zmpi2)
1830 16922 : CALL fft_dealloc(zw1)
1831 16922 : CALL fft_dealloc(zw2)
1832 16922 : CALL fft_dealloc(zt)
1833 16922 : DEALLOCATE (cosinarr)
1834 16922 : IF (nproc > 1) DEALLOCATE (zmpi1)
1835 :
1836 643036 : END SUBROUTINE F_PoissonSolver
1837 :
1838 : ! **************************************************************************************************
1839 : !> \brief ...
1840 : !> \param nfft ...
1841 : !> \param n2 ...
1842 : !> \param lot ...
1843 : !> \param n1 ...
1844 : !> \param lzt ...
1845 : !> \param zt ...
1846 : !> \param zw ...
1847 : ! **************************************************************************************************
1848 2243268 : PURE SUBROUTINE switch_upcorn(nfft, n2, lot, n1, lzt, zt, zw)
1849 : INTEGER, INTENT(IN) :: nfft, n2, lot, n1, lzt
1850 : COMPLEX(KIND=dp), INTENT(IN) :: zt(lzt, n1)
1851 : COMPLEX(KIND=dp), INTENT(INOUT) :: zw(lot, n2)
1852 :
1853 : INTEGER :: i, j
1854 :
1855 : ! WARNING: Assuming that high frequencies are in the corners
1856 : ! and that n2 is multiple of 2
1857 : ! Low frequencies
1858 :
1859 42222968 : DO j = 1, nfft
1860 2017044556 : DO i = n2/2 + 1, n2
1861 2014801288 : zw(j, i) = zt(i - n2/2, j)
1862 : END DO
1863 : END DO
1864 : ! High frequencies
1865 125580240 : DO i = 1, n2/2
1866 2100401828 : DO j = 1, nfft
1867 2098158560 : zw(j, i) = CMPLX(0.0_dp, 0.0_dp, dp)
1868 : END DO
1869 : END DO
1870 2243268 : END SUBROUTINE switch_upcorn
1871 :
1872 : ! **************************************************************************************************
1873 : !> \brief ...
1874 : !> \param j3 ...
1875 : !> \param nfft ...
1876 : !> \param Jp2stb ...
1877 : !> \param J2stb ...
1878 : !> \param lot ...
1879 : !> \param n1 ...
1880 : !> \param md2 ...
1881 : !> \param nd3 ...
1882 : !> \param nproc ...
1883 : !> \param zmpi1 ...
1884 : !> \param zw ...
1885 : ! **************************************************************************************************
1886 1187865 : PURE SUBROUTINE mpiswitch_upcorn(j3, nfft, Jp2stb, J2stb, lot, n1, md2, nd3, nproc, zmpi1, zw)
1887 : INTEGER, INTENT(IN) :: j3, nfft
1888 : INTEGER, INTENT(INOUT) :: Jp2stb, J2stb
1889 : INTEGER, INTENT(IN) :: lot, n1, md2, nd3, nproc
1890 : COMPLEX(KIND=dp), INTENT(IN) :: zmpi1(n1/2, md2/nproc, nd3/nproc, nproc)
1891 : COMPLEX(KIND=dp), INTENT(INOUT) :: zw(lot, n1)
1892 :
1893 : INTEGER :: i1, j2, jp2, mfft
1894 :
1895 : ! WARNING: Assuming that high frequencies are in the corners
1896 : ! and that n1 is multiple of 2
1897 :
1898 1187865 : mfft = 0
1899 1808671 : Main: DO Jp2 = Jp2stb, nproc
1900 21407221 : DO J2 = J2stb, md2/nproc
1901 20786415 : mfft = mfft + 1
1902 20786415 : IF (mfft > nfft) THEN
1903 796565 : Jp2stb = Jp2
1904 796565 : J2stb = J2
1905 796565 : EXIT Main
1906 : END IF
1907 1007400644 : DO I1 = 1, n1/2
1908 1007400644 : zw(mfft, I1) = CMPLX(0.0_dp, 0.0_dp, dp)
1909 : END DO
1910 1008021450 : DO I1 = n1/2 + 1, n1
1911 1007400644 : zw(mfft, I1) = zmpi1(I1 - n1/2, J2, j3, Jp2)
1912 : END DO
1913 : END DO
1914 1012106 : J2stb = 1
1915 : END DO Main
1916 1187865 : END SUBROUTINE mpiswitch_upcorn
1917 :
1918 : ! **************************************************************************************************
1919 : !> \brief ...
1920 : !> \param nfft ...
1921 : !> \param n2 ...
1922 : !> \param lot ...
1923 : !> \param n1 ...
1924 : !> \param lzt ...
1925 : !> \param zw ...
1926 : !> \param zt ...
1927 : ! **************************************************************************************************
1928 2243268 : PURE SUBROUTINE unswitch_downcorn(nfft, n2, lot, n1, lzt, zw, zt)
1929 : INTEGER, INTENT(IN) :: nfft, n2, lot, n1, lzt
1930 : COMPLEX(KIND=dp), INTENT(IN) :: zw(lot, n2)
1931 : COMPLEX(KIND=dp), INTENT(INOUT) :: zt(lzt, n1)
1932 :
1933 : INTEGER :: i, j
1934 :
1935 : ! WARNING: Assuming that high frequencies are in the corners
1936 : ! and that n2 is multiple of 2
1937 : ! Low frequencies
1938 :
1939 42222968 : DO j = 1, nfft
1940 2017044556 : DO i = 1, n2/2
1941 2014801288 : zt(i, j) = zw(j, i)
1942 : END DO
1943 : END DO
1944 2243268 : RETURN
1945 : END SUBROUTINE unswitch_downcorn
1946 :
1947 : ! **************************************************************************************************
1948 : !> \brief ...
1949 : !> \param j3 ...
1950 : !> \param nfft ...
1951 : !> \param Jp2stf ...
1952 : !> \param J2stf ...
1953 : !> \param lot ...
1954 : !> \param n1 ...
1955 : !> \param md2 ...
1956 : !> \param nd3 ...
1957 : !> \param nproc ...
1958 : !> \param zw ...
1959 : !> \param zmpi1 ...
1960 : ! **************************************************************************************************
1961 1187865 : PURE SUBROUTINE unmpiswitch_downcorn(j3, nfft, Jp2stf, J2stf, lot, n1, md2, nd3, nproc, zw, zmpi1)
1962 : INTEGER, INTENT(IN) :: j3, nfft
1963 : INTEGER, INTENT(INOUT) :: Jp2stf, J2stf
1964 : INTEGER, INTENT(IN) :: lot, n1, md2, nd3, nproc
1965 : COMPLEX(KIND=dp), INTENT(IN) :: zw(lot, n1)
1966 : COMPLEX(KIND=dp), INTENT(INOUT) :: zmpi1(n1/2, md2/nproc, nd3/nproc, nproc)
1967 :
1968 : INTEGER :: i1, j2, jp2, mfft
1969 :
1970 : ! WARNING: Assuming that high frequencies are in the corners
1971 : ! and that n1 is multiple of 2
1972 :
1973 1187865 : mfft = 0
1974 1808671 : Main: DO Jp2 = Jp2stf, nproc
1975 21407221 : DO J2 = J2stf, md2/nproc
1976 20786415 : mfft = mfft + 1
1977 20786415 : IF (mfft > nfft) THEN
1978 796565 : Jp2stf = Jp2
1979 796565 : J2stf = J2
1980 796565 : EXIT Main
1981 : END IF
1982 1008021450 : DO I1 = 1, n1/2
1983 1007400644 : zmpi1(I1, J2, j3, Jp2) = zw(mfft, I1)
1984 : END DO
1985 : END DO
1986 1012106 : J2stf = 1
1987 : END DO Main
1988 1187865 : END SUBROUTINE unmpiswitch_downcorn
1989 :
1990 : ! **************************************************************************************************
1991 : !> \brief (Based on suitable modifications of S.Goedecker routines)
1992 : !> Restore data into output array, calculating in the meanwhile
1993 : !> Hartree energy of the potential
1994 : !> \param md1 Dimensions of the undistributed part of the real grid
1995 : !> \param md3 Dimensions of the undistributed part of the real grid
1996 : !> \param lot ...
1997 : !> \param nfft number of planes
1998 : !> \param n3 (twice the) dimension of the last FFTtransform.
1999 : !> \param zw FFT work array
2000 : !> \param zf Original distributed density as well as
2001 : !> Distributed solution of the poisson equation (inout)
2002 : !> \param scal Needed to achieve unitarity and correct dimensions
2003 : !> \param ehartreetmp Hartree energy
2004 : !> \date February 2006
2005 : !> \author S. Goedecker, L. Genovese
2006 : !> \note Assuming that high frequencies are in the corners
2007 : !> and that n3 is multiple of 4
2008 : !>
2009 : !> RESTRICTIONS on USAGE
2010 : !> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
2011 : !> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
2012 : !> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
2013 : !> This file is distributed under the terms of the
2014 : !> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
2015 : ! **************************************************************************************************
2016 0 : PURE SUBROUTINE F_unfill_downcorn(md1, md3, lot, nfft, n3, zw, zf, scal, ehartreetmp)
2017 : INTEGER, INTENT(in) :: md1, md3, lot, nfft, n3
2018 : COMPLEX(KIND=dp), DIMENSION(lot, n3/2), INTENT(in) :: zw
2019 : REAL(KIND=dp), DIMENSION(md1, md3), INTENT(inout) :: zf
2020 : REAL(KIND=dp), INTENT(in) :: scal
2021 : REAL(KIND=dp), INTENT(out) :: ehartreetmp
2022 :
2023 : INTEGER :: i1, i3
2024 : REAL(KIND=dp) :: pot1
2025 :
2026 0 : ehartreetmp = 0._dp
2027 0 : DO i3 = 1, n3/4
2028 0 : DO i1 = 1, nfft
2029 0 : pot1 = scal*REAL(zw(i1, i3), dp)
2030 0 : ehartreetmp = ehartreetmp + pot1*zf(i1, 2*i3 - 1)
2031 0 : zf(i1, 2*i3 - 1) = pot1
2032 0 : pot1 = scal*AIMAG(zw(i1, i3))
2033 0 : ehartreetmp = ehartreetmp + pot1*zf(i1, 2*i3)
2034 0 : zf(i1, 2*i3) = pot1
2035 : END DO
2036 : END DO
2037 0 : END SUBROUTINE F_unfill_downcorn
2038 :
2039 : END MODULE ps_wavelet_base
|
