Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_ifft
8 : use m_juDFT
9 : CONTAINS
10 0 : INTEGER FUNCTION ifft235(ksfft, n, gmaxp)
11 : !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
12 : ! this function checks wether n can be expressed as :
13 : ! n = (2**P) * (3**Q) * (5**R) to match withe the MFFT - routines
14 : ! used in this program. If close to n there is a number n' which
15 : ! is solely expressable as n'= 2**p then this number sis choosen
16 : ! and is outputted by the function. If n is not expressable as
17 : ! n = (2**P) * (3**Q) * (5**R) a number n' is choosen close to ns
18 : ! is selected which fullfilles these requirements.
19 :
20 : ! Stefan Bl"ugel , kfa, Oct. 1993
21 : !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
22 : IMPLICIT NONE
23 :
24 : !----> declaration part
25 :
26 : INTEGER :: ksfft, n
27 : REAL :: gmaxp
28 :
29 : !----> local variabels
30 :
31 : INTEGER :: np, nn, nl, iexp, is, itrial
32 : REAL :: rl
33 : REAL :: fac_1, two, fac_2, fac_3, one, fac_4
34 : PARAMETER(two=2.0e0, fac_1=1.001e0, fac_2=1.7e0, fac_3=1.95e0)
35 : PARAMETER(one=1.0e0, fac_4=0.03e0)
36 : INTRINSIC abs, log, mod, real
37 :
38 0 : ifft235 = 0
39 0 : IF (n == 0) THEN
40 0 : CALL juDFT_error("n should not be zero", calledby="ifft235")
41 : ENDIF
42 :
43 : !====> RADIX 2 CALCULATION ONLY
44 :
45 0 : IF (ksfft == 0) THEN
46 0 : rl = log(real(n))/log(two)
47 0 : nl = rl
48 0 : IF (abs(n - 2**nl) > abs(n - 2**(nl + 1))) THEN
49 0 : np = nl + 1
50 : ELSE
51 : np = nl
52 : ENDIF
53 0 : IF ((np == nl) .AND. (gmaxp*real(nl)/rl < fac_1)) &
54 0 : np = nl + 1
55 0 : ifft235 = 2**np
56 0 : RETURN
57 : ENDIF
58 :
59 : !====> RADIX 2 , 3, 5 CALCULATION
60 :
61 : !----> since gmaxp is large enough, try also smaller n
62 :
63 0 : IF (gmaxp > fac_2 .AND. gmaxp < fac_3) THEN
64 : is = -1
65 : ELSE
66 0 : is = 1
67 : ENDIF
68 0 : np = n
69 0 : nn = n
70 :
71 0 : DO itrial = 1, 200
72 :
73 : !----> check whether there is a number close by, which is 2**NL
74 : ! ( FFT very fast )
75 :
76 0 : rl = log(real(nn))/log(two)
77 0 : nl = rl
78 0 : IF ((abs(real(nl)/rl - one) < fac_4) .AND. &
79 0 : (2**nl >= n)) THEN
80 0 : IF (gmaxp*real(nl)/rl > one) ifft235 = 2**nl
81 0 : RETURN
82 0 : ELSE IF ((abs(real(nl + 1)/rl - one) < fac_4) .AND. &
83 : (2**(nl + 1) >= n)) THEN
84 0 : ifft235 = 2**(nl + 1)
85 0 : RETURN
86 : ELSE
87 :
88 : !----> no , no binary number is arround, check wether number can be
89 : ! divided by 2,3 or 5
90 :
91 0 : DO iexp = 1, nl + 1
92 0 : IF (mod(nn, 2) == 0) THEN
93 0 : nn = nn/2
94 0 : ELSE IF (mod(nn, 3) == 0) THEN
95 0 : nn = nn/3
96 0 : ELSE IF (mod(nn, 5) == 0) THEN
97 0 : nn = nn/5
98 : ENDIF
99 0 : IF (nn == 2 .OR. nn == 3 .OR. nn == 5) THEN
100 : !---> o.k.
101 0 : ifft235 = np
102 0 : RETURN
103 : ENDIF
104 : ENDDO
105 :
106 : !----> no , n cannot be expressed as (2**P) * (3**Q) * (5**R)
107 : ! change number
108 0 : nn = n + (is**(itrial))*((itrial + 1)/2)
109 0 : np = n + (is**(itrial))*((itrial + 1)/2)
110 : ENDIF
111 : ENDDO
112 :
113 0 : call juDFT_error(' to few trials '//int2str(itrial))
114 :
115 0 : END FUNCTION ifft235
116 : !-----------------------------------------------------------------------
117 0 : INTEGER FUNCTION i2357(ii)
118 :
119 : ! simple setup to determine fft-length for ESSL calls
120 :
121 : IMPLICIT NONE
122 : INTEGER, INTENT(IN) :: ii
123 : INTEGER :: i, j, k, h, m, n, nn
124 :
125 0 : nn = 12582912
126 0 : DO m = 0, 1
127 0 : DO k = 0, 1
128 0 : DO j = 0, 1
129 0 : DO i = 0, 1
130 0 : DO h = 1, 25
131 0 : n = 2**h*3**i*5**j*7**k*11**m
132 0 : IF ((n >= ii) .AND. (n < nn)) nn = n
133 : ENDDO
134 : ENDDO
135 : ENDDO
136 : ENDDO
137 : ENDDO
138 0 : i2357 = nn
139 :
140 0 : END FUNCTION i2357
141 :
142 :
143 54573 : function next_optimal_fft_size(n) result(fft_size)
144 : implicit none
145 : integer, intent(in) :: n
146 : integer :: fft_size
147 :
148 54573 : fft_size = n
149 154097 : do while(.not. is_optimal_fft_size(fft_size))
150 49762 : fft_size = fft_size + 1
151 : enddo
152 54573 : end function next_optimal_fft_size
153 :
154 104335 : logical function is_optimal_fft_size(n)
155 : implicit none
156 : integer, intent(in) :: n
157 : integer :: i, cnt, divisor
158 : integer, parameter :: primes(4) = [2,3,5,7]
159 :
160 104335 : i = n
161 521675 : do cnt =1,4
162 417340 : divisor = primes(cnt)
163 657694 : do while(mod(i,divisor) == 0)
164 136019 : i = i / divisor
165 : enddo
166 : enddo
167 :
168 104335 : is_optimal_fft_size = (i == 1)
169 0 : end function is_optimal_fft_size
170 : END MODULE m_ifft
|
