Commit 9be0a974 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Fourier: added support for libfftw.

parent 3552faf3
......@@ -10,9 +10,68 @@
#include "cdo_output.h"
#include "cdo_options.h"
#include "cdo_fctrans.h"
#include "cimdOmp.h"
void
fc2gp(double *fc, double *gp, long nlat, long nlon, long nfc)
{
#ifdef HAVE_LIBFFTW3
struct FourierMemory
{
fftw_complex *in_fft;
double *out_fft;
fftw_plan plan;
};
std::vector<FourierMemory> ompmem(Threading::ompNumThreads);
for (int i = 0; i < Threading::ompNumThreads; i++)
{
ompmem[i].in_fft = fftw_alloc_complex(nlon / 2 + 1);
ompmem[i].out_fft = (double*) fftw_malloc(nlon*sizeof(double));
ompmem[i].plan = fftw_plan_dft_c2r_1d(nlon, ompmem[i].in_fft, ompmem[i].out_fft, FFTW_ESTIMATE);
}
if (Options::cdoVerbose) fftw_print_plan(ompmem[0].plan);
#ifdef _OPENMP
#pragma omp parallel for default(shared)
#endif
for (long ilat = 0; ilat < nlat; ++ilat)
{
const auto ompthID = cdo_omp_get_thread_num();
auto in_fft = ompmem[ompthID].in_fft;
auto out_fft = ompmem[ompthID].out_fft;
for (long ifc = 0; ifc < nfc/2; ++ifc)
{
in_fft[ifc][0] = fc[2*ifc*nlat+ilat];
in_fft[ifc][1] = fc[(2*ifc+1)*nlat+ilat];
}
for (long ifc = nfc/2; ifc < (nlon / 2 + 1); ++ifc)
{
in_fft[ifc][0] = 0;
in_fft[ifc][1] = 0;
}
fftw_execute(ompmem[ompthID].plan);
for (long ilon = 0; ilon < nlon; ++ilon) gp[ilat*nlon+ilon] = out_fft[ilon];
}
for (int i = 0; i < Threading::ompNumThreads; i++)
{
fftw_free(ompmem[i].in_fft);
fftw_free(ompmem[i].out_fft);
fftw_destroy_plan(ompmem[i].plan);
}
#else
cdoAbort("FFTW support not compiled in!");
#endif
}
void
gp2fc(const double *gp, double *fc, long nlat, long nlon, long nfc)
{
......@@ -44,6 +103,7 @@ gp2fc(const double *gp, double *fc, long nlat, long nlon, long nfc)
const auto ompthID = cdo_omp_get_thread_num();
auto in_fft = ompmem[ompthID].in_fft;
auto out_fft = ompmem[ompthID].out_fft;
for (long ilon = 0; ilon < nlon; ++ilon) in_fft[ilon] = gp[ilat*nlon+ilon];
fftw_execute(ompmem[ompthID].plan);
......
#ifndef CDO_FCTRANS_H
#define CDO_FCTRANS_H
void fc2gp(double *fc, double *gp, long nlat, long nlon, long nfc);
void gp2fc(const double *gp, double *fc, long nlat, long nlon, long nfc);
#endif
......@@ -54,7 +54,11 @@ spec2grid(SPTRANS &sptrans, int gridIDin, double *arrayIn, int gridIDout, double
std::vector<double> fpwork(nlat * nfc * nlev);
sp2fc(arrayIn, fpwork.data(), sptrans.vpoli.data(), nlev, nlat, nfc, ntr);
fc2gp(sptrans.vtrig.data(), sptrans.ifax, fpwork.data(), arrayOut, nlat, nlon, nlev, nfc);
if (sptrans.use_fftw)
fc2gp(fpwork.data(), arrayOut, nlat, nlon, nfc);
else
fc2gp(sptrans.vtrig.data(), sptrans.ifax, fpwork.data(), arrayOut, nlat, nlon, nlev, nfc);
}
void
......@@ -93,7 +97,10 @@ four2grid(SPTRANS &sptrans, int gridIDin, double *arrayIn, int gridIDout, double
const long waves = ntr + 1;
const long nfc = waves * 2;
fc2gp(sptrans.vtrig.data(), sptrans.ifax, arrayIn, arrayOut, nlat, nlon, nlev, nfc);
if (sptrans.use_fftw)
fc2gp(arrayIn, arrayOut, nlat, nlon, nfc);
else
fc2gp(sptrans.vtrig.data(), sptrans.ifax, arrayIn, arrayOut, nlat, nlon, nlev, nfc);
}
void
......@@ -189,9 +196,17 @@ trans_dv2uv(SPTRANS &sptrans, DVTRANS &dvtrans, long nlev, int gridID1, double *
sp2fc(su, fpwork.data(), sptrans.vpoli.data(), nlev, nlat, nfc, ntr);
scaluv(fpwork.data(), sptrans.vrcoslat.data(), nlat, nfc * nlev);
fc2gp(sptrans.vtrig.data(), sptrans.ifax, fpwork.data(), gu, nlat, nlon, nlev, nfc);
if (sptrans.use_fftw)
fc2gp(fpwork.data(), gu, nlat, nlon, nfc);
else
fc2gp(sptrans.vtrig.data(), sptrans.ifax, fpwork.data(), gu, nlat, nlon, nlev, nfc);
sp2fc(sv, fpwork.data(), sptrans.vpoli.data(), nlev, nlat, nfc, ntr);
scaluv(fpwork.data(), sptrans.vrcoslat.data(), nlat, nfc * nlev);
fc2gp(sptrans.vtrig.data(), sptrans.ifax, fpwork.data(), gv, nlat, nlon, nlev, nfc);
if (sptrans.use_fftw)
fc2gp(fpwork.data(), gv, nlat, nlon, nfc);
else
fc2gp(sptrans.vtrig.data(), sptrans.ifax, fpwork.data(), gv, nlat, nlon, nlev, nfc);
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment