xmdsvector.cc

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

          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;
};