(git:fdbe441)
Loading...
Searching...
No Matches
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-2026 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! **************************************************************************************************
18MODULE 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
115CONTAINS
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 cpassert(ii <= mp_comm_congruent)
1006 my_coarse_bo = coarse_coeffs_pw%pw_grid%bounds_local
1007 coarse_gbo = coarse_coeffs_pw%pw_grid%bounds
1008 fine_bo = fine_values_pw%pw_grid%bounds_local
1009 fine_gbo = fine_values_pw%pw_grid%bounds
1010 f_shift = fine_gbo(1, :) - 2*coarse_gbo(1, :)
1011 DO j = 2, 3
1012 DO i = 1, 2
1013 coarse_bo(i, j) = floor((fine_bo(i, j) - f_shift(j))/2.)
1014 END DO
1015 END DO
1016 IF (fine_bo(1, 1) <= fine_bo(2, 1)) THEN
1017 coarse_bo(1, 1) = floor((fine_bo(1, 1) - 2 - f_shift(1))/2.)
1018 coarse_bo(2, 1) = floor((fine_bo(2, 1) + 3 - f_shift(1))/2.)
1019 ELSE
1020 coarse_bo(1, 1) = coarse_gbo(2, 1)
1021 coarse_bo(2, 1) = coarse_gbo(2, 1) - 1
1022 END IF
1023 is_split = any(coarse_gbo(:, 1) /= my_coarse_bo(:, 1))
1024 IF (.NOT. is_split .OR. .NOT. pbc) THEN
1025 coarse_bo(1, 1) = max(coarse_gbo(1, 1), coarse_bo(1, 1))
1026 coarse_bo(2, 1) = min(coarse_gbo(2, 1), coarse_bo(2, 1))
1027 END IF
1028 has_i_ubound = (fine_gbo(2, 1) /= fine_bo(2, 1)) .OR. pbc .AND. is_split
1029 has_i_lbound = (fine_gbo(1, 1) /= fine_bo(1, 1)) .OR. pbc .AND. is_split
1030
1031 IF (pbc) THEN
1032 cpassert(all(fine_gbo(1, :) == 2*coarse_gbo(1, :) + f_shift))
1033 cpassert(all(fine_gbo(2, :) == 2*coarse_gbo(2, :) + 1 + f_shift))
1034 ELSE
1035 cpassert(all(fine_gbo(2, :) == 2*coarse_gbo(2, :) + f_shift))
1036 cpassert(all(fine_gbo(1, :) == 2*coarse_gbo(1, :) + f_shift))
1037 END IF
1038
1039 coarse_coeffs => coarse_coeffs_pw%array
1040 DO i = 1, 3
1041 s(i) = coarse_gbo(2, i) - coarse_gbo(1, i) + 1
1042 END DO
1043! CALL timestop(handle2)
1044 ! *** parallel case
1045 IF (is_split) THEN
1046 CALL timeset(routinen//"_comm", handle2)
1047 coarse_slice_size = (coarse_bo(2, 2) - coarse_bo(1, 2) + 1)* &
1048 (coarse_bo(2, 3) - coarse_bo(1, 3) + 1)
1049 n_procs = coarse_coeffs_pw%pw_grid%para%group%num_pe
1050 ALLOCATE (send_size(0:n_procs - 1), send_offset(0:n_procs - 1), &
1051 sent_size(0:n_procs - 1), rcv_size(0:n_procs - 1), &
1052 rcv_offset(0:n_procs - 1), real_rcv_size(0:n_procs - 1))
1053
1054 ! ** rcv size count
1055
1056 pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x
1057 p_old = pos_of_x(coarse_gbo(1, 1) &
1058 + modulo(coarse_bo(1, 1) - coarse_gbo(1, 1), s(1)))
1059 rcv_size = 0
1060 DO x = coarse_bo(1, 1), coarse_bo(2, 1)
1061 p = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
1062 rcv_size(p) = rcv_size(p) + coarse_slice_size
1063 END DO
1064
1065 ! ** send size count
1066
1067 pos_of_x => fine_values_pw%pw_grid%para%pos_of_x
1068 sf = fine_gbo(2, 1) - fine_gbo(1, 1) + 1
1069 fi_lb = 2*my_coarse_bo(1, 1) - 3 + f_shift(1)
1070 fi_ub = 2*my_coarse_bo(2, 1) + 3 + f_shift(1)
1071 IF (.NOT. pbc) THEN
1072 fi_lb = max(fi_lb, fine_gbo(1, 1))
1073 fi_ub = min(fi_ub, fine_gbo(2, 1))
1074 ELSE
1075 fi_ub = min(fi_ub, fi_lb + sf - 1)
1076 END IF
1077 p_old = pos_of_x(fine_gbo(1, 1) + modulo(fi_lb - fine_gbo(1, 1), sf))
1078 p_lb = floor((fi_lb - 2 - f_shift(1))/2.)
1079 send_size = 0
1080 DO x = fi_lb, fi_ub
1081 p = pos_of_x(fine_gbo(1, 1) + modulo(x - fine_gbo(1, 1), sf))
1082 IF (p /= p_old) THEN
1083 p_ub = floor((x - 1 + 3 - f_shift(1))/2.)
1084
1085 send_size(p_old) = send_size(p_old) + (min(p_ub, my_coarse_bo(2, 1)) &
1086 - max(p_lb, my_coarse_bo(1, 1)) + 1)*coarse_slice_size
1087
1088 IF (pbc) THEN
1089 DO xx = p_lb, coarse_gbo(1, 1) - 1
1090 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1091 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1092 send_size(p_old) = send_size(p_old) + coarse_slice_size
1093 END IF
1094 END DO
1095 DO xx = coarse_gbo(2, 1) + 1, p_ub
1096 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1097 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1098 send_size(p_old) = send_size(p_old) + coarse_slice_size
1099 END IF
1100 END DO
1101 END IF
1102
1103 p_old = p
1104 p_lb = floor((x - 2 - f_shift(1))/2.)
1105 END IF
1106 END DO
1107 p_ub = floor((fi_ub + 3 - f_shift(1))/2.)
1108
1109 send_size(p_old) = send_size(p_old) + (min(p_ub, my_coarse_bo(2, 1)) &
1110 - max(p_lb, my_coarse_bo(1, 1)) + 1)*coarse_slice_size
1111
1112 IF (pbc) THEN
1113 DO xx = p_lb, coarse_gbo(1, 1) - 1
1114 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1115 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1116 send_size(p_old) = send_size(p_old) + coarse_slice_size
1117 END IF
1118 END DO
1119 DO xx = coarse_gbo(2, 1) + 1, p_ub
1120 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1121 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1122 send_size(p_old) = send_size(p_old) + coarse_slice_size
1123 END IF
1124 END DO
1125 END IF
1126 ! ** offsets & alloc send-rcv
1127
1128 send_tot_size = 0
1129 DO ip = 0, n_procs - 1
1130 send_offset(ip) = send_tot_size
1131 send_tot_size = send_tot_size + send_size(ip)
1132 END DO
1133 ALLOCATE (send_buf(0:send_tot_size - 1))
1134
1135 rcv_tot_size = 0
1136 DO ip = 0, n_procs - 1
1137 rcv_offset(ip) = rcv_tot_size
1138 rcv_tot_size = rcv_tot_size + rcv_size(ip)
1139 END DO
1140 IF (.NOT. rcv_tot_size == (coarse_bo(2, 1) - coarse_bo(1, 1) + 1)*coarse_slice_size) THEN
1141 cpabort("Error calculating rcv_tot_size ")
1142 END IF
1143 ALLOCATE (rcv_buf(0:rcv_tot_size - 1))
1144
1145 ! ** fill send buffer
1146
1147 p_old = pos_of_x(fine_gbo(1, 1) + modulo(fi_lb - fine_gbo(1, 1), sf))
1148 p_lb = floor((fi_lb - 2 - f_shift(1))/2.)
1149 sent_size(:) = send_offset
1150 ss = my_coarse_bo(2, 1) - my_coarse_bo(1, 1) + 1
1151 DO x = fi_lb, fi_ub
1152 p = pos_of_x(fine_gbo(1, 1) + modulo(x - fine_gbo(1, 1), sf))
1153 IF (p /= p_old) THEN
1154 shift = floor((fine_gbo(1, 1) + modulo(x - 1 - fine_gbo(1, 1), sf) - f_shift(1))/2._dp) - &
1155 floor((x - 1 - f_shift(1))/2._dp)
1156 p_ub = floor((x - 1 + 3 - f_shift(1))/2._dp)
1157
1158 IF (pbc) THEN
1159 DO xx = p_lb + shift, coarse_gbo(1, 1) - 1
1160 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), sf)
1161 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1162 CALL dcopy(coarse_slice_size, &
1163 coarse_coeffs(x_att, my_coarse_bo(1, 2), &
1164 my_coarse_bo(1, 3)), ss, send_buf(sent_size(p_old)), 1)
1165 sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1166 END IF
1167 END DO
1168 END IF
1169
1170 ii = sent_size(p_old)
1171 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1172 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1173 DO i = max(p_lb + shift, my_coarse_bo(1, 1)), min(p_ub + shift, my_coarse_bo(2, 1))
1174 send_buf(ii) = coarse_coeffs(i, j, k)
1175 ii = ii + 1
1176 END DO
1177 END DO
1178 END DO
1179 sent_size(p_old) = ii
1180
1181 IF (pbc) THEN
1182 DO xx = coarse_gbo(2, 1) + 1, p_ub + shift
1183 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1184 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1185 CALL dcopy(coarse_slice_size, &
1186 coarse_coeffs(x_att, my_coarse_bo(1, 2), &
1187 my_coarse_bo(1, 3)), ss, &
1188 send_buf(sent_size(p_old)), 1)
1189 sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1190 END IF
1191 END DO
1192 END IF
1193
1194 p_old = p
1195 p_lb = floor((x - 2 - f_shift(1))/2.)
1196 END IF
1197 END DO
1198 shift = floor((fine_gbo(1, 1) + modulo(x - 1 - fine_gbo(1, 1), sf) - f_shift(1))/2._dp) - &
1199 floor((x - 1 - f_shift(1))/2._dp)
1200 p_ub = floor((fi_ub + 3 - f_shift(1))/2.)
1201
1202 IF (pbc) THEN
1203 DO xx = p_lb + shift, coarse_gbo(1, 1) - 1
1204 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1205 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1206 CALL dcopy(coarse_slice_size, &
1207 coarse_coeffs(x_att, my_coarse_bo(1, 2), &
1208 my_coarse_bo(1, 3)), ss, send_buf(sent_size(p_old)), 1)
1209 sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1210 END IF
1211 END DO
1212 END IF
1213
1214 ii = sent_size(p_old)
1215 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1216 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1217 DO i = max(p_lb + shift, my_coarse_bo(1, 1)), min(p_ub + shift, my_coarse_bo(2, 1))
1218 send_buf(ii) = coarse_coeffs(i, j, k)
1219 ii = ii + 1
1220 END DO
1221 END DO
1222 END DO
1223 sent_size(p_old) = ii
1224
1225 IF (pbc) THEN
1226 DO xx = coarse_gbo(2, 1) + 1, p_ub + shift
1227 x_att = coarse_gbo(1, 1) + modulo(xx - coarse_gbo(1, 1), s(1))
1228 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
1229 CALL dcopy(coarse_slice_size, &
1230 coarse_coeffs(x_att, my_coarse_bo(1, 2), &
1231 my_coarse_bo(1, 3)), ss, send_buf(sent_size(p_old)), 1)
1232 sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1233 END IF
1234 END DO
1235 END IF
1236
1237 cpassert(all(sent_size(:n_procs - 2) == send_offset(1:)))
1238 cpassert(sent_size(n_procs - 1) == send_tot_size)
1239 ! test send/rcv sizes
1240 CALL coarse_coeffs_pw%pw_grid%para%group%alltoall(send_size, real_rcv_size, 1)
1241 cpassert(all(real_rcv_size == rcv_size))
1242 ! all2all
1243 CALL coarse_coeffs_pw%pw_grid%para%group%alltoall(sb=send_buf, scount=send_size, sdispl=send_offset, &
1244 rb=rcv_buf, rcount=rcv_size, rdispl=rcv_offset)
1245
1246 ! ** reorder rcv buffer
1247 ! (actually reordering should be needed only with pbc)
1248
1249 ALLOCATE (coarse_coeffs(coarse_bo(1, 1):coarse_bo(2, 1), &
1250 coarse_bo(1, 2):coarse_bo(2, 2), &
1251 coarse_bo(1, 3):coarse_bo(2, 3)))
1252
1253 my_lb = max(coarse_gbo(1, 1), coarse_bo(1, 1))
1254 my_ub = min(coarse_gbo(2, 1), coarse_bo(2, 1))
1255 pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x
1256 sent_size(:) = rcv_offset
1257 ss = coarse_bo(2, 1) - coarse_bo(1, 1) + 1
1258 DO x = my_ub + 1, coarse_bo(2, 1)
1259 p_old = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
1260 CALL dcopy(coarse_slice_size, &
1261 rcv_buf(sent_size(p_old)), 1, &
1262 coarse_coeffs(x, coarse_bo(1, 2), &
1263 coarse_bo(1, 3)), ss)
1264 sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1265 END DO
1266 p_old = pos_of_x(coarse_gbo(1, 1) &
1267 + modulo(my_lb - coarse_gbo(1, 1), s(1)))
1268 p_lb = my_lb
1269 DO x = my_lb, my_ub
1270 p = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
1271 IF (p /= p_old) THEN
1272 p_ub = x - 1
1273
1274 ii = sent_size(p_old)
1275 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1276 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1277 DO i = p_lb, p_ub
1278 coarse_coeffs(i, j, k) = rcv_buf(ii)
1279 ii = ii + 1
1280 END DO
1281 END DO
1282 END DO
1283 sent_size(p_old) = ii
1284
1285 p_lb = x
1286 p_old = p
1287 END IF
1288 rcv_size(p) = rcv_size(p) + coarse_slice_size
1289 END DO
1290 p_ub = my_ub
1291 ii = sent_size(p_old)
1292 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1293 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1294 DO i = p_lb, p_ub
1295 coarse_coeffs(i, j, k) = rcv_buf(ii)
1296 ii = ii + 1
1297 END DO
1298 END DO
1299 END DO
1300 sent_size(p_old) = ii
1301 DO x = coarse_bo(1, 1), my_lb - 1
1302 p_old = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
1303 CALL dcopy(coarse_slice_size, &
1304 rcv_buf(sent_size(p_old)), 1, &
1305 coarse_coeffs(x, coarse_bo(1, 2), &
1306 coarse_bo(1, 3)), ss)
1307 sent_size(p_old) = sent_size(p_old) + coarse_slice_size
1308 END DO
1309
1310 cpassert(all(sent_size(0:n_procs - 2) == rcv_offset(1:)))
1311 cpassert(sent_size(n_procs - 1) == rcv_tot_size)
1312
1313 ! dealloc
1314 DEALLOCATE (send_size, send_offset, rcv_size, rcv_offset)
1315 DEALLOCATE (send_buf, rcv_buf, real_rcv_size)
1316 CALL timestop(handle2)
1317
1318 END IF
1319 fine_values => fine_values_pw%array
1320 w_0 = [weights_1d(3), weights_1d(1), weights_1d(3)]
1321 w_1 = [weights_1d(4), weights_1d(2), weights_1d(2), weights_1d(4)]
1322
1323 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1324 DO ik = -3, 3
1325 IF (pbc) THEN
1326 wk = weights_1d(abs(ik) + 1)
1327 fk = fine_gbo(1, 3) + modulo(2*k + ik - fine_gbo(1, 3) + f_shift(3), 2*s(3))
1328 ELSE
1329 fk = 2*k + ik + f_shift(3)
1330 IF (fk <= fine_bo(1, 3) + 1 .OR. fk >= fine_bo(2, 3) - 1) THEN
1331 IF (fk < fine_bo(1, 3) .OR. fk > fine_bo(2, 3)) cycle
1332 IF (fk == fine_bo(1, 3) .OR. fk == fine_bo(2, 3)) THEN
1333 IF (ik /= 0) cycle
1334 wk = w_border0
1335 ELSE IF (fk == 2*coarse_bo(1, 3) + 1 + f_shift(3)) THEN
1336 SELECT CASE (ik)
1337 CASE (1)
1338 wk = w_border1(1)
1339 CASE (-1)
1340 wk = w_border1(2)
1341 CASE (-3)
1342 wk = w_border1(3)
1343 CASE default
1344 cpabort("Only 1, -1, -3 are supported as the value of ik")
1345 cycle
1346 END SELECT
1347 ELSE
1348 SELECT CASE (ik)
1349 CASE (3)
1350 wk = w_border1(3)
1351 CASE (1)
1352 wk = w_border1(2)
1353 CASE (-1)
1354 wk = w_border1(1)
1355 CASE default
1356 cpabort("Only 3, 1, -1 are supported as the value of ik")
1357 cycle
1358 END SELECT
1359 END IF
1360 ELSE
1361 wk = weights_1d(abs(ik) + 1)
1362 END IF
1363 END IF
1364 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1365 DO ij = -3, 3
1366 IF (pbc) THEN
1367 wj = weights_1d(abs(ij) + 1)*wk
1368 fj = fine_gbo(1, 2) + modulo(2*j + ij - fine_gbo(1, 2) + f_shift(2), 2*s(2))
1369 ELSE
1370 fj = 2*j + ij + f_shift(2)
1371 IF (fj <= fine_bo(1, 2) + 1 .OR. fj >= fine_bo(2, 2) - 1) THEN
1372 IF (fj < fine_bo(1, 2) .OR. fj > fine_bo(2, 2)) cycle
1373 IF (fj == fine_bo(1, 2) .OR. fj == fine_bo(2, 2)) THEN
1374 IF (ij /= 0) cycle
1375 wj = w_border0*wk
1376 ELSE IF (fj == 2*coarse_bo(1, 2) + 1 + f_shift(2)) THEN
1377 SELECT CASE (ij)
1378 CASE (1)
1379 wj = w_border1(1)*wk
1380 CASE (-1)
1381 wj = w_border1(2)*wk
1382 CASE (-3)
1383 wj = w_border1(3)*wk
1384 CASE default
1385 cycle
1386 END SELECT
1387 ELSE
1388 SELECT CASE (ij)
1389 CASE (-1)
1390 wj = w_border1(1)*wk
1391 CASE (1)
1392 wj = w_border1(2)*wk
1393 CASE (3)
1394 wj = w_border1(3)*wk
1395 CASE default
1396 cycle
1397 END SELECT
1398 END IF
1399 ELSE
1400 wj = weights_1d(abs(ij) + 1)*wk
1401 END IF
1402 END IF
1403
1404 IF (fine_bo(2, 1) - fine_bo(1, 1) < 7 .OR. safe_calc) THEN
1405! CALL timeset(routineN//"_safe",handle2)
1406 DO i = coarse_bo(1, 1), coarse_bo(2, 1)
1407 DO ii = -3, 3
1408 IF (pbc .AND. .NOT. is_split) THEN
1409 wi = weights_1d(abs(ii) + 1)*wj
1410 fi = fine_gbo(1, 1) + modulo(2*i + ii - fine_gbo(1, 1) + f_shift(1), 2*s(1))
1411 ELSE
1412 fi = 2*i + ii + f_shift(1)
1413 IF (fi < fine_bo(1, 1) .OR. fi > fine_bo(2, 1)) cycle
1414 IF (.NOT. pbc .AND. (fi <= fine_gbo(1, 1) + 1 .OR. &
1415 fi >= fine_gbo(2, 1) - 1)) THEN
1416 IF (fi == fine_gbo(1, 1) .OR. fi == fine_gbo(2, 1)) THEN
1417 IF (ii /= 0) cycle
1418 wi = w_border0*wj
1419 ELSE IF (fi == fine_gbo(1, 1) + 1) THEN
1420 SELECT CASE (ii)
1421 CASE (1)
1422 wi = w_border1(1)*wj
1423 CASE (-1)
1424 wi = w_border1(2)*wj
1425 CASE (-3)
1426 wi = w_border1(3)*wj
1427 CASE default
1428 cycle
1429 END SELECT
1430 ELSE
1431 SELECT CASE (ii)
1432 CASE (-1)
1433 wi = w_border1(1)*wj
1434 CASE (1)
1435 wi = w_border1(2)*wj
1436 CASE (3)
1437 wi = w_border1(3)*wj
1438 CASE default
1439 cycle
1440 END SELECT
1441 END IF
1442 ELSE
1443 wi = weights_1d(abs(ii) + 1)*wj
1444 END IF
1445 END IF
1446 fine_values(fi, fj, fk) = &
1447 fine_values(fi, fj, fk) + &
1448 wi*coarse_coeffs(i, j, k)
1449 END DO
1450 END DO
1451! CALL timestop(handle2)
1452 ELSE
1453! CALL timeset(routineN//"_core1",handle2)
1454 ww0 = wj*w_0
1455 ww1 = wj*w_1
1456 IF (pbc .AND. .NOT. is_split) THEN
1457 v3 = coarse_coeffs(coarse_bo(2, 1), j, k)
1458 i = coarse_bo(1, 1)
1459 fi = 2*i + f_shift(1)
1460 v0 = coarse_coeffs(i, j, k)
1461 v1 = coarse_coeffs(i + 1, j, k)
1462 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1463 ww0(1)*v3 + ww0(2)*v0 + ww0(3)*v1
1464 v2 = coarse_coeffs(i + 2, j, k)
1465 fi = fi + 1
1466 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1467 ww1(1)*v3 + ww1(2)*v0 + ww1(3)*v1 + ww1(4)*v2
1468 ELSE IF (.NOT. has_i_lbound) THEN
1469 i = coarse_bo(1, 1)
1470 fi = 2*i + f_shift(1)
1471 v0 = coarse_coeffs(i, j, k)
1472 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1473 w_border0*wj*v0
1474 v1 = coarse_coeffs(i + 1, j, k)
1475 v2 = coarse_coeffs(i + 2, j, k)
1476 fi = fi + 1
1477 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1478 wj*(w_border1(1)*v0 + w_border1(2)*v1 + &
1479 w_border1(3)*v2)
1480 ELSE
1481 i = coarse_bo(1, 1)
1482 v0 = coarse_coeffs(i, j, k)
1483 v1 = coarse_coeffs(i + 1, j, k)
1484 v2 = coarse_coeffs(i + 2, j, k)
1485 fi = 2*i + f_shift(1) + 1
1486 IF (.NOT. (fi + 1 == fine_bo(1, 1) .OR. &
1487 fi + 2 == fine_bo(1, 1))) THEN
1488 CALL cp_abort(__location__, &
1489 "unexpected start index "// &
1490 trim(cp_to_string(coarse_bo(1, 1)))//" "// &
1491 trim(cp_to_string(fi)))
1492 END IF
1493 END IF
1494 fi = fi + 1
1495 IF (fi >= fine_bo(1, 1)) THEN
1496 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1497 ww0(1)*v0 + ww0(2)*v1 + &
1498 ww0(3)*v2
1499 ELSE
1500 cpassert(fi + 1 == fine_bo(1, 1))
1501 END IF
1502! CALL timestop(handle2)
1503! CALL timeset(routineN//"_core",handle2)
1504 DO i = coarse_bo(1, 1) + 3, floor((fine_bo(2, 1) - f_shift(1))/2.) - 3, 4
1505 v3 = coarse_coeffs(i, j, k)
1506 fi = fi + 1
1507 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1508 (ww1(1)*v0 + ww1(2)*v1 + &
1509 ww1(3)*v2 + ww1(4)*v3)
1510 fi = fi + 1
1511 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1512 (ww0(1)*v1 + ww0(2)*v2 + &
1513 ww0(3)*v3)
1514 v0 = coarse_coeffs(i + 1, j, k)
1515 fi = fi + 1
1516 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1517 (ww1(4)*v0 + ww1(1)*v1 + &
1518 ww1(2)*v2 + ww1(3)*v3)
1519 fi = fi + 1
1520 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1521 (ww0(1)*v2 + ww0(2)*v3 + &
1522 ww0(3)*v0)
1523 v1 = coarse_coeffs(i + 2, j, k)
1524 fi = fi + 1
1525 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1526 (ww1(3)*v0 + ww1(4)*v1 + &
1527 ww1(1)*v2 + ww1(2)*v3)
1528 fi = fi + 1
1529 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1530 (ww0(1)*v3 + ww0(2)*v0 + &
1531 ww0(3)*v1)
1532 v2 = coarse_coeffs(i + 3, j, k)
1533 fi = fi + 1
1534 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1535 (ww1(2)*v0 + ww1(3)*v1 + &
1536 ww1(4)*v2 + ww1(1)*v3)
1537 fi = fi + 1
1538 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1539 (ww0(1)*v0 + ww0(2)*v1 + &
1540 ww0(3)*v2)
1541 END DO
1542! CALL timestop(handle2)
1543! CALL timeset(routineN//"_clean",handle2)
1544 rest_b = modulo(floor((fine_bo(2, 1) - f_shift(1))/2.) - coarse_bo(1, 1) - 3 + 1, 4)
1545 IF (rest_b > 0) THEN
1546 i = floor((fine_bo(2, 1) - f_shift(1))/2.) - rest_b + 1
1547 v3 = coarse_coeffs(i, j, k)
1548 fi = fi + 1
1549 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1550 (ww1(1)*v0 + ww1(2)*v1 + &
1551 ww1(3)*v2 + ww1(4)*v3)
1552 fi = fi + 1
1553 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1554 (ww0(1)*v1 + ww0(2)*v2 + &
1555 ww0(3)*v3)
1556 IF (rest_b > 1) THEN
1557 v0 = coarse_coeffs(i + 1, j, k)
1558 fi = fi + 1
1559 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1560 (ww1(4)*v0 + ww1(1)*v1 + &
1561 ww1(2)*v2 + ww1(3)*v3)
1562 fi = fi + 1
1563 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1564 (ww0(1)*v2 + ww0(2)*v3 + &
1565 ww0(3)*v0)
1566 IF (rest_b > 2) THEN
1567 v1 = coarse_coeffs(i + 2, j, k)
1568 fi = fi + 1
1569 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1570 (ww1(3)*v0 + ww1(4)*v1 + &
1571 ww1(1)*v2 + ww1(2)*v3)
1572 fi = fi + 1
1573 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1574 (ww0(1)*v3 + ww0(2)*v0 + &
1575 ww0(3)*v1)
1576 IF (pbc .AND. .NOT. is_split) THEN
1577 v2 = coarse_coeffs(coarse_bo(1, 1), j, k)
1578 fi = fi + 1
1579 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1580 ww1(1)*v3 + ww1(2)*v0 + ww1(3)*v1 + ww1(4)*v2
1581 fi = fi + 1
1582 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1583 ww0(1)*v0 + ww0(2)*v1 + ww0(3)*v2
1584 v3 = coarse_coeffs(coarse_bo(1, 1) + 1, j, k)
1585 fi = fi + 1
1586 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1587 ww1(1)*v0 + ww1(2)*v1 + ww1(3)*v2 + ww1(4)*v3
1588 ELSE IF (has_i_ubound) THEN
1589 v2 = coarse_coeffs(i + 3, j, k)
1590 fi = fi + 1
1591 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1592 ww1(1)*v3 + ww1(2)*v0 + ww1(3)*v1 + ww1(4)*v2
1593 fi = fi + 1
1594 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1595 ww0(1)*v0 + ww0(2)*v1 + ww0(3)*v2
1596 IF (fi + 1 == fine_bo(2, 1)) THEN
1597 v3 = coarse_coeffs(i + 4, j, k)
1598 fi = fi + 1
1599 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1600 ww1(1)*v0 + ww1(2)*v1 + ww1(3)*v2 + ww1(4)*v3
1601 END IF
1602 ELSE
1603 fi = fi + 1
1604 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1605 wj*(w_border1(3)*v3 + w_border1(2)*v0 + &
1606 w_border1(1)*v1)
1607 fi = fi + 1
1608 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1609 w_border0*wj*v1
1610 END IF
1611 ELSE IF (pbc .AND. .NOT. is_split) THEN
1612 v1 = coarse_coeffs(coarse_bo(1, 1), j, k)
1613 fi = fi + 1
1614 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1615 ww1(1)*v2 + ww1(2)*v3 + ww1(3)*v0 + ww1(4)*v1
1616 fi = fi + 1
1617 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1618 ww0(1)*v3 + ww0(2)*v0 + ww0(3)*v1
1619 v2 = coarse_coeffs(coarse_bo(1, 1) + 1, j, k)
1620 fi = fi + 1
1621 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1622 ww1(1)*v3 + ww1(2)*v0 + ww1(3)*v1 + ww1(4)*v2
1623 ELSE IF (has_i_ubound) THEN
1624 v1 = coarse_coeffs(i + 2, j, k)
1625 fi = fi + 1
1626 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1627 ww1(1)*v2 + ww1(2)*v3 + ww1(3)*v0 + ww1(4)*v1
1628 fi = fi + 1
1629 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1630 ww0(1)*v3 + ww0(2)*v0 + ww0(3)*v1
1631 IF (fi + 1 == fine_bo(2, 1)) THEN
1632 v2 = coarse_coeffs(i + 3, j, k)
1633 fi = fi + 1
1634 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1635 ww1(1)*v3 + ww1(2)*v0 + ww1(3)*v1 + ww1(4)*v2
1636 END IF
1637 ELSE
1638 fi = fi + 1
1639 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1640 wj*(w_border1(3)*v2 + w_border1(2)*v3 + &
1641 w_border1(1)*v0)
1642 fi = fi + 1
1643 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1644 w_border0*wj*v0
1645 END IF
1646 ELSE IF (pbc .AND. .NOT. is_split) THEN
1647 v0 = coarse_coeffs(coarse_bo(1, 1), j, k)
1648 fi = fi + 1
1649 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1650 ww1(1)*v1 + ww1(2)*v2 + ww1(3)*v3 + ww1(4)*v0
1651 fi = fi + 1
1652 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1653 ww0(1)*v2 + ww0(2)*v3 + ww0(3)*v0
1654 v1 = coarse_coeffs(coarse_bo(1, 1) + 1, j, k)
1655 fi = fi + 1
1656 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1657 ww1(1)*v2 + ww1(2)*v3 + ww1(3)*v0 + ww1(4)*v1
1658 ELSE IF (has_i_ubound) THEN
1659 v0 = coarse_coeffs(i + 1, j, k)
1660 fi = fi + 1
1661 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1662 ww1(1)*v1 + ww1(2)*v2 + ww1(3)*v3 + ww1(4)*v0
1663 fi = fi + 1
1664 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1665 ww0(1)*v2 + ww0(2)*v3 + ww0(3)*v0
1666 IF (fi + 1 == fine_bo(2, 1)) THEN
1667 v1 = coarse_coeffs(i + 2, j, k)
1668 fi = fi + 1
1669 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1670 ww1(1)*v2 + ww1(2)*v3 + ww1(3)*v0 + ww1(4)*v1
1671 END IF
1672 ELSE
1673 fi = fi + 1
1674 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1675 wj*(w_border1(3)*v1 + w_border1(2)*v2 + &
1676 w_border1(1)*v3)
1677 fi = fi + 1
1678 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1679 w_border0*wj*v3
1680 END IF
1681 ELSE IF (pbc .AND. .NOT. is_split) THEN
1682 v3 = coarse_coeffs(coarse_bo(1, 1), j, k)
1683 fi = fi + 1
1684 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1685 ww1(1)*v0 + ww1(2)*v1 + ww1(3)*v2 + ww1(4)*v3
1686 fi = fi + 1
1687 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1688 ww0(1)*v1 + ww0(2)*v2 + ww0(3)*v3
1689 v0 = coarse_coeffs(coarse_bo(1, 1) + 1, j, k)
1690 fi = fi + 1
1691 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1692 ww1(1)*v1 + ww1(2)*v2 + ww1(3)*v3 + ww1(4)*v0
1693 ELSE IF (has_i_ubound) THEN
1694 v3 = coarse_coeffs(i, j, k)
1695 fi = fi + 1
1696 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1697 ww1(1)*v0 + ww1(2)*v1 + ww1(3)*v2 + ww1(4)*v3
1698 fi = fi + 1
1699 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1700 ww0(1)*v1 + ww0(2)*v2 + ww0(3)*v3
1701 IF (fi + 1 == fine_bo(2, 1)) THEN
1702 v0 = coarse_coeffs(i + 1, j, k)
1703 fi = fi + 1
1704 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1705 ww1(1)*v1 + ww1(2)*v2 + ww1(3)*v3 + ww1(4)*v0
1706 END IF
1707 ELSE
1708 fi = fi + 1
1709 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1710 wj*(w_border1(3)*v0 + w_border1(2)*v1 + &
1711 w_border1(1)*v2)
1712 fi = fi + 1
1713 fine_values(fi, fj, fk) = fine_values(fi, fj, fk) + &
1714 w_border0*wj*v2
1715 END IF
1716 cpassert(fi == fine_bo(2, 1))
1717 END IF
1718! CALL timestop(handle2)
1719 END DO
1720 END DO
1721 END DO
1722 END DO
1723
1724 IF (is_split) THEN
1725 DEALLOCATE (coarse_coeffs)
1726 END IF
1727 CALL timestop(handle)
1728 END SUBROUTINE add_coarse2fine
1729
1730! **************************************************************************************************
1731!> \brief low level function that adds a coarse grid (without boundary)
1732!> to a fine grid.
1733!>
1734!> It will add to
1735!>
1736!> coarse_coeffs(coarse_bounds(1,1):coarse_bounds(2,1),
1737!> coarse_bounds(1,2):coarse_bounds(2,2),
1738!> coarse_bounds(1,3):coarse_bounds(2,3))
1739!>
1740!> using
1741!>
1742!> fine_values(2*coarse_bounds(1,1):2*coarse_bounds(2,1),
1743!> 2*coarse_bounds(1,2):2*coarse_bounds(2,2),
1744!> 2*coarse_bounds(1,3):2*coarse_bounds(2,3))
1745!>
1746!> composed with the weights obtained by the direct product of the
1747!> 1d coefficients weights:
1748!>
1749!> for i,j,k in -3..3
1750!> w(i,j,k)=weights_1d(abs(i)+1)*weights_1d(abs(j)+1)*
1751!> weights_1d(abs(k)+1)
1752!> \param fine_values_pw 3d array where to add the values due to the
1753!> coarse coeffs
1754!> \param coarse_coeffs_pw 3d array with boundary of size 1 with the values of the
1755!> coefficients
1756!> \param weights_1d the weights of the 1d smearing
1757!> \param w_border0 the 1d weight at the border
1758!> \param w_border1 the 1d weights for a point one off the border
1759!> (w_border1(1) is the weight of the coefficent at the border)
1760!> \param pbc ...
1761!> \param safe_computation ...
1762!> \author fawzi
1763!> \note
1764!> see coarse2fine for some relevant notes
1765! **************************************************************************************************
1766 SUBROUTINE add_fine2coarse(fine_values_pw, coarse_coeffs_pw, &
1767 weights_1d, w_border0, w_border1, pbc, safe_computation)
1768 TYPE(pw_r3d_rs_type), INTENT(IN) :: fine_values_pw, coarse_coeffs_pw
1769 REAL(kind=dp), DIMENSION(4), INTENT(in) :: weights_1d
1770 REAL(kind=dp), INTENT(in) :: w_border0
1771 REAL(kind=dp), DIMENSION(3), INTENT(in) :: w_border1
1772 LOGICAL, INTENT(in) :: pbc
1773 LOGICAL, INTENT(in), OPTIONAL :: safe_computation
1774
1775 CHARACTER(len=*), PARAMETER :: routinen = 'add_fine2coarse'
1776
1777 INTEGER :: coarse_slice_size, f_shift(3), fi, fj, fk, handle, handle2, i, ii, ij, ik, ip, j, &
1778 k, n_procs, p, p_old, rcv_tot_size, rest_b, s(3), send_tot_size, ss, x, x_att
1779 INTEGER, ALLOCATABLE, DIMENSION(:) :: pp_lb, pp_ub, rcv_offset, rcv_size, &
1780 real_rcv_size, send_offset, send_size, &
1781 sent_size
1782 INTEGER, DIMENSION(2, 3) :: coarse_bo, coarse_gbo, fine_bo, &
1783 fine_gbo, my_coarse_bo
1784 INTEGER, DIMENSION(:), POINTER :: pos_of_x
1785 LOGICAL :: has_i_lbound, has_i_ubound, is_split, &
1786 local_data, safe_calc
1787 REAL(kind=dp) :: vv0, vv1, vv2, vv3, vv4, vv5, vv6, vv7, &
1788 wi, wj, wk
1789 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rcv_buf, send_buf
1790 REAL(kind=dp), DIMENSION(3) :: w_0, ww0
1791 REAL(kind=dp), DIMENSION(4) :: w_1, ww1
1792 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: coarse_coeffs, fine_values
1793
1794 CALL timeset(routinen, handle)
1795
1796 safe_calc = .false.
1797 IF (PRESENT(safe_computation)) safe_calc = safe_computation
1798
1799 my_coarse_bo = coarse_coeffs_pw%pw_grid%bounds_local
1800 coarse_gbo = coarse_coeffs_pw%pw_grid%bounds
1801 fine_bo = fine_values_pw%pw_grid%bounds_local
1802 fine_gbo = fine_values_pw%pw_grid%bounds
1803 f_shift = fine_gbo(1, :) - 2*coarse_gbo(1, :)
1804 is_split = any(coarse_gbo(:, 1) /= my_coarse_bo(:, 1))
1805 coarse_bo = my_coarse_bo
1806 IF (fine_bo(1, 1) <= fine_bo(2, 1)) THEN
1807 coarse_bo(1, 1) = floor(real(fine_bo(1, 1) - f_shift(1), dp)/2._dp) - 1
1808 coarse_bo(2, 1) = floor(real(fine_bo(2, 1) + 1 - f_shift(1), dp)/2._dp) + 1
1809 ELSE
1810 coarse_bo(1, 1) = coarse_gbo(2, 1)
1811 coarse_bo(2, 1) = coarse_gbo(2, 1) - 1
1812 END IF
1813 IF (.NOT. is_split .OR. .NOT. pbc) THEN
1814 coarse_bo(1, 1) = max(coarse_gbo(1, 1), coarse_bo(1, 1))
1815 coarse_bo(2, 1) = min(coarse_gbo(2, 1), coarse_bo(2, 1))
1816 END IF
1817 has_i_ubound = (fine_gbo(2, 1) /= fine_bo(2, 1)) .OR. pbc .AND. is_split
1818 has_i_lbound = (fine_gbo(1, 1) /= fine_bo(1, 1)) .OR. pbc .AND. is_split
1819
1820 IF (pbc) THEN
1821 cpassert(all(fine_gbo(1, :) == 2*coarse_gbo(1, :) + f_shift))
1822 cpassert(all(fine_gbo(2, :) == 2*coarse_gbo(2, :) + f_shift + 1))
1823 ELSE
1824 cpassert(all(fine_gbo(2, :) == 2*coarse_gbo(2, :) + f_shift))
1825 cpassert(all(fine_gbo(1, :) == 2*coarse_gbo(1, :) + f_shift))
1826 END IF
1827 cpassert(coarse_gbo(2, 1) - coarse_gbo(1, 2) > 1)
1828 local_data = is_split ! ANY(coarse_bo/=my_coarse_bo)
1829 IF (local_data) THEN
1830 ALLOCATE (coarse_coeffs(coarse_bo(1, 1):coarse_bo(2, 1), &
1831 coarse_bo(1, 2):coarse_bo(2, 2), &
1832 coarse_bo(1, 3):coarse_bo(2, 3)))
1833 coarse_coeffs = 0._dp
1834 ELSE
1835 coarse_coeffs => coarse_coeffs_pw%array
1836 END IF
1837
1838 fine_values => fine_values_pw%array
1839 w_0 = [weights_1d(3), weights_1d(1), weights_1d(3)]
1840 w_1 = [weights_1d(4), weights_1d(2), weights_1d(2), weights_1d(4)]
1841
1842 DO i = 1, 3
1843 s(i) = coarse_gbo(2, i) - coarse_gbo(1, i) + 1
1844 END DO
1845 IF (any(s < 1)) RETURN
1846
1847 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
1848 DO ik = -3, 3
1849 IF (pbc) THEN
1850 wk = weights_1d(abs(ik) + 1)
1851 fk = fine_gbo(1, 3) + modulo(2*k + ik - fine_gbo(1, 3) + f_shift(3), 2*s(3))
1852 ELSE
1853 fk = 2*k + ik + f_shift(3)
1854 IF (fk <= fine_bo(1, 3) + 1 .OR. fk >= fine_bo(2, 3) - 1) THEN
1855 IF (fk < fine_bo(1, 3) .OR. fk > fine_bo(2, 3)) cycle
1856 IF (fk == fine_bo(1, 3) .OR. fk == fine_bo(2, 3)) THEN
1857 IF (ik /= 0) cycle
1858 wk = w_border0
1859 ELSE IF (fk == fine_bo(1, 3) + 1) THEN
1860 SELECT CASE (ik)
1861 CASE (1)
1862 wk = w_border1(1)
1863 CASE (-1)
1864 wk = w_border1(2)
1865 CASE (-3)
1866 wk = w_border1(3)
1867 CASE default
1868 cpabort("Only 1, -1, -3 are supported as the value of ik")
1869 cycle
1870 END SELECT
1871 ELSE
1872 SELECT CASE (ik)
1873 CASE (3)
1874 wk = w_border1(3)
1875 CASE (1)
1876 wk = w_border1(2)
1877 CASE (-1)
1878 wk = w_border1(1)
1879 CASE default
1880 cpabort("Only 3, 1, -1 are supported as the value of ik")
1881 cycle
1882 END SELECT
1883 END IF
1884 ELSE
1885 wk = weights_1d(abs(ik) + 1)
1886 END IF
1887 END IF
1888 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
1889 DO ij = -3, 3
1890 IF (pbc) THEN
1891 fj = fine_gbo(1, 2) + modulo(2*j + ij - fine_gbo(1, 2) + f_shift(2), &
1892 2*s(2))
1893 wj = weights_1d(abs(ij) + 1)*wk
1894 ELSE
1895 fj = 2*j + ij + f_shift(2)
1896 IF (fj <= fine_bo(1, 2) + 1 .OR. fj >= fine_bo(2, 2) - 1) THEN
1897 IF (fj < fine_bo(1, 2) .OR. fj > fine_bo(2, 2)) cycle
1898 IF (fj == fine_bo(1, 2) .OR. fj == fine_bo(2, 2)) THEN
1899 IF (ij /= 0) cycle
1900 wj = w_border0*wk
1901 ELSE IF (fj == fine_bo(1, 2) + 1) THEN
1902 SELECT CASE (ij)
1903 CASE (1)
1904 wj = w_border1(1)*wk
1905 CASE (-1)
1906 wj = w_border1(2)*wk
1907 CASE (-3)
1908 wj = w_border1(3)*wk
1909 CASE default
1910 cpabort("Only 1, -1, -3 are supported as the value of ij")
1911 cycle
1912 END SELECT
1913 ELSE
1914 SELECT CASE (ij)
1915 CASE (-1)
1916 wj = w_border1(1)*wk
1917 CASE (1)
1918 wj = w_border1(2)*wk
1919 CASE (3)
1920 wj = w_border1(3)*wk
1921 CASE default
1922 cpabort("Only -1, 1, 3 are supported as the value of ij")
1923 cycle
1924 END SELECT
1925 END IF
1926 ELSE
1927 wj = weights_1d(abs(ij) + 1)*wk
1928 END IF
1929 END IF
1930
1931 IF (coarse_bo(2, 1) - coarse_bo(1, 1) < 7 .OR. safe_calc) THEN
1932 DO i = coarse_bo(1, 1), coarse_bo(2, 1)
1933 DO ii = -3, 3
1934 IF (pbc .AND. .NOT. is_split) THEN
1935 wi = weights_1d(abs(ii) + 1)*wj
1936 fi = fine_gbo(1, 1) + modulo(2*i + ii - fine_gbo(1, 1) + f_shift(1), 2*s(1))
1937 ELSE
1938 fi = 2*i + ii + f_shift(1)
1939 IF (fi < fine_bo(1, 1) .OR. fi > fine_bo(2, 1)) cycle
1940 IF (((.NOT. pbc) .AND. fi <= fine_gbo(1, 1) + 1) .OR. &
1941 ((.NOT. pbc) .AND. fi >= fine_gbo(2, 1) - 1)) THEN
1942 IF (fi == fine_gbo(1, 1) .OR. fi == fine_gbo(2, 1)) THEN
1943 IF (ii /= 0) cycle
1944 wi = w_border0*wj
1945 ELSE IF (fi == fine_gbo(1, 1) + 1) THEN
1946 SELECT CASE (ii)
1947 CASE (1)
1948 wi = w_border1(1)*wj
1949 CASE (-1)
1950 wi = w_border1(2)*wj
1951 CASE (-3)
1952 wi = w_border1(3)*wj
1953 CASE default
1954 cycle
1955 END SELECT
1956 ELSE
1957 SELECT CASE (ii)
1958 CASE (-1)
1959 wi = w_border1(1)*wj
1960 CASE (1)
1961 wi = w_border1(2)*wj
1962 CASE (3)
1963 wi = w_border1(3)*wj
1964 CASE default
1965 cycle
1966 END SELECT
1967 END IF
1968 ELSE
1969 wi = weights_1d(abs(ii) + 1)*wj
1970 END IF
1971 END IF
1972 coarse_coeffs(i, j, k) = &
1973 coarse_coeffs(i, j, k) + &
1974 wi*fine_values(fi, fj, fk)
1975 END DO
1976 END DO
1977 ELSE
1978 ww0 = wj*w_0
1979 ww1 = wj*w_1
1980 IF (pbc .AND. .NOT. is_split) THEN
1981 i = coarse_bo(1, 1) - 1
1982 vv2 = fine_values(fine_bo(2, 1) - 2, fj, fk)
1983 vv3 = fine_values(fine_bo(2, 1) - 1, fj, fk)
1984 vv4 = fine_values(fine_bo(2, 1), fj, fk)
1985 fi = fine_bo(1, 1)
1986 vv5 = fine_values(fi, fj, fk)
1987 fi = fi + 1
1988 vv6 = fine_values(fi, fj, fk)
1989 fi = fi + 1
1990 vv7 = fine_values(fi, fj, fk)
1991 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
1992 + ww1(4)*vv2 + ww0(3)*vv3 + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6 + ww0(1)*vv7
1993 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
1994 + ww1(4)*vv4 + ww0(3)*vv5 + ww1(3)*vv6 + ww0(2)*vv7
1995 coarse_coeffs(i + 3, j, k) = coarse_coeffs(i + 3, j, k) &
1996 + ww1(4)*vv6 + ww0(3)*vv7
1997 ELSE IF (has_i_lbound) THEN
1998 i = coarse_bo(1, 1)
1999 fi = fine_bo(1, 1) - 1
2000 IF (i + 1 == floor((fine_bo(1, 1) + 1 - f_shift(1))/2._dp)) THEN
2001 fi = fi + 1
2002 vv0 = fine_values(fi, fj, fk)
2003 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) + &
2004 vv0*ww0(3)
2005 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) + &
2006 vv0*ww0(2)
2007 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) + &
2008 vv0*ww0(1)
2009 END IF
2010 ELSE
2011 i = coarse_bo(1, 1)
2012 fi = 2*i + f_shift(1)
2013 vv0 = fine_values(fi, fj, fk)
2014 fi = fi + 1
2015 vv1 = fine_values(fi, fj, fk)
2016 fi = fi + 1
2017 vv2 = fine_values(fi, fj, fk)
2018 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) + &
2019 (vv0*w_border0 + vv1*w_border1(1))*wj + vv2*ww0(1)
2020 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) + &
2021 wj*w_border1(2)*vv1 + ww0(2)*vv2
2022 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) + &
2023 wj*w_border1(3)*vv1 + ww0(3)*vv2
2024 END IF
2025 DO i = coarse_bo(1, 1) + 3, floor((fine_bo(2, 1) - f_shift(1))/2._dp) - 3, 4
2026 fi = fi + 1
2027 vv0 = fine_values(fi, fj, fk)
2028 fi = fi + 1
2029 vv1 = fine_values(fi, fj, fk)
2030 fi = fi + 1
2031 vv2 = fine_values(fi, fj, fk)
2032 fi = fi + 1
2033 vv3 = fine_values(fi, fj, fk)
2034 fi = fi + 1
2035 vv4 = fine_values(fi, fj, fk)
2036 fi = fi + 1
2037 vv5 = fine_values(fi, fj, fk)
2038 fi = fi + 1
2039 vv6 = fine_values(fi, fj, fk)
2040 fi = fi + 1
2041 vv7 = fine_values(fi, fj, fk)
2042 coarse_coeffs(i - 3, j, k) = coarse_coeffs(i - 3, j, k) &
2043 + ww1(1)*vv0
2044 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2045 + ww1(2)*vv0 + ww0(1)*vv1 + ww1(1)*vv2
2046 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2047 + ww1(3)*vv0 + ww0(2)*vv1 + ww1(2)*vv2 + ww0(1)*vv3 + ww1(1)*vv4
2048 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2049 + ww1(4)*vv0 + ww0(3)*vv1 + ww1(3)*vv2 + ww0(2)*vv3 + ww1(2)*vv4 + ww0(1)*vv5 + ww1(1)*vv6
2050 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2051 + ww1(4)*vv2 + ww0(3)*vv3 + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6 + ww0(1)*vv7
2052 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2053 + ww1(4)*vv4 + ww0(3)*vv5 + ww1(3)*vv6 + ww0(2)*vv7
2054 coarse_coeffs(i + 3, j, k) = coarse_coeffs(i + 3, j, k) &
2055 + ww1(4)*vv6 + ww0(3)*vv7
2056 END DO
2057 IF (.NOT. floor((fine_bo(2, 1) - f_shift(1))/2._dp) - coarse_bo(1, 1) >= 4) THEN
2058 cpabort("FLOOR((fine_bo(2,1)-f_shift(1))/2._dp)-coarse_bo(1,1)>=4")
2059 END IF
2060 rest_b = modulo(floor((fine_bo(2, 1) - f_shift(1))/2._dp) - coarse_bo(1, 1) - 6, 4)
2061 i = floor((fine_bo(2, 1) - f_shift(1))/2._dp) - 3 - rest_b + 4
2062 cpassert(fi == (i - 2)*2 + f_shift(1))
2063 IF (rest_b > 0) THEN
2064 fi = fi + 1
2065 vv0 = fine_values(fi, fj, fk)
2066 fi = fi + 1
2067 vv1 = fine_values(fi, fj, fk)
2068 coarse_coeffs(i - 3, j, k) = coarse_coeffs(i - 3, j, k) &
2069 + ww1(1)*vv0
2070 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2071 + ww1(2)*vv0 + ww0(1)*vv1
2072 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2073 + ww1(3)*vv0 + ww0(2)*vv1
2074 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2075 + ww1(4)*vv0 + ww0(3)*vv1
2076 IF (rest_b > 1) THEN
2077 fi = fi + 1
2078 vv2 = fine_values(fi, fj, fk)
2079 fi = fi + 1
2080 vv3 = fine_values(fi, fj, fk)
2081 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2082 + ww1(1)*vv2
2083 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2084 + ww1(2)*vv2 + ww0(1)*vv3
2085 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2086 + ww1(3)*vv2 + ww0(2)*vv3
2087 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2088 + ww1(4)*vv2 + ww0(3)*vv3
2089 IF (rest_b > 2) THEN
2090 fi = fi + 1
2091 vv4 = fine_values(fi, fj, fk)
2092 fi = fi + 1
2093 vv5 = fine_values(fi, fj, fk)
2094 fi = fi + 1
2095 vv6 = fine_values(fi, fj, fk)
2096 fi = fi + 1
2097 vv7 = fine_values(fi, fj, fk)
2098 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2099 + ww1(1)*vv4
2100 IF (has_i_ubound) THEN
2101 IF (coarse_bo(2, 1) - 2 == floor((fine_bo(2, 1) - f_shift(1))/2._dp)) THEN
2102 fi = fi + 1
2103 vv0 = fine_values(fi, fj, fk)
2104 coarse_coeffs(i + 4, j, k) = coarse_coeffs(i + 4, j, k) &
2105 + vv0*ww1(4)
2106 ELSE
2107 vv0 = 0._dp
2108 END IF
2109 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2110 + ww1(2)*vv4 + ww0(1)*vv5 + ww1(1)*vv6
2111 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2112 + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6 + ww0(1)*vv7 + vv0*ww1(1)
2113 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2114 + ww1(4)*vv4 + ww0(3)*vv5 + ww1(3)*vv6 + ww0(2)*vv7 + vv0*ww1(2)
2115 coarse_coeffs(i + 3, j, k) = coarse_coeffs(i + 3, j, k) &
2116 + ww1(4)*vv6 + ww0(3)*vv7 + vv0*ww1(3)
2117 ELSE IF (pbc .AND. .NOT. is_split) THEN
2118 fi = fi + 1
2119 vv0 = fine_values(fi, fj, fk)
2120 vv1 = fine_values(fine_bo(1, 1), fj, fk)
2121 vv2 = fine_values(fine_bo(1, 1) + 1, fj, fk)
2122 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2123 + ww1(2)*vv4 + ww0(1)*vv5 + ww1(1)*vv6
2124 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2125 + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6 + ww0(1)*vv7 + vv0*ww1(1)
2126 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2127 + ww1(4)*vv4 + ww0(3)*vv5 + ww1(3)*vv6 + ww0(2)*vv7 + vv0*ww1(2) &
2128 + vv1*ww0(1) + vv2*ww1(1)
2129 ELSE
2130 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2131 + ww1(2)*vv4 + ww0(1)*vv5 + wj*w_border1(3)*vv6
2132 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2133 + ww1(3)*vv4 + ww0(2)*vv5 + wj*w_border1(2)*vv6
2134 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2135 + ww1(4)*vv4 + ww0(3)*vv5 + wj*w_border1(1)*vv6 + w_border0*wj*vv7
2136 END IF
2137 ELSE
2138 fi = fi + 1
2139 vv4 = fine_values(fi, fj, fk)
2140 fi = fi + 1
2141 vv5 = fine_values(fi, fj, fk)
2142 IF (has_i_ubound) THEN
2143 IF (coarse_bo(2, 1) - 2 == floor((fine_bo(2, 1) - f_shift(1))/2._dp)) THEN
2144 fi = fi + 1
2145 vv6 = fine_values(fi, fj, fk)
2146 coarse_coeffs(i + 3, j, k) = coarse_coeffs(i + 3, j, k) &
2147 + ww1(4)*vv6
2148 ELSE
2149 vv6 = 0._dp
2150 END IF
2151 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2152 + ww1(1)*vv4
2153 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2154 + ww1(2)*vv4 + ww0(1)*vv5 + ww1(1)*vv6
2155 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2156 + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6
2157 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2158 + ww1(4)*vv4 + ww0(3)*vv5 + ww1(3)*vv6
2159 ELSE IF (pbc .AND. .NOT. is_split) THEN
2160 fi = fi + 1
2161 vv6 = fine_values(fi, fj, fk)
2162 vv7 = fine_values(fine_bo(1, 1), fj, fk)
2163 vv0 = fine_values(fine_bo(1, 1) + 1, fj, fk)
2164 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2165 + ww1(1)*vv4
2166 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2167 + ww1(4)*vv0 + ww0(3)*vv1 + ww1(3)*vv2 + ww0(2)*vv3 + &
2168 ww1(2)*vv4 + ww0(1)*vv5 + ww1(1)*vv6
2169 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2170 + ww1(4)*vv2 + ww0(3)*vv3 + ww1(3)*vv4 + ww0(2)*vv5 + ww1(2)*vv6 &
2171 + ww0(1)*vv7 + ww1(1)*vv0
2172 ELSE
2173 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2174 + wj*w_border1(3)*vv4
2175 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2176 + wj*w_border1(2)*vv4
2177 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2178 + wj*(w_border1(1)*vv4 + w_border0*vv5)
2179 END IF
2180 END IF
2181 ELSE
2182 fi = fi + 1
2183 vv2 = fine_values(fi, fj, fk)
2184 fi = fi + 1
2185 vv3 = fine_values(fi, fj, fk)
2186 IF (has_i_ubound) THEN
2187 IF (coarse_bo(2, 1) - 2 == floor((fine_bo(2, 1) - f_shift(1))/2._dp)) THEN
2188 fi = fi + 1
2189 vv4 = fine_values(fi, fj, fk)
2190 coarse_coeffs(i + 2, j, k) = coarse_coeffs(i + 2, j, k) &
2191 + ww1(4)*vv4
2192 ELSE
2193 vv4 = 0._dp
2194 END IF
2195 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2196 + ww1(1)*vv2
2197 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2198 + ww1(2)*vv2 + ww0(1)*vv3 + ww1(1)*vv4
2199 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2200 + ww1(3)*vv2 + ww0(2)*vv3 + ww1(2)*vv4
2201 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2202 + ww1(4)*vv2 + ww0(3)*vv3 + ww1(3)*vv4
2203 ELSE IF (pbc .AND. .NOT. is_split) THEN
2204 fi = fi + 1
2205 vv4 = fine_values(fi, fj, fk)
2206 vv5 = fine_values(fine_bo(1, 1), fj, fk)
2207 vv6 = fine_values(fine_bo(1, 1) + 1, fj, fk)
2208 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2209 + ww1(1)*vv2
2210 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2211 + ww1(2)*vv2 + ww0(1)*vv3 + ww1(1)*vv4
2212 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2213 + ww1(3)*vv2 + ww0(2)*vv3 + ww1(2)*vv4 + vv5*ww0(1) + ww1(1)*vv6
2214 ELSE
2215 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2216 + wj*w_border1(3)*vv2
2217 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2218 + wj*w_border1(2)*vv2
2219 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2220 + wj*(w_border1(1)*vv2 + w_border0*vv3)
2221 END IF
2222 END IF
2223 ELSE
2224 fi = fi + 1
2225 vv0 = fine_values(fi, fj, fk)
2226 fi = fi + 1
2227 vv1 = fine_values(fi, fj, fk)
2228 IF (has_i_ubound) THEN
2229 IF (coarse_bo(2, 1) - 2 == floor((fine_bo(2, 1) - f_shift(1))/2._dp)) THEN
2230 fi = fi + 1
2231 vv2 = fine_values(fi, fj, fk)
2232 coarse_coeffs(i + 1, j, k) = coarse_coeffs(i + 1, j, k) &
2233 + ww1(4)*vv2
2234 ELSE
2235 vv2 = 0._dp
2236 END IF
2237 coarse_coeffs(i - 3, j, k) = coarse_coeffs(i - 3, j, k) &
2238 + ww1(1)*vv0
2239 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2240 + ww1(2)*vv0 + ww0(1)*vv1 + ww1(1)*vv2
2241 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2242 + ww1(3)*vv0 + ww0(2)*vv1 + ww1(2)*vv2
2243 coarse_coeffs(i, j, k) = coarse_coeffs(i, j, k) &
2244 + ww1(4)*vv0 + ww0(3)*vv1 + ww1(3)*vv2
2245 ELSE IF (pbc .AND. .NOT. is_split) THEN
2246 fi = fi + 1
2247 vv2 = fine_values(fi, fj, fk)
2248 vv3 = fine_values(fine_bo(1, 1), fk, fk)
2249 vv4 = fine_values(fine_bo(1, 1) + 1, fk, fk)
2250 coarse_coeffs(i - 3, j, k) = coarse_coeffs(i - 3, j, k) &
2251 + ww1(1)*vv0
2252 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2253 + ww1(2)*vv0 + ww0(1)*vv1 + ww1(1)*vv2
2254 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2255 + ww1(3)*vv0 + ww0(2)*vv1 + ww1(2)*vv2 + ww0(1)*vv3 + ww1(1)*vv4
2256 ELSE
2257 coarse_coeffs(i - 3, j, k) = coarse_coeffs(i - 3, j, k) &
2258 + wj*w_border1(3)*vv0
2259 coarse_coeffs(i - 2, j, k) = coarse_coeffs(i - 2, j, k) &
2260 + wj*w_border1(2)*vv0
2261 coarse_coeffs(i - 1, j, k) = coarse_coeffs(i - 1, j, k) &
2262 + wj*(w_border1(1)*vv0 + w_border0*vv1)
2263 END IF
2264 END IF
2265 cpassert(fi == fine_bo(2, 1))
2266 END IF
2267 END DO
2268 END DO
2269 END DO
2270 END DO
2271
2272 ! *** parallel case
2273 IF (is_split) THEN
2274 CALL timeset(routinen//"_comm", handle2)
2275 coarse_slice_size = (coarse_bo(2, 2) - coarse_bo(1, 2) + 1)* &
2276 (coarse_bo(2, 3) - coarse_bo(1, 3) + 1)
2277 n_procs = coarse_coeffs_pw%pw_grid%para%group%num_pe
2278 ALLOCATE (send_size(0:n_procs - 1), send_offset(0:n_procs - 1), &
2279 sent_size(0:n_procs - 1), rcv_size(0:n_procs - 1), &
2280 rcv_offset(0:n_procs - 1), pp_lb(0:n_procs - 1), &
2281 pp_ub(0:n_procs - 1), real_rcv_size(0:n_procs - 1))
2282
2283 ! ** send size count
2284
2285 pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x
2286 send_size = 0
2287 DO x = coarse_bo(1, 1), coarse_bo(2, 1)
2288 p = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
2289 send_size(p) = send_size(p) + coarse_slice_size
2290 END DO
2291
2292 ! ** rcv size count
2293
2294 pos_of_x => fine_values_pw%pw_grid%para%pos_of_x
2295 p_old = pos_of_x(fine_gbo(1, 1))
2296 pp_lb = fine_gbo(2, 1)
2297 pp_ub = fine_gbo(2, 1) - 1
2298 pp_lb(p_old) = fine_gbo(1, 1)
2299 DO x = fine_gbo(1, 1), fine_gbo(2, 1)
2300 p = pos_of_x(x)
2301 IF (p /= p_old) THEN
2302 pp_ub(p_old) = x - 1
2303 pp_lb(p) = x
2304 p_old = p
2305 END IF
2306 END DO
2307 pp_ub(p_old) = fine_gbo(2, 1)
2308
2309 DO ip = 0, n_procs - 1
2310 IF (pp_lb(ip) <= pp_ub(ip)) THEN
2311 pp_lb(ip) = floor(real(pp_lb(ip) - f_shift(1), dp)/2._dp) - 1
2312 pp_ub(ip) = floor(real(pp_ub(ip) + 1 - f_shift(1), dp)/2._dp) + 1
2313 ELSE
2314 pp_lb(ip) = coarse_gbo(2, 1)
2315 pp_ub(ip) = coarse_gbo(2, 1) - 1
2316 END IF
2317 IF (.NOT. is_split .OR. .NOT. pbc) THEN
2318 pp_lb(ip) = max(pp_lb(ip), coarse_gbo(1, 1))
2319 pp_ub(ip) = min(pp_ub(ip), coarse_gbo(2, 1))
2320 END IF
2321 END DO
2322
2323 rcv_size = 0
2324 DO ip = 0, n_procs - 1
2325 DO x = pp_lb(ip), coarse_gbo(1, 1) - 1
2326 x_att = coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1))
2327 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
2328 rcv_size(ip) = rcv_size(ip) + coarse_slice_size
2329 END IF
2330 END DO
2331 rcv_size(ip) = rcv_size(ip) + coarse_slice_size* &
2332 max(0, &
2333 min(pp_ub(ip), my_coarse_bo(2, 1)) - max(pp_lb(ip), my_coarse_bo(1, 1)) + 1)
2334 DO x = coarse_gbo(2, 1) + 1, pp_ub(ip)
2335 x_att = coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1))
2336 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
2337 rcv_size(ip) = rcv_size(ip) + coarse_slice_size
2338 END IF
2339 END DO
2340 END DO
2341
2342 ! ** offsets & alloc send-rcv
2343
2344 send_tot_size = 0
2345 DO ip = 0, n_procs - 1
2346 send_offset(ip) = send_tot_size
2347 send_tot_size = send_tot_size + send_size(ip)
2348 END DO
2349 IF (send_tot_size /= (coarse_bo(2, 1) - coarse_bo(1, 1) + 1)*coarse_slice_size) THEN
2350 cpabort("Error calculating send_tot_size")
2351 END IF
2352 ALLOCATE (send_buf(0:send_tot_size - 1))
2353
2354 rcv_tot_size = 0
2355 DO ip = 0, n_procs - 1
2356 rcv_offset(ip) = rcv_tot_size
2357 rcv_tot_size = rcv_tot_size + rcv_size(ip)
2358 END DO
2359 ALLOCATE (rcv_buf(0:rcv_tot_size - 1))
2360
2361 ! ** fill send buffer
2362
2363 pos_of_x => coarse_coeffs_pw%pw_grid%para%pos_of_x
2364 p_old = pos_of_x(coarse_gbo(1, 1) &
2365 + modulo(coarse_bo(1, 1) - coarse_gbo(1, 1), s(1)))
2366 sent_size(:) = send_offset
2367 ss = coarse_bo(2, 1) - coarse_bo(1, 1) + 1
2368 DO x = coarse_bo(1, 1), coarse_bo(2, 1)
2369 p = pos_of_x(coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1)))
2370 CALL dcopy(coarse_slice_size, &
2371 coarse_coeffs(x, coarse_bo(1, 2), &
2372 coarse_bo(1, 3)), ss, send_buf(sent_size(p)), 1)
2373 sent_size(p) = sent_size(p) + coarse_slice_size
2374 END DO
2375
2376 IF (any(sent_size(0:n_procs - 2) /= send_offset(1:n_procs - 1))) THEN
2377 cpabort("error 1 filling send buffer")
2378 END IF
2379 IF (sent_size(n_procs - 1) /= send_tot_size) THEN
2380 cpabort("error 2 filling send buffer")
2381 END IF
2382
2383 IF (local_data) THEN
2384 DEALLOCATE (coarse_coeffs)
2385 ELSE
2386 NULLIFY (coarse_coeffs)
2387 END IF
2388
2389 cpassert(all(sent_size(:n_procs - 2) == send_offset(1:)))
2390 cpassert(sent_size(n_procs - 1) == send_tot_size)
2391 ! test send/rcv sizes
2392 CALL coarse_coeffs_pw%pw_grid%para%group%alltoall(send_size, real_rcv_size, 1)
2393
2394 cpassert(all(real_rcv_size == rcv_size))
2395 ! all2all
2396 CALL coarse_coeffs_pw%pw_grid%para%group%alltoall(sb=send_buf, scount=send_size, sdispl=send_offset, &
2397 rb=rcv_buf, rcount=rcv_size, rdispl=rcv_offset)
2398
2399 ! ** sum & reorder rcv buffer
2400
2401 sent_size(:) = rcv_offset
2402 DO ip = 0, n_procs - 1
2403
2404 DO x = pp_lb(ip), coarse_gbo(1, 1) - 1
2405 x_att = coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1))
2406 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
2407 ii = sent_size(ip)
2408 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
2409 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
2410 coarse_coeffs_pw%array(x_att, j, k) = coarse_coeffs_pw%array(x_att, j, k) + rcv_buf(ii)
2411 ii = ii + 1
2412 END DO
2413 END DO
2414 sent_size(ip) = ii
2415 END IF
2416 END DO
2417
2418 ii = sent_size(ip)
2419 DO x_att = max(pp_lb(ip), my_coarse_bo(1, 1)), min(pp_ub(ip), my_coarse_bo(2, 1))
2420 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
2421 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
2422 coarse_coeffs_pw%array(x_att, j, k) = coarse_coeffs_pw%array(x_att, j, k) + rcv_buf(ii)
2423 ii = ii + 1
2424 END DO
2425 END DO
2426 END DO
2427 sent_size(ip) = ii
2428
2429 DO x = coarse_gbo(2, 1) + 1, pp_ub(ip)
2430 x_att = coarse_gbo(1, 1) + modulo(x - coarse_gbo(1, 1), s(1))
2431 IF (x_att >= my_coarse_bo(1, 1) .AND. x_att <= my_coarse_bo(2, 1)) THEN
2432 ii = sent_size(ip)
2433 DO k = coarse_bo(1, 3), coarse_bo(2, 3)
2434 DO j = coarse_bo(1, 2), coarse_bo(2, 2)
2435 coarse_coeffs_pw%array(x_att, j, k) = coarse_coeffs_pw%array(x_att, j, k) + rcv_buf(ii)
2436 ii = ii + 1
2437 END DO
2438 END DO
2439 sent_size(ip) = ii
2440 END IF
2441 END DO
2442
2443 END DO
2444
2445 IF (any(sent_size(0:n_procs - 2) /= rcv_offset(1:n_procs - 1))) THEN
2446 cpabort("error 1 handling the rcv buffer")
2447 END IF
2448 IF (sent_size(n_procs - 1) /= rcv_tot_size) THEN
2449 cpabort("error 2 handling the rcv buffer")
2450 END IF
2451
2452 ! dealloc
2453 DEALLOCATE (send_size, send_offset, rcv_size, rcv_offset)
2454 DEALLOCATE (send_buf, rcv_buf, real_rcv_size)
2455 DEALLOCATE (pp_ub, pp_lb)
2456 CALL timestop(handle2)
2457 ELSE
2458 cpassert(.NOT. local_data)
2459 END IF
2460
2461 CALL timestop(handle)
2462 END SUBROUTINE add_fine2coarse
2463
2464! **************************************************************************************************
2465!> \brief ...
2466!> \param preconditioner the preconditioner to create
2467!> \param precond_kind the kind of preconditioner to use
2468!> \param pool a pool with grids of the same type as the elements to
2469!> precondition
2470!> \param pbc if periodic boundary conditions should be applied
2471!> \param transpose ...
2472!> \author fawzi
2473! **************************************************************************************************
2474 SUBROUTINE pw_spline_precond_create(preconditioner, precond_kind, &
2475 pool, pbc, transpose)
2476 TYPE(pw_spline_precond_type), INTENT(OUT) :: preconditioner
2477 INTEGER, INTENT(in) :: precond_kind
2478 TYPE(pw_pool_type), INTENT(IN), POINTER :: pool
2479 LOGICAL, INTENT(in) :: pbc, transpose
2480
2481 preconditioner%kind = no_precond
2482 preconditioner%pool => pool
2483 preconditioner%pbc = pbc
2484 preconditioner%transpose = transpose
2485 CALL pool%retain()
2486 CALL pw_spline_precond_set_kind(preconditioner, precond_kind)
2487 END SUBROUTINE pw_spline_precond_create
2488
2489! **************************************************************************************************
2490!> \brief switches the types of preconditioner to use
2491!> \param preconditioner the preconditioner to be changed
2492!> \param precond_kind the new kind of preconditioner to use
2493!> \param pbc ...
2494!> \param transpose ...
2495!> \author fawzi
2496! **************************************************************************************************
2497 SUBROUTINE pw_spline_precond_set_kind(preconditioner, precond_kind, pbc, &
2498 transpose)
2499 TYPE(pw_spline_precond_type), INTENT(INOUT) :: preconditioner
2500 INTEGER, INTENT(in) :: precond_kind
2501 LOGICAL, INTENT(in), OPTIONAL :: pbc, transpose
2502
2503 LOGICAL :: do_3d_coeff
2504 REAL(kind=dp) :: s
2505
2506 IF (PRESENT(transpose)) preconditioner%transpose = transpose
2507 do_3d_coeff = .false.
2508 preconditioner%kind = precond_kind
2509 IF (PRESENT(pbc)) preconditioner%pbc = pbc
2510 SELECT CASE (precond_kind)
2511 CASE (no_precond)
2512 CASE (precond_spl3_aint2)
2513 preconditioner%coeffs_1d = [-1.66_dp*0.25_dp, 1.66_dp, -1.66_dp*0.25_dp]
2514 preconditioner%sharpen = .false.
2515 preconditioner%normalize = .false.
2516 do_3d_coeff = .true.
2517 CASE (precond_spl3_3)
2518 preconditioner%coeffs_1d(1) = -0.25_dp*1.6_dp
2519 preconditioner%coeffs_1d(2) = 1.6_dp
2520 preconditioner%coeffs_1d(3) = -0.25_dp*1.6_dp
2521 preconditioner%sharpen = .false.
2522 preconditioner%normalize = .false.
2523 do_3d_coeff = .true.
2524 CASE (precond_spl3_2)
2525 preconditioner%coeffs_1d(1) = -0.26_dp*1.76_dp
2526 preconditioner%coeffs_1d(2) = 1.76_dp
2527 preconditioner%coeffs_1d(3) = -0.26_dp*1.76_dp
2528 preconditioner%sharpen = .false.
2529 preconditioner%normalize = .false.
2530 do_3d_coeff = .true.
2531 CASE (precond_spl3_aint)
2532 preconditioner%coeffs_1d = spl3_1d_coeffs0
2533 preconditioner%sharpen = .true.
2534 preconditioner%normalize = .true.
2535 do_3d_coeff = .true.
2536 CASE (precond_spl3_1)
2537 preconditioner%coeffs_1d(1) = 0.5_dp/3._dp**(1._dp/3._dp)
2538 preconditioner%coeffs_1d(2) = 4._dp/3._dp**(1._dp/3._dp)
2539 preconditioner%coeffs_1d(3) = 0.5_dp/3._dp**(1._dp/3._dp)
2540 preconditioner%sharpen = .true.
2541 preconditioner%normalize = .false.
2542 do_3d_coeff = .true.
2543 CASE default
2544 cpabort("Unknown preconditioner kind for pw_spline_precond_set_kind")
2545 END SELECT
2546 IF (do_3d_coeff) THEN
2547 s = 1._dp
2548 IF (preconditioner%sharpen) s = -1._dp
2549 preconditioner%coeffs(1) = &
2550 s*preconditioner%coeffs_1d(2)* &
2551 preconditioner%coeffs_1d(2)* &
2552 preconditioner%coeffs_1d(2)
2553 preconditioner%coeffs(2) = &
2554 s*preconditioner%coeffs_1d(1)* &
2555 preconditioner%coeffs_1d(2)* &
2556 preconditioner%coeffs_1d(2)
2557 preconditioner%coeffs(3) = &
2558 s*preconditioner%coeffs_1d(1)* &
2559 preconditioner%coeffs_1d(1)* &
2560 preconditioner%coeffs_1d(2)
2561 preconditioner%coeffs(4) = &
2562 s*preconditioner%coeffs_1d(1)* &
2563 preconditioner%coeffs_1d(1)* &
2564 preconditioner%coeffs_1d(1)
2565 IF (preconditioner%sharpen) THEN
2566 IF (preconditioner%normalize) THEN
2567 preconditioner%coeffs(1) = 2._dp + &
2568 preconditioner%coeffs(1)
2569 ELSE
2570 preconditioner%coeffs(1) = -preconditioner%coeffs(1)
2571 END IF
2572 END IF
2573 END IF
2574 END SUBROUTINE pw_spline_precond_set_kind
2575
2576! **************************************************************************************************
2577!> \brief releases the preconditioner
2578!> \param preconditioner the preconditioner to release
2579!> \author fawzi
2580! **************************************************************************************************
2581 SUBROUTINE pw_spline_precond_release(preconditioner)
2582 TYPE(pw_spline_precond_type), INTENT(INOUT) :: preconditioner
2583
2584 CALL pw_pool_release(preconditioner%pool)
2585 END SUBROUTINE pw_spline_precond_release
2586
2587! **************************************************************************************************
2588!> \brief applies the preconditioner to the system of equations to find the
2589!> coefficients of the spline
2590!> \param preconditioner the preconditioner to apply
2591!> \param in_v the grid on which the preconditioner should be applied
2592!> \param out_v place to store the preconditioner applied on v_out
2593!> \author fawzi
2594! **************************************************************************************************
2595 SUBROUTINE pw_spline_do_precond(preconditioner, in_v, out_v)
2596 TYPE(pw_spline_precond_type), INTENT(IN) :: preconditioner
2597 TYPE(pw_r3d_rs_type), INTENT(IN) :: in_v
2598 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: out_v
2599
2600 SELECT CASE (preconditioner%kind)
2601 CASE (no_precond)
2602 CALL pw_copy(in_v, out_v)
2603 CASE (precond_spl3_aint, precond_spl3_1)
2604 CALL pw_zero(out_v)
2605 IF (preconditioner%pbc) THEN
2606 CALL pw_nn_smear_r(pw_in=in_v, pw_out=out_v, &
2607 coeffs=preconditioner%coeffs)
2608 ELSE
2609 CALL pw_nn_compose_r_no_pbc(weights_1d=preconditioner%coeffs_1d, &
2610 pw_in=in_v, pw_out=out_v, sharpen=preconditioner%sharpen, &
2611 normalize=preconditioner%normalize, &
2612 transpose=preconditioner%transpose)
2613 END IF
2614 CASE (precond_spl3_3, precond_spl3_2, precond_spl3_aint2)
2615 CALL pw_zero(out_v)
2616 IF (preconditioner%pbc) THEN
2617 CALL pw_nn_smear_r(pw_in=in_v, pw_out=out_v, &
2618 coeffs=preconditioner%coeffs)
2619 ELSE
2620 CALL pw_nn_compose_r_no_pbc(weights_1d=preconditioner%coeffs_1d, &
2621 pw_in=in_v, pw_out=out_v, sharpen=preconditioner%sharpen, &
2622 normalize=preconditioner%normalize, smooth_boundary=.true., &
2623 transpose=preconditioner%transpose)
2624 END IF
2625 CASE default
2626 cpabort("Unknown preconditioner kind for pw_spline_do_precond")
2627 END SELECT
2628 END SUBROUTINE pw_spline_do_precond
2629
2630! **************************************************************************************************
2631!> \brief solves iteratively (CG) a systmes of linear equations
2632!> linOp(coeffs)=values
2633!> (for example those needed to find the coefficients of a spline)
2634!> Returns true if the it succeeded to achieve the requested accuracy
2635!> \param values the right hand side of the system
2636!> \param coeffs will contain the solution of the system (and on entry
2637!> it contains the starting point)
2638!> \param linOp the linear operator to be inverted
2639!> \param preconditioner the preconditioner to apply
2640!> \param pool a pool of grids (for the temporary objects)
2641!> \param eps_r the requested precision on the residual
2642!> \param eps_x the requested precision on the solution
2643!> \param max_iter maximum number of iteration allowed
2644!> \param sumtype ...
2645!> \return ...
2646!> \author fawzi
2647! **************************************************************************************************
2648 FUNCTION find_coeffs(values, coeffs, linOp, preconditioner, pool, &
2649 eps_r, eps_x, max_iter, sumtype) RESULT(res)
2650
2651 TYPE(pw_r3d_rs_type), INTENT(IN) :: values
2652 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: coeffs
2653 INTERFACE
2654 SUBROUTINE linop(pw_in, pw_out)
2655 USE pw_types, ONLY: pw_r3d_rs_type
2656 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
2657 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
2658 END SUBROUTINE linop
2659 END INTERFACE
2660 TYPE(pw_spline_precond_type), INTENT(IN) :: preconditioner
2661 TYPE(pw_pool_type), POINTER :: pool
2662 REAL(kind=dp), INTENT(in) :: eps_r, eps_x
2663 INTEGER, INTENT(in) :: max_iter
2664 INTEGER, INTENT(in), OPTIONAL :: sumtype
2665 LOGICAL :: res
2666
2667 INTEGER :: i, iiter, iter, j, k
2668 INTEGER, DIMENSION(2, 3) :: bo
2669 LOGICAL :: last
2670 REAL(kind=dp) :: alpha, beta, eps_r_att, eps_x_att, r_z, &
2671 r_z_new
2672 TYPE(cp_logger_type), POINTER :: logger
2673 TYPE(pw_r3d_rs_type) :: ap, p, r, z
2674
2675 last = .false.
2676
2677 res = .false.
2678 logger => cp_get_default_logger()
2679 CALL pool%create_pw(r)
2680 CALL pool%create_pw(z)
2681 CALL pool%create_pw(p)
2682 CALL pool%create_pw(ap)
2683
2684 !CALL cp_add_iter_level(logger%iter_info,level_name="SPLINE_FIND_COEFFS")
2685 ext_do: DO iiter = 1, max_iter, 10
2686 CALL pw_zero(r)
2687 CALL linop(pw_in=coeffs, pw_out=r)
2688 r%array = -r%array
2689 CALL pw_axpy(values, r)
2690 CALL pw_spline_do_precond(preconditioner, in_v=r, out_v=z)
2691 CALL pw_copy(z, p)
2692 r_z = pw_integral_ab(r, z, sumtype)
2693
2694 DO iter = iiter, min(iiter + 9, max_iter)
2695 eps_r_att = sqrt(pw_integral_ab(r, r, sumtype))
2696 IF (eps_r_att == 0._dp) THEN
2697 eps_x_att = 0._dp
2698 last = .true.
2699 ELSE
2700 CALL pw_zero(ap)
2701 CALL linop(pw_in=p, pw_out=ap)
2702 alpha = r_z/pw_integral_ab(ap, p, sumtype)
2703
2704 CALL pw_axpy(p, coeffs, alpha=alpha)
2705
2706 eps_x_att = alpha*sqrt(pw_integral_ab(p, p, sumtype)) ! try to spare if unneeded?
2707 IF (eps_r_att < eps_r .AND. eps_x_att < eps_x) last = .true.
2708 END IF
2709 !CALL cp_iterate(logger%iter_info,last=last)
2710 IF (last) THEN
2711 res = .true.
2712 EXIT ext_do
2713 END IF
2714
2715 CALL pw_axpy(ap, r, alpha=-alpha)
2716
2717 CALL pw_spline_do_precond(preconditioner, in_v=r, out_v=z)
2718
2719 r_z_new = pw_integral_ab(r, z, sumtype)
2720 beta = r_z_new/r_z
2721 r_z = r_z_new
2722
2723 bo = p%pw_grid%bounds_local
2724 DO k = bo(1, 3), bo(2, 3)
2725 DO j = bo(1, 2), bo(2, 2)
2726 DO i = bo(1, 1), bo(2, 1)
2727 p%array(i, j, k) = z%array(i, j, k) + beta*p%array(i, j, k)
2728 END DO
2729 END DO
2730 END DO
2731
2732 END DO
2733 END DO ext_do
2734 !CALL cp_rm_iter_level(logger%iter_info,level_name="SPLINE_FIND_COEFFS")
2735
2736 CALL pool%give_back_pw(r)
2737 CALL pool%give_back_pw(z)
2738 CALL pool%give_back_pw(p)
2739 CALL pool%give_back_pw(ap)
2740
2741 END FUNCTION find_coeffs
2742
2743! **************************************************************************************************
2744!> \brief adds to pw_out pw_in composed with the weights
2745!> pw_out%array(i,j,k)=pw_out%array(i,j,k)+sum(pw_in%array(i+l,j+m,k+n)*
2746!> weights_1d(abs(l)+1)*weights_1d(abs(m)+1)*weights_1d(abs(n)+1),
2747!> l=-1..1,m=-1..1,n=-1..1)
2748!> \param weights_1d ...
2749!> \param pw_in ...
2750!> \param pw_out ...
2751!> \param sharpen ...
2752!> \param normalize ...
2753!> \param transpose ...
2754!> \param smooth_boundary ...
2755!> \author fawzi
2756! **************************************************************************************************
2757 SUBROUTINE pw_nn_compose_r_no_pbc(weights_1d, pw_in, pw_out, &
2758 sharpen, normalize, transpose, smooth_boundary)
2759 REAL(kind=dp), DIMENSION(-1:1) :: weights_1d
2760 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in, pw_out
2761 LOGICAL, INTENT(in), OPTIONAL :: sharpen, normalize, transpose, &
2762 smooth_boundary
2763
2764 INTEGER :: first_index, i, j, jw, k, kw, &
2765 last_index, myj, myk, n_els
2766 INTEGER, DIMENSION(2, 3) :: bo, gbo
2767 INTEGER, DIMENSION(3) :: s
2768 LOGICAL :: has_l_boundary, has_u_boundary, &
2769 is_split, my_normalize, my_sharpen, &
2770 my_smooth_boundary, my_transpose
2771 REAL(kind=dp) :: in_val_f, in_val_l, in_val_tmp, w_j, w_k
2772 REAL(kind=dp), DIMENSION(-1:1) :: w
2773 REAL(kind=dp), DIMENSION(:, :), POINTER :: l_boundary, tmp, u_boundary
2774 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: in_val, out_val
2775
2776 bo = pw_in%pw_grid%bounds_local
2777 gbo = pw_in%pw_grid%bounds
2778 in_val => pw_in%array
2779 out_val => pw_out%array
2780 my_sharpen = .false.
2781 IF (PRESENT(sharpen)) my_sharpen = sharpen
2782 my_normalize = .false.
2783 IF (PRESENT(normalize)) my_normalize = normalize
2784 my_transpose = .false.
2785 IF (PRESENT(transpose)) my_transpose = transpose
2786 my_smooth_boundary = .false.
2787 IF (PRESENT(smooth_boundary)) my_smooth_boundary = smooth_boundary
2788 cpassert(.NOT. my_normalize .OR. my_sharpen)
2789 cpassert(.NOT. my_smooth_boundary .OR. .NOT. my_sharpen)
2790 DO i = 1, 3
2791 s(i) = bo(2, i) - bo(1, i) + 1
2792 END DO
2793 IF (any(s < 1)) RETURN
2794 is_split = any(pw_in%pw_grid%bounds_local(:, 1) /= &
2795 pw_in%pw_grid%bounds(:, 1))
2796 has_l_boundary = (gbo(1, 1) == bo(1, 1))
2797 has_u_boundary = (gbo(2, 1) == bo(2, 1))
2798 IF (is_split) THEN
2799 ALLOCATE (l_boundary(bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)), &
2800 u_boundary(bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)), &
2801 tmp(bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)))
2802 tmp(:, :) = pw_in%array(bo(2, 1), :, :)
2803 CALL pw_in%pw_grid%para%group%sendrecv(tmp, pw_in%pw_grid%para%pos_of_x( &
2804 gbo(1, 1) + modulo(bo(2, 1) + 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)), &
2805 l_boundary, 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 tmp(:, :) = pw_in%array(bo(1, 1), :, :)
2808 CALL pw_in%pw_grid%para%group%sendrecv(tmp, pw_in%pw_grid%para%pos_of_x( &
2809 gbo(1, 1) + modulo(bo(1, 1) - 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)), &
2810 u_boundary, pw_in%pw_grid%para%pos_of_x( &
2811 gbo(1, 1) + modulo(bo(2, 1) + 1 - gbo(1, 1), gbo(2, 1) - gbo(1, 1) + 1)))
2812 DEALLOCATE (tmp)
2813 END IF
2814
2815 n_els = s(1)
2816 IF (has_l_boundary) THEN
2817 n_els = n_els - 1
2818 first_index = bo(1, 1) + 1
2819 ELSE
2820 first_index = bo(1, 1)
2821 END IF
2822 IF (has_u_boundary) THEN
2823 n_els = n_els - 1
2824 last_index = bo(2, 1) - 1
2825 ELSE
2826 last_index = bo(2, 1)
2827 END IF
2828!$OMP PARALLEL DO DEFAULT(NONE) &
2829!$OMP PRIVATE(k, kw, myk, j, jw, myj, in_val_f, in_val_l, w_k, w_j, in_val_tmp, w) &
2830!$OMP SHARED(bo, in_val, out_val, s, l_boundary, u_boundary, weights_1d, is_split, &
2831!$OMP my_transpose, gbo, my_smooth_boundary, has_l_boundary, has_u_boundary, &
2832!$OMP my_sharpen, last_index, first_index, my_normalize, n_els)
2833 DO k = bo(1, 3), bo(2, 3)
2834 DO kw = -1, 1
2835 myk = k + kw
2836 IF (my_transpose) THEN
2837 IF (k >= gbo(2, 3) - 1 .OR. k <= gbo(1, 3) + 1) THEN
2838 IF (k == gbo(2, 3) .OR. k == gbo(1, 3)) THEN
2839 IF (myk < gbo(2, 3) .AND. myk > gbo(1, 3)) THEN
2840 w_k = weights_1d(kw)
2841 IF (my_smooth_boundary) THEN
2842 w_k = weights_1d(kw)/weights_1d(0)
2843 END IF
2844 ELSE IF (kw == 0) THEN
2845 w_k = 1._dp
2846 ELSE
2847 cycle
2848 END IF
2849 ELSE
2850 IF (myk == gbo(2, 3) .OR. myk == gbo(1, 3)) cycle
2851 w_k = weights_1d(kw)
2852 END IF
2853 ELSE
2854 w_k = weights_1d(kw)
2855 END IF
2856 ELSE
2857 IF (k >= gbo(2, 3) - 1 .OR. k <= gbo(1, 3) + 1) THEN
2858 IF (k == gbo(2, 3) .OR. k == gbo(1, 3)) THEN
2859 IF (kw /= 0) cycle
2860 w_k = 1._dp
2861 ELSE
2862 IF (my_smooth_boundary .AND. ((k == gbo(1, 3) + 1 .AND. myk == gbo(1, 3)) .OR. &
2863 (k == gbo(2, 3) - 1 .AND. myk == gbo(2, 3)))) THEN
2864 w_k = weights_1d(kw)/weights_1d(0)
2865 ELSE
2866 w_k = weights_1d(kw)
2867 END IF
2868 END IF
2869 ELSE
2870 w_k = weights_1d(kw)
2871 END IF
2872 END IF
2873 DO j = bo(1, 2), bo(2, 2)
2874 DO jw = -1, 1
2875 myj = j + jw
2876 IF (j < gbo(2, 2) - 1 .AND. j > gbo(1, 2) + 1) THEN
2877 w_j = w_k*weights_1d(jw)
2878 ELSE
2879 IF (my_transpose) THEN
2880 IF (j == gbo(2, 2) .OR. j == gbo(1, 2)) THEN
2881 IF (myj < gbo(2, 2) .AND. myj > gbo(1, 2)) THEN
2882 w_j = weights_1d(jw)*w_k
2883 IF (my_smooth_boundary) THEN
2884 w_j = weights_1d(jw)/weights_1d(0)*w_k
2885 END IF
2886 ELSE IF (jw == 0) THEN
2887 w_j = w_k
2888 ELSE
2889 cycle
2890 END IF
2891 ELSE
2892 IF (myj == gbo(2, 2) .OR. myj == gbo(1, 2)) cycle
2893 w_j = w_k*weights_1d(jw)
2894 END IF
2895 ELSE
2896 IF (j == gbo(2, 2) .OR. j == gbo(1, 2)) THEN
2897 IF (jw /= 0) cycle
2898 w_j = w_k
2899 ELSE IF (my_smooth_boundary .AND. ((j == gbo(1, 2) + 1 .AND. myj == gbo(1, 2)) .OR. &
2900 (j == gbo(2, 2) - 1 .AND. myj == gbo(2, 2)))) THEN
2901 w_j = w_k*weights_1d(jw)/weights_1d(0)
2902 ELSE
2903 w_j = w_k*weights_1d(jw)
2904 END IF
2905 END IF
2906 END IF
2907
2908 IF (has_l_boundary) THEN
2909 IF (my_transpose) THEN
2910 IF (s(1) == 1) THEN
2911 cpassert(.NOT. has_u_boundary)
2912 in_val_tmp = u_boundary(myj, myk)
2913 ELSE
2914 in_val_tmp = in_val(bo(1, 1) + 1, myj, myk)
2915 END IF
2916 IF (my_sharpen) THEN
2917 IF (kw == 0 .AND. jw == 0) THEN
2918 IF (my_normalize) THEN
2919 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2920 (2.0_dp - w_j)*in_val(bo(1, 1), myj, myk) - &
2921 in_val_tmp*weights_1d(1)*w_j
2922 ELSE
2923 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2924 in_val(bo(1, 1), myj, myk)*w_j - &
2925 in_val_tmp*weights_1d(1)*w_j
2926 END IF
2927 ELSE
2928 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) - &
2929 in_val(bo(1, 1), myj, myk)*w_j - &
2930 in_val_tmp*weights_1d(1)*w_j
2931 END IF
2932 ELSE IF (my_smooth_boundary) THEN
2933 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2934 w_j*(in_val(bo(1, 1), myj, myk) + &
2935 in_val_tmp*weights_1d(1)/weights_1d(0))
2936 ELSE
2937 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2938 w_j*(in_val(bo(1, 1), myj, myk) + &
2939 in_val_tmp*weights_1d(1))
2940 END IF
2941 in_val_f = 0.0_dp
2942 ELSE
2943 in_val_f = in_val(bo(1, 1), myj, myk)
2944 IF (my_sharpen) THEN
2945 IF (kw == 0 .AND. jw == 0) THEN
2946 IF (my_normalize) THEN
2947 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2948 (2.0_dp - w_j)*in_val_f
2949 ELSE
2950 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2951 in_val_f*w_j
2952 END IF
2953 ELSE
2954 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) - &
2955 in_val_f*w_j
2956 END IF
2957 ELSE
2958 out_val(bo(1, 1), j, k) = out_val(bo(1, 1), j, k) + &
2959 in_val_f*w_j
2960 END IF
2961 END IF
2962 ELSE
2963 in_val_f = l_boundary(myj, myk)
2964 END IF
2965 IF (has_u_boundary) THEN
2966 IF (my_transpose) THEN
2967 in_val_l = in_val(bo(2, 1), myj, myk)
2968 IF (s(1) == 1) THEN
2969 cpassert(.NOT. has_l_boundary)
2970 in_val_tmp = l_boundary(myj, myk)
2971 ELSE
2972 in_val_tmp = in_val(bo(2, 1) - 1, myj, myk)
2973 END IF
2974 IF (my_sharpen) THEN
2975 IF (kw == 0 .AND. jw == 0) THEN
2976 IF (my_normalize) THEN
2977 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
2978 in_val_l*(2._dp - w_j) - &
2979 in_val_tmp*weights_1d(1)*w_j
2980 ELSE
2981 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
2982 in_val_l*w_j - &
2983 in_val_tmp*weights_1d(1)*w_j
2984 END IF
2985 ELSE
2986 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) - &
2987 w_j*in_val_l - &
2988 in_val_tmp*weights_1d(1)*w_j
2989 END IF
2990 ELSE IF (my_smooth_boundary) THEN
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)/weights_1d(0))
2993 ELSE
2994 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
2995 w_j*(in_val_l + in_val_tmp*weights_1d(1))
2996 END IF
2997 in_val_l = 0._dp
2998 ELSE
2999 in_val_l = in_val(bo(2, 1), myj, myk)
3000 IF (my_sharpen) THEN
3001 IF (kw == 0 .AND. jw == 0) THEN
3002 IF (my_normalize) THEN
3003 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
3004 in_val_l*(2._dp - w_j)
3005 ELSE
3006 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
3007 in_val_l*w_j
3008 END IF
3009 ELSE
3010 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) - &
3011 w_j*in_val_l
3012 END IF
3013 ELSE
3014 out_val(bo(2, 1), j, k) = out_val(bo(2, 1), j, k) + &
3015 w_j*in_val_l
3016 END IF
3017 END IF
3018 ELSE
3019 in_val_l = u_boundary(myj, myk)
3020 END IF
3021 IF (last_index >= first_index) THEN
3022 IF (my_transpose) THEN
3023 IF (bo(1, 1) - 1 == gbo(1, 1)) THEN
3024 in_val_f = 0._dp
3025 ELSE IF (bo(2, 1) + 1 == gbo(2, 1)) THEN
3026 in_val_l = 0._dp
3027 END IF
3028 END IF
3029 IF (my_sharpen) THEN
3030 w = -weights_1d*w_j
3031 IF (kw == 0 .AND. jw == 0) THEN
3032 IF (my_normalize) THEN
3033 w(0) = w(0) + 2._dp
3034 ELSE
3035 w(0) = -w(0)
3036 END IF
3037 END IF
3038 ELSE
3039 w = weights_1d*w_j
3040 END IF
3041 IF (my_smooth_boundary .AND. (.NOT. my_transpose)) THEN
3042 IF (gbo(1, 1) + 1 >= bo(1, 1) .AND. &
3043 gbo(1, 1) + 1 <= bo(2, 1) .AND. gbo(2, 1) - gbo(1, 1) > 2) THEN
3044 IF (gbo(1, 1) >= bo(1, 1)) THEN
3045 out_val(gbo(1, 1) + 1, j, k) = out_val(gbo(1, 1) + 1, j, k) + &
3046 in_val(gbo(1, 1), myj, myk)*w_j*weights_1d(-1)* &
3047 (1._dp/weights_1d(0) - 1._dp)
3048 ELSE
3049 out_val(gbo(1, 1) + 1, j, k) = out_val(gbo(1, 1) + 1, j, k) + &
3050 l_boundary(myj, myk)*w_j*weights_1d(-1)* &
3051 (1._dp/weights_1d(0) - 1._dp)
3052 END IF
3053 END IF
3054 END IF
3055 CALL pw_compose_stripe(weights=w, &
3056 in_val=in_val(first_index:last_index, myj, myk), &
3057 in_val_first=in_val_f, in_val_last=in_val_l, &
3058 out_val=out_val(first_index:last_index, j, k), &
3059 n_el=n_els)
3060!FM call pw_compose_stripe2(weights=w,&
3061!FM in_val=in_val,&
3062!FM in_val_first=in_val_f,in_val_last=in_val_l,&
3063!FM out_val=out_val,&
3064!FM first_val=first_index,last_val=last_index,&
3065!FM myj=myj,myk=myk,j=j,k=k)
3066 IF (my_smooth_boundary .AND. (.NOT. my_transpose)) THEN
3067 IF (gbo(2, 1) - 1 >= bo(1, 1) .AND. &
3068 gbo(2, 1) - 1 <= bo(2, 1) .AND. gbo(2, 1) - gbo(1, 1) > 2) THEN
3069 IF (gbo(2, 1) <= bo(2, 1)) THEN
3070 out_val(gbo(2, 1) - 1, j, k) = out_val(gbo(2, 1) - 1, j, k) + &
3071 in_val(gbo(2, 1), myj, myk)*w_j*weights_1d(1)* &
3072 (1._dp/weights_1d(0) - 1._dp)
3073 ELSE
3074 out_val(gbo(2, 1) - 1, j, k) = out_val(gbo(2, 1) - 1, j, k) + &
3075 u_boundary(myj, myk)*w_j*weights_1d(1)* &
3076 (1._dp/weights_1d(0) - 1._dp)
3077 END IF
3078 END IF
3079 END IF
3080
3081 END IF
3082 END DO
3083 END DO
3084 END DO
3085 END DO
3086
3087 IF (is_split) THEN
3088 DEALLOCATE (l_boundary, u_boundary)
3089 END IF
3090 END SUBROUTINE pw_nn_compose_r_no_pbc
3091
3092! **************************************************************************************************
3093!> \brief ...
3094!> \param pw_in ...
3095!> \param pw_out ...
3096! **************************************************************************************************
3097 SUBROUTINE spl3_nopbc(pw_in, pw_out)
3098
3099 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
3100 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
3101
3102 CALL pw_zero(pw_out)
3103 CALL pw_nn_compose_r_no_pbc(weights_1d=spl3_1d_coeffs0, pw_in=pw_in, &
3104 pw_out=pw_out, sharpen=.false., normalize=.false.)
3105
3106 END SUBROUTINE spl3_nopbc
3107
3108! **************************************************************************************************
3109!> \brief ...
3110!> \param pw_in ...
3111!> \param pw_out ...
3112! **************************************************************************************************
3113 SUBROUTINE spl3_nopbct(pw_in, pw_out)
3114
3115 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
3116 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
3117
3118 CALL pw_zero(pw_out)
3119 CALL pw_nn_compose_r_no_pbc(weights_1d=spl3_1d_coeffs0, pw_in=pw_in, &
3120 pw_out=pw_out, sharpen=.false., normalize=.false., transpose=.true.)
3121
3122 END SUBROUTINE spl3_nopbct
3123
3124! **************************************************************************************************
3125!> \brief ...
3126!> \param pw_in ...
3127!> \param pw_out ...
3128! **************************************************************************************************
3129 SUBROUTINE spl3_pbc(pw_in, pw_out)
3130
3131 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw_in
3132 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: pw_out
3133
3134 CALL pw_zero(pw_out)
3135 CALL pw_nn_smear_r(pw_in, pw_out, coeffs=spline3_coeffs)
3136
3137 END SUBROUTINE spl3_pbc
3138
3139! **************************************************************************************************
3140!> \brief Evaluates the PBC interpolated Spline (pw) function on the generic
3141!> input vector (vec)
3142!> \param vec ...
3143!> \param pw ...
3144!> \return ...
3145!> \par History
3146!> 12.2007 Adapted for use with distributed grids [rdeclerck]
3147!> \author Teodoro Laino 12/2005 [tlaino]
3148!> \note
3149!> Requires the Spline coefficients to be computed with PBC
3150! **************************************************************************************************
3151 FUNCTION eval_interp_spl3_pbc(vec, pw) RESULT(val)
3152 REAL(kind=dp), DIMENSION(3), INTENT(in) :: vec
3153 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
3154 REAL(kind=dp) :: val
3155
3156 INTEGER :: i, ivec(3), j, k, npts(3)
3157 INTEGER, DIMENSION(2, 3) :: bo, bo_l
3158 INTEGER, DIMENSION(4) :: ii, ij, ik
3159 LOGICAL :: my_mpsum
3160 REAL(kind=dp) :: a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3, dr1, dr2, dr3, e1, e2, e3, &
3161 f1, f2, f3, g1, g2, g3, h1, h2, h3, p1, p2, p3, q1, q2, q3, r1, r2, r3, s1, s2, s3, s4, &
3162 t1, t2, t3, t4, u1, u2, u3, v1, v2, v3, v4, xd1, xd2, xd3
3163 REAL(kind=dp), DIMENSION(4, 4, 4) :: box
3164 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid
3165
3166 NULLIFY (grid)
3167 my_mpsum = (pw%pw_grid%para%mode /= pw_mode_local)
3168 npts = pw%pw_grid%npts
3169 ivec = floor(vec/pw%pw_grid%dr)
3170 dr1 = pw%pw_grid%dr(1)
3171 dr2 = pw%pw_grid%dr(2)
3172 dr3 = pw%pw_grid%dr(3)
3173
3174 xd1 = (vec(1)/dr1) - real(ivec(1), kind=dp)
3175 xd2 = (vec(2)/dr2) - real(ivec(2), kind=dp)
3176 xd3 = (vec(3)/dr3) - real(ivec(3), kind=dp)
3177 grid => pw%array(:, :, :)
3178 bo = pw%pw_grid%bounds
3179 bo_l = pw%pw_grid%bounds_local
3180
3181 ik(1) = modulo(ivec(3) - 1, npts(3)) + bo(1, 3)
3182 ik(2) = modulo(ivec(3), npts(3)) + bo(1, 3)
3183 ik(3) = modulo(ivec(3) + 1, npts(3)) + bo(1, 3)
3184 ik(4) = modulo(ivec(3) + 2, npts(3)) + bo(1, 3)
3185
3186 ij(1) = modulo(ivec(2) - 1, npts(2)) + bo(1, 2)
3187 ij(2) = modulo(ivec(2), npts(2)) + bo(1, 2)
3188 ij(3) = modulo(ivec(2) + 1, npts(2)) + bo(1, 2)
3189 ij(4) = modulo(ivec(2) + 2, npts(2)) + bo(1, 2)
3190
3191 ii(1) = modulo(ivec(1) - 1, npts(1)) + bo(1, 1)
3192 ii(2) = modulo(ivec(1), npts(1)) + bo(1, 1)
3193 ii(3) = modulo(ivec(1) + 1, npts(1)) + bo(1, 1)
3194 ii(4) = modulo(ivec(1) + 2, npts(1)) + bo(1, 1)
3195
3196 DO k = 1, 4
3197 DO j = 1, 4
3198 DO i = 1, 4
3199 IF ( &
3200 ii(i) >= bo_l(1, 1) .AND. &
3201 ii(i) <= bo_l(2, 1) .AND. &
3202 ij(j) >= bo_l(1, 2) .AND. &
3203 ij(j) <= bo_l(2, 2) .AND. &
3204 ik(k) >= bo_l(1, 3) .AND. &
3205 ik(k) <= bo_l(2, 3) &
3206 ) THEN
3207 box(i, j, k) = grid(ii(i) + 1 - bo_l(1, 1), &
3208 ij(j) + 1 - bo_l(1, 2), &
3209 ik(k) + 1 - bo_l(1, 3))
3210 ELSE
3211 box(i, j, k) = 0.0_dp
3212 END IF
3213 END DO
3214 END DO
3215 END DO
3216
3217 a1 = 3.0_dp + xd1
3218 a2 = a1*a1
3219 a3 = a2*a1
3220 b1 = 2.0_dp + xd1
3221 b2 = b1*b1
3222 b3 = b2*b1
3223 c1 = 1.0_dp + xd1
3224 c2 = c1*c1
3225 c3 = c2*c1
3226 d1 = xd1
3227 d2 = d1*d1
3228 d3 = d2*d1
3229 e1 = 3.0_dp + xd2
3230 e2 = e1*e1
3231 e3 = e2*e1
3232 f1 = 2.0_dp + xd2
3233 f2 = f1*f1
3234 f3 = f2*f1
3235 g1 = 1.0_dp + xd2
3236 g2 = g1*g1
3237 g3 = g2*g1
3238 h1 = xd2
3239 h2 = h1*h1
3240 h3 = h2*h1
3241 p1 = 3.0_dp + xd3
3242 p2 = p1*p1
3243 p3 = p2*p1
3244 q1 = 2.0_dp + xd3
3245 q2 = q1*q1
3246 q3 = q2*q1
3247 r1 = 1.0_dp + xd3
3248 r2 = r1*r1
3249 r3 = r2*r1
3250 u1 = xd3
3251 u2 = u1*u1
3252 u3 = u2*u1
3253
3254 t1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3)
3255 t2 = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3
3256 t3 = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3
3257 t4 = 1.0_dp/6.0_dp*d3
3258 s1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3)
3259 s2 = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3
3260 s3 = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3
3261 s4 = 1.0_dp/6.0_dp*h3
3262 v1 = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3)
3263 v2 = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3
3264 v3 = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3
3265 v4 = 1.0_dp/6.0_dp*u3
3266
3267 val = ((box(1, 1, 1)*t1 + box(2, 1, 1)*t2 + box(3, 1, 1)*t3 + box(4, 1, 1)*t4)*s1 + &
3268 (box(1, 2, 1)*t1 + box(2, 2, 1)*t2 + box(3, 2, 1)*t3 + box(4, 2, 1)*t4)*s2 + &
3269 (box(1, 3, 1)*t1 + box(2, 3, 1)*t2 + box(3, 3, 1)*t3 + box(4, 3, 1)*t4)*s3 + &
3270 (box(1, 4, 1)*t1 + box(2, 4, 1)*t2 + box(3, 4, 1)*t3 + box(4, 4, 1)*t4)*s4)*v1 + &
3271 ((box(1, 1, 2)*t1 + box(2, 1, 2)*t2 + box(3, 1, 2)*t3 + box(4, 1, 2)*t4)*s1 + &
3272 (box(1, 2, 2)*t1 + box(2, 2, 2)*t2 + box(3, 2, 2)*t3 + box(4, 2, 2)*t4)*s2 + &
3273 (box(1, 3, 2)*t1 + box(2, 3, 2)*t2 + box(3, 3, 2)*t3 + box(4, 3, 2)*t4)*s3 + &
3274 (box(1, 4, 2)*t1 + box(2, 4, 2)*t2 + box(3, 4, 2)*t3 + box(4, 4, 2)*t4)*s4)*v2 + &
3275 ((box(1, 1, 3)*t1 + box(2, 1, 3)*t2 + box(3, 1, 3)*t3 + box(4, 1, 3)*t4)*s1 + &
3276 (box(1, 2, 3)*t1 + box(2, 2, 3)*t2 + box(3, 2, 3)*t3 + box(4, 2, 3)*t4)*s2 + &
3277 (box(1, 3, 3)*t1 + box(2, 3, 3)*t2 + box(3, 3, 3)*t3 + box(4, 3, 3)*t4)*s3 + &
3278 (box(1, 4, 3)*t1 + box(2, 4, 3)*t2 + box(3, 4, 3)*t3 + box(4, 4, 3)*t4)*s4)*v3 + &
3279 ((box(1, 1, 4)*t1 + box(2, 1, 4)*t2 + box(3, 1, 4)*t3 + box(4, 1, 4)*t4)*s1 + &
3280 (box(1, 2, 4)*t1 + box(2, 2, 4)*t2 + box(3, 2, 4)*t3 + box(4, 2, 4)*t4)*s2 + &
3281 (box(1, 3, 4)*t1 + box(2, 3, 4)*t2 + box(3, 3, 4)*t3 + box(4, 3, 4)*t4)*s3 + &
3282 (box(1, 4, 4)*t1 + box(2, 4, 4)*t2 + box(3, 4, 4)*t3 + box(4, 4, 4)*t4)*s4)*v4
3283
3284 IF (my_mpsum) CALL pw%pw_grid%para%group%sum(val)
3285
3286 END FUNCTION eval_interp_spl3_pbc
3287
3288! **************************************************************************************************
3289!> \brief Evaluates the derivatives of the PBC interpolated Spline (pw)
3290!> function on the generic input vector (vec)
3291!> \param vec ...
3292!> \param pw ...
3293!> \return ...
3294!> \par History
3295!> 12.2007 Adapted for use with distributed grids [rdeclerck]
3296!> \author Teodoro Laino 12/2005 [tlaino]
3297!> \note
3298!> Requires the Spline coefficients to be computed with PBC
3299! **************************************************************************************************
3300 FUNCTION eval_d_interp_spl3_pbc(vec, pw) RESULT(val)
3301 REAL(kind=dp), DIMENSION(3), INTENT(in) :: vec
3302 TYPE(pw_r3d_rs_type), INTENT(IN) :: pw
3303 REAL(kind=dp) :: val(3)
3304
3305 INTEGER :: i, ivec(3), j, k, npts(3)
3306 INTEGER, DIMENSION(2, 3) :: bo, bo_l
3307 INTEGER, DIMENSION(4) :: ii, ij, ik
3308 LOGICAL :: my_mpsum
3309 REAL(kind=dp) :: a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3, dr1, dr1i, dr2, dr2i, dr3, &
3310 dr3i, e1, e2, e3, f1, f2, f3, g1, g2, g3, h1, h2, h3, p1, p2, p3, q1, q2, q3, r1, r2, r3, &
3311 s1, s1d, s1o, s2, s2d, s2o, s3, s3d, s3o, s4, s4d, s4o, t1, t1d, t1o, t2, t2d, t2o, t3, &
3312 t3d, t3o, t4, t4d, t4o, u1, u2, u3, v1, v1d, v1o, v2, v2d, v2o, v3, v3d, v3o, v4, v4d, &
3313 v4o, xd1, xd2, xd3
3314 REAL(kind=dp), DIMENSION(4, 4, 4) :: box
3315 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid
3316
3317 NULLIFY (grid)
3318 my_mpsum = (pw%pw_grid%para%mode /= pw_mode_local)
3319 npts = pw%pw_grid%npts
3320 ivec = floor(vec/pw%pw_grid%dr)
3321 dr1 = pw%pw_grid%dr(1)
3322 dr2 = pw%pw_grid%dr(2)
3323 dr3 = pw%pw_grid%dr(3)
3324 dr1i = 1.0_dp/dr1
3325 dr2i = 1.0_dp/dr2
3326 dr3i = 1.0_dp/dr3
3327 xd1 = (vec(1)/dr1) - real(ivec(1), kind=dp)
3328 xd2 = (vec(2)/dr2) - real(ivec(2), kind=dp)
3329 xd3 = (vec(3)/dr3) - real(ivec(3), kind=dp)
3330 grid => pw%array(:, :, :)
3331 bo = pw%pw_grid%bounds
3332 bo_l = pw%pw_grid%bounds_local
3333
3334 ik(1) = modulo(ivec(3) - 1, npts(3)) + bo(1, 3)
3335 ik(2) = modulo(ivec(3), npts(3)) + bo(1, 3)
3336 ik(3) = modulo(ivec(3) + 1, npts(3)) + bo(1, 3)
3337 ik(4) = modulo(ivec(3) + 2, npts(3)) + bo(1, 3)
3338
3339 ij(1) = modulo(ivec(2) - 1, npts(2)) + bo(1, 2)
3340 ij(2) = modulo(ivec(2), npts(2)) + bo(1, 2)
3341 ij(3) = modulo(ivec(2) + 1, npts(2)) + bo(1, 2)
3342 ij(4) = modulo(ivec(2) + 2, npts(2)) + bo(1, 2)
3343
3344 ii(1) = modulo(ivec(1) - 1, npts(1)) + bo(1, 1)
3345 ii(2) = modulo(ivec(1), npts(1)) + bo(1, 1)
3346 ii(3) = modulo(ivec(1) + 1, npts(1)) + bo(1, 1)
3347 ii(4) = modulo(ivec(1) + 2, npts(1)) + bo(1, 1)
3348
3349 DO k = 1, 4
3350 DO j = 1, 4
3351 DO i = 1, 4
3352 IF ( &
3353 ii(i) >= bo_l(1, 1) .AND. &
3354 ii(i) <= bo_l(2, 1) .AND. &
3355 ij(j) >= bo_l(1, 2) .AND. &
3356 ij(j) <= bo_l(2, 2) .AND. &
3357 ik(k) >= bo_l(1, 3) .AND. &
3358 ik(k) <= bo_l(2, 3) &
3359 ) THEN
3360 box(i, j, k) = grid(ii(i) + 1 - bo_l(1, 1), &
3361 ij(j) + 1 - bo_l(1, 2), &
3362 ik(k) + 1 - bo_l(1, 3))
3363 ELSE
3364 box(i, j, k) = 0.0_dp
3365 END IF
3366 END DO
3367 END DO
3368 END DO
3369
3370 a1 = 3.0_dp + xd1
3371 a2 = a1*a1
3372 a3 = a2*a1
3373 b1 = 2.0_dp + xd1
3374 b2 = b1*b1
3375 b3 = b2*b1
3376 c1 = 1.0_dp + xd1
3377 c2 = c1*c1
3378 c3 = c2*c1
3379 d1 = xd1
3380 d2 = d1*d1
3381 d3 = d2*d1
3382 e1 = 3.0_dp + xd2
3383 e2 = e1*e1
3384 e3 = e2*e1
3385 f1 = 2.0_dp + xd2
3386 f2 = f1*f1
3387 f3 = f2*f1
3388 g1 = 1.0_dp + xd2
3389 g2 = g1*g1
3390 g3 = g2*g1
3391 h1 = xd2
3392 h2 = h1*h1
3393 h3 = h2*h1
3394 p1 = 3.0_dp + xd3
3395 p2 = p1*p1
3396 p3 = p2*p1
3397 q1 = 2.0_dp + xd3
3398 q2 = q1*q1
3399 q3 = q2*q1
3400 r1 = 1.0_dp + xd3
3401 r2 = r1*r1
3402 r3 = r2*r1
3403 u1 = xd3
3404 u2 = u1*u1
3405 u3 = u2*u1
3406
3407 t1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*a1 + 12.0_dp*a2 - a3)
3408 t2o = -22.0_dp/3.0_dp + 10.0_dp*b1 - 4.0_dp*b2 + 0.5_dp*b3
3409 t3o = 2.0_dp/3.0_dp - 2.0_dp*c1 + 2.0_dp*c2 - 0.5_dp*c3
3410 t4o = 1.0_dp/6.0_dp*d3
3411 s1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*e1 + 12.0_dp*e2 - e3)
3412 s2o = -22.0_dp/3.0_dp + 10.0_dp*f1 - 4.0_dp*f2 + 0.5_dp*f3
3413 s3o = 2.0_dp/3.0_dp - 2.0_dp*g1 + 2.0_dp*g2 - 0.5_dp*g3
3414 s4o = 1.0_dp/6.0_dp*h3
3415 v1o = 1.0_dp/6.0_dp*(64.0_dp - 48.0_dp*p1 + 12.0_dp*p2 - p3)
3416 v2o = -22.0_dp/3.0_dp + 10.0_dp*q1 - 4.0_dp*q2 + 0.5_dp*q3
3417 v3o = 2.0_dp/3.0_dp - 2.0_dp*r1 + 2.0_dp*r2 - 0.5_dp*r3
3418 v4o = 1.0_dp/6.0_dp*u3
3419
3420 t1d = -8.0_dp + 4.0_dp*a1 - 0.5_dp*a2
3421 t2d = 10.0_dp - 8.0_dp*b1 + 1.5_dp*b2
3422 t3d = -2.0_dp + 4.0_dp*c1 - 1.5_dp*c2
3423 t4d = 0.5_dp*d2
3424 s1d = -8.0_dp + 4.0_dp*e1 - 0.5_dp*e2
3425 s2d = 10.0_dp - 8.0_dp*f1 + 1.5_dp*f2
3426 s3d = -2.0_dp + 4.0_dp*g1 - 1.5_dp*g2
3427 s4d = 0.5_dp*h2
3428 v1d = -8.0_dp + 4.0_dp*p1 - 0.5_dp*p2
3429 v2d = 10.0_dp - 8.0_dp*q1 + 1.5_dp*q2
3430 v3d = -2.0_dp + 4.0_dp*r1 - 1.5_dp*r2
3431 v4d = 0.5_dp*u2
3432
3433 t1 = t1d*dr1i
3434 t2 = t2d*dr1i
3435 t3 = t3d*dr1i
3436 t4 = t4d*dr1i
3437 s1 = s1o
3438 s2 = s2o
3439 s3 = s3o
3440 s4 = s4o
3441 v1 = v1o
3442 v2 = v2o
3443 v3 = v3o
3444 v4 = v4o
3445 val(1) = ((box(1, 1, 1)*t1 + box(2, 1, 1)*t2 + box(3, 1, 1)*t3 + box(4, 1, 1)*t4)*s1 + &
3446 (box(1, 2, 1)*t1 + box(2, 2, 1)*t2 + box(3, 2, 1)*t3 + box(4, 2, 1)*t4)*s2 + &
3447 (box(1, 3, 1)*t1 + box(2, 3, 1)*t2 + box(3, 3, 1)*t3 + box(4, 3, 1)*t4)*s3 + &
3448 (box(1, 4, 1)*t1 + box(2, 4, 1)*t2 + box(3, 4, 1)*t3 + box(4, 4, 1)*t4)*s4)*v1 + &
3449 ((box(1, 1, 2)*t1 + box(2, 1, 2)*t2 + box(3, 1, 2)*t3 + box(4, 1, 2)*t4)*s1 + &
3450 (box(1, 2, 2)*t1 + box(2, 2, 2)*t2 + box(3, 2, 2)*t3 + box(4, 2, 2)*t4)*s2 + &
3451 (box(1, 3, 2)*t1 + box(2, 3, 2)*t2 + box(3, 3, 2)*t3 + box(4, 3, 2)*t4)*s3 + &
3452 (box(1, 4, 2)*t1 + box(2, 4, 2)*t2 + box(3, 4, 2)*t3 + box(4, 4, 2)*t4)*s4)*v2 + &
3453 ((box(1, 1, 3)*t1 + box(2, 1, 3)*t2 + box(3, 1, 3)*t3 + box(4, 1, 3)*t4)*s1 + &
3454 (box(1, 2, 3)*t1 + box(2, 2, 3)*t2 + box(3, 2, 3)*t3 + box(4, 2, 3)*t4)*s2 + &
3455 (box(1, 3, 3)*t1 + box(2, 3, 3)*t2 + box(3, 3, 3)*t3 + box(4, 3, 3)*t4)*s3 + &
3456 (box(1, 4, 3)*t1 + box(2, 4, 3)*t2 + box(3, 4, 3)*t3 + box(4, 4, 3)*t4)*s4)*v3 + &
3457 ((box(1, 1, 4)*t1 + box(2, 1, 4)*t2 + box(3, 1, 4)*t3 + box(4, 1, 4)*t4)*s1 + &
3458 (box(1, 2, 4)*t1 + box(2, 2, 4)*t2 + box(3, 2, 4)*t3 + box(4, 2, 4)*t4)*s2 + &
3459 (box(1, 3, 4)*t1 + box(2, 3, 4)*t2 + box(3, 3, 4)*t3 + box(4, 3, 4)*t4)*s3 + &
3460 (box(1, 4, 4)*t1 + box(2, 4, 4)*t2 + box(3, 4, 4)*t3 + box(4, 4, 4)*t4)*s4)*v4
3461
3462 t1 = t1o
3463 t2 = t2o
3464 t3 = t3o
3465 t4 = t4o
3466 s1 = s1d*dr2i
3467 s2 = s2d*dr2i
3468 s3 = s3d*dr2i
3469 s4 = s4d*dr2i
3470 v1 = v1o
3471 v2 = v2o
3472 v3 = v3o
3473 v4 = v4o
3474 val(2) = ((box(1, 1, 1)*t1 + box(2, 1, 1)*t2 + box(3, 1, 1)*t3 + box(4, 1, 1)*t4)*s1 + &
3475 (box(1, 2, 1)*t1 + box(2, 2, 1)*t2 + box(3, 2, 1)*t3 + box(4, 2, 1)*t4)*s2 + &
3476 (box(1, 3, 1)*t1 + box(2, 3, 1)*t2 + box(3, 3, 1)*t3 + box(4, 3, 1)*t4)*s3 + &
3477 (box(1, 4, 1)*t1 + box(2, 4, 1)*t2 + box(3, 4, 1)*t3 + box(4, 4, 1)*t4)*s4)*v1 + &
3478 ((box(1, 1, 2)*t1 + box(2, 1, 2)*t2 + box(3, 1, 2)*t3 + box(4, 1, 2)*t4)*s1 + &
3479 (box(1, 2, 2)*t1 + box(2, 2, 2)*t2 + box(3, 2, 2)*t3 + box(4, 2, 2)*t4)*s2 + &
3480 (box(1, 3, 2)*t1 + box(2, 3, 2)*t2 + box(3, 3, 2)*t3 + box(4, 3, 2)*t4)*s3 + &
3481 (box(1, 4, 2)*t1 + box(2, 4, 2)*t2 + box(3, 4, 2)*t3 + box(4, 4, 2)*t4)*s4)*v2 + &
3482 ((box(1, 1, 3)*t1 + box(2, 1, 3)*t2 + box(3, 1, 3)*t3 + box(4, 1, 3)*t4)*s1 + &
3483 (box(1, 2, 3)*t1 + box(2, 2, 3)*t2 + box(3, 2, 3)*t3 + box(4, 2, 3)*t4)*s2 + &
3484 (box(1, 3, 3)*t1 + box(2, 3, 3)*t2 + box(3, 3, 3)*t3 + box(4, 3, 3)*t4)*s3 + &
3485 (box(1, 4, 3)*t1 + box(2, 4, 3)*t2 + box(3, 4, 3)*t3 + box(4, 4, 3)*t4)*s4)*v3 + &
3486 ((box(1, 1, 4)*t1 + box(2, 1, 4)*t2 + box(3, 1, 4)*t3 + box(4, 1, 4)*t4)*s1 + &
3487 (box(1, 2, 4)*t1 + box(2, 2, 4)*t2 + box(3, 2, 4)*t3 + box(4, 2, 4)*t4)*s2 + &
3488 (box(1, 3, 4)*t1 + box(2, 3, 4)*t2 + box(3, 3, 4)*t3 + box(4, 3, 4)*t4)*s3 + &
3489 (box(1, 4, 4)*t1 + box(2, 4, 4)*t2 + box(3, 4, 4)*t3 + box(4, 4, 4)*t4)*s4)*v4
3490
3491 t1 = t1o
3492 t2 = t2o
3493 t3 = t3o
3494 t4 = t4o
3495 s1 = s1o
3496 s2 = s2o
3497 s3 = s3o
3498 s4 = s4o
3499 v1 = v1d*dr3i
3500 v2 = v2d*dr3i
3501 v3 = v3d*dr3i
3502 v4 = v4d*dr3i
3503 val(3) = ((box(1, 1, 1)*t1 + box(2, 1, 1)*t2 + box(3, 1, 1)*t3 + box(4, 1, 1)*t4)*s1 + &
3504 (box(1, 2, 1)*t1 + box(2, 2, 1)*t2 + box(3, 2, 1)*t3 + box(4, 2, 1)*t4)*s2 + &
3505 (box(1, 3, 1)*t1 + box(2, 3, 1)*t2 + box(3, 3, 1)*t3 + box(4, 3, 1)*t4)*s3 + &
3506 (box(1, 4, 1)*t1 + box(2, 4, 1)*t2 + box(3, 4, 1)*t3 + box(4, 4, 1)*t4)*s4)*v1 + &
3507 ((box(1, 1, 2)*t1 + box(2, 1, 2)*t2 + box(3, 1, 2)*t3 + box(4, 1, 2)*t4)*s1 + &
3508 (box(1, 2, 2)*t1 + box(2, 2, 2)*t2 + box(3, 2, 2)*t3 + box(4, 2, 2)*t4)*s2 + &
3509 (box(1, 3, 2)*t1 + box(2, 3, 2)*t2 + box(3, 3, 2)*t3 + box(4, 3, 2)*t4)*s3 + &
3510 (box(1, 4, 2)*t1 + box(2, 4, 2)*t2 + box(3, 4, 2)*t3 + box(4, 4, 2)*t4)*s4)*v2 + &
3511 ((box(1, 1, 3)*t1 + box(2, 1, 3)*t2 + box(3, 1, 3)*t3 + box(4, 1, 3)*t4)*s1 + &
3512 (box(1, 2, 3)*t1 + box(2, 2, 3)*t2 + box(3, 2, 3)*t3 + box(4, 2, 3)*t4)*s2 + &
3513 (box(1, 3, 3)*t1 + box(2, 3, 3)*t2 + box(3, 3, 3)*t3 + box(4, 3, 3)*t4)*s3 + &
3514 (box(1, 4, 3)*t1 + box(2, 4, 3)*t2 + box(3, 4, 3)*t3 + box(4, 4, 3)*t4)*s4)*v3 + &
3515 ((box(1, 1, 4)*t1 + box(2, 1, 4)*t2 + box(3, 1, 4)*t3 + box(4, 1, 4)*t4)*s1 + &
3516 (box(1, 2, 4)*t1 + box(2, 2, 4)*t2 + box(3, 2, 4)*t3 + box(4, 2, 4)*t4)*s2 + &
3517 (box(1, 3, 4)*t1 + box(2, 3, 4)*t2 + box(3, 3, 4)*t3 + box(4, 3, 4)*t4)*s3 + &
3518 (box(1, 4, 4)*t1 + box(2, 4, 4)*t2 + box(3, 4, 4)*t3 + box(4, 4, 4)*t4)*s4)*v4
3519
3520 IF (my_mpsum) CALL pw%pw_grid%para%group%sum(val)
3521
3522 END FUNCTION eval_d_interp_spl3_pbc
3523
3524END MODULE pw_spline_utils
3525
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:125
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_integral_ab
Definition pw_methods.F:222
pw_methods::pw_zero
Definition pw_methods.F:83
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:825
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:27
pw_grid_types::fullspace
integer, parameter, public fullspace
Definition pw_grid_types.F:26
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::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:1770
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:2584
pw_spline_utils::pw_spline_precond_create
subroutine, public pw_spline_precond_create(preconditioner, precond_kind, pool, pbc, transpose)
...
Definition pw_spline_utils.F:2478
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:2598
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 preconditioner to use
Definition pw_spline_utils.F:2501
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:3303
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:3154
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:3100
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:2652
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:3116
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:3132
pw_spline_utils::no_precond
integer, parameter, public no_precond
Definition pw_spline_utils.F:44
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
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
pw_pool_types::pw_pool_type
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:83
pw_spline_utils::pw_spline_precond_type
stores information for the preconditioner used to calculate the coeffs of splines
Definition pw_spline_utils.F:107
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62