@*
FourierTransformFFTW3.tmpl

Created by Graham Dennis on 2007-08-23.

Copyright (c) 2007-2012, Graham Dennis

This program 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 2 of the License, or
(at your option) any later version.

This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.

*@
@extends xpdeint.Features.Transforms._FourierTransformFFTW3

@def description: FFTW3

@attr $planType = "FFTW_MEASURE"
@attr $supportsInPlaceOperation = True

@def includes
  @#
  @super
  @#
#if (CFG_COMPILER == CFG_COMPILER_MSVC)
  #define FFTW_DLL
#endif

#include <fftw3.h>
#include <sys/stat.h>
#include <sys/types.h>

#define _xmds_malloc ${fftwPrefix}_malloc
#define xmds_free ${fftwPrefix}_free
@end def

@def globals
@*doc:
Return the string defining the globals needed by FFTW3.
*@
  @#
const real _inverse_sqrt_2pi = 1.0 / sqrt(2.0 * M_PI); 
string _fftwWisdomPath;
  @#
@end def

@def transformFunction(transformID, transformDict, function)
  @#
  @set runtimePrefix, prefixLattice, postfixLattice, runtimePostfix = transformDict['transformSpecifier']
// _prefix_lattice should be ${prefixLattice}${''.join([' * ' + runtimeLattice for runtimeLattice in runtimePrefix])}
// _postfix_lattice should be ${postfixLattice}${''.join([' * ' + runtimeLattice for runtimeLattice in runtimePostfix])}
static ${fftwPrefix}_plan _fftw_forward_plan = NULL;
static ${fftwPrefix}_plan _fftw_backward_plan = NULL;

if (!_fftw_forward_plan) {
  _LOG(_SIMULATION_LOG_LEVEL, "Planning for ${function.description}...");
  
  @set $transformPair = transformDict['transformPair']
  @set $dimensionsBeingTransformed = len(transformPair[0])
  @set $transformType = transformDict['transformType']
  ${fftwPrefix}_iodim _transform_sizes[${dimensionsBeingTransformed}], _loop_sizes[2];
  @if transformType == 'real'
  ${fftwPrefix}_r2r_kind _r2r_kinds[${dimensionsBeingTransformed}];
  @end if
  ${fftwPrefix}_iodim *_iodim_ptr = NULL;
  
  int _transform_sizes_index = 0, _loop_sizes_index = 0;
  
  if (_prefix_lattice > 1) {
    _iodim_ptr = &_loop_sizes[_loop_sizes_index++];
    _iodim_ptr->n = _prefix_lattice;
    _iodim_ptr->is = _iodim_ptr->os = _postfix_lattice * ${' * '.join([dimRep.globalLattice for dimRep in transformPair[0]])};
  }
  if (_postfix_lattice > 1) {
    _iodim_ptr = &_loop_sizes[_loop_sizes_index++];
    _iodim_ptr->n = _postfix_lattice;
    _iodim_ptr->is = _iodim_ptr->os = 1;
  }
  @#
  @for dimID, dimRep in enumerate(transformPair[0])
  _iodim_ptr = &_transform_sizes[_transform_sizes_index++];
  _iodim_ptr->n = ${dimRep.globalLattice};
  _iodim_ptr->is = _iodim_ptr->os = _postfix_lattice${''.join(c' * ${dr.globalLattice}' for dr in transformPair[0][dimID+1:])};
  
  @end for
  @#
  @if transformType == 'complex'
    @set $guruPlanFunction = self.createGuruDFTPlanInDirection
    @set $executeSuffix = 'dft'
    @set $reinterpretType = $fftwPrefix + '_complex'
  @else
    @set $guruPlanFunction = self.createGuruR2RPlanInDirection
    @set $executeSuffix = 'r2r'
    @set $reinterpretType = 'real'
  @end if
  @#
  @set $dataOut = '_data_out' if transformDict.get('outOfPlace', False) else '_data_in'
  @set $flags = ' | FFTW_DESTROY_INPUT' if transformDict.get('outOfPlace', False) else ''
  @#
  
  ${guruPlanFunction(transformDict, 'forward', dataOut, flags), autoIndent=True}@slurp
  ${guruPlanFunction(transformDict, 'backward', dataOut, flags), autoIndent=True}@slurp
  
  _LOG(_SIMULATION_LOG_LEVEL, " done.\n");
}

