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