(git:5e7fe52)
Loading...
Searching...
No Matches
ps_wavelet_base.F
Go to the documentation of this file.
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! **************************************************************************************************
12MODULE 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
34CONTAINS
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 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 COMPLEX(KIND=dp), POINTER, CONTIGUOUS, DIMENSION(:, :) :: zt
73 COMPLEX(KIND=dp), POINTER, CONTIGUOUS, DIMENSION(:) :: zw1, zw2
74 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: zmpi2
75 COMPLEX(KIND=dp), ALLOCATABLE, &
76 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 IF (nd1 < n1/2 + 1) cpabort("Parallel convolution:ERROR:nd1")
82 IF (nd2 < n2/2 + 1) cpabort("Parallel convolution:ERROR:nd2")
83 IF (nd3 < n3/2 + 1) cpabort("Parallel convolution:ERROR:nd3")
84 IF (md1 < n1) cpabort("Parallel convolution:ERROR:md1")
85 IF (md2 < n2) cpabort("Parallel convolution:ERROR:md2")
86 IF (md3 < n3) cpabort("Parallel convolution:ERROR:md3")
87 IF (mod(nd3, nproc) /= 0) cpabort("Parallel convolution:ERROR:nd3")
88 IF (mod(md2, nproc) /= 0) cpabort("Parallel convolution:ERROR:md2")
89
90 !defining work arrays dimensions
91 ncache = ncache_optimal
92 IF (ncache <= max(n1, n2, n3)*4) ncache = max(n1, n2, n3)*4
93
94 lzt = n2
95 IF (mod(n2, 2) == 0) lzt = lzt + 1
96
97 !Allocations
98 CALL fft_alloc(zw1, [ncache/4])
99 zw1 = cmplx(0.0_dp, 0.0_dp, kind=dp)
100 CALL fft_alloc(zw2, [ncache/4])
101 zw2 = cmplx(0.0_dp, 0.0_dp, kind=dp)
102 CALL fft_alloc(zt, [lzt, n1])
103 zt = cmplx(0.0_dp, 0.0_dp, kind=dp)
104 ALLOCATE (zmpi2(n1, md2/nproc, nd3), source=cmplx(0.0_dp, 0.0_dp, kind=dp))
105 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 lot1 = ncache/(4*n1)
110 lot2 = ncache/(4*n2)
111 lot3 = ncache/(4*n3)
112
113 ! The size of the last chunk
114 final_chunk_size1 = mod(n2, lot1)
115 final_chunk_size2 = mod(n1, lot2)
116 final_chunk_size3 = mod(n1, lot3)
117
118 ! Prevent OOB-access if n2 < lot1
119 IF (n2 >= lot1) THEN
120 CALL fft_create_plan_1d(fft_plan_bw1, bwfft, .true., .true., lot1, lzt, n1, lot1, zw1, zt)
121 CALL fft_create_plan_1d(fft_plan_fw1, fwfft, .true., .true., lzt, lot1, n1, lot1, zt, zw1)
122 END IF
123 IF (final_chunk_size1 > 0) THEN
124 CALL fft_create_plan_1d(fft_plan_bw1_last, bwfft, .true., .true., lot1, lzt, n1, &
125 final_chunk_size1, zw1, zt)
126 CALL fft_create_plan_1d(fft_plan_fw1_last, fwfft, .true., .true., lzt, lot1, n1, &
127 final_chunk_size1, zt, zw1)
128 END IF
129
130 ! Prevent OOB-access if n2 < lot1
131 IF (n1 >= lot2) THEN
132 CALL fft_create_plan_1d(fft_plan_bw2, bwfft, .true., .true., lot2, lot2, n2, lot2, zw1, zw2)
133 CALL fft_create_plan_1d(fft_plan_fw2, fwfft, .true., .true., lot2, lot2, n2, lot2, zw2, zw1)
134 END IF
135 IF (final_chunk_size2 > 0) THEN
136 CALL fft_create_plan_1d(fft_plan_bw2_last, bwfft, .true., .true., lot2, lot2, n2, &
137 final_chunk_size2, zw1, zw2)
138 CALL fft_create_plan_1d(fft_plan_fw2_last, fwfft, .true., .true., lot2, lot2, n2, &
139 final_chunk_size2, zw2, zw1)
140 END IF
141
142 ! Prevent OOB-access if n1 < lot3
143 IF (n1 >= lot3) THEN
144 CALL fft_create_plan_1d(fft_plan_fw3, fwfft, .true., .true., lot3, lot3, n3, lot3, zw1, zw2)
145 CALL fft_create_plan_1d(fft_plan_bw3, bwfft, .true., .true., lot3, lot3, n3, lot3, zw1, zw2)
146 END IF
147 IF (final_chunk_size3 > 0) THEN
148 CALL fft_create_plan_1d(fft_plan_fw3_last, fwfft, .true., .true., &
149 lot3, lot3, n3, final_chunk_size3, zw1, zw2)
150 CALL fft_create_plan_1d(fft_plan_bw3_last, bwfft, .true., .true., &
151 lot3, lot3, n3, final_chunk_size3, zw1, zw2)
152 END IF
153
154 DO j2 = 1, md2/nproc
155 !this condition ensures that we manage only the interesting part for the FFT
156 IF (iproc*(md2/nproc) + j2 <= n2) THEN
157 DO i1 = 1, n1, lot3
158 ma = i1
159 mb = min(i1 + (lot3 - 1), n1)
160 nfft = mb - ma + 1
161 !inserting real data into complex array of half length
162 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 IF (nfft == lot3) THEN
168 CALL fft_1d(fft_plan_bw3, zw1, zw2, 1.0_dp, stat)
169 ELSE
170 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 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 IF (nproc > 1) THEN
185 !communication scheduling
186 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 DO j3 = 1, nd3/nproc
192 !this condition ensures that we manage only the interesting part for the FFT
193 IF (iproc*(nd3/nproc) + j3 <= n3/2 + 1) THEN
194 jp2stb = 1
195 j2stb = 1
196 jp2stf = 1
197 j2stf = 1
198
199 ! transform along x axis
200
201 DO j = 1, n2, lot1
202 ma = j
203 mb = min(j + (lot1 - 1), n2)
204 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 IF (nproc == 1) THEN
209 CALL p_mpiswitch_upcorn(j3, nfft, jp2stb, j2stb, lot1, n1, md2, nd3, nproc, zmpi2, zw1)
210 ELSE
211 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 IF (nfft == lot1) THEN
219 CALL fft_1d(fft_plan_bw1, zw1, zt(j:, 1), 1.0_dp, stat)
220 ELSE
221 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 DO j = 1, n1, lot2
229 ma = j
230 mb = min(j + (lot2 - 1), n1)
231 nfft = mb - ma + 1
232
233 !reverse ordering
234 !input: I2,i1,j3,(jp3)
235 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 IF (nfft == lot2) THEN
242 CALL fft_1d(fft_plan_bw2, zw1, zw2, 1.0_dp, stat)
243 ELSE
244 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 i3 = iproc*(nd3/nproc) + j3
250 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 IF (nfft == lot2) THEN
258 CALL fft_1d(fft_plan_fw2, zw2, zw1, 1.0_dp, stat)
259 ELSE
260 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 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 DO j = 1, n2, lot1
272 ma = j
273 mb = min(j + (lot1 - 1), n2)
274 nfft = mb - ma + 1
275
276 !performing FFT
277 ! This is equivalent to
278 IF (nfft == lot1) THEN
279 CALL fft_1d(fft_plan_fw1, zt(j:, 1), zw2, 1.0_dp, stat)
280 ELSE
281 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 IF (nproc == 1) THEN
288 CALL p_unmpiswitch_downcorn(j3, nfft, jp2stf, j2stf, lot1, n1, md2, nd3, nproc, zw2, zmpi2)
289 ELSE
290 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 IF (nproc > 1) THEN
300 !communication scheduling
301 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 DO j2 = 1, md2/nproc
307 !this condition ensures that we manage only the interesting part for the FFT
308 IF (iproc*(md2/nproc) + j2 <= n2) THEN
309 DO i1 = 1, n1, lot3
310 ma = i1
311 mb = min(i1 + (lot3 - 1), n1)
312 nfft = mb - ma + 1
313
314 !reverse ordering
315 !input: I1,J2,i3,(Jp2)
316 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 IF (nfft == lot3) THEN
323 CALL fft_1d(fft_plan_fw3, zw1, zw2, 1.0_dp, stat)
324 ELSE
325 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 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 IF (n2 >= lot1) CALL fft_destroy_plan(fft_plan_bw1)
337 IF (final_chunk_size1 > 0) CALL fft_destroy_plan(fft_plan_bw1_last)
338 IF (n2 >= lot1) CALL fft_destroy_plan(fft_plan_fw1)
339 IF (final_chunk_size1 > 0) CALL fft_destroy_plan(fft_plan_fw1_last)
340
341 IF (n1 >= lot2) CALL fft_destroy_plan(fft_plan_bw2)
342 IF (final_chunk_size2 > 0) CALL fft_destroy_plan(fft_plan_bw2_last)
343 IF (n1 >= lot2) CALL fft_destroy_plan(fft_plan_fw2)
344 IF (final_chunk_size2 > 0) CALL fft_destroy_plan(fft_plan_fw2_last)
345
346 IF (n1 >= lot3) CALL fft_destroy_plan(fft_plan_fw3)
347 IF (final_chunk_size3 > 0) CALL fft_destroy_plan(fft_plan_fw3_last)
348 IF (n1 >= lot3) CALL fft_destroy_plan(fft_plan_bw3)
349 IF (final_chunk_size3 > 0) CALL fft_destroy_plan(fft_plan_bw3_last)
350
351 !De-allocations
352 DEALLOCATE (zmpi2)
353 CALL fft_dealloc(zw1)
354 CALL fft_dealloc(zw2)
355 CALL fft_dealloc(zt)
356 IF (nproc > 1) DEALLOCATE (zmpi1)
357 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 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 mfft = 0
386 DO jp2 = jp2stb, nproc
387 DO j2 = j2stb, md2/nproc
388 mfft = mfft + 1
389 IF (mfft > nfft) THEN
390 jp2stb = jp2
391 j2stb = j2
392 RETURN
393 END IF
394 DO i1 = 1, n1
395 zw(mfft, i1) = zmpi1(i1, j2, j3, jp2)
396 END DO
397 END DO
398 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 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 DO j = 1, nfft
421 DO i = 1, n2
422 zw(j, i) = zt(i, j)
423 END DO
424 END DO
425
426 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 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 DO j = 1, nfft
447 DO i = 1, n2
448 zt(i, j) = zw(j, i)
449 END DO
450 END DO
451
452 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 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 mfft = 0
480 DO jp2 = jp2stf, nproc
481 DO j2 = j2stf, md2/nproc
482 mfft = mfft + 1
483 IF (mfft > nfft) THEN
484 jp2stf = jp2
485 j2stf = j2
486 RETURN
487 END IF
488 DO i1 = 1, n1
489 zmpi1(i1, j2, j3, jp2) = zw(mfft, i1)
490 END DO
491 END DO
492 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 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 DO i3 = 1, n3
530 DO i1 = 1, nfft
531 pot1 = scal*real(zw(i1, i3), dp)
532 zf(i1, i3) = pot1
533 END DO
534 END DO
535
536 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 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 DO i3 = 1, n3
557 DO i1 = 1, nfft
558 zw(i1, i3) = cmplx(zf(i1, i3), 0.0_dp, dp)
559 END DO
560 END DO
561
562 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 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 DO i3 = 1, n3/2 + 1
599 DO i = 0, nfft - 1
600 zmpi2(i1 + i, j2, i3) = zw(i + 1, i3)
601 END DO
602 END DO
603
604 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 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 i3 = 1
642 DO i = 0, nfft - 1
643 zw(i + 1, i3) = zmpi2(i1 + i, j2, i3)
644 END DO
645
646 DO i3 = 2, n3/2 + 1
647 j3 = n3 + 2 - i3
648 DO i = 0, nfft - 1
649 zw(i + 1, j3) = conjg(zmpi2(i1 + i, j2, i3))
650 zw(i + 1, i3) = zmpi2(i1 + i, j2, i3)
651 END DO
652 END DO
653
654 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 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 fourpi2 = 4._dp*pi**2
694 j3 = i3 !n3/2+1-abs(n3/2+2-i3)
695 mu3 = real(j3 - 1, kind=dp)/real(n3, kind=dp)
696 mu3 = (mu3/hy)**2 !beware of the exchanged dimension
697 !Body
698 !generic case
699 DO i2 = 1, n2
700 DO i1 = 1, nfft
701 j1 = i1 + js - 1
702 j1 = j1 - (j1/(n1/2 + 2))*n1 !n1/2+1-abs(n1/2+2-jS-i1)
703 j2 = i2 - (i2/(n2/2 + 2))*n2 !n2/2+1-abs(n2/2+1-i2)
704 p1 = real(j1 - 1, kind=dp)/real(n1, kind=dp)
705 p2 = real(j2 - 1, kind=dp)/real(n2, kind=dp)
706 ker = -fourpi2*((p1/hx)**2 + (p2/hz)**2 + mu3) !beware of the exchanged dimension
707 IF (ker /= 0._dp) ker = 1._dp/ker
708 zw(i1, i2) = zw(i1, i2)*ker
709 END DO
710 END DO
711
712 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 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 DO j = 1, nfft
749 j1 = j + js - 1
750 j1 = j1 + (j1/(n1/2 + 2))*(n1 + 2 - 2*j1)
751 zw(j, 1) = zw(j, 1)*pot(j1, 1)
752 END DO
753
754 !generic case
755 DO i2 = 2, n2/2
756 DO j = 1, nfft
757 j1 = j + js - 1
758 j1 = j1 + (j1/(n1/2 + 2))*(n1 + 2 - 2*j1)
759 j2 = n2 + 2 - i2
760 zw(j, i2) = zw(j, i2)*pot(j1, i2)
761 zw(j, j2) = zw(j, j2)*pot(j1, i2)
762 END DO
763 END DO
764
765 !case i2=n2/2+1
766 DO j = 1, nfft
767 j1 = j + js - 1
768 j1 = j1 + (j1/(n1/2 + 2))*(n1 + 2 - 2*j1)
769 j2 = n2/2 + 1
770 zw(j, j2) = zw(j, j2)*pot(j1, j2)
771 END DO
772
773 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cosinarr
831 COMPLEX(KIND=dp), POINTER, CONTIGUOUS, DIMENSION(:, :) :: zt
832 COMPLEX(KIND=dp), POINTER, CONTIGUOUS, DIMENSION(:) :: zw1, zw2
833 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: zmpi2
834 COMPLEX(KIND=dp), ALLOCATABLE, &
835 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 CALL timeset(routinen, handle)
841 ! check input
842 IF (mod(n3, 2) /= 0) cpabort("Parallel convolution:ERROR:n3")
843 IF (nd1 < n1/2 + 1) cpabort("Parallel convolution:ERROR:nd1")
844 IF (nd2 < n2/2 + 1) cpabort("Parallel convolution:ERROR:nd2")
845 IF (nd3 < n3/2 + 1) cpabort("Parallel convolution:ERROR:nd3")
846 IF (md1 < n1) cpabort("Parallel convolution:ERROR:md1")
847 IF (md2 < n2) cpabort("Parallel convolution:ERROR:md2")
848 IF (md3 < n3/2) cpabort("Parallel convolution:ERROR:md3")
849 IF (mod(nd3, nproc) /= 0) cpabort("Parallel convolution:ERROR:nd3")
850 IF (mod(md2, nproc) /= 0) cpabort("Parallel convolution:ERROR:md2")
851
852 !defining work arrays dimensions
853 ncache = ncache_optimal
854 IF (ncache <= max(n1, n2, n3/2)*4) ncache = max(n1, n2, n3/2)*4
855
856 lzt = n2
857 IF (mod(n2, 2) == 0) lzt = lzt + 1
858 IF (mod(n2, 4) == 0) lzt = lzt + 1 !maybe this is useless
859
860 CALL fft_alloc(zw1, [ncache/4])
861 zw1 = cmplx(0.0_dp, 0.0_dp, kind=dp)
862 CALL fft_alloc(zw2, [ncache/4])
863 zw2 = cmplx(0.0_dp, 0.0_dp, kind=dp)
864 CALL fft_alloc(zt, [lzt, n1])
865 zt = cmplx(0.0_dp, 0.0_dp, dp)
866 ALLOCATE (zmpi2(n1, md2/nproc, nd3), source=cmplx(0.0_dp, 0.0_dp, dp))
867 ALLOCATE (cosinarr(2, n3/2))
868 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 twopion = 8._dp*atan(1._dp)/real(n3, kind=dp)
872 DO i3 = 1, n3/2
873 cosinarr(1, i3) = cos(twopion*(i3 - 1))
874 cosinarr(2, i3) = -sin(twopion*(i3 - 1))
875 END DO
876
877 !initializing integral
878
879 ! transform along z axis
880 lot1 = ncache/(4*n1)
881 lot2 = ncache/(4*n2)
882 lot3 = ncache/(2*n3)
883
884 ! The size of the last chunk
885 final_chunk_size1 = mod(n2, lot1)
886 final_chunk_size2 = mod(n1, lot2)
887 final_chunk_size3 = mod(n1, lot3)
888
889 ! Prevent OOB-access if n2 < lot1
890 IF (n2 >= lot1) THEN
891 CALL fft_create_plan_1d(fft_plan_bw1, bwfft, .true., .true., lot1, lzt, n1, lot1, zw1, zt)
892 CALL fft_create_plan_1d(fft_plan_fw1, fwfft, .true., .true., lzt, lot1, n1, lot1, zt, zw1)
893 END IF
894 IF (final_chunk_size1 > 0) THEN
895 CALL fft_create_plan_1d(fft_plan_bw1_last, bwfft, .true., .true., lot1, lzt, n1, &
896 final_chunk_size1, zw1, zt)
897 CALL fft_create_plan_1d(fft_plan_fw1_last, fwfft, .true., .true., lzt, lot1, n1, &
898 final_chunk_size1, zt, zw1)
899 END IF
900
901 ! Prevent OOB-access if n2 < lot1
902 IF (n1 >= lot2) THEN
903 CALL fft_create_plan_1d(fft_plan_bw2, bwfft, .true., .true., lot2, lot2, n2, lot2, zw1, zw2)
904 CALL fft_create_plan_1d(fft_plan_fw2, fwfft, .true., .true., lot2, lot2, n2, lot2, zw2, zw1)
905 END IF
906 IF (final_chunk_size2 > 0) THEN
907 CALL fft_create_plan_1d(fft_plan_bw2_last, bwfft, .true., .true., lot2, lot2, n2, &
908 final_chunk_size2, zw1, zw2)
909 CALL fft_create_plan_1d(fft_plan_fw2_last, fwfft, .true., .true., lot2, lot2, n2, &
910 final_chunk_size2, zw2, zw1)
911 END IF
912
913 ! Prevent OOB-access if n1 < lot3
914 IF (n1 >= lot3) THEN
915 CALL fft_create_plan_1d(fft_plan_fw3, fwfft, .true., .true., lot3, lot3, n3/2, lot3, zw1, zw2)
916 CALL fft_create_plan_1d(fft_plan_bw3, bwfft, .true., .true., lot3, lot3, n3/2, lot3, zw1, zw2)
917 END IF
918 IF (final_chunk_size3 > 0) THEN
919 CALL fft_create_plan_1d(fft_plan_fw3_last, fwfft, .true., .true., &
920 lot3, lot3, n3/2, final_chunk_size3, zw1, zw2)
921 CALL fft_create_plan_1d(fft_plan_bw3_last, bwfft, .true., .true., &
922 lot3, lot3, n3/2, final_chunk_size3, zw1, zw2)
923 END IF
924
925 DO j2 = 1, md2/nproc
926 !this condition ensures that we manage only the interesting part for the FFT
927 IF (iproc*(md2/nproc) + j2 <= n2) THEN
928 DO i1 = 1, n1, lot3
929 ma = i1
930 mb = min(i1 + (lot3 - 1), n1)
931 nfft = mb - ma + 1
932
933 !inserting real data into complex array of half length
934 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 IF (nfft == lot3) THEN
940 CALL fft_1d(fft_plan_bw3, zw1, zw2, 1.0_dp, stat)
941 ELSE
942 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 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 IF (nproc > 1) THEN
956 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 DO j3 = 1, nd3/nproc
962 !this condition ensures that we manage only the interesting part for the FFT
963 IF (iproc*(nd3/nproc) + j3 <= n3/2 + 1) THEN
964 jp2stb = 1
965 j2stb = 1
966 jp2stf = 1
967 j2stf = 1
968
969 ! transform along x axis
970
971 DO j = 1, n2, lot1
972 ma = j
973 mb = min(j + (lot1 - 1), n2)
974 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 IF (nproc == 1) THEN
979 CALL s_mpiswitch_upcorn(j3, nfft, jp2stb, j2stb, lot1, n1, md2, nd3, nproc, zmpi2, zw1)
980 ELSE
981 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 IF (nfft == lot1) THEN
989 CALL fft_1d(fft_plan_bw1, zw1, zt(j:, 1), 1.0_dp, stat)
990 ELSE
991 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 DO j = 1, n1, lot2
999 ma = j
1000 mb = min(j + (lot2 - 1), n1)
1001 nfft = mb - ma + 1
1002
1003 !reverse ordering
1004 !input: I2,i1,j3,(jp3)
1005 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 IF (nfft == lot2) THEN
1012 CALL fft_1d(fft_plan_bw2, zw1, zw2, 1.0_dp, stat)
1013 ELSE
1014 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 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 IF (nfft == lot2) THEN
1027 CALL fft_1d(fft_plan_fw2, zw2, zw1, 1.0_dp, stat)
1028 ELSE
1029 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 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 DO j = 1, n2, lot1
1041 ma = j
1042 mb = min(j + (lot1 - 1), n2)
1043 nfft = mb - ma + 1
1044
1045 !performing FFT
1046 ! This is equivalent to
1047 IF (nfft == lot1) THEN
1048 CALL fft_1d(fft_plan_fw1, zt(j:, 1), zw2, 1.0_dp, stat)
1049 ELSE
1050 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 IF (nproc == 1) THEN
1057 CALL s_unmpiswitch_downcorn(j3, nfft, jp2stf, j2stf, lot1, n1, md2, nd3, nproc, zw2, zmpi2)
1058 ELSE
1059 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 IF (nproc > 1) THEN
1069 !communication scheduling
1070 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 DO j2 = 1, md2/nproc
1078 !this condition ensures that we manage only the interesting part for the FFT
1079 IF (iproc*(md2/nproc) + j2 <= n2) THEN
1080 DO i1 = 1, n1, lot3
1081 ma = i1
1082 mb = min(i1 + (lot3 - 1), n1)
1083 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 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 IF (nfft == lot3) THEN
1094 CALL fft_1d(fft_plan_fw3, zw1, zw2, 1.0_dp, stat)
1095 ELSE
1096 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 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 IF (n2 >= lot1) CALL fft_destroy_plan(fft_plan_bw1)
1110 IF (final_chunk_size1 > 0) CALL fft_destroy_plan(fft_plan_bw1_last)
1111 IF (n2 >= lot1) CALL fft_destroy_plan(fft_plan_fw1)
1112 IF (final_chunk_size1 > 0) CALL fft_destroy_plan(fft_plan_fw1_last)
1113
1114 IF (n1 >= lot2) CALL fft_destroy_plan(fft_plan_bw2)
1115 IF (final_chunk_size2 > 0) CALL fft_destroy_plan(fft_plan_bw2_last)
1116 IF (n1 >= lot2) CALL fft_destroy_plan(fft_plan_fw2)
1117 IF (final_chunk_size2 > 0) CALL fft_destroy_plan(fft_plan_fw2_last)
1118
1119 IF (n1 >= lot3) CALL fft_destroy_plan(fft_plan_fw3)
1120 IF (final_chunk_size3 > 0) CALL fft_destroy_plan(fft_plan_fw3_last)
1121 IF (n1 >= lot3) CALL fft_destroy_plan(fft_plan_bw3)
1122 IF (final_chunk_size3 > 0) CALL fft_destroy_plan(fft_plan_bw3_last)
1123
1124 !De-allocations
1125 DEALLOCATE (zmpi2)
1126 CALL fft_dealloc(zw1)
1127 CALL fft_dealloc(zw2)
1128 CALL fft_dealloc(zt)
1129 DEALLOCATE (cosinarr)
1130 IF (nproc > 1) DEALLOCATE (zmpi1)
1131
1132 CALL timestop(handle)
1133 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 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 mfft = 0
1162 DO jp2 = jp2stb, nproc
1163 DO j2 = j2stb, md2/nproc
1164 mfft = mfft + 1
1165 IF (mfft > nfft) THEN
1166 jp2stb = jp2
1167 j2stb = j2
1168 RETURN
1169 END IF
1170 DO i1 = 1, n1
1171 zw(mfft, i1) = zmpi1(i1, j2, j3, jp2)
1172 END DO
1173 END DO
1174 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 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 DO j = 1, nfft
1197 DO i = 1, n2
1198 zw(j, i) = zt(i, j)
1199 END DO
1200 END DO
1201 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 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 DO j = 1, nfft
1222 DO i = 1, n2
1223 zt(i, j) = zw(j, i)
1224 END DO
1225 END DO
1226 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 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 mfft = 0
1254 DO jp2 = jp2stf, nproc
1255 DO j2 = j2stf, md2/nproc
1256 mfft = mfft + 1
1257 IF (mfft > nfft) THEN
1258 jp2stf = jp2
1259 j2stf = j2
1260 RETURN
1261 END IF
1262 DO i1 = 1, n1
1263 zmpi1(i1, j2, j3, jp2) = zw(mfft, i1)
1264 END DO
1265 END DO
1266 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 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 DO i3 = 1, n3/4
1305 DO i1 = 1, nfft
1306 pot1 = scal*real(zw(i1, i3), dp)
1307 !ehartreetmp =ehartreetmp + pot1* zf(i1,2*i3-1)
1308 zf(i1, 2*i3 - 1) = pot1
1309 pot1 = scal*aimag(zw(i1, i3))
1310 !ehartreetmp =ehartreetmp + pot1* zf(i1,2*i3)
1311 zf(i1, 2*i3) = pot1
1312 END DO
1313 END DO
1314 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 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 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 DO i1 = 1, nfft
1338 zw(i1, i3) = cmplx(0.0_dp, 0.0_dp, dp)
1339 END DO
1340 END DO
1341 DO i3 = n3/4 + 1, n3/2
1342 DO i1 = 1, nfft
1343 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 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 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 DO i = 0, nfft - 1
1392 a = real(zw(i + 1, 1), dp)
1393 b = aimag(zw(i + 1, 1))
1394 zmpi2(i1 + i, j2, 1) = cmplx(a + b, 0.0_dp, dp)
1395 zmpi2(i1 + i, j2, n3/2 + 1) = cmplx(a - b, 0.0_dp, dp)
1396 END DO
1397 !case 2<=i3<=n3/2
1398 DO i3 = 2, n3/2
1399 ind1 = i3
1400 ind2 = n3/2 - i3 + 2
1401 cp = cosinarr(1, i3)
1402 sp = cosinarr(2, i3)
1403 DO i = 0, nfft - 1
1404 a = real(zw(i + 1, ind1), dp)
1405 b = aimag(zw(i + 1, ind1))
1406 c = real(zw(i + 1, ind2), dp)
1407 d = aimag(zw(i + 1, ind2))
1408 fer = .5_dp*(a + c)
1409 fei = .5_dp*(b - d)
1410 for = .5_dp*(a - c)
1411 foi = .5_dp*(b + d)
1412 fr = fer + cp*foi - sp*for
1413 fi = fei - cp*for - sp*foi
1414 zmpi2(i1 + i, j2, ind1) = cmplx(fr, fi, dp)
1415 END DO
1416 END DO
1417
1418 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 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 DO i3 = 1, n3/2
1462 inda = i3
1463 indb = n3/2 + 2 - i3
1464 cp = cosinarr(1, i3)
1465 sp = cosinarr(2, i3)
1466 DO i = 0, nfft - 1
1467 a = real(zmpi2(i1 + i, j2, inda), dp)
1468 b = aimag(zmpi2(i1 + i, j2, inda))
1469 c = real(zmpi2(i1 + i, j2, indb), dp)
1470 d = -aimag(zmpi2(i1 + i, j2, indb))
1471 re = (a + c)
1472 ie = (b + d)
1473 ro = (a - c)*cp - (b - d)*sp
1474 io = (a - c)*sp + (b - d)*cp
1475 rh = re - io
1476 ih = ie + ro
1477 zw(i + 1, inda) = cmplx(rh, ih, dp)
1478 END DO
1479 END DO
1480
1481 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cosinarr
1537 COMPLEX(KIND=dp), POINTER, CONTIGUOUS, DIMENSION(:, :) :: zt
1538 COMPLEX(KIND=dp), POINTER, CONTIGUOUS, DIMENSION(:) :: zw1, zw2
1539 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: zmpi2
1540 COMPLEX(KIND=dp), ALLOCATABLE, &
1541 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 IF (mod(n1, 2) /= 0) cpabort("Parallel convolution:ERROR:n1")
1547 IF (mod(n2, 2) /= 0) cpabort("Parallel convolution:ERROR:n2")
1548 IF (mod(n3, 2) /= 0) cpabort("Parallel convolution:ERROR:n3")
1549 IF (nd1 < n1/2 + 1) cpabort("Parallel convolution:ERROR:nd1")
1550 IF (nd2 < n2/2 + 1) cpabort("Parallel convolution:ERROR:nd2")
1551 IF (nd3 < n3/2 + 1) cpabort("Parallel convolution:ERROR:nd3")
1552 IF (md1 < n1/2) cpabort("Parallel convolution:ERROR:md1")
1553 IF (md2 < n2/2) cpabort("Parallel convolution:ERROR:md2")
1554 IF (md3 < n3/2) cpabort("Parallel convolution:ERROR:md3")
1555 IF (mod(nd3, nproc) /= 0) cpabort("Parallel convolution:ERROR:nd3")
1556 IF (mod(md2, nproc) /= 0) cpabort("Parallel convolution:ERROR:md2")
1557
1558 !defining work arrays dimensions
1559
1560 ncache = ncache_optimal
1561 IF (ncache <= max(n1, n2, n3/2)*4) ncache = max(n1, n2, n3/2)*4
1562 lzt = n2/2
1563 IF (mod(n2/2, 2) == 0) lzt = lzt + 1
1564 IF (mod(n2/2, 4) == 0) lzt = lzt + 1
1565
1566 !Allocations
1567 CALL fft_alloc(zw1, [ncache/4])
1568 zw1 = cmplx(0.0_dp, 0.0_dp, kind=dp)
1569 CALL fft_alloc(zw2, [ncache/4])
1570 zw2 = cmplx(0.0_dp, 0.0_dp, kind=dp)
1571 CALL fft_alloc(zt, [lzt, n1])
1572 zt = cmplx(0.0_dp, 0.0_dp, kind=dp)
1573 ALLOCATE (zmpi2(n1, md2/nproc, nd3), source=cmplx(0.0_dp, 0.0_dp, dp))
1574 ALLOCATE (cosinarr(2, n3/2), source=0.0_dp)
1575 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 twopion = 8._dp*atan(1._dp)/real(n3, kind=dp)
1579 DO i3 = 1, n3/2
1580 cosinarr(1, i3) = cos(twopion*(i3 - 1))
1581 cosinarr(2, i3) = -sin(twopion*(i3 - 1))
1582 END DO
1583
1584 ! transform along z axis
1585 lot1 = ncache/(4*n1)
1586 lot2 = ncache/(4*n2)
1587 lot3 = ncache/(2*n3)
1588
1589 ! The size of the last chunk
1590 final_chunk_size1 = mod(n2/2, lot1)
1591 final_chunk_size2 = mod(n1, lot2)
1592 final_chunk_size3 = mod(n1/2, lot3)
1593
1594 ! Prevent OOB-access if n2 < lot1
1595 IF (n2/2 >= lot1) THEN
1596 CALL fft_create_plan_1d(fft_plan_bw1, bwfft, .true., .true., lot1, lzt, n1, lot1, zw1, zt)
1597 CALL fft_create_plan_1d(fft_plan_fw1, fwfft, .true., .true., lzt, lot1, n1, lot1, zt, zw1)
1598 END IF
1599 IF (final_chunk_size1 > 0) THEN
1600 CALL fft_create_plan_1d(fft_plan_bw1_last, bwfft, .true., .true., lot1, lzt, n1, &
1601 final_chunk_size1, zw1, zt)
1602 CALL fft_create_plan_1d(fft_plan_fw1_last, fwfft, .true., .true., lzt, lot1, n1, &
1603 final_chunk_size1, zt, zw1)
1604 END IF
1605
1606 ! Prevent OOB-access if n2 < lot1
1607 IF (n1 >= lot2) THEN
1608 CALL fft_create_plan_1d(fft_plan_bw2, bwfft, .true., .true., lot2, lot2, n2, lot2, zw1, zw2)
1609 CALL fft_create_plan_1d(fft_plan_fw2, fwfft, .true., .true., lot2, lot2, n2, lot2, zw2, zw1)
1610 END IF
1611 IF (final_chunk_size2 > 0) THEN
1612 CALL fft_create_plan_1d(fft_plan_bw2_last, bwfft, .true., .true., lot2, lot2, n2, &
1613 final_chunk_size2, zw1, zw2)
1614 CALL fft_create_plan_1d(fft_plan_fw2_last, fwfft, .true., .true., lot2, lot2, n2, &
1615 final_chunk_size2, zw2, zw1)
1616 END IF
1617
1618 ! Prevent OOB-access if n1 < lot3
1619 IF (n1/2 >= lot3) THEN
1620 CALL fft_create_plan_1d(fft_plan_fw3, fwfft, .true., .true., lot3, lot3, n3/2, lot3, zw1, zw2)
1621 CALL fft_create_plan_1d(fft_plan_bw3, bwfft, .true., .true., lot3, lot3, n3/2, lot3, zw1, zw2)
1622 END IF
1623 IF (final_chunk_size3 > 0) THEN
1624 CALL fft_create_plan_1d(fft_plan_fw3_last, fwfft, .true., .true., &
1625 lot3, lot3, n3/2, final_chunk_size3, zw1, zw2)
1626 CALL fft_create_plan_1d(fft_plan_bw3_last, bwfft, .true., .true., &
1627 lot3, lot3, n3/2, final_chunk_size3, zw1, zw2)
1628 END IF
1629
1630 DO j2 = 1, md2/nproc
1631 !this condition ensures that we manage only the interesting part for the FFT
1632 IF (iproc*(md2/nproc) + j2 <= n2/2) THEN
1633 DO i1 = 1, (n1/2), lot3
1634 ma = i1
1635 mb = min(i1 + (lot3 - 1), (n1/2))
1636 nfft = mb - ma + 1
1637
1638 !inserting real data into complex array of half length
1639 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 IF (nfft == lot3) THEN
1645 CALL fft_1d(fft_plan_bw3, zw1, zw2, 1.0_dp, stat)
1646 ELSE
1647 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 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 IF (nproc > 1) THEN
1663 !communication scheduling
1664 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 DO j3 = 1, nd3/nproc
1670 !this condition ensures that we manage only the interesting part for the FFT
1671 IF (iproc*(nd3/nproc) + j3 <= n3/2 + 1) THEN
1672 jp2stb = 1
1673 j2stb = 1
1674 jp2stf = 1
1675 j2stf = 1
1676
1677 ! transform along x axis
1678
1679 DO j = 1, n2/2, lot1
1680 ma = j
1681 mb = min(j + (lot1 - 1), n2/2)
1682 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 IF (nproc == 1) THEN
1687 CALL mpiswitch_upcorn(j3, nfft, jp2stb, j2stb, lot1, n1, md2, nd3, nproc, zmpi2, zw1)
1688 ELSE
1689 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 IF (nfft == lot1) THEN
1697 CALL fft_1d(fft_plan_bw1, zw1, zt(j:, 1), 1.0_dp, stat)
1698 ELSE
1699 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 DO j = 1, n1, lot2
1707 ma = j
1708 mb = min(j + (lot2 - 1), n1)
1709 nfft = mb - ma + 1
1710
1711 !reverse ordering
1712 !input: I2,i1,j3,(jp3)
1713 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 IF (nfft == lot2) THEN
1720 CALL fft_1d(fft_plan_bw2, zw1, zw2, 1.0_dp, stat)
1721 ELSE
1722 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 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 IF (nfft == lot2) THEN
1735 CALL fft_1d(fft_plan_fw2, zw2, zw1, 1.0_dp, stat)
1736 ELSE
1737 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 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 DO j = 1, n2/2, lot1
1749 ma = j
1750 mb = min(j + (lot1 - 1), n2/2)
1751 nfft = mb - ma + 1
1752
1753 !performing FFT
1754 ! This is equivalent to
1755 IF (nfft == lot1) THEN
1756 CALL fft_1d(fft_plan_fw1, zt(j:, 1), zw2, 1.0_dp, stat)
1757 ELSE
1758 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 IF (nproc == 1) THEN
1765 CALL unmpiswitch_downcorn(j3, nfft, jp2stf, j2stf, lot1, n1, md2, nd3, nproc, zw2, zmpi2)
1766 ELSE
1767 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 IF (nproc > 1) THEN
1777 !communication scheduling
1778 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 DO j2 = 1, md2/nproc
1785 !this condition ensures that we manage only the interesting part for the FFT
1786 IF (iproc*(md2/nproc) + j2 <= n2/2) THEN
1787 DO i1 = 1, (n1/2), lot3
1788 ma = i1
1789 mb = min(i1 + (lot3 - 1), (n1/2))
1790 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 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 IF (nfft == lot3) THEN
1801 CALL fft_1d(fft_plan_fw3, zw1, zw2, 1.0_dp, stat)
1802 ELSE
1803 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 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 IF (n2/2 >= lot1) CALL fft_destroy_plan(fft_plan_bw1)
1814 IF (final_chunk_size1 > 0) CALL fft_destroy_plan(fft_plan_bw1_last)
1815 IF (n2/2 >= lot1) CALL fft_destroy_plan(fft_plan_fw1)
1816 IF (final_chunk_size1 > 0) CALL fft_destroy_plan(fft_plan_fw1_last)
1817
1818 IF (n1 >= lot2) CALL fft_destroy_plan(fft_plan_bw2)
1819 IF (final_chunk_size2 > 0) CALL fft_destroy_plan(fft_plan_bw2_last)
1820 IF (n1 >= lot2) CALL fft_destroy_plan(fft_plan_fw2)
1821 IF (final_chunk_size2 > 0) CALL fft_destroy_plan(fft_plan_fw2_last)
1822
1823 IF (n1/2 >= lot3) CALL fft_destroy_plan(fft_plan_fw3)
1824 IF (final_chunk_size3 > 0) CALL fft_destroy_plan(fft_plan_fw3_last)
1825 IF (n1/2 >= lot3) CALL fft_destroy_plan(fft_plan_bw3)
1826 IF (final_chunk_size3 > 0) CALL fft_destroy_plan(fft_plan_bw3_last)
1827
1828 !De-allocations
1829 DEALLOCATE (zmpi2)
1830 CALL fft_dealloc(zw1)
1831 CALL fft_dealloc(zw2)
1832 CALL fft_dealloc(zt)
1833 DEALLOCATE (cosinarr)
1834 IF (nproc > 1) DEALLOCATE (zmpi1)
1835
1836 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 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 DO j = 1, nfft
1860 DO i = n2/2 + 1, n2
1861 zw(j, i) = zt(i - n2/2, j)
1862 END DO
1863 END DO
1864 ! High frequencies
1865 DO i = 1, n2/2
1866 DO j = 1, nfft
1867 zw(j, i) = cmplx(0.0_dp, 0.0_dp, dp)
1868 END DO
1869 END DO
1870 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 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 mfft = 0
1899 main: DO jp2 = jp2stb, nproc
1900 DO j2 = j2stb, md2/nproc
1901 mfft = mfft + 1
1902 IF (mfft > nfft) THEN
1903 jp2stb = jp2
1904 j2stb = j2
1905 EXIT main
1906 END IF
1907 DO i1 = 1, n1/2
1908 zw(mfft, i1) = cmplx(0.0_dp, 0.0_dp, dp)
1909 END DO
1910 DO i1 = n1/2 + 1, n1
1911 zw(mfft, i1) = zmpi1(i1 - n1/2, j2, j3, jp2)
1912 END DO
1913 END DO
1914 j2stb = 1
1915 END DO main
1916 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 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 DO j = 1, nfft
1940 DO i = 1, n2/2
1941 zt(i, j) = zw(j, i)
1942 END DO
1943 END DO
1944 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 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 mfft = 0
1974 main: DO jp2 = jp2stf, nproc
1975 DO j2 = j2stf, md2/nproc
1976 mfft = mfft + 1
1977 IF (mfft > nfft) THEN
1978 jp2stf = jp2
1979 j2stf = j2
1980 EXIT main
1981 END IF
1982 DO i1 = 1, n1/2
1983 zmpi1(i1, j2, j3, jp2) = zw(mfft, i1)
1984 END DO
1985 END DO
1986 j2stf = 1
1987 END DO main
1988 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 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 ehartreetmp = 0._dp
2027 DO i3 = 1, n3/4
2028 DO i1 = 1, nfft
2029 pot1 = scal*real(zw(i1, i3), dp)
2030 ehartreetmp = ehartreetmp + pot1*zf(i1, 2*i3 - 1)
2031 zf(i1, 2*i3 - 1) = pot1
2032 pot1 = scal*aimag(zw(i1, i3))
2033 ehartreetmp = ehartreetmp + pot1*zf(i1, 2*i3)
2034 zf(i1, 2*i3) = pot1
2035 END DO
2036 END DO
2037 END SUBROUTINE f_unfill_downcorn
2038
2039END MODULE ps_wavelet_base
main
int main(int argc, char *argv[])
Stand-alone miniapp for smoke-testing and benchmarking dbm_multiply.
Definition dbm_miniapp.c:248
for
for(int lxp=0;lxp<=lp;lxp++)
Definition grid_cpu_collint.h:66
fft_lib
Definition fft_lib.F:7
fft_lib::fft_destroy_plan
subroutine, public fft_destroy_plan(plan)
...
Definition fft_lib.F:243
fft_lib::fft_1d
subroutine, public fft_1d(plan, zin, zout, scale, stat)
...
Definition fft_lib.F:263
fft_lib::fft_create_plan_1d
subroutine, public fft_create_plan_1d(plan, fsign, trans_in, trans_out, ldx_in, ldx_out, n, m, zin, zout)
...
Definition fft_lib.F:200
fft_plan
Type to store data about a (1D or 3D) FFT, including FFTW plan.
Definition fft_plan.F:18
fft_tools
Definition fft_tools.F:29
fft_tools::bwfft
integer, parameter, public bwfft
Definition fft_tools.F:146
fft_tools::fwfft
integer, parameter, public fwfft
Definition fft_tools.F:146
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::sp
integer, parameter, public sp
Definition kinds.F:33
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
ps_wavelet_base
Creates the wavelet kernel for the wavelet based poisson solver.
Definition ps_wavelet_base.F:12
ps_wavelet_base::s_poissonsolver
subroutine, public s_poissonsolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, pot, zf, scal, mpi_group)
!HERE POT MUST BE THE KERNEL (BEWARE THE HALF DIMENSION) ****h* BigDFT/S_PoissonSolver (Based on suit...
Definition ps_wavelet_base.F:813
ps_wavelet_base::scramble_unpack
subroutine, public scramble_unpack(i1, j2, lot, nfft, n1, n3, md2, nproc, nd3, zw, zmpi2, cosinarr)
(Based on suitable modifications of S.Goedecker routines) Assign the correct planes to the work array...
Definition ps_wavelet_base.F:1378
ps_wavelet_base::f_poissonsolver
subroutine, public f_poissonsolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, pot, zf, scal, mpi_group)
(Based on suitable modifications of S.Goedecker routines) Applies the local FFT space Kernel to the d...
Definition ps_wavelet_base.F:1519
ps_wavelet_base::p_poissonsolver
subroutine, public p_poissonsolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, zf, scal, hx, hy, hz, mpi_group)
...
Definition ps_wavelet_base.F:58
fft_plan::fft_plan_type
Definition fft_plan.F:27
message_passing::mp_comm_type
Definition message_passing.F:178