Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This file contains procedure implementations of [pm_fftpack](@ref pm_fftpack).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \FatemehBagheri, Wednesday 12:20 PM, September 22, 2021, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : #if setFFTF_ENABLED
28 : real(TKC), parameter :: trifac = +2 * acos(-1._TKC)
29 : #elif setFFTR_ENABLED
30 : real(TKC), parameter :: trifac = -2 * acos(-1._TKC)
31 : #endif
32 : !%%%%%%%%%%%%%%
33 : #if setFFTI_ENABLED
34 : !%%%%%%%%%%%%%%
35 :
36 66 : call setFFTR(data)
37 : #if CK_ENABLED
38 2025 : data = data / size(data, 1, IK)
39 : #elif RK_ENABLED
40 2417 : data = data * (2._TKC / size(data, 1, IK))
41 : #else
42 : #error "Unrecognized interface."
43 : #endif
44 :
45 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
46 : #elif getFFTF_ENABLED || getFFTR_ENABLED || getFFTI_ENABLED
47 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48 :
49 11886 : fft(1 : size(data, 1, IK)) = data
50 2410 : fft(size(data, 1, IK) + 1 :) = 0._TKC
51 : #if getFFTF_ENABLED
52 90 : call setFFTF(fft)
53 : #elif getFFTI_ENABLED
54 60 : call setFFTI(fft)
55 : #elif getFFTR_ENABLED
56 30 : call setFFTR(fft)
57 : #else
58 : #error "Unrecognized interface."
59 : #endif
60 :
61 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
62 : #elif (setFFTF_ENABLED || setFFTR_ENABLED) && CK_ENABLED
63 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
64 :
65 : integer(IK) :: lendata, lenDataHalf, i, j, m, step, maxm
66 : complex(TKC) :: temp, w, wp
67 : real(TKC) :: theta, wtmp
68 108 : lendata = size(data, 1, IK)
69 108 : CHECK_ASSERTION(__LINE__, isIntPow(lenData) .and. 1 < lenData, SK_"@setFFTF/setFFTR(): The condition `isIntPow(size(data)) .and. 1 < size(data)` must hold. size(data) = "//getStr(lendata))
70 108 : lenDataHalf = lenData / 2_IK
71 : j = 1_IK
72 6780 : do i = 1, lendata
73 6672 : if (i < j) then
74 : ! swap elements.
75 2852 : temp = data(j)
76 2852 : data(j) = data(i)
77 2852 : data(i) = temp
78 : endif
79 : m = lenDataHalf
80 6456 : do
81 13128 : if (m < 2_IK .or. j <= m) exit
82 6456 : j = j - m
83 6456 : m = m / 2_IK
84 : end do
85 6780 : j = j + m
86 : end do
87 : maxm = 1_IK
88 : ! Danielson-Lanczos iteration. Repeat `log2(lenDataHalf)` times.
89 : do
90 696 : if (lendata <= maxm) exit
91 588 : step = 2_IK * maxm
92 588 : theta = trifac / step
93 588 : wp%re = -2._TKC * sin(0.5_TKC * theta)**2
94 184 : wp%im = sin(theta)
95 404 : w = (1._TKC, 0._TKC)
96 7152 : do m = 1, maxm
97 6564 : do i = m, lendata, step
98 21080 : j = i + maxm
99 21080 : temp%re = w%re * data(j)%re - w%im * data(j)%im
100 21080 : temp%im = w%re * data(j)%im + w%im * data(j)%re
101 21080 : data(j) = data(i) - temp
102 21080 : data(i) = data(i) + temp
103 : end do
104 6564 : wtmp = w%re
105 6564 : w%re = w%re * wp%re - w%im * wp%im + w%re
106 7152 : w%im = w%im * wp%re + wtmp * wp%im + w%im
107 : end do
108 108 : maxm = step
109 : end do
110 :
111 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
112 : #elif (setFFTF_ENABLED || setFFTR_ENABLED) && RK_ENABLED
113 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
114 :
115 : complex(TKC) :: h1, h2, w, wp
116 : real(TKC) :: c1, c2, theta, wtmp
117 : integer(IK) :: i, i1, i2, i3, i4, lenData, lenDataPlus3
118 108 : lenData = size(data, 1, IK)
119 108 : CHECK_ASSERTION(__LINE__, isIntPow(lenData) .and. 1 < lenData, SK_"@setFFTF/setFFTR(): The condition `isIntPow(size(data)) .and. 1 < size(data)` must hold. size(data) = "//getStr(lendata))
120 108 : theta = trifac / real(lenData, TKC)
121 36 : c1 = 0.5_TKC
122 36 : c2 = 0.5_TKC
123 : #if setFFTF_ENABLED
124 54 : call setFFTC(data)
125 18 : c2 = -c2
126 : #endif
127 108 : wp%re = -2._TKC * sin(0.5_TKC * theta)**2
128 108 : wp%im = sin(theta)
129 72 : w%im = wp%im
130 108 : w%re = 1._TKC + wp%re
131 108 : lenDataPlus3 = lenData + 3_IK
132 1888 : do i = 2_IK, lenData / 4_IK
133 1780 : i1 = 2_IK * i - 1_IK
134 : i2 = i1 + 1_IK
135 1780 : i3 = lenDataPlus3 - i2
136 1780 : i4 = i3 + 1_IK
137 1780 : h1%re = +c1 * (data(i1) + data(i3))
138 1780 : h1%im = +c1 * (data(i2) - data(i4))
139 1780 : h2%re = -c2 * (data(i2) + data(i4))
140 1780 : h2%im = +c2 * (data(i1) - data(i3))
141 1780 : data(i1) = +h1%re + w%re * h2%re - w%im * h2%im
142 1780 : data(i2) = +h1%im + w%re * h2%im + w%im * h2%re
143 1780 : data(i3) = +h1%re - w%re * h2%re + w%im * h2%im
144 1780 : data(i4) = -h1%im + w%re * h2%im + w%im * h2%re
145 604 : wtmp = w%re
146 1780 : w%re = w%re * wp%re - w%im * wp%im + w%re
147 1888 : w%im = w%im * wp%re + wtmp * wp%im + w%im
148 : end do
149 : #if setFFTF_ENABLED
150 54 : h1%re = data(1)
151 54 : data(1) = h1%re + data(2)
152 54 : data(2) = h1%re - data(2)
153 : #elif setFFTR_ENABLED
154 54 : h1%re = data(1)
155 54 : data(1) = c1 * (h1%re + data(2))
156 54 : data(2) = c1 * (h1%re - data(2))
157 54 : call setFFTC(data)
158 : #else
159 : #error "Unrecognized interface."
160 : #endif
161 : contains
162 108 : pure subroutine setFFTC(data)
163 : real(TKC), intent(inout), contiguous :: data(:)
164 : integer(IK) :: i, j, step, m, maxm, lenData, lenDataHalf
165 : real(TKC) :: tempi, tempr, theta, wi, wpi, wpr, wr, wtemp
166 108 : lenData = size(data, 1, IK)
167 108 : lenDataHalf = lenData / 2_IK
168 : j = 1_IK
169 108 : do i = 1_IK, lenData, 2_IK
170 3776 : if(i < j)then
171 1568 : tempr = data(j)
172 1568 : tempi = data(j + 1)
173 1568 : data(j) = data(i)
174 1568 : data(j + 1) = data(i + 1)
175 1568 : data(i) = tempr
176 1568 : data(i + 1) = tempi
177 : endif
178 : m = lenDataHalf
179 3668 : do
180 7444 : if (m < 2_IK .or. j <= m) exit
181 3668 : j = j - m
182 3668 : m = m / 2_IK
183 : end do
184 3776 : j = j + m
185 : end do
186 : maxm = 2_IK
187 : do
188 594 : if (lenData <= maxm) exit
189 486 : step = 2_IK * maxm
190 486 : theta = trifac / maxm
191 486 : wpr = -2._TKC * sin(0.5_TKC * theta)**2
192 486 : wpi = sin(theta)
193 164 : wr = 1._TKC
194 164 : wi = 0._TKC
195 486 : do m = 1_IK, maxm, 2_IK
196 3668 : do i = m, lenData, step
197 10480 : j = i + maxm
198 10480 : tempr = wr * data(j) - wi * data(j + 1)
199 10480 : tempi = wr * data(j + 1) + wi * data(j)
200 10480 : data(j) = data(i) - tempr
201 10480 : data(j + 1) = data(i + 1) - tempi
202 10480 : data(i) = data(i) + tempr
203 10480 : data(i + 1) = data(i + 1) + tempi
204 : end do
205 1244 : wtemp = wr
206 3668 : wr = wr * wpr - wi * wpi + wr
207 3668 : wi = wi * wpr + wtemp * wpi + wi
208 : end do
209 108 : maxm = step
210 : end do
211 108 : end subroutine
212 : #else
213 : !%%%%%%%%%%%%%%%%%%%%%%%%
214 : #error "Unrecognized interface."
215 : !%%%%%%%%%%%%%%%%%%%%%%%%
216 : #endif
|
