LCOV - code coverage report
Current view: top level - src - negf_integr_cc.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 169 171 98.8 %
Date: 2024-03-28 07:31:50 Functions: 5 7 71.4 %

          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

Generated by: LCOV version 1.15