GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-09-13 01:43:13 +08:00
Code Issues Packages Projects Releases Wiki Activity

Use std::shared_ptr for FFTW/IMKL FFT plan implementation; Fixes #2651

Browse Source
This commit is contained in:
Andrzej Ciarkowski 2023-05-05 16:58:23 +00:00 committed by Rasmus Munk Larsen
parent 1f79a6078f
commit 1698c367a0
2 changed files with 85 additions and 85 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

92
unsupported/Eigen/src/FFT/ei_fftw_impl.h
View File

@ -9,6 +9,8 @@
#include "./InternalHeaderCheck.h"
#include <memory>
namespace Eigen {
namespace internal {
@ -54,41 +56,41 @@ namespace internal {
{
typedef float scalar_type;
typedef fftwf_complex complex_type;
fftwf_plan m_plan;
fftw_plan() :m_plan(NULL) {}
~fftw_plan() {if (m_plan) fftwf_destroy_plan(m_plan);}
std::shared_ptr<fftwf_plan_s> m_plan;
fftw_plan() = default;
void set_plan(fftwf_plan p) { m_plan.reset(p, fftwf_destroy_plan); }
inline
void fwd(complex_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftwf_execute_dft( m_plan, src,dst);
if (m_plan==NULL) set_plan(fftwf_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftwf_execute_dft( m_plan.get(), src,dst);
}
inline
void inv(complex_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftwf_execute_dft( m_plan, src,dst);
if (m_plan==NULL) set_plan(fftwf_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftwf_execute_dft( m_plan.get(), src,dst);
}
inline
void fwd(complex_type * dst,scalar_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftwf_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftwf_execute_dft_r2c( m_plan,src,dst);
if (m_plan==NULL) set_plan(fftwf_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftwf_execute_dft_r2c( m_plan.get(),src,dst);
}
inline
void inv(scalar_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL)
m_plan = fftwf_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftwf_execute_dft_c2r( m_plan, src,dst);
set_plan(fftwf_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftwf_execute_dft_c2r( m_plan.get(), src,dst);
}
inline
void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftwf_execute_dft( m_plan, src,dst);
if (m_plan==NULL) set_plan(fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftwf_execute_dft( m_plan.get(), src,dst);
}
inline
void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftwf_execute_dft( m_plan, src,dst);
if (m_plan==NULL) set_plan(fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftwf_execute_dft( m_plan.get(), src,dst);
}
};
@ -97,40 +99,40 @@ namespace internal {
{
typedef double scalar_type;
typedef fftw_complex complex_type;
::fftw_plan m_plan;
fftw_plan() :m_plan(NULL) {}
~fftw_plan() {if (m_plan) fftw_destroy_plan(m_plan);}
std::shared_ptr<fftw_plan_s> m_plan;
fftw_plan() = default;
void set_plan(::fftw_plan p) { m_plan.reset(p, fftw_destroy_plan); }
inline
void fwd(complex_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftw_execute_dft( m_plan, src,dst);
if (m_plan==NULL) set_plan(fftw_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftw_execute_dft( m_plan.get(), src,dst);
}
inline
void inv(complex_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftw_execute_dft( m_plan, src,dst);
if (m_plan==NULL) set_plan(fftw_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftw_execute_dft( m_plan.get(), src,dst);
}
inline
void fwd(complex_type * dst,scalar_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftw_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftw_execute_dft_r2c( m_plan,src,dst);
if (m_plan==NULL) set_plan(fftw_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftw_execute_dft_r2c( m_plan.get(),src,dst);
}
inline
void inv(scalar_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL)
m_plan = fftw_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftw_execute_dft_c2r( m_plan, src,dst);
set_plan(fftw_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftw_execute_dft_c2r( m_plan.get(), src,dst);
}
inline
void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
if (m_plan==NULL) m_plan = fftw_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftw_execute_dft( m_plan, src,dst);
if (m_plan==NULL) set_plan(fftw_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftw_execute_dft( m_plan.get(), src,dst);
}
inline
void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
if (m_plan==NULL) m_plan = fftw_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftw_execute_dft( m_plan, src,dst);
if (m_plan==NULL) set_plan(fftw_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftw_execute_dft( m_plan.get(), src,dst);
}
};
template <>
@ -138,40 +140,40 @@ namespace internal {
{
typedef long double scalar_type;
typedef fftwl_complex complex_type;
fftwl_plan m_plan;
fftw_plan() :m_plan(NULL) {}
~fftw_plan() {if (m_plan) fftwl_destroy_plan(m_plan);}
std::shared_ptr<fftwl_plan_s> m_plan;
fftw_plan() = default;
void set_plan(fftwl_plan p) { m_plan.reset(p, fftwl_destroy_plan); }
inline
void fwd(complex_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftwl_execute_dft( m_plan, src,dst);
if (m_plan==NULL) set_plan(fftwl_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftwl_execute_dft( m_plan.get(), src,dst);
}
inline
void inv(complex_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftwl_execute_dft( m_plan, src,dst);
if (m_plan==NULL) set_plan(fftwl_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftwl_execute_dft( m_plan.get(), src,dst);
}
inline
void fwd(complex_type * dst,scalar_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftwl_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftwl_execute_dft_r2c( m_plan,src,dst);
if (m_plan==NULL) set_plan(fftwl_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftwl_execute_dft_r2c( m_plan.get(),src,dst);
}
inline
void inv(scalar_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL)
m_plan = fftwl_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftwl_execute_dft_c2r( m_plan, src,dst);
set_plan(fftwl_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftwl_execute_dft_c2r( m_plan.get(), src,dst);
}
inline
void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftwl_execute_dft( m_plan, src,dst);
if (m_plan==NULL) set_plan(fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftwl_execute_dft( m_plan.get(), src,dst);
}
inline
void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
fftwl_execute_dft( m_plan, src,dst);
if (m_plan==NULL) set_plan(fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT));
fftwl_execute_dft( m_plan.get(), src,dst);
}
};

78
unsupported/Eigen/src/FFT/ei_imklfft_impl.h
View File

@ -10,6 +10,7 @@
#include "./InternalHeaderCheck.h"
#include <complex>
#include <memory>
namespace Eigen {
namespace internal {
@ -38,7 +39,7 @@ inline MKL_Complex8* complex_cast(const std::complex<float>* p) {
* Array of type MKL_LONG otherwise. Lengths of each dimension for a
* multi-dimensional transform.
*/
inline void configure_descriptor(DFTI_DESCRIPTOR_HANDLE* handl,
inline void configure_descriptor(std::shared_ptr<DFTI_DESCRIPTOR>& handl,
enum DFTI_CONFIG_VALUE precision,
enum DFTI_CONFIG_VALUE forward_domain,
MKL_LONG dimension, MKL_LONG* sizes) {
@ -46,25 +47,28 @@ inline void configure_descriptor(DFTI_DESCRIPTOR_HANDLE* handl,
dimension == 2 &&
"Transformation dimension must be less than 3.");
DFTI_DESCRIPTOR_HANDLE res = nullptr;
if (dimension == 1) {
RUN_OR_ASSERT(DftiCreateDescriptor(handl, precision, forward_domain,
RUN_OR_ASSERT(DftiCreateDescriptor(&res, precision, forward_domain,
dimension, *sizes),
"DftiCreateDescriptor failed.")
handl.reset(res, [](DFTI_DESCRIPTOR_HANDLE handle) { DftiFreeDescriptor(&handle); });
if (forward_domain == DFTI_REAL) {
// Set CCE storage
RUN_OR_ASSERT(DftiSetValue(*handl, DFTI_CONJUGATE_EVEN_STORAGE,
RUN_OR_ASSERT(DftiSetValue(handl.get(), DFTI_CONJUGATE_EVEN_STORAGE,
DFTI_COMPLEX_COMPLEX),
"DftiSetValue failed.")
}
} else {
RUN_OR_ASSERT(
DftiCreateDescriptor(handl, precision, DFTI_COMPLEX, dimension, sizes),
DftiCreateDescriptor(&res, precision, DFTI_COMPLEX, dimension, sizes),
"DftiCreateDescriptor failed.")
handl.reset(res, [](DFTI_DESCRIPTOR_HANDLE handle) { DftiFreeDescriptor(&handle); });
}
RUN_OR_ASSERT(DftiSetValue(*handl, DFTI_PLACEMENT, DFTI_NOT_INPLACE),
RUN_OR_ASSERT(DftiSetValue(handl.get(), DFTI_PLACEMENT, DFTI_NOT_INPLACE),
"DftiSetValue failed.")
RUN_OR_ASSERT(DftiCommitDescriptor(*handl), "DftiCommitDescriptor failed.")
RUN_OR_ASSERT(DftiCommitDescriptor(handl.get()), "DftiCommitDescriptor failed.")
}
template <typename T>
@ -75,62 +79,59 @@ struct plan<float> {
typedef float scalar_type;
typedef MKL_Complex8 complex_type;
DFTI_DESCRIPTOR_HANDLE m_plan;
std::shared_ptr<DFTI_DESCRIPTOR> m_plan;
plan() : m_plan(0) {}
~plan() {
if (m_plan) DftiFreeDescriptor(&m_plan);
};
plan() = default;
enum DFTI_CONFIG_VALUE precision = DFTI_SINGLE;
inline void forward(complex_type* dst, complex_type* src, MKL_LONG nfft) {
if (m_plan == 0) {
configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 1, &nfft);
configure_descriptor(m_plan, precision, DFTI_COMPLEX, 1, &nfft);
}
RUN_OR_ASSERT(DftiComputeForward(m_plan, src, dst),
RUN_OR_ASSERT(DftiComputeForward(m_plan.get(), src, dst),
"DftiComputeForward failed.")
}
inline void inverse(complex_type* dst, complex_type* src, MKL_LONG nfft) {
if (m_plan == 0) {
configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 1, &nfft);
configure_descriptor(m_plan, precision, DFTI_COMPLEX, 1, &nfft);
}
RUN_OR_ASSERT(DftiComputeBackward(m_plan, src, dst),
RUN_OR_ASSERT(DftiComputeBackward(m_plan.get(), src, dst),
"DftiComputeBackward failed.")
}
inline void forward(complex_type* dst, scalar_type* src, MKL_LONG nfft) {
if (m_plan == 0) {
configure_descriptor(&m_plan, precision, DFTI_REAL, 1, &nfft);
configure_descriptor(m_plan, precision, DFTI_REAL, 1, &nfft);
}
RUN_OR_ASSERT(DftiComputeForward(m_plan, src, dst),
RUN_OR_ASSERT(DftiComputeForward(m_plan.get(), src, dst),
"DftiComputeForward failed.")
}
inline void inverse(scalar_type* dst, complex_type* src, MKL_LONG nfft) {
if (m_plan == 0) {
configure_descriptor(&m_plan, precision, DFTI_REAL, 1, &nfft);
configure_descriptor(m_plan, precision, DFTI_REAL, 1, &nfft);
}
RUN_OR_ASSERT(DftiComputeBackward(m_plan, src, dst),
RUN_OR_ASSERT(DftiComputeBackward(m_plan.get(), src, dst),
"DftiComputeBackward failed.")
}
inline void forward2(complex_type* dst, complex_type* src, int n0, int n1) {
if (m_plan == 0) {
MKL_LONG sizes[2] = {n0, n1};
configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 2, sizes);
configure_descriptor(m_plan, precision, DFTI_COMPLEX, 2, sizes);
}
RUN_OR_ASSERT(DftiComputeForward(m_plan, src, dst),
RUN_OR_ASSERT(DftiComputeForward(m_plan.get(), src, dst),
"DftiComputeForward failed.")
}
inline void inverse2(complex_type* dst, complex_type* src, int n0, int n1) {
if (m_plan == 0) {
MKL_LONG sizes[2] = {n0, n1};
configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 2, sizes);
configure_descriptor(m_plan, precision, DFTI_COMPLEX, 2, sizes);
}
RUN_OR_ASSERT(DftiComputeBackward(m_plan, src, dst),
RUN_OR_ASSERT(DftiComputeBackward(m_plan.get(), src, dst),
"DftiComputeBackward failed.")
}
};
@ -140,62 +141,59 @@ struct plan<double> {
typedef double scalar_type;
typedef MKL_Complex16 complex_type;
DFTI_DESCRIPTOR_HANDLE m_plan;
std::shared_ptr<DFTI_DESCRIPTOR> m_plan;
plan() : m_plan(0) {}
~plan() {
if (m_plan) DftiFreeDescriptor(&m_plan);
};
plan() = default;
enum DFTI_CONFIG_VALUE precision = DFTI_DOUBLE;
inline void forward(complex_type* dst, complex_type* src, MKL_LONG nfft) {
if (m_plan == 0) {
configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 1, &nfft);
configure_descriptor(m_plan, precision, DFTI_COMPLEX, 1, &nfft);
}
RUN_OR_ASSERT(DftiComputeForward(m_plan, src, dst),
RUN_OR_ASSERT(DftiComputeForward(m_plan.get(), src, dst),
"DftiComputeForward failed.")
}
inline void inverse(complex_type* dst, complex_type* src, MKL_LONG nfft) {
if (m_plan == 0) {
configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 1, &nfft);
configure_descriptor(m_plan, precision, DFTI_COMPLEX, 1, &nfft);
}
RUN_OR_ASSERT(DftiComputeBackward(m_plan, src, dst),
RUN_OR_ASSERT(DftiComputeBackward(m_plan.get(), src, dst),
"DftiComputeBackward failed.")
}
inline void forward(complex_type* dst, scalar_type* src, MKL_LONG nfft) {
if (m_plan == 0) {
configure_descriptor(&m_plan, precision, DFTI_REAL, 1, &nfft);
configure_descriptor(m_plan, precision, DFTI_REAL, 1, &nfft);
}
RUN_OR_ASSERT(DftiComputeForward(m_plan, src, dst),
RUN_OR_ASSERT(DftiComputeForward(m_plan.get(), src, dst),
"DftiComputeForward failed.")
}
inline void inverse(scalar_type* dst, complex_type* src, MKL_LONG nfft) {
if (m_plan == 0) {
configure_descriptor(&m_plan, precision, DFTI_REAL, 1, &nfft);
configure_descriptor(m_plan, precision, DFTI_REAL, 1, &nfft);
}
RUN_OR_ASSERT(DftiComputeBackward(m_plan, src, dst),
RUN_OR_ASSERT(DftiComputeBackward(m_plan.get(), src, dst),
"DftiComputeBackward failed.")
}
inline void forward2(complex_type* dst, complex_type* src, int n0, int n1) {
if (m_plan == 0) {
MKL_LONG sizes[2] = {n0, n1};
configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 2, sizes);
configure_descriptor(m_plan, precision, DFTI_COMPLEX, 2, sizes);
}
RUN_OR_ASSERT(DftiComputeForward(m_plan, src, dst),
RUN_OR_ASSERT(DftiComputeForward(m_plan.get(), src, dst),
"DftiComputeForward failed.")
}
inline void inverse2(complex_type* dst, complex_type* src, int n0, int n1) {
if (m_plan == 0) {
MKL_LONG sizes[2] = {n0, n1};
configure_descriptor(&m_plan, precision, DFTI_COMPLEX, 2, sizes);
configure_descriptor(m_plan, precision, DFTI_COMPLEX, 2, sizes);
}
RUN_OR_ASSERT(DftiComputeBackward(m_plan, src, dst),
RUN_OR_ASSERT(DftiComputeBackward(m_plan.get(), src, dst),
"DftiComputeBackward failed.")
}
};