(git:98357aa)
Loading...
Searching...
No Matches
ps_wavelet_scaling_function.F
Go to the documentation of this file.
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 Creates the wavelet kernel for the wavelet based poisson solver.
10!> \author Florian Schiffmann (09.2007,fschiff)
11! **************************************************************************************************
12MODULE ps_wavelet_scaling_function
13 USE kinds, ONLY: dp
14 USE lazy, ONLY: lazy_arrays
15#include "../base/base_uses.f90"
16
17 IMPLICIT NONE
18
19 PRIVATE
20
21 PUBLIC :: scaling_function, &
22 scf_recursion
23
24CONTAINS
25
26! **************************************************************************************************
27!> \brief Calculate the values of a scaling function in real uniform grid
28!> \param itype ...
29!> \param nd ...
30!> \param nrange ...
31!> \param a ...
32!> \param x ...
33! **************************************************************************************************
34 SUBROUTINE scaling_function(itype, nd, nrange, a, x)
35
36 !Type of interpolating functions
37 INTEGER, INTENT(in) :: itype, nd
38 INTEGER, INTENT(out) :: nrange
39 REAL(kind=dp), DIMENSION(0:nd), INTENT(out) :: a, x
40
41 INTEGER :: i, i_all, m, ni, nt
42 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: y
43 REAL(kind=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht
44
45!Number of points: must be 2**nex
46
47 a = 0.0_dp
48 x = 0.0_dp
49 m = itype + 2
50 CALL lazy_arrays(itype, m, ch, cg, cgt, cht)
51
52 ni = 2*itype
53 nrange = ni
54 ALLOCATE (y(0:nd), stat=i_all)
55 IF (i_all /= 0) THEN
56 cpabort("Scaling_function: problem of memory allocation")
57 END IF
58
59 ! plot scaling function
60 CALL zero(nd + 1, x)
61 CALL zero(nd + 1, y)
62 nt = ni
63 x(nt/2 - 1) = 1._dp
64 loop1: DO
65 nt = 2*nt
66
67 CALL back_trans(nd, nt, x, y, m, ch, cg)
68 CALL dcopy(nt, y, 1, x, 1)
69 IF (nt == nd) THEN
70 EXIT loop1
71 END IF
72 END DO loop1
73
74 !open (unit=1,file='scfunction',status='unknown')
75 DO i = 0, nd
76 a(i) = 1._dp*i*ni/nd - (.5_dp*ni - 1._dp)
77 END DO
78 DEALLOCATE (ch, cg, cgt, cht)
79 DEALLOCATE (y)
80 END SUBROUTINE scaling_function
81
82! **************************************************************************************************
83!> \brief Calculate the values of the wavelet function in a real uniform mesh.
84!> \param itype ...
85!> \param nd ...
86!> \param a ...
87!> \param x ...
88! **************************************************************************************************
89 SUBROUTINE wavelet_function(itype, nd, a, x)
90
91 !Type of the interpolating scaling function
92 INTEGER, INTENT(in) :: itype, nd
93 REAL(kind=dp), DIMENSION(0:nd), INTENT(out) :: a, x
94
95 INTEGER :: i, i_all, m, ni, nt
96 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: y
97 REAL(kind=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht
98
99!must be 2**nex
100
101 a = 0.0_dp
102 x = 0.0_dp
103 m = itype + 2
104 ni = 2*itype
105 CALL lazy_arrays(itype, m, ch, cg, cgt, cht)
106 ALLOCATE (y(0:nd), stat=i_all)
107 IF (i_all /= 0) THEN
108 cpabort("Wavelet_function: problem of memory allocation")
109 END IF
110
111 ! plot wavelet
112 CALL zero(nd + 1, x)
113 CALL zero(nd + 1, y)
114 nt = ni
115 x(nt + nt/2 - 1) = 1._dp
116 loop3: DO
117 nt = 2*nt
118 !WRITE(*,*) 'nd,nt',nd,nt
119 CALL back_trans(nd, nt, x, y, m, ch, cg)
120 CALL dcopy(nd, y, 1, x, 1)
121 IF (nt == nd) THEN
122 EXIT loop3
123 END IF
124 END DO loop3
125
126 !open (unit=1,file='wavelet',status='unknown')
127 DO i = 0, nd - 1
128 a(i) = 1._dp*i*ni/nd - (.5_dp*ni - .5_dp)
129 END DO
130 DEALLOCATE (ch, cg, cgt, cht)
131 DEALLOCATE (y)
132
133 END SUBROUTINE wavelet_function
134
135! **************************************************************************************************
136!> \brief Do iterations to go from p0gauss to pgauss
137!> order interpolating scaling function
138!> \param itype ...
139!> \param n_iter ...
140!> \param n_range ...
141!> \param kernel_scf ...
142!> \param kern_1_scf ...
143! **************************************************************************************************
144 SUBROUTINE scf_recursion(itype, n_iter, n_range, kernel_scf, kern_1_scf)
145 INTEGER, INTENT(in) :: itype, n_iter, n_range
146 REAL(kind=dp), INTENT(inout) :: kernel_scf(-n_range:n_range)
147 REAL(kind=dp), INTENT(out) :: kern_1_scf(-n_range:n_range)
148
149 INTEGER :: m
150 REAL(kind=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht
151
152 kern_1_scf = 0.0_dp
153 m = itype + 2
154 CALL lazy_arrays(itype, m, ch, cg, cgt, cht)
155 CALL scf_recurs(n_iter, n_range, kernel_scf, kern_1_scf, m, ch)
156 DEALLOCATE (ch, cg, cgt, cht)
157
158 END SUBROUTINE scf_recursion
159
160! **************************************************************************************************
161!> \brief Set to zero an array x(n)
162!> \param n ...
163!> \param x ...
164! **************************************************************************************************
165 PURE SUBROUTINE zero(n, x)
166 INTEGER, INTENT(in) :: n
167 REAL(kind=dp), INTENT(out) :: x(n)
168
169 INTEGER :: i
170
171 DO i = 1, n
172 x(i) = 0._dp
173 END DO
174 END SUBROUTINE zero
175
176! **************************************************************************************************
177!> \brief forward wavelet transform
178!> nd: length of data set
179!> nt length of data in data set to be transformed
180!> m filter length (m has to be even!)
181!> x input data, y output data
182!> \param nd ...
183!> \param nt ...
184!> \param x ...
185!> \param y ...
186!> \param m ...
187!> \param cgt ...
188!> \param cht ...
189! **************************************************************************************************
190 SUBROUTINE for_trans(nd, nt, x, y, m, cgt, cht)
191 INTEGER, INTENT(in) :: nd, nt
192 REAL(kind=dp), INTENT(in) :: x(0:nd - 1)
193 REAL(kind=dp), INTENT(out) :: y(0:nd - 1)
194 INTEGER :: m
195 REAL(kind=dp), DIMENSION(:), POINTER :: cgt, cht
196
197 INTEGER :: i, ind, j
198
199 y = 0.0_dp
200 DO i = 0, nt/2 - 1
201 y(i) = 0._dp
202 y(nt/2 + i) = 0._dp
203
204 DO j = -m + 1, m
205
206 ! periodically wrap index if necessary
207 ind = j + 2*i
208 loop99: DO
209 IF (ind < 0) THEN
210 ind = ind + nt
211 cycle loop99
212 END IF
213 IF (ind >= nt) THEN
214 ind = ind - nt
215 cycle loop99
216 END IF
217 EXIT loop99
218 END DO loop99
219
220 y(i) = y(i) + cht(j)*x(ind)
221 y(nt/2 + i) = y(nt/2 + i) + cgt(j)*x(ind)
222 END DO
223
224 END DO
225
226 END SUBROUTINE for_trans
227
228! **************************************************************************************************
229!> \brief ...
230!> \param nd ...
231!> \param nt ...
232!> \param x ...
233!> \param y ...
234!> \param m ...
235!> \param ch ...
236!> \param cg ...
237! **************************************************************************************************
238 SUBROUTINE back_trans(nd, nt, x, y, m, ch, cg)
239 ! backward wavelet transform
240 ! nd: length of data set
241 ! nt length of data in data set to be transformed
242 ! m filter length (m has to be even!)
243 ! x input data, y output data
244 INTEGER, INTENT(in) :: nd, nt
245 REAL(kind=dp), INTENT(in) :: x(0:nd - 1)
246 REAL(kind=dp), INTENT(out) :: y(0:nd - 1)
247 INTEGER :: m
248 REAL(kind=dp), DIMENSION(:), POINTER :: ch, cg
249
250 INTEGER :: i, ind, j
251
252 y = 0.0_dp
253
254 DO i = 0, nt/2 - 1
255 y(2*i + 0) = 0._dp
256 y(2*i + 1) = 0._dp
257
258 DO j = -m/2, m/2 - 1
259
260 ! periodically wrap index if necessary
261 ind = i - j
262 loop99: DO
263 IF (ind < 0) THEN
264 ind = ind + nt/2
265 cycle loop99
266 END IF
267 IF (ind >= nt/2) THEN
268 ind = ind - nt/2
269 cycle loop99
270 END IF
271 EXIT loop99
272 END DO loop99
273
274 y(2*i + 0) = y(2*i + 0) + ch(2*j - 0)*x(ind) + cg(2*j - 0)*x(ind + nt/2)
275 y(2*i + 1) = y(2*i + 1) + ch(2*j + 1)*x(ind) + cg(2*j + 1)*x(ind + nt/2)
276 END DO
277
278 END DO
279
280 END SUBROUTINE back_trans
281
282! **************************************************************************************************
283!> \brief Tests the 4 orthogonality relations of the filters
284!> \param m ...
285!> \param ch ...
286!> \param cg ...
287!> \param cgt ...
288!> \param cht ...
289! **************************************************************************************************
290 SUBROUTINE ftest(m, ch, cg, cgt, cht)
291 INTEGER :: m
292 REAL(kind=dp), DIMENSION(:), POINTER :: ch, cg, cgt, cht
293
294 CHARACTER(len=*), PARAMETER :: fmt22 = "(a,i3,i4,4(e17.10))"
295
296 INTEGER :: i, j, l
297 REAL(kind=dp) :: eps, t1, t2, t3, t4
298
299! do i=-m,m
300! WRITE(*,*) i,ch(i),cg(i)
301! end do
302
303 DO i = -m, m
304 DO j = -m, m
305 t1 = 0._dp
306 t2 = 0._dp
307 t3 = 0._dp
308 t4 = 0._dp
309 DO l = -3*m, 3*m
310 IF (l - 2*i >= -m .AND. l - 2*i <= m .AND. &
311 l - 2*j >= -m .AND. l - 2*j <= m) THEN
312 t1 = t1 + ch(l - 2*i)*cht(l - 2*j)
313 t2 = t2 + cg(l - 2*i)*cgt(l - 2*j)
314 t3 = t3 + ch(l - 2*i)*cgt(l - 2*j)
315 t4 = t4 + cht(l - 2*i)*cg(l - 2*j)
316 END IF
317 END DO
318 eps = 1.e-10_dp
319 IF (i == j) THEN
320 IF (abs(t1 - 1._dp) > eps .OR. abs(t2 - 1._dp) > eps .OR. &
321 abs(t3) > eps .OR. abs(t4) > eps) THEN
322 WRITE (*, fmt22) 'Orthogonality ERROR', i, j, t1, t2, t3, t4
323 END IF
324 ELSE
325 IF (abs(t1) > eps .OR. abs(t2) > eps .OR. &
326 abs(t3) > eps .OR. abs(t4) > eps) THEN
327 WRITE (*, fmt22) 'Orthogonality ERROR', i, j, t1, t2, t3, t4
328 END IF
329 END IF
330 END DO
331 END DO
332
333 WRITE (*, *) 'FILTER TEST PASSED'
334
335 END SUBROUTINE ftest
336
337! **************************************************************************************************
338!> \brief Do iterations to go from p0gauss to pgauss
339!> 8th-order interpolating scaling function
340!> \param n_iter ...
341!> \param n_range ...
342!> \param kernel_scf ...
343!> \param kern_1_scf ...
344!> \param m ...
345!> \param ch ...
346! **************************************************************************************************
347 SUBROUTINE scf_recurs(n_iter, n_range, kernel_scf, kern_1_scf, m, ch)
348 INTEGER, INTENT(in) :: n_iter, n_range
349 REAL(kind=dp), INTENT(inout) :: kernel_scf(-n_range:n_range)
350 REAL(kind=dp), INTENT(out) :: kern_1_scf(-n_range:n_range)
351 INTEGER :: m
352 REAL(kind=dp), DIMENSION(:), POINTER :: ch
353
354 INTEGER :: i, i_iter, ind, j
355 REAL(kind=dp) :: kern, kern_tot
356
357 kern_1_scf = 0.0_dp
358 !Start the iteration to go from p0gauss to pgauss
359 loop_iter_scf: DO i_iter = 1, n_iter
360 kern_1_scf(:) = kernel_scf(:)
361 kernel_scf(:) = 0._dp
362 loop_iter_i: DO i = 0, n_range
363 kern_tot = 0._dp
364 DO j = -m, m
365 ind = 2*i - j
366 IF (abs(ind) > n_range) THEN
367 kern = 0._dp
368 ELSE
369 kern = kern_1_scf(ind)
370 END IF
371 kern_tot = kern_tot + ch(j)*kern
372 END DO
373 IF (kern_tot == 0._dp) THEN
374 !zero after (be sure because strictly == 0._dp)
375 EXIT loop_iter_i
376 ELSE
377 kernel_scf(i) = 0.5_dp*kern_tot
378 kernel_scf(-i) = kernel_scf(i)
379 END IF
380 END DO loop_iter_i
381 END DO loop_iter_scf
382 END SUBROUTINE scf_recurs
383
384END MODULE ps_wavelet_scaling_function
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
lazy
Filters for interpolating scaling functions .
Definition lazy.F:12
lazy::lazy_arrays
subroutine, public lazy_arrays(itype, m, ch, cg, cgt, cht)
...
Definition lazy.F:42
ps_wavelet_scaling_function
Creates the wavelet kernel for the wavelet based poisson solver.
Definition ps_wavelet_scaling_function.F:12
ps_wavelet_scaling_function::scf_recursion
subroutine, public scf_recursion(itype, n_iter, n_range, kernel_scf, kern_1_scf)
Do iterations to go from p0gauss to pgauss order interpolating scaling function.
Definition ps_wavelet_scaling_function.F:145
ps_wavelet_scaling_function::scaling_function
subroutine, public scaling_function(itype, nd, nrange, a, x)
Calculate the values of a scaling function in real uniform grid.
Definition ps_wavelet_scaling_function.F:35