Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Definition and construction of time/frequency grids for correlation methods.
10 : !> \par History
11 : !> 05.2019 Refactored from rpa_ri_gpw [Frederick Stein]
12 : ! **************************************************************************************************
13 : MODULE time_frequency_grids
14 : USE greenx_interface, ONLY: greenx_get_minimax_grid
15 : USE kinds, ONLY: dp
16 : USE mathconstants, ONLY: pi
17 : USE minimax_exp, ONLY: get_exp_minimax_coeff
18 : USE minimax_exp_gw, ONLY: get_exp_minimax_coeff_gw
19 : USE minimax_rpa, ONLY: get_rpa_minimax_coeff,&
20 : get_rpa_minimax_coeff_larger_grid
21 : #include "./base/base_uses.f90"
22 :
23 : IMPLICIT NONE
24 :
25 : PRIVATE
26 :
27 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'time_frequency_grids'
28 :
29 : TYPE time_frequency_grid_type
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
37 : END TYPE time_frequency_grid_type
38 :
39 : INTEGER, PARAMETER :: cosine_time_to_frequency = 1, &
40 : sine_time_to_frequency = 2, &
41 : cosine_frequency_to_time = 3
42 :
43 : PUBLIC :: build_clenshaw_grid, build_minimax_time_frequency_grid, test_least_square_ft, &
44 : time_frequency_grid_release, time_frequency_grid_type
45 :
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief Release all data owned by a time_frequency_grid_type object.
50 : !> \param grid The grid object to release.
51 : ! **************************************************************************************************
52 1092 : SUBROUTINE time_frequency_grid_release(grid)
53 :
54 : TYPE(time_frequency_grid_type), INTENT(INOUT) :: grid
55 :
56 1092 : IF (ALLOCATED(grid%frequency)) DEALLOCATE (grid%frequency)
57 1092 : IF (ALLOCATED(grid%frequency_weights)) DEALLOCATE (grid%frequency_weights)
58 1092 : IF (ALLOCATED(grid%imaginary_time)) DEALLOCATE (grid%imaginary_time)
59 1092 : IF (ALLOCATED(grid%time_weights_at_zero_frequency)) DEALLOCATE (grid%time_weights_at_zero_frequency)
60 1092 : IF (ALLOCATED(grid%cosine_time_to_frequency_weights)) DEALLOCATE (grid%cosine_time_to_frequency_weights)
61 1092 : IF (ALLOCATED(grid%cosine_frequency_to_time_weights)) DEALLOCATE (grid%cosine_frequency_to_time_weights)
62 1092 : IF (ALLOCATED(grid%sine_time_to_frequency_weights)) DEALLOCATE (grid%sine_time_to_frequency_weights)
63 :
64 1092 : END SUBROUTINE time_frequency_grid_release
65 :
66 : ! **************************************************************************************************
67 : !> \brief Build a Clenshaw-Curtis frequency grid.
68 : !> \param num_points Number of integration points.
69 : !> \param grid Grid object to fill.
70 : ! **************************************************************************************************
71 116 : SUBROUTINE build_clenshaw_grid(num_points, grid)
72 :
73 : INTEGER, INTENT(IN) :: num_points
74 : TYPE(time_frequency_grid_type), INTENT(OUT) :: grid
75 :
76 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_clenshaw_grid'
77 :
78 : INTEGER :: handle, jquad
79 :
80 116 : CALL timeset(routineN, handle)
81 :
82 116 : CPASSERT(num_points > 0)
83 :
84 464 : ALLOCATE (grid%frequency(num_points), grid%frequency_weights(num_points))
85 5186 : grid%frequency = 0.0_dp
86 5186 : grid%frequency_weights = 0.0_dp
87 :
88 5070 : DO jquad = 1, num_points - 1
89 4954 : grid%frequency(jquad) = jquad*pi/(2.0_dp*num_points)
90 5070 : grid%frequency_weights(jquad) = pi/(num_points*SIN(grid%frequency(jquad))**2)
91 : END DO
92 116 : grid%frequency(num_points) = pi/2.0_dp
93 : grid%frequency_weights(num_points) = &
94 116 : pi/(2.0_dp*num_points*SIN(grid%frequency(num_points))**2)
95 :
96 116 : CALL timestop(handle)
97 :
98 116 : END SUBROUTINE build_clenshaw_grid
99 :
100 : ! **************************************************************************************************
101 : !> \brief Build a minimax time/frequency grid through the common backend boundary.
102 : !> \param num_points Number of minimax points.
103 : !> \param energy_min Lower end of the physical energy interval.
104 : !> \param energy_max Upper end of the physical energy interval.
105 : !> \param regularization Regularization used for the fitted transform weights.
106 : !> \param num_points_per_magnitude Number of fitting points per decade.
107 : !> \param grid Grid object to fill.
108 : !> \param build_frequency Whether to construct the frequency components.
109 : !> \param build_time Whether to construct the imaginary-time components.
110 : !> \param build_transforms Whether to construct the time/frequency transform weights.
111 : !> \param build_sine Whether to construct the sine time-to-frequency weights.
112 : !> \param time_scaling Scaling applied to the imaginary-time abscissas.
113 : !> \param time_weight_scaling Scaling applied to the imaginary-time weights.
114 : !> \param max_fit_error Maximum fitting error across all requested transforms.
115 : !> \param print_warning Whether the minimax coefficient routine prints warnings.
116 : !> \param unit_nr Output unit used by an external backend.
117 : !> \param prefer_external_backend Whether to try an available external backend first.
118 : !> \param used_external_backend Whether the external backend supplied the grid.
119 : ! **************************************************************************************************
120 338 : SUBROUTINE build_minimax_time_frequency_grid(num_points, energy_min, energy_max, regularization, &
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)
125 :
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
129 : TYPE(time_frequency_grid_type), INTENT(OUT) :: grid
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
138 :
139 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_minimax_time_frequency_grid'
140 :
141 : INTEGER :: external_ierr, handle, ierr
142 : REAL(KIND=dp) :: e_range, max_error
143 338 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: points_and_weights
144 :
145 338 : CALL timeset(routineN, handle)
146 :
147 338 : IF (PRESENT(used_external_backend)) used_external_backend = .FALSE.
148 338 : IF (prefer_external_backend) THEN
149 : CALL greenx_get_minimax_grid(unit_nr, num_points, energy_min, energy_max, regularization, &
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 214 : grid%sine_time_to_frequency_weights, external_ierr)
155 214 : IF (external_ierr == 0) THEN
156 78 : max_fit_error = 0.0_dp
157 78 : IF (PRESENT(used_external_backend)) used_external_backend = .TRUE.
158 78 : CALL timestop(handle)
159 : RETURN
160 : END IF
161 136 : CALL time_frequency_grid_release(grid)
162 : END IF
163 :
164 260 : CPASSERT(num_points > 0)
165 260 : CPASSERT(energy_min > 0.0_dp)
166 260 : CPASSERT(energy_max >= energy_min)
167 260 : CPASSERT(time_scaling > 0.0_dp)
168 260 : CPASSERT(time_weight_scaling > 0.0_dp)
169 260 : CPASSERT(.NOT. build_transforms .OR. (build_frequency .AND. build_time))
170 260 : CPASSERT(.NOT. build_sine .OR. build_transforms)
171 260 : max_error = 0.0_dp
172 :
173 260 : e_range = energy_max/energy_min
174 780 : ALLOCATE (points_and_weights(2*num_points))
175 :
176 260 : IF (build_frequency) THEN
177 196 : IF (num_points <= 20) THEN
178 196 : CALL get_rpa_minimax_coeff(num_points, e_range, points_and_weights, ierr, print_warning)
179 : ELSE
180 0 : CALL get_rpa_minimax_coeff_larger_grid(num_points, e_range, points_and_weights)
181 : END IF
182 :
183 588 : ALLOCATE (grid%frequency(num_points))
184 392 : ALLOCATE (grid%frequency_weights(num_points))
185 2368 : grid%frequency(:) = points_and_weights(1:num_points)*energy_min
186 2368 : grid%frequency_weights(:) = points_and_weights(num_points + 1:)*energy_min
187 196 : IF (num_points >= 26) grid%frequency_weights(:) = grid%frequency_weights(:)*4.0_dp
188 : END IF
189 :
190 260 : IF (build_time) THEN
191 230 : IF (num_points <= 20) THEN
192 230 : CALL get_exp_minimax_coeff(num_points, e_range, points_and_weights)
193 : ELSE
194 0 : CALL get_exp_minimax_coeff_gw(num_points, e_range, points_and_weights)
195 : END IF
196 :
197 690 : ALLOCATE (grid%imaginary_time(num_points))
198 460 : ALLOCATE (grid%time_weights_at_zero_frequency(num_points))
199 2526 : grid%imaginary_time(:) = points_and_weights(1:num_points)/time_scaling/energy_min
200 2526 : grid%time_weights_at_zero_frequency(:) = points_and_weights(num_points + 1:)/time_weight_scaling/energy_min
201 :
202 230 : IF (build_transforms) THEN
203 664 : ALLOCATE (grid%cosine_time_to_frequency_weights(num_points, num_points))
204 : CALL fit_l_sq_weights(cosine_time_to_frequency, num_points, grid%imaginary_time, &
205 : grid%cosine_time_to_frequency_weights, grid%frequency, energy_min, energy_max, &
206 166 : max_error, num_points_per_magnitude, regularization)
207 :
208 498 : 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 166 : max_error, num_points_per_magnitude, regularization)
212 :
213 166 : IF (build_sine) THEN
214 402 : 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 134 : max_error, num_points_per_magnitude, regularization)
218 : END IF
219 : END IF
220 : END IF
221 :
222 260 : max_fit_error = max_error
223 :
224 260 : DEALLOCATE (points_and_weights)
225 :
226 260 : CALL timestop(handle)
227 :
228 338 : END SUBROUTINE build_minimax_time_frequency_grid
229 :
230 : ! **************************************************************************************************
231 : !> \brief Calculate least-squares weights for a time/frequency transform.
232 : !> \param transform_kind Type of transform to fit.
233 : !> \param num_integ_points Number of integration points.
234 : !> \param tau_tj Imaginary-time integration points.
235 : !> \param weights Transform weights to construct.
236 : !> \param omega_tj Frequency integration points.
237 : !> \param E_min Lower end of the fitting interval.
238 : !> \param E_max Upper end of the fitting interval.
239 : !> \param max_error Maximum fitting error.
240 : !> \param num_points_per_magnitude Number of fitting points per decade.
241 : !> \param regularization Regularization parameter for the pseudoinverse.
242 : ! **************************************************************************************************
243 466 : 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)
245 :
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
254 :
255 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fit_l_sq_weights'
256 :
257 : INTEGER :: handle, iii, info, jjj, jquad, lwork, &
258 : num_x_nodes
259 466 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
260 : LOGICAL :: sine_transform, time_to_frequency
261 : REAL(KIND=dp) :: func_val, max_error_tmp, multiplicator, &
262 : omega, tau, x_value
263 466 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sing_values, vec_UTy, weight_work, work, &
264 466 : x_values, y_values
265 466 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mat_A, mat_SinvVSinvSigma, &
266 466 : mat_SinvVSinvT, mat_U
267 :
268 466 : CALL timeset(routineN, handle)
269 :
270 466 : CPASSERT(transform_kind >= cosine_time_to_frequency)
271 466 : CPASSERT(transform_kind <= cosine_frequency_to_time)
272 466 : CPASSERT(SIZE(tau_tj) == num_integ_points)
273 466 : CPASSERT(SIZE(omega_tj) == num_integ_points)
274 466 : CPASSERT(SIZE(weights, 1) == num_integ_points)
275 466 : CPASSERT(SIZE(weights, 2) == num_integ_points)
276 :
277 466 : time_to_frequency = transform_kind /= cosine_frequency_to_time
278 466 : sine_transform = transform_kind == sine_time_to_frequency
279 :
280 : ! take num_points_per_magnitude points per 10-interval
281 466 : num_x_nodes = (INT(LOG10(E_max/E_min)) + 1)*num_points_per_magnitude
282 :
283 : ! take at least as many x points as integration points to have clear
284 : ! input for the singular value decomposition
285 466 : num_x_nodes = MAX(num_x_nodes, num_integ_points)
286 :
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 10252 : mat_SinvVSinvSigma(num_integ_points, num_x_nodes), vec_UTy(num_x_nodes))
292 466 : mat_SinvVSinvSigma = 0.0_dp
293 :
294 : ! double the value nessary for 'A' to achieve good performance
295 466 : lwork = SIZE(work)
296 :
297 466 : multiplicator = (E_max/E_min)**(1.0_dp/(REAL(num_x_nodes, KIND=dp) - 1.0_dp))
298 140466 : DO iii = 1, num_x_nodes
299 140466 : x_values(iii) = E_min*multiplicator**(iii - 1)
300 : END DO
301 :
302 466 : max_error = 0.0_dp
303 :
304 6486 : DO jquad = 1, num_integ_points
305 :
306 6020 : IF (time_to_frequency) THEN
307 3984 : omega = omega_tj(jquad)
308 :
309 3984 : IF (sine_transform) THEN
310 : ! y=2*omega/(x^2+omega^2)
311 493148 : DO iii = 1, num_x_nodes
312 493148 : y_values(iii) = 2.0_dp*omega/(x_values(iii)**2 + omega**2)
313 : END DO
314 : ELSE
315 : ! y=2*x/(x^2+omega^2)
316 524436 : DO iii = 1, num_x_nodes
317 524436 : y_values(iii) = 2.0_dp*x_values(iii)/(x_values(iii)**2 + omega**2)
318 : END DO
319 : END IF
320 :
321 70928 : DO jjj = 1, num_integ_points
322 15318128 : DO iii = 1, num_x_nodes
323 15314144 : IF (sine_transform) THEN
324 7579200 : mat_A(iii, jjj) = SIN(omega*tau_tj(jjj))*EXP(-x_values(iii)*tau_tj(jjj))
325 : ELSE
326 7668000 : mat_A(iii, jjj) = COS(omega*tau_tj(jjj))*EXP(-x_values(iii)*tau_tj(jjj))
327 : END IF
328 : END DO
329 : END DO
330 : ELSE
331 2036 : tau = tau_tj(jquad)
332 :
333 : ! y=exp(-x*|tau|)
334 524436 : DO iii = 1, num_x_nodes
335 524436 : y_values(iii) = EXP(-x_values(iii)*tau)
336 : END DO
337 :
338 35632 : DO jjj = 1, num_integ_points
339 33596 : omega = omega_tj(jjj)
340 7703632 : DO iii = 1, num_x_nodes
341 7668000 : x_value = x_values(iii)
342 7701596 : mat_A(iii, jjj) = COS(tau*omega)*2.0_dp*x_value/(x_value**2 + omega**2)
343 : END DO
344 : END DO
345 : END IF
346 :
347 : ! Singular value decomposition of mat_A
348 : CALL DGESDD('A', num_x_nodes, num_integ_points, mat_A, num_x_nodes, sing_values, mat_U, num_x_nodes, &
349 6020 : mat_SinvVSinvT, num_x_nodes, work, lwork, iwork, info)
350 6020 : CPASSERT(info == 0)
351 :
352 : ! integration weights = V Sigma U^T y
353 106560 : DO jjj = 1, num_integ_points
354 1967648 : DO iii = 1, num_integ_points
355 : mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)*sing_values(jjj) &
356 1961628 : /(regularization**2 + sing_values(jjj)**2)
357 : END DO
358 : END DO
359 :
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 6020 : 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 6020 : vec_UTy, num_x_nodes, 0.0_dp, weight_work, num_integ_points)
364 :
365 106560 : weights(jquad, :) = weight_work(:)
366 :
367 6020 : max_error_tmp = 0.0_dp
368 1542020 : DO iii = 1, num_x_nodes
369 1536000 : func_val = 0.0_dp
370 1536000 : IF (time_to_frequency) THEN
371 16260800 : DO jjj = 1, num_integ_points
372 16260800 : IF (sine_transform) THEN
373 7579200 : func_val = func_val + weight_work(jjj)*SIN(omega*tau_tj(jjj))*EXP(-x_values(iii)*tau_tj(jjj))
374 : ELSE
375 7668000 : func_val = func_val + weight_work(jjj)*COS(omega*tau_tj(jjj))*EXP(-x_values(iii)*tau_tj(jjj))
376 : END IF
377 : END DO
378 : ELSE
379 8190400 : DO jjj = 1, num_integ_points
380 7668000 : omega = omega_tj(jjj)
381 : func_val = func_val + weight_work(jjj)*COS(tau*omega)*2.0_dp*x_values(iii) &
382 8190400 : /(x_values(iii)**2 + omega**2)
383 : END DO
384 : END IF
385 1542020 : max_error_tmp = MAX(max_error_tmp, ABS(y_values(iii) - func_val))
386 : END DO
387 6486 : max_error = MAX(max_error, max_error_tmp)
388 :
389 : END DO
390 :
391 0 : DEALLOCATE (x_values, y_values, mat_A, weight_work, sing_values, mat_U, mat_SinvVSinvT, work, iwork, &
392 466 : mat_SinvVSinvSigma, vec_UTy)
393 :
394 466 : CALL timestop(handle)
395 :
396 466 : END SUBROUTINE fit_l_sq_weights
397 :
398 : ! **************************************************************************************************
399 : !> \brief test the singular value decomposition for the computation of integration weights for the
400 : !> Fourier transform between time and frequency grid in cubic-scaling RPA
401 : !> \param nR ...
402 : !> \param iw ...
403 : ! **************************************************************************************************
404 0 : SUBROUTINE test_least_square_ft(nR, iw)
405 : INTEGER, INTENT(IN) :: nR, iw
406 :
407 : REAL(KIND=dp), PARAMETER :: Rc_max = 1.0E+07_dp
408 :
409 : INTEGER :: ierr, iR, num_integ_points
410 : REAL(KIND=dp) :: max_error, multiplicator, Rc
411 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tau_tj, tj, x_tw
412 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: weights_cos_tf_t_to_w
413 :
414 0 : multiplicator = Rc_max**(1.0_dp/(REAL(nR, KIND=dp) - 1.0_dp))
415 :
416 0 : DO num_integ_points = 1, 20
417 :
418 0 : ALLOCATE (x_tw(2*num_integ_points))
419 0 : x_tw = 0.0_dp
420 0 : ALLOCATE (tau_tj(num_integ_points))
421 0 : tau_tj = 0.0_dp
422 0 : ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
423 : weights_cos_tf_t_to_w = 0.0_dp
424 0 : ALLOCATE (tj(num_integ_points))
425 0 : tj = 0.0_dp
426 :
427 0 : DO iR = 0, nR - 1
428 :
429 0 : Rc = 2.0_dp*multiplicator**iR
430 :
431 0 : ierr = 0
432 0 : CALL get_rpa_minimax_coeff(num_integ_points, Rc, x_tw, ierr, print_warning=.FALSE.)
433 :
434 0 : tj(:) = x_tw(1:num_integ_points)
435 :
436 0 : x_tw = 0.0_dp
437 :
438 0 : CALL get_exp_minimax_coeff(num_integ_points, Rc, x_tw)
439 :
440 0 : tau_tj(:) = x_tw(1:num_integ_points)/2.0_dp
441 :
442 : CALL fit_l_sq_weights(cosine_time_to_frequency, num_integ_points, tau_tj, &
443 0 : weights_cos_tf_t_to_w, tj, 1.0_dp, Rc, max_error, 200, 0.0_dp)
444 :
445 0 : IF (iw > 0) THEN
446 0 : WRITE (iw, '(T2, I3, F12.1, ES12.3)') num_integ_points, Rc, max_error
447 : END IF
448 :
449 : END DO
450 :
451 0 : DEALLOCATE (x_tw, tau_tj, weights_cos_tf_t_to_w, tj)
452 :
453 : END DO
454 :
455 0 : END SUBROUTINE test_least_square_ft
456 :
457 0 : END MODULE time_frequency_grids
|
