xmdsvector.cc

XMDS is a code generator that integrates equations. You write them down in human readable form in a

/*
 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",