/*
 Copyright (C) 2000-2004

 Code contributed by Greg Collecutt, Joseph Hope and Paul Cochrane

 This file is part of xmds.
 
 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, write to the Free Software
 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/

/*
  $Id: xmdsvector.cc,v 1.23 2004/10/21 09:50:07 paultcochrane Exp $
*/

/*! @file xmdsvector.cc
  @brief Vector object classes and methods

  More detailed explanation...
*/

#include<xmlbasics.h>
#include<dom3.h>
#include<xmdsutils.h>
#include<xmdsclasses.h>

// ******************************************************************************
// ******************************************************************************
//                              xmdsVector public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsVectors=0;  //!< The number of xmds vectors

// ******************************************************************************
xmdsVector::xmdsVector(
                       const xmdsField *const yourField) :
  myField(yourField) {
  if(debugFlag) {
    nxmdsVectors++;
    printf("xmdsVector::xmdsVector\n");
    printf("nxmdsVectors=%li\n",nxmdsVectors);
  }
  myNeedsFFTWRoutines=0;
};

// ******************************************************************************
xmdsVector::~xmdsVector() {
  if(debugFlag) {
    nxmdsVectors--;
    printf("xmdsVector::~xmdsVector\n");
    printf("nxmdsVectors=%li\n",nxmdsVectors);
  }
};

// ******************************************************************************
void xmdsVector::writeDefines(
                              FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsVector::writeDefines\n");
  }

  fprintf(outfile,"// vector %s defines\n",myName.c_str());
  fprintf(outfile,"#define _%s_%s_ncomponents %li\n",field()->name()->c_str(),myName.c_str(),(long)myComponentsNamesList.size());

  unsigned long i=0;
  for(list<XMLString>::const_iterator pXMLString = myComponentsNamesList.begin();pXMLString != myComponentsNamesList.end();pXMLString++) {
    fprintf(outfile,"#define %s _active_%s_%s[_%s_%s_index_pointer + %li]\n",
            pXMLString->c_str(),field()->name()->c_str(),myName.c_str(),field()->name()->c_str(),myName.c_str(),i);
    i++;
  }

  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsVector::writeGlobals(
                              FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsVector::writeGlobals\n");
  }

  fprintf(outfile,"// vector %s globals \n",myName.c_str());
  fprintf(outfile,"\n");

  const char* typeName="";
  if(myType==COMPLEX) {
    typeName="complex";
  }
  else if(myType==DOUBLE) {
    typeName="double";
  }

  if((field()->simulation()->parameters()->usempi)&!(field()->simulation()->parameters()->stochastic)){
    fprintf(outfile,"%s *_%s_%s;\n",typeName,field()->name()->c_str(),myName.c_str());
    fprintf(outfile,"%s* _active_%s_%s;\n",typeName,field()->name()->c_str(),myName.c_str());
    fprintf(outfile,"%s *_%s_%s_work;\n",typeName,field()->name()->c_str(),myName.c_str());
    fprintf(outfile,"\n");
  }
  else {
    fprintf(outfile,"%s *_%s_%s = new %s[_%s_size*_%s_%s_ncomponents];\n",typeName,
            field()->name()->c_str(),myName.c_str(),typeName,
            field()->name()->c_str(),field()->name()->c_str(),myName.c_str());
    fprintf(outfile,"%s* _active_%s_%s;\n",typeName,field()->name()->c_str(),myName.c_str());
    fprintf(outfile,"\n");
  }

  if(!myNeedsFFTWRoutines) {
    return;
  }

  fprintf(outfile,"unsigned long _%s_%s_space;\n",field()->name()->c_str(),myName.c_str());
  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsVector::writePrototypes(
                                 FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsVector::writePrototypes\n");
  }

  fprintf(outfile,"// vector %s prototypes \n",myName.c_str());
  fprintf(outfile,"\n");

  fprintf(outfile,"void _%s_%s_initialise();\n",field()->name()->c_str(),name()->c_str());
  fprintf(outfile,"\n");

  if(!myNeedsFFTWRoutines) {
    return;
  }

  fprintf(outfile,"void _%s_%s_go_space(\n",field()->name()->c_str(),myName.c_str());
  fprintf(outfile,"     const unsigned long& _space);\n");
  fprintf(outfile,"\n");

};

