(git:58e3e09)
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 ! **************************************************************************************************
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, &
23 
24 CONTAINS
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 
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Filters for interpolating scaling functions .
Definition: lazy.F:12
subroutine, public lazy_arrays(itype, m, ch, cg, cgt, cht)
...
Definition: lazy.F:42
Creates the wavelet kernel for the wavelet based poisson solver.
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.
subroutine, public scaling_function(itype, nd, nrange, a, x)
Calculate the values of a scaling function in real uniform grid.