(git:08825b9)
Loading...
Searching...
No Matches
ps_wavelet_kernel.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_kernel
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 fft_alloc,&
19 fft_dealloc
20 USE kinds, ONLY: dp
21 USE mathconstants, ONLY: pi
22 USE message_passing, ONLY: mp_comm_type
23 USE ps_wavelet_base, ONLY: scramble_unpack
24 USE ps_wavelet_scaling_function, ONLY: scaling_function,&
25 scf_recursion
26 USE ps_wavelet_util, ONLY: f_fft_dimensions,&
27 s_fft_dimensions
28#include "../base/base_uses.f90"
29
30 IMPLICIT NONE
31
32 PRIVATE
33
34 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ps_wavelet_kernel'
35
36! *** Public data types ***
37
38 PUBLIC :: createkernel
39
40CONTAINS
41
42! **************************************************************************************************
43!> \brief Allocate a pointer which corresponds to the zero-padded FFT slice needed for
44!> calculating the convolution with the kernel expressed in the interpolating scaling
45!> function basis. The kernel pointer is unallocated on input, allocated on output.
46!> \param geocode Indicates the boundary conditions (BC) of the problem:
47!> 'F' free BC, isolated systems.
48!> The program calculates the solution as if the given density is
49!> "alone" in R^3 space.
50!> 'S' surface BC, isolated in y direction, periodic in xz plane
51!> The given density is supposed to be periodic in the xz plane,
52!> so the dimensions in these direction mus be compatible with the FFT
53!> Beware of the fact that the isolated direction is y!
54!> 'P' periodic BC.
55!> The density is supposed to be periodic in all the three directions,
56!> then all the dimensions must be compatible with the FFT.
57!> No need for setting up the kernel.
58!> \param n01 dimensions of the real space grid to be hit with the Poisson Solver
59!> \param n02 dimensions of the real space grid to be hit with the Poisson Solver
60!> \param n03 dimensions of the real space grid to be hit with the Poisson Solver
61!> \param hx grid spacings. For the isolated BC case for the moment they are supposed to
62!> be equal in the three directions
63!> \param hy grid spacings. For the isolated BC case for the moment they are supposed to
64!> be equal in the three directions
65!> \param hz grid spacings. For the isolated BC case for the moment they are supposed to
66!> be equal in the three directions
67!> \param itype_scf order of the interpolating scaling functions used in the decomposition
68!> \param iproc ,nproc: number of process, number of processes
69!> \param nproc ...
70!> \param kernel pointer for the kernel FFT. Unallocated on input, allocated on output.
71!> Its dimensions are equivalent to the region of the FFT space for which the
72!> kernel is injective. This will divide by two each direction,
73!> since the kernel for the zero-padded convolution is real and symmetric.
74!> \param mpi_group ...
75!> \date February 2007
76!> \author Luigi Genovese
77!> \note Due to the fact that the kernel dimensions are unknown before the calling, the kernel
78!> must be declared as pointer in input of this routine.
79!> To avoid that, one can properly define the kernel dimensions by adding
80!> the nd1,nd2,nd3 arguments to the PS_dim4allocation routine, then eliminating the pointer
81!> declaration.
82! **************************************************************************************************
83 SUBROUTINE createkernel(geocode, n01, n02, n03, hx, hy, hz, itype_scf, iproc, nproc, kernel, mpi_group)
84
85 CHARACTER(len=1), INTENT(in) :: geocode
86 INTEGER, INTENT(in) :: n01, n02, n03
87 REAL(kind=dp), INTENT(in) :: hx, hy, hz
88 INTEGER, INTENT(in) :: itype_scf, iproc, nproc
89 REAL(kind=dp), POINTER :: kernel(:)
90
91 CLASS(mp_comm_type), INTENT(in) :: mpi_group
92
93 INTEGER :: m1, m2, m3, md1, md2, md3, n1, n2, n3, &
94 nd1, nd2, nd3, nlimd, nlimk
95 REAL(kind=dp) :: hgrid
96
97 hgrid = max(hx, hy, hz)
98
99 IF (geocode == 'P') THEN
100
101 CALL f_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, &
102 md1, md2, md3, nd1, nd2, nd3, nproc)
103
104 ALLOCATE (kernel(1))
105 nlimd = n2
106 nlimk = 0
107
108 ELSE IF (geocode == 'S') THEN
109
110 CALL s_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, &
111 md1, md2, md3, nd1, nd2, nd3, nproc)
112
113 ALLOCATE (kernel(nd1*nd2*nd3/nproc))
114
115 !the kernel must be built and scattered to all the processes
116
117 CALL surfaces_kernel(n1, n2, n3, m3, nd1, nd2, nd3, hx, hz, hy, &
118 itype_scf, kernel, iproc, nproc, mpi_group)
119 !last plane calculated for the density and the kernel
120
121 nlimd = n2
122 nlimk = n3/2 + 1
123 ELSE IF (geocode == 'F') THEN
124
125 !Build the Kernel
126
127 CALL f_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, &
128 md1, md2, md3, nd1, nd2, nd3, nproc)
129 ALLOCATE (kernel(nd1*nd2*nd3/nproc))
130
131 !the kernel must be built and scattered to all the processes
132 CALL free_kernel(n01, n02, n03, n1, n2, n3, nd1, nd2, nd3, hgrid, &
133 itype_scf, iproc, nproc, kernel, mpi_group)
134
135 !last plane calculated for the density and the kernel
136 nlimd = n2/2
137 nlimk = n3/2 + 1
138
139 ELSE
140
141 cpabort("No wavelet based poisson solver for given geometry")
142
143 END IF
144 END SUBROUTINE createkernel
145
146! **************************************************************************************************
147!> \brief Build the kernel of the Poisson operator with
148!> surfaces Boundary conditions
149!> in an interpolating scaling functions basis.
150!> Beware of the fact that the nonperiodic direction is y!
151!> \param n1 Dimensions for the FFT
152!> \param n2 Dimensions for the FFT
153!> \param n3 Dimensions for the FFT
154!> \param m3 Actual dimension in non-periodic direction
155!> \param nker1 Dimensions of the kernel (nker3=n3/2+1) nker(1,2)=n(1,2)/2+1
156!> \param nker2 Dimensions of the kernel (nker3=n3/2+1) nker(1,2)=n(1,2)/2+1
157!> \param nker3 Dimensions of the kernel (nker3=n3/2+1) nker(1,2)=n(1,2)/2+1
158!> \param h1 Mesh steps in the three dimensions
159!> \param h2 Mesh steps in the three dimensions
160!> \param h3 Mesh steps in the three dimensions
161!> \param itype_scf Order of the scaling function
162!> \param karray output array
163!> \param iproc Number of process
164!> \param nproc number of processes
165!> \param mpi_group ...
166!> \date October 2006
167!> \author L. Genovese
168! **************************************************************************************************
169 SUBROUTINE surfaces_kernel(n1, n2, n3, m3, nker1, nker2, nker3, h1, h2, h3, &
170 itype_scf, karray, iproc, nproc, mpi_group)
171
172 INTEGER, INTENT(in) :: n1, n2, n3, m3, nker1, nker2, nker3
173 REAL(kind=dp), INTENT(in) :: h1, h2, h3
174 INTEGER, INTENT(in) :: itype_scf, nproc, iproc
175 REAL(kind=dp), &
176 DIMENSION(nker1, nker2, nker3/nproc), &
177 INTENT(out) :: karray
178 TYPE(mp_comm_type), INTENT(in) :: mpi_group
179
180 INTEGER, PARAMETER :: n_points = 2**6, ncache_optimal = 8*1024
181
182 COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:), &
183 POINTER :: halfft_cache1, halfft_cache2
184 INTEGER :: final_chunk_size, i1, i2, i3, iend, imu, ind1, ind2, ipolyord, ireim, istart, j2, &
185 j2nd, j2st, jnd1, jp2, jreim, n_cell, n_range, n_scf, nact2, ncache, nfft, num_of_mus, &
186 shift, stat
187 REAL(kind=dp) :: a, b, c, cp, d, diff, dx, fei, fer, foi, &
188 for, fr, mu1, pion, ponx, pony, sp, &
189 value, x
190 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: kernel_scf, x_scf, y_scf
191 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cossinarr
192 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: kernel
193 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kernel_mpi
194 REAL(kind=dp), DIMENSION(9, 8) :: cpol
195 TYPE(fft_plan_type) :: fft_plan_bw, fft_plan_bw_last
196
197!Better if higher (1024 points are enough 10^{-14}: 2*itype_scf*n_points)
198!FFT arrays
199!coefficients for the polynomial interpolation
200!assign the values of the coefficients
201
202 karray = 0.0_dp
203 cpol(:, :) = 1._dp
204 cpol(1, 2) = .25_dp
205 cpol(1, 3) = 1._dp/3._dp
206 cpol(1, 4) = 7._dp/12._dp
207 cpol(2, 4) = 8._dp/3._dp
208 cpol(1, 5) = 19._dp/50._dp
209 cpol(2, 5) = 3._dp/2._dp
210 cpol(1, 6) = 41._dp/272._dp
211 cpol(2, 6) = 27._dp/34._dp
212 cpol(3, 6) = 27._dp/272._dp
213 cpol(1, 7) = 751._dp/2989._dp
214 cpol(2, 7) = 73._dp/61._dp
215 cpol(3, 7) = 27._dp/61._dp
216 cpol(1, 8) = -989._dp/4540._dp
217 cpol(2, 8) = -1472._dp/1135._dp
218 cpol(3, 8) = 232._dp/1135._dp
219 cpol(4, 8) = -2624._dp/1135._dp
220
221 !renormalize values
222 cpol(1, 1) = .5_dp*cpol(1, 1)
223 cpol(1:2, 2) = 2._dp/3._dp*cpol(1:2, 2)
224 cpol(1:2, 3) = 3._dp/8._dp*cpol(1:2, 3)
225 cpol(1:3, 4) = 2._dp/15._dp*cpol(1:3, 4)
226 cpol(1:3, 5) = 25._dp/144._dp*cpol(1:3, 5)
227 cpol(1:4, 6) = 34._dp/105._dp*cpol(1:4, 6)
228 cpol(1:4, 7) = 2989._dp/17280._dp*cpol(1:4, 7)
229 cpol(1:5, 8) = -454._dp/2835._dp*cpol(1:5, 8)
230
231 !assign the complete values
232 cpol(2, 1) = cpol(1, 1)
233 cpol(3, 2) = cpol(1, 2)
234 cpol(3, 3) = cpol(2, 3)
235 cpol(4, 3) = cpol(1, 3)
236 cpol(4, 4) = cpol(2, 4)
237 cpol(5, 4) = cpol(1, 4)
238 cpol(4, 5) = cpol(3, 5)
239 cpol(5, 5) = cpol(2, 5)
240 cpol(6, 5) = cpol(1, 5)
241 cpol(5, 6) = cpol(3, 6)
242 cpol(6, 6) = cpol(2, 6)
243 cpol(7, 6) = cpol(1, 6)
244 cpol(5, 7) = cpol(4, 7)
245 cpol(6, 7) = cpol(3, 7)
246 cpol(7, 7) = cpol(2, 7)
247 cpol(8, 7) = cpol(1, 7)
248 cpol(6, 8) = cpol(4, 8)
249 cpol(7, 8) = cpol(3, 8)
250 cpol(8, 8) = cpol(2, 8)
251 cpol(9, 8) = cpol(1, 8)
252
253 !Number of integration points : 2*itype_scf*n_points
254 n_scf = 2*itype_scf*n_points
255 !Allocations
256 ALLOCATE (x_scf(0:n_scf))
257 ALLOCATE (y_scf(0:n_scf))
258
259 !Build the scaling function
260 CALL scaling_function(itype_scf, n_scf, n_range, x_scf, y_scf)
261 !Step grid for the integration
262 dx = real(n_range, kind=dp)/real(n_scf, kind=dp)
263 !Extend the range (no more calculations because fill in by 0._dp)
264 n_cell = m3
265 n_range = max(n_cell, n_range)
266
267 !Allocations
268 ncache = ncache_optimal
269 !the HalFFT must be performed only in the third dimension,
270 !and nker3=n3/2+1, hence
271 IF (ncache <= (nker3 - 1)*4) ncache = nker3 - 1*4
272
273 !enlarge the second dimension of the kernel to be compatible with nproc
274 nact2 = nker2
275 enlarge_ydim: DO
276 IF (nproc*(nact2/nproc) /= nact2) THEN
277 nact2 = nact2 + 1
278 ELSE
279 EXIT enlarge_ydim
280 END IF
281 END DO enlarge_ydim
282
283 !array for the MPI procedure
284 ALLOCATE (kernel(nker1, nact2/nproc, nker3))
285 ALLOCATE (kernel_mpi(nker1, nact2/nproc, nker3/nproc, nproc))
286 ALLOCATE (kernel_scf(n_range))
287 CALL fft_alloc(halfft_cache1, [max(ncache/4, n3/2)])
288 halfft_cache1 = cmplx(0.0_dp, 0.0_dp, dp)
289 CALL fft_alloc(halfft_cache2, [max(ncache/4, n3/2)])
290 halfft_cache2 = cmplx(0.0_dp, 0.0_dp, dp)
291 ALLOCATE (cossinarr(2, n3/2 - 1))
292
293 !build the phases for the HalFFT reconstruction
294 pion = 2._dp*pi/real(n3, kind=dp)
295 DO i3 = 2, n3/2
296 x = real(i3 - 1, kind=dp)*pion
297 cossinarr(1, i3 - 1) = cos(x)
298 cossinarr(2, i3 - 1) = -sin(x)
299 END DO
300
301 ! satisfy valgrind, init arrays to large value, even if the offending bit is (likely?) padding
302 kernel = huge(0._dp)
303 kernel_mpi = huge(0._dp)
304
305 !calculate the limits of the FFT calculations
306 !that can be performed in a row remaining inside the cache
307 num_of_mus = max(1, ncache/(2*n3))
308
309 diff = 0._dp
310 !order of the polynomial to be used for integration (must be a power of two)
311 ipolyord = 8 !this part should be incorporated inside the numerical integration
312 !here we have to choice the piece of the x-y grid to cover
313
314 !let us now calculate the fraction of mu that will be considered
315 j2st = iproc*(nact2/nproc)
316 j2nd = min((iproc + 1)*(nact2/nproc), n2/2 + 1)
317
318 ! The size of the last chunk
319 final_chunk_size = mod((n1/2 + 1)*(j2nd - j2st), num_of_mus)
320
321 ! Prevent OOB-access if n1 < num_of_mus
322 IF ((n1/2 + 1)*(j2nd - j2st) >= num_of_mus) THEN
323 CALL fft_create_plan_1d(fft_plan_bw, bwfft, .true., .true., num_of_mus, num_of_mus, &
324 n3/2, num_of_mus, halfft_cache1, halfft_cache2)
325 END IF
326 IF (final_chunk_size > 0) THEN
327 CALL fft_create_plan_1d(fft_plan_bw_last, bwfft, .true., .true., num_of_mus, num_of_mus, &
328 n3/2, final_chunk_size, halfft_cache1, halfft_cache2)
329 END IF
330
331 DO ind2 = (n1/2 + 1)*j2st + 1, (n1/2 + 1)*j2nd, num_of_mus
332 istart = ind2
333 iend = min(ind2 + (num_of_mus - 1), (n1/2 + 1)*j2nd)
334 nfft = iend - istart + 1
335 shift = 0
336
337 !initialization of the interesting part of the cache array
338 ! THis could be replaced using a R2C-FFT (not yet available in CP2K)
339 halfft_cache1 = cmplx(0.0_dp, 0.0_dp, dp)
340 halfft_cache1 = cmplx(0.0_dp, 0.0_dp, dp)
341
342 IF (istart == 1) THEN
343 !i2=1
344 shift = 1
345
346 CALL calculates_green_opt_muzero(n_range, n_scf, ipolyord, x_scf, y_scf, &
347 cpol(1, ipolyord), dx, kernel_scf)
348
349 !copy of the first zero value
350 halfft_cache1(1) = cmplx(0._dp, aimag(halfft_cache1(1)), dp)
351
352 DO i3 = 1, m3
353
354 value = 0.5_dp*h3*kernel_scf(i3)
355 !index in where to copy the value of the kernel
356 CALL indices(ireim, num_of_mus, n3/2 + i3, 1, ind1)
357 !index in where to copy the symmetric value
358 CALL indices(jreim, num_of_mus, n3/2 + 2 - i3, 1, jnd1)
359 IF (ireim == 1) THEN
360 halfft_cache1(ind1) = cmplx(value, aimag(halfft_cache1(ind1)), dp)
361 ELSE
362 halfft_cache1(ind1) = cmplx(real(halfft_cache1(ind1), dp), value, dp)
363 END IF
364 IF (jreim == 1) THEN
365 halfft_cache1(jnd1) = cmplx(value, aimag(halfft_cache1(jnd1)), dp)
366 ELSE
367 halfft_cache1(jnd1) = cmplx(real(halfft_cache1(jnd1), dp), value, dp)
368 END IF
369
370 END DO
371
372 END IF
373
374 loopimpulses: DO imu = istart + shift, iend
375
376 !here there is the value of mu associated to hgrid
377 !note that we have multiplicated mu for hgrid to be comparable
378 !with mu0ref
379
380 !calculate the proper value of mu taking into account the periodic dimensions
381 !corresponding value of i1 and i2
382 i1 = mod(imu, n1/2 + 1)
383 IF (i1 == 0) i1 = n1/2 + 1
384 i2 = (imu - i1)/(n1/2 + 1) + 1
385 ponx = real(i1 - 1, kind=dp)/real(n1, kind=dp)
386 pony = real(i2 - 1, kind=dp)/real(n2, kind=dp)
387
388 mu1 = 2._dp*pi*sqrt((ponx/h1)**2 + (pony/h2)**2)*h3
389
390 CALL calculates_green_opt(n_range, n_scf, itype_scf, ipolyord, x_scf, y_scf, &
391 cpol(1, ipolyord), mu1, dx, kernel_scf)
392
393 !readjust the coefficient and define the final kernel
394
395 !copy of the first zero value
396 halfft_cache1(imu - istart + 1) = cmplx(0._dp, aimag(halfft_cache1(imu - istart + 1)), dp)
397 DO i3 = 1, m3
398 value = -0.5_dp*h3/mu1*kernel_scf(i3)
399 !index in where to copy the value of the kernel
400 CALL indices(ireim, num_of_mus, n3/2 + i3, imu - istart + 1, ind1)
401 !index in where to copy the symmetric value
402 CALL indices(jreim, num_of_mus, n3/2 + 2 - i3, imu - istart + 1, jnd1)
403 IF (ireim == 1) THEN
404 halfft_cache1(ind1) = cmplx(value, aimag(halfft_cache1(ind1)), dp)
405 ELSE
406 halfft_cache1(ind1) = cmplx(real(halfft_cache1(ind1), dp), value, dp)
407 END IF
408 IF (jreim == 1) THEN
409 halfft_cache1(jnd1) = cmplx(value, aimag(halfft_cache1(jnd1)), dp)
410 ELSE
411 halfft_cache1(jnd1) = cmplx(real(halfft_cache1(jnd1), dp), value, dp)
412 END IF
413 END DO
414
415 END DO loopimpulses
416
417 !now perform the FFT of the array in cache
418 ! This is equivalent to
419 IF (nfft == num_of_mus) THEN
420 CALL fft_1d(fft_plan_bw, halfft_cache1, halfft_cache2, 1.0_dp, stat)
421 ELSE
422 CALL fft_1d(fft_plan_bw_last, halfft_cache1, halfft_cache2, 1.0_dp, stat)
423 END IF
424
425 !assign the values of the FFT array
426 !and compare with the good results
427 DO imu = istart, iend
428
429 !corresponding value of i1 and i2
430 i1 = mod(imu, n1/2 + 1)
431 IF (i1 == 0) i1 = n1/2 + 1
432 i2 = (imu - i1)/(n1/2 + 1) + 1
433
434 j2 = i2 - j2st
435
436 a = real(halfft_cache2(imu - istart + 1), dp)
437 b = aimag(halfft_cache2(imu - istart + 1))
438 kernel(i1, j2, 1) = a + b
439 kernel(i1, j2, n3/2 + 1) = a - b
440
441 DO i3 = 2, n3/2
442 ind1 = imu - istart + 1 + num_of_mus*(i3 - 1)
443 jnd1 = imu - istart + 1 + num_of_mus*(n3/2 + 2 - i3 - 1)
444 cp = cossinarr(1, i3 - 1)
445 sp = cossinarr(2, i3 - 1)
446 a = real(halfft_cache2(ind1), dp)
447 b = aimag(halfft_cache2(ind1))
448 c = real(halfft_cache2(jnd1), dp)
449 d = aimag(halfft_cache2(jnd1))
450 fer = .5_dp*(a + c)
451 fei = .5_dp*(b - d)
452 for = .5_dp*(a - c)
453 foi = .5_dp*(b + d)
454 fr = fer + cp*foi - sp*for
455 kernel(i1, j2, i3) = fr
456 END DO
457 END DO
458
459 END DO
460
461 IF ((n1/2 + 1)*(j2nd - j2st) >= num_of_mus) CALL fft_destroy_plan(fft_plan_bw)
462 IF (final_chunk_size > 0) CALL fft_destroy_plan(fft_plan_bw_last)
463
464 !give to each processor a slice of the third dimension
465 IF (nproc > 1) THEN
466 CALL mpi_group%alltoall(kernel, &!nker1*(nact2/nproc)*(nker3/nproc), &
467 kernel_mpi, nker1*(nact2/nproc)*(nker3/nproc))
468
469 DO jp2 = 1, nproc
470 DO i3 = 1, nker3/nproc
471 DO i2 = 1, nact2/nproc
472 j2 = i2 + (jp2 - 1)*(nact2/nproc)
473 IF (j2 <= nker2) THEN
474 DO i1 = 1, nker1
475 karray(i1, j2, i3) = &
476 kernel_mpi(i1, i2, i3, jp2)
477 END DO
478 END IF
479 END DO
480 END DO
481 END DO
482
483 ELSE
484 karray(1:nker1, 1:nker2, 1:nker3) = kernel(1:nker1, 1:nker2, 1:nker3)
485 END IF
486
487 !De-allocations
488 DEALLOCATE (kernel)
489 DEALLOCATE (kernel_mpi)
490 CALL fft_dealloc(halfft_cache1)
491 CALL fft_dealloc(halfft_cache2)
492 DEALLOCATE (kernel_scf)
493 DEALLOCATE (x_scf)
494 DEALLOCATE (y_scf)
495
496 END SUBROUTINE surfaces_kernel
497
498! **************************************************************************************************
499!> \brief ...
500!> \param n ...
501!> \param n_scf ...
502!> \param itype_scf ...
503!> \param intorder ...
504!> \param xval ...
505!> \param yval ...
506!> \param c ...
507!> \param mu ...
508!> \param hres ...
509!> \param g_mu ...
510! **************************************************************************************************
511 SUBROUTINE calculates_green_opt(n, n_scf, itype_scf, intorder, xval, yval, c, mu, hres, g_mu)
512 INTEGER, INTENT(in) :: n, n_scf, itype_scf, intorder
513 REAL(kind=dp), DIMENSION(0:n_scf), INTENT(in) :: xval, yval
514 REAL(kind=dp), DIMENSION(intorder+1), INTENT(in) :: c
515 REAL(kind=dp), INTENT(in) :: mu, hres
516 REAL(kind=dp), DIMENSION(n), INTENT(out) :: g_mu
517
518 REAL(kind=dp), PARAMETER :: mu_max = 0.2_dp
519
520 INTEGER :: i, iend, ikern, ivalue, izero, n_iter, &
521 nrec
522 REAL(kind=dp) :: f, filter, fl, fr, gleft, gltmp, gright, &
523 grtmp, mu0, ratio, x, x0, x1
524 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: green, green1
525
526 g_mu = 0.0_dp
527 !We calculate the number of iterations to go from mu0 to mu0_ref
528 IF (mu <= mu_max) THEN
529 n_iter = 0
530 mu0 = mu
531 ELSE
532 n_iter = 1
533 loop_iter: DO
534 ratio = real(2**n_iter, kind=dp)
535 mu0 = mu/ratio
536 IF (mu0 <= mu_max) THEN
537 EXIT loop_iter
538 END IF
539 n_iter = n_iter + 1
540 END DO loop_iter
541 END IF
542
543 !dimension needed for the correct calculation of the recursion
544 nrec = 2**n_iter*n
545
546 ALLOCATE (green(-nrec:nrec))
547
548 !initialization of the branching value
549 ikern = 0
550 izero = 0
551 initialization: DO
552 IF (xval(izero) >= real(ikern, kind=dp) .OR. izero == n_scf) EXIT initialization
553 izero = izero + 1
554 END DO initialization
555 green = 0._dp
556 !now perform the interpolation in right direction
557 ivalue = izero
558 gright = 0._dp
559 loop_right: DO
560 IF (ivalue >= n_scf - intorder - 1) EXIT loop_right
561 DO i = 1, intorder + 1
562 x = xval(ivalue) - real(ikern, kind=dp)
563 f = yval(ivalue)*exp(-mu0*x)
564 filter = intorder*c(i)
565 gright = gright + filter*f
566 ivalue = ivalue + 1
567 END DO
568 ivalue = ivalue - 1
569 END DO loop_right
570 iend = n_scf - ivalue
571 DO i = 1, iend
572 x = xval(ivalue) - real(ikern, kind=dp)
573 f = yval(ivalue)*exp(-mu0*x)
574 filter = intorder*c(i)
575 gright = gright + filter*f
576 ivalue = ivalue + 1
577 END DO
578 gright = hres*gright
579
580 !the scaling function is symmetric, so the same for the other direction
581 gleft = gright
582
583 green(ikern) = gleft + gright
584
585 !now the loop until the last value
586 DO ikern = 1, nrec
587 gltmp = 0._dp
588 grtmp = 0._dp
589 ivalue = izero
590 x0 = xval(izero)
591 loop_integration: DO
592 IF (izero == n_scf) EXIT loop_integration
593 DO i = 1, intorder + 1
594 x = xval(ivalue)
595 fl = yval(ivalue)*exp(mu0*x)
596 fr = yval(ivalue)*exp(-mu0*x)
597 filter = intorder*c(i)
598 gltmp = gltmp + filter*fl
599 grtmp = grtmp + filter*fr
600 ivalue = ivalue + 1
601 IF (xval(izero) >= real(ikern, kind=dp) .OR. izero == n_scf) THEN
602 x1 = xval(izero)
603 EXIT loop_integration
604 END IF
605 izero = izero + 1
606 END DO
607 ivalue = ivalue - 1
608 izero = izero - 1
609 END DO loop_integration
610 gleft = exp(-mu0)*(gleft + hres*exp(-mu0*real(ikern - 1, kind=dp))*gltmp)
611 IF (izero == n_scf) THEN
612 gright = 0._dp
613 ELSE
614 gright = exp(mu0)*(gright - hres*exp(mu0*real(ikern - 1, kind=dp))*grtmp)
615 END IF
616 green(ikern) = gleft + gright
617 green(-ikern) = gleft + gright
618 IF (abs(green(ikern)) <= 1.e-20_dp) THEN
619 nrec = ikern
620 EXIT
621 END IF
622 END DO
623 !now we must calculate the recursion
624 ALLOCATE (green1(-nrec:nrec))
625
626 !Start the iteration to go from mu0 to mu
627 CALL scf_recursion(itype_scf, n_iter, nrec, green(-nrec), green1(-nrec))
628
629 DO i = 1, min(n, nrec)
630 g_mu(i) = green(i - 1)
631 END DO
632 DO i = min(n, nrec) + 1, n
633 g_mu(i) = 0._dp
634 END DO
635
636 DEALLOCATE (green, green1)
637
638 END SUBROUTINE calculates_green_opt
639
640! **************************************************************************************************
641!> \brief ...
642!> \param n ...
643!> \param n_scf ...
644!> \param intorder ...
645!> \param xval ...
646!> \param yval ...
647!> \param c ...
648!> \param hres ...
649!> \param green ...
650! **************************************************************************************************
651 PURE SUBROUTINE calculates_green_opt_muzero(n, n_scf, intorder, xval, yval, c, hres, green)
652 INTEGER, INTENT(in) :: n, n_scf, intorder
653 REAL(kind=dp), DIMENSION(0:n_scf), INTENT(in) :: xval, yval
654 REAL(kind=dp), DIMENSION(intorder+1), INTENT(in) :: c
655 REAL(kind=dp), INTENT(in) :: hres
656 REAL(kind=dp), DIMENSION(n), INTENT(out) :: green
657
658 INTEGER :: i, iend, ikern, ivalue, izero
659 REAL(kind=dp) :: c0, c1, filter, gl0, gl1, gr0, gr1, x, y
660
661!initialization of the branching value
662
663 ikern = 0
664 izero = 0
665 initialization: DO
666 IF (xval(izero) >= real(ikern, kind=dp) .OR. izero == n_scf) EXIT initialization
667 izero = izero + 1
668 END DO initialization
669 green = 0._dp
670 !first case, ikern=0
671 !now perform the interpolation in right direction
672 ivalue = izero
673 gr1 = 0._dp
674 loop_right: DO
675 IF (ivalue >= n_scf - intorder - 1) EXIT loop_right
676 DO i = 1, intorder + 1
677 x = xval(ivalue)
678 y = yval(ivalue)
679 filter = intorder*c(i)
680 gr1 = gr1 + filter*x*y
681 ivalue = ivalue + 1
682 END DO
683 ivalue = ivalue - 1
684 END DO loop_right
685 iend = n_scf - ivalue
686 DO i = 1, iend
687 x = xval(ivalue)
688 y = yval(ivalue)
689 filter = intorder*c(i)
690 gr1 = gr1 + filter*x*y
691 ivalue = ivalue + 1
692 END DO
693 gr1 = hres*gr1
694 !the scaling function is symmetric
695 gl1 = -gr1
696 gl0 = 0.5_dp
697 gr0 = 0.5_dp
698
699 green(1) = 2._dp*gr1
700
701 !now the loop until the last value
702 DO ikern = 1, n - 1
703 c0 = 0._dp
704 c1 = 0._dp
705 ivalue = izero
706 loop_integration: DO
707 IF (izero == n_scf) EXIT loop_integration
708 DO i = 1, intorder + 1
709 x = xval(ivalue)
710 y = yval(ivalue)
711 filter = intorder*c(i)
712 c0 = c0 + filter*y
713 c1 = c1 + filter*y*x
714 ivalue = ivalue + 1
715 IF (xval(izero) >= real(ikern, kind=dp) .OR. izero == n_scf) THEN
716 EXIT loop_integration
717 END IF
718 izero = izero + 1
719 END DO
720 ivalue = ivalue - 1
721 izero = izero - 1
722 END DO loop_integration
723 c0 = hres*c0
724 c1 = hres*c1
725
726 gl0 = gl0 + c0
727 gl1 = gl1 + c1
728 gr0 = gr0 - c0
729 gr1 = gr1 - c1
730 !general case
731 green(ikern + 1) = real(ikern, kind=dp)*(gl0 - gr0) + gr1 - gl1
732 END DO
733
734 END SUBROUTINE calculates_green_opt_muzero
735
736! **************************************************************************************************
737!> \brief ...
738!> \param var_realimag ...
739!> \param nelem ...
740!> \param intrn ...
741!> \param extrn ...
742!> \param index ...
743! **************************************************************************************************
744 PURE ELEMENTAL SUBROUTINE indices(var_realimag, nelem, intrn, extrn, index)
745
746 INTEGER, INTENT(out) :: var_realimag
747 INTEGER, INTENT(in) :: nelem, intrn, extrn
748 INTEGER, INTENT(out) :: index
749
750 INTEGER :: i
751
752!real or imaginary part
753
754 var_realimag = 2 - mod(intrn, 2)
755!actual index over half the length
756
757 i = (intrn + 1)/2
758 !complete index to be assigned
759 index = extrn + nelem*(i - 1)
760
761 END SUBROUTINE indices
762
763! **************************************************************************************************
764!> \brief Build the kernel of a gaussian function
765!> for interpolating scaling functions.
766!> Do the parallel HalFFT of the symmetrized function and stores into
767!> memory only 1/8 of the grid divided by the number of processes nproc
768!>
769!> Build the kernel (karray) of a gaussian function
770!> for interpolating scaling functions
771!> $$ K(j) = \sum_k \omega_k \int \int \phi(x) g_k(x'-x) \delta(x'- j) dx dx' $$
772!> \param n01 Mesh dimensions of the density
773!> \param n02 Mesh dimensions of the density
774!> \param n03 Mesh dimensions of the density
775!> \param nfft1 Dimensions of the FFT grid (HalFFT in the third direction)
776!> \param nfft2 Dimensions of the FFT grid (HalFFT in the third direction)
777!> \param nfft3 Dimensions of the FFT grid (HalFFT in the third direction)
778!> \param n1k Dimensions of the kernel FFT
779!> \param n2k Dimensions of the kernel FFT
780!> \param n3k Dimensions of the kernel FFT
781!> \param hgrid Mesh step
782!> \param itype_scf Order of the scaling function (8,14,16)
783!> \param iproc ...
784!> \param nproc ...
785!> \param karray ...
786!> \param mpi_group ...
787!> \date February 2006
788!> \author T. Deutsch, L. Genovese
789! **************************************************************************************************
790 SUBROUTINE free_kernel(n01, n02, n03, nfft1, nfft2, nfft3, n1k, n2k, n3k, &
791 hgrid, itype_scf, iproc, nproc, karray, mpi_group)
792
793 INTEGER, INTENT(in) :: n01, n02, n03, nfft1, nfft2, nfft3, n1k, &
794 n2k, n3k
795 REAL(kind=dp), INTENT(in) :: hgrid
796 INTEGER, INTENT(in) :: itype_scf, iproc, nproc
797 REAL(kind=dp), DIMENSION(n1k, n2k, n3k/nproc), &
798 INTENT(out) :: karray
799 TYPE(mp_comm_type), INTENT(in) :: mpi_group
800
801 INTEGER, PARAMETER :: n_gauss = 89, n_points = 2**6
802 REAL(kind=dp), PARAMETER :: p0_ref = 1._dp
803
804 INTEGER :: i, i01, i02, i03, i1, i2, i3, i_gauss, &
805 i_kern, iend, istart, istart1, n1h, &
806 n2h, n3h, n_cell, n_iter, n_range, &
807 n_scf, nker1, nker2, nker3
808 REAL(kind=dp) :: a1, a2, a3, a_range, absci, acc_gauss, &
809 dr_gauss, dx, factor, factor2, kern, &
810 p0_cell, p0gauss, pgauss, ur_gauss
811 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: kern_1_scf, kernel_scf, x_scf, y_scf
812 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: kp
813 REAL(kind=dp), DIMENSION(n_gauss) :: p_gauss, w_gauss
814
815!Do not touch !!!!
816!Better if higher (1024 points are enough 10^{-14}: 2*itype_scf*n_points)
817!Better p_gauss for calculation
818!(the support of the exponential should be inside [-n_range/2,n_range/2])
819!Number of integration points : 2*itype_scf*n_points
820
821 n_scf = 2*itype_scf*n_points
822 !Set karray
823 karray = 0.0_dp
824 !here we must set the dimensions for the fft part, starting from the nfft
825 !remember that actually nfft2 is associated to n03 and viceversa
826
827 !dimensions that define the center of symmetry
828 n1h = nfft1/2
829 n2h = nfft2/2
830 n3h = nfft3/2
831
832 !Auxiliary dimensions only for building the FFT part
833 nker1 = nfft1
834 nker2 = nfft2
835 nker3 = nfft3/2 + 1
836
837 !adjusting the last two dimensions to be multiples of nproc
838 DO
839 IF (modulo(nker2, nproc) == 0) EXIT
840 nker2 = nker2 + 1
841 END DO
842 DO
843 IF (modulo(nker3, nproc) == 0) EXIT
844 nker3 = nker3 + 1
845 END DO
846
847 !this will be the array of the kernel in the real space
848 ALLOCATE (kp(n1h + 1, n3h + 1, nker2/nproc))
849
850 !defining proper extremes for the calculation of the
851 !local part of the kernel
852
853 istart = iproc*nker2/nproc + 1
854 iend = min((iproc + 1)*nker2/nproc, n2h + n03)
855
856 istart1 = istart
857 IF (iproc == 0) istart1 = n2h - n03 + 2
858
859 !Allocations
860 ALLOCATE (x_scf(0:n_scf))
861 ALLOCATE (y_scf(0:n_scf))
862
863 !Build the scaling function
864 CALL scaling_function(itype_scf, n_scf, n_range, x_scf, y_scf)
865 !Step grid for the integration
866 dx = real(n_range, kind=dp)/real(n_scf, kind=dp)
867 !Extend the range (no more calculations because fill in by 0._dp)
868 n_cell = max(n01, n02, n03)
869 n_range = max(n_cell, n_range)
870
871 !Allocations
872 ALLOCATE (kernel_scf(-n_range:n_range))
873 ALLOCATE (kern_1_scf(-n_range:n_range))
874
875 !Lengthes of the box (use FFT dimension)
876 a1 = hgrid*real(n01, kind=dp)
877 a2 = hgrid*real(n02, kind=dp)
878 a3 = hgrid*real(n03, kind=dp)
879
880 x_scf(:) = hgrid*x_scf(:)
881 y_scf(:) = 1._dp/hgrid*y_scf(:)
882 dx = hgrid*dx
883 !To have a correct integration
884 p0_cell = p0_ref/(hgrid*hgrid)
885
886 !Initialization of the gaussian (Beylkin)
887 CALL gequad(p_gauss, w_gauss, ur_gauss, dr_gauss, acc_gauss)
888 !In order to have a range from a_range=sqrt(a1*a1+a2*a2+a3*a3)
889 !(biggest length in the cube)
890 !We divide the p_gauss by a_range**2 and a_gauss by a_range
891 a_range = sqrt(a1*a1 + a2*a2 + a3*a3)
892 factor = 1._dp/a_range
893 !factor2 = factor*factor
894 factor2 = 1._dp/(a1*a1 + a2*a2 + a3*a3)
895 DO i_gauss = 1, n_gauss
896 p_gauss(i_gauss) = factor2*p_gauss(i_gauss)
897 END DO
898 DO i_gauss = 1, n_gauss
899 w_gauss(i_gauss) = factor*w_gauss(i_gauss)
900 END DO
901
902 kp(:, :, :) = 0._dp
903 !Use in this order (better for accuracy).
904 loop_gauss: DO i_gauss = n_gauss, 1, -1
905 !Gaussian
906 pgauss = p_gauss(i_gauss)
907
908 !We calculate the number of iterations to go from pgauss to p0_ref
909 n_iter = nint((log(pgauss) - log(p0_cell))/log(4._dp))
910 IF (n_iter <= 0) THEN
911 n_iter = 0
912 p0gauss = pgauss
913 ELSE
914 p0gauss = pgauss/4._dp**n_iter
915 END IF
916
917 !Stupid integration
918 !Do the integration with the exponential centered in i_kern
919 kernel_scf(:) = 0._dp
920 DO i_kern = 0, n_range
921 kern = 0._dp
922 DO i = 0, n_scf
923 absci = x_scf(i) - real(i_kern, kind=dp)*hgrid
924 absci = absci*absci
925 kern = kern + y_scf(i)*exp(-p0gauss*absci)*dx
926 END DO
927 kernel_scf(i_kern) = kern
928 kernel_scf(-i_kern) = kern
929 IF (abs(kern) < 1.e-18_dp) THEN
930 !Too small not useful to calculate
931 EXIT
932 END IF
933 END DO
934
935 !Start the iteration to go from p0gauss to pgauss
936 CALL scf_recursion(itype_scf, n_iter, n_range, kernel_scf, kern_1_scf)
937
938 !Add to the kernel (only the local part)
939
940 DO i3 = istart1, iend
941 i03 = i3 - n2h - 1
942 ! Crash if index out of range
943 ! Without compiler bounds checking, the calculation might run successfully but
944 ! it is also possible that the Hartree energy will blow up
945 ! This seems to happen with large MPI processor counts if the size of the
946 ! RS grid is not directly compatible with the allowed FFT dimensions
947 IF (i03 < -n_range .OR. i03 > n_range) THEN
948 CALL cp_abort(__location__, "Index out of range in wavelet solver. "// &
949 "Try decreasing the number of MPI processors, or adjust the PW_CUTOFF or cell size "// &
950 "so that 2*(number of RS grid points) matches the allowed FFT dimensions "// &
951 "(see ps_wavelet_fft3d.F) exactly.")
952 END IF
953 DO i2 = 1, n02
954 i02 = i2 - 1
955 DO i1 = 1, n01
956 i01 = i1 - 1
957 kp(i1, i2, i3 - istart + 1) = kp(i1, i2, i3 - istart + 1) + w_gauss(i_gauss)* &
958 kernel_scf(i01)*kernel_scf(i02)*kernel_scf(i03)
959 END DO
960 END DO
961 END DO
962
963 END DO loop_gauss
964
965 !De-allocations
966 DEALLOCATE (kernel_scf)
967 DEALLOCATE (kern_1_scf)
968 DEALLOCATE (x_scf)
969 DEALLOCATE (y_scf)
970
971!!!!END KERNEL CONSTRUCTION
972
973!!$ if(iproc .eq. 0) print *,"Do a 3D PHalFFT for the kernel"
974
975 CALL kernelfft(nfft1, nfft2, nfft3, nker1, nker2, nker3, n1k, n2k, n3k, nproc, iproc, &
976 kp, karray, mpi_group)
977
978 !De-allocations
979 DEALLOCATE (kp)
980
981 END SUBROUTINE free_kernel
982
983! **************************************************************************************************
984!> \brief ...
985!> \param n1 ...
986!> \param n3 ...
987!> \param lot ...
988!> \param nfft ...
989!> \param i1 ...
990!> \param zf ...
991!> \param zw ...
992! **************************************************************************************************
993 PURE SUBROUTINE inserthalf(n1, n3, lot, nfft, i1, zf, zw)
994 INTEGER, INTENT(in) :: n1, n3, lot, nfft, i1
995 REAL(kind=dp), DIMENSION(n1/2+1, n3/2+1), &
996 INTENT(in) :: zf
997 COMPLEX(KIND=dp), DIMENSION(lot, n3/2), &
998 INTENT(out) :: zw
999
1000 INTEGER :: i01, i03i, i03r, i3, l1, l3
1001
1002 zw = 0.0_dp
1003 i3 = 0
1004 DO l3 = 1, n3, 2
1005 i3 = i3 + 1
1006 i03r = abs(l3 - n3/2 - 1) + 1
1007 i03i = abs(l3 - n3/2) + 1
1008 DO l1 = 1, nfft
1009 i01 = abs(l1 - 1 + i1 - n1/2 - 1) + 1
1010 zw(l1, i3) = cmplx(zf(i01, i03r), zf(i01, i03i), dp)
1011 END DO
1012 END DO
1013
1014 END SUBROUTINE inserthalf
1015
1016! **************************************************************************************************
1017!> \brief (Based on suitable modifications of S.Goedecker routines)
1018!> Calculates the FFT of the distributed kernel
1019!> \param n1 logical dimension of the transform.
1020!> \param n2 logical dimension of the transform.
1021!> \param n3 logical dimension of the transform.
1022!> \param nd1 Dimensions of work arrays
1023!> \param nd2 Dimensions of work arrays
1024!> \param nd3 Dimensions of work arrays
1025!> \param nk1 ...
1026!> \param nk2 ...
1027!> \param nk3 ...
1028!> \param nproc number of processors used as returned by MPI_COMM_SIZE
1029!> \param iproc [0:nproc-1] number of processor as returned by MPI_COMM_RANK
1030!> \param zf Real kernel (input)
1031!> zf(i1,i2,i3)
1032!> \param zr Distributed Kernel FFT
1033!> zr(2,i1,i2,i3)
1034!> \param mpi_group ...
1035!> \date February 2006
1036!> \par Restrictions
1037!> Copyright (C) Stefan Goedecker, Cornell University, Ithaca, USA, 1994
1038!> Copyright (C) Stefan Goedecker, MPI Stuttgart, Germany, 1999
1039!> Copyright (C) 2002 Stefan Goedecker, CEA Grenoble
1040!> This file is distributed under the terms of the
1041!> GNU General Public License, see http://www.gnu.org/copyleft/gpl.txt .
1042!> \author S. Goedecker, L. Genovese
1043! **************************************************************************************************
1044 SUBROUTINE kernelfft(n1, n2, n3, nd1, nd2, nd3, nk1, nk2, nk3, nproc, iproc, zf, zr, mpi_group)
1045
1046 INTEGER, INTENT(in) :: n1, n2, n3, nd1, nd2, nd3, nk1, nk2, &
1047 nk3, nproc, iproc
1048 REAL(kind=dp), &
1049 DIMENSION(n1/2+1, n3/2+1, nd2/nproc), &
1050 INTENT(in) :: zf
1051 REAL(kind=dp), DIMENSION(nk1, nk2, nk3/nproc), &
1052 INTENT(inout) :: zr
1053 TYPE(mp_comm_type), INTENT(in) :: mpi_group
1054
1055 INTEGER, PARAMETER :: ncache_optimal = 8*1024
1056
1057 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: zmpi2
1058 COMPLEX(KIND=dp), ALLOCATABLE, &
1059 DIMENSION(:, :, :, :) :: zmpi1
1060 COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:), &
1061 POINTER :: zw1, zw2
1062 COMPLEX(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1063 POINTER :: zt
1064 INTEGER :: final_chunk_size1, final_chunk_size2, final_chunk_size3, i1, i3, j, j2, j2st, j3, &
1065 jp2st, lot1, lot2, lot3, lzt, ma, mb, ncache, nfft, stat
1066 REAL(kind=dp) :: twopion
1067 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cosinarr
1068 TYPE(fft_plan_type) :: fft_plan_bw1, fft_plan_bw1_last, &
1069 fft_plan_bw2, fft_plan_bw2_last, &
1070 fft_plan_bw3, fft_plan_bw3_last
1071
1072!work arrays for transpositions
1073!work arrays for MPI
1074!cache work array
1075!FFT work arrays
1076!Body
1077!check input
1078
1079 cpassert(nd1 >= n1)
1080 cpassert(nd2 >= n2)
1081 cpassert(nd3 >= n3/2 + 1)
1082 cpassert(mod(nd3, nproc) == 0)
1083 cpassert(mod(nd2, nproc) == 0)
1084 mark_used(nd1)
1085
1086 !defining work arrays dimensions
1087 ncache = ncache_optimal
1088 IF (ncache <= max(n1, n2, n3/2)*4) ncache = max(n1, n2, n3/2)*4
1089 lzt = n2
1090 IF (mod(n2, 2) == 0) lzt = lzt + 1
1091
1092 !Allocations
1093 CALL fft_alloc(zw1, [ncache/4])
1094 zw1 = cmplx(0.0_dp, 0.0_dp, kind=dp)
1095 CALL fft_alloc(zw2, [ncache/4])
1096 zw2 = cmplx(0.0_dp, 0.0_dp, kind=dp)
1097 CALL fft_alloc(zt, [lzt, n1])
1098 zt = cmplx(0.0_dp, 0.0_dp, kind=dp)
1099 ALLOCATE (zmpi2(n1, nd2/nproc, nd3), source=cmplx(0.0_dp, 0.0_dp, dp))
1100 ALLOCATE (cosinarr(2, n3/2))
1101 IF (nproc > 1) ALLOCATE (zmpi1(n1, nd2/nproc, nd3/nproc, nproc), source=cmplx(0.0_dp, 0.0_dp, dp))
1102
1103 !Calculating array of phases for HalFFT decoding
1104 twopion = 8._dp*atan(1._dp)/real(n3, kind=dp)
1105 DO i3 = 1, n3/2
1106 cosinarr(1, i3) = cos(twopion*(i3 - 1))
1107 cosinarr(2, i3) = -sin(twopion*(i3 - 1))
1108 END DO
1109
1110 !transform along z axis
1111
1112 lot1 = ncache/(4*n1)
1113 lot2 = ncache/(4*n2)
1114 lot3 = ncache/(2*n3)
1115
1116 ! The size of the last chunk
1117 final_chunk_size1 = mod(n2, lot1)
1118 final_chunk_size2 = mod(nk1, lot2)
1119 final_chunk_size3 = mod(n1, lot3)
1120
1121 ! Prevent OOB-access if n2 < lot1
1122 IF (n2 >= lot1) THEN
1123 CALL fft_create_plan_1d(fft_plan_bw1, bwfft, .true., .true., lot1, lzt, n1, lot1, zw1, zt)
1124 END IF
1125 IF (final_chunk_size1 > 0) THEN
1126 CALL fft_create_plan_1d(fft_plan_bw1_last, bwfft, .true., .true., lot1, lzt, n1, &
1127 final_chunk_size1, zw1, zt)
1128 END IF
1129
1130 ! Prevent OOB-access if n2 < lot1
1131 IF (nk1 >= lot2) THEN
1132 CALL fft_create_plan_1d(fft_plan_bw2, bwfft, .true., .true., lot2, lot2, n2, lot2, zw1, zw2)
1133 END IF
1134 IF (final_chunk_size2 > 0) THEN
1135 CALL fft_create_plan_1d(fft_plan_bw2_last, bwfft, .true., .true., lot2, lot2, n2, &
1136 final_chunk_size2, zw1, zw2)
1137 END IF
1138
1139 ! Prevent OOB-access if n1 < lot3
1140 IF (n1 >= lot3) THEN
1141 CALL fft_create_plan_1d(fft_plan_bw3, bwfft, .true., .true., lot3, lot3, n3/2, lot3, zw1, zw2)
1142 END IF
1143 IF (final_chunk_size3 > 0) THEN
1144 CALL fft_create_plan_1d(fft_plan_bw3_last, bwfft, .true., .true., &
1145 lot3, lot3, n3/2, final_chunk_size3, zw1, zw2)
1146 END IF
1147
1148 DO j2 = 1, nd2/nproc
1149 !this condition ensures that we manage only the interesting part for the FFT
1150 IF (iproc*(nd2/nproc) + j2 <= n2) THEN
1151 DO i1 = 1, n1, lot3
1152 ma = i1
1153 mb = min(i1 + (lot3 - 1), n1)
1154 nfft = mb - ma + 1
1155
1156 !inserting real data into complex array of half length
1157 !input: I1,I3,J2,(Jp2)
1158
1159 CALL inserthalf(n1, n3, lot3, nfft, i1, zf(1, 1, j2), zw1)
1160
1161 !performing FFT
1162 ! This is equivalent to
1163 IF (nfft == lot3) THEN
1164 CALL fft_1d(fft_plan_bw3, zw1, zw2, 1.0_dp, stat)
1165 ELSE
1166 CALL fft_1d(fft_plan_bw3_last, zw1, zw2, 1.0_dp, stat)
1167 END IF
1168 !output: I1,i3,J2,(Jp2)
1169
1170 !unpacking FFT in order to restore correct result,
1171 !while exchanging components
1172 !input: I1,i3,J2,(Jp2)
1173 CALL scramble_unpack(i1, j2, lot3, nfft, n1, n3, nd2, nproc, nd3, zw2, zmpi2, cosinarr)
1174 !output: I1,J2,i3,(Jp2)
1175 END DO
1176 END IF
1177 END DO
1178
1179 !Interprocessor data transposition
1180 !input: I1,J2,j3,jp3,(Jp2)
1181 IF (nproc > 1) THEN
1182 !communication scheduling
1183 CALL mpi_group%alltoall(zmpi2, zmpi1, n1*(nd2/nproc)*(nd3/nproc))
1184 ! output: I1,J2,j3,Jp2,(jp3)
1185 END IF
1186
1187 DO j3 = 1, nd3/nproc
1188 !this condition ensures that we manage only the interesting part for the FFT
1189 IF (iproc*(nd3/nproc) + j3 <= n3/2 + 1) THEN
1190 jp2st = 1
1191 j2st = 1
1192
1193 !transform along x axis
1194
1195 DO j = 1, n2, lot1
1196 ma = j
1197 mb = min(j + (lot1 - 1), n2)
1198 nfft = mb - ma + 1
1199
1200 !reverse ordering
1201 !input: I1,J2,j3,Jp2,(jp3)
1202 IF (nproc == 1) THEN
1203 CALL mpiswitch(j3, nfft, jp2st, j2st, lot1, n1, nd2, nd3, nproc, zmpi2, zw1)
1204 ELSE
1205 CALL mpiswitch(j3, nfft, jp2st, j2st, lot1, n1, nd2, nd3, nproc, zmpi1, zw1)
1206 END IF
1207 !output: J2,Jp2,I1,j3,(jp3)
1208
1209 !performing FFT
1210 !input: I2,I1,j3,(jp3)
1211 ! This is equivalent to
1212 IF (nfft == lot1) THEN
1213 CALL fft_1d(fft_plan_bw1, zw1, zt(j:, 1), 1.0_dp, stat)
1214 ELSE
1215 CALL fft_1d(fft_plan_bw1_last, zw1, zt(j:, 1), 1.0_dp, stat)
1216 END IF
1217 !output: I2,i1,j3,(jp3)
1218 END DO
1219
1220 !transform along y axis, and taking only the first half
1221
1222 DO j = 1, nk1, lot2
1223 ma = j
1224 mb = min(j + (lot2 - 1), nk1)
1225 nfft = mb - ma + 1
1226
1227 !reverse ordering
1228 !input: I2,i1,j3,(jp3)
1229 CALL switch(nfft, n2, lot2, n1, lzt, zt(:, j), zw1)
1230 !output: i1,I2,j3,(jp3)
1231
1232 !performing FFT
1233 !input: i1,I2,j3,(jp3)
1234 ! This is equivalent to
1235 IF (nfft == lot2) THEN
1236 CALL fft_1d(fft_plan_bw2, zw1, zw2, 1.0_dp, stat)
1237 ELSE
1238 CALL fft_1d(fft_plan_bw2_last, zw1, zw2, 1.0_dp, stat)
1239 END IF
1240
1241 CALL realcopy(lot2, nfft, n2, nk1, nk2, zw2, zr(j, 1, j3))
1242
1243 END DO
1244 !output: i1,i2,j3,(jp3)
1245 END IF
1246 END DO
1247
1248 IF (n2 >= lot1) CALL fft_destroy_plan(fft_plan_bw1)
1249 IF (final_chunk_size1 > 0) CALL fft_destroy_plan(fft_plan_bw1_last)
1250
1251 IF (nk1 >= lot2) CALL fft_destroy_plan(fft_plan_bw2)
1252 IF (final_chunk_size2 > 0) CALL fft_destroy_plan(fft_plan_bw2_last)
1253
1254 IF (n1 >= lot3) CALL fft_destroy_plan(fft_plan_bw3)
1255 IF (final_chunk_size3 > 0) CALL fft_destroy_plan(fft_plan_bw3_last)
1256
1257 !De-allocations
1258 DEALLOCATE (zmpi2)
1259 CALL fft_dealloc(zw1)
1260 CALL fft_dealloc(zw2)
1261 CALL fft_dealloc(zt)
1262 DEALLOCATE (cosinarr)
1263 IF (nproc > 1) DEALLOCATE (zmpi1)
1264
1265 END SUBROUTINE kernelfft
1266
1267! **************************************************************************************************
1268!> \brief ...
1269!> \param lot ...
1270!> \param nfft ...
1271!> \param n2 ...
1272!> \param nk1 ...
1273!> \param nk2 ...
1274!> \param zin ...
1275!> \param zout ...
1276! **************************************************************************************************
1277 PURE SUBROUTINE realcopy(lot, nfft, n2, nk1, nk2, zin, zout)
1278 INTEGER, INTENT(in) :: lot, nfft, n2, nk1, nk2
1279 COMPLEX(KIND=dp), DIMENSION(lot, n2), INTENT(in) :: zin
1280 REAL(kind=dp), DIMENSION(nk1, nk2), INTENT(inout) :: zout
1281
1282 INTEGER :: i, j
1283
1284 DO i = 1, nk2
1285 DO j = 1, nfft
1286 zout(j, i) = real(zin(j, i), dp)
1287 END DO
1288 END DO
1289
1290 END SUBROUTINE realcopy
1291
1292! **************************************************************************************************
1293!> \brief ...
1294!> \param nfft ...
1295!> \param n2 ...
1296!> \param lot ...
1297!> \param n1 ...
1298!> \param lzt ...
1299!> \param zt ...
1300!> \param zw ...
1301! **************************************************************************************************
1302 PURE SUBROUTINE switch(nfft, n2, lot, n1, lzt, zt, zw)
1303 INTEGER, INTENT(IN) :: nfft, n2, lot, n1, lzt
1304 COMPLEX(KIND=dp), INTENT(IN) :: zt(lzt, n1)
1305 COMPLEX(KIND=dp), INTENT(INOUT) :: zw(lot, n2)
1306
1307 INTEGER :: i, j
1308
1309 DO j = 1, nfft
1310 DO i = 1, n2
1311 zw(j, i) = zt(i, j)
1312 END DO
1313 END DO
1314 END SUBROUTINE switch
1315
1316! **************************************************************************************************
1317!> \brief ...
1318!> \param j3 ...
1319!> \param nfft ...
1320!> \param Jp2st ...
1321!> \param J2st ...
1322!> \param lot ...
1323!> \param n1 ...
1324!> \param nd2 ...
1325!> \param nd3 ...
1326!> \param nproc ...
1327!> \param zmpi1 ...
1328!> \param zw ...
1329! **************************************************************************************************
1330 PURE SUBROUTINE mpiswitch(j3, nfft, Jp2st, J2st, lot, n1, nd2, nd3, nproc, zmpi1, zw)
1331 INTEGER, INTENT(IN) :: j3, nfft
1332 INTEGER, INTENT(INOUT) :: jp2st, j2st
1333 INTEGER, INTENT(IN) :: lot, n1, nd2, nd3, nproc
1334 COMPLEX(KIND=dp), INTENT(IN) :: zmpi1(n1, nd2/nproc, nd3/nproc, nproc)
1335 COMPLEX(KIND=dp), INTENT(INOUT) :: zw(lot, n1)
1336
1337 INTEGER :: i1, j2, jp2, mfft
1338
1339 mfft = 0
1340 DO jp2 = jp2st, nproc
1341 DO j2 = j2st, nd2/nproc
1342 mfft = mfft + 1
1343 IF (mfft > nfft) THEN
1344 jp2st = jp2
1345 j2st = j2
1346 RETURN
1347 END IF
1348 DO i1 = 1, n1
1349 zw(mfft, i1) = zmpi1(i1, j2, j3, jp2)
1350 END DO
1351 END DO
1352 j2st = 1
1353 END DO
1354 END SUBROUTINE mpiswitch
1355
1356! **************************************************************************************************
1357!> \brief ...
1358!> \param p ...
1359!> \param w ...
1360!> \param urange ...
1361!> \param drange ...
1362!> \param acc ...
1363! **************************************************************************************************
1364 PURE SUBROUTINE gequad(p, w, urange, drange, acc)
1365!
1366 REAL(kind=dp), INTENT(OUT) :: p(89), w(89), urange, drange, acc
1367
1368!
1369!
1370! range [10^(-9),1] and accuracy ~10^(-8);
1371!
1372!
1373
1374 p(1) = 4.96142640560223544e19_dp
1375 p(2) = 1.37454269147978052e19_dp
1376 p(3) = 7.58610013441204679e18_dp
1377 p(4) = 4.42040691347806996e18_dp
1378 p(5) = 2.61986077948367892e18_dp
1379 p(6) = 1.56320138155496681e18_dp
1380 p(7) = 9.35645215863028402e17_dp
1381 p(8) = 5.60962910452691703e17_dp
1382 p(9) = 3.3666225119686761e17_dp
1383 p(10) = 2.0218253197947866e17_dp
1384 p(11) = 1.21477756091902017e17_dp
1385 p(12) = 7.3012982513608503e16_dp
1386 p(13) = 4.38951893556421099e16_dp
1387 p(14) = 2.63949482512262325e16_dp
1388 p(15) = 1.58742054072786174e16_dp
1389 p(16) = 9.54806587737665531e15_dp
1390 p(17) = 5.74353712364571709e15_dp
1391 p(18) = 3.455214877389445e15_dp
1392 p(19) = 2.07871658520326804e15_dp
1393 p(20) = 1.25064667315629928e15_dp
1394 p(21) = 7.52469429541933745e14_dp
1395 p(22) = 4.5274603337253175e14_dp
1396 p(23) = 2.72414006900059548e14_dp
1397 p(24) = 1.63912168349216752e14_dp
1398 p(25) = 9.86275802590865738e13_dp
1399 p(26) = 5.93457701624974985e13_dp
1400 p(27) = 3.5709554322296296e13_dp
1401 p(28) = 2.14872890367310454e13_dp
1402 p(29) = 1.29294719957726902e13_dp
1403 p(30) = 7.78003375426361016e12_dp
1404 p(31) = 4.68148199759876704e12_dp
1405 p(32) = 2.8169955024829868e12_dp
1406 p(33) = 1.69507790481958464e12_dp
1407 p(34) = 1.01998486064607581e12_dp
1408 p(35) = 6.13759486539856459e11_dp
1409 p(36) = 3.69320183828682544e11_dp
1410 p(37) = 2.22232783898905102e11_dp
1411 p(38) = 1.33725247623668682e11_dp
1412 p(39) = 8.0467192739036288e10_dp
1413 p(40) = 4.84199582415144143e10_dp
1414 p(41) = 2.91360091170559564e10_dp
1415 p(42) = 1.75321747475309216e10_dp
1416 p(43) = 1.0549735552210995e10_dp
1417 p(44) = 6.34815321079006586e9_dp
1418 p(45) = 3.81991113733594231e9_dp
1419 p(46) = 2.29857747533101109e9_dp
1420 p(47) = 1.38313653595483694e9_dp
1421 p(48) = 8.32282908580025358e8_dp
1422 p(49) = 5.00814519374587467e8_dp
1423 p(50) = 3.01358090773319025e8_dp
1424 p(51) = 1.81337994217503535e8_dp
1425 p(52) = 1.09117589961086823e8_dp
1426 p(53) = 6.56599771718640323e7_dp
1427 p(54) = 3.95099693638497164e7_dp
1428 p(55) = 2.37745694710665991e7_dp
1429 p(56) = 1.43060135285912813e7_dp
1430 p(57) = 8.60844290313506695e6_dp
1431 p(58) = 5.18000974075383424e6_dp
1432 p(59) = 3.116998193057466e6_dp
1433 p(60) = 1.87560993870024029e6_dp
1434 p(61) = 1.12862197183979562e6_dp
1435 p(62) = 679132.441326077231_dp
1436 p(63) = 408658.421279877969_dp
1437 p(64) = 245904.473450669789_dp
1438 p(65) = 147969.568088321005_dp
1439 p(66) = 89038.612357311147_dp
1440 p(67) = 53577.7362552358895_dp
1441 p(68) = 32239.6513926914668_dp
1442 p(69) = 19399.7580852362791_dp
1443 p(70) = 11673.5323603058634_dp
1444 p(71) = 7024.38438577707758_dp
1445 p(72) = 4226.82479307685999_dp
1446 p(73) = 2543.43254175354295_dp
1447 p(74) = 1530.47486269122675_dp
1448 p(75) = 920.941785160749482_dp
1449 p(76) = 554.163803906291646_dp
1450 p(77) = 333.46029740785694_dp
1451 p(78) = 200.6550575335041_dp
1452 p(79) = 120.741366914147284_dp
1453 p(80) = 72.6544243200329916_dp
1454 p(81) = 43.7187810415471025_dp
1455 p(82) = 26.3071631447061043_dp
1456 p(83) = 15.8299486353816329_dp
1457 p(84) = 9.52493152341244004_dp
1458 p(85) = 5.72200417067776041_dp
1459 p(86) = 3.36242234070940928_dp
1460 p(87) = 1.75371394604499472_dp
1461 p(88) = 0.64705932650658966_dp
1462 p(89) = 0.072765905943708247_dp
1463 !
1464 w(1) = 47.67445484528304247e10_dp
1465 w(2) = 11.37485774750442175e9_dp
1466 w(3) = 78.64340976880190239e8_dp
1467 w(4) = 46.27335788759590498e8_dp
1468 w(5) = 24.7380464827152951e8_dp
1469 w(6) = 13.62904116438987719e8_dp
1470 w(7) = 92.79560029045882433e8_dp
1471 w(8) = 52.15931216254660251e8_dp
1472 w(9) = 31.67018011061666244e8_dp
1473 w(10) = 1.29291036801493046e8_dp
1474 w(11) = 1.00139319988015862e8_dp
1475 w(12) = 7.75892350510188341e7_dp
1476 w(13) = 6.01333567950731271e7_dp
1477 w(14) = 4.66141178654796875e7_dp
1478 w(15) = 3.61398903394911448e7_dp
1479 w(16) = 2.80225846672956389e7_dp
1480 w(17) = 2.1730509180930247e7_dp
1481 w(18) = 1.68524482625876965e7_dp
1482 w(19) = 1.30701489345870338e7_dp
1483 w(20) = 1.01371784832269282e7_dp
1484 w(21) = 7.86264116300379329e6_dp
1485 w(22) = 6.09861667912273717e6_dp
1486 w(23) = 4.73045784039455683e6_dp
1487 w(24) = 3.66928949951594161e6_dp
1488 w(25) = 2.8462050836230259e6_dp
1489 w(26) = 2.20777394798527011e6_dp
1490 w(27) = 1.71256191589205524e6_dp
1491 w(28) = 1.32843556197737076e6_dp
1492 w(29) = 1.0304731275955989e6_dp
1493 w(30) = 799345.206572271448_dp
1494 w(31) = 620059.354143595343_dp
1495 w(32) = 480986.704107449333_dp
1496 w(33) = 373107.167700228515_dp
1497 w(34) = 289424.08337412132_dp
1498 w(35) = 224510.248231581788_dp
1499 w(36) = 174155.825690028966_dp
1500 w(37) = 135095.256919654065_dp
1501 w(38) = 104795.442776800312_dp
1502 w(39) = 81291.4458222430418_dp
1503 w(40) = 63059.0493649328682_dp
1504 w(41) = 48915.9040455329689_dp
1505 w(42) = 37944.8484018048756_dp
1506 w(43) = 29434.4290473253969_dp
1507 w(44) = 22832.7622054490044_dp
1508 w(45) = 17711.743950151233_dp
1509 w(46) = 13739.287867104177_dp
1510 w(47) = 10657.7895710752585_dp
1511 w(48) = 8267.42141053961834_dp
1512 w(49) = 6413.17397520136448_dp
1513 w(50) = 4974.80402838654277_dp
1514 w(51) = 3859.03698188553047_dp
1515 w(52) = 2993.51824493299154_dp
1516 w(53) = 2322.1211966811754_dp
1517 w(54) = 1801.30750964719641_dp
1518 w(55) = 1397.30379659817038_dp
1519 w(56) = 1083.91149143250697_dp
1520 w(57) = 840.807939169209188_dp
1521 w(58) = 652.228524366749422_dp
1522 w(59) = 505.944376983506128_dp
1523 w(60) = 392.469362317941064_dp
1524 w(61) = 304.444930257324312_dp
1525 w(62) = 236.162932842453601_dp
1526 w(63) = 183.195466078603525_dp
1527 w(64) = 142.107732186551471_dp
1528 w(65) = 110.23530215723992_dp
1529 w(66) = 85.5113346705382257_dp
1530 w(67) = 66.3325469806696621_dp
1531 w(68) = 51.4552463353841373_dp
1532 w(69) = 39.9146798429449273_dp
1533 w(70) = 30.9624728409162095_dp
1534 w(71) = 24.018098812215013_dp
1535 w(72) = 18.6312338024296588_dp
1536 w(73) = 14.4525541233150501_dp
1537 w(74) = 11.2110836519105938_dp
1538 w(75) = 8.69662175848497178_dp
1539 w(76) = 6.74611236165731961_dp
1540 w(77) = 5.23307018057529994_dp
1541 w(78) = 4.05937850501539556_dp
1542 w(79) = 3.14892659076635714_dp
1543 w(80) = 2.44267408211071604_dp
1544 w(81) = 1.89482240522855261_dp
1545 w(82) = 1.46984505907050079_dp
1546 w(83) = 1.14019261330527007_dp
1547 w(84) = 0.884791217422925293_dp
1548 w(85) = 0.692686387080616483_dp
1549 w(86) = 0.585244576897023282_dp
1550 w(87) = 0.576182522545327589_dp
1551 w(88) = 0.596688817388997178_dp
1552 w(89) = 0.607879901151108771_dp
1553 !
1554 !
1555 urange = 1._dp
1556 drange = 1e-08_dp
1557 acc = 1e-08_dp
1558 END SUBROUTINE gequad
1559
1560END MODULE ps_wavelet_kernel
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:126
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:256
fft_lib::fft_1d
subroutine, public fft_1d(plan, zin, zout, scale, stat)
...
Definition fft_lib.F:276
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:213
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
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::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_kernel
Creates the wavelet kernel for the wavelet based poisson solver.
Definition ps_wavelet_kernel.F:12
ps_wavelet_kernel::createkernel
subroutine, public createkernel(geocode, n01, n02, n03, hx, hy, hz, itype_scf, iproc, nproc, kernel, mpi_group)
Allocate a pointer which corresponds to the zero-padded FFT slice needed for calculating the convolut...
Definition ps_wavelet_kernel.F:84
ps_wavelet_scaling_function
Creates the wavelet kernel for the wavelet based poisson solver.
Definition ps_wavelet_scaling_function.F:12
ps_wavelet_scaling_function::scaling_function
subroutine, public scaling_function(itype, nd, nrange, a, x)
Calculate the values of a scaling function in real uniform grid.
Definition ps_wavelet_scaling_function.F:35
ps_wavelet_scaling_function::scf_recursion
pure subroutine, public scf_recursion(itype, n_iter, n_range, kernel_scf, kern_1_scf)
Do iterations to go from p0gauss to pgauss order interpolating scaling function.
Definition ps_wavelet_scaling_function.F:144
ps_wavelet_util
Performs a wavelet based solution of the Poisson equation.
Definition ps_wavelet_util.F:12
ps_wavelet_util::s_fft_dimensions
subroutine, public s_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
Calculate four sets of dimension needed for the calculation of the convolution for the surface system...
Definition ps_wavelet_util.F:326
ps_wavelet_util::f_fft_dimensions
subroutine, public f_fft_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
Calculate four sets of dimension needed for the calculation of the zero-padded convolution.
Definition ps_wavelet_util.F:432
fft_plan::fft_plan_type
Definition fft_plan.F:27
message_passing::mp_comm_type
Definition message_passing.F:179