37 equidistant_nodes_a_b,&
40 #include "./base/base_uses.f90"
45 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'negf_integr_cc'
46 LOGICAL,
PARAMETER,
PRIVATE :: debug_this_module = .true.
67 COMPLEX(kind=dp) :: a, b
73 INTEGER :: interval_id
77 REAL(kind=
dp) :: error
79 TYPE(cp_cfm_type),
POINTER :: integral
81 TYPE(cp_fm_type),
POINTER :: error_fm
83 TYPE(cp_fm_type),
POINTER :: weights
86 TYPE(cp_cfm_type),
ALLOCATABLE,
DIMENSION(:) :: zdata_cache
87 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: tnodes
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
122 CHARACTER(len=*),
PARAMETER :: routinen =
'ccquad_init'
124 INTEGER :: handle, icol, ipoint, irow, ncols, &
126 REAL(kind=
dp),
CONTIGUOUS,
DIMENSION(:, :), &
127 POINTER :: w_data, w_data_my
128 TYPE(cp_fm_struct_type),
POINTER :: fm_struct
130 CALL timeset(routinen, handle)
135 nnodes = 2*((nnodes - 1)/2) + 1
137 cc_env%interval_id = interval_id
138 cc_env%shape_id = shape_id
141 cc_env%error = huge(0.0_dp)
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)
152 w_data_my(irow, icol) = abs(w_data(irow, icol))
156 SELECT CASE (interval_id)
158 nnodes_half = nnodes/2 + 1
162 cpabort(
"Unimplemented interval type")
165 ALLOCATE (cc_env%tnodes(nnodes))
167 IF (
PRESENT(tnodes_restart))
THEN
168 cc_env%tnodes(1:nnodes) = tnodes_restart(1:nnodes)
170 CALL equidistant_nodes_a_b(-1.0_dp, 0.0_dp, nnodes_half, cc_env%tnodes)
175 IF (nnodes_half > 2) &
178 SELECT CASE (interval_id)
181 DO ipoint = nnodes_half - 1, 1, -1
182 cc_env%tnodes(nnodes_half + ipoint) = -cc_env%tnodes(nnodes_half - ipoint)
186 cc_env%tnodes(1:nnodes_half) = 2.0_dp*cc_env%tnodes(1:nnodes_half) + 1.0_dp
192 CALL timestop(handle)
202 TYPE(ccquad_type),
INTENT(inout) :: cc_env
204 CHARACTER(len=*),
PARAMETER :: routinen =
'ccquad_release'
206 INTEGER :: handle, ipoint
208 CALL timeset(routinen, handle)
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)
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)
222 IF (
ASSOCIATED(cc_env%integral))
THEN
224 DEALLOCATE (cc_env%integral)
225 NULLIFY (cc_env%integral)
228 IF (
ALLOCATED(cc_env%zdata_cache))
THEN
229 DO ipoint =
SIZE(cc_env%zdata_cache), 1, -1
233 DEALLOCATE (cc_env%zdata_cache)
236 IF (
ALLOCATED(cc_env%tnodes))
DEALLOCATE (cc_env%tnodes)
238 CALL timestop(handle)
250 TYPE(ccquad_type),
INTENT(inout) :: cc_env
251 COMPLEX(kind=dp),
ALLOCATABLE,
DIMENSION(:), &
252 INTENT(inout) :: xnodes_next
254 CHARACTER(len=*),
PARAMETER :: routinen =
'ccquad_double_number_of_points'
256 INTEGER :: handle, ipoint, nnodes_exist, &
257 nnodes_half, nnodes_next
258 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: tnodes, tnodes_old
260 CALL timeset(routinen, handle)
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))
268 nnodes_exist =
SIZE(cc_env%zdata_cache)
271 nnodes_half = nnodes_exist - 1
273 SELECT CASE (cc_env%interval_id)
276 nnodes_next = 2*nnodes_half
278 nnodes_next = nnodes_half
280 cpabort(
"Unimplemented interval type")
283 ALLOCATE (xnodes_next(nnodes_next))
284 ALLOCATE (tnodes(nnodes_next))
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), &
292 SELECT CASE (cc_env%interval_id)
295 DO ipoint = 1, nnodes_half
296 tnodes(nnodes_half + ipoint) = -tnodes(nnodes_half - ipoint + 1)
300 tnodes(1:nnodes_half) = 2.0_dp*tnodes(1:nnodes_half) + 1.0_dp
304 CALL move_alloc(cc_env%tnodes, tnodes_old)
305 nnodes_exist =
SIZE(tnodes_old)
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)
316 CALL timestop(handle)
333 TYPE(ccquad_type),
INTENT(inout) :: cc_env
334 TYPE(cp_cfm_type),
DIMENSION(:),
INTENT(inout) :: zdata_next
336 CHARACTER(len=*),
PARAMETER :: routinen =
'ccquad_reduce_and_append_zdata'
337 TYPE(cp_cfm_type),
PARAMETER :: cfm_null = cp_cfm_type()
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
344 CALL timeset(routinen, handle)
346 nnodes_next =
SIZE(zdata_next)
347 cpassert(nnodes_next > 0)
350 nnodes_exist =
SIZE(cc_env%tnodes)
351 cpassert(nnodes_exist >= nnodes_next)
353 ALLOCATE (zscale(nnodes_next))
355 cc_env%a, cc_env%b, cc_env%shape_id, weights=zscale)
360 DO ipoint = 1, nnodes_next
361 CALL cp_cfm_scale(zscale(ipoint), zdata_next(ipoint))
366 IF (
ALLOCATED(cc_env%zdata_cache))
THEN
367 nnodes_exist =
SIZE(cc_env%zdata_cache)
372 SELECT CASE (cc_env%interval_id)
374 IF (
ALLOCATED(cc_env%zdata_cache))
THEN
375 cpassert(nnodes_exist == nnodes_next/2 + 1)
376 nnodes_half = nnodes_exist - 1
378 cpassert(mod(nnodes_next, 2) == 1)
379 nnodes_half = nnodes_next/2 + 1
382 IF (
ALLOCATED(cc_env%zdata_cache))
THEN
383 cpassert(nnodes_exist == nnodes_next + 1)
386 nnodes_half = nnodes_next
390 DO ipoint = nnodes_next/2, 1, -1
395 IF (
ALLOCATED(cc_env%zdata_cache))
THEN
397 ALLOCATE (zdata_tmp(nnodes_half + nnodes_exist))
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
404 zdata_tmp(nnodes_half + nnodes_exist) = cc_env%zdata_cache(nnodes_exist)
406 CALL move_alloc(zdata_tmp, cc_env%zdata_cache)
408 CALL cp_cfm_scale(2.0_dp, zdata_next(nnodes_half))
410 ALLOCATE (cc_env%zdata_cache(nnodes_half))
412 DO ipoint = 1, nnodes_half
413 cc_env%zdata_cache(ipoint) = zdata_next(ipoint)
414 zdata_next(ipoint) = cfm_null
418 CALL timestop(handle)
428 TYPE(ccquad_type),
INTENT(inout) :: cc_env
430 CHARACTER(len=*),
PARAMETER :: routinen =
'ccquad_refine_integral'
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
437 REAL(kind=
dp) :: rscale
438 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: weights
439 TYPE(cp_fm_struct_type),
POINTER :: fm_struct
444 CALL timeset(routinen, handle)
446 cpassert(
ALLOCATED(cc_env%zdata_cache))
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)
455 IF (.NOT.
ASSOCIATED(cc_env%integral))
THEN
460 ALLOCATE (cc_env%integral)
462 NULLIFY (cc_env%error_fm)
463 ALLOCATE (cc_env%error_fm)
467 IF (debug_this_module)
THEN
468 DO ipoint = 1, nintervals_half_plus_1
474 CALL cp_cfm_get_info(cc_env%integral, nrow_local=nrows_local, ncol_local=ncols_local)
476 ALLOCATE (weights(nintervals_half))
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)
484 rscale = real(nintervals, kind=
dp)
485 weights(1) = 1.0_dp/(1.0_dp - rscale*rscale)
488 rscale = 1.0_dp/rscale
490 CALL fft_alloc(ztmp, [nintervals, nrows_local, ncols_local])
491 CALL fft_alloc(ztmp_dct, [nintervals, nrows_local, ncols_local])
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)
501 DO ipoint = 2, nintervals_half
502 ztmp(nintervals_half + ipoint, irow, icol) = ztmp(nintervals_half_plus_2 - ipoint, irow, icol)
508 CALL fft_fw1d(nintervals, nrows_local*ncols_local, .false., ztmp, ztmp_dct, 1.0_dp, stat)
510 CALL cp_abort(__location__, &
511 "An FFT library is required for Clenshaw-Curtis quadrature. "// &
512 "You can use an alternative integration method instead.")
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))
525 ztmp_dct(nintervals_half_plus_1, irow, icol) = weights(1)*ztmp_dct(nintervals_half_plus_1, irow, icol)
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))
533 CALL fft_dealloc(ztmp)
534 CALL fft_dealloc(ztmp_dct)
536 CALL cp_fm_trace(cc_env%error_fm, cc_env%weights, cc_env%error)
539 CALL timestop(handle)
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.
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
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.
basic linear algebra operations for full matrices
represent the structure of a full matrix
logical function, public cp_fm_struct_equivalent(fmstruct1, fmstruct2)
returns true if the two matrix structures are equivalent, false otherwise.
represent a full matrix distributed on many processors
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
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Defines the basic variable types.
integer, parameter, public int_8
integer, parameter, public dp
Definition of mathematical constants and functions.
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