Commit 2d760561 authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

FFTW3 version compiles

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/branches/libenv201008
SVN Revision: 3196
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 8a129987
......@@ -25,6 +25,7 @@
# - new doxygen definitions
# - package creation not yet implemented
# - set TF_REMCMMNT to cat if not defined
# 07/10/2010 V1.4 - migrate to FFTW3
#
# ============================================================================
#
......@@ -60,8 +61,17 @@ AR=ar
AS=as
RANLIB=ranlib
#----------------------------------------------------------------------
# uncomment this to activate the FFTW fallback solution
# i.e. to use FFTW instead of FFTW3
# this option will disappear in the near future! Better not use it...
#FFTWFALLBACK=-DFFTWFALLBACK
#----------------------------------------------------------------------
CFLAGS += -O2 -I${SERVERINCLUDEDIR} -I${LOCINCLUDEDIR}
CPPFLAGS += -O2 -I${SERVERINCLUDEDIR} -I${LOCINCLUDEDIR}
CPPFLAGS += -O2 -I${SERVERINCLUDEDIR} -I${LOCINCLUDEDIR} \
$(FFTWFALLBACK)
FFLAGS += -ff2c -Wall -ffixed-line-length-0 -fno-backslash $(FLAGS)
LIBSRC=$(wildcard *.f)
......@@ -98,9 +108,6 @@ libfourier.a: $(LIBOBS)
# http://www.eti.pg.gda.pl/KATEDRY/kecs/lab-cpp/snippets/
# If it is not available to you, you should set TF_REMCMMNT=cat (see below)
#
# REVISIONS and CHANGES
# 30/12/2002 V1.0 Thomas Forbriger
#
# ============================================================================
#
# environment variables
......
......@@ -14,18 +14,19 @@
* REVISIONS and CHANGES
* - 11/07/2006 V1.0 Thomas Forbriger
* - 12/09/2007 V1.1 first running version
*
* - 07/10/2010 V1.2
* - migrate to FFTW3
*
* ============================================================================
*/
#define TF_FFTWAFF_CC_VERSION \
"TF_FFTWAFF_CC V1.1"
"TF_FFTWAFF_CC V1.2"
#define TF_FFTWAFF_CC_CVSID \
"$Id$"
#include <fourier/fftwaff.h>
#include <fourier/error.h>
#include <tfxx/error.h>
#include <tfxx/misc.h>
#include <fourier/error.h>
#include <aff/seriesoperators.h>
#include <aff/dump.h>
#include <cmath>
......@@ -42,9 +43,16 @@ namespace fourier {
if (Mplan_forward==0)
{
//Mplan_forward=rfftw_create_plan(Msize, FFTW_FORWARD, 0);
#ifdef FFTWFALLBACK
Mplan_forward=rfftw_create_plan(Msize, FFTW_REAL_TO_COMPLEX,
FFTW_ESTIMATE);
TFXX_assert(Mplan_forward!=0,
#else
this->create_arrays();
Mplan_forward=fftw_plan_dft_r2c_1d(Msize, Mseriesarray,
Mspectrumarray,
FFTW_ESTIMATE);
#endif
FOURIER_assert(Mplan_forward!=0,
"Error (DRFFTWAFF::create_plan_forward): "
"could not create plan!")
}
......@@ -57,8 +65,15 @@ namespace fourier {
{
if (Mplan_backward==0)
{
#ifdef FFTWFALLBACK
Mplan_backward=rfftw_create_plan(Msize, FFTW_BACKWARD, 0);
TFXX_assert(Mplan_backward!=0,
#else
this->create_arrays();
Mplan_backward=fftw_plan_dft_c2r_1d(Msize, Mspectrumarray,
Mseriesarray,
FFTW_ESTIMATE);
#endif
FOURIER_assert(Mplan_backward!=0,
"Error (DRFFTWAFF::create_plan_backward): "
"could not create plan!")
}
......@@ -70,6 +85,9 @@ namespace fourier {
DRFFTWAFF::~DRFFTWAFF()
{
this->delete_plans();
#ifndef FFTWFALLBACK
fftw_cleanup();
#endif
} // DRFFTWAFF::~DRFFTWAFF()
/*----------------------------------------------------------------------*/
......@@ -81,6 +99,9 @@ namespace fourier {
{
Msize=n;
this->delete_plans();
#ifdef FFTWFALLBACK
this->delete_arrays();
#endif
}
} // DRFFTWAFF::~DRFFTWAFF()
......@@ -89,11 +110,52 @@ namespace fourier {
//! delete plans.
void DRFFTWAFF::delete_plans() const
{
#ifdef FFTWFALLBACK
if (Mplan_forward != 0) { rfftw_destroy_plan(Mplan_forward); }
if (Mplan_backward != 0) { rfftw_destroy_plan(Mplan_backward); }
#else
if (Mplan_forward != 0) { fftw_destroy_plan(Mplan_forward); }
if (Mplan_backward != 0) { fftw_destroy_plan(Mplan_backward); }
#endif
Mplan_forward=0;
Mplan_backward=0;
} // DRFFTWAFF::~DRFFTWAFF()
} // DRFFTWAFF::delete_plans()
/*----------------------------------------------------------------------*/
#ifndef FFTWFALLBACK
//! create plans.
void DRFFTWAFF::create_arrays() const
{
this->delete_arrays();
Mseriesarray = (double *) fftw_malloc(sizeof(double)*Msize);
FOURIER_assert(Mseriesarray!=0,
"Error (DRFFTWAFF::create_plan_forward): "
"could not create series array!")
Mseries=Tseries(Tseries::Trepresentation(Mseriesarray, Msize));
Mspectrumarray = (fftw_complex *)
fftw_malloc(sizeof(double)*this->ssize());
FOURIER_assert(Mspectrumarray!=0,
"Error (DRFFTWAFF::create_plan_forward): "
"could not create spectrum array!")
Mspectrum=Tspectrum(Tspectrum::Trepresentation(
reinterpret_cast<Tcoeff*>(Mspectrumarray),
this->ssize()));
} // DRFFTWAFF::create_arrays()
/*----------------------------------------------------------------------*/
//! delete arrays.
void DRFFTWAFF::delete_arrays() const
{
if (Mseriesarray != 0) { fftw_free(Mseriesarray); }
if (Mspectrumarray != 0) { fftw_free(Mspectrumarray); }
Mplan_forward=0;
Mplan_backward=0;
Mseries=Tseries(0);
Mspectrum=Tspectrum(0);
} // DRFFTWAFF::delete_arrays()
#endif
/*----------------------------------------------------------------------*/
......@@ -105,22 +167,23 @@ namespace fourier {
const bool& debug) const
{
this->set_size(s.size());
#ifdef FFTWFALLBACK
Tspectrum retval(this->Msize/2+1);
aff::Series<fftw_real> out(this->Msize);
aff::Series<fftw_real> in(this->Msize);
fftw_real* pout=out.pointer();
fftw_real* pin=in.pointer();
TFXX_debug(debug, "DRFFTWAFF::operator()",
FOURIER_debug(debug, "DRFFTWAFF::operator()",
"processing arrays are created; copy in series");
in.copyin(s);
TFXX_debug(debug, "DRFFTWAFF::operator()",
FOURIER_debug(debug, "DRFFTWAFF::operator()",
"create plan forward");
this->create_plan_forward();
TFXX_debug(debug, "DRFFTWAFF::operator()",
FOURIER_debug(debug, "DRFFTWAFF::operator()",
"call rfftw_one");
rfftw_one(Mplan_forward, pin, pout);
TFXX_debug(debug, "DRFFTWAFF::operator()",
FOURIER_debug(debug, "DRFFTWAFF::operator()",
"copy results to output");
retval(0)=out(0);
for (int i=1; i<((Msize+1)/2); ++i)
......@@ -134,13 +197,33 @@ namespace fourier {
/*
if (debug)
{
DUMP(s);
DUMP(in);
DUMP(out);
}
*/
#else
Tspectrum retval(this->ssize());
this->create_plan_forward();
Mseries.copyin(s);
/*
for (int i=0; i<this->size(); ++i)
{ Mseries(Mseries.first()+i)=s(s.first()+i); }
*/
fftw_execute(Mplan_forward);
retval.copyin(Mspectrum);
/*
for (int i=0; <this->ssize(); ++i)
{ retval(retval.first()+i)=Mspectrum(Mspectrum.first()+i); }
*/
#endif
/*
if (debug)
{
DUMP(s);
DUMP(retval);
}
*/
TFXX_debug(debug, "DRFFTWAFF::operator()",
FOURIER_debug(debug, "DRFFTWAFF::operator()",
"finished; return");
return(retval);
} // Tspectrum DRFFTWAFF::operator()(const Tseries::Tcoc& s) const
......@@ -161,6 +244,7 @@ namespace fourier {
{ --seriessize; }
this->set_size(seriessize);
#ifdef FFTWFALLBACK
Tseries retval(Msize);
aff::Series<fftw_real> out(Msize);
aff::Series<fftw_real> in(Msize);
......@@ -179,6 +263,13 @@ namespace fourier {
this->create_plan_backward();
rfftw_one(Mplan_backward, pin, pout);
retval.copyin(out);
#else
Tseries retval(this->size());
this->create_plan_backward();
Mspectrum.copyin(s);
fftw_execute(Mplan_backward);
retval.copyin(Mseries);
#endif
return(retval);
} // Tseries DRFFTWAFF::operator()(const Tspectrum::Tcoc& s) const
......
......@@ -16,6 +16,37 @@
* REVISIONS and CHANGES
* - 11/07/2006 V1.0 Thomas Forbriger
* - 12/09/2007 V1.1 first running version
* - 07/10/2010 V1.2
* - migrate to FFTW3:
* - use different fftw header file (fftw3.h)
* - type of FFTW plan has changed
*
* Migration to FFTW3:
* - The location of the arrays in memory are part of the plan.
* Consequently we have to either create a new plan for each transform or to
* allocate the working array together with the plan and keep it available in
* the background. The solution of choice is to keep arrays. We will
* introduce a control parameter to the arguments of the constructor which
* optionally allows to get rid of the arrays after each transformation.
* - FFTW3 uses wisdom by default.
* \code
* void fftw_forget_wisdom(void);
* \endcode
* should be called by the destructor of the class. Alternatively
* \code
* void fftw_cleanup(void);
* \endcode
* could be called which apparently is even more rigorous. Plans, however,
* must be destroyed prior to calling \c fftw_cleanup.
* - The FFTW documentation states:
* To the extent that this is true, if you have a variable
* \code
* complex<double> *x,
* \endcode
* you can pass it directly to FFTW via
* \code
* reinterpret_cast<fftw_complex*>(x).
* \endcode
*
* ============================================================================
*/
......@@ -24,12 +55,16 @@
#ifndef TF_FFTWAFF_H_VERSION
#define TF_FFTWAFF_H_VERSION \
"TF_FFTWAFF_H V1.1"
"TF_FFTWAFF_H V1.2"
#define TF_FFTWAFF_H_CVSID \
"$Id$"
#include<complex>
#ifdef FFTWFALLBACK
#include<drfftw.h>
#else
#include<fftw3.h>
#endif
#include<aff/series.h>
namespace fourier {
......@@ -60,10 +95,21 @@ namespace fourier {
typedef std::complex<Tsample> Tcoeff;
typedef aff::Series<Tsample> Tseries;
typedef aff::Series<Tcoeff> Tspectrum;
#ifdef FFTWFALLBACK
DRFFTWAFF(const int& n):
Msize(n), Mplan_forward(0), Mplan_backward(0) { }
DRFFTWAFF():
Msize(0), Mplan_forward(0), Mplan_backward(0) { }
#else
DRFFTWAFF(const int& n, const bool& deletearrays=false):
Msize(n), Mplan_forward(0), Mplan_backward(0),
Mseriesarray(0), Mspectrumarray(0),
Mdeletearrays(deletearrays) { }
DRFFTWAFF(const bool& deletearrays=false):
Msize(0), Mplan_forward(0), Mplan_backward(0),
Mseriesarray(0), Mspectrumarray(0),
Mdeletearrays(deletearrays) { }
#endif
~DRFFTWAFF();
Tspectrum operator()(const Tseries::Tcoc& s,
const bool& debug=false) const;
......@@ -82,10 +128,25 @@ namespace fourier {
void create_plan_forward() const;
void create_plan_backward() const;
void delete_plans() const;
#ifndef FFTWFALLBACK
void create_arrays() const;
void delete_arrays() const;
int ssize() const { return(this->size()/2+1); }
#endif
void set_size(const int& n) const;
mutable int Msize;
#ifdef FFTWFALLBACK
mutable rfftw_plan Mplan_forward;
mutable rfftw_plan Mplan_backward;
#else
mutable fftw_plan Mplan_forward;
mutable fftw_plan Mplan_backward;
mutable double *Mseriesarray;
mutable fftw_complex *Mspectrumarray;
mutable Tspectrum Mspectrum;
mutable Tseries Mseries;
bool Mdeletearrays;
#endif
}; // class DRFFTWAFF
} // namespace ftt
......
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