Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Adaptive Clenshaw-Curtis quadrature algorithm to integrate a complex-valued function in
10 : !> a complex plane
11 : !> \par History
12 : !> * 05.2017 created [Sergey Chulkov]
13 : ! **************************************************************************************************
14 : MODULE negf_integr_cc
15 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale,&
16 : cp_cfm_scale_and_add
17 : USE cp_cfm_types, ONLY: cp_cfm_create,&
18 : cp_cfm_get_info,&
19 : cp_cfm_release,&
20 : cp_cfm_type
21 : USE cp_fm_basic_linalg, ONLY: cp_fm_trace
22 : USE cp_fm_struct, ONLY: cp_fm_struct_equivalent,&
23 : cp_fm_struct_type
24 : USE cp_fm_types, ONLY: cp_fm_create,&
25 : cp_fm_get_info,&
26 : cp_fm_release,&
27 : cp_fm_type
28 : USE fft_tools, ONLY: fft_alloc,&
29 : fft_dealloc,&
30 : fft_fw1d
31 : USE kahan_sum, ONLY: accurate_sum
32 : USE kinds, ONLY: dp,&
33 : int_8
34 : USE mathconstants, ONLY: z_one,&
35 : z_zero
36 : USE negf_integr_utils, ONLY: contour_shape_arc,&
37 : contour_shape_linear,&
38 : equidistant_nodes_a_b,&
39 : rescale_nodes_cos,&
40 : rescale_normalised_nodes
41 : #include "./base/base_uses.f90"
42 :
43 : IMPLICIT NONE
44 : PRIVATE
45 :
46 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_integr_cc'
47 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
48 :
49 : INTEGER, PARAMETER, PUBLIC :: cc_interval_full = 0, &
50 : cc_interval_half = 1
51 :
52 : INTEGER, PARAMETER, PUBLIC :: cc_shape_linear = contour_shape_linear, &
53 : cc_shape_arc = contour_shape_arc
54 :
55 : PUBLIC :: ccquad_type
56 :
57 : PUBLIC :: ccquad_init, &
58 : ccquad_release, &
59 : ccquad_double_number_of_points, &
60 : ccquad_reduce_and_append_zdata, &
61 : ccquad_refine_integral
62 :
63 : ! **************************************************************************************************
64 : !> \brief Adaptive Clenshaw-Curtis environment.
65 : ! **************************************************************************************************
66 : TYPE ccquad_type
67 : !> integration lower and upper bounds
68 : COMPLEX(kind=dp) :: a = z_zero, b = z_zero
69 : !> integration interval:
70 : !> cc_interval_full -- [a .. b],
71 : !> grid density: 'a' .. . . . . . .. 'b';
72 : !> cc_interval_half -- [a .. 2b-a], assuming int_{b}^{2b-a} f(x) dx = 0,
73 : !> grid density: 'a' .. . . . 'b'
74 : INTEGER :: interval_id = -1
75 : !> integration shape
76 : INTEGER :: shape_id = -1
77 : !> estimated error
78 : REAL(kind=dp) :: error = -1.0_dp
79 : !> approximate integral value
80 : TYPE(cp_cfm_type), POINTER :: integral => NULL()
81 : !> error estimate for every element of the 'integral' matrix
82 : TYPE(cp_fm_type), POINTER :: error_fm => NULL()
83 : !> weights associated with matrix elements; the 'error' variable contains the value Trace(error_fm * weights)
84 : TYPE(cp_fm_type), POINTER :: weights => NULL()
85 : !> integrand value at grid points. Due to symmetry of Clenshaw-Curtis quadratures,
86 : !> we only need to keep the left half-interval
87 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: zdata_cache
88 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tnodes
89 : END TYPE ccquad_type
90 :
91 : CONTAINS
92 :
93 : ! **************************************************************************************************
94 : !> \brief Initialise a Clenshaw-Curtis quadrature environment variable.
95 : !> \param cc_env environment variable to initialise
96 : !> \param xnodes points at which an integrand needs to be computed (initialised on exit)
97 : !> \param nnodes initial number of points to compute (initialised on exit)
98 : !> \param a integral lower bound
99 : !> \param b integral upper bound
100 : !> \param interval_id full [-1 .. 1] or half [-1 .. 0] interval
101 : !> \param shape_id shape of a curve along which the integral will be evaluated
102 : !> \param weights weights associated with matrix elements; used to compute cumulative error
103 : !> \param tnodes_restart list of nodes over the interval [-1 .. 1] from a previous integral evaluation.
104 : !> If present, the same set of 'xnodes' will be used to compute this integral.
105 : !> \par History
106 : !> * 05.2017 created [Sergey Chulkov]
107 : !> \note Clenshaw-Curtis quadratures are defined on the interval [-1 .. 1] and have non-uniforms node
108 : !> distribution which is symmetric and much sparse about 0. When the half-interval [-1 .. 0]
109 : !> is requested, the integrand value on another subinterval (0 .. 1] is assumed to be zero.
110 : !> Half interval mode is typically useful for rapidly decaying integrands (e.g. multiplied by
111 : !> Fermi function), so we do not actually need a fine grid spacing on this tail.
112 : ! **************************************************************************************************
113 160 : SUBROUTINE ccquad_init(cc_env, xnodes, nnodes, a, b, interval_id, shape_id, weights, tnodes_restart)
114 : TYPE(ccquad_type), INTENT(out) :: cc_env
115 : INTEGER, INTENT(inout) :: nnodes
116 : COMPLEX(kind=dp), DIMENSION(nnodes), INTENT(out) :: xnodes
117 : COMPLEX(kind=dp), INTENT(in) :: a, b
118 : INTEGER, INTENT(in) :: interval_id, shape_id
119 : TYPE(cp_fm_type), INTENT(IN) :: weights
120 : REAL(kind=dp), DIMENSION(nnodes), INTENT(in), &
121 : OPTIONAL :: tnodes_restart
122 :
123 : CHARACTER(len=*), PARAMETER :: routineN = 'ccquad_init'
124 :
125 : INTEGER :: handle, icol, ipoint, irow, ncols, &
126 : nnodes_half, nrows
127 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
128 112 : POINTER :: w_data, w_data_my
129 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
130 :
131 112 : CALL timeset(routineN, handle)
132 :
133 112 : CPASSERT(nnodes > 2)
134 :
135 : ! ensure that MOD(nnodes-1, 2) == 0
136 112 : nnodes = 2*((nnodes - 1)/2) + 1
137 :
138 112 : cc_env%interval_id = interval_id
139 112 : cc_env%shape_id = shape_id
140 112 : cc_env%a = a
141 112 : cc_env%b = b
142 112 : cc_env%error = HUGE(0.0_dp)
143 :
144 112 : NULLIFY (cc_env%integral, cc_env%error_fm, cc_env%weights)
145 112 : ALLOCATE (cc_env%weights)
146 112 : CALL cp_fm_get_info(weights, local_data=w_data, nrow_local=nrows, ncol_local=ncols, matrix_struct=fm_struct)
147 112 : CALL cp_fm_create(cc_env%weights, fm_struct)
148 112 : CALL cp_fm_get_info(cc_env%weights, local_data=w_data_my)
149 :
150 : ! use the explicit loop to avoid temporary arrays
151 1904 : DO icol = 1, ncols
152 19824 : DO irow = 1, nrows
153 19712 : w_data_my(irow, icol) = ABS(w_data(irow, icol))
154 : END DO
155 : END DO
156 :
157 56 : SELECT CASE (interval_id)
158 : CASE (cc_interval_full)
159 56 : nnodes_half = nnodes/2 + 1
160 : CASE (cc_interval_half)
161 56 : nnodes_half = nnodes
162 : CASE DEFAULT
163 112 : CPABORT("Unimplemented interval type")
164 : END SELECT
165 :
166 336 : ALLOCATE (cc_env%tnodes(nnodes))
167 :
168 112 : IF (PRESENT(tnodes_restart)) THEN
169 2112 : cc_env%tnodes(1:nnodes) = tnodes_restart(1:nnodes)
170 : ELSE
171 64 : CALL equidistant_nodes_a_b(-1.0_dp, 0.0_dp, nnodes_half, cc_env%tnodes)
172 :
173 : ! rescale all but the end-points, as they are transformed into themselves (-1.0 -> -1.0; 0.0 -> 0.0).
174 : ! Moreover, by applying this rescaling transformation to the end-points we cannot guarantee the exact
175 : ! result due to rounding errors in evaluation of COS function.
176 64 : IF (nnodes_half > 2) &
177 64 : CALL rescale_nodes_cos(nnodes_half - 2, cc_env%tnodes(2:))
178 :
179 32 : SELECT CASE (interval_id)
180 : CASE (cc_interval_full)
181 : ! reflect symmetric nodes
182 256 : DO ipoint = nnodes_half - 1, 1, -1
183 256 : cc_env%tnodes(nnodes_half + ipoint) = -cc_env%tnodes(nnodes_half - ipoint)
184 : END DO
185 : CASE (cc_interval_half)
186 : ! rescale half-interval : [-1 .. 0] -> [-1 .. 1]
187 544 : cc_env%tnodes(1:nnodes_half) = 2.0_dp*cc_env%tnodes(1:nnodes_half) + 1.0_dp
188 : END SELECT
189 : END IF
190 :
191 112 : CALL rescale_normalised_nodes(nnodes, cc_env%tnodes, a, b, shape_id, xnodes)
192 :
193 112 : CALL timestop(handle)
194 384 : END SUBROUTINE ccquad_init
195 :
196 : ! **************************************************************************************************
197 : !> \brief Release a Clenshaw-Curtis quadrature environment variable.
198 : !> \param cc_env environment variable to release (modified on exit)
199 : !> \par History
200 : !> * 05.2017 created [Sergey Chulkov]
201 : ! **************************************************************************************************
202 112 : SUBROUTINE ccquad_release(cc_env)
203 : TYPE(ccquad_type), INTENT(inout) :: cc_env
204 :
205 : CHARACTER(len=*), PARAMETER :: routineN = 'ccquad_release'
206 :
207 : INTEGER :: handle, ipoint
208 :
209 112 : CALL timeset(routineN, handle)
210 :
211 112 : IF (ASSOCIATED(cc_env%error_fm)) THEN
212 112 : CALL cp_fm_release(cc_env%error_fm)
213 112 : DEALLOCATE (cc_env%error_fm)
214 : NULLIFY (cc_env%error_fm)
215 : END IF
216 :
217 112 : IF (ASSOCIATED(cc_env%weights)) THEN
218 112 : CALL cp_fm_release(cc_env%weights)
219 112 : DEALLOCATE (cc_env%weights)
220 : NULLIFY (cc_env%weights)
221 : END IF
222 :
223 112 : IF (ASSOCIATED(cc_env%integral)) THEN
224 112 : CALL cp_cfm_release(cc_env%integral)
225 112 : DEALLOCATE (cc_env%integral)
226 : NULLIFY (cc_env%integral)
227 : END IF
228 :
229 112 : IF (ALLOCATED(cc_env%zdata_cache)) THEN
230 3360 : DO ipoint = SIZE(cc_env%zdata_cache), 1, -1
231 3360 : CALL cp_cfm_release(cc_env%zdata_cache(ipoint))
232 : END DO
233 :
234 112 : DEALLOCATE (cc_env%zdata_cache)
235 : END IF
236 :
237 112 : IF (ALLOCATED(cc_env%tnodes)) DEALLOCATE (cc_env%tnodes)
238 :
239 112 : CALL timestop(handle)
240 112 : END SUBROUTINE ccquad_release
241 :
242 : ! **************************************************************************************************
243 : !> \brief Get the next set of points at which the integrand needs to be computed. These points are
244 : !> then can be used to refine the integral approximation.
245 : !> \param cc_env environment variable (modified on exit)
246 : !> \param xnodes_next set of additional points (allocated and initialised on exit)
247 : !> \par History
248 : !> * 05.2017 created [Sergey Chulkov]
249 : ! **************************************************************************************************
250 96 : SUBROUTINE ccquad_double_number_of_points(cc_env, xnodes_next)
251 : TYPE(ccquad_type), INTENT(inout) :: cc_env
252 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:), &
253 : INTENT(inout) :: xnodes_next
254 :
255 : CHARACTER(len=*), PARAMETER :: routineN = 'ccquad_double_number_of_points'
256 :
257 : INTEGER :: handle, ipoint, nnodes_exist, &
258 : nnodes_half, nnodes_next
259 96 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tnodes, tnodes_old
260 :
261 96 : CALL timeset(routineN, handle)
262 :
263 96 : CPASSERT(.NOT. ALLOCATED(xnodes_next))
264 96 : CPASSERT(ASSOCIATED(cc_env%integral))
265 96 : CPASSERT(ASSOCIATED(cc_env%error_fm))
266 96 : CPASSERT(ALLOCATED(cc_env%zdata_cache))
267 :
268 : ! due to symmetry of Clenshaw-Curtis quadratures, we only need to keep the left half-interval [-1 .. 0]
269 96 : nnodes_exist = SIZE(cc_env%zdata_cache)
270 : ! new nodes will be placed between the existed ones, so the number of nodes
271 : ! on the left half-interval [-1 .. 0] is equal to nnodes_exist - 1
272 96 : nnodes_half = nnodes_exist - 1
273 :
274 160 : SELECT CASE (cc_env%interval_id)
275 : CASE (cc_interval_full)
276 : ! double number of nodes as we have 2 half-intervals [-1 .. 0] and [0 .. 1]
277 64 : nnodes_next = 2*nnodes_half
278 : CASE (cc_interval_half)
279 32 : nnodes_next = nnodes_half
280 : CASE DEFAULT
281 96 : CPABORT("Unimplemented interval type")
282 : END SELECT
283 :
284 288 : ALLOCATE (xnodes_next(nnodes_next))
285 288 : ALLOCATE (tnodes(nnodes_next))
286 :
287 : CALL equidistant_nodes_a_b(0.5_dp/REAL(nnodes_half, kind=dp) - 1.0_dp, &
288 : -0.5_dp/REAL(nnodes_half, kind=dp), &
289 96 : nnodes_half, tnodes)
290 :
291 96 : CALL rescale_nodes_cos(nnodes_half, tnodes)
292 :
293 96 : SELECT CASE (cc_env%interval_id)
294 : CASE (cc_interval_full)
295 : ! reflect symmetric nodes
296 736 : DO ipoint = 1, nnodes_half
297 736 : tnodes(nnodes_half + ipoint) = -tnodes(nnodes_half - ipoint + 1)
298 : END DO
299 : CASE (cc_interval_half)
300 : ! rescale half-interval : [-1 .. 0] -> [-1 .. 1]
301 544 : tnodes(1:nnodes_half) = 2.0_dp*tnodes(1:nnodes_half) + 1.0_dp
302 : END SELECT
303 :
304 : ! append new tnodes to the cache
305 96 : CALL MOVE_ALLOC(cc_env%tnodes, tnodes_old)
306 96 : nnodes_exist = SIZE(tnodes_old)
307 :
308 288 : ALLOCATE (cc_env%tnodes(nnodes_exist + nnodes_next))
309 1984 : cc_env%tnodes(1:nnodes_exist) = tnodes_old(1:nnodes_exist)
310 1888 : cc_env%tnodes(nnodes_exist + 1:nnodes_exist + nnodes_next) = tnodes(1:nnodes_next)
311 96 : DEALLOCATE (tnodes_old)
312 :
313 : ! rescale nodes [-1 .. 1] -> [a .. b] according to the shape
314 96 : CALL rescale_normalised_nodes(nnodes_next, tnodes, cc_env%a, cc_env%b, cc_env%shape_id, xnodes_next)
315 :
316 96 : DEALLOCATE (tnodes)
317 96 : CALL timestop(handle)
318 96 : END SUBROUTINE ccquad_double_number_of_points
319 :
320 : ! **************************************************************************************************
321 : !> \brief Prepare Clenshaw-Curtis environment for the subsequent refinement of the integral.
322 : !> \param cc_env environment variable (modified on exit)
323 : !> \param zdata_next additional integrand value at additional points (modified on exit)
324 : !> \par History
325 : !> * 05.2017 created [Sergey Chulkov]
326 : !> \note Due to symmetry of Clenshaw-Curtis quadratures (weight(x) == weight(-x)), we do not need to
327 : !> keep all the matrices from 'zdata_next', only 'zdata_next(x) + zdata_next(-x)' is needed.
328 : !> In order to reduce the number of matrix allocations, we move some of the matrices from the
329 : !> end of the 'zdata_new' array to the 'cc_env%zdata_cache' array, and nullify the corresponding
330 : !> pointers at 'zdata_next' array. So the calling subroutine need to release the remained
331 : !> matrices or reuse them but taking into account the missed ones.
332 : ! **************************************************************************************************
333 208 : SUBROUTINE ccquad_reduce_and_append_zdata(cc_env, zdata_next)
334 : TYPE(ccquad_type), INTENT(inout) :: cc_env
335 : TYPE(cp_cfm_type), DIMENSION(:), INTENT(inout) :: zdata_next
336 :
337 : CHARACTER(len=*), PARAMETER :: routineN = 'ccquad_reduce_and_append_zdata'
338 : TYPE(cp_cfm_type), PARAMETER :: cfm_null = cp_cfm_type()
339 :
340 208 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: zscale
341 : INTEGER :: handle, ipoint, nnodes_exist, &
342 : nnodes_half, nnodes_next
343 208 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: zdata_tmp
344 :
345 208 : CALL timeset(routineN, handle)
346 :
347 208 : nnodes_next = SIZE(zdata_next)
348 208 : CPASSERT(nnodes_next > 0)
349 :
350 : ! compute weights of new points on a complex contour according to their values of the 't' parameter
351 208 : nnodes_exist = SIZE(cc_env%tnodes)
352 208 : CPASSERT(nnodes_exist >= nnodes_next)
353 :
354 624 : ALLOCATE (zscale(nnodes_next))
355 : CALL rescale_normalised_nodes(nnodes_next, cc_env%tnodes(nnodes_exist - nnodes_next + 1:nnodes_exist), &
356 208 : cc_env%a, cc_env%b, cc_env%shape_id, weights=zscale)
357 :
358 1920 : IF (cc_env%interval_id == cc_interval_half) zscale(:) = 2.0_dp*zscale(:)
359 :
360 : ! rescale integrand values
361 5024 : DO ipoint = 1, nnodes_next
362 5024 : CALL cp_cfm_scale(zscale(ipoint), zdata_next(ipoint))
363 : END DO
364 208 : DEALLOCATE (zscale)
365 :
366 : ! squash points with the same clenshaw-curtis weights together
367 208 : IF (ALLOCATED(cc_env%zdata_cache)) THEN
368 96 : nnodes_exist = SIZE(cc_env%zdata_cache)
369 : ELSE
370 : nnodes_exist = 0
371 : END IF
372 :
373 328 : SELECT CASE (cc_env%interval_id)
374 : CASE (cc_interval_full)
375 120 : IF (ALLOCATED(cc_env%zdata_cache)) THEN
376 64 : CPASSERT(nnodes_exist == nnodes_next/2 + 1)
377 64 : nnodes_half = nnodes_exist - 1
378 : ELSE
379 56 : CPASSERT(MOD(nnodes_next, 2) == 1)
380 56 : nnodes_half = nnodes_next/2 + 1
381 : END IF
382 : CASE (cc_interval_half)
383 88 : IF (ALLOCATED(cc_env%zdata_cache)) THEN
384 32 : CPASSERT(nnodes_exist == nnodes_next + 1)
385 : END IF
386 :
387 208 : nnodes_half = nnodes_next
388 : END SELECT
389 :
390 208 : IF (cc_env%interval_id == cc_interval_full) THEN
391 1688 : DO ipoint = nnodes_next/2, 1, -1
392 1688 : CALL cp_cfm_scale_and_add(z_one, zdata_next(ipoint), z_one, zdata_next(nnodes_next - ipoint + 1))
393 : END DO
394 : END IF
395 :
396 208 : IF (ALLOCATED(cc_env%zdata_cache)) THEN
397 : ! note that nnodes_half+1 == nnodes_exist for both half- and full-intervals
398 2624 : ALLOCATE (zdata_tmp(nnodes_half + nnodes_exist))
399 :
400 1216 : DO ipoint = 1, nnodes_half
401 1120 : zdata_tmp(2*ipoint - 1) = cc_env%zdata_cache(ipoint)
402 1120 : zdata_tmp(2*ipoint) = zdata_next(ipoint)
403 1216 : zdata_next(ipoint) = cfm_null
404 : END DO
405 96 : zdata_tmp(nnodes_half + nnodes_exist) = cc_env%zdata_cache(nnodes_exist)
406 :
407 96 : CALL MOVE_ALLOC(zdata_tmp, cc_env%zdata_cache)
408 : ELSE
409 112 : CALL cp_cfm_scale(2.0_dp, zdata_next(nnodes_half))
410 :
411 2464 : ALLOCATE (cc_env%zdata_cache(nnodes_half))
412 :
413 2240 : DO ipoint = 1, nnodes_half
414 2128 : cc_env%zdata_cache(ipoint) = zdata_next(ipoint)
415 2240 : zdata_next(ipoint) = cfm_null
416 : END DO
417 : END IF
418 :
419 208 : CALL timestop(handle)
420 208 : END SUBROUTINE ccquad_reduce_and_append_zdata
421 :
422 : ! **************************************************************************************************
423 : !> \brief Refine approximated integral.
424 : !> \param cc_env environment variable (modified on exit)
425 : !> \par History
426 : !> * 05.2017 created [Sergey Chulkov]
427 : ! **************************************************************************************************
428 208 : SUBROUTINE ccquad_refine_integral(cc_env)
429 : TYPE(ccquad_type), INTENT(inout) :: cc_env
430 :
431 : CHARACTER(len=*), PARAMETER :: routineN = 'ccquad_refine_integral'
432 :
433 : COMPLEX(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
434 208 : POINTER :: ztmp, ztmp_dct
435 : INTEGER :: handle, icol, ipoint, irow, ncols_local, nintervals, nintervals_half, &
436 : nintervals_half_plus_1, nintervals_half_plus_2, nintervals_plus_2, nrows_local, stat
437 : LOGICAL :: equiv
438 : REAL(kind=dp) :: rscale
439 208 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: weights
440 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
441 :
442 : ! TYPE(fft_plan_type) :: fft_plan
443 : ! INTEGER(kind=int_8) :: plan
444 :
445 208 : CALL timeset(routineN, handle)
446 :
447 208 : CPASSERT(ALLOCATED(cc_env%zdata_cache))
448 :
449 208 : nintervals_half_plus_1 = SIZE(cc_env%zdata_cache)
450 208 : nintervals_half = nintervals_half_plus_1 - 1
451 208 : nintervals_half_plus_2 = nintervals_half_plus_1 + 1
452 208 : nintervals = 2*nintervals_half
453 208 : nintervals_plus_2 = nintervals + 2
454 208 : CPASSERT(nintervals_half > 1)
455 :
456 208 : IF (.NOT. ASSOCIATED(cc_env%integral)) THEN
457 112 : CALL cp_cfm_get_info(cc_env%zdata_cache(1), matrix_struct=fm_struct)
458 112 : equiv = cp_fm_struct_equivalent(fm_struct, cc_env%weights%matrix_struct)
459 112 : CPASSERT(equiv)
460 :
461 112 : ALLOCATE (cc_env%integral)
462 112 : CALL cp_cfm_create(cc_env%integral, fm_struct)
463 : NULLIFY (cc_env%error_fm)
464 112 : ALLOCATE (cc_env%error_fm)
465 112 : CALL cp_fm_create(cc_env%error_fm, fm_struct)
466 : END IF
467 :
468 : IF (debug_this_module) THEN
469 4672 : DO ipoint = 1, nintervals_half_plus_1
470 4464 : equiv = cp_fm_struct_equivalent(cc_env%zdata_cache(ipoint)%matrix_struct, cc_env%integral%matrix_struct)
471 4672 : CPASSERT(equiv)
472 : END DO
473 : END IF
474 :
475 208 : CALL cp_cfm_get_info(cc_env%integral, nrow_local=nrows_local, ncol_local=ncols_local)
476 :
477 624 : ALLOCATE (weights(nintervals_half))
478 :
479 : ! omit the trivial weights(1) = 0.5
480 4256 : DO ipoint = 2, nintervals_half
481 4048 : rscale = REAL(2*(ipoint - 1), kind=dp)
482 4256 : weights(ipoint) = 1.0_dp/(1.0_dp - rscale*rscale)
483 : END DO
484 : ! weights(1) <- weights(intervals_half + 1)
485 208 : rscale = REAL(nintervals, kind=dp)
486 208 : weights(1) = 1.0_dp/(1.0_dp - rscale*rscale)
487 :
488 : ! 1.0 / nintervals
489 208 : rscale = 1.0_dp/rscale
490 :
491 832 : CALL fft_alloc(ztmp, [nintervals, nrows_local, ncols_local])
492 832 : CALL fft_alloc(ztmp_dct, [nintervals, nrows_local, ncols_local])
493 :
494 : !$OMP PARALLEL DO DEFAULT(NONE), PRIVATE(icol, ipoint, irow), &
495 208 : !$OMP SHARED(cc_env, ncols_local, nintervals_half, nintervals_half_plus_1, nintervals_half_plus_2, nrows_local, ztmp)
496 : DO icol = 1, ncols_local
497 : DO irow = 1, nrows_local
498 : DO ipoint = 1, nintervals_half_plus_1
499 : ztmp(ipoint, irow, icol) = cc_env%zdata_cache(ipoint)%local_data(irow, icol)
500 : END DO
501 :
502 : DO ipoint = 2, nintervals_half
503 : ztmp(nintervals_half + ipoint, irow, icol) = ztmp(nintervals_half_plus_2 - ipoint, irow, icol)
504 : END DO
505 : END DO
506 : END DO
507 : !$OMP END PARALLEL DO
508 :
509 208 : CALL fft_fw1d(nintervals, nrows_local*ncols_local, .FALSE., ztmp, ztmp_dct, 1.0_dp, stat)
510 208 : IF (stat /= 0) THEN
511 : CALL cp_abort(__LOCATION__, &
512 : "An FFT library is required for Clenshaw-Curtis quadrature. "// &
513 0 : "You can use an alternative integration method instead.")
514 : END IF
515 :
516 : !$OMP PARALLEL DO DEFAULT(NONE), PRIVATE(icol, ipoint, irow), &
517 : !$OMP SHARED(cc_env, rscale, ncols_local, nintervals_half, nintervals_half_plus_1, nintervals_plus_2), &
518 208 : !$OMP SHARED(nrows_local, weights, ztmp_dct)
519 : DO icol = 1, ncols_local
520 : DO irow = 1, nrows_local
521 : ztmp_dct(1, irow, icol) = 0.5_dp*ztmp_dct(1, irow, icol)
522 : DO ipoint = 2, nintervals_half
523 : ztmp_dct(ipoint, irow, icol) = 0.5_dp*weights(ipoint)*(ztmp_dct(ipoint, irow, icol) + &
524 : ztmp_dct(nintervals_plus_2 - ipoint, irow, icol))
525 : END DO
526 : ztmp_dct(nintervals_half_plus_1, irow, icol) = weights(1)*ztmp_dct(nintervals_half_plus_1, irow, icol)
527 :
528 : cc_env%integral%local_data(irow, icol) = rscale*accurate_sum(ztmp_dct(1:nintervals_half_plus_1, irow, icol))
529 : cc_env%error_fm%local_data(irow, icol) = rscale*ABS(ztmp_dct(nintervals_half_plus_1, irow, icol))
530 : END DO
531 : END DO
532 : !$OMP END PARALLEL DO
533 :
534 208 : CALL fft_dealloc(ztmp)
535 208 : CALL fft_dealloc(ztmp_dct)
536 :
537 208 : CALL cp_fm_trace(cc_env%error_fm, cc_env%weights, cc_env%error)
538 :
539 208 : DEALLOCATE (weights)
540 208 : CALL timestop(handle)
541 624 : END SUBROUTINE ccquad_refine_integral
542 :
543 0 : END MODULE negf_integr_cc
|
