/*
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;
};
syntax highlighted by Code2HTML, v. 0.9.1