Line data Source code
1 : ! **************************************************************************************************
2 : ! Copyright (C) 2020-2023 Green-X library
3 : ! This file is distributed under the terms of the APACHE2 License.
4 : !
5 : ! **************************************************************************************************
6 : !> \brief Routines to calculate frequency and time grids (integration points and weights)
7 : !> for correlation methods as well as weights for the inhomogeneous cosine/sine transform.
8 : !>
9 : !> NB: When dealing with runtime exceptions, we set ierr to a non-zero value and return immediately
10 : ! to the caller so we don't need to goto to a cleanup section at the end of the procedure.
11 : ! Assume -std=f2008: i.e. allocatable arrays are automatically deallocated when going out of scope.
12 : !> reference: [https://doi.org/10.1021/ct5001268](https://doi.org/10.1021/ct5001268)
13 : !> reference: [https://doi.org/10.1103/PhysRevB.94.165109](https://doi.org/10.1103/PhysRevB.94.165109)
14 : ! **************************************************************************************************
15 :
16 : module minimax_grids
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 : #include "abi_common.h"
21 :
22 : #include "gx_common.h"
23 : use defs_basis, only: dp, pi
24 : use m_errors
25 : !use kinds, only: dp
26 : !use error_handling, only: register_exc
27 : !use constants, only: pi
28 : use minimax_tau, only: get_points_weights_tau
29 : use minimax_omega, only: get_points_weights_omega
30 : use minimax_utils, only: cosine_wt, cosine_tw, sine_tw
31 : !use lapack_interfaces, only: dgemm, dgesdd
32 :
33 : implicit none
34 :
35 : private
36 :
37 : !> Main entry point for client code.
38 : public :: gx_minimax_grid, gx_minimax_grid_frequency
39 :
40 : contains
41 :
42 : !> \brief Compute minimax grid for GW calculation on imaginary time/frequency domain.
43 : !! @param[in] num_points: Number of mesh points.
44 : !! @param[in] e_min: Minimum transition energy. Arbitrary units as we only need e_max/e_min
45 : !! @param[in] e_max: Maximum transition energy.
46 : !! @param[out] tau_points: imaginary time grid points
47 : !! @param[out] tau_weights: weights for imaginary time grid points weights
48 : !! @param[out] omega_points: imaginary frequency grid points
49 : !! @param[out] omega_weights: weights for imaginary frequency grid points
50 : !! @param[out] cosft_wt: weights for tau -> omega cosine transform. cos(w*t) factor is included.
51 : !! @param[out] cosft_tw: weights for omega -> tau cosine transform. cos(w*t) factor is included.
52 : !! @param[out] sinft_wt: weights for tau -> omega sine transform. sin(w*t) factor is included.
53 : !! @param[out] max_errors: Max error for the three kind of transforms (same order as previous args)
54 : !! @param[out] cosft_duality_error. Max_{ij} |AB - I| where A and B are the cosft_wt and cosft_tw matrices.
55 : !! @param[out] ierr: Exit status
56 : !! @param[in] bare_cos_sin_weights: whether the cosine and sine weights are multiplied by cos and sin term, optional
57 2 : subroutine gx_minimax_grid(num_points, e_min, e_max, &
58 : tau_points, tau_weights, omega_points, omega_weights, &
59 : cosft_wt, cosft_tw, sinft_wt, &
60 : max_errors, cosft_duality_error, ierr, bare_cos_sin_weights, regterm)
61 :
62 : integer, intent(in) :: num_points
63 : real(kind=dp), intent(in) :: e_min, e_max
64 : real(kind=dp), allocatable, dimension(:), &
65 : intent(out) :: tau_points, tau_weights
66 : real(kind=dp), allocatable, dimension(:), &
67 : intent(out) :: omega_points, omega_weights
68 : real(kind=dp), allocatable, dimension(:, :), &
69 : intent(out) :: cosft_wt, cosft_tw, sinft_wt
70 : real(kind=dp), intent(out) :: max_errors(3), cosft_duality_error
71 : logical, intent(in), optional :: bare_cos_sin_weights
72 : integer, intent(out) :: ierr
73 : real(kind=dp),optional,intent(in) :: regterm
74 :
75 : ! Internal variables
76 : logical :: my_bare_cos_sin_weights
77 : integer, parameter :: cos_t_to_cos_w = 1
78 : integer, parameter :: cos_w_to_cos_t = 2
79 : integer, parameter :: sin_t_to_sin_w = 3
80 : integer :: i_point, j_point
81 : real(kind=dp) :: e_range, scaling, regterm__
82 2 : real(kind=dp), dimension(:), allocatable :: x_tw
83 2 : real(kind=dp), dimension(:, :), allocatable :: mat
84 2 : real(kind=dp), dimension(:, :), allocatable :: tmp_cosft_wt, tmp_cosft_tw
85 :
86 2 : my_bare_cos_sin_weights = .false.
87 2 : if (present(bare_cos_sin_weights)) then
88 0 : my_bare_cos_sin_weights = bare_cos_sin_weights
89 : endif
90 :
91 2 : regterm__ = 0.0_dp; if (present(regterm)) regterm__ = regterm
92 :
93 : ! Begin work
94 2 : e_range = e_max/e_min
95 2 : ierr = 0
96 :
97 : ! Allocations
98 6 : allocate (x_tw(2*num_points))
99 2 : if (.not. allocated(omega_points)) then
100 6 : allocate (omega_points(num_points))
101 : end if
102 2 : if (.not. allocated(omega_weights)) then
103 6 : allocate (omega_weights(num_points))
104 : end if
105 2 : if (.not. allocated(tau_points)) then
106 6 : allocate (tau_points(num_points))
107 : end if
108 2 : if (.not. allocated(tau_weights)) then
109 6 : allocate (tau_weights(num_points))
110 : end if
111 :
112 : ! Get the frequency grid points and weights
113 2 : call get_points_weights_omega(num_points, e_range, x_tw, ierr)
114 2 : if (ierr /= 0) return
115 :
116 : ! Scale the frequency grid points and weights from [1,R] to [e_min,e_max]
117 : ! Note: the frequency grid points and weights include a factor of two
118 2 : scaling = e_min
119 32 : omega_points(:) = x_tw(1: num_points) *scaling
120 32 : omega_weights(:) = x_tw(num_points+1: 2* num_points) *scaling
121 :
122 : ! Get the time grid points and weights
123 2 : call get_points_weights_tau(num_points, e_range, x_tw, ierr)
124 2 : if (ierr /= 0) return
125 :
126 : ! Scale the time grid points and weights from [1,R] to [e_min,e_max]
127 2 : scaling = 2.0_dp *e_min
128 32 : tau_points(:) = x_tw(1:num_points) / scaling
129 32 : tau_weights(:) = x_tw(num_points+1:2*num_points) / scaling
130 :
131 8 : allocate (cosft_wt(num_points, num_points))
132 6 : allocate (cosft_tw(num_points, num_points))
133 6 : allocate (sinft_wt(num_points, num_points))
134 6 : allocate (tmp_cosft_wt(num_points, num_points))
135 6 : allocate (tmp_cosft_tw(num_points, num_points))
136 :
137 : ! get the weights for the cosine transform W^c(it) -> W^c(iw)
138 : call get_transformation_weights(num_points, tau_points, omega_points, cosft_wt, e_min, e_max, &
139 2 : max_errors(1), cos_t_to_cos_w, regterm__, ierr)
140 2 : if (ierr /= 0) return
141 :
142 : ! get the weights for the cosine transform W^c(iw) -> W^c(it)
143 : call get_transformation_weights(num_points, tau_points, omega_points, cosft_tw, e_min, e_max, &
144 2 : max_errors(2), cos_w_to_cos_t, regterm__, ierr)
145 2 : if (ierr /= 0) return
146 :
147 : ! get the weights for the sine transform Sigma^sin(it) -> Sigma^sin(iw) (PRB 94, 165109 (2016), Eq. 71)
148 : call get_transformation_weights(num_points, tau_points, omega_points, sinft_wt, e_min, e_max, &
149 2 : max_errors(3), sin_t_to_sin_w, regterm__, ierr)
150 2 : if (ierr /= 0) return
151 :
152 : ! Compute the actual weights used for the inhomogeneous cosine/ FT and check whether
153 : ! the two matrices for the forward/backward transform are the inverse of each other.
154 2 : if(.not.my_bare_cos_sin_weights) then
155 32 : do j_point = 1, num_points
156 532 : do i_point = 1, num_points
157 500 : cosft_wt(j_point, i_point) = cosft_wt(j_point, i_point)*cos(tau_points(i_point)*omega_points(j_point))
158 500 : cosft_tw(i_point, j_point) = cosft_tw(i_point, j_point)*cos(tau_points(i_point)*omega_points(j_point))
159 530 : sinft_wt(j_point, i_point) = sinft_wt(j_point, i_point)*sin(tau_points(i_point)*omega_points(j_point))
160 : end do
161 : end do
162 : else
163 0 : do j_point = 1, num_points
164 0 : do i_point = 1, num_points
165 0 : tmp_cosft_wt(j_point, i_point) = cosft_wt(j_point, i_point)*cos(tau_points(i_point)*omega_points(j_point))
166 0 : tmp_cosft_tw(i_point, j_point) = cosft_tw(i_point, j_point)*cos(tau_points(i_point)*omega_points(j_point))
167 : end do
168 : end do
169 : end if
170 :
171 6 : allocate (mat(num_points, num_points))
172 2 : if(.not.my_bare_cos_sin_weights) then
173 37062 : mat(:,:) = matmul(cosft_wt, cosft_tw)
174 : else
175 0 : mat(:,:) = matmul(tmp_cosft_wt, tmp_cosft_tw)
176 : endif
177 :
178 32 : do i_point = 1, num_points
179 32 : mat(i_point, i_point) = mat(i_point, i_point) - 1.0_dp
180 : end do
181 532 : cosft_duality_error = maxval(abs(mat))
182 :
183 2 : deallocate (mat)
184 2 : deallocate (x_tw)
185 2 : deallocate (tmp_cosft_wt,tmp_cosft_tw)
186 :
187 4 : end subroutine gx_minimax_grid
188 :
189 : !> \brief Retrieves the frequency grid for a canonical GW/RPA calculation
190 : !! @param[in] num_points: Number of mesh points.
191 : !! @param[in] e_min: Minimum transition energy. Arbitrary units as we only need e_max/e_min
192 : !! @param[in] e_max: Maximum transition energy.
193 : !! @param[out] omega_points: imaginary frequency grid points
194 : !! @param[out] omega_weights: weights for imaginary frequency grid points
195 : !! @param[out] ierr: Exit status
196 0 : subroutine gx_minimax_grid_frequency (num_points, e_min, e_max, omega_points, omega_weights, ierr)
197 : integer, intent(in) :: num_points
198 : real(kind=dp), intent(in) :: e_min, e_max
199 : real(kind=dp), allocatable, dimension(:), &
200 : intent(out) :: omega_points(:), omega_weights(:)
201 : integer, intent(out) :: ierr
202 :
203 : ! Internal variables
204 : real(kind=dp) :: e_range, scaling
205 0 : real(kind=dp), dimension(:), allocatable :: x_tw
206 :
207 : ! Begin work
208 0 : e_range = e_max/e_min
209 0 : ierr = 0
210 :
211 : ! Allocations
212 0 : allocate (x_tw(2*num_points))
213 0 : if (.not. allocated(omega_points)) then
214 0 : allocate (omega_points(num_points))
215 : end if
216 0 : if (.not. allocated(omega_weights)) then
217 0 : allocate (omega_weights(num_points))
218 : end if
219 :
220 : ! Get the frequency grid points and weights
221 0 : call get_points_weights_omega(num_points, e_range, x_tw, ierr)
222 0 : if (ierr /= 0) return
223 :
224 : ! Scale the frequency grid points and weights from [1,R] to [e_min,e_max]
225 : ! Note: the frequency grid points and weights include a factor of two
226 0 : scaling = e_min
227 0 : omega_points(:) = x_tw(1: num_points) *scaling
228 0 : omega_weights(:) = x_tw(num_points+1: 2* num_points) *scaling
229 :
230 0 : deallocate (x_tw)
231 :
232 0 : end subroutine gx_minimax_grid_frequency
233 :
234 :
235 : !> \brief Get the weights eiter for the cosine/sin transformation for tau to omega or viceversa
236 : !! @param[in] num_points: Number of mesh points.
237 : !! @param[in] tau_points: imaginary time grid points
238 : !! @param[in] omega_points: imaginary frequency grid points
239 : !! @param[inout] weights: corresponding tranformation weights
240 : !! @param[in] e_min: Minimum transition energy.
241 : !! @param[in] e_max: Maximum transition energy.
242 : !! @param[inout] max_error: Max error for the transform
243 : !! @param[in] transformation type : 1 the cosine transform cos(it) -> cos(iw)
244 : !! : 2 the cosine transform cos(iw) -> cos(it)
245 : !! : 3 the sine transform sin(it) -> sin(iw)
246 : !! @param[in] ierr: exit status
247 6 : subroutine get_transformation_weights(num_points, tau_points, omega_points, weights, e_min, e_max, &
248 : max_error, transformation_type, regterm, ierr)
249 :
250 : integer, intent(in) :: num_points
251 : real(kind=dp), allocatable, dimension(:), &
252 : intent(in) :: tau_points
253 : real(kind=dp), allocatable, dimension(:), &
254 : intent(in) :: omega_points
255 : real(kind=dp), allocatable, dimension(:, :), &
256 : intent(inout) :: weights
257 : real(kind=dp), intent(in) :: e_min, e_max
258 : real(kind=dp), intent(inout) :: max_error
259 : integer, intent(in) :: transformation_type
260 : integer, intent(out) :: ierr
261 : real(kind=dp),intent(in) :: regterm
262 :
263 : ! Internal variables
264 : integer :: i_node, i_point, j_point, k_point, &
265 : num_x_nodes
266 : integer, parameter :: nodes_factor = 200
267 : real(kind=dp) :: current_point, x_factor
268 : real(kind=dp), allocatable, dimension(:) :: weights_work, x_mu, psi
269 : real(kind=dp), allocatable, dimension(:, :) :: mat_A
270 :
271 : integer :: lwork
272 6 : integer, allocatable, dimension(:) :: iwork
273 6 : real(kind=dp), allocatable, dimension(:) :: vec_S, vec_UT_psi, work
274 6 : real(kind=dp), allocatable, dimension(:, :) :: mat_U, mat_VT, mat_VT_s
275 :
276 : ! Begin work
277 6 : ierr = 0
278 :
279 108 : allocate (weights_work(num_points), source=0.0_dp)
280 :
281 : ! compute the number of x nodes per magnitude points per 10-interval
282 6 : num_x_nodes = (int(log10(e_max/e_min)) + 1)*nodes_factor
283 :
284 : ! make sure that the number of x nodes are at least as many integration points
285 6 : num_x_nodes = max(num_x_nodes, num_points)
286 :
287 3018 : allocate (x_mu(num_x_nodes), source=0.0_dp)
288 3012 : allocate (psi(num_x_nodes), source=0.0_dp)
289 :
290 : ! Allocations for the BLAS routines
291 : ! double the value nessary for 'A' to achieve good performance
292 6 : lwork = 8*num_points*num_points + 12*num_points + 2*num_x_nodes
293 738 : allocate (iwork(8*num_points), source=0)
294 19098 : allocate (work(lwork), source=0.0_dp)
295 :
296 42114 : allocate (mat_A(num_x_nodes, num_points), source=0.0_dp)
297 1563024 : allocate (mat_U(num_x_nodes, num_x_nodes), source=0.0_dp)
298 42108 : allocate (mat_VT(num_x_nodes, num_points), source=0.0_dp)
299 45018 : allocate (mat_VT_s(num_points, num_x_nodes), source=0.0_dp)
300 102 : allocate (vec_S(num_points), source=0.0_dp)
301 3012 : allocate (vec_UT_psi(num_x_nodes), source=0.0_dp)
302 :
303 : ! set the x-mu logarithmically in the interval [e_min,e_max]
304 6 : x_factor = (e_max/e_min)**(1.0_dp/(real(num_x_nodes, kind=dp) - 1.0_dp))
305 3006 : do i_node = 1, num_x_nodes
306 3006 : x_mu(i_node) = e_min*x_factor**(i_node - 1)
307 : end do
308 :
309 6 : current_point = 0.0_dp
310 6 : max_error = 0.0_dp
311 :
312 : ! loop over all grid points
313 96 : do i_point = 1, num_points
314 : ! calculate psi and mat_A
315 : call calculate_psi_and_mat_A(num_points, tau_points, omega_points, num_x_nodes, x_mu, psi, &
316 90 : mat_A, i_point, current_point, transformation_type)
317 :
318 : ! Singular value decomposition of mat_A = U*Sigma*V^T
319 : call dgesdd('A', num_x_nodes, num_points, mat_A, num_x_nodes, vec_S, mat_U, num_x_nodes, &
320 90 : mat_VT, num_x_nodes, work, lwork, iwork, ierr)
321 :
322 90 : if (ierr /= 0) then
323 0 : _REGISTER_EXC("DGESDD returned ierr != 0")
324 0 : return
325 : end if
326 :
327 : ! integration weights = (V Sigma^-1 U^T)*psi
328 :
329 : ! 1) V * Sigma^-1
330 : !regterm = 0.01_dp
331 1590 : do j_point = 1, num_points
332 28590 : do k_point = 1, num_points
333 28500 : if (regterm > tiny(regterm)) then
334 : mat_VT_s(k_point, j_point) = mat_VT(j_point, k_point)*vec_S(j_point) / &
335 0 : (vec_S(j_point)**2+regterm**2)
336 : else
337 27000 : mat_VT_s(k_point, j_point) = mat_VT(j_point, k_point) / vec_S(j_point)
338 : end if
339 : end do ! k_point
340 : end do ! j_point
341 :
342 : ! 2) (U^T)*psi
343 : call dgemm('T', 'N', num_x_nodes, 1, num_x_nodes, 1.0_dp, mat_U, num_x_nodes, psi, num_x_nodes, &
344 90 : 0.0_dp, vec_UT_psi, num_x_nodes)
345 :
346 : ! 3) (V*Sigma^-1) * (U^T*psi)
347 : call dgemm('N', 'N', num_points, 1, num_x_nodes, 1.0_dp, mat_VT_s, num_points, vec_UT_psi, &
348 90 : num_x_nodes, 0.0_dp, weights_work, num_points)
349 :
350 1590 : weights(i_point, :) = weights_work(:)
351 :
352 : ! calculate the maximum error of the fitting
353 : call calculate_max_error(num_points, tau_points, omega_points, weights_work, num_x_nodes, x_mu, &
354 96 : psi, current_point, max_error, transformation_type)
355 : end do ! i_point
356 :
357 6 : deallocate (x_mu, psi, mat_A, weights_work, vec_S, mat_U, mat_VT, work, iwork, mat_VT_s, vec_UT_psi)
358 :
359 6 : end subroutine get_transformation_weights
360 :
361 : !> \brief Calculate the auxiliary matrix for cosine/sin transformation for tau to omega or viceversa
362 : !! @param[in] num_points: Number of mesh points.
363 : !! @param[in] tau_points: imaginary time grid points
364 : !! @param[in] omega_points: imaginary frequency grid points
365 : !! @param[in] num_x_nodes: Number of node in the interval [e_min,e_max]
366 : !! @param[in] x_mu : Transition energy (nodes in the interval [e_min,e_max])
367 : !! @param[inout] psi: corresponding auxiliary function (see transformation type definition)
368 : !! @param[inout] mat_A: auxiliary matrix (see transformation type definition)
369 : !! @param[in] i_point: pointer for the current grid point
370 : !! @param[inout] current_point: current grid point ether omega(i_point) or tau_(i_point)
371 : !! @param[in] transformation type :
372 : !! (1) the cosine transform cos(it) -> cos(iw): psi(omega,x), mat_A = cos(omega*tau)*psi(tau,x)
373 : !! (2) the cosine transform cos(iw) -> cos(it): psi(tau,x) , mat_A = cos(omega*tau)*psi(omega,x)
374 : !! (3) the sine transform sin(it) -> sin(iw): psi(omega,x), mat_A = sin(omega*tau)*psi(tau,x)
375 90 : subroutine calculate_psi_and_mat_A(num_points, tau_points, omega_points, num_x_nodes, x_mu, psi, &
376 : mat_A, i_point, current_point, transformation_type)
377 :
378 : integer, intent(in) :: num_points, num_x_nodes, i_point
379 : real(kind=dp), allocatable, dimension(:), &
380 : intent(in) :: tau_points, omega_points, x_mu
381 : real(kind=dp), allocatable, dimension(:), &
382 : intent(inout) :: psi
383 : real(kind=dp), allocatable, dimension(:, :), &
384 : intent(inout) :: mat_A
385 : real(kind=dp), intent(inout) :: current_point
386 : integer, intent(in) :: transformation_type
387 :
388 : ! Internal variables
389 : integer :: i_node, j_point
390 : real(kind=dp) :: tau, omega
391 :
392 : ! Begin work
393 :
394 : ! the cosine transform cos(it) -> cos(iw)
395 90 : if (transformation_type == cosine_tw) then
396 30 : omega = omega_points(i_point)
397 30 : current_point = omega
398 :
399 : ! psi(omega_k,x) = 2x/(x^2+omega_k^2)
400 14030 : do i_node = 1, num_x_nodes
401 14030 : psi(i_node) = 2.0_dp*x_mu(i_node)/((x_mu(i_node))**2 + omega**2)
402 : end do
403 :
404 : ! mat_A = cos(omega_k * tau) psi(tau,x)
405 530 : do j_point = 1, num_points
406 500 : tau = tau_points(j_point)
407 220530 : do i_node = 1, num_x_nodes
408 220500 : mat_A(i_node, j_point) = cos(omega*tau)*exp(-x_mu(i_node)*tau)
409 : end do
410 : end do
411 :
412 : ! the cosine transform cos(iw) -> cos(it)
413 60 : else if (transformation_type == cosine_wt) then
414 30 : tau = tau_points(i_point)
415 30 : current_point = tau
416 :
417 : ! psi(tau_k,x) = =exp(-x*|tau_k|)
418 14030 : do i_node = 1, num_x_nodes
419 14030 : psi(i_node) = exp(-x_mu(i_node)*tau)
420 : end do
421 :
422 : ! mat_A = cos(tau_k,omega) psi(omega,x)
423 530 : do j_point = 1, num_points
424 500 : omega = omega_points(j_point)
425 220530 : do i_node = 1, num_x_nodes
426 220500 : mat_A(i_node, j_point) = cos(tau*omega)*2.0_dp*x_mu(i_node)/(x_mu(i_node)**2 + omega**2)
427 : end do
428 : end do
429 :
430 : ! the sine transform sin(it) -> sin(iw)
431 30 : else if (transformation_type == sine_tw) then
432 30 : omega = omega_points(i_point)
433 30 : current_point = omega
434 :
435 : ! psi(omega_k,x) = 2*omega_k/(x^2+omega_k^2)
436 14030 : do i_node = 1, num_x_nodes
437 14030 : psi(i_node) = 2.0_dp*omega/(x_mu(i_node)**2 + omega**2)
438 : end do
439 :
440 : ! mat_A = sin(omega_k,tau)*phi(tau,x)
441 530 : do j_point = 1, num_points
442 500 : tau = tau_points(j_point)
443 220530 : do i_node = 1, num_x_nodes
444 220500 : mat_A(i_node, j_point) = sin(omega*tau)*exp(-x_mu(i_node)*tau)
445 : end do
446 : end do
447 : end if
448 :
449 90 : end subroutine calculate_psi_and_mat_A
450 :
451 : !> /brief calculate the error of the fit function for the cosine/sin transformation for tau to omega or viceversa
452 : !! @param[in] num_points: Number of mesh points.
453 : !! @param[in] tau_points: imaginary time grid points
454 : !! @param[in] omega_points: imaginary frequency grid points
455 : !! @param[in] weights_work: work vector for the transformation weights
456 : !! @param[in] num_x_nodes: Number of node in the interval [e_min,e_max]
457 : !! @param[in] x_mu : Transition energy (nodes in the interval [e_min,e_max])
458 : !! @param[in] psi: corresponding auxiliary function (see transformation type definition)
459 : !! @param[in] current_point: current grid point ether omega(i_point) or tau_(i_point)
460 : !! @param[out] max_error: Max error for the transform.
461 : !! @param[in] transformation type : 1 fit function for the cosine transform cos(it) -> cos(iw), psi(omeaga,x)
462 : !! : 2 fit function for the cosine transform cos(iw) -> cos(it), psi(tau,x)
463 : !! : 3 fit function for the sine transform sin(it) -> sin(iw), psi(omega,x)
464 90 : subroutine calculate_max_error(num_points, tau_points, omega_points, weights_work, num_x_nodes, x_mu, &
465 : psi, current_point, max_error, transformation_type)
466 :
467 : real(kind=dp), intent(out) :: max_error
468 : real(kind=dp), intent(in) :: current_point
469 : real(kind=dp), allocatable, dimension(:), &
470 : intent(in) :: tau_points, omega_points, x_mu, psi, &
471 : weights_work
472 : integer, intent(in) :: num_points, num_x_nodes
473 : integer, intent(in) :: transformation_type
474 :
475 : ! Internal variables
476 : integer :: i_node,i_point
477 : real(kind=dp) :: func_val, func_val_temp, max_error_tmp, &
478 : tau, omega, x_val
479 :
480 : ! Begin work
481 90 : max_error_tmp = 0.0_dp
482 :
483 : ! the cosine transform cos(it) -> cos(iw)
484 90 : if (transformation_type == cosine_tw) then
485 30 : omega=current_point
486 :
487 14030 : do i_node = 1, num_x_nodes
488 14000 : func_val = 0.0_dp
489 : ! calculate value of the fit function f(x) = f(x) + weights(omega)cos(omega*tau)psi(tau.x)
490 234000 : do i_point = 1, num_points
491 220000 : tau = tau_points(i_point)
492 234000 : func_val = func_val + weights_work(i_point)*cos(omega*tau)*exp(-x_mu(i_node)*tau)
493 : end do
494 :
495 14030 : if (abs(psi(i_node) - func_val) > max_error_tmp) then
496 : max_error_tmp = abs(psi(i_node) - func_val)
497 : func_val_temp = func_val
498 : end if
499 : end do
500 :
501 : ! the cosine transform cos(iw) -> cos(it)
502 60 : else if (transformation_type == cosine_wt) then
503 30 : tau = current_point
504 :
505 14030 : do i_node = 1, num_x_nodes
506 14000 : func_val = 0.0_dp
507 14000 : x_val=x_mu(i_node)
508 : ! calculate value of the fit function f(x) = f(x) + weights(tau)cos(omega*tau)psi(omega.x)
509 234000 : do i_point = 1, num_points
510 220000 : omega = omega_points(i_point)
511 234000 : func_val = func_val + weights_work(i_point)*cos(tau*omega)*2.0_dp*x_val/(x_val**2 + omega**2)
512 : end do
513 :
514 14030 : if (abs(psi(i_node) - func_val) > max_error_tmp) then
515 : max_error_tmp = abs(psi(i_node) - func_val)
516 : func_val_temp = func_val
517 : end if
518 : end do
519 :
520 : ! the sine transform sin(it) -> sin(iw)
521 30 : else if (transformation_type == sine_tw) then
522 30 : omega = current_point
523 :
524 14030 : do i_node = 1, num_x_nodes
525 14000 : func_val = 0.0_dp
526 : ! calculate value of the fit function f(x) = f(x) + weights(omega)sin(omega*tau)psi(tau.x)
527 234000 : do i_point = 1, num_points
528 220000 : tau = tau_points(i_point)
529 234000 : func_val = func_val + weights_work(i_point)*sin(omega*tau)*exp(-x_mu(i_node)*tau)
530 : end do
531 :
532 14030 : if (abs(psi(i_node) - func_val) > max_error_tmp) then
533 : max_error_tmp = abs(psi(i_node) - func_val)
534 : func_val_temp = func_val
535 : end if
536 : end do
537 : end if
538 :
539 90 : if (max_error_tmp > max_error) then
540 13 : max_error = max_error_tmp
541 : end if
542 :
543 90 : end subroutine calculate_max_error
544 :
545 : end module minimax_grids
|
