21#include "./base/base_uses.f90"
27 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'time_frequency_grids'
30 REAL(kind=
dp),
DIMENSION(:),
ALLOCATABLE :: frequency
31 REAL(kind=
dp),
DIMENSION(:),
ALLOCATABLE :: frequency_weights
32 REAL(kind=
dp),
DIMENSION(:),
ALLOCATABLE :: imaginary_time
33 REAL(kind=
dp),
DIMENSION(:),
ALLOCATABLE :: time_weights_at_zero_frequency
34 REAL(kind=
dp),
DIMENSION(:, :),
ALLOCATABLE :: cosine_time_to_frequency_weights
35 REAL(kind=
dp),
DIMENSION(:, :),
ALLOCATABLE :: cosine_frequency_to_time_weights
36 REAL(kind=
dp),
DIMENSION(:, :),
ALLOCATABLE :: sine_time_to_frequency_weights
40 sine_time_to_frequency = 2, &
41 cosine_frequency_to_time = 3
56 IF (
ALLOCATED(grid%frequency))
DEALLOCATE (grid%frequency)
57 IF (
ALLOCATED(grid%frequency_weights))
DEALLOCATE (grid%frequency_weights)
58 IF (
ALLOCATED(grid%imaginary_time))
DEALLOCATE (grid%imaginary_time)
59 IF (
ALLOCATED(grid%time_weights_at_zero_frequency))
DEALLOCATE (grid%time_weights_at_zero_frequency)
60 IF (
ALLOCATED(grid%cosine_time_to_frequency_weights))
DEALLOCATE (grid%cosine_time_to_frequency_weights)
61 IF (
ALLOCATED(grid%cosine_frequency_to_time_weights))
DEALLOCATE (grid%cosine_frequency_to_time_weights)
62 IF (
ALLOCATED(grid%sine_time_to_frequency_weights))
DEALLOCATE (grid%sine_time_to_frequency_weights)
73 INTEGER,
INTENT(IN) :: num_points
76 CHARACTER(LEN=*),
PARAMETER :: routinen =
'build_clenshaw_grid'
78 INTEGER :: handle, jquad
80 CALL timeset(routinen, handle)
82 cpassert(num_points > 0)
84 ALLOCATE (grid%frequency(num_points), grid%frequency_weights(num_points))
85 grid%frequency = 0.0_dp
86 grid%frequency_weights = 0.0_dp
88 DO jquad = 1, num_points - 1
89 grid%frequency(jquad) = jquad*
pi/(2.0_dp*num_points)
90 grid%frequency_weights(jquad) =
pi/(num_points*sin(grid%frequency(jquad))**2)
92 grid%frequency(num_points) =
pi/2.0_dp
93 grid%frequency_weights(num_points) = &
94 pi/(2.0_dp*num_points*sin(grid%frequency(num_points))**2)
121 num_points_per_magnitude, grid, build_frequency, build_time, &
122 build_transforms, build_sine, time_scaling, time_weight_scaling, &
123 max_fit_error, print_warning, unit_nr, prefer_external_backend, &
124 used_external_backend)
126 INTEGER,
INTENT(IN) :: num_points
127 REAL(kind=
dp),
INTENT(IN) :: energy_min, energy_max, regularization
128 INTEGER,
INTENT(IN) :: num_points_per_magnitude
130 LOGICAL,
INTENT(IN) :: build_frequency, build_time, &
131 build_transforms, build_sine
132 REAL(kind=
dp),
INTENT(IN) :: time_scaling, time_weight_scaling
133 REAL(kind=
dp),
INTENT(OUT) :: max_fit_error
134 LOGICAL,
INTENT(IN) :: print_warning
135 INTEGER,
INTENT(IN) :: unit_nr
136 LOGICAL,
INTENT(IN) :: prefer_external_backend
137 LOGICAL,
INTENT(OUT),
OPTIONAL :: used_external_backend
139 CHARACTER(LEN=*),
PARAMETER :: routinen =
'build_minimax_time_frequency_grid'
141 INTEGER :: external_ierr, handle, ierr
142 REAL(kind=
dp) :: e_range, max_error
143 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: points_and_weights
145 CALL timeset(routinen, handle)
147 IF (
PRESENT(used_external_backend)) used_external_backend = .false.
148 IF (prefer_external_backend)
THEN
150 grid%imaginary_time, grid%time_weights_at_zero_frequency, &
151 grid%frequency, grid%frequency_weights, &
152 grid%cosine_time_to_frequency_weights, &
153 grid%cosine_frequency_to_time_weights, &
154 grid%sine_time_to_frequency_weights, external_ierr)
155 IF (external_ierr == 0)
THEN
156 max_fit_error = 0.0_dp
157 IF (
PRESENT(used_external_backend)) used_external_backend = .true.
158 CALL timestop(handle)
164 cpassert(num_points > 0)
165 cpassert(energy_min > 0.0_dp)
166 cpassert(energy_max >= energy_min)
167 cpassert(time_scaling > 0.0_dp)
168 cpassert(time_weight_scaling > 0.0_dp)
169 cpassert(.NOT. build_transforms .OR. (build_frequency .AND. build_time))
170 cpassert(.NOT. build_sine .OR. build_transforms)
173 e_range = energy_max/energy_min
174 ALLOCATE (points_and_weights(2*num_points))
176 IF (build_frequency)
THEN
177 IF (num_points <= 20)
THEN
183 ALLOCATE (grid%frequency(num_points))
184 ALLOCATE (grid%frequency_weights(num_points))
185 grid%frequency(:) = points_and_weights(1:num_points)*energy_min
186 grid%frequency_weights(:) = points_and_weights(num_points + 1:)*energy_min
187 IF (num_points >= 26) grid%frequency_weights(:) = grid%frequency_weights(:)*4.0_dp
191 IF (num_points <= 20)
THEN
197 ALLOCATE (grid%imaginary_time(num_points))
198 ALLOCATE (grid%time_weights_at_zero_frequency(num_points))
199 grid%imaginary_time(:) = points_and_weights(1:num_points)/time_scaling/energy_min
200 grid%time_weights_at_zero_frequency(:) = points_and_weights(num_points + 1:)/time_weight_scaling/energy_min
202 IF (build_transforms)
THEN
203 ALLOCATE (grid%cosine_time_to_frequency_weights(num_points, num_points))
205 grid%cosine_time_to_frequency_weights, grid%frequency, energy_min, energy_max, &
206 max_error, num_points_per_magnitude, regularization)
208 ALLOCATE (grid%cosine_frequency_to_time_weights(num_points, num_points))
209 CALL fit_l_sq_weights(cosine_frequency_to_time, num_points, grid%imaginary_time, &
210 grid%cosine_frequency_to_time_weights, grid%frequency, energy_min, energy_max, &
211 max_error, num_points_per_magnitude, regularization)
214 ALLOCATE (grid%sine_time_to_frequency_weights(num_points, num_points))
215 CALL fit_l_sq_weights(sine_time_to_frequency, num_points, grid%imaginary_time, &
216 grid%sine_time_to_frequency_weights, grid%frequency, energy_min, energy_max, &
217 max_error, num_points_per_magnitude, regularization)
222 max_fit_error = max_error
224 DEALLOCATE (points_and_weights)
226 CALL timestop(handle)
243 SUBROUTINE fit_l_sq_weights(transform_kind, num_integ_points, tau_tj, weights, omega_tj, E_min, E_max, &
244 max_error, num_points_per_magnitude, regularization)
246 INTEGER,
INTENT(IN) :: transform_kind, num_integ_points
247 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: tau_tj
248 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT) :: weights
249 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: omega_tj
250 REAL(kind=
dp),
INTENT(IN) :: e_min, e_max
251 REAL(kind=
dp),
INTENT(OUT) :: max_error
252 INTEGER,
INTENT(IN) :: num_points_per_magnitude
253 REAL(kind=
dp),
INTENT(IN) :: regularization
255 CHARACTER(LEN=*),
PARAMETER :: routinen =
'fit_l_sq_weights'
257 INTEGER :: handle, iii, info, jjj, jquad, lwork, &
259 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: iwork
260 LOGICAL :: sine_transform, time_to_frequency
261 REAL(kind=
dp) :: func_val, max_error_tmp, multiplicator, &
263 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: sing_values, vec_uty, weight_work, work, &
265 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: mat_a, mat_sinvvsinvsigma, &
266 mat_sinvvsinvt, mat_u
268 CALL timeset(routinen, handle)
271 cpassert(transform_kind <= cosine_frequency_to_time)
272 cpassert(
SIZE(tau_tj) == num_integ_points)
273 cpassert(
SIZE(omega_tj) == num_integ_points)
274 cpassert(
SIZE(weights, 1) == num_integ_points)
275 cpassert(
SIZE(weights, 2) == num_integ_points)
277 time_to_frequency = transform_kind /= cosine_frequency_to_time
278 sine_transform = transform_kind == sine_time_to_frequency
281 num_x_nodes = (int(log10(e_max/e_min)) + 1)*num_points_per_magnitude
285 num_x_nodes = max(num_x_nodes, num_integ_points)
287 ALLOCATE (x_values(num_x_nodes), y_values(num_x_nodes), mat_a(num_x_nodes, num_integ_points), &
288 weight_work(num_integ_points), sing_values(num_integ_points), mat_u(num_x_nodes, num_x_nodes), &
289 mat_sinvvsinvt(num_x_nodes, num_integ_points), work(8*num_integ_points*num_integ_points + &
290 12*num_integ_points + 2*num_x_nodes), iwork(8*num_integ_points), &
291 mat_sinvvsinvsigma(num_integ_points, num_x_nodes), vec_uty(num_x_nodes))
292 mat_sinvvsinvsigma = 0.0_dp
297 multiplicator = (e_max/e_min)**(1.0_dp/(real(num_x_nodes, kind=
dp) - 1.0_dp))
298 DO iii = 1, num_x_nodes
299 x_values(iii) = e_min*multiplicator**(iii - 1)
304 DO jquad = 1, num_integ_points
306 IF (time_to_frequency)
THEN
307 omega = omega_tj(jquad)
309 IF (sine_transform)
THEN
311 DO iii = 1, num_x_nodes
312 y_values(iii) = 2.0_dp*omega/(x_values(iii)**2 + omega**2)
316 DO iii = 1, num_x_nodes
317 y_values(iii) = 2.0_dp*x_values(iii)/(x_values(iii)**2 + omega**2)
321 DO jjj = 1, num_integ_points
322 DO iii = 1, num_x_nodes
323 IF (sine_transform)
THEN
324 mat_a(iii, jjj) = sin(omega*tau_tj(jjj))*exp(-x_values(iii)*tau_tj(jjj))
326 mat_a(iii, jjj) = cos(omega*tau_tj(jjj))*exp(-x_values(iii)*tau_tj(jjj))
334 DO iii = 1, num_x_nodes
335 y_values(iii) = exp(-x_values(iii)*tau)
338 DO jjj = 1, num_integ_points
339 omega = omega_tj(jjj)
340 DO iii = 1, num_x_nodes
341 x_value = x_values(iii)
342 mat_a(iii, jjj) = cos(tau*omega)*2.0_dp*x_value/(x_value**2 + omega**2)
348 CALL dgesdd(
'A', num_x_nodes, num_integ_points, mat_a, num_x_nodes, sing_values, mat_u, num_x_nodes, &
349 mat_sinvvsinvt, num_x_nodes, work, lwork, iwork, info)
353 DO jjj = 1, num_integ_points
354 DO iii = 1, num_integ_points
355 mat_sinvvsinvsigma(iii, jjj) = mat_sinvvsinvt(jjj, iii)*sing_values(jjj) &
356 /(regularization**2 + sing_values(jjj)**2)
360 CALL dgemm(
'T',
'N', num_x_nodes, 1, num_x_nodes, 1.0_dp, mat_u, num_x_nodes, y_values, num_x_nodes, &
361 0.0_dp, vec_uty, num_x_nodes)
362 CALL dgemm(
'N',
'N', num_integ_points, 1, num_x_nodes, 1.0_dp, mat_sinvvsinvsigma, num_integ_points, &
363 vec_uty, num_x_nodes, 0.0_dp, weight_work, num_integ_points)
365 weights(jquad, :) = weight_work(:)
367 max_error_tmp = 0.0_dp
368 DO iii = 1, num_x_nodes
370 IF (time_to_frequency)
THEN
371 DO jjj = 1, num_integ_points
372 IF (sine_transform)
THEN
373 func_val = func_val + weight_work(jjj)*sin(omega*tau_tj(jjj))*exp(-x_values(iii)*tau_tj(jjj))
375 func_val = func_val + weight_work(jjj)*cos(omega*tau_tj(jjj))*exp(-x_values(iii)*tau_tj(jjj))
379 DO jjj = 1, num_integ_points
380 omega = omega_tj(jjj)
381 func_val = func_val + weight_work(jjj)*cos(tau*omega)*2.0_dp*x_values(iii) &
382 /(x_values(iii)**2 + omega**2)
385 max_error_tmp = max(max_error_tmp, abs(y_values(iii) - func_val))
387 max_error = max(max_error, max_error_tmp)
391 DEALLOCATE (x_values, y_values, mat_a, weight_work, sing_values, mat_u, mat_sinvvsinvt, work, iwork, &
392 mat_sinvvsinvsigma, vec_uty)
394 CALL timestop(handle)
396 END SUBROUTINE fit_l_sq_weights
405 INTEGER,
INTENT(IN) :: nr, iw
407 INTEGER :: ierr, ir, num_integ_points
408 REAL(kind=
dp) :: max_error, multiplicator, rc, rc_max
409 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: tau_tj, tj, x_tw
410 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: weights_cos_tf_t_to_w
414 multiplicator = rc_max**(1.0_dp/(real(nr, kind=
dp) - 1.0_dp))
416 DO num_integ_points = 1, 20
418 ALLOCATE (x_tw(2*num_integ_points))
420 ALLOCATE (tau_tj(num_integ_points))
422 ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
423 weights_cos_tf_t_to_w = 0.0_dp
424 ALLOCATE (tj(num_integ_points))
429 rc = 2.0_dp*multiplicator**ir
434 tj(:) = x_tw(1:num_integ_points)
440 tau_tj(:) = x_tw(1:num_integ_points)/2.0_dp
443 weights_cos_tf_t_to_w, tj, 1.0_dp, rc, max_error, 200, 0.0_dp)
446 WRITE (iw,
'(T2, I3, F12.1, ES12.3)') num_integ_points, rc, max_error
451 DEALLOCATE (x_tw, tau_tj, weights_cos_tf_t_to_w, tj)
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Interface to the Greenx library.
subroutine, public greenx_get_minimax_grid(unit_nr, num_integ_points, emin, emax, regularization_minimax, imaginary_time, time_weights_at_zero_frequency, frequency, frequency_weights, cosine_time_to_frequency_weights, cosine_frequency_to_time_weights, sine_time_to_frequency_weights, ierr)
Get a minimax grid from GreenX when it is available.
Defines the basic variable types.
integer, parameter, public dp
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Routines to calculate the minimax coefficients in order to approximate 1/x as a sum over exponential ...
subroutine, public get_exp_minimax_coeff_gw(k, e_range, aw)
...
Routines to calculate the minimax coefficients in order to approximate 1/x as a sum over exponential ...
subroutine, public get_exp_minimax_coeff(k, rc, aw, mm_error, which_coeffs)
Get best minimax approximation for given input parameters. Automatically chooses the most exact set o...
Routines to calculate the minimax coefficients for approximating 1/x as 1/x ~ 1/pi SUM_{i}...
subroutine, public get_rpa_minimax_coeff_larger_grid(k, e_range, aw)
...
subroutine, public get_rpa_minimax_coeff(k, e_range, aw, ierr, print_warning)
The a_i and w_i coefficient are stored in aw such that the first 1:K elements correspond to a_i and t...
Definition and construction of time/frequency grids for correlation methods.
integer, parameter cosine_time_to_frequency
subroutine, public build_minimax_time_frequency_grid(num_points, energy_min, energy_max, regularization, num_points_per_magnitude, grid, build_frequency, build_time, build_transforms, build_sine, time_scaling, time_weight_scaling, max_fit_error, print_warning, unit_nr, prefer_external_backend, used_external_backend)
Build a minimax time/frequency grid through the common backend boundary.
subroutine, public build_clenshaw_grid(num_points, grid)
Build a Clenshaw-Curtis frequency grid.
subroutine, public test_least_square_ft(nr, iw)
test the singular value decomposition for the computation of integration weights for the Fourier tran...
subroutine, public time_frequency_grid_release(grid)
Release all data owned by a time_frequency_grid_type object.