if (_forward) {
  ${fftwPrefix}_execute_${executeSuffix}(
    _fftw_forward_plan,
    reinterpret_cast<${reinterpretType}*>(_data_in),
    reinterpret_cast<${reinterpretType}*>(${dataOut})
  );
} else {
  ${fftwPrefix}_execute_${executeSuffix}(
    _fftw_backward_plan,
    reinterpret_cast<${reinterpretType}*>(_data_in),
    reinterpret_cast<${reinterpretType}*>(${dataOut})
  );
}
  @#
@end def

@def createGuruDFTPlanInDirection($transformDict, $direction, $dataOut, $flags)
  @#
_fftw_${direction}_plan = ${fftwPrefix}_plan_guru_dft(
  _transform_sizes_index, _transform_sizes,
  _loop_sizes_index, _loop_sizes,
  reinterpret_cast<${fftwPrefix}_complex*>(_data_in), reinterpret_cast<${fftwPrefix}_complex*>($dataOut),
  FFTW_${direction.upper()}, ${planType}${flags}
);
if (!_fftw_${direction}_plan)
  _LOG(_ERROR_LOG_LEVEL, "(%s: %i) Unable to create ${direction} dft plan.\n", __FILE__, __LINE__);

  @#
@end def

@def createGuruR2RPlanInDirection($transformDict, $direction, $dataOut, $flags)
  @#
  @for idx, dimRep in enumerate(transformDict['transformPair'][0])
_r2r_kinds[${idx}] = ${r2rKindForDimensionAndDirection(dimRep.name, direction)};
  @end for

_fftw_${direction}_plan = ${fftwPrefix}_plan_guru_r2r(
  _transform_sizes_index, _transform_sizes,
  _loop_sizes_index, _loop_sizes,
  reinterpret_cast<real*>(_data_in), reinterpret_cast<real*>($dataOut),
  _r2r_kinds, ${planType}${flags}
);
if (!_fftw_${direction}_plan)
  _LOG(_ERROR_LOG_LEVEL, "(%s: %i) Unable to create ${direction} r2r plan.\n", __FILE__, __LINE__);

  @#
@end def

@def mainBegin($dict)
// load wisdom
#if CFG_OSAPI == CFG_OSAPI_POSIX // Don't load wisdom on windows
${loadWisdom}@slurp
#endif // POSIX
@end def


@def loadWisdom
  @#
{
  char _hostName[256];
  gethostname(_hostName, 256);
  _hostName[255] = '\0'; // just in case
  
  string _pathToWisdom = getenv("HOME");
  _pathToWisdom += "/.xmds/wisdom/";
  
  string _wisdomFileName = _hostName;
  _wisdomFileName += ".wisdom";
  _wisdomFileName += "${wisdomExtension}";
  
  FILE *_fp = NULL;
  
  _fp = fopen(_pathToWisdom.c_str(), "r");
  if (_fp) {
    fclose(_fp);
  } else {
    int _result = mkdir((string(getenv("HOME")) + "/.xmds").c_str(), S_IRWXU);
    if (mkdir(_pathToWisdom.c_str(), S_IRWXU)) {
      // We failed to create the ~/.xmds/wisdom directory
      _LOG(_WARNING_LOG_LEVEL, "Warning: Cannot find enlightenment, the path to wisdom ~/.xmds/wisdom doesn't seem to exist and we couldn't create it.\n"
                               "         I'll use the current path instead.\n");
      _pathToWisdom = ""; // present directory
    }
    
  }
  
  _fftwWisdomPath = _pathToWisdom + _wisdomFileName;
  
  FILE *_wisdomFile = NULL;
  if ( (_wisdomFile = fopen(_fftwWisdomPath.c_str(), "r")) != NULL) {
    _LOG(_SIMULATION_LOG_LEVEL, "Found enlightenment... (Importing wisdom)\n");
    ${fftwPrefix}_import_wisdom_from_file(_wisdomFile);
    fclose(_wisdomFile);
  }
}
  @#
@end def


@def saveWisdom
  @#
{
  FILE *_wisdomFile = NULL;
  if ( (_wisdomFile = fopen(_fftwWisdomPath.c_str(), "w")) != NULL) {
    ${fftwPrefix}_export_wisdom_to_file(_wisdomFile);
    fclose(_wisdomFile);
  }
}
  @#
@end def


@def mainEnd($dict)
  @#

// Save wisdom
#if CFG_OSAPI == CFG_OSAPI_POSIX
${saveWisdom, autoIndent=True}@slurp
#endif // POSIX

${fftwPrefix}_cleanup();
  @#
@end def