minimax_exp.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Routines to calculate the minimax coefficients in order to
!>        approximate 1/x as a sum over exponential functions
!>        1/x ~ SUM_{i}^{K} w_i EXP(-a_i * x) for x belonging to [1:Rc].
!>
!>        This module is an extension of original minimax module minimax_exp_k15
!>        (up to K = 15) to provide minimax approximations for larger
!>        ranges Rc (up to K = 53).
!>
!>        k53 implementation is based on directly tabulated coefficients from
!>        D. Braess and W. Hackbusch, IMA Journal of Numerical Analysis 25.4 (2005): 685-697
!>        http://www.mis.mpg.de/scicomp/EXP_SUM/1_x
!>
!>        Note: Due to discrete Rc values, the k53 implementation does not yield
!>        optimal approximations for arbitrary Rc. If optimal minimax coefficients
!>        are needed, the minimax_exp_k15 module should be extended by interpolating
!>        k53 coefficients.
!> \par History
!>      02.2016 created [Patrick Seewald]
! **************************************************************************************************
 
MODULE minimax_exp
   USE cp_log_handling,                 ONLY: cp_to_string
   USE kinds,                           ONLY: dp
   USE minimax_exp_k15,                 ONLY: check_range_k15,&
                                              get_minimax_coeff_k15,&
                                              get_minimax_numerical_error
   USE minimax_exp_k53,                 ONLY: r_max,&
                                              r_mm,&
                                              err_mm,&
                                              get_minimax_coeff_low,&
                                              k_max,&
                                              k_mm,&
                                              k_p,&
                                              n_approx
#include "../base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'minimax_exp'
 
   INTEGER, PARAMETER :: mm_k15 = 0, mm_k53 = 1
 
   PUBLIC :: get_exp_minimax_coeff, validate_exp_minimax, check_exp_minimax_range
 
   ! Imported from minimax_k53:
 
   ! Number of tabulated minimax approximations:
   ! INTEGER, PARAMETER :: n_approx
 
   ! Number of K values:
   ! INTEGER, PARAMETER :: n_k
 
   ! Maximum K value:
   ! INTEGER, PARAMETER :: k_max
 
   ! Maximum range Rc:
   ! REAL(KIND=dp), PARAMETER :: R_max
 
   ! Values of K:
   ! INTEGER, PARAMETER, DIMENSION(n_approx) :: k_mm
 
   ! Values of Rc:
   ! REAL(KIND=dp), PARAMETER, DIMENSION(n_approx) :: R_mm
 
   ! Values of minimax error:
   ! REAL(KIND=dp), PARAMETER, DIMENSION(n_approx) :: err_mm
 
   ! Note: the coefficients (k_mm, R_mm, err_mm) are sorted w.r.t. 1) k_mm, 2) R_mm
 
   ! Given the ith value of K, k_p(i) points to the first minimax
   ! approximation with K terms:
   ! INTEGER, PARAMETER, DIMENSION(n_k+1) :: k_p
 
   ! Minimax coefficients aw of the ith minimax approximation:
   ! SUBROUTINE get_minimax_coeff_low(i, aw)
 
CONTAINS
 
! **************************************************************************************************
!> \brief Check that a minimax approximation is available for given input k, Rc.
!>        ierr ==  0: everything ok
!>        ierr ==  1: Rc too small
!>        ierr == -1: k too large
!> \param k ...
!> \param Rc ...
!> \param ierr ...
!> \note: ierr ==  1 is not a fatal error since get_exp_minimax_coeff will return
!>        k53 minimax coefficients with smallest possible range.
! **************************************************************************************************
   SUBROUTINE check_exp_minimax_range(k, Rc, ierr)
      INTEGER, INTENT(IN)                                :: k
      REAL(kind=dp), INTENT(IN)                          :: rc
      INTEGER, INTENT(OUT)                               :: ierr
 
      ierr = 0
      IF (k .LE. 15) THEN
         CALL check_range_k15(k, rc, ierr)
      ELSE
         IF (k .GT. k_max) ierr = -1
      END IF
 
   END SUBROUTINE check_exp_minimax_range
 
