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 :
14 : USE kinds, ONLY: dp
15 : USE mathconstants, ONLY: fourpi
16 : USE ps_wavelet_base, ONLY: f_poissonsolver,&
17 : p_poissonsolver,&
18 : s_poissonsolver
19 : USE ps_wavelet_fft3d, ONLY: fourier_dim
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 33295 : 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 33295 : 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 33295 : 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 15972 : 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 15918 : ELSE IF (geocode == 'F') THEN
105 15918 : 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 166475 : ALLOCATE (zf(md1, md3, md2/nproc))
111 :
112 : ! CALL timing(iproc,'Exchangecorr ','ON')
113 : !dimension for exchange-correlation (different in the global or distributed case)
114 : !let us calculate the dimension of the portion of the rhopot array to be passed
115 : !to the xc routine
116 : !this portion will depend on the need of calculating the gradient or not,
117 : !and whether the White-Bird correction must be inserted or not
118 : !(absent only in the LB ixc=13 case)
119 :
120 : !nxc is the effective part of the third dimension that is being processed
121 : !nxt is the dimension of the part of rhopot that must be passed to the gradient routine
122 : !nwb is the dimension of the part of rhopot in the wb-postprocessing routine
123 : !note: nxc <= nwb <= nxt
124 : !the dimension are related by the values of nwbl and nwbr
125 : ! nxc+nxcl+nxcr-2 = nwb
126 : ! nwb+nwbl+nwbr = nxt
127 33295 : istart = iproc*(md2/nproc)
128 33295 : iend = MIN((iproc + 1)*md2/nproc, m2)
129 :
130 33295 : nxc = iend - istart
131 33295 : nwbl = 0
132 33295 : nwbr = 0
133 33295 : nxcl = 1
134 33295 : nxcr = 1
135 :
136 33295 : nwb = nxcl + nxc + nxcr - 2
137 33295 : nxt = nwbr + nwb + nwbl
138 :
139 : !calculate the actual limit of the array for the zero padded FFT
140 33295 : IF (geocode == 'P') THEN
141 17323 : nlim = n2
142 15972 : ELSE IF (geocode == 'S') THEN
143 54 : nlim = n2
144 15918 : ELSE IF (geocode == 'F') THEN
145 15918 : nlim = n2/2
146 : END IF
147 :
148 : !!$ print *,'density must go from',min(istart+1,m2),'to',iend,'with n2/2=',n2/2
149 : !!$ print *,' it goes from',i3start+nwbl+nxcl-1,'to',i3start+nxc-1
150 :
151 33295 : IF (istart + 1 <= m2) THEN
152 33295 : red_fact = 1._dp
153 33295 : CALL scale_and_distribute(m1, m3, md1, md2, md3, nxc, rhopot, zf, nproc, red_fact)
154 0 : ELSE IF (istart + 1 <= nlim) THEN !this condition assures that we have perform good zero padding
155 0 : DO i2 = istart + 1, MIN(nlim, istart + md2/nproc)
156 0 : j2 = i2 - istart
157 0 : DO i3 = 1, md3
158 0 : DO i1 = 1, md1
159 0 : zf(i1, i3, j2) = 0._dp
160 : END DO
161 : END DO
162 : END DO
163 : END IF
164 :
165 : !this routine builds the values for each process of the potential (zf), multiplying by scal
166 33295 : IF (geocode == 'P') THEN
167 : !no powers of hgrid because they are incorporated in the plane wave treatment
168 17323 : scal = 1._dp/REAL(n1*n2*n3, KIND=dp)
169 : CALL P_PoissonSolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, zf, &
170 17323 : scal, hx, hy, hz, pw_grid%para%group)
171 15972 : ELSE IF (geocode == 'S') THEN
172 : !only one power of hgrid
173 54 : scal = hy/REAL(n1*n2*n3, KIND=dp)
174 : CALL S_PoissonSolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, karray, zf, &
175 54 : scal, pw_grid%para%group)
176 15918 : ELSE IF (geocode == 'F') THEN
177 15918 : hgrid = MAX(hx, hy, hz)
178 15918 : scal = hgrid**3/REAL(n1*n2*n3, KIND=dp)
179 : CALL F_PoissonSolver(n1, n2, n3, nd1, nd2, nd3, md1, md2, md3, nproc, iproc, karray, zf, &
180 15918 : scal, pw_grid%para%group)
181 15918 : factor = 0.5_dp*hgrid**3
182 : END IF
183 :
184 : ! call timing(iproc,'PSolv_comput ','ON')
185 :
186 : !the value of the shift depends on the distributed i/o or not
187 33295 : IF (geocode == 'F') THEN
188 15918 : red_fact = 1._dp
189 : ELSE
190 17377 : red_fact = -fourpi
191 : END IF
192 :
193 33295 : CALL scale_and_distribute(m1, m3, md1, md2, md3, nxc, zf, rhopot, nproc, red_fact)
194 :
195 33295 : DEALLOCATE (zf)
196 :
197 33295 : END SUBROUTINE PSolver
198 :
199 : ! **************************************************************************************************
200 : !> \brief Calculate four sets of dimension needed for the calculation of the
201 : !> convolution for the periodic system
202 : !> \param n01 original real dimensions (input)
203 : !> \param n02 original real dimensions (input)
204 : !> \param n03 original real dimensions (input)
205 : !> \param m1 original real dimension, with m2 and m3 exchanged
206 : !> \param m2 original real dimension, with m2 and m3 exchanged
207 : !> \param m3 original real dimension, with m2 and m3 exchanged
208 : !> \param n1 the first FFT dimensions, for the moment supposed to be even
209 : !> \param n2 the first FFT dimensions, for the moment supposed to be even
210 : !> \param n3 the first FFT dimensions, for the moment supposed to be even
211 : !> \param md1 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 md2 the n1,n2,n3 dimensions. They contain the real unpadded space,
215 : !> properly enlarged to be compatible with the FFT dimensions n_i.
216 : !> md2 is further enlarged to be a multiple of nproc
217 : !> \param md3 the n1,n2,n3 dimensions. They contain the real unpadded space,
218 : !> properly enlarged to be compatible with the FFT dimensions n_i.
219 : !> md2 is further enlarged to be a multiple of nproc
220 : !> \param nd1 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 nd2 fourier dimensions for which the kernel is injective,
224 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
225 : !> enlarged to be a multiple of nproc
226 : !> \param nd3 fourier dimensions for which the kernel is injective,
227 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
228 : !> enlarged to be a multiple of nproc
229 : !> \param nproc ...
230 : !> \date October 2006
231 : !> \author Luigi Genovese
232 : !> \note This four sets of dimensions are actually redundant (mi=n0i),
233 : !> due to the backward-compatibility
234 : !> with the other geometries of the Poisson Solver.
235 : !> The dimensions 2 and 3 are exchanged.
236 : ! **************************************************************************************************
237 52947 : SUBROUTINE P_FFT_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
238 : INTEGER, INTENT(in) :: n01, n02, n03
239 : INTEGER, INTENT(out) :: m1, m2, m3, n1, n2, n3, md1, md2, md3, &
240 : nd1, nd2, nd3
241 : INTEGER, INTENT(in) :: nproc
242 :
243 : CHARACTER(len=80) :: err
244 : INTEGER :: l1, l2, l3
245 :
246 : !dimensions of the density in the real space
247 :
248 17649 : m1 = n01
249 17649 : m2 = n03
250 17649 : m3 = n02
251 :
252 : ! real space grid dimension (suitable for number of processors)
253 17649 : l1 = m1
254 17649 : l2 = m2
255 17649 : l3 = m3 !beware of the half dimension
256 17649 : CALL fourier_dim(l1, n1)
257 17649 : IF (n1 == m1) THEN
258 : ELSE
259 0 : WRITE (err, *) 'the FFT in the x direction is not allowed; n01 dimension ', n01
260 0 : CPABORT(TRIM(err))
261 : END IF
262 : l1 = l1 + 1
263 17649 : CALL fourier_dim(l2, n2)
264 17649 : IF (n2 == m2) THEN
265 : ELSE
266 0 : WRITE (err, *) 'the FFT in the z direction is not allowed; n03 dimension ', n03
267 0 : CPABORT(TRIM(err))
268 : END IF
269 17649 : CALL fourier_dim(l3, n3)
270 17649 : IF (n3 == m3) THEN
271 : ELSE
272 0 : WRITE (err, *) 'the FFT in the y direction is not allowed; n02 dimension ', n02
273 0 : CPABORT(TRIM(err))
274 : END IF
275 :
276 : !dimensions that contain the unpadded real space,
277 : ! compatible with the number of processes
278 17649 : md1 = n1
279 17649 : md2 = n2
280 17649 : md3 = n3
281 19029 : DO WHILE (nproc*(md2/nproc) < n2)
282 1380 : md2 = md2 + 1
283 : END DO
284 :
285 : !dimensions of the kernel, 1/8 of the total volume,
286 : !compatible with nproc
287 17649 : nd1 = n1/2 + 1
288 17649 : nd2 = n2/2 + 1
289 17649 : nd3 = n3/2 + 1
290 19579 : DO WHILE (MODULO(nd3, nproc) /= 0)
291 1930 : nd3 = nd3 + 1
292 : END DO
293 :
294 17649 : END SUBROUTINE P_FFT_dimensions
295 :
296 : ! **************************************************************************************************
297 : !> \brief Calculate four sets of dimension needed for the calculation of the
298 : !> convolution for the surface system
299 : !> \param n01 original real dimensions (input)
300 : !> \param n02 original real dimensions (input)
301 : !> \param n03 original real dimensions (input)
302 : !> \param m1 original real dimension, with 2 and 3 exchanged
303 : !> \param m2 original real dimension, with 2 and 3 exchanged
304 : !> \param m3 original real dimension, with 2 and 3 exchanged
305 : !> \param n1 the first FFT dimensions, for the moment supposed to be even
306 : !> \param n2 the first FFT dimensions, for the moment supposed to be even
307 : !> \param n3 the double of the first FFT even dimension greater than m3
308 : !> (improved for the HalFFT procedure)
309 : !> \param md1 the n1,n2 dimensions.
310 : !> \param md2 the n1,n2,n3 dimensions.
311 : !> \param md3 the half of n3 dimension. They contain the real unpadded space,
312 : !> properly enlarged to be compatible with the FFT dimensions n_i.
313 : !> md2 is further enlarged to be a multiple of nproc
314 : !> \param nd1 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 nd2 fourier dimensions for which the kernel is injective,
318 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
319 : !> enlarged to be a multiple of nproc
320 : !> \param nd3 fourier dimensions for which the kernel is injective,
321 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
322 : !> enlarged to be a multiple of nproc
323 : !> \param nproc ...
324 : !> \date October 2006
325 : !> \author Luigi Genovese
326 : !> \note This four sets of dimensions are actually redundant (mi=n0i),
327 : !> due to the backward-compatibility
328 : !> with the Poisson Solver with other geometries.
329 : !> Dimensions n02 and n03 were exchanged
330 : ! **************************************************************************************************
331 198 : SUBROUTINE S_FFT_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
332 : INTEGER, INTENT(in) :: n01, n02, n03
333 : INTEGER, INTENT(out) :: m1, m2, m3, n1, n2, n3, md1, md2, md3, &
334 : nd1, nd2, nd3
335 : INTEGER, INTENT(in) :: nproc
336 :
337 : CHARACTER(len=*), PARAMETER :: routineN = 'S_FFT_dimensions'
338 :
339 : CHARACTER(len=80) :: err
340 : INTEGER :: handle, l1, l2, l3
341 :
342 : !dimensions of the density in the real space
343 :
344 66 : CALL timeset(routineN, handle)
345 66 : m1 = n01
346 66 : m2 = n03
347 66 : m3 = n02
348 :
349 : ! real space grid dimension (suitable for number of processors)
350 66 : l1 = m1
351 66 : l2 = m2
352 66 : l3 = m3 !beware of the half dimension
353 66 : CALL fourier_dim(l1, n1)
354 66 : IF (n1 == m1) THEN
355 : ELSE
356 0 : WRITE (err, *) 'the FFT in the x direction is not allowed; n01 dimension', n01
357 0 : CPABORT(TRIM(err))
358 : END IF
359 : l1 = l1 + 1
360 66 : CALL fourier_dim(l2, n2)
361 66 : IF (n2 == m2) THEN
362 : ELSE
363 0 : WRITE (err, *) 'the FFT in the z direction is not allowed; n03 dimension', n03
364 0 : CPABORT(TRIM(err))
365 : END IF
366 66 : DO
367 66 : CALL fourier_dim(l3, n3)
368 66 : IF (MODULO(n3, 2) == 0) THEN
369 : EXIT
370 : END IF
371 0 : l3 = l3 + 1
372 : END DO
373 66 : n3 = 2*n3
374 :
375 : !dimensions that contain the unpadded real space,
376 : ! compatible with the number of processes
377 66 : md1 = n1
378 66 : md2 = n2
379 66 : md3 = n3/2
380 66 : DO WHILE (nproc*(md2/nproc) < n2)
381 0 : md2 = md2 + 1
382 : END DO
383 :
384 : !dimensions of the kernel, 1/8 of the total volume,
385 : !compatible with nproc
386 :
387 : !these two dimensions are like that since they are even
388 66 : nd1 = n1/2 + 1
389 66 : nd2 = n2/2 + 1
390 :
391 66 : nd3 = n3/2 + 1
392 132 : DO WHILE (MODULO(nd3, nproc) /= 0)
393 66 : nd3 = nd3 + 1
394 : END DO
395 66 : CALL timestop(handle)
396 :
397 66 : END SUBROUTINE S_FFT_dimensions
398 :
399 : ! **************************************************************************************************
400 : !> \brief Calculate four sets of dimension needed for the calculation of the
401 : !> zero-padded convolution
402 : !> \param n01 original real dimensions (input)
403 : !> \param n02 original real dimensions (input)
404 : !> \param n03 original real dimensions (input)
405 : !> \param m1 original real dimension with the dimension 2 and 3 exchanged
406 : !> \param m2 original real dimension with the dimension 2 and 3 exchanged
407 : !> \param m3 original real dimension with the dimension 2 and 3 exchanged
408 : !> \param n1 ...
409 : !> \param n2 ...
410 : !> \param n3 the double of the first FFT even dimension greater than m3
411 : !> (improved for the HalFFT procedure)
412 : !> \param md1 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 md2 half of n1,n2,n3 dimension. They contain the real unpadded space,
416 : !> properly enlarged to be compatible with the FFT dimensions n_i.
417 : !> md2 is further enlarged to be a multiple of nproc
418 : !> \param md3 half of n1,n2,n3 dimension. They contain the real unpadded space,
419 : !> properly enlarged to be compatible with the FFT dimensions n_i.
420 : !> md2 is further enlarged to be a multiple of nproc
421 : !> \param nd1 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 nd2 fourier dimensions for which the kernel FFT is injective,
425 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
426 : !> enlarged to be a multiple of nproc
427 : !> \param nd3 fourier dimensions for which the kernel FFT is injective,
428 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
429 : !> enlarged to be a multiple of nproc
430 : !> \param nproc ...
431 : !> \date February 2006
432 : !> \author Luigi Genovese
433 : !> \note The dimension m2 and m3 correspond to n03 and n02 respectively
434 : !> this is needed since the convolution routine manage arrays of dimension
435 : !> (md1,md3,md2/nproc)
436 : ! **************************************************************************************************
437 17292 : SUBROUTINE F_FFT_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
438 : INTEGER, INTENT(in) :: n01, n02, n03
439 : INTEGER, INTENT(out) :: m1, m2, m3, n1, n2, n3, md1, md2, md3, &
440 : nd1, nd2, nd3
441 : INTEGER, INTENT(in) :: nproc
442 :
443 : INTEGER :: l1, l2, l3
444 :
445 : !dimensions of the density in the real space, inverted for convenience
446 :
447 17292 : m1 = n01
448 17292 : m2 = n03
449 17292 : m3 = n02
450 : ! real space grid dimension (suitable for number of processors)
451 17292 : l1 = 2*m1
452 17292 : l2 = 2*m2
453 17292 : l3 = m3 !beware of the half dimension
454 17292 : DO
455 17292 : CALL fourier_dim(l1, n1)
456 17292 : IF (MODULO(n1, 2) == 0) THEN
457 : EXIT
458 : END IF
459 17292 : l1 = l1 + 1
460 : END DO
461 17292 : DO
462 17292 : CALL fourier_dim(l2, n2)
463 17292 : IF (MODULO(n2, 2) == 0) THEN
464 : EXIT
465 : END IF
466 17292 : l2 = l2 + 1
467 : END DO
468 27082 : DO
469 27082 : CALL fourier_dim(l3, n3)
470 27082 : IF (MODULO(n3, 2) == 0) THEN
471 : EXIT
472 : END IF
473 9790 : l3 = l3 + 1
474 : END DO
475 17292 : n3 = 2*n3
476 :
477 : !dimensions that contain the unpadded real space,
478 : ! compatible with the number of processes
479 17292 : md1 = n1/2
480 17292 : md2 = n2/2
481 17292 : md3 = n3/2
482 19640 : DO WHILE (nproc*(md2/nproc) < n2/2)
483 2348 : md2 = md2 + 1
484 : END DO
485 :
486 : !dimensions of the kernel, 1/8 of the total volume,
487 : !compatible with nproc
488 17292 : nd1 = n1/2 + 1
489 17292 : nd2 = n2/2 + 1
490 17292 : nd3 = n3/2 + 1
491 :
492 26456 : DO WHILE (MODULO(nd3, nproc) /= 0)
493 9164 : nd3 = nd3 + 1
494 : END DO
495 :
496 17292 : END SUBROUTINE F_FFT_dimensions
497 :
498 : ! **************************************************************************************************
499 : !> \brief ...
500 : !> \param m1 ...
501 : !> \param m3 ...
502 : !> \param md1 ...
503 : !> \param md2 ...
504 : !> \param md3 ...
505 : !> \param nxc ...
506 : !> \param rhopot ...
507 : !> \param zf ...
508 : !> \param nproc ...
509 : !> \param factor ...
510 : ! **************************************************************************************************
511 66590 : SUBROUTINE scale_and_distribute(m1, m3, md1, md2, md3, nxc, &
512 66590 : rhopot, zf, nproc, factor)
513 :
514 : !Arguments----------------------
515 : INTEGER, INTENT(in) :: m1, m3, md1, md2, md3, nxc, nproc
516 : REAL(KIND=dp), DIMENSION(md1, md3, md2/nproc), &
517 : INTENT(inout) :: zf, rhopot
518 : REAL(KIND=dp), INTENT(in) :: factor
519 :
520 : CHARACTER(len=*), PARAMETER :: routineN = 'scale_and_distribute'
521 :
522 : INTEGER :: handle, j1, j3, jp2
523 :
524 66590 : CALL timeset(routineN, handle)
525 :
526 66590 : IF (nxc >= 1) THEN
527 1489542 : DO jp2 = 1, nxc
528 48002796 : DO j3 = 1, m3
529 1919562780 : DO j1 = 1, m1
530 1919562780 : zf(j1, j3, jp2) = factor*rhopot(j1, j3, jp2)
531 : END DO
532 54001546 : DO j1 = m1 + 1, md1
533 52578594 : zf(j1, j3, jp2) = 0._dp
534 : END DO
535 : END DO
536 2220754 : DO j3 = m3 + 1, md3
537 29985686 : DO j1 = 1, md1
538 28562734 : zf(j1, j3, jp2) = 0._dp
539 : END DO
540 : END DO
541 : END DO
542 77300 : DO jp2 = nxc + 1, md2/nproc
543 419162 : DO j3 = 1, md3
544 11845654 : DO j1 = 1, md1
545 11834944 : zf(j1, j3, jp2) = 0._dp
546 : END DO
547 : END DO
548 : END DO
549 : ELSE
550 0 : zf = 0._dp
551 : END IF
552 66590 : CALL timestop(handle)
553 :
554 66590 : END SUBROUTINE scale_and_distribute
555 : END MODULE ps_wavelet_util
|
