Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F82200135
fftw_engine_c2c.cc
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
Tue, Sep 10, 03:31
Size
5 KB
Mime Type
text/x-c
Expires
Thu, Sep 12, 03:31 (2 d)
Engine
blob
Format
Raw Data
Handle
20657683
Attached To
rMSPPROTO µSpectre prototype implementation
fftw_engine_c2c.cc
View Options
/**
* file fftw_engine_c2c.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 11 May 2017
*
* @brief fft_engine usinge FFTW
*
* @section LICENCE
*
* Copyright (C) 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre 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
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include <complex>
#include <array>
#include "system/fftw_engine_c2c.hh"
#include <unsupported/Eigen/CXX11/Tensor>
#include <memory>
namespace
muSpectre
{
template
<
Dim_t
DimS
,
Dim_t
DimM
>
FFTW_EngineC2C
<
DimS
,
DimM
>::
FFTW_EngineC2C
(
Index
nb_pixels
)
:
FFT_Engine
<
DimS
,
DimM
>
(
nb_pixels
)
{
this
->
i_work_space
=
fftw_alloc_complex
(
this
->
tot_nb_pixels
*
DimM
*
DimM
);
this
->
o_work_space
=
fftw_alloc_complex
(
this
->
tot_nb_pixels
*
DimM
*
DimM
);
}
//----------------------------------------------------------------------------//
template
<
Dim_t
DimS
,
Dim_t
DimM
>
FFTW_EngineC2C
<
DimS
,
DimM
>::~
FFTW_EngineC2C
()
{
fftw_destroy_plan
(
this
->
plan_fft
);
fftw_destroy_plan
(
this
->
plan_ifft
);
fftw_free
(
this
->
i_work_space
);
fftw_free
(
this
->
o_work_space
);
fftw_cleanup
();
}
//----------------------------------------------------------------------------//
template
<
Dim_t
DimS
,
Dim_t
DimM
>
void
FFTW_EngineC2C
<
DimS
,
DimM
>::
init
(
FFT_PlanFlags
inflags
)
{
const
int
&
rank
=
DimS
;
std
::
array
<
int
,
DimS
>
narr
;
const
int
*
const
n
=
&
narr
[
0
];
std
::
copy
(
this
->
nb_pixels
.
begin
(),
this
->
nb_pixels
.
end
(),
narr
.
begin
());
int
howmany
=
DimM
*
DimM
;
fftw_complex
*
in
=
reinterpret_cast
<
fftw_complex
*>
(
this
->
i_work_space
);
const
int
*
const
inembed
=
n
;
int
istride
=
howmany
;
int
idist
=
1
;
fftw_complex
*
out
=
in
;
const
int
*
const
onembed
=
n
;
int
ostride
=
howmany
;
int
odist
=
idist
;
int
sign
=
FFTW_FORWARD
;
unsigned
int
flags
;
if
(
inflags
==
FFT_PlanFlags
::
estimate
)
{
flags
=
FFTW_ESTIMATE
;
}
if
(
inflags
==
FFT_PlanFlags
::
measure
)
{
flags
=
FFTW_MEASURE
;
}
if
(
inflags
==
FFT_PlanFlags
::
patient
)
{
flags
=
FFTW_PATIENT
;
}
this
->
plan_fft
=
fftw_plan_many_dft
(
rank
,
n
,
howmany
,
in
,
inembed
,
istride
,
idist
,
out
,
onembed
,
ostride
,
odist
,
sign
,
flags
);
sign
=
FFTW_BACKWARD
;
in
=
reinterpret_cast
<
fftw_complex
*>
(
this
->
o_work_space
);
out
=
reinterpret_cast
<
fftw_complex
*>
(
this
->
o_work_space
);
this
->
plan_ifft
=
fftw_plan_many_dft
(
rank
,
n
,
howmany
,
in
,
inembed
,
istride
,
idist
,
out
,
onembed
,
ostride
,
odist
,
sign
,
flags
);
}
//----------------------------------------------------------------------------//
template
<
Dim_t
DimS
,
Dim_t
DimM
>
Real
*
FFTW_EngineC2C
<
DimS
,
DimM
>::
convolve
(
const
Real
*
const
Ghat
,
const
Real
*
const
arr
)
{
std
::
copy
(
arr
,
arr
+
this
->
tot_nb_pixels
*
DimM
*
DimM
,
reinterpret_cast
<
Complex
*>
(
this
->
i_work_space
));
fftw_execute
(
this
->
plan_fft
);
//----------------------------------------------------------------------------//
const
std
::
array
<
Eigen
::
IndexPair
<
int
>
,
2
>
prod_dims
{
Eigen
::
IndexPair
<
int
>
{
2
,
1
},
Eigen
::
IndexPair
<
int
>
{
3
,
0
}};
//set output to zero
Eigen
::
Map
<
Eigen
::
Matrix
<
Complex
,
Eigen
::
Dynamic
,
1
>>
o_space_handle
(
reinterpret_cast
<
Complex
*>
(
this
->
o_work_space
),
this
->
tot_nb_pixels
*
DimM
*
DimM
)
;
o_space_handle
.
setZero
();
for
(
size_t
pix_id
=
0
;
pix_id
<
this
->
tot_nb_pixels
;
++
pix_id
)
{
using
T4shape
=
Eigen
::
Sizes
<
DimM
,
DimM
,
DimM
,
DimM
>
;
using
T2shape
=
Eigen
::
Sizes
<
DimM
,
DimM
>
;
using
T4type
=
Eigen
::
TensorFixedSize
<
Real
,
T4shape
>
;
using
T2type
=
Eigen
::
TensorFixedSize
<
Complex
,
T2shape
>
;
using
T4map
=
Eigen
::
TensorMap
<
T4type
>
;
using
T2map
=
Eigen
::
TensorMap
<
T2type
>
;
const
Real
*
const
G_ptr
=
&
Ghat
[
pix_id
*
DimM
*
DimM
*
DimM
*
DimM
];
const
T4map
G
(
const_cast
<
Real
*>
(
G_ptr
),
DimM
,
DimM
,
DimM
,
DimM
);
Complex
*
P_ptr
=
reinterpret_cast
<
Complex
*>
(
&
this
->
i_work_space
[
pix_id
*
DimM
*
DimM
]);
const
T2map
P
(
P_ptr
,
DimM
,
DimM
);
Complex
*
Pout_ptr
=
reinterpret_cast
<
Complex
*>
(
&
this
->
o_work_space
[
pix_id
*
DimM
*
DimM
]);
T2map
Pout
(
Pout_ptr
,
DimM
,
DimM
);
////TODO manual contraction to work around complex-
//for (Dim_t i = 0; i < DimM; ++i) {
// for (Dim_t j = 0; j < DimM; ++j) {
// for (Dim_t k = 0; k < DimM; ++k) {
// for (Dim_t l = 0; l < DimM; ++l) {
// Pout(i,j) = G(i,j,k,l)*P(l,k);
// }
// }
// }
//}
Eigen
::
TensorFixedSize
<
Complex
,
T4shape
>
Gc
=
G
.
template
cast
<
Complex
>
();
Pout
=
Gc
.
contract
(
P
,
prod_dims
);
}
//----------------------------------------------------------------------------//
fftw_execute
(
this
->
plan_ifft
);
return
reinterpret_cast
<
Real
*>
(
this
->
o_work_space
);
}
template
class
FFTW_EngineC2C
<
1
,
1
>
;
template
class
FFTW_EngineC2C
<
2
,
2
>
;
template
class
FFTW_EngineC2C
<
3
,
3
>
;
}
// muSpectre
Event Timeline
Log In to Comment