1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines to calculate the minimax coefficients in order to 8!> approximate 1/x as a sum over exponential functions 9!> 1/x ~ SUM_{i}^{K} w_i EXP(-a_i * x) for x belonging to [1:Rc]. 10!> 11!> This module is an extension of original minimax module minimax_exp_k15 12!> (up to K = 15) to provide minimax approximations for larger 13!> ranges Rc (up to K = 53). 14!> 15!> k53 implementation is based on directly tabulated coefficients from 16!> D. Braess and W. Hackbusch, IMA Journal of Numerical Analysis 25.4 (2005): 685-697 17!> http://www.mis.mpg.de/scicomp/EXP_SUM/1_x 18!> 19!> Note: Due to discrete Rc values, the k53 implementation does not yield 20!> optimal approximations for arbitrary Rc. If optimal minimax coefficients 21!> are needed, the minimax_exp_k15 module should be extended by interpolating 22!> k53 coefficients. 23!> \par History 24!> 02.2016 created [Patrick Seewald] 25! ************************************************************************************************** 26 27MODULE minimax_exp 28 USE cp_log_handling, ONLY: cp_to_string 29 USE kinds, ONLY: dp 30 USE minimax_exp_k15, ONLY: check_range_k15,& 31 get_minimax_coeff_k15,& 32 get_minimax_numerical_error 33 USE minimax_exp_k53, ONLY: R_max,& 34 R_mm,& 35 err_mm,& 36 get_minimax_coeff_low,& 37 k_max,& 38 k_mm,& 39 k_p,& 40 n_approx 41#include "../base/base_uses.f90" 42 43 IMPLICIT NONE 44 45 PRIVATE 46 47 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'minimax_exp' 48 49 INTEGER, PARAMETER :: mm_k15 = 0, mm_k53 = 1 50 51 PUBLIC :: get_exp_minimax_coeff, validate_exp_minimax, check_exp_minimax_range 52 53 ! Imported from minimax_k53: 54 55 ! Number of tabulated minimax approximations: 56 ! INTEGER, PARAMETER :: n_approx 57 58 ! Number of K values: 59 ! INTEGER, PARAMETER :: n_k 60 61 ! Maximum K value: 62 ! INTEGER, PARAMETER :: k_max 63 64 ! Maximum range Rc: 65 ! REAL(KIND=dp), PARAMETER :: R_max 66 67 ! Values of K: 68 ! INTEGER, PARAMETER, DIMENSION(n_approx) :: k_mm 69 70 ! Values of Rc: 71 ! REAL(KIND=dp), PARAMETER, DIMENSION(n_approx) :: R_mm 72 73 ! Values of minimax error: 74 ! REAL(KIND=dp), PARAMETER, DIMENSION(n_approx) :: err_mm 75 76 ! Note: the coefficients (k_mm, R_mm, err_mm) are sorted w.r.t. 1) k_mm, 2) R_mm 77 78 ! Given the ith value of K, k_p(i) points to the first minimax 79 ! approximation with K terms: 80 ! INTEGER, PARAMETER, DIMENSION(n_k+1) :: k_p 81 82 ! Minimax coefficients aw of the ith minimax approximation: 83 ! SUBROUTINE get_minimax_coeff_low(i, aw) 84 85CONTAINS 86 87! ************************************************************************************************** 88!> \brief Check that a minimax approximation is available for given input k, Rc. 89!> ierr == 0: everything ok 90!> ierr == 1: Rc too small 91!> ierr == -1: k too large 92!> \param k ... 93!> \param Rc ... 94!> \param ierr ... 95!> \note: ierr == 1 is not a fatal error since get_exp_minimax_coeff will return 96!> k53 minimax coefficients with smallest possible range. 97! ************************************************************************************************** 98 SUBROUTINE check_exp_minimax_range(k, Rc, ierr) 99 INTEGER, INTENT(IN) :: k 100 REAL(KIND=dp), INTENT(IN) :: Rc 101 INTEGER, INTENT(OUT) :: ierr 102 103 ierr = 0 104 IF (k .LE. 15) THEN 105 CALL check_range_k15(k, Rc, ierr) 106 ELSE 107 IF (k .GT. k_max) ierr = -1 108 ENDIF 109 110 END SUBROUTINE check_exp_minimax_range 111 112! ************************************************************************************************** 113!> \brief Get best minimax approximation for given input parameters. Automatically 114!> chooses the most exact set of minimax coefficients (k15 or k53) for 115!> given k, Rc. 116!> \param k Number of minimax terms 117!> \param Rc Minimax range 118!> \param aw The a_i and w_i coefficient are stored in aw such that the first 1:K 119!> elements correspond to a_i and the K+1:2k correspond to w_i. 120!> \param mm_error Numerical error of minimax approximation for given k, Rc 121!> \param which_coeffs Whether the coefficients returned have been generated from 122!> k15 or k53 coefficients (mm_k15 or mm_k53). 123! ************************************************************************************************** 124 SUBROUTINE get_exp_minimax_coeff(k, Rc, aw, mm_error, which_coeffs) 125 INTEGER, INTENT(IN) :: k 126 REAL(KIND=dp), INTENT(IN) :: Rc 127 REAL(KIND=dp), DIMENSION(2*k), INTENT(OUT) :: aw 128 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: mm_error 129 INTEGER, INTENT(OUT), OPTIONAL :: which_coeffs 130 131 INTEGER :: ierr 132 133 IF (k .LE. 15) THEN 134 CALL check_range_k15(k, Rc, ierr) 135 IF (ierr .EQ. 1) THEN ! Rc too small for k15 coeffs --> use k53 136 CALL get_minimax_coeff_k53(k, Rc, aw, mm_error) 137 IF (PRESENT(which_coeffs)) which_coeffs = mm_k53 138 ELSE 139 CPASSERT(ierr .EQ. 0) 140 CALL get_minimax_coeff_k15(k, Rc, aw, mm_error) 141 IF (PRESENT(which_coeffs)) which_coeffs = mm_k15 142 ENDIF 143 ELSEIF (k .LE. 53) THEN 144 CALL get_minimax_coeff_k53(k, Rc, aw, mm_error) 145 IF (PRESENT(which_coeffs)) which_coeffs = mm_k53 146 ELSE 147 CPABORT("No minimax approximations available for k = "//cp_to_string(k)) 148 ENDIF 149 END SUBROUTINE get_exp_minimax_coeff 150 151! ************************************************************************************************** 152!> \brief Get minimax coefficients: k53 implementation (coefficients up to k=53 terms). 153!> All a_i and w_i for a set of discrete values Rc, k are tabulated and 154!> the most accurate coefficients for given input k, Rc are returned. 155!> \param k ... 156!> \param Rc ... 157!> \param aw ... 158!> \param mm_error ... 159! ************************************************************************************************** 160 SUBROUTINE get_minimax_coeff_k53(k, Rc, aw, mm_error) 161 INTEGER, INTENT(IN) :: k 162 REAL(KIND=dp), INTENT(IN) :: Rc 163 REAL(KIND=dp), DIMENSION(2*k), INTENT(OUT) :: aw 164 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: mm_error 165 166 INTEGER :: i_mm 167 168 CALL get_best_minimax_approx_k53(k, Rc, i_mm) 169 CALL get_minimax_coeff_low(i_mm, aw) 170 IF (PRESENT(mm_error)) mm_error = get_minimax_numerical_error(Rc, aw) 171 172 END SUBROUTINE get_minimax_coeff_k53 173 174! ************************************************************************************************** 175!> \brief find minimax approx. with k terms that is most accurate for range Rc. 176!> \param k ... 177!> \param Rc ... 178!> \param i_mm ... 179!> \param ge_Rc Whether the tabulated range of the returned minimax approximation 180!> must be strictly greater than or equal to Rc. Default is .FALSE. 181! ************************************************************************************************** 182 SUBROUTINE get_best_minimax_approx_k53(k, Rc, i_mm, ge_Rc) 183 INTEGER, INTENT(IN) :: k 184 REAL(KIND=dp), INTENT(IN) :: Rc 185 INTEGER, INTENT(OUT) :: i_mm 186 LOGICAL, INTENT(IN), OPTIONAL :: ge_Rc 187 188 INTEGER :: i, i_k, i_l, i_r 189 REAL(KIND=dp) :: error_l, error_r, R_k_max, R_k_min 190 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: aw 191 192 CPASSERT(k .LE. k_max) 193 194 ! find k pointer and smallest and largest R_mm value for this k 195 i_k = 1 196 DO WHILE (k_mm(k_p(i_k)) .LT. k) 197 i_k = i_k + 1 198 ENDDO 199 CPASSERT(k_mm(k_p(i_k)) .EQ. k) 200 201 R_k_min = R_mm(k_p(i_k)) 202 R_k_max = R_mm(k_p(i_k + 1) - 1) 203 204 IF (Rc .GE. R_k_max) THEN 205 i_mm = k_p(i_k + 1) - 1 ! pointer to largest Rc for current k 206 ELSE IF (Rc .LE. R_k_min) THEN 207 i_mm = k_p(i_k) ! pointer to smallest Rc for current k 208 ELSE 209 i = k_p(i_k) 210 DO WHILE (Rc .GT. R_mm(i)) 211 i = i + 1 212 ENDDO 213 i_r = i ! pointer to closest R_mm >= Rc 214 i_l = i - 1 ! pointer to closest R_mm < Rc 215 216 IF (PRESENT(ge_Rc)) THEN 217 IF (ge_Rc) THEN 218 i_mm = i_r 219 RETURN 220 ENDIF 221 ENDIF 222 223 ALLOCATE (aw(2*k_mm(i_r))) 224 CALL get_minimax_coeff_low(i_r, aw) 225 error_l = get_minimax_numerical_error(Rc, aw) 226 DEALLOCATE (aw) 227 ALLOCATE (aw(2*k_mm(i_l))) 228 CALL get_minimax_coeff_low(i_l, aw) 229 error_r = get_minimax_numerical_error(Rc, aw) 230 DEALLOCATE (aw) 231 i_mm = MERGE(i_r, i_l, error_l .LE. error_r) 232 ENDIF 233 234 END SUBROUTINE get_best_minimax_approx_k53 235 236! ************************************************************************************************** 237!> \brief Unit test checking that numerical error of minimax approximations 238!> generated using any k15 or k53 coefficients is consistent with 239!> tabulated error. 240!> \param n_R Number of Rc values to be tested. 241!> \param iw ... 242! ************************************************************************************************** 243 SUBROUTINE validate_exp_minimax(n_R, iw) 244 INTEGER, INTENT(IN) :: n_R, iw 245 246 INTEGER :: i_mm, i_R, ierr, k, which_coeffs 247 LOGICAL :: do_exit 248 REAL(KIND=dp) :: dR, mm_error, R, ref_error 249 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: aw 250 251 IF (iw > 0) THEN 252 WRITE (iw, '(//T2,A)') & 253 "Unit tests for minimax 1/x ~ SUM_{i}^{K} w_i EXP(-a_i * x) for x belonging to [1:Rc]" 254 WRITE (iw, '(T2,84("*"))') 255 ENDIF 256 257 IF (iw > 0) THEN 258 WRITE (iw, '(/T2,A)') & 259 "1) checking numerical error against tabulated error at tabulated values Rc" 260 WRITE (iw, '(/T2,A)') & 261 "which coeffs, K, Rc, num. error, ref. error, rel. diff (num. error - ref. error)/(ref. error)" 262 WRITE (iw, '(T2,54("-"))') 263 ENDIF 264 DO i_mm = 1, n_approx 265 R = R_mm(i_mm) 266 k = k_mm(i_mm) 267 CALL check_exp_minimax_range(k, R, ierr) 268 IF (ierr .EQ. 0) THEN 269 ALLOCATE (aw(2*k)) 270 CALL get_exp_minimax_coeff(k, R, aw, mm_error, which_coeffs) 271 ref_error = err_mm(i_mm) 272 DEALLOCATE (aw) 273 IF (iw > 0) WRITE (iw, '(T2,A4, I3, ES10.1, ES12.3, ES12.3, ES12.3)') & 274 MERGE("k15", "k53", which_coeffs .EQ. mm_k15), k, R, & 275 mm_error, ref_error, (mm_error - ref_error)/ref_error 276 CPASSERT(mm_error .LE. ref_error*1.05_dp + 1.0E-15_dp) 277 ELSE 278 IF (iw > 0) WRITE (iw, '(T2,A4, I3, ES10.1, 3X, A)') "k15", k, R, "missing" 279 ENDIF 280 281 IF (k .LE. 15) THEN 282 ALLOCATE (aw(2*k)) 283 CALL get_minimax_coeff_k53(k, R, aw, mm_error) 284 ref_error = err_mm(i_mm) 285 DEALLOCATE (aw) 286 IF (iw > 0) WRITE (iw, '(T2,A4,I3, ES10.1, ES12.3, ES12.3, ES12.3)') & 287 "k53", k, R, mm_error, ref_error, (mm_error - ref_error)/ref_error 288 IF (mm_error .GT. ref_error*1.05_dp + 1.0E-15_dp) THEN 289 CPABORT("Test 1 failed: numerical error is larger than tabulated error") 290 ENDIF 291 ENDIF 292 ENDDO 293 294 IF (iw > 0 .AND. n_R .GT. 0) THEN 295 WRITE (iw, '(T2,54("-"))') 296 WRITE (iw, '(/T2,A)') "Test 1 OK" 297 WRITE (iw, '(/T2,A)') & 298 "2) checking numerical error against tabulated error at arbitrary values Rc" 299 WRITE (iw, '(/T2,A)') & 300 "which coeffs, K, Rc, num. error, ref. error, rel. diff (num. error - ref. error)/(ref. error)" 301 WRITE (iw, '(T2,54("-"))') 302 ENDIF 303 dR = R_max**(1.0_dp/n_R) 304 305 DO k = 1, k_max 306 ALLOCATE (aw(2*k)) 307 do_exit = .FALSE. 308 DO i_R = 1, n_R 309 R = dR**i_R 310 CALL get_exp_minimax_coeff(k, R, aw, mm_error, which_coeffs) 311 CALL get_best_minimax_approx_k53(k, R, i_mm, ge_Rc=.TRUE.) 312 IF (R .GT. R_mm(i_mm)) THEN 313 R = R_max 314 do_exit = .TRUE. 315 ENDIF 316 ref_error = err_mm(i_mm) 317 IF (iw > 0) WRITE (iw, '(T2, A4, I3, ES10.1, ES12.3, ES12.3, ES12.3)') & 318 MERGE("k15", "k53", which_coeffs .EQ. mm_k15), k, R, & 319 mm_error, ref_error, (mm_error - ref_error)/ref_error 320 IF (mm_error .GT. ref_error*1.05_dp + 1.0E-15_dp) THEN 321 CPABORT("Test 2 failed: numerical error is larger than tabulated error") 322 ENDIF 323 IF (do_exit) EXIT 324 ENDDO 325 DEALLOCATE (aw) 326 ENDDO 327 IF (iw > 0) THEN 328 WRITE (iw, '(T2,54("-"))') 329 WRITE (iw, '(/T2,A)') "Test 2 OK" 330 ENDIF 331 END SUBROUTINE validate_exp_minimax 332 333END MODULE minimax_exp 334