Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 32597 : 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 32597 : 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 32597 : 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 15274 : ELSE IF (geocode == 'S') THEN
103 18 : CALL S_FFT_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
104 15256 : ELSE IF (geocode == 'F') THEN
105 15256 : 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 162985 : 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 32597 : istart = iproc*(md2/nproc)
128 32597 : iend = MIN((iproc + 1)*md2/nproc, m2)
129 :
130 32597 : nxc = iend - istart
131 32597 : nwbl = 0
132 32597 : nwbr = 0
133 32597 : nxcl = 1
134 32597 : nxcr = 1
135 :
136 32597 : nwb = nxcl + nxc + nxcr - 2
137 32597 : nxt = nwbr + nwb + nwbl
138 :
139 : !calculate the actual limit of the array for the zero padded FFT
140 32597 : IF (geocode == 'P') THEN
141 17323 : nlim = n2
142 15274 : ELSE IF (geocode == 'S') THEN
143 18 : nlim = n2
144 15256 : ELSE IF (geocode == 'F') THEN
145 15256 : 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 32597 : IF (istart + 1 <= m2) THEN
152 32597 : red_fact = 1._dp
153 32597 : 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 32597 : 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 15274 : ELSE IF (geocode == 'S') THEN
172 : !only one power of hgrid
173 18 : 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 18 : scal, pw_grid%para%group)
176 15256 : ELSE IF (geocode == 'F') THEN
177 15256 : hgrid = MAX(hx, hy, hz)
178 15256 : 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 15256 : scal, pw_grid%para%group)
181 15256 : 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 32597 : IF (geocode == 'F') THEN
188 15256 : red_fact = 1._dp
189 : ELSE
190 17341 : red_fact = -fourpi
191 : END IF
192 :
193 32597 : CALL scale_and_distribute(m1, m3, md1, md2, md3, nxc, zf, rhopot, nproc, red_fact)
194 :
195 32597 : DEALLOCATE (zf)
196 :
197 32597 : 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 : INTEGER :: l1, l2, l3
244 :
245 : !dimensions of the density in the real space
246 :
247 17649 : m1 = n01
248 17649 : m2 = n03
249 17649 : m3 = n02
250 :
251 : ! real space grid dimension (suitable for number of processors)
252 17649 : l1 = m1
253 17649 : l2 = m2
254 17649 : l3 = m3 !beware of the half dimension
255 17649 : CALL fourier_dim(l1, n1)
256 17649 : IF (n1 == m1) THEN
257 : ELSE
258 0 : PRINT *, 'the FFT in the x direction is not allowed'
259 0 : PRINT *, 'n01 dimension', n01
260 0 : CPABORT("")
261 : END IF
262 17649 : l1 = l1 + 1
263 17649 : CALL fourier_dim(l2, n2)
264 17649 : IF (n2 == m2) THEN
265 : ELSE
266 0 : PRINT *, 'the FFT in the z direction is not allowed'
267 0 : PRINT *, 'n03 dimension', n03
268 0 : CPABORT("")
269 : END IF
270 17649 : CALL fourier_dim(l3, n3)
271 17649 : IF (n3 == m3) THEN
272 : ELSE
273 0 : PRINT *, 'the FFT in the y direction is not allowed'
274 0 : PRINT *, 'n02 dimension', n02
275 0 : CPABORT("")
276 : END IF
277 :
278 : !dimensions that contain the unpadded real space,
279 : ! compatible with the number of processes
280 17649 : md1 = n1
281 17649 : md2 = n2
282 17649 : md3 = n3
283 19029 : DO WHILE (nproc*(md2/nproc) < n2)
284 1380 : md2 = md2 + 1
285 : END DO
286 :
287 : !dimensions of the kernel, 1/8 of the total volume,
288 : !compatible with nproc
289 17649 : nd1 = n1/2 + 1
290 17649 : nd2 = n2/2 + 1
291 17649 : nd3 = n3/2 + 1
292 19579 : DO WHILE (MODULO(nd3, nproc) /= 0)
293 1930 : nd3 = nd3 + 1
294 : END DO
295 :
296 17649 : END SUBROUTINE P_FFT_dimensions
297 :
298 : ! **************************************************************************************************
299 : !> \brief Calculate four sets of dimension needed for the calculation of the
300 : !> convolution for the surface system
301 : !> \param n01 original real dimensions (input)
302 : !> \param n02 original real dimensions (input)
303 : !> \param n03 original real dimensions (input)
304 : !> \param m1 original real dimension, with 2 and 3 exchanged
305 : !> \param m2 original real dimension, with 2 and 3 exchanged
306 : !> \param m3 original real dimension, with 2 and 3 exchanged
307 : !> \param n1 the first FFT dimensions, for the moment supposed to be even
308 : !> \param n2 the first FFT dimensions, for the moment supposed to be even
309 : !> \param n3 the double of the first FFT even dimension greater than m3
310 : !> (improved for the HalFFT procedure)
311 : !> \param md1 the n1,n2 dimensions.
312 : !> \param md2 the n1,n2,n3 dimensions.
313 : !> \param md3 the half of n3 dimension. They contain the real unpadded space,
314 : !> properly enlarged to be compatible with the FFT dimensions n_i.
315 : !> md2 is further enlarged to be a multiple of nproc
316 : !> \param nd1 fourier dimensions for which the kernel is injective,
317 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
318 : !> enlarged to be a multiple of nproc
319 : !> \param nd2 fourier dimensions for which the kernel is injective,
320 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
321 : !> enlarged to be a multiple of nproc
322 : !> \param nd3 fourier dimensions for which the kernel is injective,
323 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
324 : !> enlarged to be a multiple of nproc
325 : !> \param nproc ...
326 : !> \date October 2006
327 : !> \author Luigi Genovese
328 : !> \note This four sets of dimensions are actually redundant (mi=n0i),
329 : !> due to the backward-compatibility
330 : !> with the Poisson Solver with other geometries.
331 : !> Dimensions n02 and n03 were exchanged
332 : ! **************************************************************************************************
333 66 : SUBROUTINE S_FFT_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
334 : INTEGER, INTENT(in) :: n01, n02, n03
335 : INTEGER, INTENT(out) :: m1, m2, m3, n1, n2, n3, md1, md2, md3, &
336 : nd1, nd2, nd3
337 : INTEGER, INTENT(in) :: nproc
338 :
339 : CHARACTER(len=*), PARAMETER :: routineN = 'S_FFT_dimensions'
340 :
341 : INTEGER :: handle, l1, l2, l3
342 :
343 : !dimensions of the density in the real space
344 :
345 22 : CALL timeset(routineN, handle)
346 22 : m1 = n01
347 22 : m2 = n03
348 22 : m3 = n02
349 :
350 : ! real space grid dimension (suitable for number of processors)
351 22 : l1 = m1
352 22 : l2 = m2
353 22 : l3 = m3 !beware of the half dimension
354 22 : CALL fourier_dim(l1, n1)
355 22 : IF (n1 == m1) THEN
356 : ELSE
357 0 : PRINT *, 'the FFT in the x direction is not allowed'
358 0 : PRINT *, 'n01 dimension', n01
359 0 : CPABORT("")
360 : END IF
361 22 : l1 = l1 + 1
362 22 : CALL fourier_dim(l2, n2)
363 22 : IF (n2 == m2) THEN
364 : ELSE
365 0 : PRINT *, 'the FFT in the z direction is not allowed'
366 0 : PRINT *, 'n03 dimension', n03
367 0 : CPABORT("")
368 : END IF
369 22 : DO
370 22 : CALL fourier_dim(l3, n3)
371 22 : IF (MODULO(n3, 2) == 0) THEN
372 : EXIT
373 : END IF
374 0 : l3 = l3 + 1
375 : END DO
376 22 : n3 = 2*n3
377 :
378 : !dimensions that contain the unpadded real space,
379 : ! compatible with the number of processes
380 22 : md1 = n1
381 22 : md2 = n2
382 22 : md3 = n3/2
383 22 : DO WHILE (nproc*(md2/nproc) < n2)
384 0 : md2 = md2 + 1
385 : END DO
386 :
387 : !dimensions of the kernel, 1/8 of the total volume,
388 : !compatible with nproc
389 :
390 : !these two dimensions are like that since they are even
391 22 : nd1 = n1/2 + 1
392 22 : nd2 = n2/2 + 1
393 :
394 22 : nd3 = n3/2 + 1
395 44 : DO WHILE (MODULO(nd3, nproc) /= 0)
396 22 : nd3 = nd3 + 1
397 : END DO
398 22 : CALL timestop(handle)
399 :
400 22 : END SUBROUTINE S_FFT_dimensions
401 :
402 : ! **************************************************************************************************
403 : !> \brief Calculate four sets of dimension needed for the calculation of the
404 : !> zero-padded convolution
405 : !> \param n01 original real dimensions (input)
406 : !> \param n02 original real dimensions (input)
407 : !> \param n03 original real dimensions (input)
408 : !> \param m1 original real dimension with the dimension 2 and 3 exchanged
409 : !> \param m2 original real dimension with the dimension 2 and 3 exchanged
410 : !> \param m3 original real dimension with the dimension 2 and 3 exchanged
411 : !> \param n1 ...
412 : !> \param n2 ...
413 : !> \param n3 the double of the first FFT even dimension greater than m3
414 : !> (improved for the HalFFT procedure)
415 : !> \param md1 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 md2 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 md3 half of n1,n2,n3 dimension. They contain the real unpadded space,
422 : !> properly enlarged to be compatible with the FFT dimensions n_i.
423 : !> md2 is further enlarged to be a multiple of nproc
424 : !> \param nd1 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 nd2 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 nd3 fourier dimensions for which the kernel FFT is injective,
431 : !> formally 1/8 of the fourier grid. Here the dimension nd3 is
432 : !> enlarged to be a multiple of nproc
433 : !> \param nproc ...
434 : !> \date February 2006
435 : !> \author Luigi Genovese
436 : !> \note The dimension m2 and m3 correspond to n03 and n02 respectively
437 : !> this is needed since the convolution routine manage arrays of dimension
438 : !> (md1,md3,md2/nproc)
439 : ! **************************************************************************************************
440 16598 : SUBROUTINE F_FFT_dimensions(n01, n02, n03, m1, m2, m3, n1, n2, n3, md1, md2, md3, nd1, nd2, nd3, nproc)
441 : INTEGER, INTENT(in) :: n01, n02, n03
442 : INTEGER, INTENT(out) :: m1, m2, m3, n1, n2, n3, md1, md2, md3, &
443 : nd1, nd2, nd3
444 : INTEGER, INTENT(in) :: nproc
445 :
446 : INTEGER :: l1, l2, l3
447 :
448 : !dimensions of the density in the real space, inverted for convenience
449 :
450 16598 : m1 = n01
451 16598 : m2 = n03
452 16598 : m3 = n02
453 : ! real space grid dimension (suitable for number of processors)
454 16598 : l1 = 2*m1
455 16598 : l2 = 2*m2
456 16598 : l3 = m3 !beware of the half dimension
457 16598 : DO
458 16598 : CALL fourier_dim(l1, n1)
459 16598 : IF (MODULO(n1, 2) == 0) THEN
460 : EXIT
461 : END IF
462 16598 : l1 = l1 + 1
463 : END DO
464 16598 : DO
465 16598 : CALL fourier_dim(l2, n2)
466 16598 : IF (MODULO(n2, 2) == 0) THEN
467 : EXIT
468 : END IF
469 16598 : l2 = l2 + 1
470 : END DO
471 26194 : DO
472 26194 : CALL fourier_dim(l3, n3)
473 26194 : IF (MODULO(n3, 2) == 0) THEN
474 : EXIT
475 : END IF
476 9596 : l3 = l3 + 1
477 : END DO
478 16598 : n3 = 2*n3
479 :
480 : !dimensions that contain the unpadded real space,
481 : ! compatible with the number of processes
482 16598 : md1 = n1/2
483 16598 : md2 = n2/2
484 16598 : md3 = n3/2
485 18780 : DO WHILE (nproc*(md2/nproc) < n2/2)
486 2182 : md2 = md2 + 1
487 : END DO
488 :
489 : !dimensions of the kernel, 1/8 of the total volume,
490 : !compatible with nproc
491 16598 : nd1 = n1/2 + 1
492 16598 : nd2 = n2/2 + 1
493 16598 : nd3 = n3/2 + 1
494 :
495 25068 : DO WHILE (MODULO(nd3, nproc) /= 0)
496 8470 : nd3 = nd3 + 1
497 : END DO
498 :
499 16598 : END SUBROUTINE F_FFT_dimensions
500 :
501 : ! **************************************************************************************************
502 : !> \brief ...
503 : !> \param m1 ...
504 : !> \param m3 ...
505 : !> \param md1 ...
506 : !> \param md2 ...
507 : !> \param md3 ...
508 : !> \param nxc ...
509 : !> \param rhopot ...
510 : !> \param zf ...
511 : !> \param nproc ...
512 : !> \param factor ...
513 : ! **************************************************************************************************
514 65194 : SUBROUTINE scale_and_distribute(m1, m3, md1, md2, md3, nxc, &
515 65194 : rhopot, zf, nproc, factor)
516 :
517 : !Arguments----------------------
518 : INTEGER, INTENT(in) :: m1, m3, md1, md2, md3, nxc, nproc
519 : REAL(KIND=dp), DIMENSION(md1, md3, md2/nproc), &
520 : INTENT(inout) :: zf, rhopot
521 : REAL(KIND=dp), INTENT(in) :: factor
522 :
523 : CHARACTER(len=*), PARAMETER :: routineN = 'scale_and_distribute'
524 :
525 : INTEGER :: handle, j1, j3, jp2
526 :
527 65194 : CALL timeset(routineN, handle)
528 :
529 65194 : IF (nxc >= 1) THEN
530 1456540 : DO jp2 = 1, nxc
531 46555924 : DO j3 = 1, m3
532 1870036468 : DO j1 = 1, m1
533 1870036468 : zf(j1, j3, jp2) = factor*rhopot(j1, j3, jp2)
534 : END DO
535 51429674 : DO j1 = m1 + 1, md1
536 50038328 : zf(j1, j3, jp2) = 0._dp
537 : END DO
538 : END DO
539 2137238 : DO j3 = m3 + 1, md3
540 25953592 : DO j1 = 1, md1
541 24562246 : zf(j1, j3, jp2) = 0._dp
542 : END DO
543 : END DO
544 : END DO
545 75678 : DO jp2 = nxc + 1, md2/nproc
546 399136 : DO j3 = 1, md3
547 10096876 : DO j1 = 1, md1
548 10086392 : zf(j1, j3, jp2) = 0._dp
549 : END DO
550 : END DO
551 : END DO
552 : ELSE
553 0 : zf = 0._dp
554 : END IF
555 65194 : CALL timestop(handle)
556 :
557 65194 : END SUBROUTINE scale_and_distribute
558 : END MODULE ps_wavelet_util
|
