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