(git:92574dc)
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 DO i = 0, nd
75 a(i) = 1._dp*i*ni/nd - (.5_dp*ni - 1._dp)
76 END DO
77 DEALLOCATE (ch, cg, cgt, cht)
78 DEALLOCATE (y)
79 END SUBROUTINE scaling_function
80
81! **************************************************************************************************
82!> \brief Calculate the values of the wavelet function in a real uniform mesh.
83!> \param itype ...
84!> \param nd ...
85!> \param a ...
86!> \param x ...
87! **************************************************************************************************
88 SUBROUTINE wavelet_function(itype, nd, a, x)
89
90 !Type of the interpolating scaling function
91 INTEGER, INTENT(in) :: itype, nd
92 REAL(kind=dp), DIMENSION(0:nd), INTENT(out) :: a, x
93
94 INTEGER :: i, i_all, m, ni, nt
95 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: y
96 REAL(kind=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht
97
98!must be 2**nex
99
100 a = 0.0_dp
101 x = 0.0_dp
102 m = itype + 2
103 ni = 2*itype
104 CALL lazy_arrays(itype, m, ch, cg, cgt, cht)
105 ALLOCATE (y(0:nd), stat=i_all)
106 IF (i_all /= 0) THEN
107 cpabort("Wavelet_function: problem of memory allocation")
108 END IF
109
110 ! plot wavelet
111 CALL zero(nd + 1, x)
112 CALL zero(nd + 1, y)
113 nt = ni
114 x(nt + nt/2 - 1) = 1._dp
115 loop3: DO
116 nt = 2*nt
117 !WRITE(*,*) 'nd,nt',nd,nt
118 CALL back_trans(nd, nt, x, y, m, ch, cg)
119 CALL dcopy(nd, y, 1, x, 1)
120 IF (nt == nd) THEN
121 EXIT loop3
122 END IF
123 END DO loop3
124
125 !open (unit=1,file='wavelet',status='unknown')
126 DO i = 0, nd - 1
127 a(i) = 1._dp*i*ni/nd - (.5_dp*ni - .5_dp)
128 END DO
129 DEALLOCATE (ch, cg, cgt, cht)
130 DEALLOCATE (y)
131
132 END SUBROUTINE wavelet_function
133
134! **************************************************************************************************
135!> \brief Do iterations to go from p0gauss to pgauss
136!> order interpolating scaling function
137!> \param itype ...
138!> \param n_iter ...
139!> \param n_range ...
140!> \param kernel_scf ...
141!> \param kern_1_scf ...
142! **************************************************************************************************
143 PURE SUBROUTINE scf_recursion(itype, n_iter, n_range, kernel_scf, kern_1_scf)
144 INTEGER, INTENT(in) :: itype, n_iter, n_range
145 REAL(kind=dp), INTENT(inout) :: kernel_scf(-n_range:n_range)
146 REAL(kind=dp), INTENT(out) :: kern_1_scf(-n_range:n_range)
147
148 INTEGER :: m
149 REAL(kind=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht
150
151 kern_1_scf = 0.0_dp
152 m = itype + 2
153 CALL lazy_arrays(itype, m, ch, cg, cgt, cht)
154 CALL scf_recurs(n_iter, n_range, kernel_scf, kern_1_scf, m, ch)
155 DEALLOCATE (ch, cg, cgt, cht)
156
157 END SUBROUTINE scf_recursion
158
159! **************************************************************************************************
160!> \brief Set to zero an array x(n)
161!> \param n ...
162!> \param x ...
163! **************************************************************************************************
164 PURE SUBROUTINE zero(n, x)
165 INTEGER, INTENT(in) :: n
166 REAL(kind=dp), INTENT(out) :: x(n)
167
168 INTEGER :: i
169
170 DO i = 1, n
171 x(i) = 0._dp
172 END DO
173 END SUBROUTINE zero
174
175! **************************************************************************************************
176!> \brief forward wavelet transform
177!> nd: length of data set
178!> nt length of data in data set to be transformed
179!> m filter length (m has to be even!)
180!> x input data, y output data
181!> \param nd ...
182!> \param nt ...
183!> \param x ...
184!> \param y ...
185!> \param m ...
186!> \param cgt ...
187!> \param cht ...
188! **************************************************************************************************
189 PURE SUBROUTINE for_trans(nd, nt, x, y, m, cgt, cht)
190 INTEGER, INTENT(in) :: nd, nt
191 REAL(kind=dp), INTENT(in) :: x(0:nd - 1)
192 REAL(kind=dp), INTENT(out) :: y(0:nd - 1)
193 INTEGER, INTENT(in) :: m
194 REAL(kind=dp), DIMENSION(:), POINTER :: cgt, cht
195
196 INTEGER :: i, ind, j
197
198 y = 0.0_dp
199 DO i = 0, nt/2 - 1
200 y(i) = 0._dp
201 y(nt/2 + i) = 0._dp
202
203 DO j = -m + 1, m
204
205 ! periodically wrap index if necessary
206 ind = j + 2*i
207 loop99: DO
208 IF (ind < 0) THEN
209 ind = ind + nt
210 cycle loop99
211 END IF
212 IF (ind >= nt) THEN
213 ind = ind - nt
214 cycle loop99
215 END IF
216 EXIT loop99
217 END DO loop99
218
219 y(i) = y(i) + cht(j)*x(ind)
220 y(nt/2 + i) = y(nt/2 + i) + cgt(j)*x(ind)
221 END DO
222
223 END DO
224
225 END SUBROUTINE for_trans
226
227! **************************************************************************************************
228!> \brief ...
229!> \param nd ...
230!> \param nt ...
231!> \param x ...
232!> \param y ...
233!> \param m ...
234!> \param ch ...
235!> \param cg ...
236! **************************************************************************************************
237 PURE SUBROUTINE back_trans(nd, nt, x, y, m, ch, cg)
238 ! backward wavelet transform
239 ! nd: length of data set
240 ! nt length of data in data set to be transformed
241 ! m filter length (m has to be even!)
242 ! x input data, y output data
243 INTEGER, INTENT(in) :: nd, nt
244 REAL(kind=dp), INTENT(in) :: x(0:nd - 1)
245 REAL(kind=dp), INTENT(out) :: y(0:nd - 1)
246 INTEGER, INTENT(in) :: m
247 REAL(kind=dp), DIMENSION(:), POINTER :: ch, cg
248
249 INTEGER :: i, ind, j
250
251 y = 0.0_dp
252
253 DO i = 0, nt/2 - 1
254 y(2*i + 0) = 0._dp
255 y(2*i + 1) = 0._dp
256
257 DO j = -m/2, m/2 - 1
258
259 ! periodically wrap index if necessary
260 ind = i - j
261 loop99: DO
262 IF (ind < 0) THEN
263 ind = ind + nt/2
264 cycle loop99
265 END IF
266 IF (ind >= nt/2) THEN
267 ind = ind - nt/2
268 cycle loop99
269 END IF
270 EXIT loop99
271 END DO loop99
272
273 y(2*i + 0) = y(2*i + 0) + ch(2*j - 0)*x(ind) + cg(2*j - 0)*x(ind + nt/2)
274 y(2*i + 1) = y(2*i + 1) + ch(2*j + 1)*x(ind) + cg(2*j + 1)*x(ind + nt/2)
275 END DO
276
277 END DO
278
279 END SUBROUTINE back_trans
280
281! **************************************************************************************************
282!> \brief Do iterations to go from p0gauss to pgauss
283!> 8th-order interpolating scaling function
284!> \param n_iter ...
285!> \param n_range ...
286!> \param kernel_scf ...
287!> \param kern_1_scf ...
288!> \param m ...
289!> \param ch ...
290! **************************************************************************************************
291 PURE SUBROUTINE scf_recurs(n_iter, n_range, kernel_scf, kern_1_scf, m, ch)
292 INTEGER, INTENT(in) :: n_iter, n_range
293 REAL(kind=dp), INTENT(inout) :: kernel_scf(-n_range:n_range)
294 REAL(kind=dp), INTENT(out) :: kern_1_scf(-n_range:n_range)
295 INTEGER, INTENT(in) :: m
296 REAL(kind=dp), DIMENSION(:), POINTER :: ch
297
298 INTEGER :: i, i_iter, ind, j
299 REAL(kind=dp) :: kern, kern_tot
300
301 kern_1_scf = 0.0_dp
302 !Start the iteration to go from p0gauss to pgauss
303 loop_iter_scf: DO i_iter = 1, n_iter
304 kern_1_scf(:) = kernel_scf(:)
305 kernel_scf(:) = 0._dp
306 loop_iter_i: DO i = 0, n_range
307 kern_tot = 0._dp
308 DO j = -m, m
309 ind = 2*i - j
310 IF (abs(ind) > n_range) THEN
311 kern = 0._dp
312 ELSE
313 kern = kern_1_scf(ind)
314 END IF
315 kern_tot = kern_tot + ch(j)*kern
316 END DO
317 IF (kern_tot == 0._dp) THEN
318 !zero after (be sure because strictly == 0._dp)
319 EXIT loop_iter_i
320 ELSE
321 kernel_scf(i) = 0.5_dp*kern_tot
322 kernel_scf(-i) = kernel_scf(i)
323 END IF
324 END DO loop_iter_i
325 END DO loop_iter_scf
326 END SUBROUTINE scf_recurs
327
328END 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
pure 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::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
ps_wavelet_scaling_function::scf_recursion
pure 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:144