// ******************************************************************************
void xmdsVector::writeRoutines(
                               FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsVector::writeRoutines\n");
  }

  const char *const fieldName = field()->name()->c_str();

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"//           field %s vector %s routines\n",fieldName,myName.c_str());
  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"\n");

  this->writeInitialiseRoutine(outfile);

  if(!myNeedsFFTWRoutines) {
    return;
  }

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _%s_%s_go_space(\n",fieldName,myName.c_str());
  fprintf(outfile,"     const unsigned long& _newspace) {\n");
  fprintf(outfile,"\n");

  fprintf(outfile,"if(_%s_%s_space==_newspace)\n",fieldName,myName.c_str());
  fprintf(outfile,"     return;\n");
  fprintf(outfile,"\n");

  unsigned long fullSpace = myField->geometry()->fullSpace();

  fprintf(outfile,"double _c=1.0;\n");
  fprintf(outfile,"\n");

  fprintf(outfile,"if((_%s_%s_space==0)&(_newspace==%li)) {\n",fieldName,myName.c_str(),fullSpace);
  fprintf(outfile,"\n");

  if(field()->simulation()->parameters()->nThreads>1) {
    fprintf(outfile,"   fftwnd_threads(_num_threads,_%s_forward_plan,_%s_%s_ncomponents,_active_%s_%s,_%s_%s_ncomponents,1,0,0,0);\n",
            fieldName,fieldName,myName.c_str(),
            fieldName,myName.c_str(),fieldName,myName.c_str());
  }
  else if((field()->simulation()->parameters()->usempi)&!(field()->simulation()->parameters()->stochastic)) {
    fprintf(outfile,"   fftwnd_mpi(_%s_forward_plan,_%s_%s_ncomponents,_active_%s_%s,_%s_%s_work,FFTW_TRANSPOSED_ORDER);\n",
            fieldName,fieldName,myName.c_str(),fieldName,myName.c_str(),
            fieldName,myName.c_str());
  }
  else {
    fprintf(outfile,"   fftwnd(_%s_forward_plan,_%s_%s_ncomponents,_active_%s_%s,_%s_%s_ncomponents,1,0,0,0);\n",
            fieldName,fieldName,myName.c_str(),
            fieldName,myName.c_str(),fieldName,myName.c_str());
  }
  fprintf(outfile,"\n");
  for(unsigned long i=0;i<myField->geometry()->nDims();i++) {
    fprintf(outfile,"   _c *= _%s_dx%li/sqrt(2*M_PI);\n",fieldName,i);
  }
  fprintf(outfile,"\n");
  fprintf(outfile,"     _%s_%s_space=_newspace;\n",fieldName,myName.c_str());
  fprintf(outfile,"     }\n");

  fprintf(outfile,"else if((_%s_%s_space==%li)&(_newspace==0)) {\n",fieldName,myName.c_str(),fullSpace);
  fprintf(outfile,"\n");
  if(field()->simulation()->parameters()->nThreads>1) {
    fprintf(outfile,"   fftwnd_threads(_num_threads,_%s_backward_plan,_%s_%s_ncomponents,_active_%s_%s,_%s_%s_ncomponents,1,0,0,0);\n",
            fieldName,fieldName,myName.c_str(),
            fieldName,myName.c_str(),fieldName,myName.c_str());
  }
  else if((field()->simulation()->parameters()->usempi)&!(field()->simulation()->parameters()->stochastic)) {
    fprintf(outfile,"   fftwnd_mpi(_%s_backward_plan,_%s_%s_ncomponents,_active_%s_%s,_%s_%s_work,FFTW_TRANSPOSED_ORDER);\n",
            fieldName,fieldName,myName.c_str(),fieldName,myName.c_str(),
            fieldName,myName.c_str());
  }
  else {
    fprintf(outfile,"   fftwnd(_%s_backward_plan,_%s_%s_ncomponents,_active_%s_%s,_%s_%s_ncomponents,1,0,0,0);\n",
            fieldName,fieldName,myName.c_str(),
            fieldName,myName.c_str(),fieldName,myName.c_str());
  }
  fprintf(outfile,"\n");
  for(unsigned long i=0;i<myField->geometry()->nDims();i++) {
    fprintf(outfile,"   _c *= _%s_dk%li/sqrt(2*M_PI);\n",fieldName,i);
  }
  fprintf(outfile,"\n");
  fprintf(outfile,"     _%s_%s_space=_newspace;\n",fieldName,myName.c_str());
  fprintf(outfile,"     }\n");

  if(myField->geometry()->nDims()>1) {
    if((field()->simulation()->parameters()->usempi)&!(field()->simulation()->parameters()->stochastic)) {
        
        // mixed transforms stuff
        // note this does not need to be ultra efficient since it is probably only
        // ever used in filters, moment stuff, and after initialisation
    
        fprintf(outfile,"else {\n");
        fprintf(outfile,"  long _howmany;\n");
        fprintf(outfile,"  int _lattice;\n");
        fprintf(outfile,"\n");
 
        fprintf(outfile,"  if((_newspace&1)&((_newspace>>1)&1)) {\n");
        fprintf(outfile,"    if(!((_%s_%s_space&1)&&((_%s_%s_space>>1)&1))) {\n",fieldName,myName.c_str(),fieldName,myName.c_str());
        // 1+ -> _%s_%s_space(0) (not swapped)
                                        
        unsigned long two2n=2;
            
        for(unsigned long i=1; i<myField->geometry()->nDims(); i++) {
        
          fprintf(outfile,"      _howmany = _%s_%s_ncomponents;\n",fieldName,myName.c_str());
          for(unsigned long j=myField->geometry()->nDims();j>i+1;j--) {
            fprintf(outfile,"      _howmany *= _%s_lattice%li;\n",fieldName,j-1);
	  }
          fprintf(outfile,"\n");
          fprintf(outfile,"      _lattice = _%s_lattice%li;\n",       fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"      if((!(_%s_%s_space&%li))&&(_%s_%s_space&1)) {\n",fieldName,myName.c_str(),two2n,fieldName,myName.c_str());
          fprintf(outfile,"\n");
          fprintf(outfile,"        fftwnd_plan _plan = fftwnd_create_plan(1,&_lattice,FFTW_FORWARD,FFTW_IN_PLACE);\n");
          fprintf(outfile,"\n");
                            
          fprintf(outfile,"        for(long _i0=0; _i0<local_nx");
          for(unsigned long j=1; j<i; j++) {
            fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	  }
          fprintf(outfile,"; _i0++)\n");
          if(field()->simulation()->parameters()->nThreads>1) {
            fprintf(outfile,"          fftwnd_threads(_num_threads,_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          else {
            fprintf(outfile,"          fftwnd(_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          fprintf(outfile,"\n");
          fprintf(outfile,"          _c *= _%s_dx%li/sqrt(2*M_PI);\n",fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"          delete _plan;\n");
          fprintf(outfile,"        }\n");
          fprintf(outfile,"        else if((_%s_%s_space&%li)&&!(_%s_%s_space&1)) {\n",fieldName,myName.c_str(),two2n,fieldName,myName.c_str());
          fprintf(outfile,"\n");
          fprintf(outfile,"          fftwnd_plan _plan = fftwnd_create_plan(1,&_lattice,FFTW_BACKWARD,FFTW_IN_PLACE);\n");
          fprintf(outfile,"\n");
        
          fprintf(outfile,"          for(long _i0=0; _i0<local_nx");
          for(unsigned long j=1; j<i; j++) {
            fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	  }
          fprintf(outfile,"; _i0++)\n");
          if(field()->simulation()->parameters()->nThreads>1) {
            fprintf(outfile,"            fftwnd_threads(_num_threads,_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          else {
            fprintf(outfile,"            fftwnd(_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
                                    
          fprintf(outfile,"\n");
          fprintf(outfile,"            _c *= _%s_dk%li/sqrt(2*M_PI);\n",fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"            delete _plan;\n");
          fprintf(outfile,"        }\n");
          fprintf(outfile,"\n");
            
          two2n *= 2;
        }

        // if _%s_%s_space&1  BNTFR
        fprintf(outfile,"          if(_%s_%s_space&1) {\n",fieldName,myName.c_str());
        fprintf(outfile,"             fftwnd_mpi(_%s_backward_plan,_%s_%s_ncomponents,_active_%s_%s,_%s_%s_work,FFTW_NORMAL_ORDER);\n\n",
                fieldName,fieldName,myName.c_str(),fieldName,myName.c_str(),
                fieldName,myName.c_str());
                            
        for(unsigned long i=0;i<myField->geometry()->nDims();i++) {
          fprintf(outfile,"              _c *= _%s_dk%li/sqrt(2*M_PI);\n",fieldName,i);
	}
        fprintf(outfile,"              }\n\n");

                                                    
        // BF
        fprintf(outfile,"          fftwnd_mpi(_%s_forward_plan,_%s_%s_ncomponents,_active_%s_%s,_%s_%s_work,FFTW_TRANSPOSED_ORDER);\n\n",
                fieldName,fieldName,myName.c_str(),fieldName,myName.c_str(),
                fieldName,myName.c_str());

        for(unsigned long i=0;i<myField->geometry()->nDims();i++) {
          fprintf(outfile,"              _c *= _%s_dx%li/sqrt(2*M_PI);\n",fieldName,i);
	}
        fprintf(outfile,"\n");

        fprintf(outfile,"                  _%s_%s_space=%li;\n",fieldName,myName.c_str(),two2n-1);

                                        
        fprintf(outfile,"         }\n");
        fprintf(outfile,"      }\n");
        fprintf(outfile,"   else {\n");
        fprintf(outfile,"      if((_%s_%s_space&1)&&((_%s_%s_space>>1)&1)) {\n",fieldName,myName.c_str(),fieldName,myName.c_str());
        // 2+ -> k (swapped)
                                        
        two2n=4;
                
        for(unsigned long i=2; i<myField->geometry()->nDims(); i++) {
            
          fprintf(outfile,"        _howmany = _%s_%s_ncomponents;\n",fieldName,myName.c_str());
          for(unsigned long j=myField->geometry()->nDims();j>i+1;j--) {
            fprintf(outfile,"      _howmany *= _%s_lattice%li;\n",fieldName,j-1);
	  }
          fprintf(outfile,"\n");
          fprintf(outfile,"        _lattice = _%s_lattice%li;\n",       fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"        if(!(_%s_%s_space&%li)) {\n",fieldName,myName.c_str(),two2n);
          fprintf(outfile,"\n");
          fprintf(outfile,"               fftwnd_plan _plan = fftwnd_create_plan(1,&_lattice,FFTW_FORWARD,FFTW_IN_PLACE);\n");
          fprintf(outfile,"\n");
                                
          fprintf(outfile,"                for(long _i0=0; _i0<local_ny_after_transpose*_%s_lattice0",fieldName);
          for(unsigned long j=2; j<i; j++) {
            fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	  }
          fprintf(outfile,"; _i0++)\n");
          if(field()->simulation()->parameters()->nThreads>1) {
            fprintf(outfile,"                       fftwnd_threads(_num_threads,_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          else {
            fprintf(outfile,"                       fftwnd(_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          fprintf(outfile,"\n");
          fprintf(outfile,"                _c *= _%s_dx%li/sqrt(2*M_PI);\n",fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"                delete _plan;\n");
          fprintf(outfile,"        }\n");
          fprintf(outfile,"\n");
                
          two2n *= 2;
        }

        // BFR
        fprintf(outfile,"           fftwnd_mpi(_%s_backward_plan,_%s_%s_ncomponents,_active_%s_%s,_%s_%s_work,FFTW_TRANSPOSED_ORDER);\n\n",
                fieldName,fieldName,myName.c_str(),fieldName,myName.c_str(),
                fieldName,myName.c_str());

        for(unsigned long i=0;i<myField->geometry()->nDims();i++) {
          fprintf(outfile,"             _c *= _%s_dk%li/sqrt(2*M_PI);\n",fieldName,i);
	}
        fprintf(outfile,"\n");

        fprintf(outfile,"                  _%s_%s_space=0;\n",fieldName,myName.c_str());
        fprintf(outfile,"         }\n");
        fprintf(outfile,"      else if((_%s_%s_space+_newspace)&1) {  // if first index changes\n",fieldName,myName.c_str());
        // 1+ -> _%s_%s_space(0)  (not swapped)

                                        
        two2n=2;
                
        for(unsigned long i=1; i<myField->geometry()->nDims(); i++) {
            
          fprintf(outfile,"   _howmany = _%s_%s_ncomponents;\n",fieldName,myName.c_str());
          for(unsigned long j=myField->geometry()->nDims();j>i+1;j--) {
            fprintf(outfile,"   _howmany *= _%s_lattice%li;\n",fieldName,j-1);
	  }
          fprintf(outfile,"\n");
          fprintf(outfile,"   _lattice = _%s_lattice%li;\n",    fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"   if((!(_%s_%s_space&%li))&&(_%s_%s_space&1)) {\n",fieldName,myName.c_str(),two2n,fieldName,myName.c_str());
          fprintf(outfile,"\n");
          fprintf(outfile,"       fftwnd_plan _plan = fftwnd_create_plan(1,&_lattice,FFTW_FORWARD,FFTW_IN_PLACE);\n");
          fprintf(outfile,"\n");
                                
          fprintf(outfile,"                for(long _i0=0; _i0<local_nx");
          for(unsigned long j=1; j<i; j++) {
            fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	  }
          fprintf(outfile,"; _i0++)\n");
          if(field()->simulation()->parameters()->nThreads>1) {
            fprintf(outfile,"               fftwnd_threads(_num_threads,_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          else {
            fprintf(outfile,"               fftwnd(_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          fprintf(outfile,"\n");
          fprintf(outfile,"        _c *= _%s_dx%li/sqrt(2*M_PI);\n",fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"        delete _plan;\n");
          fprintf(outfile,"        }\n");
          fprintf(outfile,"   else if((_%s_%s_space&%li)&&!(_%s_%s_space&1)) {\n",fieldName,myName.c_str(),two2n,fieldName,myName.c_str());
          fprintf(outfile,"\n");
          fprintf(outfile,"        fftwnd_plan _plan = fftwnd_create_plan(1,&_lattice,FFTW_BACKWARD,FFTW_IN_PLACE);\n");
          fprintf(outfile,"\n");
            
          fprintf(outfile,"                for(long _i0=0; _i0<local_nx");
          for(unsigned long j=1; j<i; j++) {
            fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	  }
          fprintf(outfile,"; _i0++)\n");
          if(field()->simulation()->parameters()->nThreads>1) {
            fprintf(outfile,"                  fftwnd_threads(_num_threads,_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          else {
            fprintf(outfile,"                  fftwnd(_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
                                        
          fprintf(outfile,"\n");
          fprintf(outfile,"        _c *= _%s_dk%li/sqrt(2*M_PI);\n",fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"        delete _plan;\n");
          fprintf(outfile,"        }\n");
          fprintf(outfile,"\n");
                
          two2n *= 2;
        }


        // BNTF(R) (to _newspace(0))
        fprintf(outfile,"          if(_newspace&1) {\n");
        fprintf(outfile,"               fftwnd_mpi(_%s_forward_plan,_%s_%s_ncomponents,_active_%s_%s,_%s_%s_work,FFTW_NORMAL_ORDER);\n\n",
                fieldName,fieldName,myName.c_str(),fieldName,myName.c_str(),
                fieldName,myName.c_str());
        for(unsigned long i=0;i<myField->geometry()->nDims();i++) {
          fprintf(outfile,"            _c *= _%s_dx%li/sqrt(2*M_PI);\n",fieldName,i);
	}
        fprintf(outfile,"                  _%s_%s_space=%li;\n",fieldName,myName.c_str(),two2n-1);
        fprintf(outfile,"              }\n");
        fprintf(outfile,"          else {\n");
        fprintf(outfile,"               fftwnd_mpi(_%s_backward_plan,_%s_%s_ncomponents,_active_%s_%s,_%s_%s_work,FFTW_NORMAL_ORDER);\n\n",
                fieldName,fieldName,myName.c_str(),fieldName,myName.c_str(),
                fieldName,myName.c_str());
        for(unsigned long i=0;i<myField->geometry()->nDims();i++) {
          fprintf(outfile,"           _c *= _%s_dk%li/sqrt(2*M_PI);\n",fieldName,i);
	}
        fprintf(outfile,"                  _%s_%s_space=0;\n",fieldName,myName.c_str());
        fprintf(outfile,"              }\n");
                                        
        fprintf(outfile,"         }\n");
        // 1 -> _newspace(1) (not swapped)
                                        
        two2n=2;
                
        for(unsigned long i=1; i<2; i++) {
            
          fprintf(outfile,"   _howmany = _%s_%s_ncomponents;\n",fieldName,myName.c_str());
          for(unsigned long j=myField->geometry()->nDims();j>i+1;j--) {
            fprintf(outfile,"   _howmany *= _%s_lattice%li;\n",fieldName,j-1);
	  }
          fprintf(outfile,"\n");
          fprintf(outfile,"   _lattice = _%s_lattice%li;\n",    fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"   if((!(_%s_%s_space&%li))&&(_newspace&%li)) {\n",fieldName,myName.c_str(),two2n,two2n);
          fprintf(outfile,"\n");
          fprintf(outfile,"       fftwnd_plan _plan = fftwnd_create_plan(1,&_lattice,FFTW_FORWARD,FFTW_IN_PLACE);\n");
          fprintf(outfile,"\n");
                                
          fprintf(outfile,"                for(long _i0=0; _i0<local_nx; _i0++)\n");
          if(field()->simulation()->parameters()->nThreads>1) {
            fprintf(outfile,"               fftwnd_threads(_num_threads,_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          else {
            fprintf(outfile,"               fftwnd(_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          fprintf(outfile,"\n");
          fprintf(outfile,"        _c *= _%s_dx%li/sqrt(2*M_PI);\n",fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"        delete _plan;\n");
          fprintf(outfile,"        }\n");
          fprintf(outfile,"   else if((_%s_%s_space&%li)&&!(_newspace&%li)) {\n",fieldName,myName.c_str(),two2n,two2n);
          fprintf(outfile,"\n");
          fprintf(outfile,"        fftwnd_plan _plan = fftwnd_create_plan(1,&_lattice,FFTW_BACKWARD,FFTW_IN_PLACE);\n");
          fprintf(outfile,"\n");
            
          fprintf(outfile,"                for(long _i0=0; _i0<local_nx; _i0++)\n");
          if(field()->simulation()->parameters()->nThreads>1) {
            fprintf(outfile,"                  fftwnd_threads(_num_threads,_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          else {
            fprintf(outfile,"                  fftwnd(_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
                                        
          fprintf(outfile,"\n");
          fprintf(outfile,"        _c *= _%s_dk%li/sqrt(2*M_PI);\n",fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"        delete _plan;\n");
          fprintf(outfile,"        }\n");
          fprintf(outfile,"\n");
                
          two2n *= 2;
        }


        fprintf(outfile,"      }\n\n\n");
       
        two2n=4;
    
        for(unsigned long i=2; i<myField->geometry()->nDims(); i++) {

          fprintf(outfile,"   _howmany = _%s_%s_ncomponents;\n",fieldName,myName.c_str());
          for(unsigned long j=myField->geometry()->nDims();j>i+1;j--) {
            fprintf(outfile,"   _howmany *= _%s_lattice%li;\n",fieldName,j-1);
	  }
          fprintf(outfile,"\n");
          fprintf(outfile,"   _lattice = _%s_lattice%li;\n",    fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"   if((!(_%s_%s_space&%li))&&(_newspace&%li)) {\n",fieldName,myName.c_str(),two2n,two2n);
          fprintf(outfile,"\n");
          fprintf(outfile,"       fftwnd_plan _plan = fftwnd_create_plan(1,&_lattice,FFTW_FORWARD,FFTW_IN_PLACE);\n");
          fprintf(outfile,"\n");
                    
          fprintf(outfile,"             if((_newspace&1)&&((_newspace>>1)&1))\n");
          fprintf(outfile,"                for(long _i0=0; _i0<local_ny_after_transpose*_%s_lattice0",fieldName);
          for(unsigned long j=2; j<i; j++) {
            fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	  }
          fprintf(outfile,"; _i0++)\n");
          if(field()->simulation()->parameters()->nThreads>1) {
            fprintf(outfile,"               fftwnd_threads(_num_threads,_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          else {
            fprintf(outfile,"               fftwnd(_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          fprintf(outfile,"             else\n");
          fprintf(outfile,"                for(long _i0=0; _i0<local_nx*_%s_lattice1",fieldName);
          for(unsigned long j=2; j<i; j++) {
            fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	  }
          fprintf(outfile,"; _i0++)\n");
          if(field()->simulation()->parameters()->nThreads>1) {
            fprintf(outfile,"               fftwnd_threads(_num_threads,_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          else {
            fprintf(outfile,"               fftwnd(_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          fprintf(outfile,"\n");
          fprintf(outfile,"        _c *= _%s_dx%li/sqrt(2*M_PI);\n",fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"        delete _plan;\n");
          fprintf(outfile,"        }\n");
          fprintf(outfile,"   else if((_%s_%s_space&%li)&&!(_newspace&%li)) {\n",fieldName,myName.c_str(),two2n,two2n);
          fprintf(outfile,"\n");
          fprintf(outfile,"        fftwnd_plan _plan = fftwnd_create_plan(1,&_lattice,FFTW_BACKWARD,FFTW_IN_PLACE);\n");
          fprintf(outfile,"\n");

          fprintf(outfile,"             if((_newspace&1)&&((_newspace>>1)&1))\n");
          fprintf(outfile,"                for(long _i0=0; _i0<local_ny_after_transpose*_%s_lattice0",fieldName);
          for(unsigned long j=2; j<i; j++) {
            fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	  }
          fprintf(outfile,"; _i0++)\n");
          if(field()->simulation()->parameters()->nThreads>1) {
            fprintf(outfile,"                  fftwnd_threads(_num_threads,_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          else {
            fprintf(outfile,"                  fftwnd(_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          fprintf(outfile,"             else\n");
          fprintf(outfile,"                for(long _i0=0; _i0<local_nx*_%s_lattice1",fieldName);
          for(unsigned long j=2; j<i; j++) {
            fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	  }
          fprintf(outfile,"; _i0++)\n");
          if(field()->simulation()->parameters()->nThreads>1) {
            fprintf(outfile,"                  fftwnd_threads(_num_threads,_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
          else {
            fprintf(outfile,"                  fftwnd(_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                    fieldName,myName.c_str());
	  }
                            
          fprintf(outfile,"\n");
          fprintf(outfile,"        _c *= _%s_dk%li/sqrt(2*M_PI);\n",fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"        delete _plan;\n");
          fprintf(outfile,"        }\n");
          fprintf(outfile,"\n");
    
          two2n *= 2;
        }
    
        fprintf(outfile,"   _%s_%s_space=_newspace;\n",fieldName,myName.c_str());
        fprintf(outfile,"   }\n\n");
      }
    else {
        
        // mixed transforms stuff
        // note this does not need to be ultra efficient since it is probably only
        // ever used in filters, moment stuff, and after initialisation
    
        fprintf(outfile,"else {\n");
        fprintf(outfile,"   unsigned long _howmany;\n");
        fprintf(outfile,"   int _lattice;\n");
        fprintf(outfile,"\n");
    
        unsigned long two2n=1;
    
        for(unsigned long i=0; i<myField->geometry()->nDims(); i++) {
    
          fprintf(outfile,"   _howmany = _%s_%s_ncomponents;\n",fieldName,myName.c_str());
          for(unsigned long j=myField->geometry()->nDims();j>i+1;j--) {
            fprintf(outfile,"   _howmany *= _%s_lattice%li;\n",fieldName,j-1);
	  }
          fprintf(outfile,"\n");
          fprintf(outfile,"   _lattice = _%s_lattice%li;\n",    fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"   if((!(_%s_%s_space&%li))&&(_newspace&%li)) {\n",fieldName,myName.c_str(),two2n,two2n);
          fprintf(outfile,"\n");
          fprintf(outfile,"     fftwnd_plan _plan = fftwnd_create_plan(1,&_lattice,FFTW_FORWARD,FFTW_IN_PLACE);\n");
          fprintf(outfile,"\n");
          if(i>0) {
            fprintf(outfile,"           for(unsigned long _i0=0; _i0<_%s_lattice0",fieldName);
            for(unsigned long j=1; j<i; j++) {
              fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	    }
            fprintf(outfile,"; _i0++)\n");
            if(field()->simulation()->parameters()->nThreads>1) {
              fprintf(outfile,"         fftwnd_threads(_num_threads,_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                      fieldName,myName.c_str());
	    }
            else {
              fprintf(outfile,"         fftwnd(_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                      fieldName,myName.c_str());
	    }
          }
          else {
            if(field()->simulation()->parameters()->nThreads>1) {
              fprintf(outfile,"    fftwnd_threads(_num_threads,_plan,_howmany,_active_%s_%s,_howmany,1,0,0,0);\n",
                      fieldName,myName.c_str());
	    }
            else {
              fprintf(outfile,"    fftwnd(_plan,_howmany,_active_%s_%s,_howmany,1,0,0,0);\n",
                      fieldName,myName.c_str());
	    }
	  }
          fprintf(outfile,"\n");
          fprintf(outfile,"        _c *= _%s_dx%li/sqrt(2*M_PI);\n",fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"        delete _plan;\n");
          fprintf(outfile,"        }\n");
          fprintf(outfile,"   else if((_%s_%s_space&%li)&&!(_newspace&%li)) {\n",fieldName,myName.c_str(),two2n,two2n);
          fprintf(outfile,"\n");
          fprintf(outfile,"        fftwnd_plan _plan = fftwnd_create_plan(1,&_lattice,FFTW_BACKWARD,FFTW_IN_PLACE);\n");
          fprintf(outfile,"\n");
          if(i>0) {
            fprintf(outfile,"           for(unsigned long _i0=0; _i0<_%s_lattice0",fieldName);
            for(unsigned long j=1; j<i; j++) {
              fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	    }
            fprintf(outfile,"; _i0++)\n");
            if(field()->simulation()->parameters()->nThreads>1) {
              fprintf(outfile,"            fftwnd_threads(_num_threads,_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                      fieldName,myName.c_str());
	    }
            else {
              fprintf(outfile,"            fftwnd(_plan,_howmany,_active_%s_%s + _i0*_lattice*_howmany,_howmany,1,0,0,0);\n",
                      fieldName,myName.c_str());
	    }
          }
          else {
            if(field()->simulation()->parameters()->nThreads>1) {
              fprintf(outfile,"         fftwnd_threads(_num_threads,_plan,_howmany,_active_%s_%s,_howmany,1,0,0,0);\n",
                      fieldName,myName.c_str());
	    }
            else {
              fprintf(outfile,"    fftwnd(_plan,_howmany,_active_%s_%s,_howmany,1,0,0,0);\n",
                      fieldName,myName.c_str());
	    }
	  }
          fprintf(outfile,"\n");
          fprintf(outfile,"        _c *= _%s_dk%li/sqrt(2*M_PI);\n",fieldName,i);
          fprintf(outfile,"\n");
          fprintf(outfile,"        delete _plan;\n");
          fprintf(outfile,"        }\n");
          fprintf(outfile,"\n");
    
          two2n *= 2;
        }
    
        fprintf(outfile,"   _%s_%s_space=_newspace;\n",fieldName,myName.c_str());
        fprintf(outfile,"   }\n\n");
    }
  }

  if((field()->simulation()->parameters()->usempi)&!(field()->simulation()->parameters()->stochastic)) {
    fprintf(outfile,"for(long _i0=0;_i0<total_local_size*_%s_%s_ncomponents;_i0++) {\n",fieldName,myName.c_str());
  }
  else {
    fprintf(outfile,"for(unsigned long _i0=0;_i0<_%s_size*_%s_%s_ncomponents;_i0++) {\n",fieldName,fieldName,myName.c_str());
  }
  fprintf(outfile,"     _active_%s_%s[_i0].re *= _c;\n",fieldName,myName.c_str());
  fprintf(outfile,"     _active_%s_%s[_i0].im *= _c;\n",fieldName,myName.c_str());
  fprintf(outfile,"     }\n");
  fprintf(outfile,"\n");

  fprintf(outfile,"}\n");
  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsVector::writeInitialisationCall(
                                         FILE *const outfile,
                                         const char *const indent) const {
  if(debugFlag) {
    printf("xmdsVector::writeInitialisationCall\n");
  }

  fprintf(outfile,"%s_%s_%s_initialise();\n",indent,field()->name()->c_str(),name()->c_str());
};

// ******************************************************************************
const XMLString* xmdsVector::name() const {
  if(debugFlag) {
    printf("xmdsVector::name\n");
  }

  return &myName;
};

// ******************************************************************************
void xmdsVector::setName(
                         const XMLString& yourName) {
  if(debugFlag) {
    printf("xmdsVector::setName\n");
  }

  myName = yourName;
};

// ******************************************************************************
xmdsVectorType xmdsVector::vectorType() const {
  if(debugFlag) {
    printf("xmdsVector::vectorType\n");
  }

  return myType;
};

// ******************************************************************************
void xmdsVector::setVectorType(
                               const xmdsVectorType& yourType) {
  if(debugFlag) {
    printf("xmdsVector::setVectorType\n");
  }

  myType = yourType;
};

// ******************************************************************************
unsigned long xmdsVector::nComponents() const {
  if(debugFlag) {
    printf("xmdsVector::nComponents\n");
  }

  return myComponentsNamesList.size();
};

// ******************************************************************************
const XMLString* xmdsVector::componentName(
                                           const unsigned long& index) const {
  if(debugFlag) {
    printf("xmdsVector::componentName\n");
  }

  if(index>=myComponentsNamesList.size()) {
    return 0;
  }

  list<XMLString>::const_iterator pXMLString = myComponentsNamesList.begin();
  for(unsigned long i=0;i<index;i++) {
    pXMLString++;
  }

  return &*pXMLString;
};

// ******************************************************************************
bool xmdsVector::getComponent(
                              const XMLString& ofName,
                              unsigned long& index) const {
  if(debugFlag) {
    printf("xmdsVector::getComponent\n");
  }

  index=0;
  for(list<XMLString>::const_iterator pXMLString = myComponentsNamesList.begin(); pXMLString != myComponentsNamesList.end(); pXMLString++) {
    if(*pXMLString==ofName) {
      return 1;
    }
    index++;
  }

  return 0;
};

// ******************************************************************************
void xmdsVector::setComponents(
                               const list<XMLString>& yourComponentsNamesList) {
  if(debugFlag) {
    printf("xmdsVector::setComponents\n");
  }

  myComponentsNamesList = yourComponentsNamesList;
};

// ******************************************************************************
bool xmdsVector::needsFFTWRoutines() const {
  if(debugFlag) {
    printf("xmdsVector::needsFFTWRoutines\n");
  }

  return myNeedsFFTWRoutines;
};

// ******************************************************************************
void xmdsVector::setNeedsFFTWRoutines() const {
  if(debugFlag) {
    printf("xmdsVector::setNeedsFFTWRoutines\n");
  }

  myNeedsFFTWRoutines=1;
};

// ******************************************************************************
unsigned long xmdsVector::initialSpace() const {
  if(debugFlag) {
    printf("xmdsVector::initialSpace\n");
  }

  return myInitialSpace;
};

// ******************************************************************************
void xmdsVector::setInitialSpace(
                                 const unsigned long& yourInitialSpace) {
  if(debugFlag) {
    printf("xmdsVector::setInitialSpace\n");
  }

  myInitialSpace=yourInitialSpace;
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsVector protected
// ******************************************************************************
// ******************************************************************************

// ******************************************************************************
void xmdsVector::writeInitialiseRoutine(
                                        FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsVector::writeInitialise\n");
  }

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _%s_%s_initialise() {\n",field()->name()->c_str(),myName.c_str());
  fprintf(outfile,"\n");
  fprintf(outfile,"for(unsigned long _i0=0;_i0<_%s_size*_%s_%s_ncomponents;_i0++)\n",
          field()->name()->c_str(),field()->name()->c_str(),name()->c_str());
  fprintf(outfile,"     _%s_%s[_i0] = 0;\n",field()->name()->c_str(),myName.c_str());

  if(myNeedsFFTWRoutines) {
    fprintf(outfile,"\n");
    fprintf(outfile,"_%s_%s_space=%li;\n",field()->name()->c_str(),myName.c_str(),myInitialSpace);
  }

  fprintf(outfile,"}\n");
  fprintf(outfile,"\n");
};

// ******************************************************************************
const xmdsField* xmdsVector::field() const {
  if(debugFlag) {
    printf("xmdsVector::field\n");
  }

  return myField;
};

// ******************************************************************************
bool xmdsVector::space(
                       const long unsigned int& index) const {
  if(debugFlag) {
    printf("xmdsVector::space\n");
  }

  if(index>=field()->geometry()->nDims()) {
    throw xmdsException("Internal range error in xmdsVector::space()");
  }

  return (myInitialSpace >> index)&1;
};