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