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 600 : 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 600 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: y
43 600 : REAL(KIND=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht
44 :
45 : !Number of points: must be 2**nex
46 :
47 3039920 : a = 0.0_dp
48 3039920 : x = 0.0_dp
49 600 : m = itype + 2
50 600 : CALL lazy_arrays(itype, m, ch, cg, cgt, cht)
51 :
52 600 : ni = 2*itype
53 600 : nrange = ni
54 1800 : ALLOCATE (y(0:nd), stat=i_all)
55 600 : IF (i_all /= 0) THEN
56 0 : CPABORT("Scaling_function: problem of memory allocation")
57 : END IF
58 :
59 : ! plot scaling function
60 600 : CALL zero(nd + 1, x)
61 600 : CALL zero(nd + 1, y)
62 600 : nt = ni
63 600 : x(nt/2 - 1) = 1._dp
64 : loop1: DO
65 3600 : nt = 2*nt
66 :
67 3600 : CALL back_trans(nd, nt, x, y, m, ch, cg)
68 3600 : CALL dcopy(nt, y, 1, x, 1)
69 3600 : IF (nt == nd) THEN
70 : EXIT loop1
71 : END IF
72 : END DO loop1
73 :
74 3039920 : DO i = 0, nd
75 3039920 : a(i) = 1._dp*i*ni/nd - (.5_dp*ni - 1._dp)
76 : END DO
77 600 : DEALLOCATE (ch, cg, cgt, cht)
78 600 : DEALLOCATE (y)
79 600 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: y
96 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht
97 :
98 : !must be 2**nex
99 :
100 0 : a = 0.0_dp
101 0 : x = 0.0_dp
102 0 : m = itype + 2
103 0 : ni = 2*itype
104 0 : CALL lazy_arrays(itype, m, ch, cg, cgt, cht)
105 0 : ALLOCATE (y(0:nd), stat=i_all)
106 0 : IF (i_all /= 0) THEN
107 0 : CPABORT("Wavelet_function: problem of memory allocation")
108 : END IF
109 :
110 : ! plot wavelet
111 0 : CALL zero(nd + 1, x)
112 0 : CALL zero(nd + 1, y)
113 0 : nt = ni
114 0 : x(nt + nt/2 - 1) = 1._dp
115 : loop3: DO
116 0 : nt = 2*nt
117 : !WRITE(*,*) 'nd,nt',nd,nt
118 0 : CALL back_trans(nd, nt, x, y, m, ch, cg)
119 0 : CALL dcopy(nd, y, 1, x, 1)
120 0 : IF (nt == nd) THEN
121 : EXIT loop3
122 : END IF
123 : END DO loop3
124 :
125 : !open (unit=1,file='wavelet',status='unknown')
126 0 : DO i = 0, nd - 1
127 0 : a(i) = 1._dp*i*ni/nd - (.5_dp*ni - .5_dp)
128 : END DO
129 0 : DEALLOCATE (ch, cg, cgt, cht)
130 0 : DEALLOCATE (y)
131 :
132 0 : 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 55215 : 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 55215 : REAL(KIND=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht
150 :
151 10106380 : kern_1_scf = 0.0_dp
152 55215 : m = itype + 2
153 55215 : CALL lazy_arrays(itype, m, ch, cg, cgt, cht)
154 55215 : CALL scf_recurs(n_iter, n_range, kernel_scf, kern_1_scf, m, ch)
155 55215 : DEALLOCATE (ch, cg, cgt, cht)
156 :
157 55215 : END SUBROUTINE scf_recursion
158 :
159 : ! **************************************************************************************************
160 : !> \brief Set to zero an array x(n)
161 : !> \param n ...
162 : !> \param x ...
163 : ! **************************************************************************************************
164 1200 : PURE SUBROUTINE zero(n, x)
165 : INTEGER, INTENT(in) :: n
166 : REAL(KIND=dp), INTENT(out) :: x(n)
167 :
168 : INTEGER :: i
169 :
170 6079840 : DO i = 1, n
171 6079840 : x(i) = 0._dp
172 : END DO
173 1200 : 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 0 : 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 0 : y = 0.0_dp
199 0 : DO i = 0, nt/2 - 1
200 0 : y(i) = 0._dp
201 0 : y(nt/2 + i) = 0._dp
202 :
203 0 : DO j = -m + 1, m
204 :
205 : ! periodically wrap index if necessary
206 0 : ind = j + 2*i
207 : loop99: DO
208 0 : IF (ind < 0) THEN
209 0 : ind = ind + nt
210 0 : CYCLE loop99
211 : END IF
212 0 : IF (ind >= nt) THEN
213 0 : ind = ind - nt
214 0 : CYCLE loop99
215 : END IF
216 : EXIT loop99
217 : END DO loop99
218 :
219 0 : y(i) = y(i) + cht(j)*x(ind)
220 0 : y(nt/2 + i) = y(nt/2 + i) + cgt(j)*x(ind)
221 : END DO
222 :
223 : END DO
224 :
225 0 : 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 3600 : 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 18235920 : y = 0.0_dp
252 :
253 2994840 : DO i = 0, nt/2 - 1
254 2991240 : y(2*i + 0) = 0._dp
255 2991240 : y(2*i + 1) = 0._dp
256 :
257 128168280 : DO j = -m/2, m/2 - 1
258 :
259 : ! periodically wrap index if necessary
260 125173440 : ind = i - j
261 : loop99: DO
262 126738420 : IF (ind < 0) THEN
263 745080 : ind = ind + nt/2
264 745080 : CYCLE loop99
265 : END IF
266 125993340 : IF (ind >= nt/2) THEN
267 819900 : ind = ind - nt/2
268 819900 : CYCLE loop99
269 : END IF
270 : EXIT loop99
271 : END DO loop99
272 :
273 125173440 : y(2*i + 0) = y(2*i + 0) + ch(2*j - 0)*x(ind) + cg(2*j - 0)*x(ind + nt/2)
274 128164680 : 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 3600 : 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 55215 : 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 10106380 : kern_1_scf = 0.0_dp
302 : !Start the iteration to go from p0gauss to pgauss
303 614391 : loop_iter_scf: DO i_iter = 1, n_iter
304 95428584 : kern_1_scf(:) = kernel_scf(:)
305 95428584 : kernel_scf(:) = 0._dp
306 22345873 : loop_iter_i: DO i = 0, n_range
307 22290658 : kern_tot = 0._dp
308 1871962508 : DO j = -m, m
309 1849671850 : ind = 2*i - j
310 1849671850 : IF (ABS(ind) > n_range) THEN
311 : kern = 0._dp
312 : ELSE
313 1636785514 : kern = kern_1_scf(ind)
314 : END IF
315 1871962508 : kern_tot = kern_tot + ch(j)*kern
316 : END DO
317 22290658 : IF (kern_tot == 0._dp) THEN
318 : !zero after (be sure because strictly == 0._dp)
319 : EXIT loop_iter_i
320 : ELSE
321 21731482 : kernel_scf(i) = 0.5_dp*kern_tot
322 21731482 : kernel_scf(-i) = kernel_scf(i)
323 : END IF
324 : END DO loop_iter_i
325 : END DO loop_iter_scf
326 55215 : END SUBROUTINE scf_recurs
327 :
328 : END MODULE ps_wavelet_scaling_function
|
