(git:0de0cc2)
negf_integr_cc.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief 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 ! **************************************************************************************************
15  USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale,&
17  USE cp_cfm_types, ONLY: cp_cfm_create,&
20  cp_cfm_type
21  USE cp_fm_basic_linalg, ONLY: cp_fm_trace
23  cp_fm_struct_type
24  USE cp_fm_types, ONLY: cp_fm_create,&
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
37  equidistant_nodes_a_b,&
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, &
50 
51  INTEGER, PARAMETER, PUBLIC :: cc_shape_linear = contour_shape_linear, &
53 
54  PUBLIC :: ccquad_type
55 
56  PUBLIC :: ccquad_init, &
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  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  POINTER :: w_data, w_data_my
128  TYPE(cp_fm_struct_type), POINTER :: fm_struct
129 
130  CALL timeset(routinen, handle)
131 
132  cpassert(nnodes > 2)
133 
134  ! ensure that MOD(nnodes-1, 2) == 0
135  nnodes = 2*((nnodes - 1)/2) + 1
136 
137  cc_env%interval_id = interval_id
138  cc_env%shape_id = shape_id
139  cc_env%a = a
140  cc_env%b = b
141  cc_env%error = huge(0.0_dp)
142 
143  NULLIFY (cc_env%integral, cc_env%error_fm, cc_env%weights)
144  ALLOCATE (cc_env%weights)
145  CALL cp_fm_get_info(weights, local_data=w_data, nrow_local=nrows, ncol_local=ncols, matrix_struct=fm_struct)
146  CALL cp_fm_create(cc_env%weights, fm_struct)
147  CALL cp_fm_get_info(cc_env%weights, local_data=w_data_my)
148 
149  ! use the explicit loop to avoid temporary arrays
150  DO icol = 1, ncols
151  DO irow = 1, nrows
152  w_data_my(irow, icol) = abs(w_data(irow, icol))
153  END DO
154  END DO
155 
156  SELECT CASE (interval_id)
157  CASE (cc_interval_full)
158  nnodes_half = nnodes/2 + 1
159  CASE (cc_interval_half)
160  nnodes_half = nnodes
161  CASE DEFAULT
162  cpabort("Unimplemented interval type")
163  END SELECT
164 
165  ALLOCATE (cc_env%tnodes(nnodes))
166 
167  IF (PRESENT(tnodes_restart)) THEN
168  cc_env%tnodes(1:nnodes) = tnodes_restart(1:nnodes)
169  ELSE
170  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  IF (nnodes_half > 2) &
176  CALL rescale_nodes_cos(nnodes_half - 2, cc_env%tnodes(2:))
177 
178  SELECT CASE (interval_id)
179  CASE (cc_interval_full)
180  ! reflect symmetric nodes
181  DO ipoint = nnodes_half - 1, 1, -1
182  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  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  CALL rescale_normalised_nodes(nnodes, cc_env%tnodes, a, b, shape_id, xnodes)
191 
192  CALL timestop(handle)
193  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  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  CALL timeset(routinen, handle)
209 
210  IF (ASSOCIATED(cc_env%error_fm)) THEN
211  CALL cp_fm_release(cc_env%error_fm)
212  DEALLOCATE (cc_env%error_fm)
213  NULLIFY (cc_env%error_fm)
214  END IF
215 
216  IF (ASSOCIATED(cc_env%weights)) THEN
217  CALL cp_fm_release(cc_env%weights)
218  DEALLOCATE (cc_env%weights)
219  NULLIFY (cc_env%weights)
220  END IF
221 
222  IF (ASSOCIATED(cc_env%integral)) THEN
223  CALL cp_cfm_release(cc_env%integral)
224  DEALLOCATE (cc_env%integral)
225  NULLIFY (cc_env%integral)
226  END IF
227 
228  IF (ALLOCATED(cc_env%zdata_cache)) THEN
229  DO ipoint = SIZE(cc_env%zdata_cache), 1, -1
230  CALL cp_cfm_release(cc_env%zdata_cache(ipoint))
231  END DO
232 
233  DEALLOCATE (cc_env%zdata_cache)
234  END IF
235 
236  IF (ALLOCATED(cc_env%tnodes)) DEALLOCATE (cc_env%tnodes)
237 
238  CALL timestop(handle)
239  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tnodes, tnodes_old
259 
260  CALL timeset(routinen, handle)
261 
262  cpassert(.NOT. ALLOCATED(xnodes_next))
263  cpassert(ASSOCIATED(cc_env%integral))
264  cpassert(ASSOCIATED(cc_env%error_fm))
265  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  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  nnodes_half = nnodes_exist - 1
272 
273  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  nnodes_next = 2*nnodes_half
277  CASE (cc_interval_half)
278  nnodes_next = nnodes_half
279  CASE DEFAULT
280  cpabort("Unimplemented interval type")
281  END SELECT
282 
283  ALLOCATE (xnodes_next(nnodes_next))
284  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  nnodes_half, tnodes)
289 
290  CALL rescale_nodes_cos(nnodes_half, tnodes)
291 
292  SELECT CASE (cc_env%interval_id)
293  CASE (cc_interval_full)
294  ! reflect symmetric nodes
295  DO ipoint = 1, nnodes_half
296  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  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  CALL move_alloc(cc_env%tnodes, tnodes_old)
305  nnodes_exist = SIZE(tnodes_old)
306 
307  ALLOCATE (cc_env%tnodes(nnodes_exist + nnodes_next))
308  cc_env%tnodes(1:nnodes_exist) = tnodes_old(1:nnodes_exist)
309  cc_env%tnodes(nnodes_exist + 1:nnodes_exist + nnodes_next) = tnodes(1:nnodes_next)
310  DEALLOCATE (tnodes_old)
311 
312  ! rescale nodes [-1 .. 1] -> [a .. b] according to the shape
313  CALL rescale_normalised_nodes(nnodes_next, tnodes, cc_env%a, cc_env%b, cc_env%shape_id, xnodes_next)
314 
315  DEALLOCATE (tnodes)
316  CALL timestop(handle)
317  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  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  COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: zscale
340  INTEGER :: handle, ipoint, nnodes_exist, &
341  nnodes_half, nnodes_next
342  TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: zdata_tmp
343 
344  CALL timeset(routinen, handle)
345 
346  nnodes_next = SIZE(zdata_next)
347  cpassert(nnodes_next > 0)
348 
349  ! compute weights of new points on a complex contour according to their values of the 't' parameter
350  nnodes_exist = SIZE(cc_env%tnodes)
351  cpassert(nnodes_exist >= nnodes_next)
352 
353  ALLOCATE (zscale(nnodes_next))
354  CALL rescale_normalised_nodes(nnodes_next, cc_env%tnodes(nnodes_exist - nnodes_next + 1:nnodes_exist), &
355  cc_env%a, cc_env%b, cc_env%shape_id, weights=zscale)
356 
357  IF (cc_env%interval_id == cc_interval_half) zscale(:) = 2.0_dp*zscale(:)
358 
359  ! rescale integrand values
360  DO ipoint = 1, nnodes_next
361  CALL cp_cfm_scale(zscale(ipoint), zdata_next(ipoint))
362  END DO
363  DEALLOCATE (zscale)
364 
365  ! squash points with the same clenshaw-curtis weights together
366  IF (ALLOCATED(cc_env%zdata_cache)) THEN
367  nnodes_exist = SIZE(cc_env%zdata_cache)
368  ELSE
369  nnodes_exist = 0
370  END IF
371 
372  SELECT CASE (cc_env%interval_id)
373  CASE (cc_interval_full)
374  IF (ALLOCATED(cc_env%zdata_cache)) THEN
375  cpassert(nnodes_exist == nnodes_next/2 + 1)
376  nnodes_half = nnodes_exist - 1
377  ELSE
378  cpassert(mod(nnodes_next, 2) == 1)
379  nnodes_half = nnodes_next/2 + 1
380  END IF
381  CASE (cc_interval_half)
382  IF (ALLOCATED(cc_env%zdata_cache)) THEN
383  cpassert(nnodes_exist == nnodes_next + 1)
384  END IF
385 
386  nnodes_half = nnodes_next
387  END SELECT
388 
389  IF (cc_env%interval_id == cc_interval_full) THEN
390  DO ipoint = nnodes_next/2, 1, -1
391  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  IF (ALLOCATED(cc_env%zdata_cache)) THEN
396  ! note that nnodes_half+1 == nnodes_exist for both half- and full-intervals
397  ALLOCATE (zdata_tmp(nnodes_half + nnodes_exist))
398 
399  DO ipoint = 1, nnodes_half
400  zdata_tmp(2*ipoint - 1) = cc_env%zdata_cache(ipoint)
401  zdata_tmp(2*ipoint) = zdata_next(ipoint)
402  zdata_next(ipoint) = cfm_null
403  END DO
404  zdata_tmp(nnodes_half + nnodes_exist) = cc_env%zdata_cache(nnodes_exist)
405 
406  CALL move_alloc(zdata_tmp, cc_env%zdata_cache)
407  ELSE
408  CALL cp_cfm_scale(2.0_dp, zdata_next(nnodes_half))
409 
410  ALLOCATE (cc_env%zdata_cache(nnodes_half))
411 
412  DO ipoint = 1, nnodes_half
413  cc_env%zdata_cache(ipoint) = zdata_next(ipoint)
414  zdata_next(ipoint) = cfm_null
415  END DO
416  END IF
417 
418  CALL timestop(handle)
419  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  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  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  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  CALL timeset(routinen, handle)
445 
446  cpassert(ALLOCATED(cc_env%zdata_cache))
447 
448  nintervals_half_plus_1 = SIZE(cc_env%zdata_cache)
449  nintervals_half = nintervals_half_plus_1 - 1
450  nintervals_half_plus_2 = nintervals_half_plus_1 + 1
451  nintervals = 2*nintervals_half
452  nintervals_plus_2 = nintervals + 2
453  cpassert(nintervals_half > 1)
454 
455  IF (.NOT. ASSOCIATED(cc_env%integral)) THEN
456  CALL cp_cfm_get_info(cc_env%zdata_cache(1), matrix_struct=fm_struct)
457  equiv = cp_fm_struct_equivalent(fm_struct, cc_env%weights%matrix_struct)
458  cpassert(equiv)
459 
460  ALLOCATE (cc_env%integral)
461  CALL cp_cfm_create(cc_env%integral, fm_struct)
462  NULLIFY (cc_env%error_fm)
463  ALLOCATE (cc_env%error_fm)
464  CALL cp_fm_create(cc_env%error_fm, fm_struct)
465  END IF
466 
467  IF (debug_this_module) THEN
468  DO ipoint = 1, nintervals_half_plus_1
469  equiv = cp_fm_struct_equivalent(cc_env%zdata_cache(ipoint)%matrix_struct, cc_env%integral%matrix_struct)
470  cpassert(equiv)
471  END DO
472  END IF
473 
474  CALL cp_cfm_get_info(cc_env%integral, nrow_local=nrows_local, ncol_local=ncols_local)
475 
476  ALLOCATE (weights(nintervals_half))
477 
478  ! omit the trivial weights(1) = 0.5
479  DO ipoint = 2, nintervals_half
480  rscale = real(2*(ipoint - 1), kind=dp)
481  weights(ipoint) = 1.0_dp/(1.0_dp - rscale*rscale)
482  END DO
483  ! weights(1) <- weights(intervals_half + 1)
484  rscale = real(nintervals, kind=dp)
485  weights(1) = 1.0_dp/(1.0_dp - rscale*rscale)
486 
487  ! 1.0 / nintervals
488  rscale = 1.0_dp/rscale
489 
490  CALL fft_alloc(ztmp, [nintervals, nrows_local, ncols_local])
491  CALL fft_alloc(ztmp_dct, [nintervals, nrows_local, ncols_local])
492 
493 !$OMP PARALLEL DO DEFAULT(NONE), PRIVATE(icol, ipoint, irow), &
494 !$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  CALL fft_fw1d(nintervals, nrows_local*ncols_local, .false., ztmp, ztmp_dct, 1.0_dp, stat)
509  IF (stat /= 0) THEN
510  CALL cp_abort(__location__, &
511  "An FFT library is required for Clenshaw-Curtis quadrature. "// &
512  "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 !$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  CALL fft_dealloc(ztmp)
534  CALL fft_dealloc(ztmp_dct)
535 
536  CALL cp_fm_trace(cc_env%error_fm, cc_env%weights, cc_env%error)
537 
538  DEALLOCATE (weights)
539  CALL timestop(handle)
540  END SUBROUTINE ccquad_refine_integral
541 
542 END MODULE negf_integr_cc
Basic linear algebra operations for complex full matrices.
subroutine, public cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b).
Represents a complex full matrix distributed on many processors.
Definition: cp_cfm_types.F:12
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
Definition: cp_cfm_types.F:121
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition: cp_cfm_types.F:159
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition: cp_cfm_types.F:607
basic linear algebra operations for full matrices
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
logical function, public cp_fm_struct_equivalent(fmstruct1, fmstruct2)
returns true if the two matrix structures are equivalent, false otherwise.
Definition: cp_fm_struct.F:357
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
subroutine, public fft_fw1d(n, m, trans, zin, zout, scale, stat)
Performs m 1-D forward FFT-s of size n.
Definition: fft_tools.F:331
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition: kahan_sum.F:29
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public int_8
Definition: kinds.F:54
integer, parameter, public dp
Definition: kinds.F:34
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
complex(kind=dp), parameter, public z_one
Adaptive Clenshaw-Curtis quadrature algorithm to integrate a complex-valued function in a complex pla...
integer, parameter, public cc_shape_linear
subroutine, public ccquad_refine_integral(cc_env)
Refine approximated integral.
integer, parameter, public cc_interval_full
subroutine, public ccquad_double_number_of_points(cc_env, xnodes_next)
Get the next set of points at which the integrand needs to be computed. These points are then can be ...
subroutine, public ccquad_release(cc_env)
Release a Clenshaw-Curtis quadrature environment variable.
integer, parameter, public cc_interval_half
subroutine, public ccquad_init(cc_env, xnodes, nnodes, a, b, interval_id, shape_id, weights, tnodes_restart)
Initialise a Clenshaw-Curtis quadrature environment variable.
integer, parameter, public cc_shape_arc
subroutine, public ccquad_reduce_and_append_zdata(cc_env, zdata_next)
Prepare Clenshaw-Curtis environment for the subsequent refinement of the integral.
Helper functions for integration routines.
subroutine, public rescale_normalised_nodes(nnodes, tnodes, a, b, shape_id, xnodes, weights)
subroutine, public rescale_nodes_cos(nnodes, tnodes)
Rescale nodes tnodes(i) = cos(pi/2 * (1-tnodes(i))); tnodes \in [-1 .. 1] .
integer, parameter, public contour_shape_arc
integer, parameter, public contour_shape_linear