! **************************************************************************************************
!> \brief Get best minimax approximation for given input parameters. Automatically
!>        chooses the most exact set of minimax coefficients (k15 or k53) for
!>        given k, Rc.
!> \param k Number of minimax terms
!> \param Rc Minimax range
!> \param aw The a_i and w_i coefficient are stored in aw such that the first 1:K
!>        elements correspond to a_i and the K+1:2k correspond to w_i.
!> \param mm_error Numerical error of minimax approximation for given k, Rc
!> \param which_coeffs Whether the coefficients returned have been generated from
!>        k15 or k53 coefficients (mm_k15 or mm_k53).
! **************************************************************************************************
   SUBROUTINE get_exp_minimax_coeff(k, Rc, aw, mm_error, which_coeffs)
      INTEGER, INTENT(IN)                                :: k
      REAL(kind=dp), INTENT(IN)                          :: rc
      REAL(kind=dp), DIMENSION(2*k), INTENT(OUT)         :: aw
      REAL(kind=dp), INTENT(OUT), OPTIONAL               :: mm_error
      INTEGER, INTENT(OUT), OPTIONAL                     :: which_coeffs
 
      INTEGER                                            :: ierr
 
      IF (k .LE. 15) THEN
         CALL check_range_k15(k, rc, ierr)
         IF (ierr .EQ. 1) THEN ! Rc too small for k15 coeffs --> use k53
            CALL get_minimax_coeff_k53(k, rc, aw, mm_error)
            IF (PRESENT(which_coeffs)) which_coeffs = mm_k53
         ELSE
            cpassert(ierr .EQ. 0)
            CALL get_minimax_coeff_k15(k, rc, aw, mm_error)
            IF (PRESENT(which_coeffs)) which_coeffs = mm_k15
         END IF
      ELSEIF (k .LE. 53) THEN
         CALL get_minimax_coeff_k53(k, rc, aw, mm_error)
         IF (PRESENT(which_coeffs)) which_coeffs = mm_k53
      ELSE
         cpabort("No minimax approximations available for k = "//cp_to_string(k))
      END IF
   END SUBROUTINE get_exp_minimax_coeff
 
! **************************************************************************************************
!> \brief Get minimax coefficients: k53 implementation (coefficients up to k=53 terms).
!>        All a_i and w_i for a set of discrete values Rc, k are tabulated and
!>        the most accurate coefficients for given input k, Rc are returned.
!> \param k ...
!> \param Rc ...
!> \param aw ...
!> \param mm_error ...
! **************************************************************************************************
   SUBROUTINE get_minimax_coeff_k53(k, Rc, aw, mm_error)
      INTEGER, INTENT(IN)                                :: k
      REAL(kind=dp), INTENT(IN)                          :: rc
      REAL(kind=dp), DIMENSION(2*k), INTENT(OUT)         :: aw
      REAL(kind=dp), INTENT(OUT), OPTIONAL               :: mm_error
 
      INTEGER                                            :: i_mm
 
      CALL get_best_minimax_approx_k53(k, rc, i_mm)
      CALL get_minimax_coeff_low(i_mm, aw)
      IF (PRESENT(mm_error)) mm_error = get_minimax_numerical_error(rc, aw)
 
   END SUBROUTINE get_minimax_coeff_k53
 
! **************************************************************************************************
!> \brief find minimax approx. with k terms that is most accurate for range Rc.
!> \param k ...
!> \param Rc ...
!> \param i_mm ...
!> \param ge_Rc Whether the tabulated range of the returned minimax approximation
!>              must be strictly greater than or equal to Rc. Default is .FALSE.
! **************************************************************************************************
   SUBROUTINE get_best_minimax_approx_k53(k, Rc, i_mm, ge_Rc)
      INTEGER, INTENT(IN)                                :: k
      REAL(kind=dp), INTENT(IN)                          :: rc
      INTEGER, INTENT(OUT)                               :: i_mm
      LOGICAL, INTENT(IN), OPTIONAL                      :: ge_rc
 
      INTEGER                                            :: i, i_k, i_l, i_r
      REAL(kind=dp)                                      :: error_l, error_r, r_k_max, r_k_min
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: aw
 
      cpassert(k .LE. k_max)
 
      ! find k pointer and smallest and largest R_mm value for this k
      i_k = 1
      DO WHILE (k_mm(k_p(i_k)) .LT. k)
         i_k = i_k + 1
      END DO
      cpassert(k_mm(k_p(i_k)) .EQ. k)
 
      r_k_min = r_mm(k_p(i_k))
      r_k_max = r_mm(k_p(i_k + 1) - 1)
 
      IF (rc .GE. r_k_max) THEN
         i_mm = k_p(i_k + 1) - 1 ! pointer to largest Rc for current k
      ELSE IF (rc .LE. r_k_min) THEN
         i_mm = k_p(i_k) ! pointer to smallest Rc for current k
      ELSE
         i = k_p(i_k)
         DO WHILE (rc .GT. r_mm(i))
            i = i + 1
         END DO
         i_r = i ! pointer to closest R_mm >= Rc
         i_l = i - 1 ! pointer to closest R_mm < Rc
 
         IF (PRESENT(ge_rc)) THEN
            IF (ge_rc) THEN
               i_mm = i_r
               RETURN
            END IF
         END IF
 
         ALLOCATE (aw(2*k_mm(i_r)))
         CALL get_minimax_coeff_low(i_r, aw)
         error_l = get_minimax_numerical_error(rc, aw)
         DEALLOCATE (aw)
         ALLOCATE (aw(2*k_mm(i_l)))
         CALL get_minimax_coeff_low(i_l, aw)
         error_r = get_minimax_numerical_error(rc, aw)
         DEALLOCATE (aw)
         i_mm = merge(i_r, i_l, error_l .LE. error_r)
      END IF
 
   END SUBROUTINE get_best_minimax_approx_k53
 
! **************************************************************************************************
!> \brief Unit test checking that numerical error of minimax approximations
!>        generated using any k15 or k53 coefficients is consistent with
!>        tabulated error.
!> \param n_R Number of Rc values to be tested.
!> \param iw ...
! **************************************************************************************************
   SUBROUTINE validate_exp_minimax(n_R, iw)
      INTEGER, INTENT(IN)                                :: n_r, iw
 
      INTEGER                                            :: i_mm, i_r, ierr, k, which_coeffs
      LOGICAL                                            :: do_exit
      REAL(kind=dp)                                      :: dr, mm_error, r, ref_error
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: aw
 
      IF (iw > 0) THEN
         WRITE (iw, '(//T2,A)') &
            "Unit tests for minimax 1/x ~ SUM_{i}^{K} w_i EXP(-a_i * x) for x belonging to [1:Rc]"
         WRITE (iw, '(T2,84("*"))')
      END IF
 
      IF (iw > 0) THEN
         WRITE (iw, '(/T2,A)') &
            "1) checking numerical error against tabulated error at tabulated values Rc"
         WRITE (iw, '(/T2,A)') &
            "which coeffs, K, Rc, num. error, ref. error, rel. diff (num. error - ref. error)/(ref. error)"
         WRITE (iw, '(T2,54("-"))')
      END IF
      DO i_mm = 1, n_approx
         r = r_mm(i_mm)
         k = k_mm(i_mm)
         CALL check_exp_minimax_range(k, r, ierr)
         IF (ierr .EQ. 0) THEN
            ALLOCATE (aw(2*k))
            CALL get_exp_minimax_coeff(k, r, aw, mm_error, which_coeffs)
            ref_error = err_mm(i_mm)
            DEALLOCATE (aw)
            IF (iw > 0) WRITE (iw, '(T2,A4, I3, ES10.1, ES12.3, ES12.3, ES12.3)') &
               merge("k15", "k53", which_coeffs .EQ. mm_k15), k, r, &
               mm_error, ref_error, (mm_error - ref_error)/ref_error
            cpassert(mm_error .LE. ref_error*1.05_dp + 1.0e-15_dp)
         ELSE
            IF (iw > 0) WRITE (iw, '(T2,A4, I3, ES10.1, 3X, A)') "k15", k, r, "missing"
         END IF
 
         IF (k .LE. 15) THEN
            ALLOCATE (aw(2*k))
            CALL get_minimax_coeff_k53(k, r, aw, mm_error)
            ref_error = err_mm(i_mm)
            DEALLOCATE (aw)
            IF (iw > 0) WRITE (iw, '(T2,A4,I3, ES10.1, ES12.3, ES12.3, ES12.3)') &
               "k53", k, r, mm_error, ref_error, (mm_error - ref_error)/ref_error
            IF (mm_error .GT. ref_error*1.05_dp + 1.0e-15_dp) THEN
               cpabort("Test 1 failed: numerical error is larger than tabulated error")
            END IF
         END IF
      END DO
 
      IF (iw > 0 .AND. n_r .GT. 0) THEN
         WRITE (iw, '(T2,54("-"))')
         WRITE (iw, '(/T2,A)') "Test 1 OK"
         WRITE (iw, '(/T2,A)') &
            "2) checking numerical error against tabulated error at arbitrary values Rc"
         WRITE (iw, '(/T2,A)') &
            "which coeffs, K, Rc, num. error, ref. error, rel. diff (num. error - ref. error)/(ref. error)"
         WRITE (iw, '(T2,54("-"))')
      END IF
      dr = r_max**(1.0_dp/n_r)
 
      DO k = 1, k_max
         ALLOCATE (aw(2*k))
         do_exit = .false.
         DO i_r = 1, n_r
            r = dr**i_r
            CALL get_exp_minimax_coeff(k, r, aw, mm_error, which_coeffs)
            CALL get_best_minimax_approx_k53(k, r, i_mm, ge_rc=.true.)
            IF (r .GT. r_mm(i_mm)) THEN
               r = r_max
               do_exit = .true.
            END IF
            ref_error = err_mm(i_mm)
            IF (iw > 0) WRITE (iw, '(T2, A4, I3, ES10.1, ES12.3, ES12.3, ES12.3)') &
               merge("k15", "k53", which_coeffs .EQ. mm_k15), k, r, &
               mm_error, ref_error, (mm_error - ref_error)/ref_error
            IF (mm_error .GT. ref_error*1.05_dp + 1.0e-15_dp) THEN
               cpabort("Test 2 failed: numerical error is larger than tabulated error")
            END IF
            IF (do_exit) EXIT
         END DO
         DEALLOCATE (aw)
      END DO
      IF (iw > 0) THEN
         WRITE (iw, '(T2,54("-"))')
         WRITE (iw, '(/T2,A)') "Test 2 OK"
      END IF
   END SUBROUTINE validate_exp_minimax
 
END MODULE minimax_exp