1#ifndef PRDC_CUDA_FFT_TPP
2#define PRDC_CUDA_FFT_TPP
83 cufftDestroy(rcfPlan_);
86 cufftDestroy(criPlan_);
89 cufftDestroy(ccPlan_);
97 void FFT<D>::setup(IntVec<D>
const & meshDimensions)
105 for (
int i = 0; i < D; ++i) {
107 meshDimensions_[i] = meshDimensions[i];
108 rSize_ *= meshDimensions[i];
110 kSize_ *= meshDimensions[i];
112 kSize_ *= (meshDimensions[i]/2 + 1);
117 if (kFieldCopy_.isAllocated()) {
118 if (kFieldCopy_.capacity() != kSize_) {
119 kFieldCopy_.deallocate();
120 kFieldCopy_.allocate(meshDimensions);
134 void FFT<D>::forwardTransform(RField<D>
const & rField,
135 RFieldDft<D>& kField)
const
145 #ifdef SINGLE_PRECISION
146 result = cufftExecR2C(rcfPlan_,
const_cast<cudaReal*
>(rField.cArray()),
149 result = cufftExecD2Z(rcfPlan_,
const_cast<cudaReal*
>(rField.cArray()),
152 if (result != CUFFT_SUCCESS) {
153 UTIL_THROW(
"Failure in cufft real-to-complex forward transform");
157 cudaReal scale = 1.0/cudaReal(rSize_);
158 VecOp::mulEqS(kField, scale);
165 void FFT<D>::inverseTransformUnsafe(RFieldDft<D>& kField,
166 RField<D>& rField)
const
172 UTIL_CHECK(rField.meshDimensions() == meshDimensions_);
173 UTIL_CHECK(kField.meshDimensions() == meshDimensions_);
176 #ifdef SINGLE_PRECISION
177 result = cufftExecC2R(criPlan_, kField.cArray(), rField.cArray());
179 result = cufftExecZ2D(criPlan_, kField.cArray(), rField.cArray());
181 if (result != CUFFT_SUCCESS) {
182 UTIL_THROW(
"Failure in cufft complex-to-real inverse transform");
191 void FFT<D>::inverseTransformSafe(RFieldDft<D>
const & kField,
192 RField<D>& rField)
const
195 if (kFieldCopy_.isAllocated()) {
197 UTIL_CHECK(kFieldCopy_.meshDimensions() == meshDimensions_);
201 kFieldCopy_ = kField;
204 inverseTransformUnsafe(kFieldCopy_, rField);
213 void FFT<D>::forwardTransform(CField<D>
const & rField, CField<D>& kField)
220 UTIL_CHECK(rField.meshDimensions() == meshDimensions_);
221 UTIL_CHECK(kField.meshDimensions() == meshDimensions_);
226 #ifdef SINGLE_PRECISION
227 result = cufftExecC2C(ccPlan_,
const_cast<cudaComplex*
>(rField.cArray()),
228 kField.cArray(), CUFFT_FORWARD);
230 result = cufftExecZ2Z(ccPlan_,
const_cast<cudaComplex*
>(rField.cArray()),
231 kField.cArray(), CUFFT_FORWARD);
233 if (result != CUFFT_SUCCESS) {
234 UTIL_THROW(
"Failure in cufft complex-to-complex forward transform");
238 cudaReal scale = 1.0/cudaReal(rSize_);
239 VecOp::mulEqS(kField, scale);
246 void FFT<D>::inverseTransform(CField<D>
const & kField, CField<D>& rField)
253 UTIL_CHECK(rField.meshDimensions() == meshDimensions_);
254 UTIL_CHECK(kField.meshDimensions() == meshDimensions_);
259 #ifdef SINGLE_PRECISION
260 result = cufftExecC2C(ccPlan_,
const_cast<cudaComplex*
>(kField.cArray()),
261 rField.cArray(), CUFFT_INVERSE);
263 result = cufftExecZ2Z(ccPlan_,
const_cast<cudaComplex*
>(kField.cArray()),
264 rField.cArray(), CUFFT_INVERSE);
266 if (result != CUFFT_SUCCESS) {
267 UTIL_THROW(
"Failure in cufft complex-to-complex inverse transform");
FFT()
Default constructor.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
PSCF package top-level namespace.
Utility classes for scientific computation.