Line data Source code
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 : ! **************************************************************************************************
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 528 : 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 528 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: y
43 528 : REAL(KIND=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht
44 :
45 : !Number of points: must be 2**nex
46 :
47 2684448 : a = 0.0_dp
48 2684448 : x = 0.0_dp
49 528 : m = itype + 2
50 528 : CALL lazy_arrays(itype, m, ch, cg, cgt, cht)
51 :
52 528 : ni = 2*itype
53 528 : nrange = ni
54 1584 : ALLOCATE (y(0:nd), stat=i_all)
55 528 : IF (i_all /= 0) THEN
56 0 : WRITE (*, *) ' scaling_function: problem of memory allocation'
57 0 : CPABORT("")
58 : END IF
59 :
60 : ! plot scaling function
61 528 : CALL zero(nd + 1, x)
62 528 : CALL zero(nd + 1, y)
63 528 : nt = ni
64 528 : x(nt/2 - 1) = 1._dp
65 : loop1: DO
66 3168 : nt = 2*nt
67 :
68 3168 : CALL back_trans(nd, nt, x, y, m, ch, cg)
69 3168 : CALL dcopy(nt, y, 1, x, 1)
70 3168 : 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 2684448 : DO i = 0, nd
77 2684448 : 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 528 : DEALLOCATE (ch, cg, cgt, cht)
82 528 : DEALLOCATE (y)
83 528 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: y
100 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht
101 :
102 : !must be 2**nex
103 :
104 0 : a = 0.0_dp
105 0 : x = 0.0_dp
106 0 : m = itype + 2
107 0 : ni = 2*itype
108 0 : CALL lazy_arrays(itype, m, ch, cg, cgt, cht)
109 0 : ALLOCATE (y(0:nd), stat=i_all)
110 0 : IF (i_all /= 0) THEN
111 0 : WRITE (*, *) ' wavelet_function: problem of memory allocation'
112 0 : CPABORT("")
113 : END IF
114 :
115 : ! plot wavelet
116 0 : CALL zero(nd + 1, x)
117 0 : CALL zero(nd + 1, y)
118 0 : nt = ni
119 0 : x(nt + nt/2 - 1) = 1._dp
120 : loop3: DO
121 0 : nt = 2*nt
122 : !WRITE(*,*) 'nd,nt',nd,nt
123 0 : CALL back_trans(nd, nt, x, y, m, ch, cg)
124 0 : CALL dcopy(nd, y, 1, x, 1)
125 0 : 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 0 : DO i = 0, nd - 1
132 0 : 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 0 : DEALLOCATE (ch, cg, cgt, cht)
137 0 : DEALLOCATE (y)
138 :
139 0 : 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 47597 : 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 47597 : REAL(KIND=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht
157 :
158 8089116 : kern_1_scf = 0.0_dp
159 47597 : m = itype + 2
160 47597 : CALL lazy_arrays(itype, m, ch, cg, cgt, cht)
161 47597 : CALL scf_recurs(n_iter, n_range, kernel_scf, kern_1_scf, m, ch)
162 47597 : DEALLOCATE (ch, cg, cgt, cht)
163 :
164 47597 : END SUBROUTINE scf_recursion
165 :
166 : ! **************************************************************************************************
167 : !> \brief Set to zero an array x(n)
168 : !> \param n ...
169 : !> \param x ...
170 : ! **************************************************************************************************
171 1056 : SUBROUTINE zero(n, x)
172 : INTEGER, INTENT(in) :: n
173 : REAL(KIND=dp), INTENT(out) :: x(n)
174 :
175 : INTEGER :: i
176 :
177 5368896 : DO i = 1, n
178 5368896 : x(i) = 0._dp
179 : END DO
180 1056 : 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 0 : 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 0 : y = 0.0_dp
206 0 : DO i = 0, nt/2 - 1
207 0 : y(i) = 0._dp
208 0 : y(nt/2 + i) = 0._dp
209 :
210 0 : DO j = -m + 1, m
211 :
212 : ! periodically wrap index if necessary
213 0 : ind = j + 2*i
214 : loop99: DO
215 0 : IF (ind .LT. 0) THEN
216 0 : ind = ind + nt
217 0 : CYCLE loop99
218 : END IF
219 0 : IF (ind .GE. nt) THEN
220 0 : ind = ind - nt
221 0 : CYCLE loop99
222 : END IF
223 : EXIT loop99
224 : END DO loop99
225 :
226 0 : y(i) = y(i) + cht(j)*x(ind)
227 0 : y(nt/2 + i) = y(nt/2 + i) + cgt(j)*x(ind)
228 : END DO
229 :
230 : END DO
231 :
232 0 : 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 3168 : 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 16103520 : y = 0.0_dp
259 :
260 2644632 : DO i = 0, nt/2 - 1
261 2641464 : y(2*i + 0) = 0._dp
262 2641464 : y(2*i + 1) = 0._dp
263 :
264 113310936 : DO j = -m/2, m/2 - 1
265 :
266 : ! periodically wrap index if necessary
267 110666304 : ind = i - j
268 : loop99: DO
269 112049820 : IF (ind .LT. 0) THEN
270 658728 : ind = ind + nt/2
271 658728 : CYCLE loop99
272 : END IF
273 111391092 : IF (ind .GE. nt/2) THEN
274 724788 : ind = ind - nt/2
275 724788 : CYCLE loop99
276 : END IF
277 : EXIT loop99
278 : END DO loop99
279 :
280 110666304 : y(2*i + 0) = y(2*i + 0) + ch(2*j - 0)*x(ind) + cg(2*j - 0)*x(ind + nt/2)
281 113307768 : 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 3168 : 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 0 : 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 0 : DO i = -m, m
310 0 : DO j = -m, m
311 0 : t1 = 0._dp
312 0 : t2 = 0._dp
313 0 : t3 = 0._dp
314 0 : t4 = 0._dp
315 0 : DO l = -3*m, 3*m
316 : IF (l - 2*i .GE. -m .AND. l - 2*i .LE. m .AND. &
317 0 : l - 2*j .GE. -m .AND. l - 2*j .LE. m) THEN
318 0 : t1 = t1 + ch(l - 2*i)*cht(l - 2*j)
319 0 : t2 = t2 + cg(l - 2*i)*cgt(l - 2*j)
320 0 : t3 = t3 + ch(l - 2*i)*cgt(l - 2*j)
321 0 : t4 = t4 + cht(l - 2*i)*cg(l - 2*j)
322 : END IF
323 : END DO
324 0 : eps = 1.e-10_dp
325 0 : IF (i .EQ. j) THEN
326 : IF (ABS(t1 - 1._dp) .GT. eps .OR. ABS(t2 - 1._dp) .GT. eps .OR. &
327 0 : ABS(t3) .GT. eps .OR. ABS(t4) .GT. eps) THEN
328 0 : 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 0 : ABS(t3) .GT. eps .OR. ABS(t4) .GT. eps) THEN
333 0 : WRITE (*, fmt22) 'Orthogonality ERROR', i, j, t1, t2, t3, t4
334 : END IF
335 : END IF
336 : END DO
337 : END DO
338 :
339 0 : WRITE (*, *) 'FILTER TEST PASSED'
340 :
341 0 : 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 47597 : 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 8089116 : kern_1_scf = 0.0_dp
364 : !Start the iteration to go from p0gauss to pgauss
365 541295 : loop_iter_scf: DO i_iter = 1, n_iter
366 81500504 : kern_1_scf(:) = kernel_scf(:)
367 81500504 : kernel_scf(:) = 0._dp
368 19471599 : loop_iter_i: DO i = 0, n_range
369 19424002 : kern_tot = 0._dp
370 1653575820 : DO j = -m, m
371 1634151818 : ind = 2*i - j
372 1634151818 : IF (ABS(ind) > n_range) THEN
373 : kern = 0._dp
374 : ELSE
375 1443764078 : kern = kern_1_scf(ind)
376 : END IF
377 1653575820 : kern_tot = kern_tot + ch(j)*kern
378 : END DO
379 19424002 : IF (kern_tot == 0._dp) THEN
380 : !zero after (be sure because strictly == 0._dp)
381 : EXIT loop_iter_i
382 : ELSE
383 18930304 : kernel_scf(i) = 0.5_dp*kern_tot
384 18930304 : kernel_scf(-i) = kernel_scf(i)
385 : END IF
386 : END DO loop_iter_i
387 : END DO loop_iter_scf
388 47597 : END SUBROUTINE scf_recurs
389 :
390 : END MODULE ps_wavelet_scaling_function
|
