xmdsvector.cc

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

                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");