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 : END DO
78 530 : DEALLOCATE (ch, cg, cgt, cht)
79 530 : DEALLOCATE (y)
80 530 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: y
97 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht
98 :
99 : !must be 2**nex
100 :
101 0 : a = 0.0_dp
102 0 : x = 0.0_dp
103 0 : m = itype + 2
104 0 : ni = 2*itype
105 0 : CALL lazy_arrays(itype, m, ch, cg, cgt, cht)
106 0 : ALLOCATE (y(0:nd), stat=i_all)
107 0 : IF (i_all /= 0) THEN
108 0 : CPABORT("Wavelet_function: problem of memory allocation")
109 : END IF
110 :
111 : ! plot wavelet
112 0 : CALL zero(nd + 1, x)
113 0 : CALL zero(nd + 1, y)
114 0 : nt = ni
115 0 : x(nt + nt/2 - 1) = 1._dp
116 : loop3: DO
117 0 : nt = 2*nt
118 : !WRITE(*,*) 'nd,nt',nd,nt
119 0 : CALL back_trans(nd, nt, x, y, m, ch, cg)
120 0 : CALL dcopy(nd, y, 1, x, 1)
121 0 : IF (nt == nd) THEN
122 : EXIT loop3
123 : END IF
124 : END DO loop3
125 :
126 : !open (unit=1,file='wavelet',status='unknown')
127 0 : DO i = 0, nd - 1
128 0 : a(i) = 1._dp*i*ni/nd - (.5_dp*ni - .5_dp)
129 : END DO
130 0 : DEALLOCATE (ch, cg, cgt, cht)
131 0 : DEALLOCATE (y)
132 :
133 0 : 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 48985 : 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 48985 : REAL(KIND=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht
151 :
152 9097120 : kern_1_scf = 0.0_dp
153 48985 : m = itype + 2
154 48985 : CALL lazy_arrays(itype, m, ch, cg, cgt, cht)
155 48985 : CALL scf_recurs(n_iter, n_range, kernel_scf, kern_1_scf, m, ch)
156 48985 : DEALLOCATE (ch, cg, cgt, cht)
157 :
158 48985 : END SUBROUTINE scf_recursion
159 :
160 : ! **************************************************************************************************
161 : !> \brief Set to zero an array x(n)
162 : !> \param n ...
163 : !> \param x ...
164 : ! **************************************************************************************************
165 1060 : PURE SUBROUTINE zero(n, x)
166 : INTEGER, INTENT(in) :: n
167 : REAL(KIND=dp), INTENT(out) :: x(n)
168 :
169 : INTEGER :: i
170 :
171 5362760 : DO i = 1, n
172 5362760 : x(i) = 0._dp
173 : END DO
174 1060 : 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 0 : 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 0 : y = 0.0_dp
200 0 : DO i = 0, nt/2 - 1
201 0 : y(i) = 0._dp
202 0 : y(nt/2 + i) = 0._dp
203 :
204 0 : DO j = -m + 1, m
205 :
206 : ! periodically wrap index if necessary
207 0 : ind = j + 2*i
208 : loop99: DO
209 0 : IF (ind < 0) THEN
210 0 : ind = ind + nt
211 0 : CYCLE loop99
212 : END IF
213 0 : IF (ind >= nt) THEN
214 0 : ind = ind - nt
215 0 : CYCLE loop99
216 : END IF
217 : EXIT loop99
218 : END DO loop99
219 :
220 0 : y(i) = y(i) + cht(j)*x(ind)
221 0 : y(nt/2 + i) = y(nt/2 + i) + cgt(j)*x(ind)
222 : END DO
223 :
224 : END DO
225 :
226 0 : 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 3180 : 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 16085100 : y = 0.0_dp
253 :
254 2641620 : DO i = 0, nt/2 - 1
255 2638440 : y(2*i + 0) = 0._dp
256 2638440 : y(2*i + 1) = 0._dp
257 :
258 112997460 : DO j = -m/2, m/2 - 1
259 :
260 : ! periodically wrap index if necessary
261 110355840 : ind = i - j
262 : loop99: DO
263 111735600 : IF (ind < 0) THEN
264 656880 : ind = ind + nt/2
265 656880 : CYCLE loop99
266 : END IF
267 111078720 : IF (ind >= nt/2) THEN
268 722880 : ind = ind - nt/2
269 722880 : CYCLE loop99
270 : END IF
271 : EXIT loop99
272 : END DO loop99
273 :
274 110355840 : y(2*i + 0) = y(2*i + 0) + ch(2*j - 0)*x(ind) + cg(2*j - 0)*x(ind + nt/2)
275 112994280 : 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 3180 : 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 0 : 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 0 : DO i = -m, m
304 0 : DO j = -m, m
305 0 : t1 = 0._dp
306 0 : t2 = 0._dp
307 0 : t3 = 0._dp
308 0 : t4 = 0._dp
309 0 : DO l = -3*m, 3*m
310 : IF (l - 2*i >= -m .AND. l - 2*i <= m .AND. &
311 0 : l - 2*j >= -m .AND. l - 2*j <= m) THEN
312 0 : t1 = t1 + ch(l - 2*i)*cht(l - 2*j)
313 0 : t2 = t2 + cg(l - 2*i)*cgt(l - 2*j)
314 0 : t3 = t3 + ch(l - 2*i)*cgt(l - 2*j)
315 0 : t4 = t4 + cht(l - 2*i)*cg(l - 2*j)
316 : END IF
317 : END DO
318 0 : eps = 1.e-10_dp
319 0 : IF (i == j) THEN
320 : IF (ABS(t1 - 1._dp) > eps .OR. ABS(t2 - 1._dp) > eps .OR. &
321 0 : ABS(t3) > eps .OR. ABS(t4) > eps) THEN
322 0 : 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 0 : ABS(t3) > eps .OR. ABS(t4) > eps) THEN
327 0 : WRITE (*, fmt22) 'Orthogonality ERROR', i, j, t1, t2, t3, t4
328 : END IF
329 : END IF
330 : END DO
331 : END DO
332 :
333 0 : WRITE (*, *) 'FILTER TEST PASSED'
334 :
335 0 : 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 48985 : 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 9097120 : kern_1_scf = 0.0_dp
358 : !Start the iteration to go from p0gauss to pgauss
359 546697 : loop_iter_scf: DO i_iter = 1, n_iter
360 85471416 : kern_1_scf(:) = kernel_scf(:)
361 85471416 : kernel_scf(:) = 0._dp
362 19930365 : loop_iter_i: DO i = 0, n_range
363 19881380 : kern_tot = 0._dp
364 1664764600 : DO j = -m, m
365 1644883220 : ind = 2*i - j
366 1644883220 : IF (ABS(ind) > n_range) THEN
367 : kern = 0._dp
368 : ELSE
369 1456043460 : kern = kern_1_scf(ind)
370 : END IF
371 1664764600 : kern_tot = kern_tot + ch(j)*kern
372 : END DO
373 19881380 : IF (kern_tot == 0._dp) THEN
374 : !zero after (be sure because strictly == 0._dp)
375 : EXIT loop_iter_i
376 : ELSE
377 19383668 : kernel_scf(i) = 0.5_dp*kern_tot
378 19383668 : kernel_scf(-i) = kernel_scf(i)
379 : END IF
380 : END DO loop_iter_i
381 : END DO loop_iter_scf
382 48985 : END SUBROUTINE scf_recurs
383 :
384 : END MODULE ps_wavelet_scaling_function
|
