(git:b279b6b)
pw_spline_utils.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief different utils that are useful to manipulate splines on the regular
10 !> grid of a pw
11 !> \par History
12 !> 05.2003 created [fawzi]
13 !> 08.2004 removed spline evaluation method using more than 2 read streams
14 !> (pw_compose_stripe_rs3), added linear solver based spline
15 !> inversion [fawzi]
16 !> \author Fawzi Mohamed
17 ! **************************************************************************************************
18 MODULE pw_spline_utils
19 
20  USE cp_log_handling, ONLY: cp_get_default_logger,&
21  cp_logger_type,&
22  cp_to_string
23  USE kinds, ONLY: dp
24  USE mathconstants, ONLY: twopi
25  USE message_passing, ONLY: mp_comm_congruent
26  USE pw_grid_types, ONLY: fullspace,&
27  pw_mode_local
28  USE pw_methods, ONLY: pw_axpy,&
29  pw_copy,&
30  pw_integral_ab,&
31  pw_zero
32  USE pw_pool_types, ONLY: pw_pool_release,&
33  pw_pool_type
34  USE pw_types, ONLY: pw_c1d_gs_type,&
35  pw_r3d_rs_type
36 #include "../base/base_uses.f90"
37 
38  IMPLICIT NONE
39  PRIVATE
40 
41  LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
42  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_spline_utils'
43 
44  INTEGER, PARAMETER, PUBLIC :: no_precond = 0, &
45  precond_spl3_aint = 1, &
46  precond_spl3_1 = 2, &
47  precond_spl3_aint2 = 3, &
48  precond_spl3_2 = 4, &
49  precond_spl3_3 = 5
50 
51  REAL(kind=dp), PUBLIC, PARAMETER, DIMENSION(4) :: nn10_coeffs = &
52  (/125._dp/216._dp, 25._dp/432._dp, 5._dp/864._dp, 1._dp/1728._dp/), &
53  spline3_coeffs = &
54  (/8._dp/(27._dp), 2._dp/(27._dp), 1._dp/(27._dp*2._dp), &
55  1._dp/(27._dp*8._dp)/), &
56  spline2_coeffs = &
57  (/27._dp/(64._dp), 9._dp/(64._dp*2_dp), 3._dp/(64._dp*4._dp), &
58  1._dp/(64._dp*8._dp)/), &
59  nn50_coeffs = &
60  (/15625._dp/17576._dp, 625._dp/35152._dp, 25._dp/70304._dp, &
61  1._dp/140608._dp/), &
62  spl3_aint_coeff = &
63  (/46._dp/27._dp, -2._dp/(27._dp), -1._dp/(27._dp*2._dp), &
64  -1._dp/(27._dp*8._dp)/), &
65  spl3_precond1_coeff = &
66  (/64._dp/3._dp, -8._dp/3._dp, -1._dp/3._dp, -1._dp/24._dp/), &
67  spl3_1d_transf_coeffs = &
68  (/2._dp/3._dp, 23._dp/48._dp, 1._dp/6._dp, 1._dp/48._dp/)
69 
70  REAL(kind=dp), PUBLIC, PARAMETER, DIMENSION(3) :: spline3_deriv_coeffs = &
71  (/2.0_dp/9.0_dp, 1.0_dp/18.0_dp, 1.0_dp/72.0_dp/), &
72  spline2_deriv_coeffs = &
73  (/9.0_dp/32.0_dp, 3.0_dp/64.0_dp, 1.0_dp/128.0_dp/), &
74  nn10_deriv_coeffs = &
75  (/25._dp/72._dp, 5._dp/144, 1._dp/288._dp/), &
76  nn50_deriv_coeffs = &
77  (/625._dp/1352._dp, 25._dp/2704._dp, 1._dp/5408._dp/), &
78  spl3_1d_coeffs0 = &
79  (/1._dp/6_dp, 2._dp/3._dp, 1._dp/6._dp/), &
80  spl3_1d_transf_border1 = &
81  (/0.517977704_dp, 0.464044595_dp, 0.17977701e-1_dp/)
82 
83  PUBLIC :: pw_spline3_interpolate_values_g, &
84  pw_spline3_deriv_g
85  PUBLIC :: pw_spline_scale_deriv
86  PUBLIC :: pw_spline2_interpolate_values_g, &
87  pw_spline2_deriv_g
88  PUBLIC :: pw_nn_smear_r, pw_nn_deriv_r, &
89  spl3_nopbc, spl3_pbc, spl3_nopbct
90  PUBLIC :: add_fine2coarse, add_coarse2fine
91  PUBLIC :: pw_spline_precond_create, &
92  pw_spline_do_precond, &
93  pw_spline_precond_set_kind, &
94  find_coeffs, &
95  pw_spline_precond_release, &
96  pw_spline_precond_type, &
97  eval_interp_spl3_pbc, &
98  eval_d_interp_spl3_pbc
99 
100 !***
101 
102 ! **************************************************************************************************
103 !> \brief stores information for the preconditioner used to calculate the
104 !> coeffs of splines
105 !> \author fawzi
106 ! **************************************************************************************************
107  TYPE pw_spline_precond_type
108  INTEGER :: kind = no_precond
109  REAL(kind=dp), DIMENSION(4) :: coeffs = 0.0_dp
110  REAL(kind=dp), DIMENSION(3) :: coeffs_1d = 0.0_dp
111  LOGICAL :: sharpen = .false., normalize = .false., pbc = .false., transpose = .false.
112  TYPE(pw_pool_type), POINTER :: pool => null()
113  END TYPE pw_spline_precond_type
114 
115 CONTAINS
116 
117 ! **************************************************************************************************
118 !> \brief calculates the FFT of the coefficients of the quadratic spline that
119 !> interpolates the given values
120 !> \param spline_g on entry the FFT of the values to interpolate as cc,
121 !> will contain the FFT of the coefficients of the spline
122 !> \par History
123 !> 06.2003 created [fawzi]
124 !> \author Fawzi Mohamed
125 !> \note
126 !> does not work with spherical cutoff
127 ! **************************************************************************************************
128  SUBROUTINE pw_spline2_interpolate_values_g(spline_g)
129  TYPE(pw_c1d_gs_type), INTENT(IN) :: spline_g
130 
131  CHARACTER(len=*), PARAMETER :: routinen = 'pw_spline2_interpolate_values_g'
132 
133  INTEGER :: handle, i, ii, j, k
134  INTEGER, DIMENSION(2, 3) :: gbo
135  INTEGER, DIMENSION(3) :: n_tot
136  REAL(kind=dp) :: c23, coeff
137  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cosivals, cosjvals, coskvals
138 
139  CALL timeset(routinen, handle)
140 
141  n_tot(1:3) = spline_g%pw_grid%npts(1:3)
142  gbo = spline_g%pw_grid%bounds
143 
144  cpassert(.NOT. spline_g%pw_grid%spherical)
145  cpassert(spline_g%pw_grid%grid_span == fullspace)
146 
147  ALLOCATE (cosivals(gbo(1, 1):gbo(2, 1)), cosjvals(gbo(1, 2):gbo(2, 2)), &
148  coskvals(gbo(1, 3):gbo(2, 3)))
149 
150  coeff = twopi/n_tot(1)
151 !$OMP PARALLEL DO PRIVATE(i) DEFAULT(NONE) SHARED(cosIVals,coeff,gbo)
152  DO i = gbo(1, 1), gbo(2, 1)
153  cosivals(i) = cos(coeff*real(i, dp))
154  END DO
155  coeff = twopi/n_tot(2)
156 !$OMP PARALLEL DO PRIVATE(j) DEFAULT(NONE) SHARED(cosJVals,coeff,gbo)
157  DO j = gbo(1, 2), gbo(2, 2)
158  cosjvals(j) = cos(coeff*real(j, dp))
159  END DO
160  coeff = twopi/n_tot(3)
161 !$OMP PARALLEL DO PRIVATE(k) DEFAULT(NONE) SHARED(cosKVals,coeff,gbo)
162  DO k = gbo(1, 3), gbo(2, 3)
163  coskvals(k) = cos(coeff*real(k, dp))
164  END DO
165 
166 !$OMP PARALLEL DO PRIVATE(i,j,k,ii,coeff,c23) DEFAULT(NONE) SHARED(spline_g,cosIVals,cosJVals,cosKVals)
167  DO ii = 1, SIZE(spline_g%array)
168  i = spline_g%pw_grid%g_hat(1, ii)
169  j = spline_g%pw_grid%g_hat(2, ii)
170  k = spline_g%pw_grid%g_hat(3, ii)
171 
172  c23 = cosjvals(j)*coskvals(k)
173  coeff = 64.0_dp/(cosivals(i)*c23 + &
174  (cosivals(i)*cosjvals(j) + cosivals(i)*coskvals(k) + c23)*3.0_dp + &
175  (cosivals(i) + cosjvals(j) + coskvals(k))*9.0_dp + &
176  27.0_dp)
177 
178  spline_g%array(ii) = spline_g%array(ii)*coeff
179 
180  END DO
181  DEALLOCATE (cosivals, cosjvals, coskvals)
182 
183  CALL timestop(handle)
184  END SUBROUTINE pw_spline2_interpolate_values_g
185 
186 ! **************************************************************************************************
187 !> \brief calculates the FFT of the coefficients of the2 cubic spline that
188 !> interpolates the given values
189 !> \param spline_g on entry the FFT of the values to interpolate as cc,
190 !> will contain the FFT of the coefficients of the spline
191 !> \par History
192 !> 06.2003 created [fawzi]
193 !> \author Fawzi Mohamed
194 !> \note
195 !> does not work with spherical cutoff
196 !> stupid distribution for cos calculation, it should calculate only the
197 !> needed cos, and avoid the mpi_allreduce
198 ! **************************************************************************************************
199  SUBROUTINE pw_spline3_interpolate_values_g(spline_g)
200  TYPE(pw_c1d_gs_type), INTENT(IN) :: spline_g
201 
202  CHARACTER(len=*), PARAMETER :: routinen = 'pw_spline3_interpolate_values_g'
203 
204  INTEGER :: handle, i, ii, j, k
205  INTEGER, DIMENSION(2, 3) :: gbo
206  INTEGER, DIMENSION(3) :: n_tot
207  REAL(kind=dp) :: c23, coeff
208  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cosivals, cosjvals, coskvals
209 
210  CALL timeset(routinen, handle)
211 
212  n_tot(1:3) = spline_g%pw_grid%npts(1:3)
213  gbo = spline_g%pw_grid%bounds
214 
215  cpassert(.NOT. spline_g%pw_grid%spherical)
216  cpassert(spline_g%pw_grid%grid_span == fullspace)
217 
218  ALLOCATE (cosivals(gbo(1, 1):gbo(2, 1)), &
219  cosjvals(gbo(1, 2):gbo(2, 2)), &
220  coskvals(gbo(1, 3):gbo(2, 3)))
221 
222  coeff = twopi/n_tot(1)
223 !$OMP PARALLEL DO PRIVATE(i) DEFAULT(NONE) SHARED(cosIVals,coeff,gbo)
224  DO i = gbo(1, 1), gbo(2, 1)
225  cosivals(i) = cos(coeff*real(i, dp))
226  END DO
227  coeff = twopi/n_tot(2)
228 !$OMP PARALLEL DO PRIVATE(j) DEFAULT(NONE) SHARED(cosJVals,coeff,gbo)
229  DO j = gbo(1, 2), gbo(2, 2)
230  cosjvals(j) = cos(coeff*real(j, dp))
231  END DO
232  coeff = twopi/n_tot(3)
233 !$OMP PARALLEL DO PRIVATE(k) DEFAULT(NONE) SHARED(cosKVals,coeff,gbo)
234  DO k = gbo(1, 3), gbo(2, 3)
235  coskvals(k) = cos(coeff*real(k, dp))
236  END DO
237 
238 !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i,j,k,ii,coeff,c23) SHARED(spline_g,cosIVals,cosJVals,cosKVals)
239  DO ii = 1, SIZE(spline_g%array)
240  i = spline_g%pw_grid%g_hat(1, ii)
241  j = spline_g%pw_grid%g_hat(2, ii)
242  k = spline_g%pw_grid%g_hat(3, ii)
243  ! no opt
244 !FM coeff=1.0/((cosVal(1)*cosVal(2)*cosVal(3))/27.0_dp+&
245 !FM (cosVal(1)*cosVal(2)+cosVal(1)*cosVal(3)+&
246 !FM cosVal(2)*cosVal(3))*2.0_dp/27.0_dp+&
247 !FM (cosVal(1)+cosVal(2)+cosVal(3))*4.0_dp/27.0_dp+&
248 !FM 8.0_dp/27.0_dp)
249  ! opt
250  c23 = cosjvals(j)*coskvals(k)
251  coeff = 27.0_dp/(cosivals(i)*c23 + &
252  (cosivals(i)*cosjvals(j) + cosivals(i)*coskvals(k) + c23)*2.0_dp + &
253  (cosivals(i) + cosjvals(j) + coskvals(k))*4.0_dp + &
254  8.0_dp)
255 
256  spline_g%array(ii) = spline_g%array(ii)*coeff
257 
258  END DO
259  DEALLOCATE (cosivals, cosjvals, coskvals)
260 
261  CALL timestop(handle)
262  END SUBROUTINE pw_spline3_interpolate_values_g
263 
264 ! **************************************************************************************************
265 !> \brief rescales the derivatives from gridspacing=1 to the real derivatives
266 !> \param deriv_vals_r an array of x,y,z derivatives
267 !> \param transpose if true applies the transpose of the map (defaults to
268 !> false)
269 !> \param scale a scaling factor (defaults to 1.0)
270 !> \par History
271 !> 06.2003 created [fawzi]
272 !> \author Fawzi Mohamed
273 ! **************************************************************************************************
274  SUBROUTINE pw_spline_scale_deriv(deriv_vals_r, transpose, scale)
275  TYPE(pw_r3d_rs_type), DIMENSION(3), INTENT(IN) :: deriv_vals_r
276  LOGICAL, INTENT(in), OPTIONAL :: transpose
277  REAL(kind=dp), INTENT(in), OPTIONAL :: scale
278 
279  CHARACTER(len=*), PARAMETER :: routinen = 'pw_spline_scale_deriv'
280 
281  INTEGER :: handle, i, idir, j, k
282  INTEGER, DIMENSION(2, 3) :: bo
283  INTEGER, DIMENSION(3) :: n_tot
284  LOGICAL :: diag, my_transpose
285  REAL(kind=dp) :: dval1, dval2, dval3, my_scale, scalef
286  REAL(kind=dp), DIMENSION(3, 3) :: dh_inv, h_grid
287  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: ddata, ddata2, ddata3
288 
289  CALL timeset(routinen, handle)
290 
291  my_transpose = .false.
292  IF (PRESENT(transpose)) my_transpose = transpose
293  my_scale = 1.0_dp
294  IF (PRESENT(scale)) my_scale = scale
295  n_tot(1:3) = deriv_vals_r(1)%pw_grid%npts(1:3)
296  bo = deriv_vals_r(1)%pw_grid%bounds_local
297  dh_inv = deriv_vals_r(1)%pw_grid%dh_inv
298 
299  ! map grid to real derivative
300  diag = .true.
301  IF (my_transpose) THEN
302  DO j = 1, 3
303  DO i = 1, 3
304  h_grid(j, i) = my_scale*dh_inv(i, j) ! REAL(n_tot(i),dp)*cell_h_inv(i,j)
305  IF (i /= j .AND. h_grid(j, i) /= 0.0_dp) diag = .false.
306  END DO
307  END DO
308  ELSE
309  DO j = 1, 3
310  DO i = 1, 3
311  h_grid(i, j) = my_scale*dh_inv(i, j) ! REAL(n_tot(i),dp)*cell_h_inv(i,j)
312  IF (i /= j .AND. h_grid(i, j) /= 0.0_dp) diag = .false.
313  END DO
314  END DO
315  END IF
316 
317  IF (diag) THEN
318  DO idir = 1, 3
319  ddata => deriv_vals_r(idir)%array
320  scalef = h_grid(idir, idir)
321  CALL dscal((bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1), &
322  scalef, ddata, 1)
323  END DO
324  ELSE
325  ddata => deriv_vals_r(1)%array
326  ddata2 => deriv_vals_r(2)%array
327  ddata3 => deriv_vals_r(3)%array
328 !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(k,j,i,dVal1,dVal2,dVal3) &
329 !$OMP SHARED(ddata,ddata2,ddata3,h_grid,bo)
330  DO k = bo(1, 3), bo(2, 3)
331  DO j = bo(1, 2), bo(2, 2)
332  DO i = bo(1, 1), bo(2, 1)
333 
334  dval1 = ddata(i, j, k)
335  dval2 = ddata2(i, j, k)
336  dval3 = ddata3(i, j, k)
337 
338  ddata(i, j, k) = h_grid(1, 1)*dval1 + &
339  h_grid(2, 1)*dval2 + h_grid(3, 1)*dval3
340  ddata2(i, j, k) = h_grid(1, 2)*dval1 + &
341  h_grid(2, 2)*dval2 + h_grid(3, 2)*dval3
342  ddata3(i, j, k) = h_grid(1, 3)*dval1 + &
343  h_grid(2, 3)*dval2 + h_grid(3, 3)*dval3
344 
345  END DO
346  END DO
347  END DO
348  END IF
349 
350  CALL timestop(handle)
351  END SUBROUTINE pw_spline_scale_deriv
352 
353 ! **************************************************************************************************
354 !> \brief calculates the FFT of the values of the x,y,z (idir=1,2,3)
355 !> derivative of the cubic spline
356 !> \param spline_g on entry the FFT of the coefficients of the spline
357 !> will contain the FFT of the derivative
358 !> \param idir direction of the derivative
359 !> \par History
360 !> 06.2003 created [fawzi]
361 !> \author Fawzi Mohamed
362 !> \note
363 !> the distance between gridpoints is assumed to be 1
364 ! **************************************************************************************************
365  SUBROUTINE pw_spline3_deriv_g(spline_g, idir)
366  TYPE(pw_c1d_gs_type), INTENT(IN) :: spline_g
367  INTEGER, INTENT(in) :: idir
368 
369  CHARACTER(len=*), PARAMETER :: routinen = 'pw_spline3_deriv_g'
370  REAL(kind=dp), PARAMETER :: inv9 = 1.0_dp/9.0_dp
371 
372  INTEGER :: handle, i, ii, j, k
373  INTEGER, DIMENSION(2, 3) :: bo, gbo
374  INTEGER, DIMENSION(3) :: n, n_tot
375  REAL(kind=dp) :: coeff, tmp
376  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: csivals, csjvals, cskvals
377 
378  CALL timeset(routinen, handle)
379 
380  n(1:3) = spline_g%pw_grid%npts_local(1:3)
381  n_tot(1:3) = spline_g%pw_grid%npts(1:3)
382  bo = spline_g%pw_grid%bounds_local
383  gbo = spline_g%pw_grid%bounds
384 
385  cpassert(.NOT. spline_g%pw_grid%spherical)
386  cpassert(spline_g%pw_grid%grid_span == fullspace)
387 
388  ALLOCATE (csivals(gbo(1, 1):gbo(2, 1)), &
389  csjvals(gbo(1, 2):gbo(2, 2)), &
390  cskvals(gbo(1, 3):gbo(2, 3)))
391 
392  coeff = twopi/n_tot(1)
393  IF (idir == 1) THEN
394 !$OMP PARALLEL DO PRIVATE(i) DEFAULT(NONE) SHARED(gbo,csIVals,coeff)
395  DO i = gbo(1, 1), gbo(2, 1)
396  csivals(i) = sin(coeff*real(i, dp))
397  END DO
398  ELSE
399 !$OMP PARALLEL DO PRIVATE(i) DEFAULT(NONE) SHARED(gbo,csIVals,coeff)
400  DO i = gbo(1, 1), gbo(2, 1)
401  csivals(i) = cos(coeff*real(i, dp))
402  END DO
403  END IF
404  coeff = twopi/n_tot(2)
405  IF (idir == 2) THEN
406 !$OMP PARALLEL DO PRIVATE(j) DEFAULT(NONE) SHARED(gbo,csJVals,coeff)
407  DO j = gbo(1, 2), gbo(2, 2)
408  csjvals(j) = sin(coeff*real(j, dp))
409  END DO
410  ELSE
411 !$OMP PARALLEL DO PRIVATE(j) DEFAULT(NONE) SHARED(gbo,csJVals,coeff)
412  DO j = gbo(1, 2), gbo(2, 2)
413  csjvals(j) = cos(coeff*real(j, dp))
414  END DO
415  END IF
416  coeff = twopi/n_tot(3)
417  IF (idir == 3) THEN
418 !$OMP PARALLEL DO PRIVATE(k) DEFAULT(NONE) SHARED(gbo,csKVals,coeff)
419  DO k = gbo(1, 3), gbo(2, 3)
420  cskvals(k) = sin(coeff*real(k, dp))
421  END DO
422  ELSE
423 !$OMP PARALLEL DO PRIVATE(k) DEFAULT(NONE) SHARED(gbo,csKVals,coeff)
424  DO k = gbo(1, 3), gbo(2, 3)
425  cskvals(k) = cos(coeff*real(k, dp))
426  END DO
427  END IF
428 
429  SELECT CASE (idir)
430  CASE (1)
431  ! x deriv
432 !$OMP PARALLEL DO PRIVATE(ii,k,j,i,coeff,tmp) DEFAULT(NONE) SHARED(spline_g,csIVals,csJVals,csKVals)
433  DO ii = 1, SIZE(spline_g%array)
434  i = spline_g%pw_grid%g_hat(1, ii)
435  j = spline_g%pw_grid%g_hat(2, ii)
436  k = spline_g%pw_grid%g_hat(3, ii)
437 !FM ! formula
438 !FM coeff=(sinVal(1)*cosVal(2)*cosVal(3))/9.0_dp+&
439 !FM (sinVal(1)*cosVal(2)+sinVal(1)*cosVal(3))*2.0_dp/9.0_dp+&
440 !FM sinVal(1)*4.0_dp/9.0_dp
441  tmp = csivals(i)*csjvals(j)
442  coeff = (tmp*cskvals(k) + &
443  (tmp + csivals(i)*cskvals(k))*2.0_dp + &
444  csivals(i)*4.0_dp)*inv9
445 
446  spline_g%array(ii) = spline_g%array(ii)* &
447  cmplx(0.0_dp, coeff, dp)
448  END DO
449  CASE (2)
450  ! y deriv
451 !$OMP PARALLEL DO PRIVATE(ii,k,j,i,coeff,tmp) DEFAULT(NONE) SHARED(spline_g,csIVals,csJVals,csKVals)
452  DO ii = 1, SIZE(spline_g%array)
453  i = spline_g%pw_grid%g_hat(1, ii)
454  j = spline_g%pw_grid%g_hat(2, ii)
455  k = spline_g%pw_grid%g_hat(3, ii)
456 
457  tmp = csivals(i)*csjvals(j)
458  coeff = (tmp*cskvals(k) + &
459  (tmp + csjvals(j)*cskvals(k))*2.0_dp + &
460  csjvals(j)*4.0_dp)*inv9
461 
462  spline_g%array(ii) = spline_g%array(ii)* &
463  cmplx(0.0_dp, coeff, dp)
464  END DO
465  CASE (3)
466  ! z deriv
467 !$OMP PARALLEL DO PRIVATE(ii,k,j,i,coeff,tmp) DEFAULT(NONE) SHARED(spline_g,csIVals,csJVals,csKVals)
468  DO ii = 1, SIZE(spline_g%array)
469  i = spline_g%pw_grid%g_hat(1, ii)
470  j = spline_g%pw_grid%g_hat(2, ii)
471  k = spline_g%pw_grid%g_hat(3, ii)
472 
473  tmp = csivals(i)*cskvals(k)
474  coeff = (tmp*csjvals(j) + &
475  (tmp + csjvals(j)*cskvals(k))*2.0_dp + &
476  cskvals(k)*4.0_dp)*inv9
477 
478  spline_g%array(ii) = spline_g%array(ii)* &
479  cmplx(0.0_dp, coeff, dp)
480  END DO
481  END SELECT
482 
483  DEALLOCATE (csivals, csjvals, cskvals)
484 
485  CALL timestop(handle)
486  END SUBROUTINE pw_spline3_deriv_g
487 
488 ! **************************************************************************************************
489 !> \brief calculates the FFT of the values of the x,y,z (idir=1,2,3)
490 !> derivative of the quadratic spline
491 !> \param spline_g on entry the FFT of the coefficients of the spline
492 !> will contain the FFT of the derivative
493 !> \param idir direction of the derivative
494 !> \par History
495 !> 06.2003 created [fawzi]
496 !> \author Fawzi Mohamed
497 !> \note
498 !> the distance between gridpoints is assumed to be 1
499 ! **************************************************************************************************
500  SUBROUTINE pw_spline2_deriv_g(spline_g, idir)
501  TYPE(pw_c1d_gs_type), INTENT(IN) :: spline_g
502  INTEGER, INTENT(in) :: idir
503 
504  CHARACTER(len=*), PARAMETER :: routinen = 'pw_spline2_deriv_g'
505  REAL(kind=dp), PARAMETER :: inv16 = 1.0_dp/16.0_dp
506 
507  INTEGER :: handle, i, ii, j, k
508  INTEGER, DIMENSION(2, 3) :: bo
509  INTEGER, DIMENSION(3) :: n, n_tot
510  REAL(kind=dp) :: coeff, tmp
511  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: csivals, csjvals, cskvals
512 
513  CALL timeset(routinen, handle)
514 
515  n(1:3) = spline_g%pw_grid%npts_local(1:3)
516  n_tot(1:3) = spline_g%pw_grid%npts(1:3)
517  bo = spline_g%pw_grid%bounds
518 
519  cpassert(.NOT. spline_g%pw_grid%spherical)
520  cpassert(spline_g%pw_grid%grid_span == fullspace)
521 
522  ALLOCATE (csivals(bo(1, 1):bo(2, 1)), csjvals(bo(1, 2):bo(2, 2)), &
523  cskvals(bo(1, 3):bo(2, 3)))
524 
525  coeff = twopi/n_tot(1)
526  IF (idir == 1) THEN
527 !$OMP PARALLEL DO PRIVATE(i) DEFAULT(NONE) SHARED(bo,coeff,csIVals)
528  DO i = bo(1, 1), bo(2, 1)
529  csivals(i) = sin(coeff*real(i, dp))
530  END DO
531  ELSE
532 !$OMP PARALLEL DO PRIVATE(i) DEFAULT(NONE) SHARED(bo,coeff,csIVals)
533  DO i = bo(1, 1), bo(2, 1)
534  csivals(i) = cos(coeff*real(i, dp))
535  END DO
536  END IF
537  coeff = twopi/n_tot(2)
538  IF (idir == 2) THEN
539 !$OMP PARALLEL DO PRIVATE(j) DEFAULT(NONE) SHARED(bo,coeff,csJVals)
540  DO j = bo(1, 2), bo(2, 2)
541  csjvals(j) = sin(coeff*real(j, dp))
542  END DO
543  ELSE
544 !$OMP PARALLEL DO PRIVATE(j) DEFAULT(NONE) SHARED(bo,coeff,csJVals)
545  DO j = bo(1, 2), bo(2, 2)
546  csjvals(j) = cos(coeff*real(j, dp))
547  END DO
548  END IF
549  coeff = twopi/n_tot(3)
550  IF (idir == 3) THEN
551 !$OMP PARALLEL DO PRIVATE(k) DEFAULT(NONE) SHARED(bo,coeff,csKVals)
552  DO k = bo(1, 3), bo(2, 3)
553  cskvals(k) = sin(coeff*real(k, dp))
554  END DO
555  ELSE
556 !$OMP PARALLEL DO PRIVATE(k) DEFAULT(NONE) SHARED(bo,coeff,csKVals)
557  DO k = bo(1, 3), bo(2, 3)
558  cskvals(k) = cos(coeff*real(k, dp))
559  END DO
560  END IF
561 
562  SELECT CASE (idir)
563  CASE (1)
564  ! x deriv
565 !$OMP PARALLEL DO PRIVATE(ii,k,j,i,coeff,tmp) SHARED(spline_g,csIVals,csJVals,csKVals) DEFAULT(NONE)
566  DO ii = 1, SIZE(spline_g%array)
567  i = spline_g%pw_grid%g_hat(1, ii)
568  j = spline_g%pw_grid%g_hat(2, ii)
569  k = spline_g%pw_grid%g_hat(3, ii)
570 !FM ! formula
571 !FM coeff=(sinVal(1)*cosVal(2)*cosVal(3))/16.0_dp+&
572 !FM (sinVal(1)*cosVal(2)+sinVal(1)*cosVal(3))*3.0_dp/16.0_dp+&
573 !FM sinVal(1)*9.0_dp/16.0_dp
574  tmp = csivals(i)*csjvals(j)
575  coeff = (tmp*cskvals(k) + &
576  (tmp + csivals(i)*cskvals(k))*3.0_dp + &
577  csivals(i)*9.0_dp)*inv16
578 
579  spline_g%array(ii) = spline_g%array(ii)* &
580  cmplx(0.0_dp, coeff, dp)
581  END DO
582  CASE (2)
583  ! y deriv
584 !$OMP PARALLEL DO PRIVATE(ii,k,j,i,coeff,tmp) DEFAULT(NONE) SHARED(spline_g,csIVals,csJVals,csKVals)
585  DO ii = 1, SIZE(spline_g%array)
586  i = spline_g%pw_grid%g_hat(1, ii)
587  j = spline_g%pw_grid%g_hat(2, ii)
588  k = spline_g%pw_grid%g_hat(3, ii)
589 
590  tmp = csivals(i)*csjvals(j)
591  coeff = (tmp*cskvals(k) + &
592  (tmp + csjvals(j)*cskvals(k))*3.0_dp + &
593  csjvals(j)*9.0_dp)*inv16
594 
595  spline_g%array(ii) = spline_g%array(ii)* &
596  cmplx(0.0_dp, coeff, dp)
597  END DO
598  CASE (3)
599  ! z deriv
600 !$OMP PARALLEL DO PRIVATE(ii,k,j,i,coeff,tmp) DEFAULT(NONE) SHARED(spline_g,csIVals,csJVals,csKVals)
601  DO ii = 1, SIZE(spline_g%array)
602  i = spline_g%pw_grid%g_hat(1, ii)
603  j = spline_g%pw_grid%g_hat(2, ii)
604  k = spline_g%pw_grid%g_hat(3, ii)
605 
606  tmp = csivals(i)*cskvals(k)
607  coeff = (tmp*csjvals(j) + &
608  (tmp + csjvals(j)*cskvals(k))*3.0_dp + &
609  cskvals(k)*9.0_dp)*inv16
610 
611  spline_g%array(ii) = spline_g%array(ii)* &
612  cmplx(0.0_dp, coeff, dp)
613  END DO
614  END SELECT
615 
616  DEALLOCATE (csivals, csjvals, cskvals)
617 
618  CALL timestop(handle)
619  END SUBROUTINE pw_spline2_deriv_g
620 
621 ! **************************************************************************************************
622 !> \brief applies a nearest neighbor linear operator to a stripe in x direction:
623 !> out_val(i)=sum(weight(j)*in_val(i+j-1),j=0..2)
624 !> \param weights the weights of the linear operator
625 !> \param in_val the argument to the operator
626 !> \param in_val_first the first argument (needed to calculate out_val(1))
627 !> \param in_val_last the last argument (needed to calculate out_val(n_el))
628 !> \param out_val the place where the result is accumulated
629 !> \param n_el the number of elements in in_v and out_v
630 !> \par History
631 !> 04.2004 created [fawzi]
632 !> \author fawzi
633 !> \note
634 !> uses 2 read streams and 1 write stream
635 ! **************************************************************************************************
636  SUBROUTINE pw_compose_stripe(weights, in_val, in_val_first, in_val_last, &
637  out_val, n_el)
638  REAL(kind=dp), DIMENSION(0:2), INTENT(in) :: weights
639  REAL(kind=dp), DIMENSION(*), INTENT(in) :: in_val
640  REAL(kind=dp), INTENT(in) :: in_val_first, in_val_last
641  REAL(kind=dp), DIMENSION(*), INTENT(inout) :: out_val
642  INTEGER :: n_el
643 
644  INTEGER :: i
645  REAL(kind=dp) :: v0, v1, v2
646 
647 !1:n_el), &
648 !1:n_el), &
649 
650  IF (n_el < 1) RETURN
651  v0 = in_val_first
652  v1 = in_val(1)
653  IF (weights(1) == 0.0_dp) THEN
654  ! optimized version for x deriv
655  DO i = 1, n_el - 3, 3
656  v2 = in_val(i + 1)
657  out_val(i) = out_val(i) + &
658  weights(0)*v0 + &
659  weights(2)*v2
660  v0 = in_val(i + 2)
661  out_val(i + 1) = out_val(i + 1) + &
662  weights(0)*v1 + &
663  weights(2)*v0
664  v1 = in_val(i + 3)
665  out_val(i + 2) = out_val(i + 2) + &
666  weights(0)*v2 + &
667  weights(2)*v1
668  END DO
669  ELSE
670  ! generic version
671  DO i = 1, n_el - 3, 3
672  v2 = in_val(i + 1)
673  out_val(i) = out_val(i) + &
674  weights(0)*v0 + &
675  weights(1)*v1 + &
676  weights(2)*v2
677  v0 = in_val(i + 2)
678  out_val(i + 1) = out_val(i + 1) + &
679  weights(0)*v1 + &
680  weights(1)*v2 + &
681  weights(2)*v0
682  v1 = in_val(i + 3)
683  out_val(i + 2) = out_val(i + 2) + &
684  weights(0)*v2 + &
685  weights(1)*v0 + &
686  weights(2)*v1
687  END DO
688  END IF
689  SELECT CASE (modulo(n_el - 1, 3))
690  CASE (0)
691  v2 = in_val_last
692  out_val(n_el) = out_val(n_el) + &
693  weights(0)*v0 + &
694  weights(1)*v1 + &
695  weights(2)*v2
696  CASE (1)
697  v2 = in_val(n_el)
698  out_val(n_el - 1) = out_val(n_el - 1) + &
699  weights(0)*v0 + &
700  weights(1)*v1 + &
701  weights(2)*v2
702  v0 = in_val_last
703  out_val(n_el) = out_val(n_el) + &
704  weights(0)*v1 + &
705  weights(1)*v2 + &
706  weights(2)*v0
707  CASE (2)
708  v2 = in_val(n_el - 1)
709  out_val(n_el - 2) = out_val(n_el - 2) + &
710  weights(0)*v0 + &
711  weights(1)*v1 + &
712  weights(2)*v2
713  v0 = in_val(n_el)
714  out_val(n_el - 1) = out_val(n_el - 1) + &
715  weights(0)*v1 + &
716  weights(1)*v2 + &
717  weights(2)*v0
718  v1 = in_val_last
719  out_val(n_el) = out_val(n_el) + &
720  weights(0)*v2 + &
721  weights(1)*v0 + &
722  weights(2)*v1
723  END SELECT
724 
725  END SUBROUTINE pw_compose_stripe
726 
727 ! **************************************************************************************************
728 !> \brief private routine that computes pw_nn_compose_r (it seems that without
729 !> passing arrays in this way either some compiler do a copyin/out (xlf)
730 !> or by inlining suboptimal code is produced (nag))
731 !> \param weights a 3x3x3 array with the linear operator
732 !> \param in_val the argument for the linear operator
733 !> \param out_val place where the value of the linear oprator should be added
734 !> \param pw_in pw to be able to get the needed meta data about in_val and
735 !> out_val
736 !> \param bo boundaries of in_val and out_val
737 !> \author fawzi
738 ! **************************************************************************************************
739  SUBROUTINE pw_nn_compose_r_work(weights, in_val, out_val, pw_in, bo)
740  REAL(kind=dp), DIMENSION(0:2, 0:2, 0:2) :: weights
741  INTEGER, DIMENSION(2, 3) :: bo
742  TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
743  REAL(kind=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, & 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(inout) :: out_val
744  REAL(kind=dp), DIMENSION(bo(1, 1):bo(2, 1), bo(1, & 2):bo(2, 2), bo(1, 3):bo(2, 3)), INTENT(in) :: in_val
745 
746  INTEGER :: i, j, jw, k, kw, myj, myk
747  INTEGER, DIMENSION(2, 3) :: gbo
748  INTEGER, DIMENSION(3) :: s
749  LOGICAL :: has_boundary, yderiv, zderiv
750  REAL(kind=dp) :: in_val_f, in_val_l
751  REAL(kind=dp), DIMENSION(:, :), POINTER :: l_boundary, tmp, u_boundary
752 
753  zderiv = all(weights(:, :, 1) == 0.0_dp)
754  yderiv = all(weights(:, 1, :) == 0.0_dp)
755  bo = pw_in%pw_grid%bounds_local
756  gbo = pw_in%pw_grid%bounds
757  DO i = 1, 3
758  s(i) = bo(2, i) - bo(1, i) + 1
759  END DO
760  IF (any(s < 1)) RETURN
761  has_boundary = any(pw_in%pw_grid%bounds_local(:, 1) /= &
762  pw_in%pw_grid%bounds(:, 1))
763  IF (has_boundary) THEN
764  ALLOCATE (l_boundary(bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)), &
765  u_boundary(bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)), &
766  tmp(bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
767  tmp(:, :) = pw_in%array(bo(2, 1), :, :)
768  CALL pw_in%pw_grid%para%group%sendrecv(tmp, pw_in%pw_grid%para%pos_of_x( &
769  gbo(1, 1) + modulo(bo(2, 1) + 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)), &
770  l_boundary, pw_in%pw_grid%para%pos_of_x( &
771  gbo(1, 1) + modulo(bo(1, 1) - 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)))
772  tmp(:, :) = pw_in%array(bo(1, 1), :, :)
773  CALL pw_in%pw_grid%para%group%sendrecv(tmp, pw_in%pw_grid%para%pos_of_x( &
774  gbo(1, 1) + modulo(bo(1, 1) - 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)), &
775  u_boundary, pw_in%pw_grid%para%pos_of_x( &
776  gbo(1, 1) + modulo(bo(2, 1) + 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)))
777  DEALLOCATE (tmp)
778  END IF
779 
780 !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(k,kw,myk,j,jw,myj,in_val_f,&
781 !$OMP in_val_l) SHARED(zderiv,yderiv,bo,in_val,out_val,s,l_boundary,&
782 !$OMP u_boundary,weights,has_boundary)
783  DO k = 0, s(3) - 1
784  DO kw = 0, 2
785  myk = bo(1, 3) + modulo(k + kw - 1, s(3))
786  IF (zderiv .AND. kw == 1) cycle
787  DO j = 0, s(2) - 1
788  DO jw = 0, 2
789  myj = bo(1, 2) + modulo(j + jw - 1, s(2))
790  IF (yderiv .AND. jw == 1) cycle
791  IF (has_boundary) THEN
792  in_val_f = l_boundary(myj, myk)
793  in_val_l = u_boundary(myj, myk)
794  ELSE
795  in_val_f = in_val(bo(2, 1), myj, myk)
796  in_val_l = in_val(bo(1, 1), myj, myk)
797  END IF
798  CALL pw_compose_stripe(weights=weights(:, jw, kw), &
799  in_val=in_val(:, myj, myk), &
800  in_val_first=in_val_f, in_val_last=in_val_l, &
801  out_val=out_val(:, bo(1, 2) + j, bo(1, 3) + k), n_el=s(1))
802  END DO
803  END DO
804  END DO
805  END DO
806  IF (has_boundary) THEN
807  DEALLOCATE (l_boundary, u_boundary)
808  END IF
809  END SUBROUTINE pw_nn_compose_r_work
810 
811 ! **************************************************************************************************
812 !> \brief applies a nearest neighbor linear operator to a pw in real space
813 !> \param weights a 3x3x3 array with the linear operator
814 !> \param pw_in the argument for the linear operator
815 !> \param pw_out place where the value of the linear oprator should be added
816 !> \author fawzi
817 !> \note
818 !> has specialized versions for derivative operator (with central values==0)
819 ! **************************************************************************************************
820  SUBROUTINE pw_nn_compose_r(weights, pw_in, pw_out)
821  REAL(kind=dp), DIMENSION(0:2, 0:2, 0:2) :: weights
822  TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in, pw_out
823 
824  CHARACTER(len=*), PARAMETER :: routinen = 'pw_nn_compose_r'
825 
826  INTEGER :: handle
827 
828  CALL timeset(routinen, handle)
829  IF (.NOT. all(pw_in%pw_grid%bounds_local(:, 2:3) == pw_in%pw_grid%bounds(:, 2:3))) THEN
830  cpabort("wrong pw distribution")
831  END IF
832  CALL pw_nn_compose_r_work(weights=weights, in_val=pw_in%array, &
833  out_val=pw_out%array, pw_in=pw_in, bo=pw_in%pw_grid%bounds_local)
834  CALL timestop(handle)
835  END SUBROUTINE pw_nn_compose_r
836 
837 ! **************************************************************************************************
838 !> \brief calculates the values of a nearest neighbor smearing
839 !> \param pw_in the argument for the linear operator
840 !> \param pw_out place where the smeared values should be added
841 !> \param coeffs array with the coefficent of the smearing, ordered with
842 !> the distance from the center: coeffs(1) the coeff of the central
843 !> element, coeffs(2) the coeff of the 6 element with distance 1,
844 !> coeff(3) the coeff of the 12 elements at distance sqrt(2),
845 !> coeff(4) the coeff of the 8 elements at distance sqrt(3).
846 !> \author Fawzi Mohamed
847 !> \note
848 !> does not normalize the smear to 1.
849 !> with coeff=(/ 8._dp/27._dp, 2._dp/27._dp, 1._dp/54._dp, 1._dp/216._dp /)
850 !> is equivalent to pw_spline3_evaluate_values_g, with
851 !> coeff=(/ 27._dp/64._dp, 9._dp/128._dp, 3._dp/256._dp, 1._dp/512._dp /)
852 ! **************************************************************************************************
853  SUBROUTINE pw_nn_smear_r(pw_in, pw_out, coeffs)
854  TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in, pw_out
855  REAL(kind=dp), DIMENSION(4), INTENT(in) :: coeffs
856 
857  INTEGER :: i, j, k
858  REAL(kind=dp), DIMENSION(-1:1, -1:1, -1:1) :: weights
859 
860  DO k = -1, 1
861  DO j = -1, 1
862  DO i = -1, 1
863  weights(i, j, k) = coeffs(abs(i) + abs(j) + abs(k) + 1)
864  END DO
865  END DO
866  END DO
867 
868  CALL pw_nn_compose_r(weights=weights, pw_in=pw_in, pw_out=pw_out)
869  END SUBROUTINE pw_nn_smear_r
870 
871 ! **************************************************************************************************
872 !> \brief calculates a nearest neighbor central derivative.
873 !> for the x dir:
874 !> pw_out%array(i,j,k)=( pw_in(i+1,j,k)-pw_in(i-1,j,k) )*coeff(1)+
875 !> ( pw_in(i+1,j(+-)1,k)-pw_in(i-1,j(+-)1,k)+
876 !> pw_in(i+1,j,k(+-)1)-pw_in(i-1,j,k(+-)1) )*coeff(2)+
877 !> ( pw_in(i+1,j(+-)1,k(+-)1)-pw_in(i-1,j(+-)1,k(+-)1)+
878 !> pw_in(i+1,j(+-)1,k(-+)1)-pw_in(i-1,j(+-)1,k(-+)1) )*coeff(3)
879 !> periodic boundary conditions are applied
880 !> \param pw_in the argument for the linear operator
881 !> \param pw_out place where the smeared values should be added
882 !> \param coeffs array with the coefficent of the front (positive) plane
883 !> of the central derivative, ordered with
884 !> the distance from the center: coeffs(1) the coeff of the central
885 !> element, coeffs(2) the coeff of the 4 element with distance 1,
886 !> coeff(3) the coeff of the 4 elements at distance sqrt(2)
887 !> \param idir ...
888 !> \author Fawzi Mohamed
889 !> \note
890 !> with coeff=(/ 2.0_dp/9.0_dp, 1.0_dp/18.0_dp, 1.0_dp/72.0_dp /)
891 !> is equivalent to pw_spline3_deriv_r, with
892 !> coeff=(/ 9.0_dp/32.0_dp, 3.0_dp/64.0_dp, 1.0_dp/128.0_dp /)
893 !> to pw_spline2_deriv_r
894 !> coeff=(/ 25._dp/72._dp, 5._dp/144, 1._dp/288._dp /)
895 ! **************************************************************************************************
896  SUBROUTINE pw_nn_deriv_r(pw_in, pw_out, coeffs, idir)
897  TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in, pw_out
898  REAL(kind=dp), DIMENSION(3), INTENT(in) :: coeffs
899  INTEGER :: idir
900 
901  INTEGER :: i, idirval, j, k
902  REAL(kind=dp), DIMENSION(-1:1, -1:1, -1:1) :: weights
903 
904  DO k = -1, 1
905  DO j = -1, 1
906  DO i = -1, 1
907  SELECT CASE (idir)
908  CASE (1)
909  idirval = i
910  CASE (2)
911  idirval = j
912  CASE (3)
913  idirval = k
914  CASE default
915  cpabort("invalid idir ("//trim(cp_to_string(idir))//")")
916  END SELECT
917  IF (idirval == 0) THEN
918  weights(i, j, k) = 0.0_dp
919  ELSE
920  weights(i, j, k) = real(idirval, dp)*coeffs(abs(i) + abs(j) + abs(k))
921  END IF
922  END DO
923  END DO
924  END DO
925 
926  CALL pw_nn_compose_r(weights=weights, pw_in=pw_in, pw_out=pw_out)
927  END SUBROUTINE pw_nn_deriv_r
928 
929 ! **************************************************************************************************
930 !> \brief low level function that adds a coarse grid
931 !> to a fine grid.
932 !> If pbc is true periodic boundary conditions are applied
933 !>
934 !> It will add to
935 !>
936 !> fine_values(2*coarse_bounds(1,1):2*coarse_bounds(2,1),
937 !> 2*coarse_bounds(1,2):2*coarse_bounds(2,2),
938 !> 2*coarse_bounds(1,3):2*coarse_bounds(2,3))
939 !>
940 !> using
941 !>
942 !> coarse_coeffs(coarse_bounds(1,1):coarse_bounds(2,1),
943 !> coarse_bounds(1,2):coarse_bounds(2,2),
944 !> coarse_bounds(1,3):coarse_bounds(2,3))
945 !>
946 !> composed with the weights obtained by the direct product of the
947 !> 1d coefficients weights:
948 !>
949 !> for i,j,k in -3..3
950 !> w(i,j,k)=weights_1d(abs(i)+1)*weights_1d(abs(j)+1)*
951 !> weights_1d(abs(k)+1)
952 !> \param coarse_coeffs_pw the values of the coefficients
953 !> \param fine_values_pw where to add the values due to the
954 !> coarse coeffs
955 !> \param weights_1d the weights of the 1d smearing
956 !> \param w_border0 the 1d weight at the border (when pbc is false)
957 !> \param w_border1 the 1d weights for a point one off the border
958 !> (w_border1(1) is the weight of the coefficent at the border)
959 !> (used if pbc is false)
960 !> \param pbc if periodic boundary conditions should be applied
961 !> \param safe_computation ...
962 !> \author fawzi
963 !> \note
964 !> coarse looping is continuos, I did not check if keeping the fine looping
965 !> contiguous is better.
966 !> And I ask myself yet again why, why we use x-slice distribution,
967 !> z-slice distribution would be much better performancewise
968 !> (and would semplify this code enormously).
969 !> fine2coarse has much more understandable parallel part (build up of
970 !> send/rcv sizes,... but worse if you have really a lot of processors,
971 !> probabily irrelevant because it is not critical) [fawzi].
972 ! **************************************************************************************************
973  SUBROUTINE add_coarse2fine(coarse_coeffs_pw, fine_values_pw, &
974  weights_1d, w_border0, w_border1, pbc, safe_computation)
975  TYPE(pw_r3d_rs_type), INTENT(IN) :: coarse_coeffs_pw, fine_values_pw
976  REAL(kind=dp), DIMENSION(4), INTENT(in) :: weights_1d
977  REAL(kind=dp), INTENT(in) :: w_border0
978  REAL(kind=dp), DIMENSION(3), INTENT(in) :: w_border1
979  LOGICAL, INTENT(in) :: pbc
980  LOGICAL, INTENT(in), OPTIONAL :: safe_computation
981 
982  CHARACTER(len=*), PARAMETER :: routinen = 'add_coarse2fine'
983 
984  INTEGER :: coarse_slice_size, f_shift(3), fi, fi_lb, fi_ub, fj, fk, handle, handle2, i, ii, &
985  ij, ik, ip, j, k, my_lb, my_ub, n_procs, p, p_lb, p_old, p_ub, rcv_tot_size, rest_b, &
986  s(3), send_tot_size, sf, shift, ss, x, x_att, xx
987  INTEGER, ALLOCATABLE, DIMENSION(:) :: rcv_offset, rcv_size, real_rcv_size, &
988  send_offset, send_size, sent_size
989  INTEGER, DIMENSION(2, 3) :: coarse_bo, coarse_gbo, fine_bo, &
990  fine_gbo, my_coarse_bo
991  INTEGER, DIMENSION(:), POINTER :: pos_of_x
992  LOGICAL :: has_i_lbound, has_i_ubound, is_split, &
993  safe_calc
994  REAL(kind=dp) :: v0, v1, v2, v3, wi, wj, wk
995  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rcv_buf, send_buf
996  REAL(kind=dp), DIMENSION(3) :: w_0, ww0
997  REAL(kind=dp), DIMENSION(4) :: w_1, ww1
998  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: coarse_coeffs, fine_values
999 
1000  CALL timeset(routinen, handle)
1001 ! CALL timeset(routineN//"_pre",handle2)
1002  safe_calc = .false.
1003  IF (PRESENT(safe_computation)) safe_calc = safe_computation
1004  ii = coarse_coeffs_pw%pw_grid%para%group%compare(fine_values_pw%pw_grid%para%group)
1005  IF (ii > mp_comm_congruent) THEN
1006  cpabort("")
1007  END IF
1008  my_coarse_bo = coarse_coeffs_pw%pw_grid%bounds_local
1009  coarse_gbo = coarse_coeffs_pw%pw_grid%bounds
1010  fine_bo = fine_values_pw%pw_grid%bounds_local
1011  fine_gbo = fine_values_pw%pw_grid%bounds
1012  f_shift = fine_gbo(1, :) - 2*coarse_gbo(1, :)
1013  DO j = 2, 3
1014  DO i = 1, 2
1015  coarse_bo(i, j) = floor((fine_bo(i, j) - f_shift(j))/2.)
1016  END DO
1017  END DO
1018  IF (fine_bo(1, 1) <= fine_bo(2, 1)) THEN
1019  coarse_bo(1, 1) = floor((fine_bo(1, 1) - 2 - f_shift(1))/2.)
1020  coarse_bo(2, 1) = floor((fine_bo(2, 1) + 3 - f_shift(1))/2.)
1021  ELSE
1022  coarse_bo(1, 1) = coarse_gbo(2, 1)
1023  coarse_bo(2, 1) = coarse_gbo(2, 1) - 1
1024  END IF
1025  is_split = any(coarse_gbo(:, 1) /= my_coarse_bo(:, 1))
1026  IF (.NOT. is_split .OR. .NOT. pbc) THEN
1027  coarse_bo(1, 1) = max(coarse_gbo(1, 1), coarse_bo(1, 1))
1028  coarse_bo(2, 1) = min(coarse_gbo(2, 1), coarse_bo(2, 1))
1029  END IF
1030  has_i_ubound = (fine_gbo(2, 1) /= fine_bo(2, 1)) .OR. pbc .AND. is_split
1031  has_i_lbound = (fine_gbo(1, 1) /= fine_bo(1, 1)) .OR. pbc .AND. is_split
1032 
1033  IF (pbc) THEN
1034  cpassert(all(fine_gbo(1, :) == 2*coarse_gbo(1, :) + f_shift))
1035  cpassert(all(fine_gbo(2, :) == 2*coarse_gbo(2, :) + 1 + f_shift))
1036  ELSE
1037  cpassert(all(fine_gbo(2, :) == 2*coarse_gbo(2, :) + f_shift))
1038  cpassert(all(fine_gbo(1, :) == 2*coarse_gbo(1, :) + f_shift))
1039  END IF
1040 
1041  coarse_coeffs => coarse_coeffs_pw%array
1042  DO i = 1, 3
1043  s(i) = coarse_gbo(2, i) - coarse_gbo(1, i) + 1
1044  END DO
1045 ! CALL timestop(handle2)
1046  ! *** parallel case
1047  IF (is_split) THEN
1048  CALL timeset(routinen//"_comm", handle2)
1049  coarse_slice_size = (coarse_bo(2, 2) - coarse_bo(1, 2) + 1)* &
1050  (coarse_bo(2, 3) - coarse_bo(1, 3) + 1)
1051  n_procs = coarse_coeffs_pw%pw_grid%para%group_size
1052  ALLOCATE (send_size(0:n_procs - 1), send_offset(0:n_procs - 1), &
1053  sent_size(0:n_procs - 1), rcv_size(0:n_procs - 1), &
1054  rcv_offset(0:n_procs - 1), real_rcv_size(0:n_procs - 1))
1055 
1056  ! ** rcv size count
1057 
1058  pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x
1059  p_old = pos_of_x(coarse_gbo(1, 1) &
1060  + modulo(coarse_bo(1, 1) - coarse_gbo(1, 1), s(1)))
1061  rcv_size = 0
1062  DO x = coarse_bo(1, 1), coarse_bo(2, 1)
1063  p = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
1064  rcv_size(p) = rcv_size(p) + coarse_slice_size
1065  END DO
1066 
1067  ! ** send size count
1068 
1069  pos_of_x => fine_values_pw%pw_grid%para%pos_of_x
1070  sf = fine_gbo(2, 1) - fine_gbo(1, 1) + 1
1071  fi_lb = 2*my_coarse_bo(1, 1) - 3 + f_shift(1)
1072  fi_ub = 2*my_coarse_bo(2, 1) + 3 + f_shift(1)
1073  IF (.NOT. pbc) THEN
1074  fi_lb = max(fi_lb, fine_gbo(1, 1))
1075  fi_ub = min(fi_ub, fine_gbo(2, 1))
1076  ELSE
1077  fi_ub = min(fi_ub, fi_lb + sf - 1)
1078  END IF
1079  p_old = pos_of_x(fine_gbo(1, 1) + modulo(fi_lb - fine_gbo(1, 1), sf))
1080  p_lb = floor((fi_lb - 2 - f_shift(1))/2.)
1081  send_size = 0
1082  DO x = fi_lb, fi_ub
1083  p = pos_of_x(fine_gbo(1, 1) + modulo(x - fine_gbo(1, 1), sf))
1084  IF (p /= p_old) THEN
1085  p_ub = floor((x - 1 + 3 - f_shift(1))/2.)
1086 
1087  send_size(p_old) = send_size(p_old) + (min(p_ub, my_coarse_bo(2, 1)) &
1088  - max(p_lb, my_coarse_bo(1, 1)) + 1)*coarse_slice_size
1089 
1090  IF (pbc) THEN
1091  DO xx = p_lb, coarse_gbo(1, 1) - 1
1092  x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1093  IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1094  send_size(p_old) = send_size(p_old) + coarse_slice_size
1095  END IF
1096  END DO
1097  DO xx = coarse_gbo(2, 1) + 1, p_ub
1098  x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1099  IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1100  send_size(p_old) = send_size(p_old) + coarse_slice_size
1101  END IF
1102  END DO
1103  END IF
1104 
1105  p_old = p
1106  p_lb = floor((x - 2 - f_shift(1))/2.)
1107  END IF
1108  END DO
1109  p_ub = floor((fi_ub + 3 - f_shift(1))/2.)
1110 
1111  send_size(p_old) = send_size(p_old) + (min(p_ub, my_coarse_bo(2, 1)) &
1112  - max(p_lb, my_coarse_bo(1, 1)) + 1)*coarse_slice_size
1113 
1114  IF (pbc) THEN
1115  DO xx = p_lb, coarse_gbo(1, 1) - 1
1116  x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1117  IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1118  send_size(p_old) = send_size(p_old) + coarse_slice_size
1119  END IF
1120  END DO
1121  DO xx = coarse_gbo(2, 1) + 1, p_ub
1122  x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1123  IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1124  send_size(p_old) = send_size(p_old) + coarse_slice_size
1125  END IF
1126  END DO
1127  END IF
1128  ! ** offsets & alloc send-rcv
1129 
1130  send_tot_size = 0
1131  DO ip = 0, n_procs - 1
1132  send_offset(ip) = send_tot_size
1133  send_tot_size = send_tot_size + send_size(ip)
1134  END DO
1135  ALLOCATE (send_buf(0:send_tot_size - 1))
1136 
1137  rcv_tot_size = 0
1138  DO ip = 0, n_procs - 1
1139  rcv_offset(ip) = rcv_tot_size
1140  rcv_tot_size = rcv_tot_size + rcv_size(ip)
1141  END DO
1142  IF (.NOT. rcv_tot_size == (coarse_bo(2, 1) - coarse_bo(1, 1) + 1)*coarse_slice_size) THEN
1143  cpabort("Error calculating rcv_tot_size ")
1144  END IF
1145  ALLOCATE (rcv_buf(0:rcv_tot_size - 1))
1146 
1147  ! ** fill send buffer
1148 
1149  p_old = pos_of_x(fine_gbo(1, 1) + modulo(fi_lb - fine_gbo(1, 1), sf))
1150  p_lb = floor((fi_lb - 2 - f_shift(1))/2.)
1151  sent_size(:) = send_offset
1152  ss = my_coarse_bo(2, 1) - my_coarse_bo(1, 1) + 1
1153  DO x = fi_lb, fi_ub
1154  p = pos_of_x(fine_gbo(1, 1) + modulo(x - fine_gbo(1, 1), sf))
1155  IF (p /= p_old) THEN
1156  shift = floor((fine_gbo(1, 1) + modulo(x - 1 - fine_gbo(1, 1), sf) - f_shift(1))/2._dp) - &
1157  floor((x - 1 - f_shift(1))/2._dp)
1158  p_ub = floor((x - 1 + 3 - f_shift(1))/2._dp)
1159 
1160  IF (pbc) THEN
1161  DO xx = p_lb + shift, coarse_gbo(1, 1) - 1
1162  x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), sf)
1163  IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1164  CALL dcopy(coarse_slice_size, &
1165  coarse_coeffs(x_att, my_coarse_bo(1, 2), &
1166  my_coarse_bo(1, 3)), ss, send_buf(sent_size(p_old)), 1)
1167  sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1168  END IF
1169  END DO
1170  END IF
1171 
1172  ii = sent_size(p_old)
1173  DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1174  DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1175  DO i = max(p_lb + shift, my_coarse_bo(1, 1)), min(p_ub + shift, my_coarse_bo(2, 1))
1176  send_buf(ii) = coarse_coeffs(i, j, k)
1177  ii = ii + 1
1178  END DO
1179  END DO
1180  END DO
1181  sent_size(p_old) = ii
1182 
1183  IF (pbc) THEN
1184  DO xx = coarse_gbo(2, 1) + 1, p_ub + shift
1185  x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1186  IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1187  CALL dcopy(coarse_slice_size, &
1188  coarse_coeffs(x_att, my_coarse_bo(1, 2), &
1189  my_coarse_bo(1, 3)), ss, &
1190  send_buf(sent_size(p_old)), 1)
1191  sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1192  END IF
1193  END DO
1194  END IF
1195 
1196  p_old = p
1197  p_lb = floor((x - 2 - f_shift(1))/2.)
1198  END IF
1199  END DO
1200  shift = floor((fine_gbo(1, 1) + modulo(x - 1 - fine_gbo(1, 1), sf) - f_shift(1))/2._dp) - &
1201  floor((x - 1 - f_shift(1))/2._dp)
1202  p_ub = floor((fi_ub + 3 - f_shift(1))/2.)
1203 
1204  IF (pbc) THEN
1205  DO xx = p_lb + shift, coarse_gbo(1, 1) - 1
1206  x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1207  IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1208  CALL dcopy(coarse_slice_size, &
1209  coarse_coeffs(x_att, my_coarse_bo(1, 2), &
1210  my_coarse_bo(1, 3)), ss, send_buf(sent_size(p_old)), 1)
1211  sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1212  END IF
1213  END DO
1214  END IF
1215 
1216  ii = sent_size(p_old)
1217  DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1218  DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1219  DO i = max(p_lb + shift, my_coarse_bo(1, 1)), min(p_ub + shift, my_coarse_bo(2, 1))
1220  send_buf(ii) = coarse_coeffs(i, j, k)
1221  ii = ii + 1
1222  END DO
1223  END DO
1224  END DO
1225  sent_size(p_old) = ii
1226 
1227  IF (pbc) THEN
1228  DO xx = coarse_gbo(2, 1) + 1, p_ub + shift
1229  x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1230  IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1231  CALL dcopy(coarse_slice_size, &
1232  coarse_coeffs(x_att, my_coarse_bo(1, 2), &
1233  my_coarse_bo(1, 3)), ss, send_buf(sent_size(p_old)), 1)
1234  sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1235  END IF
1236  END DO
1237  END IF
1238 
1239  cpassert(all(sent_size(:n_procs - 2) == send_offset(1:)))
1240  cpassert(sent_size(n_procs - 1) == send_tot_size)
1241  ! test send/rcv sizes
1242  CALL coarse_coeffs_pw%pw_grid%para%group%alltoall(send_size, real_rcv_size, 1)
1243  cpassert(all(real_rcv_size == rcv_size))
1244  ! all2all
1245  CALL coarse_coeffs_pw%pw_grid%para%group%alltoall(sb=send_buf, scount=send_size, sdispl=send_offset, &
1246  rb=rcv_buf, rcount=rcv_size, rdispl=rcv_offset)
1247 
1248  ! ** reorder rcv buffer
1249  ! (actually reordering should be needed only with pbc)
1250 
1251  ALLOCATE (coarse_coeffs(coarse_bo(1, 1):coarse_bo(2, 1), &
1252  coarse_bo(1, 2):coarse_bo(2, 2), &
1253  coarse_bo(1, 3):coarse_bo(2, 3)))
1254 
1255  my_lb = max(coarse_gbo(1, 1), coarse_bo(1, 1))
1256  my_ub = min(coarse_gbo(2, 1), coarse_bo(2, 1))
1257  pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x
1258  sent_size(:) = rcv_offset
1259  ss = coarse_bo(2, 1) - coarse_bo(1, 1) + 1
1260  DO x = my_ub + 1, coarse_bo(2, 1)
1261  p_old = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
1262  CALL dcopy(coarse_slice_size, &
1263  rcv_buf(sent_size(p_old)), 1, &
1264  coarse_coeffs(x, coarse_bo(1, 2), &
1265  coarse_bo(1, 3)), ss)
1266  sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1267  END DO
1268  p_old = pos_of_x(coarse_gbo(1, 1) &
1269  + modulo(my_lb - coarse_gbo(1, 1), s(1)))
1270  p_lb = my_lb
1271  DO x = my_lb, my_ub
1272  p = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
1273  IF (p /= p_old) THEN
1274  p_ub = x - 1
1275 
1276  ii = sent_size(p_old)
1277  DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1278  DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1279  DO i = p_lb, p_ub
1280  coarse_coeffs(i, j, k) = rcv_buf(ii)
1281  ii = ii + 1
1282  END DO
1283  END DO
1284  END DO
1285  sent_size(p_old) = ii
1286 
1287  p_lb = x
1288  p_old = p
1289  END IF
1290  rcv_size(p) = rcv_size(p) + coarse_slice_size
1291  END DO
1292  p_ub = my_ub
1293  ii = sent_size(p_old)
1294  DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1295  DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1296  DO i = p_lb, p_ub
1297  coarse_coeffs(i, j, k) = rcv_buf(ii)
1298  ii = ii + 1
1299  END DO
1300  END DO
1301  END DO
1302  sent_size(p_old) = ii
1303  DO x = coarse_bo(1, 1), my_lb - 1
1304  p_old = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
1305  CALL dcopy(coarse_slice_size, &
1306  rcv_buf(sent_size(p_old)), 1, &
1307  coarse_coeffs(x, coarse_bo(1, 2), &
1308  coarse_bo(1, 3)), ss)
1309  sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1310  END DO
1311 
1312  cpassert(all(sent_size(0:n_procs - 2) == rcv_offset(1:)))
1313  cpassert(sent_size(n_procs - 1) == rcv_tot_size)
1314 
1315  ! dealloc
1316  DEALLOCATE (send_size, send_offset, rcv_size, rcv_offset)
1317  DEALLOCATE (send_buf, rcv_buf, real_rcv_size)
1318  CALL timestop(handle2)
1319 
1320  END IF
1321  fine_values => fine_values_pw%array
1322  w_0 = (/weights_1d(3), weights_1d(1), weights_1d(3)/)
1323  w_1 = (/weights_1d(4), weights_1d(2), weights_1d(2), weights_1d(4)/)
1324 
1325  DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1326  DO ik = -3, 3
1327  IF (pbc) THEN
1328  wk = weights_1d(abs(ik) + 1)
1329  fk = fine_gbo(1, 3) + modulo(2*k + ik - fine_gbo(1, 3) + f_shift(3), 2*s(3))
1330  ELSE
1331  fk = 2*k + ik + f_shift(3)
1332  IF (fk <= fine_bo(1, 3) + 1 .OR. fk >= fine_bo(2, 3) - 1) THEN
1333  IF (fk < fine_bo(1, 3) .OR. fk > fine_bo(2, 3)) cycle
1334  IF (fk == fine_bo(1, 3) .OR. fk == fine_bo(2, 3)) THEN
1335  IF (ik /= 0) cycle
1336  wk = w_border0
1337  ELSE IF (fk == 2*coarse_bo(1, 3) + 1 + f_shift(3)) THEN
1338  SELECT CASE (ik)
1339  CASE (1)
1340  wk = w_border1(1)
1341  CASE (-1)
1342  wk = w_border1(2)
1343  CASE (-3)
1344  wk = w_border1(3)
1345  CASE default
1346  cpabort("")
1347  cycle
1348  END SELECT
1349  ELSE
1350  SELECT CASE (ik)
1351  CASE (3)
1352  wk = w_border1(3)
1353  CASE (1)
1354  wk = w_border1(2)
1355  CASE (-1)
1356  wk = w_border1(1)
1357  CASE default
1358  cpabort("")
1359  cycle
1360  END SELECT
1361  END IF
1362  ELSE
1363  wk = weights_1d(abs(ik) + 1)
1364  END IF
1365  END IF
1366  DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1367  DO ij = -3, 3
1368  IF (pbc) THEN
1369  wj = weights_1d(abs(ij) + 1)*wk
1370  fj = fine_gbo(1, 2) + modulo(2*j + ij - fine_gbo(1, 2) + f_shift(2), 2*s(2))
1371  ELSE
1372  fj = 2*j + ij + f_shift(2)
1373  IF (fj <= fine_bo(1, 2) + 1 .OR. fj >= fine_bo(2, 2) - 1) THEN
1374  IF (fj < fine_bo(1, 2) .OR. fj > fine_bo(2, 2)) cycle
1375  IF (fj == fine_bo(1, 2) .OR. fj == fine_bo(2, 2)) THEN
1376  IF (ij /= 0) cycle
1377  wj = w_border0*wk
1378  ELSE IF (fj == 2*coarse_bo(1, 2) + 1 + f_shift(2)) THEN
1379  SELECT CASE (ij)
1380  CASE (1)
1381  wj = w_border1(1)*wk
1382  CASE (-1)
1383  wj = w_border1(2)*wk
1384  CASE (-3)
1385  wj = w_border1(3)*wk
1386  CASE default
1387  cycle
1388  END SELECT
1389  ELSE
1390  SELECT CASE (ij)
1391  CASE (-1)
1392  wj = w_border1(1)*wk
1393  CASE (1)
1394  wj = w_border1(2)*wk
1395  CASE (3)
1396  wj = w_border1(3)*wk
1397  CASE default
1398  cycle
1399  END SELECT
1400  END IF
1401  ELSE
1402  wj = weights_1d(abs(ij) + 1)*wk
1403  END IF
1404  END IF
1405 
1406  IF (fine_bo(2, 1) - fine_bo(1, 1) < 7 .OR. safe_calc) THEN
1407 ! CALL timeset(routineN//"_safe",handle2)
1408  DO i = coarse_bo(1, 1), coarse_bo(2, 1)
1409  DO ii = -3, 3
1410  IF (pbc .AND. .NOT. is_split) THEN
1411  wi = weights_1d(abs(ii) + 1)*wj
1412  fi = fine_gbo(1, 1) + modulo(2*i + ii - fine_gbo(1, 1) + f_shift(1), 2*s(1))
1413  ELSE
1414  fi = 2*i + ii + f_shift(1)
1415  IF (fi < fine_bo(1, 1) .OR. fi > fine_bo(2, 1)) cycle
1416  IF (.NOT. pbc .AND. (fi <= fine_gbo(1, 1) + 1 .OR. &
1417  fi >= fine_gbo(2, 1) - 1)) THEN
1418  IF (fi == fine_gbo(1, 1) .OR. fi == fine_gbo(2, 1)) THEN
1419  IF (ii /= 0) cycle
1420  wi = w_border0*wj
1421  ELSE IF (fi == fine_gbo(1, 1) + 1) THEN
1422  SELECT CASE (ii)
1423  CASE (1)
1424  wi = w_border1(1)*wj
1425  CASE (-1)
1426  wi = w_border1(2)*wj
1427  CASE (-3)
1428  wi = w_border1(3)*wj
1429  CASE default
1430  cycle
1431  END SELECT
1432  ELSE
1433  SELECT CASE (ii)
1434  CASE (-1)
1435  wi = w_border1(1)*wj
1436  CASE (1)
1437  wi = w_border1(2)*wj
1438  CASE (3)
1439  wi = w_border1(3)*wj
1440  CASE default
1441  cycle
1442  END SELECT
1443  END IF
1444  ELSE
1445  wi = weights_1d(abs(ii) + 1)*wj
1446  END IF
1447  END IF
1448  fine_values(fi, fj, fk) = &
1449  fine_values(fi, fj, fk) + &
1450  wi*coarse_coeffs(i, j, k)
1451  END DO
1452  END DO
1453 ! CALL timestop(handle2)
1454  ELSE
1455 ! CALL timeset(routineN//"_core1",handle2)
1456  ww0 = wj*w_0
1457  ww1 = wj*w_1
1458  IF (pbc .AND. .NOT. is_split) THEN
1459  v3 = coarse_coeffs(coarse_bo(2, 1), j, k)
1460  i = coarse_bo(1, 1)
1461  fi = 2*i + f_shift(1)
1462  v0 = coarse_coeffs(i, j, k)
1463  v1 = coarse_coeffs(i + 1, j, k)
1464  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1465  ww0(1)*v3 + ww0(2)*v0 + ww0(3)*v1
1466  v2 = coarse_coeffs(i + 2, j, k)
1467  fi = fi + 1
1468  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1469  ww1(1)*v3 + ww1(2)*v0 + ww1(3)*v1 + ww1(4)*v2
1470  ELSE IF (.NOT. has_i_lbound) THEN
1471  i = coarse_bo(1, 1)
1472  fi = 2*i + f_shift(1)
1473  v0 = coarse_coeffs(i, j, k)
1474  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1475  w_border0*wj*v0
1476  v1 = coarse_coeffs(i + 1, j, k)
1477  v2 = coarse_coeffs(i + 2, j, k)
1478  fi = fi + 1
1479  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1480  wj*(w_border1(1)*v0 + w_border1(2)*v1 + &
1481  w_border1(3)*v2)
1482  ELSE
1483  i = coarse_bo(1, 1)
1484  v0 = coarse_coeffs(i, j, k)
1485  v1 = coarse_coeffs(i + 1, j, k)
1486  v2 = coarse_coeffs(i + 2, j, k)
1487  fi = 2*i + f_shift(1) + 1
1488  IF (.NOT. (fi + 1 == fine_bo(1, 1) .OR. &
1489  fi + 2 == fine_bo(1, 1))) THEN
1490  CALL cp_abort(__location__, &
1491  "unexpected start index "// &
1492  trim(cp_to_string(coarse_bo(1, 1)))//" "// &
1493  trim(cp_to_string(fi)))
1494  END IF
1495  END IF
1496  fi = fi + 1
1497  IF (fi >= fine_bo(1, 1)) THEN
1498  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1499  ww0(1)*v0 + ww0(2)*v1 + &
1500  ww0(3)*v2
1501  ELSE
1502  cpassert(fi + 1 == fine_bo(1, 1))
1503  END IF
1504 ! CALL timestop(handle2)
1505 ! CALL timeset(routineN//"_core",handle2)
1506  DO i = coarse_bo(1, 1) + 3, floor((fine_bo(2, 1) - f_shift(1))/2.) - 3, 4
1507  v3 = coarse_coeffs(i, j, k)
1508  fi = fi + 1
1509  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1510  (ww1(1)*v0 + ww1(2)*v1 + &
1511  ww1(3)*v2 + ww1(4)*v3)
1512  fi = fi + 1
1513  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1514  (ww0(1)*v1 + ww0(2)*v2 + &
1515  ww0(3)*v3)
1516  v0 = coarse_coeffs(i + 1, j, k)
1517  fi = fi + 1
1518  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1519  (ww1(4)*v0 + ww1(1)*v1 + &
1520  ww1(2)*v2 + ww1(3)*v3)
1521  fi = fi + 1
1522  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1523  (ww0(1)*v2 + ww0(2)*v3 + &
1524  ww0(3)*v0)
1525  v1 = coarse_coeffs(i + 2, j, k)
1526  fi = fi + 1
1527  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1528  (ww1(3)*v0 + ww1(4)*v1 + &
1529  ww1(1)*v2 + ww1(2)*v3)
1530  fi = fi + 1
1531  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1532  (ww0(1)*v3 + ww0(2)*v0 + &
1533  ww0(3)*v1)
1534  v2 = coarse_coeffs(i + 3, j, k)
1535  fi = fi + 1
1536  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1537  (ww1(2)*v0 + ww1(3)*v1 + &
1538  ww1(4)*v2 + ww1(1)*v3)
1539  fi = fi + 1
1540  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1541  (ww0(1)*v0 + ww0(2)*v1 + &
1542  ww0(3)*v2)
1543  END DO
1544 ! CALL timestop(handle2)
1545 ! CALL timeset(routineN//"_clean",handle2)
1546  rest_b = modulo(floor((fine_bo(2, 1) - f_shift(1))/2.) - coarse_bo(1, 1) - 3 + 1, 4)
1547  IF (rest_b > 0) THEN
1548  i = floor((fine_bo(2, 1) - f_shift(1))/2.) - rest_b + 1
1549  v3 = coarse_coeffs(i, j, k)
1550  fi = fi + 1
1551  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1552  (ww1(1)*v0 + ww1(2)*v1 + &
1553  ww1(3)*v2 + ww1(4)*v3)
1554  fi = fi + 1
1555  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1556  (ww0(1)*v1 + ww0(2)*v2 + &
1557  ww0(3)*v3)
1558  IF (rest_b > 1) THEN
1559  v0 = coarse_coeffs(i + 1, j, k)
1560  fi = fi + 1
1561  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1562  (ww1(4)*v0 + ww1(1)*v1 + &
1563  ww1(2)*v2 + ww1(3)*v3)
1564  fi = fi + 1
1565  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1566  (ww0(1)*v2 + ww0(2)*v3 + &
1567  ww0(3)*v0)
1568  IF (rest_b > 2) THEN
1569  v1 = coarse_coeffs(i + 2, j, k)
1570  fi = fi + 1
1571  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1572  (ww1(3)*v0 + ww1(4)*v1 + &
1573  ww1(1)*v2 + ww1(2)*v3)
1574  fi = fi + 1
1575  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1576  (ww0(1)*v3 + ww0(2)*v0 + &
1577  ww0(3)*v1)
1578  IF (pbc .AND. .NOT. is_split) THEN
1579  v2 = coarse_coeffs(coarse_bo(1, 1), j, k)
1580  fi = fi + 1
1581  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1582  ww1(1)*v3 + ww1(2)*v0 + ww1(3)*v1 + ww1(4)*v2
1583  fi = fi + 1
1584  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1585  ww0(1)*v0 + ww0(2)*v1 + ww0(3)*v2
1586  v3 = coarse_coeffs(coarse_bo(1, 1) + 1, j, k)
1587  fi = fi + 1
1588  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1589  ww1(1)*v0 + ww1(2)*v1 + ww1(3)*v2 + ww1(4)*v3
1590  ELSE IF (has_i_ubound) THEN
1591  v2 = coarse_coeffs(i + 3, j, k)
1592  fi = fi + 1
1593  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1594  ww1(1)*v3 + ww1(2)*v0 + ww1(3)*v1 + ww1(4)*v2
1595  fi = fi + 1
1596  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1597  ww0(1)*v0 + ww0(2)*v1 + ww0(3)*v2
1598  IF (fi + 1 == fine_bo(2, 1)) THEN
1599  v3 = coarse_coeffs(i + 4, j, k)
1600  fi = fi + 1
1601  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1602  ww1(1)*v0 + ww1(2)*v1 + ww1(3)*v2 + ww1(4)*v3
1603  END IF
1604  ELSE
1605  fi = fi + 1
1606  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1607  wj*(w_border1(3)*v3 + w_border1(2)*v0 + &
1608  w_border1(1)*v1)
1609  fi = fi + 1
1610  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1611  w_border0*wj*v1
1612  END IF
1613  ELSE IF (pbc .AND. .NOT. is_split) THEN
1614  v1 = coarse_coeffs(coarse_bo(1, 1), j, k)
1615  fi = fi + 1
1616  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1617  ww1(1)*v2 + ww1(2)*v3 + ww1(3)*v0 + ww1(4)*v1
1618  fi = fi + 1
1619  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1620  ww0(1)*v3 + ww0(2)*v0 + ww0(3)*v1
1621  v2 = coarse_coeffs(coarse_bo(1, 1) + 1, j, k)
1622  fi = fi + 1
1623  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1624  ww1(1)*v3 + ww1(2)*v0 + ww1(3)*v1 + ww1(4)*v2
1625  ELSE IF (has_i_ubound) THEN
1626  v1 = coarse_coeffs(i + 2, j, k)
1627  fi = fi + 1
1628  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1629  ww1(1)*v2 + ww1(2)*v3 + ww1(3)*v0 + ww1(4)*v1
1630  fi = fi + 1
1631  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1632  ww0(1)*v3 + ww0(2)*v0 + ww0(3)*v1
1633  IF (fi + 1 == fine_bo(2, 1)) THEN
1634  v2 = coarse_coeffs(i + 3, j, k)
1635  fi = fi + 1
1636  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1637  ww1(1)*v3 + ww1(2)*v0 + ww1(3)*v1 + ww1(4)*v2
1638  END IF
1639  ELSE
1640  fi = fi + 1
1641  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1642  wj*(w_border1(3)*v2 + w_border1(2)*v3 + &
1643  w_border1(1)*v0)
1644  fi = fi + 1
1645  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1646  w_border0*wj*v0
1647  END IF
1648  ELSE IF (pbc .AND. .NOT. is_split) THEN
1649  v0 = coarse_coeffs(coarse_bo(1, 1), j, k)
1650  fi = fi + 1
1651  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1652  ww1(1)*v1 + ww1(2)*v2 + ww1(3)*v3 + ww1(4)*v0
1653  fi = fi + 1
1654  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1655  ww0(1)*v2 + ww0(2)*v3 + ww0(3)*v0
1656  v1 = coarse_coeffs(coarse_bo(1, 1) + 1, j, k)
1657  fi = fi + 1
1658  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1659  ww1(1)*v2 + ww1(2)*v3 + ww1(3)*v0 + ww1(4)*v1
1660  ELSE IF (has_i_ubound) THEN
1661  v0 = coarse_coeffs(i + 1, j, k)
1662  fi = fi + 1
1663  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1664  ww1(1)*v1 + ww1(2)*v2 + ww1(3)*v3 + ww1(4)*v0
1665  fi = fi + 1
1666  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1667  ww0(1)*v2 + ww0(2)*v3 + ww0(3)*v0
1668  IF (fi + 1 == fine_bo(2, 1)) THEN
1669  v1 = coarse_coeffs(i + 2, j, k)
1670  fi = fi + 1
1671  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1672  ww1(1)*v2 + ww1(2)*v3 + ww1(3)*v0 + ww1(4)*v1
1673  END IF
1674  ELSE
1675  fi = fi + 1
1676  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1677  wj*(w_border1(3)*v1 + w_border1(2)*v2 + &
1678  w_border1(1)*v3)
1679  fi = fi + 1
1680  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1681  w_border0*wj*v3
1682  END IF
1683  ELSE IF (pbc .AND. .NOT. is_split) THEN
1684  v3 = coarse_coeffs(coarse_bo(1, 1), j, k)
1685  fi = fi + 1
1686  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1687  ww1(1)*v0 + ww1(2)*v1 + ww1(3)*v2 + ww1(4)*v3
1688  fi = fi + 1
1689  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1690  ww0(1)*v1 + ww0(2)*v2 + ww0(3)*v3
1691  v0 = coarse_coeffs(coarse_bo(1, 1) + 1, j, k)
1692  fi = fi + 1
1693  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1694  ww1(1)*v1 + ww1(2)*v2 + ww1(3)*v3 + ww1(4)*v0
1695  ELSE IF (has_i_ubound) THEN
1696  v3 = coarse_coeffs(i, j, k)
1697  fi = fi + 1
1698  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1699  ww1(1)*v0 + ww1(2)*v1 + ww1(3)*v2 + ww1(4)*v3
1700  fi = fi + 1
1701  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1702  ww0(1)*v1 + ww0(2)*v2 + ww0(3)*v3
1703  IF (fi + 1 == fine_bo(2, 1)) THEN
1704  v0 = coarse_coeffs(i + 1, j, k)
1705  fi = fi + 1
1706  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1707  ww1(1)*v1 + ww1(2)*v2 + ww1(3)*v3 + ww1(4)*v0
1708  END IF
1709  ELSE
1710  fi = fi + 1
1711  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1712  wj*(w_border1(3)*v0 + w_border1(2)*v1 + &
1713  w_border1(1)*v2)
1714  fi = fi + 1
1715  fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1716  w_border0*wj*v2
1717  END IF
1718  cpassert(fi == fine_bo(2, 1))
1719  END IF
1720 ! CALL timestop(handle2)
1721  END DO
1722  END DO
1723  END DO
1724  END DO
1725 
1726  IF (is_split) THEN
1727  DEALLOCATE (coarse_coeffs)
1728  END IF
1729  CALL timestop(handle)
1730  END SUBROUTINE add_coarse2fine
1731 
1732 ! **************************************************************************************************
1733 !> \brief low level function that adds a coarse grid (without boundary)
1734 !> to a fine grid.
1735 !>
1736 !> It will add to
1737 !>
1738 !> coarse_coeffs(coarse_bounds(1,1):coarse_bounds(2,1),
1739 !> coarse_bounds(1,2):coarse_bounds(2,2),
1740 !> coarse_bounds(1,3):coarse_bounds(2,3))
1741 !>
1742 !> using
1743 !>
1744 !> fine_values(2*coarse_bounds(1,1):2*coarse_bounds(2,1),
1745 !> 2*coarse_bounds(1,2):2*coarse_bounds(2,2),
1746 !> 2*coarse_bounds(1,3):2*coarse_bounds(2,3))
1747 !>
1748 !> composed with the weights obtained by the direct product of the
1749 !> 1d coefficients weights:
1750 !>
1751 !> for i,j,k in -3..3
1752 !> w(i,j,k)=weights_1d(abs(i)+1)*weights_1d(abs(j)+1)*
1753 !> weights_1d(abs(k)+1)
1754 !> \param fine_values_pw 3d array where to add the values due to the
1755 !> coarse coeffs
1756 !> \param coarse_coeffs_pw 3d array with boundary of size 1 with the values of the
1757 !> coefficients
1758 !> \param weights_1d the weights of the 1d smearing
1759 !> \param w_border0 the 1d weight at the border
1760 !> \param w_border1 the 1d weights for a point one off the border
1761 !> (w_border1(1) is the weight of the coefficent at the border)
1762 !> \param pbc ...
1763 !> \param safe_computation ...
1764 !> \author fawzi
1765 !> \note
1766 !> see coarse2fine for some relevant notes
1767 ! **************************************************************************************************
1768  SUBROUTINE add_fine2coarse(fine_values_pw, coarse_coeffs_pw, &
1769  weights_1d, w_border0, w_border1, pbc, safe_computation)
1770  TYPE(pw_r3d_rs_type), INTENT(IN) :: fine_values_pw, coarse_coeffs_pw
1771  REAL(kind=dp), DIMENSION(4), INTENT(in) :: weights_1d
1772  REAL(kind=dp), INTENT(in) :: w_border0
1773  REAL(kind=dp), DIMENSION(3), INTENT(in) :: w_border1
1774  LOGICAL, INTENT(in) :: pbc
1775  LOGICAL, INTENT(in), OPTIONAL :: safe_computation
1776 
1777  CHARACTER(len=*), PARAMETER :: routinen = 'add_fine2coarse'
1778 
1779  INTEGER :: coarse_slice_size, f_shift(3), fi, fj, fk, handle, handle2, i, ii, ij, ik, ip, j, &
1780  k, n_procs, p, p_old, rcv_tot_size, rest_b, s(3), send_tot_size, ss, x, x_att
1781  INTEGER, ALLOCATABLE, DIMENSION(:) :: pp_lb, pp_ub, rcv_offset, rcv_size, &
1782  real_rcv_size, send_offset, send_size, &
1783  sent_size
1784  INTEGER, DIMENSION(2, 3) :: coarse_bo, coarse_gbo, fine_bo, &
1785  fine_gbo, my_coarse_bo
1786  INTEGER, DIMENSION(:), POINTER :: pos_of_x
1787  LOGICAL :: has_i_lbound, has_i_ubound, is_split, &
1788  local_data, safe_calc
1789  REAL(kind=dp) :: vv0, vv1, vv2, vv3, vv4, vv5, vv6, vv7, &
1790  wi, wj, wk
1791  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rcv_buf, send_buf
1792  REAL(kind=dp), DIMENSION(3) :: w_0, ww0
1793  REAL(kind=dp), DIMENSION(4) :: w_1, ww1
1794  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: coarse_coeffs, fine_values
1795 
1796  CALL timeset(routinen, handle)
1797 
1798  safe_calc = .false.
1799  IF (PRESENT(safe_computation)) safe_calc = safe_computation
1800 
1801  my_coarse_bo = coarse_coeffs_pw%pw_grid%bounds_local
1802  coarse_gbo = coarse_coeffs_pw%pw_grid%bounds
1803  fine_bo = fine_values_pw%pw_grid%bounds_local
1804  fine_gbo = fine_values_pw%pw_grid%bounds
1805  f_shift = fine_gbo(1, :) - 2*coarse_gbo(1, :)
1806  is_split = any(coarse_gbo(:, 1) /= my_coarse_bo(:, 1))
1807  coarse_bo = my_coarse_bo
1808  IF (fine_bo(1, 1) <= fine_bo(2, 1)) THEN
1809  coarse_bo(1, 1) = floor(real(fine_bo(1, 1) - f_shift(1), dp)/2._dp) - 1
1810  coarse_bo(2, 1) = floor(real(fine_bo(2, 1) + 1 - f_shift(1), dp)/2._dp) + 1
1811  ELSE
1812  coarse_bo(1, 1) = coarse_gbo(2, 1)
1813  coarse_bo(2, 1) = coarse_gbo(2, 1) - 1
1814  END IF
1815  IF (.NOT. is_split .OR. .NOT. pbc) THEN
1816  coarse_bo(1, 1) = max(coarse_gbo(1, 1), coarse_bo(1, 1))
1817  coarse_bo(2, 1) = min(coarse_gbo(2, 1), coarse_bo(2, 1))
1818  END IF
1819  has_i_ubound = (fine_gbo(2, 1) /= fine_bo(2, 1)) .OR. pbc .AND. is_split
1820  has_i_lbound = (fine_gbo(1, 1) /= fine_bo(1, 1)) .OR. pbc .AND. is_split
1821 
1822  IF (pbc) THEN
1823  cpassert(all(fine_gbo(1, :) == 2*coarse_gbo(1, :) + f_shift))
1824  cpassert(all(fine_gbo(2, :) == 2*coarse_gbo(2, :) + f_shift + 1))
1825  ELSE
1826  cpassert(all(fine_gbo(2, :) == 2*coarse_gbo(2, :) + f_shift))
1827  cpassert(all(fine_gbo(1, :) == 2*coarse_gbo(1, :) + f_shift))
1828  END IF
1829  cpassert(coarse_gbo(2, 1) - coarse_gbo(1, 2) > 1)
1830  local_data = is_split ! ANY(coarse_bo/=my_coarse_bo)
1831  IF (local_data) THEN
1832  ALLOCATE (coarse_coeffs(coarse_bo(1, 1):coarse_bo(2, 1), &
1833  coarse_bo(1, 2):coarse_bo(2, 2), &
1834  coarse_bo(1, 3):coarse_bo(2, 3)))
1835  coarse_coeffs = 0._dp
1836  ELSE
1837  coarse_coeffs => coarse_coeffs_pw%array
1838  END IF
1839 
1840  fine_values => fine_values_pw%array
1841  w_0 = (/weights_1d(3), weights_1d(1), weights_1d(3)/)
1842  w_1 = (/weights_1d(4), weights_1d(2), weights_1d(2), weights_1d(4)/)
1843 
1844  DO i = 1, 3
1845  s(i) = coarse_gbo(2, i) - coarse_gbo(1, i) + 1
1846  END DO
1847  IF (any(s < 1)) RETURN
1848 
1849  DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1850  DO ik = -3, 3
1851  IF (pbc) THEN
1852  wk = weights_1d(abs(ik) + 1)
1853  fk = fine_gbo(1, 3) + modulo(2*k + ik - fine_gbo(1, 3) + f_shift(3), 2*s(3))
1854  ELSE
1855  fk = 2*k + ik + f_shift(3)
1856  IF (fk <= fine_bo(1, 3) + 1 .OR. fk >= fine_bo(2, 3) - 1) THEN
1857  IF (fk < fine_bo(1, 3) .OR. fk > fine_bo(2, 3)) cycle
1858  IF (fk == fine_bo(1, 3) .OR. fk == fine_bo(2, 3)) THEN
1859  IF (ik /= 0) cycle
1860  wk = w_border0
1861  ELSE IF (fk == fine_bo(1, 3) + 1) THEN
1862  SELECT CASE (ik)
1863  CASE (1)
1864  wk = w_border1(1)
1865  CASE (-1)
1866  wk = w_border1(2)
1867  CASE (-3)
1868  wk = w_border1(3)
1869  CASE default
1870  cpabort("")
1871  cycle
1872  END SELECT
1873  ELSE
1874  SELECT CASE (ik)
1875  CASE (3)
1876  wk = w_border1(3)
1877  CASE (1)
1878  wk = w_border1(2)
1879  CASE (-1)
1880  wk = w_border1(1)
1881  CASE default
1882  cpabort("")
1883  cycle
1884  END SELECT
1885  END IF
1886  ELSE
1887  wk = weights_1d(abs(ik) + 1)
1888  END IF
1889  END IF
1890  DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1891  DO ij = -3, 3
1892  IF (pbc) THEN
1893  fj = fine_gbo(1, 2) + modulo(2*j + ij - fine_gbo(1, 2) + f_shift(2), &
1894  2*s(2))
1895  wj = weights_1d(abs(ij) + 1)*wk
1896  ELSE
1897  fj = 2*j + ij + f_shift(2)
1898  IF (fj <= fine_bo(1, 2) + 1 .OR. fj >= fine_bo(2, 2) - 1) THEN
1899  IF (fj < fine_bo(1, 2) .OR. fj > fine_bo(2, 2)) cycle
1900  IF (fj == fine_bo(1, 2) .OR. fj == fine_bo(2, 2)) THEN
1901  IF (ij /= 0) cycle
1902  wj = w_border0*wk
1903  ELSE IF (fj == fine_bo(1, 2) + 1) THEN
1904  SELECT CASE (ij)
1905  CASE (1)
1906  wj = w_border1(1)*wk
1907  CASE (-1)
1908  wj = w_border1(2)*wk
1909  CASE (-3)
1910  wj = w_border1(3)*wk
1911  CASE default
1912  cpabort("")
1913  cycle
1914  END SELECT
1915  ELSE
1916  SELECT CASE (ij)
1917  CASE (-1)
1918  wj = w_border1(1)*wk
1919  CASE (1)
1920  wj = w_border1(2)*wk
1921  CASE (3)
1922  wj = w_border1(3)*wk
1923  CASE default
1924  cpabort("")
1925  cycle
1926  END SELECT
1927  END IF
1928  ELSE
1929  wj = weights_1d(abs(ij) + 1)*wk
1930  END IF
1931  END IF
1932 
1933  IF (coarse_bo(2, 1) - coarse_bo(1, 1) < 7 .OR. safe_calc) THEN
1934  DO i = coarse_bo(1, 1), coarse_bo(2, 1)
1935  DO ii = -3, 3
1936  IF (pbc .AND. .NOT. is_split) THEN
1937  wi = weights_1d(abs(ii) + 1)*wj
1938  fi = fine_gbo(1, 1) + modulo(2*i + ii - fine_gbo(1, 1) + f_shift(1), 2*s(1))
1939  ELSE
1940  fi = 2*i + ii + f_shift(1)
1941  IF (fi < fine_bo(1, 1) .OR. fi > fine_bo(2, 1)) cycle
1942  IF (((.NOT. pbc) .AND. fi <= fine_gbo(1, 1) + 1) .OR. &
1943  ((.NOT. pbc) .AND. fi >= fine_gbo(2, 1) - 1)) THEN
1944  IF (fi == fine_gbo(1, 1) .OR. fi == fine_gbo(2, 1)) THEN
1945  IF (ii /= 0) cycle
1946  wi = w_border0*wj
1947  ELSE IF (fi == fine_gbo(1, 1) + 1) THEN
1948  SELECT CASE (ii)
1949  CASE (1)
1950  wi = w_border1(1)*wj
1951  CASE (-1)
1952  wi = w_border1(2)*wj
1953  CASE (-3)
1954  wi = w_border1(3)*wj
1955  CASE default
1956  cycle
1957  END SELECT
1958  ELSE
1959  SELECT CASE (ii)
1960  CASE (-1)
1961  wi = w_border1(1)*wj
1962  CASE (1)
1963  wi = w_border1(2)*wj
1964  CASE (3)
1965  wi = w_border1(3)*wj
1966  CASE default
1967  cycle
1968  END SELECT
1969  END IF
1970  ELSE
1971  wi = weights_1d(abs(ii) + 1)*wj
1972  END IF
1973  END IF
1974  coarse_coeffs(i, j, k) = &
1975  coarse_coeffs(i, j, k) + &
1976  wi*fine_values(fi, fj, fk)
1977  END DO
1978  END DO
1979  ELSE
1980  ww0 = wj*w_0
1981  ww1 = wj*w_1
1982  IF (pbc .AND. .NOT. is_split) THEN
1983  i = coarse_bo(1, 1) - 1
1984  vv2 = fine_values(fine_bo(2, 1) - 2, fj, fk)
1985  vv3 = fine_values(fine_bo(2, 1) - 1, fj, fk)
1986  vv4 = fine_values(fine_bo(2, 1), fj, fk)
1987  fi = fine_bo(1, 1)
1988  vv5 = fine_values(fi, fj, fk)
1989  fi = fi + 1
1990  vv6 = fine_values(fi, fj, fk)
1991  fi = fi + 1
1992  vv7 = fine_values(fi, fj, fk)
1993  coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
1994  + ww1(4)*vv2 + ww0(3)*vv3 + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6 + ww0(1)*vv7
1995  coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
1996  + ww1(4)*vv4 + ww0(3)*vv5 + ww1(3)*vv6 + ww0(2)*vv7
1997  coarse_coeffs(i + 3, j, k) = coarse_coeffs(i + 3, j, k) &
1998  + ww1(4)*vv6 + ww0(3)*vv7
1999  ELSE IF (has_i_lbound) THEN
2000  i = coarse_bo(1, 1)
2001  fi = fine_bo(1, 1) - 1
2002  IF (i + 1 == floor((fine_bo(1, 1) + 1 - f_shift(1))/2._dp)) THEN
2003  fi = fi + 1
2004  vv0 = fine_values(fi, fj, fk)
2005  coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) + &
2006  vv0*ww0(3)
2007  coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) + &
2008  vv0*ww0(2)
2009  coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) + &
2010  vv0*ww0(1)
2011  END IF
2012  ELSE
2013  i = coarse_bo(1, 1)
2014  fi = 2*i + f_shift(1)
2015  vv0 = fine_values(fi, fj, fk)
2016  fi = fi + 1
2017  vv1 = fine_values(fi, fj, fk)
2018  fi = fi + 1
2019  vv2 = fine_values(fi, fj, fk)
2020  coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) + &
2021  (vv0*w_border0 + vv1*w_border1(1))*wj + vv2*ww0(1)
2022  coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) + &
2023  wj*w_border1(2)*vv1 + ww0(2)*vv2
2024  coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) + &
2025  wj*w_border1(3)*vv1 + ww0(3)*vv2
2026  END IF
2027  DO i = coarse_bo(1, 1) + 3, floor((fine_bo(2, 1) - f_shift(1))/2._dp) - 3, 4
2028  fi = fi + 1
2029  vv0 = fine_values(fi, fj, fk)
2030  fi = fi + 1
2031  vv1 = fine_values(fi, fj, fk)
2032  fi = fi + 1
2033  vv2 = fine_values(fi, fj, fk)
2034  fi = fi + 1
2035  vv3 = fine_values(fi, fj, fk)
2036  fi = fi + 1
2037  vv4 = fine_values(fi, fj, fk)
2038  fi = fi + 1
2039  vv5 = fine_values(fi, fj, fk)
2040  fi = fi + 1
2041  vv6 = fine_values(fi, fj, fk)
2042  fi = fi + 1
2043  vv7 = fine_values(fi, fj, fk)
2044  coarse_coeffs(i - 3, j, k) = coarse_coeffs(i - 3, j, k) &
2045  + ww1(1)*vv0
2046  coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2047  + ww1(2)*vv0 + ww0(1)*vv1 + ww1(1)*vv2
2048  coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2049  + ww1(3)*vv0 + ww0(2)*vv1 + ww1(2)*vv2 + ww0(1)*vv3 + ww1(1)*vv4
2050  coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2051  + ww1(4)*vv0 + ww0(3)*vv1 + ww1(3)*vv2 + ww0(2)*vv3 + ww1(2)*vv4 + ww0(1)*vv5 + ww1(1)*vv6
2052  coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2053  + ww1(4)*vv2 + ww0(3)*vv3 + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6 + ww0(1)*vv7
2054  coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2055  + ww1(4)*vv4 + ww0(3)*vv5 + ww1(3)*vv6 + ww0(2)*vv7
2056  coarse_coeffs(i + 3, j, k) = coarse_coeffs(i + 3, j, k) &
2057  + ww1(4)*vv6 + ww0(3)*vv7
2058  END DO
2059  IF (.NOT. floor((fine_bo(2, 1) - f_shift(1))/2._dp) - coarse_bo(1, 1) >= 4) THEN
2060  cpabort("FLOOR((fine_bo(2,1)-f_shift(1))/2._dp)-coarse_bo(1,1)>=4")
2061  END IF
2062  rest_b = modulo(floor((fine_bo(2, 1) - f_shift(1))/2._dp) - coarse_bo(1, 1) - 6, 4)
2063  i = floor((fine_bo(2, 1) - f_shift(1))/2._dp) - 3 - rest_b + 4
2064  cpassert(fi == (i - 2)*2 + f_shift(1))
2065  IF (rest_b > 0) THEN
2066  fi = fi + 1
2067  vv0 = fine_values(fi, fj, fk)
2068  fi = fi + 1
2069  vv1 = fine_values(fi, fj, fk)
2070  coarse_coeffs(i - 3, j, k) = coarse_coeffs(i - 3, j, k) &
2071  + ww1(1)*vv0
2072  coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2073  + ww1(2)*vv0 + ww0(1)*vv1
2074  coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2075  + ww1(3)*vv0 + ww0(2)*vv1
2076  coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2077  + ww1(4)*vv0 + ww0(3)*vv1
2078  IF (rest_b > 1) THEN
2079  fi = fi + 1
2080  vv2 = fine_values(fi, fj, fk)
2081  fi = fi + 1
2082  vv3 = fine_values(fi, fj, fk)
2083  coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2084  + ww1(1)*vv2
2085  coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2086  + ww1(2)*vv2 + ww0(1)*vv3
2087  coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2088  + ww1(3)*vv2 + ww0(2)*vv3
2089  coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2090  + ww1(4)*vv2 + ww0(3)*vv3
2091  IF (rest_b > 2) THEN
2092  fi = fi + 1
2093  vv4 = fine_values(fi, fj, fk)
2094  fi = fi + 1
2095  vv5 = fine_values(fi, fj, fk)
2096  fi = fi + 1
2097  vv6 = fine_values(fi, fj, fk)
2098  fi = fi + 1
2099  vv7 = fine_values(fi, fj, fk)
2100  coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2101  + ww1(1)*vv4
2102  IF (has_i_ubound) THEN
2103  IF (coarse_bo(2, 1) - 2 == floor((fine_bo(2, 1) - f_shift(1))/2._dp)) THEN
2104  fi = fi + 1
2105  vv0 = fine_values(fi, fj, fk)
2106  coarse_coeffs(i + 4, j, k) = coarse_coeffs(i + 4, j, k) &
2107  + vv0*ww1(4)
2108  ELSE
2109  vv0 = 0._dp
2110  END IF
2111  coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2112  + ww1(2)*vv4 + ww0(1)*vv5 + ww1(1)*vv6
2113  coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2114  + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6 + ww0(1)*vv7 + vv0*ww1(1)
2115  coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2116  + ww1(4)*vv4 + ww0(3)*vv5 + ww1(3)*vv6 + ww0(2)*vv7 + vv0*ww1(2)
2117  coarse_coeffs(i + 3, j, k) = coarse_coeffs(i + 3, j, k) &
2118  + ww1(4)*vv6 + ww0(3)*vv7 + vv0*ww1(3)
2119  ELSEIF (pbc .AND. .NOT. is_split) THEN
2120  fi = fi + 1
2121  vv0 = fine_values(fi, fj, fk)
2122  vv1 = fine_values(fine_bo(1, 1), fj, fk)
2123  vv2 = fine_values(fine_bo(1, 1) + 1, fj, fk)
2124  coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2125  + ww1(2)*vv4 + ww0(1)*vv5 + ww1(1)*vv6
2126  coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2127  + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6 + ww0(1)*vv7 + vv0*ww1(1)
2128  coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2129  + ww1(4)*vv4 + ww0(3)*vv5 + ww1(3)*vv6 + ww0(2)*vv7 + vv0*ww1(2) &
2130  + vv1*ww0(1) + vv2*ww1(1)
2131  ELSE
2132  coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2133  + ww1(2)*vv4 + ww0(1)*vv5 + wj*w_border1(3)*vv6
2134  coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2135  + ww1(3)*vv4 + ww0(2)*vv5 + wj*w_border1(2)*vv6
2136  coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2137  + ww1(4)*vv4 + ww0(3)*vv5 + wj*w_border1(1)*vv6 + w_border0*wj*vv7
2138  END IF
2139  ELSE
2140  fi = fi + 1
2141  vv4 = fine_values(fi, fj, fk)
2142  fi = fi + 1
2143  vv5 = fine_values(fi, fj, fk)
2144  IF (has_i_ubound) THEN
2145  IF (coarse_bo(2, 1) - 2 == floor((fine_bo(2, 1) - f_shift(1))/2._dp)) THEN
2146  fi = fi + 1
2147  vv6 = fine_values(fi, fj, fk)
2148  coarse_coeffs(i + 3, j, k) = coarse_coeffs(i + 3, j, k) &
2149  + ww1(4)*vv6
2150  ELSE
2151  vv6 = 0._dp
2152  END IF
2153  coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2154  + ww1(1)*vv4
2155  coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2156  + ww1(2)*vv4 + ww0(1)*vv5 + ww1(1)*vv6
2157  coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2158  + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6
2159  coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2160  + ww1(4)*vv4 + ww0(3)*vv5 + ww1(3)*vv6
2161  ELSEIF (pbc .AND. .NOT. is_split) THEN
2162  fi = fi + 1
2163  vv6 = fine_values(fi, fj, fk)
2164  vv7 = fine_values(fine_bo(1, 1), fj, fk)
2165  vv0 = fine_values(fine_bo(1, 1) + 1, fj, fk)
2166  coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2167  + ww1(1)*vv4
2168  coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2169  + ww1(4)*vv0 + ww0(3)*vv1 + ww1(3)*vv2 + ww0(2)*vv3 + &
2170  ww1(2)*vv4 + ww0(1)*vv5 + ww1(1)*vv6
2171  coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2172  + ww1(4)*vv2 + ww0(3)*vv3 + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6 &
2173  + ww0(1)*vv7 + ww1(1)*vv0
2174  ELSE
2175  coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2176  + wj*w_border1(3)*vv4
2177  coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2178  + wj*w_border1(2)*vv4
2179  coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2180  + wj*(w_border1(1)*vv4 + w_border0*vv5)
2181  END IF
2182  END IF
2183  ELSE
2184  fi = fi + 1
2185  vv2 = fine_values(fi, fj, fk)
2186  fi = fi + 1
2187  vv3 = fine_values(fi, fj, fk)
2188  IF (has_i_ubound) THEN
2189  IF (coarse_bo(2, 1) - 2 == floor((fine_bo(2, 1) - f_shift(1))/2._dp)) THEN
2190  fi = fi + 1
2191  vv4 = fine_values(fi, fj, fk)
2192  coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2193  + ww1(4)*vv4
2194  ELSE
2195  vv4 = 0._dp
2196  END IF
2197  coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2198  + ww1(1)*vv2
2199  coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2200  + ww1(2)*vv2 + ww0(1)*vv3 + ww1(1)*vv4
2201  coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2202  + ww1(3)*vv2 + ww0(2)*vv3 + ww1(2)*vv4
2203  coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2204  + ww1(4)*vv2 + ww0(3)*vv3 + ww1(3)*vv4
2205  ELSEIF (pbc .AND. .NOT. is_split) THEN
2206  fi = fi + 1
2207  vv4 = fine_values(fi, fj, fk)
2208  vv5 = fine_values(fine_bo(1, 1), fj, fk)
2209  vv6 = fine_values(fine_bo(1, 1) + 1, fj, fk)
2210  coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2211  + ww1(1)*vv2
2212  coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2213  + ww1(2)*vv2 + ww0(1)*vv3 + ww1(1)*vv4
2214  coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2215  + ww1(3)*vv2 + ww0(2)*vv3 + ww1(2)*vv4 + vv5*ww0(1) + ww1(1)*vv6
2216  ELSE
2217  coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2218  + wj*w_border1(3)*vv2
2219  coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2220  + wj*w_border1(2)*vv2
2221  coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2222  + wj*(w_border1(1)*vv2 + w_border0*vv3)
2223  END IF
2224  END IF
2225  ELSE
2226  fi = fi + 1
2227  vv0 = fine_values(fi, fj, fk)
2228  fi = fi + 1
2229  vv1 = fine_values(fi, fj, fk)
2230  IF (has_i_ubound) THEN
2231  IF (coarse_bo(2, 1) - 2 == floor((fine_bo(2, 1) - f_shift(1))/2._dp)) THEN
2232  fi = fi + 1
2233  vv2 = fine_values(fi, fj, fk)
2234  coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2235  + ww1(4)*vv2
2236  ELSE
2237  vv2 = 0._dp
2238  END IF
2239  coarse_coeffs(i - 3, j, k) = coarse_coeffs(i - 3, j, k) &
2240  + ww1(1)*vv0
2241  coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2242  + ww1(2)*vv0 + ww0(1)*vv1 + ww1(1)*vv2
2243  coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2244  + ww1(3)*vv0 + ww0(2)*vv1 + ww1(2)*vv2
2245  coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2246  + ww1(4)*vv0 + ww0(3)*vv1 + ww1(3)*vv2
2247  ELSEIF (pbc .AND. .NOT. is_split) THEN
2248  fi = fi + 1
2249  vv2 = fine_values(fi, fj, fk)
2250  vv3 = fine_values(fine_bo(1, 1), fk, fk)
2251  vv4 = fine_values(fine_bo(1, 1) + 1, fk, fk)
2252  coarse_coeffs(i - 3, j, k) = coarse_coeffs(i - 3, j, k) &
2253  + ww1(1)*vv0
2254  coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2255  + ww1(2)*vv0 + ww0(1)*vv1 + ww1(1)*vv2
2256  coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2257  + ww1(3)*vv0 + ww0(2)*vv1 + ww1(2)*vv2 + ww0(1)*vv3 + ww1(1)*vv4
2258  ELSE
2259  coarse_coeffs(i - 3, j, k) = coarse_coeffs(i - 3, j, k) &
2260  + wj*w_border1(3)*vv0
2261  coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2262  + wj*w_border1(2)*vv0
2263  coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2264  + wj*(w_border1(1)*vv0 + w_border0*vv1)
2265  END IF
2266  END IF
2267  cpassert(fi == fine_bo(2, 1))
2268  END IF
2269  END DO
2270  END DO
2271  END DO
2272  END DO
2273 
2274  ! *** parallel case
2275  IF (is_split) THEN
2276  CALL timeset(routinen//"_comm", handle2)
2277  coarse_slice_size = (coarse_bo(2, 2) - coarse_bo(1, 2) + 1)* &
2278  (coarse_bo(2, 3) - coarse_bo(1, 3) + 1)
2279  n_procs = coarse_coeffs_pw%pw_grid%para%group_size
2280  ALLOCATE (send_size(0:n_procs - 1), send_offset(0:n_procs - 1), &
2281  sent_size(0:n_procs - 1), rcv_size(0:n_procs - 1), &
2282  rcv_offset(0:n_procs - 1), pp_lb(0:n_procs - 1), &
2283  pp_ub(0:n_procs - 1), real_rcv_size(0:n_procs - 1))
2284 
2285  ! ** send size count
2286 
2287  pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x
2288  send_size = 0
2289  DO x = coarse_bo(1, 1), coarse_bo(2, 1)
2290  p = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
2291  send_size(p) = send_size(p) + coarse_slice_size
2292  END DO
2293 
2294  ! ** rcv size count
2295 
2296  pos_of_x => fine_values_pw%pw_grid%para%pos_of_x
2297  p_old = pos_of_x(fine_gbo(1, 1))
2298  pp_lb = fine_gbo(2, 1)
2299  pp_ub = fine_gbo(2, 1) - 1
2300  pp_lb(p_old) = fine_gbo(1, 1)
2301  DO x = fine_gbo(1, 1), fine_gbo(2, 1)
2302  p = pos_of_x(x)
2303  IF (p /= p_old) THEN
2304  pp_ub(p_old) = x - 1
2305  pp_lb(p) = x
2306  p_old = p
2307  END IF
2308  END DO
2309  pp_ub(p_old) = fine_gbo(2, 1)
2310 
2311  DO ip = 0, n_procs - 1
2312  IF (pp_lb(ip) <= pp_ub(ip)) THEN
2313  pp_lb(ip) = floor(real(pp_lb(ip) - f_shift(1), dp)/2._dp) - 1
2314  pp_ub(ip) = floor(real(pp_ub(ip) + 1 - f_shift(1), dp)/2._dp) + 1
2315  ELSE
2316  pp_lb(ip) = coarse_gbo(2, 1)
2317  pp_ub(ip) = coarse_gbo(2, 1) - 1
2318  END IF
2319  IF (.NOT. is_split .OR. .NOT. pbc) THEN
2320  pp_lb(ip) = max(pp_lb(ip), coarse_gbo(1, 1))
2321  pp_ub(ip) = min(pp_ub(ip), coarse_gbo(2, 1))
2322  END IF
2323  END DO
2324 
2325  rcv_size = 0
2326  DO ip = 0, n_procs - 1
2327  DO x = pp_lb(ip), coarse_gbo(1, 1) - 1
2328  x_att = coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1))
2329  IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
2330  rcv_size(ip) = rcv_size(ip) + coarse_slice_size
2331  END IF
2332  END DO
2333  rcv_size(ip) = rcv_size(ip) + coarse_slice_size* &
2334  max(0, &
2335  min(pp_ub(ip), my_coarse_bo(2, 1)) - max(pp_lb(ip), my_coarse_bo(1, 1)) + 1)
2336  DO x = coarse_gbo(2, 1) + 1, pp_ub(ip)
2337  x_att = coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1))
2338  IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
2339  rcv_size(ip) = rcv_size(ip) + coarse_slice_size
2340  END IF
2341  END DO
2342  END DO
2343 
2344  ! ** offsets & alloc send-rcv
2345 
2346  send_tot_size = 0
2347  DO ip = 0, n_procs - 1
2348  send_offset(ip) = send_tot_size
2349  send_tot_size = send_tot_size + send_size(ip)
2350  END DO
2351  IF (send_tot_size /= (coarse_bo(2, 1) - coarse_bo(1, 1) + 1)*coarse_slice_size) &
2352  cpabort("Error calculating send_tot_size")
2353  ALLOCATE (send_buf(0:send_tot_size - 1))
2354 
2355  rcv_tot_size = 0
2356  DO ip = 0, n_procs - 1
2357  rcv_offset(ip) = rcv_tot_size
2358  rcv_tot_size = rcv_tot_size + rcv_size(ip)
2359  END DO
2360  ALLOCATE (rcv_buf(0:rcv_tot_size - 1))
2361 
2362  ! ** fill send buffer
2363 
2364  pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x
2365  p_old = pos_of_x(coarse_gbo(1, 1) &
2366  + modulo(coarse_bo(1, 1) - coarse_gbo(1, 1), s(1)))
2367  sent_size(:) = send_offset
2368  ss = coarse_bo(2, 1) - coarse_bo(1, 1) + 1
2369  DO x = coarse_bo(1, 1), coarse_bo(2, 1)
2370  p = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
2371  CALL dcopy(coarse_slice_size, &
2372  coarse_coeffs(x, coarse_bo(1, 2), &
2373  coarse_bo(1, 3)), ss, send_buf(sent_size(p)), 1)
2374  sent_size(p) = sent_size(p) + coarse_slice_size
2375  END DO
2376 
2377  IF (any(sent_size(0:n_procs - 2) /= send_offset(1:n_procs - 1))) &
2378  cpabort("error 1 filling send buffer")
2379  IF (sent_size(n_procs - 1) /= send_tot_size) &
2380  cpabort("error 2 filling send buffer")
2381 
2382  IF (local_data) THEN
2383  DEALLOCATE (coarse_coeffs)
2384  ELSE
2385  NULLIFY (coarse_coeffs)
2386  END IF
2387 
2388  cpassert(all(sent_size(:n_procs - 2) == send_offset(1:)))
2389  cpassert(sent_size(n_procs - 1) == send_tot_size)
2390  ! test send/rcv sizes
2391  CALL coarse_coeffs_pw%pw_grid%para%group%alltoall(send_size, real_rcv_size, 1)
2392 
2393  cpassert(all(real_rcv_size == rcv_size))
2394  ! all2all
2395  CALL coarse_coeffs_pw%pw_grid%para%group%alltoall(sb=send_buf, scount=send_size, sdispl=send_offset, &
2396  rb=rcv_buf, rcount=rcv_size, rdispl=rcv_offset)
2397 
2398  ! ** sum & reorder rcv buffer
2399 
2400  sent_size(:) = rcv_offset
2401  DO ip = 0, n_procs - 1
2402 
2403  DO x = pp_lb(ip), coarse_gbo(1, 1) - 1
2404  x_att = coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1))
2405  IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
2406  ii = sent_size(ip)
2407  DO k = coarse_bo(1, 3), coarse_bo(2, 3)
2408  DO j = coarse_bo(1, 2), coarse_bo(2, 2)
2409  coarse_coeffs_pw%array(x_att, j, k) = coarse_coeffs_pw%array(x_att, j, k) + rcv_buf(ii)
2410  ii = ii + 1
2411  END DO
2412  END DO
2413  sent_size(ip) = ii
2414  END IF
2415  END DO
2416 
2417  ii = sent_size(ip)
2418  DO x_att = max(pp_lb(ip), my_coarse_bo(1, 1)), min(pp_ub(ip), my_coarse_bo(2, 1))
2419  DO k = coarse_bo(1, 3), coarse_bo(2, 3)
2420  DO j = coarse_bo(1, 2), coarse_bo(2, 2)
2421  coarse_coeffs_pw%array(x_att, j, k) = coarse_coeffs_pw%array(x_att, j, k) + rcv_buf(ii)
2422  ii = ii + 1
2423  END DO
2424  END DO
2425  END DO
2426  sent_size(ip) = ii
2427 
2428  DO x = coarse_gbo(2, 1) + 1, pp_ub(ip)
2429  x_att = coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1))
2430  IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
2431  ii = sent_size(ip)
2432  DO k = coarse_bo(1, 3), coarse_bo(2, 3)
2433  DO j = coarse_bo(1, 2), coarse_bo(2, 2)
2434  coarse_coeffs_pw%array(x_att, j, k) = coarse_coeffs_pw%array(x_att, j, k) + rcv_buf(ii)
2435  ii = ii + 1
2436  END DO
2437  END DO
2438  sent_size(ip) = ii
2439  END IF
2440  END DO
2441 
2442  END DO
2443 
2444  IF (any(sent_size(0:n_procs - 2) /= rcv_offset(1:n_procs - 1))) &
2445  cpabort("error 1 handling the rcv buffer")
2446  IF (sent_size(n_procs - 1) /= rcv_tot_size) &
2447  cpabort("error 2 handling the rcv buffer")
2448 
2449  ! dealloc
2450  DEALLOCATE (send_size, send_offset, rcv_size, rcv_offset)
2451  DEALLOCATE (send_buf, rcv_buf, real_rcv_size)
2452  DEALLOCATE (pp_ub, pp_lb)
2453  CALL timestop(handle2)
2454  ELSE
2455  cpassert(.NOT. local_data)
2456  END IF
2457 
2458  CALL timestop(handle)
2459  END SUBROUTINE add_fine2coarse
2460 
2461 ! **************************************************************************************************
2462 !> \brief ...
2463 !> \param preconditioner the preconditioner to create
2464 !> \param precond_kind the kind of preconditioner to use
2465 !> \param pool a pool with grids of the same type as the elements to
2466 !> precondition
2467 !> \param pbc if periodic boundary conditions should be applied
2468 !> \param transpose ...
2469 !> \author fawzi
2470 ! **************************************************************************************************
2471  SUBROUTINE pw_spline_precond_create(preconditioner, precond_kind, &
2472  pool, pbc, transpose)
2473  TYPE(pw_spline_precond_type), INTENT(OUT) :: preconditioner
2474  INTEGER, INTENT(in) :: precond_kind
2475  TYPE(pw_pool_type), INTENT(IN), POINTER :: pool
2476  LOGICAL, INTENT(in) :: pbc, transpose
2477 
2478  preconditioner%kind = no_precond
2479  preconditioner%pool => pool
2480  preconditioner%pbc = pbc
2481  preconditioner%transpose = transpose
2482  CALL pool%retain()
2483  CALL pw_spline_precond_set_kind(preconditioner, precond_kind)
2484  END SUBROUTINE pw_spline_precond_create
2485 
2486 ! **************************************************************************************************
2487 !> \brief switches the types of precoditioner to use
2488 !> \param preconditioner the preconditioner to be changed
2489 !> \param precond_kind the new kind of preconditioner to use
2490 !> \param pbc ...
2491 !> \param transpose ...
2492 !> \author fawzi
2493 ! **************************************************************************************************
2494  SUBROUTINE pw_spline_precond_set_kind(preconditioner, precond_kind, pbc, &
2495  transpose)
2496  TYPE(pw_spline_precond_type), INTENT(INOUT) :: preconditioner
2497  INTEGER, INTENT(in) :: precond_kind
2498  LOGICAL, INTENT(in), OPTIONAL :: pbc, transpose
2499 
2500  LOGICAL :: do_3d_coeff
2501  REAL(kind=dp) :: s
2502 
2503  IF (PRESENT(transpose)) preconditioner%transpose = transpose
2504  do_3d_coeff = .false.
2505  preconditioner%kind = precond_kind
2506  IF (PRESENT(pbc)) preconditioner%pbc = pbc
2507  SELECT CASE (precond_kind)
2508  CASE (no_precond)
2509  CASE (precond_spl3_aint2)
2510  preconditioner%coeffs_1d = (/-1.66_dp*0.25_dp, 1.66_dp, -1.66_dp*0.25_dp/)
2511  preconditioner%sharpen = .false.
2512  preconditioner%normalize = .false.
2513  do_3d_coeff = .true.
2514  CASE (precond_spl3_3)
2515  preconditioner%coeffs_1d(1) = -0.25_dp*1.6_dp
2516  preconditioner%coeffs_1d(2) = 1.6_dp
2517  preconditioner%coeffs_1d(3) = -0.25_dp*1.6_dp
2518  preconditioner%sharpen = .false.
2519  preconditioner%normalize = .false.
2520  do_3d_coeff = .true.
2521  CASE (precond_spl3_2)
2522  preconditioner%coeffs_1d(1) = -0.26_dp*1.76_dp
2523  preconditioner%coeffs_1d(2) = 1.76_dp
2524  preconditioner%coeffs_1d(3) = -0.26_dp*1.76_dp
2525  preconditioner%sharpen = .false.
2526  preconditioner%normalize = .false.
2527  do_3d_coeff = .true.
2528  CASE (precond_spl3_aint)
2529  preconditioner%coeffs_1d = spl3_1d_coeffs0
2530  preconditioner%sharpen = .true.
2531  preconditioner%normalize = .true.
2532  do_3d_coeff = .true.
2533  CASE (precond_spl3_1)
2534  preconditioner%coeffs_1d(1) = 0.5_dp/3._dp**(1._dp/3._dp)
2535  preconditioner%coeffs_1d(2) = 4._dp/3._dp**(1._dp/3._dp)
2536  preconditioner%coeffs_1d(3) = 0.5_dp/3._dp**(1._dp/3._dp)
2537  preconditioner%sharpen = .true.
2538  preconditioner%normalize = .false.
2539  do_3d_coeff = .true.
2540  CASE default
2541  cpabort("")
2542  END SELECT
2543  IF (do_3d_coeff) THEN
2544  s = 1._dp
2545  IF (preconditioner%sharpen) s = -1._dp
2546  preconditioner%coeffs(1) = &
2547  s*preconditioner%coeffs_1d(2)* &
2548  preconditioner%coeffs_1d(2)* &
2549  preconditioner%coeffs_1d(2)
2550  preconditioner%coeffs(2) = &
2551  s*preconditioner%coeffs_1d(1)* &
2552  preconditioner%coeffs_1d(2)* &
2553  preconditioner%coeffs_1d(2)
2554  preconditioner%coeffs(3) = &
2555  s*preconditioner%coeffs_1d(1)* &
2556  preconditioner%coeffs_1d(1)* &
2557  preconditioner%coeffs_1d(2)
2558  preconditioner%coeffs(4) = &
2559  s*preconditioner%coeffs_1d(1)* &
2560  preconditioner%coeffs_1d(1)* &
2561  preconditioner%coeffs_1d(1)
2562  IF (preconditioner%sharpen) THEN
2563  IF (preconditioner%normalize) THEN
2564  preconditioner%coeffs(1) = 2._dp + &
2565  preconditioner%coeffs(1)
2566  ELSE
2567  preconditioner%coeffs(1) = -preconditioner%coeffs(1)
2568  END IF
2569  END IF
2570  END IF
2571  END SUBROUTINE pw_spline_precond_set_kind
2572 
2573 ! **************************************************************************************************
2574 !> \brief releases the preconditioner
2575 !> \param preconditioner the preconditioner to release
2576 !> \author fawzi
2577 ! **************************************************************************************************
2578  SUBROUTINE pw_spline_precond_release(preconditioner)
2579  TYPE(pw_spline_precond_type), INTENT(INOUT) :: preconditioner
2580 
2581  CALL pw_pool_release(preconditioner%pool)
2582  END SUBROUTINE pw_spline_precond_release
2583 
2584 ! **************************************************************************************************
2585 !> \brief applies the preconditioner to the system of equations to find the
2586 !> coefficients of the spline
2587 !> \param preconditioner the preconditioner to apply
2588 !> \param in_v the grid on which the preconditioner should be applied
2589 !> \param out_v place to store the preconditioner applied on v_out
2590 !> \author fawzi
2591 ! **************************************************************************************************
2592  SUBROUTINE pw_spline_do_precond(preconditioner, in_v, out_v)
2593  TYPE(pw_spline_precond_type), INTENT(IN) :: preconditioner
2594  TYPE(pw_r3d_rs_type), INTENT(IN) :: in_v
2595  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: out_v
2596 
2597  SELECT CASE (preconditioner%kind)
2598  CASE (no_precond)
2599  CALL pw_copy(in_v, out_v)
2600  CASE (precond_spl3_aint, precond_spl3_1)
2601  CALL pw_zero(out_v)
2602  IF (preconditioner%pbc) THEN
2603  CALL pw_nn_smear_r(pw_in=in_v, pw_out=out_v, &
2604  coeffs=preconditioner%coeffs)
2605  ELSE
2606  CALL pw_nn_compose_r_no_pbc(weights_1d=preconditioner%coeffs_1d, &
2607  pw_in=in_v, pw_out=out_v, sharpen=preconditioner%sharpen, &
2608  normalize=preconditioner%normalize, &
2609  transpose=preconditioner%transpose)
2610  END IF
2611  CASE (precond_spl3_3, precond_spl3_2, precond_spl3_aint2)
2612  CALL pw_zero(out_v)
2613  IF (preconditioner%pbc) THEN
2614  CALL pw_nn_smear_r(pw_in=in_v, pw_out=out_v, &
2615  coeffs=preconditioner%coeffs)
2616  ELSE
2617  CALL pw_nn_compose_r_no_pbc(weights_1d=preconditioner%coeffs_1d, &
2618  pw_in=in_v, pw_out=out_v, sharpen=preconditioner%sharpen, &
2619  normalize=preconditioner%normalize, smooth_boundary=.true., &
2620  transpose=preconditioner%transpose)
2621  END IF
2622  CASE default
2623  cpabort("")
2624  END SELECT
2625  END SUBROUTINE pw_spline_do_precond
2626 
2627 ! **************************************************************************************************
2628 !> \brief solves iteratively (CG) a systmes of linear equations
2629 !> linOp(coeffs)=values
2630 !> (for example those needed to find the coefficients of a spline)
2631 !> Returns true if the it succeeded to achieve the requested accuracy
2632 !> \param values the right hand side of the system
2633 !> \param coeffs will contain the solution of the system (and on entry
2634 !> it contains the starting point)
2635 !> \param linOp the linear operator to be inverted
2636 !> \param preconditioner the preconditioner to apply
2637 !> \param pool a pool of grids (for the temporary objects)
2638 !> \param eps_r the requested precision on the residual
2639 !> \param eps_x the requested precision on the solution
2640 !> \param max_iter maximum number of iteration allowed
2641 !> \param sumtype ...
2642 !> \return ...
2643 !> \author fawzi
2644 ! **************************************************************************************************
2645  FUNCTION find_coeffs(values, coeffs, linOp, preconditioner, pool, &
2646  eps_r, eps_x, max_iter, sumtype) RESULT(res)
2647 
2648  TYPE(pw_r3d_rs_type), INTENT(IN) :: values
2649  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: coeffs
2650  INTERFACE
2651  SUBROUTINE linop(pw_in, pw_out)
2652  USE pw_types, ONLY: pw_r3d_rs_type
2653  TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
2654  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
2655  END SUBROUTINE linop
2656  END INTERFACE
2657  TYPE(pw_spline_precond_type), INTENT(IN) :: preconditioner
2658  TYPE(pw_pool_type), POINTER :: pool
2659  REAL(kind=dp), INTENT(in) :: eps_r, eps_x
2660  INTEGER, INTENT(in) :: max_iter
2661  INTEGER, INTENT(in), OPTIONAL :: sumtype
2662  LOGICAL :: res
2663 
2664  INTEGER :: i, iiter, iter, j, k
2665  INTEGER, DIMENSION(2, 3) :: bo
2666  LOGICAL :: last
2667  REAL(kind=dp) :: alpha, beta, eps_r_att, eps_x_att, r_z, &
2668  r_z_new
2669  TYPE(cp_logger_type), POINTER :: logger
2670  TYPE(pw_r3d_rs_type) :: ap, p, r, z
2671 
2672  last = .false.
2673 
2674  res = .false.
2675  logger => cp_get_default_logger()
2676  CALL pool%create_pw(r)
2677  CALL pool%create_pw(z)
2678  CALL pool%create_pw(p)
2679  CALL pool%create_pw(ap)
2680 
2681  !CALL cp_add_iter_level(logger%iter_info,level_name="SPLINE_FIND_COEFFS")
2682  ext_do: DO iiter = 1, max_iter, 10
2683  CALL pw_zero(r)
2684  CALL linop(pw_in=coeffs, pw_out=r)
2685  r%array = -r%array
2686  CALL pw_axpy(values, r)
2687  CALL pw_spline_do_precond(preconditioner, in_v=r, out_v=z)
2688  CALL pw_copy(z, p)
2689  r_z = pw_integral_ab(r, z, sumtype)
2690 
2691  DO iter = iiter, min(iiter + 9, max_iter)
2692  eps_r_att = sqrt(pw_integral_ab(r, r, sumtype))
2693  IF (eps_r_att == 0._dp) THEN
2694  eps_x_att = 0._dp
2695  last = .true.
2696  ELSE
2697  CALL pw_zero(ap)
2698  CALL linop(pw_in=p, pw_out=ap)
2699  alpha = r_z/pw_integral_ab(ap, p, sumtype)
2700 
2701  CALL pw_axpy(p, coeffs, alpha=alpha)
2702 
2703  eps_x_att = alpha*sqrt(pw_integral_ab(p, p, sumtype)) ! try to spare if unneeded?
2704  IF (eps_r_att < eps_r .AND. eps_x_att < eps_x) last = .true.
2705  END IF
2706  !CALL cp_iterate(logger%iter_info,last=last)
2707  IF (last) THEN
2708  res = .true.
2709  EXIT ext_do
2710  END IF
2711 
2712  CALL pw_axpy(ap, r, alpha=-alpha)
2713 
2714  CALL pw_spline_do_precond(preconditioner, in_v=r, out_v=z)
2715 
2716  r_z_new = pw_integral_ab(r, z, sumtype)
2717  beta = r_z_new/r_z
2718  r_z = r_z_new
2719 
2720  bo = p%pw_grid%bounds_local
2721  DO k = bo(1, 3), bo(2, 3)
2722  DO j = bo(1, 2), bo(2, 2)
2723  DO i = bo(1, 1), bo(2, 1)
2724  p%array(i, j, k) = z%array(i, j, k) + beta*p%array(i, j, k)
2725  END DO
2726  END DO
2727  END DO
2728 
2729  END DO
2730  END DO ext_do
2731  !CALL cp_rm_iter_level(logger%iter_info,level_name="SPLINE_FIND_COEFFS")
2732 
2733  CALL pool%give_back_pw(r)
2734  CALL pool%give_back_pw(z)
2735  CALL pool%give_back_pw(p)
2736  CALL pool%give_back_pw(ap)
2737 
2738  END FUNCTION find_coeffs
2739 
2740 ! **************************************************************************************************
2741 !> \brief adds to pw_out pw_in composed with the weights
2742 !> pw_out%array(i,j,k)=pw_out%array(i,j,k)+sum(pw_in%array(i+l,j+m,k+n)*
2743 !> weights_1d(abs(l)+1)*weights_1d(abs(m)+1)*weights_1d(abs(n)+1),
2744 !> l=-1..1,m=-1..1,n=-1..1)
2745 !> \param weights_1d ...
2746 !> \param pw_in ...
2747 !> \param pw_out ...
2748 !> \param sharpen ...
2749 !> \param normalize ...
2750 !> \param transpose ...
2751 !> \param smooth_boundary ...
2752 !> \author fawzi
2753 ! **************************************************************************************************
2754  SUBROUTINE pw_nn_compose_r_no_pbc(weights_1d, pw_in, pw_out, &
2755  sharpen, normalize, transpose, smooth_boundary)
2756  REAL(kind=dp), DIMENSION(-1:1) :: weights_1d
2757  TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in, pw_out
2758  LOGICAL, INTENT(in), OPTIONAL :: sharpen, normalize, transpose, &
2759  smooth_boundary
2760 
2761  INTEGER :: first_index, i, j, jw, k, kw, &
2762  last_index, myj, myk, n_els
2763  INTEGER, DIMENSION(2, 3) :: bo, gbo
2764  INTEGER, DIMENSION(3) :: s
2765  LOGICAL :: has_l_boundary, has_u_boundary, &
2766  is_split, my_normalize, my_sharpen, &
2767  my_smooth_boundary, my_transpose
2768  REAL(kind=dp) :: in_val_f, in_val_l, in_val_tmp, w_j, w_k
2769  REAL(kind=dp), DIMENSION(-1:1) :: w
2770  REAL(kind=dp), DIMENSION(:, :), POINTER :: l_boundary, tmp, u_boundary
2771  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: in_val, out_val
2772 
2773  bo = pw_in%pw_grid%bounds_local
2774  gbo = pw_in%pw_grid%bounds
2775  in_val => pw_in%array
2776  out_val => pw_out%array
2777  my_sharpen = .false.
2778  IF (PRESENT(sharpen)) my_sharpen = sharpen
2779  my_normalize = .false.
2780  IF (PRESENT(normalize)) my_normalize = normalize
2781  my_transpose = .false.
2782  IF (PRESENT(transpose)) my_transpose = transpose
2783  my_smooth_boundary = .false.
2784  IF (PRESENT(smooth_boundary)) my_smooth_boundary = smooth_boundary
2785  cpassert(.NOT. my_normalize .OR. my_sharpen)
2786  cpassert(.NOT. my_smooth_boundary .OR. .NOT. my_sharpen)
2787  DO i = 1, 3
2788  s(i) = bo(2, i) - bo(1, i) + 1
2789  END DO
2790  IF (any(s < 1)) RETURN
2791  is_split = any(pw_in%pw_grid%bounds_local(:, 1) /= &
2792  pw_in%pw_grid%bounds(:, 1))
2793  has_l_boundary = (gbo(1, 1) == bo(1, 1))
2794  has_u_boundary = (gbo(2, 1) == bo(2, 1))
2795  IF (is_split) THEN
2796  ALLOCATE (l_boundary(bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)), &
2797  u_boundary(bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)), &
2798  tmp(bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
2799  tmp(:, :) = pw_in%array(bo(2, 1), :, :)
2800  CALL pw_in%pw_grid%para%group%sendrecv(tmp, pw_in%pw_grid%para%pos_of_x( &
2801  gbo(1, 1) + modulo(bo(2, 1) + 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)), &
2802  l_boundary, pw_in%pw_grid%para%pos_of_x( &
2803  gbo(1, 1) + modulo(bo(1, 1) - 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)))
2804  tmp(:, :) = pw_in%array(bo(1, 1), :, :)
2805  CALL pw_in%pw_grid%para%group%sendrecv(tmp, pw_in%pw_grid%para%pos_of_x( &
2806  gbo(1, 1) + modulo(bo(1, 1) - 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)), &
2807  u_boundary, pw_in%pw_grid%para%pos_of_x( &
2808  gbo(1, 1) + modulo(bo(2, 1) + 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)))
2809  DEALLOCATE (tmp)
2810  END IF
2811 
2812  n_els = s(1)
2813  IF (has_l_boundary) THEN
2814  n_els = n_els - 1
2815  first_index = bo(1, 1) + 1
2816  ELSE
2817  first_index = bo(1, 1)
2818  END IF
2819  IF (has_u_boundary) THEN
2820  n_els = n_els - 1
2821  last_index = bo(2, 1) - 1
2822  ELSE
2823  last_index = bo(2, 1)
2824  END IF
2825 !$OMP PARALLEL DO DEFAULT(NONE) &
2826 !$OMP PRIVATE(k, kw, myk, j, jw, myj, in_val_f, in_val_l, w_k, w_j, in_val_tmp, w) &
2827 !$OMP SHARED(bo, in_val, out_val, s, l_boundary, u_boundary, weights_1d, is_split, &
2828 !$OMP my_transpose, gbo, my_smooth_boundary, has_l_boundary, has_u_boundary, &
2829 !$OMP my_sharpen, last_index, first_index, my_normalize, n_els)
2830  DO k = bo(1, 3), bo(2, 3)
2831  DO kw = -1, 1
2832  myk = k + kw
2833  IF (my_transpose) THEN
2834  IF (k >= gbo(2, 3) - 1 .OR. k <= gbo(1, 3) + 1) THEN
2835  IF (k == gbo(2, 3) .OR. k == gbo(1, 3)) THEN
2836  IF (myk < gbo(2, 3) .AND. myk > gbo(1, 3)) THEN
2837  w_k = weights_1d(kw)
2838  IF (my_smooth_boundary) THEN
2839  w_k = weights_1d(kw)/weights_1d(0)
2840  END IF
2841  ELSE IF (kw == 0) THEN
2842  w_k = 1._dp
2843  ELSE
2844  cycle
2845  END IF
2846  ELSE
2847  IF (myk == gbo(2, 3) .OR. myk == gbo(1, 3)) cycle
2848  w_k = weights_1d(kw)
2849  END IF
2850  ELSE
2851  w_k = weights_1d(kw)
2852  END IF
2853  ELSE
2854  IF (k >= gbo(2, 3) - 1 .OR. k <= gbo(1, 3) + 1) THEN
2855  IF (k == gbo(2, 3) .OR. k == gbo(1, 3)) THEN
2856  IF (kw /= 0) cycle
2857  w_k = 1._dp
2858  ELSE
2859  IF (my_smooth_boundary .AND. ((k == gbo(1, 3) + 1 .AND. myk == gbo(1, 3)) .OR. &
2860  (k == gbo(2, 3) - 1 .AND. myk == gbo(2, 3)))) THEN
2861  w_k = weights_1d(kw)/weights_1d(0)
2862  ELSE
2863  w_k = weights_1d(kw)
2864  END IF
2865  END IF
2866  ELSE
2867  w_k = weights_1d(kw)
2868  END IF
2869  END IF
2870  DO j = bo(1, 2), bo(2, 2)
2871  DO jw = -1, 1
2872  myj = j + jw
2873  IF (j < gbo(2, 2) - 1 .AND. j > gbo(1, 2) + 1) THEN
2874  w_j = w_k*weights_1d(jw)
2875  ELSE
2876  IF (my_transpose) THEN
2877  IF (j == gbo(2, 2) .OR. j == gbo(1, 2)) THEN
2878  IF (myj < gbo(2, 2) .AND. myj > gbo(1, 2)) THEN
2879  w_j = weights_1d(jw)*w_k
2880  IF (my_smooth_boundary) THEN
2881  w_j = weights_1d(jw)/weights_1d(0)*w_k
2882  END IF
2883  ELSE IF (jw == 0) THEN
2884  w_j = w_k
2885  ELSE
2886  cycle
2887  END IF
2888  ELSE
2889  IF (myj == gbo(2, 2) .OR. myj == gbo(1, 2)) cycle
2890  w_j = w_k*weights_1d(jw)
2891  END IF
2892  ELSE
2893  IF (j == gbo(2, 2) .OR. j == gbo(1, 2)) THEN
2894  IF (jw /= 0) cycle
2895  w_j = w_k
2896  ELSE IF (my_smooth_boundary .AND. ((j == gbo(1, 2) + 1 .AND. myj == gbo(1, 2)) .OR. &
2897  (j == gbo(2, 2) - 1 .AND. myj == gbo(2, 2)))) THEN
2898  w_j = w_k*weights_1d(jw)/weights_1d(0)
2899  ELSE
2900  w_j = w_k*weights_1d(jw)
2901  END IF
2902  END IF
2903  END IF
2904 
2905  IF (has_l_boundary) THEN
2906  IF (my_transpose) THEN
2907  IF (s(1) == 1) THEN
2908  cpassert(.NOT. has_u_boundary)
2909  in_val_tmp = u_boundary(myj, myk)
2910  ELSE
2911  in_val_tmp = in_val(bo(1, 1) + 1, myj, myk)
2912  END IF
2913  IF (my_sharpen) THEN
2914  IF (kw == 0 .AND. jw == 0) THEN
2915  IF (my_normalize) THEN
2916  out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2917  (2.0_dp - w_j)*in_val(bo(1, 1), myj, myk) - &
2918  in_val_tmp*weights_1d(1)*w_j
2919  ELSE
2920  out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2921  in_val(bo(1, 1), myj, myk)*w_j - &
2922  in_val_tmp*weights_1d(1)*w_j
2923  END IF
2924  ELSE
2925  out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) - &
2926  in_val(bo(1, 1), myj, myk)*w_j - &
2927  in_val_tmp*weights_1d(1)*w_j
2928  END IF
2929  ELSE IF (my_smooth_boundary) THEN
2930  out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2931  w_j*(in_val(bo(1, 1), myj, myk) + &
2932  in_val_tmp*weights_1d(1)/weights_1d(0))
2933  ELSE
2934  out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2935  w_j*(in_val(bo(1, 1), myj, myk) + &
2936  in_val_tmp*weights_1d(1))
2937  END IF
2938  in_val_f = 0.0_dp
2939  ELSE
2940  in_val_f = in_val(bo(1, 1), myj, myk)
2941  IF (my_sharpen) THEN
2942  IF (kw == 0 .AND. jw == 0) THEN
2943  IF (my_normalize) THEN
2944  out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2945  (2.0_dp - w_j)*in_val_f
2946  ELSE
2947  out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2948  in_val_f*w_j
2949  END IF
2950  ELSE
2951  out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) - &
2952  in_val_f*w_j
2953  END IF
2954  ELSE
2955  out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2956  in_val_f*w_j
2957  END IF
2958  END IF
2959  ELSE
2960  in_val_f = l_boundary(myj, myk)
2961  END IF
2962  IF (has_u_boundary) THEN
2963  IF (my_transpose) THEN
2964  in_val_l = in_val(bo(2, 1), myj, myk)
2965  IF (s(1) == 1) THEN
2966  cpassert(.NOT. has_l_boundary)
2967  in_val_tmp = l_boundary(myj, myk)
2968  ELSE
2969  in_val_tmp = in_val(bo(2, 1) - 1, myj, myk)
2970  END IF
2971  IF (my_sharpen) THEN
2972  IF (kw == 0 .AND. jw == 0) THEN
2973  IF (my_normalize) THEN
2974  out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
2975  in_val_l*(2._dp - w_j) - &
2976  in_val_tmp*weights_1d(1)*w_j
2977  ELSE
2978  out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
2979  in_val_l*w_j - &
2980  in_val_tmp*weights_1d(1)*w_j
2981  END IF
2982  ELSE
2983  out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) - &
2984  w_j*in_val_l - &
2985  in_val_tmp*weights_1d(1)*w_j
2986  END IF
2987  ELSE IF (my_smooth_boundary) THEN
2988  out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
2989  w_j*(in_val_l + in_val_tmp*weights_1d(1)/weights_1d(0))
2990  ELSE
2991  out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
2992  w_j*(in_val_l + in_val_tmp*weights_1d(1))
2993  END IF
2994  in_val_l = 0._dp
2995  ELSE
2996  in_val_l = in_val(bo(2, 1), myj, myk)
2997  IF (my_sharpen) THEN
2998  IF (kw == 0 .AND. jw == 0) THEN
2999  IF (my_normalize) THEN
3000  out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
3001  in_val_l*(2._dp - w_j)
3002  ELSE
3003  out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
3004  in_val_l*w_j
3005  END IF
3006  ELSE
3007  out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) - &
3008  w_j*in_val_l
3009  END IF
3010  ELSE
3011  out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
3012  w_j*in_val_l
3013  END IF
3014  END IF
3015  ELSE
3016  in_val_l = u_boundary(myj, myk)
3017  END IF
3018  IF (last_index >= first_index) THEN
3019  IF (my_transpose) THEN
3020  IF (bo(1, 1) - 1 == gbo(1, 1)) THEN
3021  in_val_f = 0._dp
3022  ELSE IF (bo(2, 1) + 1 == gbo(2, 1)) THEN
3023  in_val_l = 0._dp
3024  END IF
3025  END IF
3026  IF (my_sharpen) THEN
3027  w = -weights_1d*w_j
3028  IF (kw == 0 .AND. jw == 0) THEN
3029  IF (my_normalize) THEN
3030  w(0) = w(0) + 2._dp
3031  ELSE
3032  w(0) = -w(0)
3033  END IF
3034  END IF
3035  ELSE
3036  w = weights_1d*w_j
3037  END IF
3038  IF (my_smooth_boundary .AND. (.NOT. my_transpose)) THEN
3039  IF (gbo(1, 1) + 1 >= bo(1, 1) .AND. &
3040  gbo(1, 1) + 1 <= bo(2, 1) .AND. gbo(2, 1) - gbo(1, 1) > 2) THEN
3041  IF (gbo(1, 1) >= bo(1, 1)) THEN
3042  out_val(gbo(1, 1) + 1, j, k) = out_val(gbo(1, 1) + 1, j, k) + &
3043  in_val(gbo(1, 1), myj, myk)*w_j*weights_1d(-1)* &
3044  (1._dp/weights_1d(0) - 1._dp)
3045  ELSE
3046  out_val(gbo(1, 1) + 1, j, k) = out_val(gbo(1, 1) + 1, j, k) + &
3047  l_boundary(myj, myk)*w_j*weights_1d(-1)* &
3048  (1._dp/weights_1d(0) - 1._dp)
3049  END IF
3050  END IF
3051  END IF
3052  CALL pw_compose_stripe(weights=w, &
3053  in_val=in_val(first_index:last_index, myj, myk), &
3054  in_val_first=in_val_f, in_val_last=in_val_l, &
3055  out_val=out_val(first_index:last_index, j, k), &
3056  n_el=n_els)
3057 !FM call pw_compose_stripe2(weights=w,&
3058 !FM in_val=in_val,&
3059 !FM in_val_first=in_val_f,in_val_last=in_val_l,&
3060 !FM out_val=out_val,&
3061 !FM first_val=first_index,last_val=last_index,&
3062 !FM myj=myj,myk=myk,j=j,k=k)
3063  IF (my_smooth_boundary .AND. (.NOT. my_transpose)) THEN
3064  IF (gbo(2, 1) - 1 >= bo(1, 1) .AND. &
3065  gbo(2, 1) - 1 <= bo(2, 1) .AND. gbo(2, 1) - gbo(1, 1) > 2) THEN
3066  IF (gbo(2, 1) <= bo(2, 1)) THEN
3067  out_val(gbo(2, 1) - 1, j, k) = out_val(gbo(2, 1) - 1, j, k) + &
3068  in_val(gbo(2, 1), myj, myk)*w_j*weights_1d(1)* &
3069  (1._dp/weights_1d(0) - 1._dp)
3070  ELSE
3071  out_val(gbo(2, 1) - 1, j, k) = out_val(gbo(2, 1) - 1, j, k) + &
3072  u_boundary(myj, myk)*w_j*weights_1d(1)* &
3073  (1._dp/weights_1d(0) - 1._dp)
3074  END IF
3075  END IF
3076  END IF
3077 
3078  END IF
3079  END DO
3080  END DO
3081  END DO
3082  END DO
3083 
3084  IF (is_split) THEN
3085  DEALLOCATE (l_boundary, u_boundary)
3086  END IF
3087  END SUBROUTINE pw_nn_compose_r_no_pbc
3088 
3089 ! **************************************************************************************************
3090 !> \brief ...
3091 !> \param pw_in ...
3092 !> \param pw_out ...
3093 ! **************************************************************************************************
3094  SUBROUTINE spl3_nopbc(pw_in, pw_out)
3095 
3096  TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
3097  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
3098 
3099  CALL pw_zero(pw_out)
3100  CALL pw_nn_compose_r_no_pbc(weights_1d=spl3_1d_coeffs0, pw_in=pw_in, &
3101  pw_out=pw_out, sharpen=.false., normalize=.false.)
3102 
3103  END SUBROUTINE spl3_nopbc
3104 
3105 ! **************************************************************************************************
3106 !> \brief ...
3107 !> \param pw_in ...
3108 !> \param pw_out ...
3109 ! **************************************************************************************************
3110  SUBROUTINE spl3_nopbct(pw_in, pw_out)
3111 
3112  TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
3113  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
3114 
3115  CALL pw_zero(pw_out)
3116  CALL pw_nn_compose_r_no_pbc(weights_1d=spl3_1d_coeffs0, pw_in=pw_in, &
3117  pw_out=pw_out, sharpen=.false., normalize=.false., transpose=.true.)
3118 
3119  END SUBROUTINE spl3_nopbct
3120 
3121 ! **************************************************************************************************
3122 !> \brief ...
3123 !> \param pw_in ...
3124 !> \param pw_out ...
3125 ! **************************************************************************************************
3126  SUBROUTINE spl3_pbc(pw_in, pw_out)
3127 
3128  TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
3129  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
3130 
3131  CALL pw_zero(pw_out)
3132  CALL pw_nn_smear_r(pw_in, pw_out, coeffs=spline3_coeffs)
3133 
3134  END SUBROUTINE spl3_pbc
3135 
3136 ! **************************************************************************************************
3137 !> \brief Evaluates the PBC interpolated Spline (pw) function on the generic
3138 !> input vector (vec)
3139 !> \param vec ...
3140 !> \param pw ...
3141 !> \return ...
3142 !> \par History
3143 !> 12.2007 Adapted for use with distributed grids [rdeclerck]
3144 !> \author Teodoro Laino 12/2005 [tlaino]
3145 !> \note
3146 !> Requires the Spline coefficients to be computed with PBC
3147 ! **************************************************************************************************
3148  FUNCTION eval_interp_spl3_pbc(vec, pw) RESULT(val)
3149  REAL(kind=dp), DIMENSION(3), INTENT(in) :: vec
3150  TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
3151  REAL(kind=dp) :: val
3152 
3153  INTEGER :: i, ivec(3), j, k, npts(3)
3154  INTEGER, DIMENSION(2, 3) :: bo, bo_l
3155  INTEGER, DIMENSION(4) :: ii, ij, ik
3156  LOGICAL :: my_mpsum
3157  REAL(kind=dp) :: a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3, dr1, dr2, dr3, e1, e2, e3, &
3158  f1, f2, f3, g1, g2, g3, h1, h2, h3, p1, p2, p3, q1, q2, q3, r1, r2, r3, s1, s2, s3, s4, &
3159  t1, t2, t3, t4, u1, u2, u3, v1, v2, v3, v4, xd1, xd2, xd3
3160  REAL(kind=dp), DIMENSION(4, 4, 4) :: box
3161  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid
3162 
3163  NULLIFY (grid)
3164  my_mpsum = (pw%pw_grid%para%mode /= pw_mode_local)
3165  npts = pw%pw_grid%npts
3166  ivec = floor(vec/pw%pw_grid%dr)
3167  dr1 = pw%pw_grid%dr(1)
3168  dr2 = pw%pw_grid%dr(2)
3169  dr3 = pw%pw_grid%dr(3)
3170 
3171  xd1 = (vec(1)/dr1) - real(ivec(1), kind=dp)
3172  xd2 = (vec(2)/dr2) - real(ivec(2), kind=dp)
3173  xd3 = (vec(3)/dr3) - real(ivec(3), kind=dp)
3174  grid => pw%array(:, :, :)
3175  bo = pw%pw_grid%bounds
3176  bo_l = pw%pw_grid%bounds_local
3177 
3178  ik(1) = modulo(ivec(3) - 1, npts(3)) + bo(1, 3)
3179  ik(2) = modulo(ivec(3), npts(3)) + bo(1, 3)
3180  ik(3) = modulo(ivec(3) + 1, npts(3)) + bo(1, 3)
3181  ik(4) = modulo(ivec(3) + 2, npts(3)) + bo(1, 3)
3182 
3183  ij(1) = modulo(ivec(2) - 1, npts(2)) + bo(1, 2)
3184  ij(2) = modulo(ivec(2), npts(2)) + bo(1, 2)
3185  ij(3) = modulo(ivec(2) + 1, npts(2)) + bo(1, 2)
3186  ij(4) = modulo(ivec(2) + 2, npts(2)) + bo(1, 2)
3187 
3188  ii(1) = modulo(ivec(1) - 1, npts(1)) + bo(1, 1)
3189  ii(2) = modulo(ivec(1), npts(1)) + bo(1, 1)
3190  ii(3) = modulo(ivec(1) + 1, npts(1)) + bo(1, 1)
3191  ii(4) = modulo(ivec(1) + 2, npts(1)) + bo(1, 1)
3192 
3193  DO k = 1, 4
3194  DO j = 1, 4
3195  DO i = 1, 4
3196  IF ( &
3197  ii(i) >= bo_l(1, 1) .AND. &
3198  ii(i) <= bo_l(2, 1) .AND. &
3199  ij(j) >= bo_l(1, 2) .AND. &
3200  ij(j) <= bo_l(2, 2) .AND. &
3201  ik(k) >= bo_l(1, 3) .AND. &
3202  ik(k) <= bo_l(2, 3) &
3203  ) THEN
3204  box(i, j, k) = grid(ii(i) + 1 - bo_l(1, 1), &
3205  ij(j) + 1 - bo_l(1, 2), &
3206  ik(k) + 1 - bo_l(1, 3))
3207  ELSE
3208  box(i, j, k) = 0.0_dp
3209  END IF
3210  END DO
3211  END DO
3212  END DO
3213 
3214  a1 = 3.0_dp + xd1
3215  a2 = a1*a1
3216  a3 = a2*a1
3217  b1 = 2.0_dp + xd1
3218  b2 = b1*b1
3219  b3 = b2*b1
3220  c1 = 1.0_dp + xd1
3221  c2 = c1*c1
3222  c3 = c2*c1
3223  d1 = xd1
3224  d2 = d1*d1
3225  d3 = d2*d1
3226  e1 = 3.0_dp + xd2
3227  e2 = e1*e1
3228  e3 = e2*e1
3229  f1 = 2.0_dp + xd2
3230  f2 = f1*f1
3231  f3 = f2*f1
3232  g1 = 1.0_dp + xd2
3233  g2 = g1*g1
3234  g3 = g2*g1
3235  h1 = xd2
3236  h2 = h1*h1
3237  h3 = h2*h1
3238  p1 = 3.0_dp + xd3
3239  p2 = p1*p1
3240  p3 = p2*p1
3241  q1 = 2.0_dp + xd3
3242  q2 = q1*q1
3243  q3 = q2*q1
3244  r1 = 1.0_dp + xd3
3245  r2 = r1*r1
3246  r3 = r2*r1
3247  u1 = xd3
3248  u2 = u1*u1
3249  u3 = u2*u1
3250 
3251  t1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3)
3252  t2 = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3
3253  t3 = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3
3254  t4 = 1.0_dp/6.0_dp*d3
3255  s1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3)
3256  s2 = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3
3257  s3 = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3
3258  s4 = 1.0_dp/6.0_dp*h3
3259  v1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3)
3260  v2 = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3
3261  v3 = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3
3262  v4 = 1.0_dp/6.0_dp*u3
3263 
3264  val = ((box(1, 1, 1)*t1 + box(2, 1, 1)*t2 + box(3, 1, 1)*t3 + box(4, 1, 1)*t4)*s1 + &
3265  (box(1, 2, 1)*t1 + box(2, 2, 1)*t2 + box(3, 2, 1)*t3 + box(4, 2, 1)*t4)*s2 + &
3266  (box(1, 3, 1)*t1 + box(2, 3, 1)*t2 + box(3, 3, 1)*t3 + box(4, 3, 1)*t4)*s3 + &
3267  (box(1, 4, 1)*t1 + box(2, 4, 1)*t2 + box(3, 4, 1)*t3 + box(4, 4, 1)*t4)*s4)*v1 + &
3268  ((box(1, 1, 2)*t1 + box(2, 1, 2)*t2 + box(3, 1, 2)*t3 + box(4, 1, 2)*t4)*s1 + &
3269  (box(1, 2, 2)*t1 + box(2, 2, 2)*t2 + box(3, 2, 2)*t3 + box(4, 2, 2)*t4)*s2 + &
3270  (box(1, 3, 2)*t1 + box(2, 3, 2)*t2 + box(3, 3, 2)*t3 + box(4, 3, 2)*t4)*s3 + &
3271  (box(1, 4, 2)*t1 + box(2, 4, 2)*t2 + box(3, 4, 2)*t3 + box(4, 4, 2)*t4)*s4)*v2 + &
3272  ((box(1, 1, 3)*t1 + box(2, 1, 3)*t2 + box(3, 1, 3)*t3 + box(4, 1, 3)*t4)*s1 + &
3273  (box(1, 2, 3)*t1 + box(2, 2, 3)*t2 + box(3, 2, 3)*t3 + box(4, 2, 3)*t4)*s2 + &
3274  (box(1, 3, 3)*t1 + box(2, 3, 3)*t2 + box(3, 3, 3)*t3 + box(4, 3, 3)*t4)*s3 + &
3275  (box(1, 4, 3)*t1 + box(2, 4, 3)*t2 + box(3, 4, 3)*t3 + box(4, 4, 3)*t4)*s4)*v3 + &
3276  ((box(1, 1, 4)*t1 + box(2, 1, 4)*t2 + box(3, 1, 4)*t3 + box(4, 1, 4)*t4)*s1 + &
3277  (box(1, 2, 4)*t1 + box(2, 2, 4)*t2 + box(3, 2, 4)*t3 + box(4, 2, 4)*t4)*s2 + &
3278  (box(1, 3, 4)*t1 + box(2, 3, 4)*t2 + box(3, 3, 4)*t3 + box(4, 3, 4)*t4)*s3 + &
3279  (box(1, 4, 4)*t1 + box(2, 4, 4)*t2 + box(3, 4, 4)*t3 + box(4, 4, 4)*t4)*s4)*v4
3280 
3281  IF (my_mpsum) CALL pw%pw_grid%para%group%sum(val)
3282 
3283  END FUNCTION eval_interp_spl3_pbc
3284 
3285 ! **************************************************************************************************
3286 !> \brief Evaluates the derivatives of the PBC interpolated Spline (pw)
3287 !> function on the generic input vector (vec)
3288 !> \param vec ...
3289 !> \param pw ...
3290 !> \return ...
3291 !> \par History
3292 !> 12.2007 Adapted for use with distributed grids [rdeclerck]
3293 !> \author Teodoro Laino 12/2005 [tlaino]
3294 !> \note
3295 !> Requires the Spline coefficients to be computed with PBC
3296 ! **************************************************************************************************
3297  FUNCTION eval_d_interp_spl3_pbc(vec, pw) RESULT(val)
3298  REAL(kind=dp), DIMENSION(3), INTENT(in) :: vec
3299  TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
3300  REAL(kind=dp) :: val(3)
3301 
3302  INTEGER :: i, ivec(3), j, k, npts(3)
3303  INTEGER, DIMENSION(2, 3) :: bo, bo_l
3304  INTEGER, DIMENSION(4) :: ii, ij, ik
3305  LOGICAL :: my_mpsum
3306  REAL(kind=dp) :: a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3, dr1, dr1i, dr2, dr2i, dr3, &
3307  dr3i, e1, e2, e3, f1, f2, f3, g1, g2, g3, h1, h2, h3, p1, p2, p3, q1, q2, q3, r1, r2, r3, &
3308  s1, s1d, s1o, s2, s2d, s2o, s3, s3d, s3o, s4, s4d, s4o, t1, t1d, t1o, t2, t2d, t2o, t3, &
3309  t3d, t3o, t4, t4d, t4o, u1, u2, u3, v1, v1d, v1o, v2, v2d, v2o, v3, v3d, v3o, v4, v4d, &
3310  v4o, xd1, xd2, xd3
3311  REAL(kind=dp), DIMENSION(4, 4, 4) :: box
3312  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid
3313 
3314  NULLIFY (grid)
3315  my_mpsum = (pw%pw_grid%para%mode /= pw_mode_local)
3316  npts = pw%pw_grid%npts
3317  ivec = floor(vec/pw%pw_grid%dr)
3318  dr1 = pw%pw_grid%dr(1)
3319  dr2 = pw%pw_grid%dr(2)
3320  dr3 = pw%pw_grid%dr(3)
3321  dr1i = 1.0_dp/dr1
3322  dr2i = 1.0_dp/dr2
3323  dr3i = 1.0_dp/dr3
3324  xd1 = (vec(1)/dr1) - real(ivec(1), kind=dp)
3325  xd2 = (vec(2)/dr2) - real(ivec(2), kind=dp)
3326  xd3 = (vec(3)/dr3) - real(ivec(3), kind=dp)
3327  grid => pw%array(:, :, :)
3328  bo = pw%pw_grid%bounds
3329  bo_l = pw%pw_grid%bounds_local
3330 
3331  ik(1) = modulo(ivec(3) - 1, npts(3)) + bo(1, 3)
3332  ik(2) = modulo(ivec(3), npts(3)) + bo(1, 3)
3333  ik(3) = modulo(ivec(3) + 1, npts(3)) + bo(1, 3)
3334  ik(4) = modulo(ivec(3) + 2, npts(3)) + bo(1, 3)
3335 
3336  ij(1) = modulo(ivec(2) - 1, npts(2)) + bo(1, 2)
3337  ij(2) = modulo(ivec(2), npts(2)) + bo(1, 2)
3338  ij(3) = modulo(ivec(2) + 1, npts(2)) + bo(1, 2)
3339  ij(4) = modulo(ivec(2) + 2, npts(2)) + bo(1, 2)
3340 
3341  ii(1) = modulo(ivec(1) - 1, npts(1)) + bo(1, 1)
3342  ii(2) = modulo(ivec(1), npts(1)) + bo(1, 1)
3343  ii(3) = modulo(ivec(1) + 1, npts(1)) + bo(1, 1)
3344  ii(4) = modulo(ivec(1) + 2, npts(1)) + bo(1, 1)
3345 
3346  DO k = 1, 4
3347  DO j = 1, 4
3348  DO i = 1, 4
3349  IF ( &
3350  ii(i) >= bo_l(1, 1) .AND. &
3351  ii(i) <= bo_l(2, 1) .AND. &
3352  ij(j) >= bo_l(1, 2) .AND. &
3353  ij(j) <= bo_l(2, 2) .AND. &
3354  ik(k) >= bo_l(1, 3) .AND. &
3355  ik(k) <= bo_l(2, 3) &
3356  ) THEN
3357  box(i, j, k) = grid(ii(i) + 1 - bo_l(1, 1), &
3358  ij(j) + 1 - bo_l(1, 2), &
3359  ik(k) + 1 - bo_l(1, 3))
3360  ELSE
3361  box(i, j, k) = 0.0_dp
3362  END IF
3363  END DO
3364  END DO
3365  END DO
3366 
3367  a1 = 3.0_dp + xd1
3368  a2 = a1*a1
3369  a3 = a2*a1
3370  b1 = 2.0_dp + xd1
3371  b2 = b1*b1
3372  b3 = b2*b1
3373  c1 = 1.0_dp + xd1
3374  c2 = c1*c1
3375  c3 = c2*c1
3376  d1 = xd1
3377  d2 = d1*d1
3378  d3 = d2*d1
3379  e1 = 3.0_dp + xd2
3380  e2 = e1*e1
3381  e3 = e2*e1
3382  f1 = 2.0_dp + xd2
3383  f2 = f1*f1
3384  f3 = f2*f1
3385  g1 = 1.0_dp + xd2
3386  g2 = g1*g1
3387  g3 = g2*g1
3388  h1 = xd2
3389  h2 = h1*h1
3390  h3 = h2*h1
3391  p1 = 3.0_dp + xd3
3392  p2 = p1*p1
3393  p3 = p2*p1
3394  q1 = 2.0_dp + xd3
3395  q2 = q1*q1
3396  q3 = q2*q1
3397  r1 = 1.0_dp + xd3
3398  r2 = r1*r1
3399  r3 = r2*r1
3400  u1 = xd3
3401  u2 = u1*u1
3402  u3 = u2*u1
3403 
3404  t1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3)
3405  t2o = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3
3406  t3o = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3
3407  t4o = 1.0_dp/6.0_dp*d3
3408  s1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3)
3409  s2o = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3
3410  s3o = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3
3411  s4o = 1.0_dp/6.0_dp*h3
3412  v1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3)
3413  v2o = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3
3414  v3o = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3
3415  v4o = 1.0_dp/6.0_dp*u3
3416 
3417  t1d = -8.0_dp + 4.0_dp*a1 - 0.5_dp*a2
3418  t2d = 10.0_dp - 8.0_dp*b1 + 1.5_dp*b2
3419  t3d = -2.0_dp + 4.0_dp*c1 - 1.5_dp*c2
3420  t4d = 0.5_dp*d2
3421  s1d = -8.0_dp + 4.0_dp*e1 - 0.5_dp*e2
3422  s2d = 10.0_dp - 8.0_dp*f1 + 1.5_dp*f2
3423  s3d = -2.0_dp + 4.0_dp*g1 - 1.5_dp*g2
3424  s4d = 0.5_dp*h2
3425  v1d = -8.0_dp + 4.0_dp*p1 - 0.5_dp*p2
3426  v2d = 10.0_dp - 8.0_dp*q1 + 1.5_dp*q2
3427  v3d = -2.0_dp + 4.0_dp*r1 - 1.5_dp*r2
3428  v4d = 0.5_dp*u2
3429 
3430  t1 = t1d*dr1i
3431  t2 = t2d*dr1i
3432  t3 = t3d*dr1i
3433  t4 = t4d*dr1i
3434  s1 = s1o
3435  s2 = s2o
3436  s3 = s3o
3437  s4 = s4o
3438  v1 = v1o
3439  v2 = v2o
3440  v3 = v3o
3441  v4 = v4o
3442  val(1) = ((box(1, 1, 1)*t1 + box(2, 1, 1)*t2 + box(3, 1, 1)*t3 + box(4, 1, 1)*t4)*s1 + &
3443  (box(1, 2, 1)*t1 + box(2, 2, 1)*t2 + box(3, 2, 1)*t3 + box(4, 2, 1)*t4)*s2 + &
3444  (box(1, 3, 1)*t1 + box(2, 3, 1)*t2 + box(3, 3, 1)*t3 + box(4, 3, 1)*t4)*s3 + &
3445  (box(1, 4, 1)*t1 + box(2, 4, 1)*t2 + box(3, 4, 1)*t3 + box(4, 4, 1)*t4)*s4)*v1 + &
3446  ((box(1, 1, 2)*t1 + box(2, 1, 2)*t2 + box(3, 1, 2)*t3 + box(4, 1, 2)*t4)*s1 + &
3447  (box(1, 2, 2)*t1 + box(2, 2, 2)*t2 + box(3, 2, 2)*t3 + box(4, 2, 2)*t4)*s2 + &
3448  (box(1, 3, 2)*t1 + box(2, 3, 2)*t2 + box(3, 3, 2)*t3 + box(4, 3, 2)*t4)*s3 + &
3449  (box(1, 4, 2)*t1 + box(2, 4, 2)*t2 + box(3, 4, 2)*t3 + box(4, 4, 2)*t4)*s4)*v2 + &
3450  ((box(1, 1, 3)*t1 + box(2, 1, 3)*t2 + box(3, 1, 3)*t3 + box(4, 1, 3)*t4)*s1 + &
3451  (box(1, 2, 3)*t1 + box(2, 2, 3)*t2 + box(3, 2, 3)*t3 + box(4, 2, 3)*t4)*s2 + &
3452  (box(1, 3, 3)*t1 + box(2, 3, 3)*t2 + box(3, 3, 3)*t3 + box(4, 3, 3)*t4)*s3 + &
3453  (box(1, 4, 3)*t1 + box(2, 4, 3)*t2 + box(3, 4, 3)*t3 + box(4, 4, 3)*t4)*s4)*v3 + &
3454  ((box(1, 1, 4)*t1 + box(2, 1, 4)*t2 + box(3, 1, 4)*t3 + box(4, 1, 4)*t4)*s1 + &
3455  (box(1, 2, 4)*t1 + box(2, 2, 4)*t2 + box(3, 2, 4)*t3 + box(4, 2, 4)*t4)*s2 + &
3456  (box(1, 3, 4)*t1 + box(2, 3, 4)*t2 + box(3, 3, 4)*t3 + box(4, 3, 4)*t4)*s3 + &
3457  (box(1, 4, 4)*t1 + box(2, 4, 4)*t2 + box(3, 4, 4)*t3 + box(4, 4, 4)*t4)*s4)*v4
3458 
3459  t1 = t1o
3460  t2 = t2o
3461  t3 = t3o
3462  t4 = t4o
3463  s1 = s1d*dr2i
3464  s2 = s2d*dr2i
3465  s3 = s3d*dr2i
3466  s4 = s4d*dr2i
3467  v1 = v1o
3468  v2 = v2o
3469  v3 = v3o
3470  v4 = v4o
3471  val(2) = ((box(1, 1, 1)*t1 + box(2, 1, 1)*t2 + box(3, 1, 1)*t3 + box(4, 1, 1)*t4)*s1 + &
3472  (box(1, 2, 1)*t1 + box(2, 2, 1)*t2 + box(3, 2, 1)*t3 + box(4, 2, 1)*t4)*s2 + &
3473  (box(1, 3, 1)*t1 + box(2, 3, 1)*t2 + box(3, 3, 1)*t3 + box(4, 3, 1)*t4)*s3 + &
3474  (box(1, 4, 1)*t1 + box(2, 4, 1)*t2 + box(3, 4, 1)*t3 + box(4, 4, 1)*t4)*s4)*v1 + &
3475  ((box(1, 1, 2)*t1 + box(2, 1, 2)*t2 + box(3, 1, 2)*t3 + box(4, 1, 2)*t4)*s1 + &
3476  (box(1, 2, 2)*t1 + box(2, 2, 2)*t2 + box(3, 2, 2)*t3 + box(4, 2, 2)*t4)*s2 + &
3477  (box(1, 3, 2)*t1 + box(2, 3, 2)*t2 + box(3, 3, 2)*t3 + box(4, 3, 2)*t4)*s3 + &
3478  (box(1, 4, 2)*t1 + box(2, 4, 2)*t2 + box(3, 4, 2)*t3 + box(4, 4, 2)*t4)*s4)*v2 + &
3479  ((box(1, 1, 3)*t1 + box(2, 1, 3)*t2 + box(3, 1, 3)*t3 + box(4, 1, 3)*t4)*s1 + &
3480  (box(1, 2, 3)*t1 + box(2, 2, 3)*t2 + box(3, 2, 3)*t3 + box(4, 2, 3)*t4)*s2 + &
3481  (box(1, 3, 3)*t1 + box(2, 3, 3)*t2 + box(3, 3, 3)*t3 + box(4, 3, 3)*t4)*s3 + &
3482  (box(1, 4, 3)*t1 + box(2, 4, 3)*t2 + box(3, 4, 3)*t3 + box(4, 4, 3)*t4)*s4)*v3 + &
3483  ((box(1, 1, 4)*t1 + box(2, 1, 4)*t2 + box(3, 1, 4)*t3 + box(4, 1, 4)*t4)*s1 + &
3484  (box(1, 2, 4)*t1 + box(2, 2, 4)*t2 + box(3, 2, 4)*t3 + box(4, 2, 4)*t4)*s2 + &
3485  (box(1, 3, 4)*t1 + box(2, 3, 4)*t2 + box(3, 3, 4)*t3 + box(4, 3, 4)*t4)*s3 + &
3486  (box(1, 4, 4)*t1 + box(2, 4, 4)*t2 + box(3, 4, 4)*t3 + box(4, 4, 4)*t4)*s4)*v4
3487 
3488  t1 = t1o
3489  t2 = t2o
3490  t3 = t3o
3491  t4 = t4o
3492  s1 = s1o
3493  s2 = s2o
3494  s3 = s3o
3495  s4 = s4o
3496  v1 = v1d*dr3i
3497  v2 = v2d*dr3i
3498  v3 = v3d*dr3i
3499  v4 = v4d*dr3i
3500  val(3) = ((box(1, 1, 1)*t1 + box(2, 1, 1)*t2 + box(3, 1, 1)*t3 + box(4, 1, 1)*t4)*s1 + &
3501  (box(1, 2, 1)*t1 + box(2, 2, 1)*t2 + box(3, 2, 1)*t3 + box(4, 2, 1)*t4)*s2 + &
3502  (box(1, 3, 1)*t1 + box(2, 3, 1)*t2 + box(3, 3, 1)*t3 + box(4, 3, 1)*t4)*s3 + &
3503  (box(1, 4, 1)*t1 + box(2, 4, 1)*t2 + box(3, 4, 1)*t3 + box(4, 4, 1)*t4)*s4)*v1 + &
3504  ((box(1, 1, 2)*t1 + box(2, 1, 2)*t2 + box(3, 1, 2)*t3 + box(4, 1, 2)*t4)*s1 + &
3505  (box(1, 2, 2)*t1 + box(2, 2, 2)*t2 + box(3, 2, 2)*t3 + box(4, 2, 2)*t4)*s2 + &
3506  (box(1, 3, 2)*t1 + box(2, 3, 2)*t2 + box(3, 3, 2)*t3 + box(4, 3, 2)*t4)*s3 + &
3507  (box(1, 4, 2)*t1 + box(2, 4, 2)*t2 + box(3, 4, 2)*t3 + box(4, 4, 2)*t4)*s4)*v2 + &
3508  ((box(1, 1, 3)*t1 + box(2, 1, 3)*t2 + box(3, 1, 3)*t3 + box(4, 1, 3)*t4)*s1 + &
3509  (box(1, 2, 3)*t1 + box(2, 2, 3)*t2 + box(3, 2, 3)*t3 + box(4, 2, 3)*t4)*s2 + &
3510  (box(1, 3, 3)*t1 + box(2, 3, 3)*t2 + box(3, 3, 3)*t3 + box(4, 3, 3)*t4)*s3 + &
3511  (box(1, 4, 3)*t1 + box(2, 4, 3)*t2 + box(3, 4, 3)*t3 + box(4, 4, 3)*t4)*s4)*v3 + &
3512  ((box(1, 1, 4)*t1 + box(2, 1, 4)*t2 + box(3, 1, 4)*t3 + box(4, 1, 4)*t4)*s1 + &
3513  (box(1, 2, 4)*t1 + box(2, 2, 4)*t2 + box(3, 2, 4)*t3 + box(4, 2, 4)*t4)*s2 + &
3514  (box(1, 3, 4)*t1 + box(2, 3, 4)*t2 + box(3, 3, 4)*t3 + box(4, 3, 4)*t4)*s3 + &
3515  (box(1, 4, 4)*t1 + box(2, 4, 4)*t2 + box(3, 4, 4)*t3 + box(4, 4, 4)*t4)*s4)*v4
3516 
3517  IF (my_mpsum) CALL pw%pw_grid%para%group%sum(val)
3518 
3519  END FUNCTION eval_d_interp_spl3_pbc
3520 
3521 END MODULE pw_spline_utils
3522 
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
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:117
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition: mathconstants.F:112
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
message_passing::mp_comm_congruent
integer, parameter, public mp_comm_congruent
Definition: message_passing.F:827
preconditioner
computes preconditioners, and implements methods to apply them currently used in qs_ot
Definition: preconditioner.F:15
pw_grid_types
Definition: pw_grid_types.F:13
pw_grid_types::pw_mode_local
integer, parameter, public pw_mode_local
Definition: pw_grid_types.F:29
pw_grid_types::fullspace
integer, parameter, public fullspace
Definition: pw_grid_types.F:28
pw_methods
Definition: pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
pw_pool_types::pw_pool_release
subroutine, public pw_pool_release(pool)
releases the given pool (see cp2k/doc/ReferenceCounting.html)
Definition: pw_pool_types.F:259
pw_spline_utils
different utils that are useful to manipulate splines on the regular grid of a pw
Definition: pw_spline_utils.F:18
pw_spline_utils::find_coeffs
logical function, public find_coeffs(values, coeffs, linOp, preconditioner, pool, eps_r, eps_x, max_iter, sumtype)
solves iteratively (CG) a systmes of linear equations linOp(coeffs)=values (for example those needed ...
Definition: pw_spline_utils.F:2649
pw_spline_utils::add_fine2coarse
subroutine, public add_fine2coarse(fine_values_pw, coarse_coeffs_pw, weights_1d, w_border0, w_border1, pbc, safe_computation)
low level function that adds a coarse grid (without boundary) to a fine grid.
Definition: pw_spline_utils.F:1772
pw_spline_utils::precond_spl3_3
integer, parameter, public precond_spl3_3
Definition: pw_spline_utils.F:44
pw_spline_utils::pw_spline_precond_release
subroutine, public pw_spline_precond_release(preconditioner)
releases the preconditioner
Definition: pw_spline_utils.F:2581
pw_spline_utils::pw_spline_precond_create
subroutine, public pw_spline_precond_create(preconditioner, precond_kind, pool, pbc, transpose)
...
Definition: pw_spline_utils.F:2475
pw_spline_utils::pw_nn_deriv_r
subroutine, public pw_nn_deriv_r(pw_in, pw_out, coeffs, idir)
calculates a nearest neighbor central derivative. for the x dir: pw_outarray(i,j,k)=( pw_in(i+1,...
Definition: pw_spline_utils.F:899
pw_spline_utils::pw_spline_do_precond
subroutine, public pw_spline_do_precond(preconditioner, in_v, out_v)
applies the preconditioner to the system of equations to find the coefficients of the spline
Definition: pw_spline_utils.F:2595
pw_spline_utils::pw_spline3_deriv_g
subroutine, public pw_spline3_deriv_g(spline_g, idir)
calculates the FFT of the values of the x,y,z (idir=1,2,3) derivative of the cubic spline
Definition: pw_spline_utils.F:366
pw_spline_utils::pw_spline_precond_set_kind
subroutine, public pw_spline_precond_set_kind(preconditioner, precond_kind, pbc, transpose)
switches the types of precoditioner to use
Definition: pw_spline_utils.F:2498
pw_spline_utils::spl3_1d_transf_coeffs
real(kind=dp), dimension(4), parameter, public spl3_1d_transf_coeffs
Definition: pw_spline_utils.F:51
pw_spline_utils::spl3_1d_transf_border1
real(kind=dp), dimension(3), parameter, public spl3_1d_transf_border1
Definition: pw_spline_utils.F:70
pw_spline_utils::spline2_coeffs
real(kind=dp), dimension(4), parameter, public spline2_coeffs
Definition: pw_spline_utils.F:51
pw_spline_utils::add_coarse2fine
subroutine, public add_coarse2fine(coarse_coeffs_pw, fine_values_pw, weights_1d, w_border0, w_border1, pbc, safe_computation)
low level function that adds a coarse grid to a fine grid. If pbc is true periodic boundary condition...
Definition: pw_spline_utils.F:977
pw_spline_utils::eval_d_interp_spl3_pbc
real(kind=dp) function, dimension(3), public eval_d_interp_spl3_pbc(vec, pw)
Evaluates the derivatives of the PBC interpolated Spline (pw) function on the generic input vector (v...
Definition: pw_spline_utils.F:3300
pw_spline_utils::eval_interp_spl3_pbc
real(kind=dp) function, public eval_interp_spl3_pbc(vec, pw)
Evaluates the PBC interpolated Spline (pw) function on the generic input vector (vec)
Definition: pw_spline_utils.F:3151
pw_spline_utils::spline3_deriv_coeffs
real(kind=dp), dimension(3), parameter, public spline3_deriv_coeffs
Definition: pw_spline_utils.F:70
pw_spline_utils::precond_spl3_aint
integer, parameter, public precond_spl3_aint
Definition: pw_spline_utils.F:44
pw_spline_utils::pw_spline_scale_deriv
subroutine, public pw_spline_scale_deriv(deriv_vals_r, transpose, scale)
rescales the derivatives from gridspacing=1 to the real derivatives
Definition: pw_spline_utils.F:275
pw_spline_utils::spline2_deriv_coeffs
real(kind=dp), dimension(3), parameter, public spline2_deriv_coeffs
Definition: pw_spline_utils.F:70
pw_spline_utils::pw_nn_smear_r
subroutine, public pw_nn_smear_r(pw_in, pw_out, coeffs)
calculates the values of a nearest neighbor smearing
Definition: pw_spline_utils.F:856
pw_spline_utils::spline3_coeffs
real(kind=dp), dimension(4), parameter, public spline3_coeffs
Definition: pw_spline_utils.F:51
pw_spline_utils::spl3_nopbc
subroutine, public spl3_nopbc(pw_in, pw_out)
...
Definition: pw_spline_utils.F:3097
pw_spline_utils::nn50_deriv_coeffs
real(kind=dp), dimension(3), parameter, public nn50_deriv_coeffs
Definition: pw_spline_utils.F:70
pw_spline_utils::pw_spline3_interpolate_values_g
subroutine, public pw_spline3_interpolate_values_g(spline_g)
calculates the FFT of the coefficients of the2 cubic spline that interpolates the given values
Definition: pw_spline_utils.F:200
pw_spline_utils::pw_spline2_interpolate_values_g
subroutine, public pw_spline2_interpolate_values_g(spline_g)
calculates the FFT of the coefficients of the quadratic spline that interpolates the given values
Definition: pw_spline_utils.F:129
pw_spline_utils::spl3_nopbct
subroutine, public spl3_nopbct(pw_in, pw_out)
...
Definition: pw_spline_utils.F:3113
pw_spline_utils::nn10_deriv_coeffs
real(kind=dp), dimension(3), parameter, public nn10_deriv_coeffs
Definition: pw_spline_utils.F:70
pw_spline_utils::spl3_precond1_coeff
real(kind=dp), dimension(4), parameter, public spl3_precond1_coeff
Definition: pw_spline_utils.F:51
pw_spline_utils::spl3_1d_coeffs0
real(kind=dp), dimension(3), parameter, public spl3_1d_coeffs0
Definition: pw_spline_utils.F:70
pw_spline_utils::pw_spline2_deriv_g
subroutine, public pw_spline2_deriv_g(spline_g, idir)
calculates the FFT of the values of the x,y,z (idir=1,2,3) derivative of the quadratic spline
Definition: pw_spline_utils.F:501
pw_spline_utils::spl3_aint_coeff
real(kind=dp), dimension(4), parameter, public spl3_aint_coeff
Definition: pw_spline_utils.F:51
pw_spline_utils::spl3_pbc
subroutine, public spl3_pbc(pw_in, pw_out)
...
Definition: pw_spline_utils.F:3129
pw_spline_utils::no_precond
integer, parameter, public no_precond
Definition: pw_spline_utils.F:44
pw_spline_utils::pw_nn_compose_r_no_pbc
subroutine pw_nn_compose_r_no_pbc(weights_1d, pw_in, pw_out, sharpen, normalize, transpose, smooth_boundary)
adds to pw_out pw_in composed with the weights pw_outarray(i,j,k)=pw_outarray(i,j,...
Definition: pw_spline_utils.F:2758
pw_spline_utils::precond_spl3_2
integer, parameter, public precond_spl3_2
Definition: pw_spline_utils.F:44
pw_spline_utils::nn10_coeffs
real(kind=dp), dimension(4), parameter, public nn10_coeffs
Definition: pw_spline_utils.F:51
pw_spline_utils::precond_spl3_aint2
integer, parameter, public precond_spl3_aint2
Definition: pw_spline_utils.F:44
pw_spline_utils::nn50_coeffs
real(kind=dp), dimension(4), parameter, public nn50_coeffs
Definition: pw_spline_utils.F:51
pw_spline_utils::precond_spl3_1
integer, parameter, public precond_spl3_1
Definition: pw_spline_utils.F:44
pw_types
Definition: pw_types.F:24