-
Notifications
You must be signed in to change notification settings - Fork 4
/
fft_support.f90
177 lines (136 loc) · 4.61 KB
/
fft_support.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
!
! Copyright (C) Quantum ESPRESSO group
!
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
#include "fft_defs.h"
!=----------------------------------------------------------------------=!
MODULE fft_support
!=----------------------------------------------------------------------=!
USE, intrinsic :: iso_c_binding
IMPLICIT NONE
SAVE
INTEGER, PARAMETER :: DP = selected_real_kind(14,200)
PRIVATE
PUBLIC :: good_fft_dimension, allowed, good_fft_order
!=----------------------------------------------------------------------=!
CONTAINS
!=----------------------------------------------------------------------=!
!
! FFT support Functions/Subroutines
!
!=----------------------------------------------------------------------=!
!
!
integer function good_fft_dimension (n)
!
! Determines the optimal maximum dimensions of fft arrays
! Useful on some machines to avoid memory conflicts
!
IMPLICIT NONE
INTEGER :: n, nx
REAL(DP) :: log2n
!
! this is the default: max dimension = fft dimension
nx = n
!
#if defined(__LINUX_ESSL)
log2n = LOG ( dble (n) ) / LOG ( 2.0_DP )
! log2n is the logarithm of n in base 2
IF ( ABS (NINT(log2n) - log2n) < 1.0d-8 ) nx = n + 1
! if n is a power of 2 (log2n is integer) increase dimension by 1
#elif defined(__SX6)
!
if (mod (n, 2) ==0) nx = n + 1
! for nec vector machines: if n is even increase dimension by 1
!
#endif
!
good_fft_dimension = nx
return
end function good_fft_dimension
!=----------------------------------------------------------------------=!
function allowed (nr)
! find if the fft dimension is a good one
! a "bad one" is either not implemented (as on IBM with ESSL)
! or implemented but with awful performances (most other cases)
implicit none
integer :: nr
logical :: allowed
integer :: pwr (5)
integer :: mr, i, fac, p, maxpwr
integer :: factors( 5 ) = (/ 2, 3, 5, 7, 11 /)
! find the factors of the fft dimension
mr = nr
pwr = 0
factors_loop: do i = 1, 5
fac = factors (i)
maxpwr = NINT ( LOG( DBLE (mr) ) / LOG( DBLE (fac) ) ) + 1
do p = 1, maxpwr
if ( mr == 1 ) EXIT factors_loop
if ( MOD (mr, fac) == 0 ) then
mr = mr / fac
pwr (i) = pwr (i) + 1
endif
enddo
end do factors_loop
IF ( nr /= ( mr * 2**pwr (1) * 3**pwr (2) * 5**pwr (3) * 7**pwr (4) * 11**pwr (5) ) ) &
CALL fftx_error__ (' allowed ', ' what ?!? ', 1 )
if ( mr /= 1 ) then
! fft dimension contains factors > 11 : no good in any case
allowed = .false.
else
#if defined(__LINUX_ESSL)
! IBM machines with essl libraries
allowed = ( pwr(1) >= 1 ) .and. ( pwr(2) <= 2 ) .and. ( pwr(3) <= 1 ) .and. &
( pwr(4) <= 1 ) .and. ( pwr(5) <= 1 ) .and. &
( ( (pwr(2) == 0 ) .and. ( pwr(3) + pwr(4) + pwr(5) ) <= 2 ) .or. &
( (pwr(2) /= 0 ) .and. ( pwr(3) + pwr(4) + pwr(5) ) <= 1 ) )
#else
! fftw and all other cases: no factors 7 and 11
allowed = ( ( pwr(4) == 0 ) .and. ( pwr(5) == 0 ) )
#endif
endif
return
end function allowed
!=----------------------------------------------------------------------=!
INTEGER FUNCTION good_fft_order( nr, np )
!
! This function find a "good" fft order value greater or equal to "nr"
!
! nr (input) tentative order n of a fft
!
! np (optional input) if present restrict the search of the order
! in the ensemble of multiples of np
!
! Output: the same if n is a good number
! the closest higher number that is good
! an fft order is not good if not implemented (as on IBM with ESSL)
! or implemented but with awful performances (most other cases)
!
IMPLICIT NONE
#include "fft_param.f90"
INTEGER, INTENT(IN) :: nr
INTEGER, OPTIONAL, INTENT(IN) :: np
INTEGER :: new
new = nr
IF( PRESENT( np ) ) THEN
DO WHILE( ( ( .NOT. allowed( new ) ) .OR. ( MOD( new, np ) /= 0 ) ) .AND. ( new <= nfftx ) )
new = new + 1
END DO
ELSE
DO WHILE( ( .NOT. allowed( new ) ) .AND. ( new <= nfftx ) )
new = new + 1
END DO
END IF
IF( new > nfftx ) &
CALL fftx_error__( ' good_fft_order ', ' fft order too large ', new )
good_fft_order = new
RETURN
END FUNCTION good_fft_order
!=----------------------------------------------------------------------=!
END MODULE fft_support
!=----------------------------------------------------------------------=!