Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78647164
fft_engine.hh
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Aug 22, 04:07
Size
5 KB
Mime Type
text/x-c++
Expires
Sat, Aug 24, 04:07 (2 d)
Engine
blob
Format
Raw Data
Handle
20080471
Attached To
rTAMAAS tamaas
fft_engine.hh
View Options
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef FFT_ENGINE_H
#define FFT_ENGINE_H
/* -------------------------------------------------------------------------- */
#include "fftw_interface.hh"
#include "grid.hh"
#include "grid_base.hh"
#include "grid_hermitian.hh"
#include "partitioner.hh"
#include <algorithm>
#include <array>
#include <map>
#include <memory>
#include <string>
/* -------------------------------------------------------------------------- */
namespace
tamaas
{
class
FFTEngine
{
protected
:
using
plan_t
=
std
::
pair
<
fftw
::
plan
<
Real
>
,
fftw
::
plan
<
Real
>>
;
using
key_t
=
std
::
basic_string
<
UInt
>
;
using
complex_t
=
fftw
::
helper
<
Real
>::
complex
;
public
:
virtual
~
FFTEngine
()
noexcept
=
default
;
/// Execute a forward plan on real and spectral 1D
virtual
void
forward
(
const
Grid
<
Real
,
1
>&
real
,
GridHermitian
<
Real
,
1
>&
spectral
)
=
0
;
/// Execute a forward plan on real and spectral 2D
virtual
void
forward
(
const
Grid
<
Real
,
2
>&
real
,
GridHermitian
<
Real
,
2
>&
spectral
)
=
0
;
/// Execute a backward plan on real and spectral 1D
virtual
void
backward
(
Grid
<
Real
,
1
>&
real
,
GridHermitian
<
Real
,
1
>&
spectral
)
=
0
;
/// Execute a backward plan on real and spectral 2D
virtual
void
backward
(
Grid
<
Real
,
2
>&
real
,
GridHermitian
<
Real
,
2
>&
spectral
)
=
0
;
/// Fill a grid with wavevector values in appropriate ordering
template
<
typename
T
,
UInt
dim
,
bool
hermitian
>
static
Grid
<
T
,
dim
>
computeFrequencies
(
const
std
::
array
<
UInt
,
dim
>&
sizes
);
/// Instanciate an appropriate FFTEngine subclass
static
std
::
unique_ptr
<
FFTEngine
>
makeEngine
(
unsigned
int
flags
=
FFTW_ESTIMATE
);
protected
:
/// Make a transform signature from a pair of grids
template
<
UInt
dim
>
static
key_t
make_key
(
const
Grid
<
Real
,
dim
>&
real
,
const
GridHermitian
<
Real
,
dim
>&
spectral
);
/// Cast to FFTW complex type
static
auto
cast
(
Complex
*
data
)
{
return
reinterpret_cast
<
complex_t
*>
(
data
);
}
static
auto
cast
(
const
Complex
*
data
)
{
return
const_cast
<
complex_t
*>
(
reinterpret_cast
<
const
complex_t
*>
(
data
));
}
};
template
<
UInt
dim
>
FFTEngine
::
key_t
FFTEngine
::
make_key
(
const
Grid
<
Real
,
dim
>&
real
,
const
GridHermitian
<
Real
,
dim
>&
spectral
)
{
if
(
real
.
getNbComponents
()
!=
spectral
.
getNbComponents
())
TAMAAS_EXCEPTION
(
"Components do not match"
);
auto
hermitian_dims
=
GridHermitian
<
Real
,
dim
>::
hermitianDimensions
(
real
.
sizes
());
if
(
not
std
::
equal
(
hermitian_dims
.
begin
(),
hermitian_dims
.
end
(),
spectral
.
sizes
().
begin
()))
TAMAAS_EXCEPTION
(
"Spectral grid does not have hermitian size"
);
// Reserve space for dimensions + components + both last strides
key_t
key
(
real
.
getDimension
()
+
3
,
0
);
std
::
copy_n
(
real
.
sizes
().
begin
(),
dim
,
key
.
begin
());
key
[
dim
]
=
real
.
getNbComponents
();
key
[
dim
+
1
]
=
real
.
getStrides
().
back
();
key
[
dim
+
2
]
=
spectral
.
getStrides
().
back
();
return
key
;
}
template
<
typename
T
,
UInt
dim
,
bool
hermitian
>
Grid
<
T
,
dim
>
FFTEngine
::
computeFrequencies
(
const
std
::
array
<
UInt
,
dim
>&
local_sizes
)
{
// If hermitian is true, we suppose the dimensions of freq are
// reduced based on hermitian symetry and that it has dim components
Grid
<
T
,
dim
>
freq
(
local_sizes
,
dim
);
constexpr
UInt
dmax
=
dim
-
static_cast
<
UInt
>
(
hermitian
);
// For MPI
const
auto
&
n
=
Partitioner
<
dim
>::
global_size
(
local_sizes
);
const
auto
offset
=
Partitioner
<
dim
>::
local_offset
(
freq
)
/
dim
;
#pragma omp parallel for
// to get rid of a compilation warning from nvcc
for
(
UInt
i
=
0
;
i
<
freq
.
getNbPoints
();
++
i
)
{
UInt
index
=
i
+
offset
;
VectorProxy
<
T
,
dim
>
wavevector
(
freq
(
i
*
freq
.
getNbComponents
()));
std
::
array
<
UInt
,
dim
>
tuple
{{
0
}};
/// Computing tuple from index
for
(
Int
d
=
dim
-
1
;
d
>=
0
;
d
--
)
{
tuple
[
d
]
=
index
%
n
[
d
];
index
-=
tuple
[
d
];
index
/=
n
[
d
];
}
if
(
hermitian
)
wavevector
(
dim
-
1
)
=
tuple
[
dim
-
1
];
for
(
UInt
d
=
0
;
d
<
dmax
;
d
++
)
{
// Type conversion
T
td
=
tuple
[
d
];
T
nd
=
n
[
d
];
T
q
=
(
tuple
[
d
]
<
n
[
d
]
/
2
)
?
td
:
td
-
nd
;
wavevector
(
d
)
=
q
;
}
}
return
freq
;
}
}
// namespace tamaas
#endif
Event Timeline
Log In to Comment