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