Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Performs a wavelet based solution of the Poisson equation.
10 : !> \author Florian Schiffmann (09.2007,fschiff)
11 : ! **************************************************************************************************
12 : MODULE ps_wavelet_util
13 : USE fft_tools, ONLY: FFT_RADIX_NEXT,&
14 : fft_radix_operations
15 : USE kinds, ONLY: dp
16 : USE mathconstants, ONLY: fourpi
17 : USE ps_wavelet_base, ONLY: f_poissonsolver,&
18 : p_poissonsolver,&
19 : s_poissonsolver
20 : USE pw_grid_types, ONLY: pw_grid_type
21 : #include "../base/base_uses.f90"
22 :
23 : IMPLICIT NONE
24 :
25 : PRIVATE
26 :
27 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ps_wavelet_util'
28 :
29 : ! *** Public data types ***
30 :
31 : PUBLIC :: PSolver, &
32 : P_FFT_dimensions, &
33 : S_FFT_dimensions, &
34 : F_FFT_dimensions
35 :
36 : CONTAINS
37 :
38 : ! **************************************************************************************************
39 : !> \brief Calculate the Poisson equation $\nabla^2 V(x,y,z)=-4 \pi \rho(x,y,z)$
40 : !> from a given $\rho$, for different boundary conditions an for different data distributions.
41 : !> Following the boundary conditions, it applies the Poisson Kernel previously calculated.
42 : !> \param geocode Indicates the boundary conditions (BC) of the problem:
43 : !> 'F' free BC, isolated systems.
44 : !> The program calculates the solution as if the given density is
45 : !> "alone" in R^3 space.
46 : !> 'S' surface BC, isolated in y direction, periodic in xz plane
47 : !> The given density is supposed to be periodic in the xz plane,
48 : !> so the dimensions in these direction mus be compatible with the FFT
49 : !> Beware of the fact that the isolated direction is y!
50 : !> 'P' periodic BC.
51 : !> The density is supposed to be periodic in all the three directions,
52 : !> then all the dimensions must be compatible with the FFT.
53 : !> No need for setting up the kernel.
54 : !> \param iproc label of the process,from 0 to nproc-1
55 : !> \param nproc number of processors
56 : !> \param n01 global dimension in the three directions.
57 : !> \param n02 global dimension in the three directions.
58 : !> \param n03 global dimension in the three directions.
59 : !> \param hx grid spacings. For the isolated BC case for the moment they are supposed to
60 : !> be equal in the three directions
61 : !> \param hy grid spacings. For the isolated BC case for the moment they are supposed to
62 : !> be equal in the three directions
63 : !> \param hz grid spacings. For the isolated BC case for the moment they are supposed to
64 : !> be equal in the three directions
65 : !> \param rhopot main input/output array.
66 : !> On input, it represents the density values on the grid points
67 : !> On output, it is the Hartree potential, namely the solution of the Poisson
68 : !> equation PLUS (when ixc/=0) the XC potential PLUS (again for ixc/=0) the
69 : !> pot_ion array. The output is non overlapping, in the sense that it does not
70 : !> consider the points that are related to gradient and WB calculation
71 : !> \param karray kernel of the poisson equation. It is provided in distributed case, with
72 : !> dimensions that are related to the output of the PS_dim4allocation routine
73 : !> it MUST be created by following the same geocode as the Poisson Solver.
74 : !> \param pw_grid ...
75 : !> \date February 2007
76 : !> \author Luigi Genovese
77 : !> \note The dimensions of the arrays must be compatible with geocode, nproc,
78 : !> ixc and iproc. Since the arguments of these routines are indicated with the *, it
79 : !> is IMPERATIVE to use the PS_dim4allocation routine for calculation arrays sizes.
80 : ! **************************************************************************************************
81 34299 : SUBROUTINE PSolver(geocode, iproc, nproc, n01, n02, n03, hx, hy, hz, &
82 : rhopot, karray, pw_grid)
83 : CHARACTER(len=1), INTENT(in) :: geocode
84 : INTEGER, INTENT(in) :: iproc, nproc, n01, n02, n03
85 : REAL(KIND=dp), INTENT(in) :: hx, hy, hz
86 : REAL(KIND=dp), DIMENSION(*), INTENT(inout) :: rhopot
87 : REAL(KIND=dp), DIMENSION(*), INTENT(in) :: karray
88 : TYPE(pw_grid_type), POINTER :: pw_grid
89 :
90 : INTEGER :: i1, i2, i3, iend, istart, j2, m1, m2, &
91 : m3, md1, md2, md3, n1, n2, n3, nd1, &
92 : nd2, nd3, nlim, nwb, nwbl, nwbr, nxc, &
93 : nxcl, nxcr, nxt
94 : REAL(KIND=dp) :: factor, hgrid, red_fact, scal
95 34299 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: zf
96 :
97 : !the order of the finite-difference gradient (fixed)
98 : !calculate the dimensions wrt the geocode
99 :
100 34299 : IF (geocode == 'P') THEN
101 17323 : CALL P_FFT_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
102 16976 : ELSE IF (geocode == 'S') THEN
103 54 : CALL S_FFT_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
104 16922 : ELSE IF (geocode == 'F') THEN
105 16922 : CALL F_FFT_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
106 : ELSE
107 0 : CPABORT("PSolver: geometry code not admitted")
108 : END IF
109 : !array allocations
110 171495 : ALLOCATE (zf(md1, md3, md2/nproc))
111 :
112 : !dimension for exchange-correlation (different in the global or distributed case)
113 : !let us calculate the dimension of the portion of the rhopot array to be passed
114 : !to the xc routine
115 : !this portion will depend on the need of calculating the gradient or not,
116 : !and whether the White-Bird correction must be inserted or not
117 : !(absent only in the LB ixc=13 case)
118 :
119 : !nxc is the effective part of the third dimension that is being processed
120 : !nxt is the dimension of the part of rhopot that must be passed to the gradient routine
121 : !nwb is the dimension of the part of rhopot in the wb-postprocessing routine
122 : !note: nxc <= nwb <= nxt
123 : !the dimension are related by the values of nwbl and nwbr
124 : ! nxc+nxcl+nxcr-2 = nwb
125 : ! nwb+nwbl+nwbr = nxt
126 34299 : istart = iproc*(md2/nproc)
127 34299 : iend = MIN((iproc + 1)*md2/nproc, m2)
128 :
129 34299 : nxc = iend - istart
130 34299 : nwbl = 0
131 34299 : nwbr = 0
132 34299 : nxcl = 1
133 34299 : nxcr = 1
134 :
135 34299 : nwb = nxcl + nxc + nxcr - 2
136 34299 : nxt = nwbr + nwb + nwbl
137 :
138 : !calculate the actual limit of the array for the zero padded FFT
139 34299 : IF (geocode == 'P') THEN
140 17323 : nlim = n2
141 16976 : ELSE IF (geocode == 'S') THEN
142 54 : nlim = n2
143 16922 : ELSE IF (geocode == 'F') THEN
144 16922 : nlim = n2/2
145 : END IF
146 :
147 34299 : IF (istart + 1 <= m2) THEN
148 34299 : red_fact = 1._dp
149 34299 : CALL scale_and_distribute(m1, m3, md1, md2, md3, nxc, rhopot, zf, nproc, red_fact)
150 0 : ELSE IF (istart + 1 <= nlim) THEN !this condition assures that we have perform good zero padding
151 0 : DO i2 = istart + 1, MIN(nlim, istart + md2/nproc)
152 0 : j2 = i2 - istart
153 0 : DO i3 = 1, md3
154 0 : DO i1 = 1, md1
155 0 : zf(i1, i3, j2) = 0._dp
156 : END DO
157 : END DO
158 : END DO
159 : END IF
160 :
161 : !this routine builds the values for each process of the potential (zf), multiplying by scal
162 34299 : IF (geocode == 'P') THEN
163 : !no powers of hgrid because they are incorporated in the plane wave treatment
164 17323 : scal = 1._dp/REAL(n1*n2*n3, KIND=dp)
165 : CALL P_PoissonSolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, zf, &
166 17323 : scal, hx, hy, hz, pw_grid%para%group)
167 16976 : ELSE IF (geocode == 'S') THEN
168 : !only one power of hgrid
169 54 : scal = hy/REAL(n1*n2*n3, KIND=dp)
170 : CALL S_PoissonSolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, karray, zf, &
171 54 : scal, pw_grid%para%group)
172 16922 : ELSE IF (geocode == 'F') THEN
173 16922 : hgrid = MAX(hx, hy, hz)
174 16922 : scal = hgrid**3/REAL(n1*n2*n3, KIND=dp)
175 : CALL F_PoissonSolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, karray, zf, &
176 16922 : scal, pw_grid%para%group)
177 16922 : factor = 0.5_dp*hgrid**3
178 : END IF
179 :
180 : !the value of the shift depends on the distributed i/o or not
181 34299 : IF (geocode == 'F') THEN
182 16922 : red_fact = 1._dp
183 : ELSE
184 17377 : red_fact = -fourpi
185 : END IF
186 :
187 34299 : CALL scale_and_distribute(m1, m3, md1, md2, md3, nxc, zf, rhopot, nproc, red_fact)
188 :
189 34299 : DEALLOCATE (zf)
190 :
191 34299 : END SUBROUTINE PSolver
192 :
193 : ! **************************************************************************************************
194 : !> \brief Calculate four sets of dimension needed for the calculation of the
195 : !> convolution for the periodic system
196 : !> \param n01 original real dimensions (input)
197 : !> \param n02 original real dimensions (input)
198 : !> \param n03 original real dimensions (input)
199 : !> \param m1 original real dimension, with m2 and m3 exchanged
200 : !> \param m2 original real dimension, with m2 and m3 exchanged
201 : !> \param m3 original real dimension, with m2 and m3 exchanged
202 : !> \param n1 the first FFT dimensions, for the moment supposed to be even
203 : !> \param n2 the first FFT dimensions, for the moment supposed to be even
204 : !> \param n3 the first FFT dimensions, for the moment supposed to be even
205 : !> \param md1 the n1,n2,n3 dimensions. They contain the real unpadded space,
206 : !> properly enlarged to be compatible with the FFT dimensions n_i.
207 : !> md2 is further enlarged to be a multiple of nproc
208 : !> \param md2 the n1,n2,n3 dimensions. They contain the real unpadded space,
209 : !> properly enlarged to be compatible with the FFT dimensions n_i.
210 : !> md2 is further enlarged to be a multiple of nproc
211 : !> \param md3 the n1,n2,n3 dimensions. They contain the real unpadded space,
212 : !> properly enlarged to be compatible with the FFT dimensions n_i.
213 : !> md2 is further enlarged to be a multiple of nproc
214 : !> \param nd1 fourier dimensions for which the kernel is injective,
215 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
216 : !> enlarged to be a multiple of nproc
217 : !> \param nd2 fourier dimensions for which the kernel is injective,
218 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
219 : !> enlarged to be a multiple of nproc
220 : !> \param nd3 fourier dimensions for which the kernel is injective,
221 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
222 : !> enlarged to be a multiple of nproc
223 : !> \param nproc ...
224 : !> \date October 2006
225 : !> \author Luigi Genovese
226 : !> \note This four sets of dimensions are actually redundant (mi=n0i),
227 : !> due to the backward-compatibility
228 : !> with the other geometries of the Poisson Solver.
229 : !> The dimensions 2 and 3 are exchanged.
230 : ! **************************************************************************************************
231 17649 : SUBROUTINE P_FFT_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
232 : INTEGER, INTENT(in) :: n01, n02, n03
233 : INTEGER, INTENT(out) :: m1, m2, m3, n1, n2, n3, md1, md2, md3, &
234 : nd1, nd2, nd3
235 : INTEGER, INTENT(in) :: nproc
236 :
237 : CHARACTER(len=80) :: err
238 : INTEGER :: l1, l2, l3
239 :
240 : !dimensions of the density in the real space
241 :
242 17649 : m1 = n01
243 17649 : m2 = n03
244 17649 : m3 = n02
245 :
246 : ! real space grid dimension (suitable for number of processors)
247 17649 : l1 = m1
248 17649 : l2 = m2
249 17649 : l3 = m3 !beware of the half dimension
250 : ! Find suitable and fast radices (in principal, every radix is possible, but not every radix is reasonably fast)
251 17649 : CALL fft_radix_operations(l1, n1, FFT_RADIX_NEXT)
252 17649 : IF (n1 == m1) THEN
253 : ELSE
254 0 : WRITE (err, *) 'the FFT in the x direction is not allowed; n01 dimension ', n01
255 0 : CPABORT(TRIM(err))
256 : END IF
257 17649 : CALL fft_radix_operations(l2, n2, FFT_RADIX_NEXT)
258 17649 : IF (n2 == m2) THEN
259 : ELSE
260 0 : WRITE (err, *) 'the FFT in the z direction is not allowed; n03 dimension ', n03
261 0 : CPABORT(TRIM(err))
262 : END IF
263 17649 : CALL fft_radix_operations(l3, n3, FFT_RADIX_NEXT)
264 17649 : IF (n3 == m3) THEN
265 : ELSE
266 0 : WRITE (err, *) 'the FFT in the y direction is not allowed; n02 dimension ', n02
267 0 : CPABORT(TRIM(err))
268 : END IF
269 :
270 : !dimensions that contain the unpadded real space,
271 : ! compatible with the number of processes
272 17649 : md1 = n1
273 17649 : md2 = n2
274 17649 : md3 = n3
275 19029 : DO WHILE (nproc*(md2/nproc) < n2)
276 1380 : md2 = md2 + 1
277 : END DO
278 :
279 : !dimensions of the kernel, 1/8 of the total volume,
280 : !compatible with nproc
281 17649 : nd1 = n1/2 + 1
282 17649 : nd2 = n2/2 + 1
283 17649 : nd3 = n3/2 + 1
284 19579 : DO WHILE (MODULO(nd3, nproc) /= 0)
285 1930 : nd3 = nd3 + 1
286 : END DO
287 :
288 17649 : END SUBROUTINE P_FFT_dimensions
289 :
290 : ! **************************************************************************************************
291 : !> \brief Calculate four sets of dimension needed for the calculation of the
292 : !> convolution for the surface system
293 : !> \param n01 original real dimensions (input)
294 : !> \param n02 original real dimensions (input)
295 : !> \param n03 original real dimensions (input)
296 : !> \param m1 original real dimension, with 2 and 3 exchanged
297 : !> \param m2 original real dimension, with 2 and 3 exchanged
298 : !> \param m3 original real dimension, with 2 and 3 exchanged
299 : !> \param n1 the first FFT dimensions, for the moment supposed to be even
300 : !> \param n2 the first FFT dimensions, for the moment supposed to be even
301 : !> \param n3 the double of the first FFT even dimension greater than m3
302 : !> (improved for the HalFFT procedure)
303 : !> \param md1 the n1,n2 dimensions.
304 : !> \param md2 the n1,n2,n3 dimensions.
305 : !> \param md3 the half of n3 dimension. They contain the real unpadded space,
306 : !> properly enlarged to be compatible with the FFT dimensions n_i.
307 : !> md2 is further enlarged to be a multiple of nproc
308 : !> \param nd1 fourier dimensions for which the kernel is injective,
309 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
310 : !> enlarged to be a multiple of nproc
311 : !> \param nd2 fourier dimensions for which the kernel is injective,
312 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
313 : !> enlarged to be a multiple of nproc
314 : !> \param nd3 fourier dimensions for which the kernel is injective,
315 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
316 : !> enlarged to be a multiple of nproc
317 : !> \param nproc ...
318 : !> \date October 2006
319 : !> \author Luigi Genovese
320 : !> \note This four sets of dimensions are actually redundant (mi=n0i),
321 : !> due to the backward-compatibility
322 : !> with the Poisson Solver with other geometries.
323 : !> Dimensions n02 and n03 were exchanged
324 : ! **************************************************************************************************
325 66 : SUBROUTINE S_FFT_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
326 : INTEGER, INTENT(in) :: n01, n02, n03
327 : INTEGER, INTENT(out) :: m1, m2, m3, n1, n2, n3, md1, md2, md3, &
328 : nd1, nd2, nd3
329 : INTEGER, INTENT(in) :: nproc
330 :
331 : CHARACTER(len=*), PARAMETER :: routineN = 'S_FFT_dimensions'
332 :
333 : CHARACTER(len=80) :: err
334 : INTEGER :: handle, l1, l2, l3
335 :
336 : !dimensions of the density in the real space
337 :
338 66 : CALL timeset(routineN, handle)
339 66 : m1 = n01
340 66 : m2 = n03
341 66 : m3 = n02
342 :
343 : ! real space grid dimension (suitable for number of processors)
344 66 : l1 = m1
345 66 : l2 = m2
346 66 : l3 = m3 !beware of the half dimension
347 : ! Find suitable and fast radices (in principal, every radix is possible, but not every radix is reasonably fast)
348 66 : CALL fft_radix_operations(l1, n1, FFT_RADIX_NEXT)
349 66 : IF (n1 == m1) THEN
350 : ELSE
351 0 : WRITE (err, *) 'the FFT in the x direction is not allowed; n01 dimension', n01
352 0 : CPABORT(TRIM(err))
353 : END IF
354 66 : CALL fft_radix_operations(l2, n2, FFT_RADIX_NEXT)
355 66 : IF (n2 == m2) THEN
356 : ELSE
357 0 : WRITE (err, *) 'the FFT in the z direction is not allowed; n03 dimension', n03
358 0 : CPABORT(TRIM(err))
359 : END IF
360 0 : DO
361 66 : CALL fft_radix_operations(l3, n3, FFT_RADIX_NEXT)
362 66 : IF (MODULO(n3, 2) == 0) THEN
363 : EXIT
364 : END IF
365 0 : l3 = l3 + 1
366 : END DO
367 66 : n3 = 2*n3
368 :
369 : !dimensions that contain the unpadded real space,
370 : ! compatible with the number of processes
371 66 : md1 = n1
372 66 : md2 = n2
373 66 : md3 = n3/2
374 66 : DO WHILE (nproc*(md2/nproc) < n2)
375 0 : md2 = md2 + 1
376 : END DO
377 :
378 : !dimensions of the kernel, 1/8 of the total volume,
379 : !compatible with nproc
380 :
381 : !these two dimensions are like that since they are even
382 66 : nd1 = n1/2 + 1
383 66 : nd2 = n2/2 + 1
384 :
385 66 : nd3 = n3/2 + 1
386 132 : DO WHILE (MODULO(nd3, nproc) /= 0)
387 66 : nd3 = nd3 + 1
388 : END DO
389 66 : CALL timestop(handle)
390 :
391 66 : END SUBROUTINE S_FFT_dimensions
392 :
393 : ! **************************************************************************************************
394 : !> \brief Calculate four sets of dimension needed for the calculation of the
395 : !> zero-padded convolution
396 : !> \param n01 original real dimensions (input)
397 : !> \param n02 original real dimensions (input)
398 : !> \param n03 original real dimensions (input)
399 : !> \param m1 original real dimension with the dimension 2 and 3 exchanged
400 : !> \param m2 original real dimension with the dimension 2 and 3 exchanged
401 : !> \param m3 original real dimension with the dimension 2 and 3 exchanged
402 : !> \param n1 ...
403 : !> \param n2 ...
404 : !> \param n3 the double of the first FFT even dimension greater than m3
405 : !> (improved for the HalFFT procedure)
406 : !> \param md1 half of n1,n2,n3 dimension. They contain the real unpadded space,
407 : !> properly enlarged to be compatible with the FFT dimensions n_i.
408 : !> md2 is further enlarged to be a multiple of nproc
409 : !> \param md2 half of n1,n2,n3 dimension. They contain the real unpadded space,
410 : !> properly enlarged to be compatible with the FFT dimensions n_i.
411 : !> md2 is further enlarged to be a multiple of nproc
412 : !> \param md3 half of n1,n2,n3 dimension. They contain the real unpadded space,
413 : !> properly enlarged to be compatible with the FFT dimensions n_i.
414 : !> md2 is further enlarged to be a multiple of nproc
415 : !> \param nd1 fourier dimensions for which the kernel FFT is injective,
416 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
417 : !> enlarged to be a multiple of nproc
418 : !> \param nd2 fourier dimensions for which the kernel FFT is injective,
419 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
420 : !> enlarged to be a multiple of nproc
421 : !> \param nd3 fourier dimensions for which the kernel FFT is injective,
422 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
423 : !> enlarged to be a multiple of nproc
424 : !> \param nproc ...
425 : !> \date February 2006
426 : !> \author Luigi Genovese
427 : !> \note The dimension m2 and m3 correspond to n03 and n02 respectively
428 : !> this is needed since the convolution routine manage arrays of dimension
429 : !> (md1,md3,md2/nproc)
430 : ! **************************************************************************************************
431 18436 : SUBROUTINE F_FFT_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
432 : INTEGER, INTENT(in) :: n01, n02, n03
433 : INTEGER, INTENT(out) :: m1, m2, m3, n1, n2, n3, md1, md2, md3, &
434 : nd1, nd2, nd3
435 : INTEGER, INTENT(in) :: nproc
436 :
437 : INTEGER :: l1, l2, l3
438 :
439 : !dimensions of the density in the real space, inverted for convenience
440 :
441 18436 : m1 = n01
442 18436 : m2 = n03
443 18436 : m3 = n02
444 : ! real space grid dimension (suitable for number of processors)
445 18436 : l1 = 2*m1
446 18436 : l2 = 2*m2
447 18436 : l3 = m3 !beware of the half dimension
448 : ! Find suitable and fast radices (in principal, every radix is possible, but not every radix is reasonably fast)
449 0 : DO
450 18436 : CALL fft_radix_operations(l1, n1, FFT_RADIX_NEXT)
451 18436 : IF (MODULO(n1, 2) == 0) THEN
452 : EXIT
453 : END IF
454 18436 : l1 = l1 + 1
455 : END DO
456 0 : DO
457 18436 : CALL fft_radix_operations(l2, n2, FFT_RADIX_NEXT)
458 18436 : IF (MODULO(n2, 2) == 0) THEN
459 : EXIT
460 : END IF
461 18436 : l2 = l2 + 1
462 : END DO
463 10800 : DO
464 29236 : CALL fft_radix_operations(l3, n3, FFT_RADIX_NEXT)
465 29236 : IF (MODULO(n3, 2) == 0) THEN
466 : EXIT
467 : END IF
468 10800 : l3 = l3 + 1
469 : END DO
470 18436 : n3 = 2*n3
471 :
472 : !dimensions that contain the unpadded real space,
473 : ! compatible with the number of processes
474 18436 : md1 = n1/2
475 18436 : md2 = n2/2
476 18436 : md3 = n3/2
477 21794 : DO WHILE (nproc*(md2/nproc) < n2/2)
478 3358 : md2 = md2 + 1
479 : END DO
480 :
481 : !dimensions of the kernel, 1/8 of the total volume,
482 : !compatible with nproc
483 18436 : nd1 = n1/2 + 1
484 18436 : nd2 = n2/2 + 1
485 18436 : nd3 = n3/2 + 1
486 :
487 28744 : DO WHILE (MODULO(nd3, nproc) /= 0)
488 10308 : nd3 = nd3 + 1
489 : END DO
490 :
491 18436 : END SUBROUTINE F_FFT_dimensions
492 :
493 : ! **************************************************************************************************
494 : !> \brief ...
495 : !> \param m1 ...
496 : !> \param m3 ...
497 : !> \param md1 ...
498 : !> \param md2 ...
499 : !> \param md3 ...
500 : !> \param nxc ...
501 : !> \param rhopot ...
502 : !> \param zf ...
503 : !> \param nproc ...
504 : !> \param factor ...
505 : ! **************************************************************************************************
506 68598 : SUBROUTINE scale_and_distribute(m1, m3, md1, md2, md3, nxc, &
507 68598 : rhopot, zf, nproc, factor)
508 :
509 : !Arguments----------------------
510 : INTEGER, INTENT(in) :: m1, m3, md1, md2, md3, nxc, nproc
511 : REAL(KIND=dp), DIMENSION(md1, md3, md2/nproc), &
512 : INTENT(inout) :: zf, rhopot
513 : REAL(KIND=dp), INTENT(in) :: factor
514 :
515 : CHARACTER(len=*), PARAMETER :: routineN = 'scale_and_distribute'
516 :
517 : INTEGER :: handle, j1, j3, jp2
518 :
519 68598 : CALL timeset(routineN, handle)
520 :
521 68598 : IF (nxc >= 1) THEN
522 1560472 : DO jp2 = 1, nxc
523 52947816 : DO j3 = 1, m3
524 2276481120 : DO j1 = 1, m1
525 2276481120 : zf(j1, j3, jp2) = factor*rhopot(j1, j3, jp2)
526 : END DO
527 58946566 : DO j1 = m1 + 1, md1
528 57454692 : zf(j1, j3, jp2) = 0._dp
529 : END DO
530 : END DO
531 2598194 : DO j3 = m3 + 1, md3
532 52839068 : DO j1 = 1, md1
533 51347194 : zf(j1, j3, jp2) = 0._dp
534 : END DO
535 : END DO
536 : END DO
537 80206 : DO jp2 = nxc + 1, md2/nproc
538 489876 : DO j3 = 1, md3
539 16818520 : DO j1 = 1, md1
540 16806912 : zf(j1, j3, jp2) = 0._dp
541 : END DO
542 : END DO
543 : END DO
544 : ELSE
545 0 : zf = 0._dp
546 : END IF
547 68598 : CALL timestop(handle)
548 :
549 68598 : END SUBROUTINE scale_and_distribute
550 : END MODULE ps_wavelet_util
|
