Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 0 : 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 0 : POINTER :: w_data, w_data_my
129 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
130 :
131 0 : CALL timeset(routineN, handle)
132 :
133 0 : CPASSERT(nnodes > 2)
134 :
135 : ! ensure that MOD(nnodes-1, 2) == 0
136 0 : nnodes = 2*((nnodes - 1)/2) + 1
137 :
138 0 : cc_env%interval_id = interval_id
139 0 : cc_env%shape_id = shape_id
140 0 : cc_env%a = a
141 0 : cc_env%b = b
142 0 : cc_env%error = HUGE(0.0_dp)
143 :
144 0 : NULLIFY (cc_env%integral, cc_env%error_fm, cc_env%weights)
145 0 : ALLOCATE (cc_env%weights)
146 0 : CALL cp_fm_get_info(weights, local_data=w_data, nrow_local=nrows, ncol_local=ncols, matrix_struct=fm_struct)
147 0 : CALL cp_fm_create(cc_env%weights, fm_struct)
148 0 : CALL cp_fm_get_info(cc_env%weights, local_data=w_data_my)
149 :
150 : ! use the explicit loop to avoid temporary arrays
151 0 : DO icol = 1, ncols
152 0 : DO irow = 1, nrows
153 0 : w_data_my(irow, icol) = ABS(w_data(irow, icol))
154 : END DO
155 : END DO
156 :
157 0 : SELECT CASE (interval_id)
158 : CASE (cc_interval_full)
159 0 : nnodes_half = nnodes/2 + 1
160 : CASE (cc_interval_half)
161 0 : nnodes_half = nnodes
162 : CASE DEFAULT
163 0 : CPABORT("Unimplemented interval type")
164 : END SELECT
165 :
166 0 : ALLOCATE (cc_env%tnodes(nnodes))
167 :
168 0 : IF (PRESENT(tnodes_restart)) THEN
169 0 : cc_env%tnodes(1:nnodes) = tnodes_restart(1:nnodes)
170 : ELSE
171 0 : 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 0 : IF (nnodes_half > 2) THEN
177 0 : CALL rescale_nodes_cos(nnodes_half - 2, cc_env%tnodes(2:))
178 : END IF
179 :
180 0 : SELECT CASE (interval_id)
181 : CASE (cc_interval_full)
182 : ! reflect symmetric nodes
183 0 : DO ipoint = nnodes_half - 1, 1, -1
184 0 : cc_env%tnodes(nnodes_half + ipoint) = -cc_env%tnodes(nnodes_half - ipoint)
185 : END DO
186 : CASE (cc_interval_half)
187 : ! rescale half-interval : [-1 .. 0] -> [-1 .. 1]
188 0 : cc_env%tnodes(1:nnodes_half) = 2.0_dp*cc_env%tnodes(1:nnodes_half) + 1.0_dp
189 : END SELECT
190 : END IF
191 :
192 0 : CALL rescale_normalised_nodes(nnodes, cc_env%tnodes, a, b, shape_id, xnodes)
193 :
194 0 : CALL timestop(handle)
195 0 : END SUBROUTINE ccquad_init
196 :
197 : ! **************************************************************************************************
198 : !> \brief Release a Clenshaw-Curtis quadrature environment variable.
199 : !> \param cc_env environment variable to release (modified on exit)
200 : !> \par History
201 : !> * 05.2017 created [Sergey Chulkov]
202 : ! **************************************************************************************************
203 0 : SUBROUTINE ccquad_release(cc_env)
204 : TYPE(ccquad_type), INTENT(inout) :: cc_env
205 :
206 : CHARACTER(len=*), PARAMETER :: routineN = 'ccquad_release'
207 :
208 : INTEGER :: handle, ipoint
209 :
210 0 : CALL timeset(routineN, handle)
211 :
212 0 : IF (ASSOCIATED(cc_env%error_fm)) THEN
213 0 : CALL cp_fm_release(cc_env%error_fm)
214 0 : DEALLOCATE (cc_env%error_fm)
215 : NULLIFY (cc_env%error_fm)
216 : END IF
217 :
218 0 : IF (ASSOCIATED(cc_env%weights)) THEN
219 0 : CALL cp_fm_release(cc_env%weights)
220 0 : DEALLOCATE (cc_env%weights)
221 : NULLIFY (cc_env%weights)
222 : END IF
223 :
224 0 : IF (ASSOCIATED(cc_env%integral)) THEN
225 0 : CALL cp_cfm_release(cc_env%integral)
226 0 : DEALLOCATE (cc_env%integral)
227 : NULLIFY (cc_env%integral)
228 : END IF
229 :
230 0 : IF (ALLOCATED(cc_env%zdata_cache)) THEN
231 0 : DO ipoint = SIZE(cc_env%zdata_cache), 1, -1
232 0 : CALL cp_cfm_release(cc_env%zdata_cache(ipoint))
233 : END DO
234 :
235 0 : DEALLOCATE (cc_env%zdata_cache)
236 : END IF
237 :
238 0 : IF (ALLOCATED(cc_env%tnodes)) DEALLOCATE (cc_env%tnodes)
239 :
240 0 : CALL timestop(handle)
241 0 : END SUBROUTINE ccquad_release
242 :
243 : ! **************************************************************************************************
244 : !> \brief Get the next set of points at which the integrand needs to be computed. These points are
245 : !> then can be used to refine the integral approximation.
246 : !> \param cc_env environment variable (modified on exit)
247 : !> \param xnodes_next set of additional points (allocated and initialised on exit)
248 : !> \par History
249 : !> * 05.2017 created [Sergey Chulkov]
250 : ! **************************************************************************************************
251 0 : SUBROUTINE ccquad_double_number_of_points(cc_env, xnodes_next)
252 : TYPE(ccquad_type), INTENT(inout) :: cc_env
253 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:), &
254 : INTENT(inout) :: xnodes_next
255 :
256 : CHARACTER(len=*), PARAMETER :: routineN = 'ccquad_double_number_of_points'
257 :
258 : INTEGER :: handle, ipoint, nnodes_exist, &
259 : nnodes_half, nnodes_next
260 0 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tnodes, tnodes_old
261 :
262 0 : CALL timeset(routineN, handle)
263 :
264 0 : CPASSERT(.NOT. ALLOCATED(xnodes_next))
265 0 : CPASSERT(ASSOCIATED(cc_env%integral))
266 0 : CPASSERT(ASSOCIATED(cc_env%error_fm))
267 0 : CPASSERT(ALLOCATED(cc_env%zdata_cache))
268 :
269 : ! due to symmetry of Clenshaw-Curtis quadratures, we only need to keep the left half-interval [-1 .. 0]
270 0 : nnodes_exist = SIZE(cc_env%zdata_cache)
271 : ! new nodes will be placed between the existed ones, so the number of nodes
272 : ! on the left half-interval [-1 .. 0] is equal to nnodes_exist - 1
273 0 : nnodes_half = nnodes_exist - 1
274 :
275 0 : SELECT CASE (cc_env%interval_id)
276 : CASE (cc_interval_full)
277 : ! double number of nodes as we have 2 half-intervals [-1 .. 0] and [0 .. 1]
278 0 : nnodes_next = 2*nnodes_half
279 : CASE (cc_interval_half)
280 0 : nnodes_next = nnodes_half
281 : CASE DEFAULT
282 0 : CPABORT("Unimplemented interval type")
283 : END SELECT
284 :
285 0 : ALLOCATE (xnodes_next(nnodes_next))
286 0 : ALLOCATE (tnodes(nnodes_next))
287 :
288 : CALL equidistant_nodes_a_b(0.5_dp/REAL(nnodes_half, kind=dp) - 1.0_dp, &
289 : -0.5_dp/REAL(nnodes_half, kind=dp), &
290 0 : nnodes_half, tnodes)
291 :
292 0 : CALL rescale_nodes_cos(nnodes_half, tnodes)
293 :
294 0 : SELECT CASE (cc_env%interval_id)
295 : CASE (cc_interval_full)
296 : ! reflect symmetric nodes
297 0 : DO ipoint = 1, nnodes_half
298 0 : tnodes(nnodes_half + ipoint) = -tnodes(nnodes_half - ipoint + 1)
299 : END DO
300 : CASE (cc_interval_half)
301 : ! rescale half-interval : [-1 .. 0] -> [-1 .. 1]
302 0 : tnodes(1:nnodes_half) = 2.0_dp*tnodes(1:nnodes_half) + 1.0_dp
303 : END SELECT
304 :
305 : ! append new tnodes to the cache
306 0 : CALL MOVE_ALLOC(cc_env%tnodes, tnodes_old)
307 0 : nnodes_exist = SIZE(tnodes_old)
308 :
309 0 : ALLOCATE (cc_env%tnodes(nnodes_exist + nnodes_next))
310 0 : cc_env%tnodes(1:nnodes_exist) = tnodes_old(1:nnodes_exist)
311 0 : cc_env%tnodes(nnodes_exist + 1:nnodes_exist + nnodes_next) = tnodes(1:nnodes_next)
312 0 : DEALLOCATE (tnodes_old)
313 :
314 : ! rescale nodes [-1 .. 1] -> [a .. b] according to the shape
315 0 : CALL rescale_normalised_nodes(nnodes_next, tnodes, cc_env%a, cc_env%b, cc_env%shape_id, xnodes_next)
316 :
317 0 : DEALLOCATE (tnodes)
318 0 : CALL timestop(handle)
319 0 : END SUBROUTINE ccquad_double_number_of_points
320 :
321 : ! **************************************************************************************************
322 : !> \brief Prepare Clenshaw-Curtis environment for the subsequent refinement of the integral.
323 : !> \param cc_env environment variable (modified on exit)
324 : !> \param zdata_next additional integrand value at additional points (modified on exit)
325 : !> \par History
326 : !> * 05.2017 created [Sergey Chulkov]
327 : !> \note Due to symmetry of Clenshaw-Curtis quadratures (weight(x) == weight(-x)), we do not need to
328 : !> keep all the matrices from 'zdata_next', only 'zdata_next(x) + zdata_next(-x)' is needed.
329 : !> In order to reduce the number of matrix allocations, we move some of the matrices from the
330 : !> end of the 'zdata_new' array to the 'cc_env%zdata_cache' array, and nullify the corresponding
331 : !> pointers at 'zdata_next' array. So the calling subroutine need to release the remained
332 : !> matrices or reuse them but taking into account the missed ones.
333 : ! **************************************************************************************************
334 0 : SUBROUTINE ccquad_reduce_and_append_zdata(cc_env, zdata_next)
335 : TYPE(ccquad_type), INTENT(inout) :: cc_env
336 : TYPE(cp_cfm_type), DIMENSION(:), INTENT(inout) :: zdata_next
337 :
338 : CHARACTER(len=*), PARAMETER :: routineN = 'ccquad_reduce_and_append_zdata'
339 : TYPE(cp_cfm_type), PARAMETER :: cfm_null = cp_cfm_type()
340 :
341 0 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: zscale
342 : INTEGER :: handle, ipoint, nnodes_exist, &
343 : nnodes_half, nnodes_next
344 0 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: zdata_tmp
345 :
346 0 : CALL timeset(routineN, handle)
347 :
348 0 : nnodes_next = SIZE(zdata_next)
349 0 : CPASSERT(nnodes_next > 0)
350 :
351 : ! compute weights of new points on a complex contour according to their values of the 't' parameter
352 0 : nnodes_exist = SIZE(cc_env%tnodes)
353 0 : CPASSERT(nnodes_exist >= nnodes_next)
354 :
355 0 : ALLOCATE (zscale(nnodes_next))
356 : CALL rescale_normalised_nodes(nnodes_next, cc_env%tnodes(nnodes_exist - nnodes_next + 1:nnodes_exist), &
357 0 : cc_env%a, cc_env%b, cc_env%shape_id, weights=zscale)
358 :
359 0 : IF (cc_env%interval_id == cc_interval_half) zscale(:) = 2.0_dp*zscale(:)
360 :
361 : ! rescale integrand values
362 0 : DO ipoint = 1, nnodes_next
363 0 : CALL cp_cfm_scale(zscale(ipoint), zdata_next(ipoint))
364 : END DO
365 0 : DEALLOCATE (zscale)
366 :
367 : ! squash points with the same clenshaw-curtis weights together
368 0 : IF (ALLOCATED(cc_env%zdata_cache)) THEN
369 0 : nnodes_exist = SIZE(cc_env%zdata_cache)
370 : ELSE
371 : nnodes_exist = 0
372 : END IF
373 :
374 0 : SELECT CASE (cc_env%interval_id)
375 : CASE (cc_interval_full)
376 0 : IF (ALLOCATED(cc_env%zdata_cache)) THEN
377 0 : CPASSERT(nnodes_exist == nnodes_next/2 + 1)
378 0 : nnodes_half = nnodes_exist - 1
379 : ELSE
380 0 : CPASSERT(MOD(nnodes_next, 2) == 1)
381 0 : nnodes_half = nnodes_next/2 + 1
382 : END IF
383 : CASE (cc_interval_half)
384 0 : IF (ALLOCATED(cc_env%zdata_cache)) THEN
385 0 : CPASSERT(nnodes_exist == nnodes_next + 1)
386 : END IF
387 :
388 0 : nnodes_half = nnodes_next
389 : END SELECT
390 :
391 0 : IF (cc_env%interval_id == cc_interval_full) THEN
392 0 : DO ipoint = nnodes_next/2, 1, -1
393 0 : CALL cp_cfm_scale_and_add(z_one, zdata_next(ipoint), z_one, zdata_next(nnodes_next - ipoint + 1))
394 : END DO
395 : END IF
396 :
397 0 : IF (ALLOCATED(cc_env%zdata_cache)) THEN
398 : ! note that nnodes_half+1 == nnodes_exist for both half- and full-intervals
399 0 : ALLOCATE (zdata_tmp(nnodes_half + nnodes_exist))
400 :
401 0 : DO ipoint = 1, nnodes_half
402 0 : zdata_tmp(2*ipoint - 1) = cc_env%zdata_cache(ipoint)
403 0 : zdata_tmp(2*ipoint) = zdata_next(ipoint)
404 0 : zdata_next(ipoint) = cfm_null
405 : END DO
406 0 : zdata_tmp(nnodes_half + nnodes_exist) = cc_env%zdata_cache(nnodes_exist)
407 :
408 0 : CALL MOVE_ALLOC(zdata_tmp, cc_env%zdata_cache)
409 : ELSE
410 0 : CALL cp_cfm_scale(2.0_dp, zdata_next(nnodes_half))
411 :
412 0 : ALLOCATE (cc_env%zdata_cache(nnodes_half))
413 :
414 0 : DO ipoint = 1, nnodes_half
415 0 : cc_env%zdata_cache(ipoint) = zdata_next(ipoint)
416 0 : zdata_next(ipoint) = cfm_null
417 : END DO
418 : END IF
419 :
420 0 : CALL timestop(handle)
421 0 : END SUBROUTINE ccquad_reduce_and_append_zdata
422 :
423 : ! **************************************************************************************************
424 : !> \brief Refine approximated integral.
425 : !> \param cc_env environment variable (modified on exit)
426 : !> \par History
427 : !> * 05.2017 created [Sergey Chulkov]
428 : ! **************************************************************************************************
429 0 : SUBROUTINE ccquad_refine_integral(cc_env)
430 : TYPE(ccquad_type), INTENT(inout) :: cc_env
431 :
432 : CHARACTER(len=*), PARAMETER :: routineN = 'ccquad_refine_integral'
433 :
434 : COMPLEX(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
435 0 : POINTER :: ztmp, ztmp_dct
436 : INTEGER :: handle, icol, ipoint, irow, ncols_local, nintervals, nintervals_half, &
437 : nintervals_half_plus_1, nintervals_half_plus_2, nintervals_plus_2, nrows_local, stat
438 : LOGICAL :: equiv
439 : REAL(kind=dp) :: rscale
440 0 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: weights
441 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
442 :
443 : ! TYPE(fft_plan_type) :: fft_plan
444 : ! INTEGER(kind=int_8) :: plan
445 :
446 0 : CALL timeset(routineN, handle)
447 :
448 0 : CPASSERT(ALLOCATED(cc_env%zdata_cache))
449 :
450 0 : nintervals_half_plus_1 = SIZE(cc_env%zdata_cache)
451 0 : nintervals_half = nintervals_half_plus_1 - 1
452 0 : nintervals_half_plus_2 = nintervals_half_plus_1 + 1
453 0 : nintervals = 2*nintervals_half
454 0 : nintervals_plus_2 = nintervals + 2
455 0 : CPASSERT(nintervals_half > 1)
456 :
457 0 : IF (.NOT. ASSOCIATED(cc_env%integral)) THEN
458 0 : CALL cp_cfm_get_info(cc_env%zdata_cache(1), matrix_struct=fm_struct)
459 0 : equiv = cp_fm_struct_equivalent(fm_struct, cc_env%weights%matrix_struct)
460 0 : CPASSERT(equiv)
461 :
462 0 : ALLOCATE (cc_env%integral)
463 0 : CALL cp_cfm_create(cc_env%integral, fm_struct)
464 : NULLIFY (cc_env%error_fm)
465 0 : ALLOCATE (cc_env%error_fm)
466 0 : CALL cp_fm_create(cc_env%error_fm, fm_struct)
467 : END IF
468 :
469 : IF (debug_this_module) THEN
470 0 : DO ipoint = 1, nintervals_half_plus_1
471 0 : equiv = cp_fm_struct_equivalent(cc_env%zdata_cache(ipoint)%matrix_struct, cc_env%integral%matrix_struct)
472 0 : CPASSERT(equiv)
473 : END DO
474 : END IF
475 :
476 0 : CALL cp_cfm_get_info(cc_env%integral, nrow_local=nrows_local, ncol_local=ncols_local)
477 :
478 0 : ALLOCATE (weights(nintervals_half))
479 :
480 : ! omit the trivial weights(1) = 0.5
481 0 : DO ipoint = 2, nintervals_half
482 0 : rscale = REAL(2*(ipoint - 1), kind=dp)
483 0 : weights(ipoint) = 1.0_dp/(1.0_dp - rscale*rscale)
484 : END DO
485 : ! weights(1) <- weights(intervals_half + 1)
486 0 : rscale = REAL(nintervals, kind=dp)
487 0 : weights(1) = 1.0_dp/(1.0_dp - rscale*rscale)
488 :
489 : ! 1.0 / nintervals
490 0 : rscale = 1.0_dp/rscale
491 :
492 0 : CALL fft_alloc(ztmp, [nintervals, nrows_local, ncols_local])
493 0 : CALL fft_alloc(ztmp_dct, [nintervals, nrows_local, ncols_local])
494 :
495 : !$OMP PARALLEL DO DEFAULT(NONE), PRIVATE(icol, ipoint, irow), &
496 0 : !$OMP SHARED(cc_env, ncols_local, nintervals_half, nintervals_half_plus_1, nintervals_half_plus_2, nrows_local, ztmp)
497 : DO icol = 1, ncols_local
498 : DO irow = 1, nrows_local
499 : DO ipoint = 1, nintervals_half_plus_1
500 : ztmp(ipoint, irow, icol) = cc_env%zdata_cache(ipoint)%local_data(irow, icol)
501 : END DO
502 :
503 : DO ipoint = 2, nintervals_half
504 : ztmp(nintervals_half + ipoint, irow, icol) = ztmp(nintervals_half_plus_2 - ipoint, irow, icol)
505 : END DO
506 : END DO
507 : END DO
508 : !$OMP END PARALLEL DO
509 :
510 0 : CALL fft_fw1d(nintervals, nrows_local*ncols_local, .FALSE., ztmp, ztmp_dct, 1.0_dp, stat)
511 0 : IF (stat /= 0) THEN
512 : CALL cp_abort(__LOCATION__, &
513 : "An FFT library is required for Clenshaw-Curtis quadrature. "// &
514 0 : "You can use an alternative integration method instead.")
515 : END IF
516 :
517 : !$OMP PARALLEL DO DEFAULT(NONE), PRIVATE(icol, ipoint, irow), &
518 : !$OMP SHARED(cc_env, rscale, ncols_local, nintervals_half, nintervals_half_plus_1, nintervals_plus_2), &
519 0 : !$OMP SHARED(nrows_local, weights, ztmp_dct)
520 : DO icol = 1, ncols_local
521 : DO irow = 1, nrows_local
522 : ztmp_dct(1, irow, icol) = 0.5_dp*ztmp_dct(1, irow, icol)
523 : DO ipoint = 2, nintervals_half
524 : ztmp_dct(ipoint, irow, icol) = 0.5_dp*weights(ipoint)*(ztmp_dct(ipoint, irow, icol) + &
525 : ztmp_dct(nintervals_plus_2 - ipoint, irow, icol))
526 : END DO
527 : ztmp_dct(nintervals_half_plus_1, irow, icol) = weights(1)*ztmp_dct(nintervals_half_plus_1, irow, icol)
528 :
529 : cc_env%integral%local_data(irow, icol) = rscale*accurate_sum(ztmp_dct(1:nintervals_half_plus_1, irow, icol))
530 : cc_env%error_fm%local_data(irow, icol) = rscale*ABS(ztmp_dct(nintervals_half_plus_1, irow, icol))
531 : END DO
532 : END DO
533 : !$OMP END PARALLEL DO
534 :
535 0 : CALL fft_dealloc(ztmp)
536 0 : CALL fft_dealloc(ztmp_dct)
537 :
538 0 : CALL cp_fm_trace(cc_env%error_fm, cc_env%weights, cc_env%error)
539 :
540 0 : DEALLOCATE (weights)
541 0 : CALL timestop(handle)
542 0 : END SUBROUTINE ccquad_refine_integral
543 :
544 0 : END MODULE negf_integr_cc
|
