--- /dev/null
+# Try to find gnu scientific library GSL
+# See
+# http://www.gnu.org/software/gsl/ and
+# http://gnuwin32.sourceforge.net/packages/gsl.htm
+#
+# Based on a script of Felix Woelk and Jan Woetzel
+# (www.mip.informatik.uni-kiel.de)
+#
+# It defines the following variables:
+# GSL_FOUND - system has GSL lib
+# GSL_INCLUDE_DIRS - where to find headers
+# GSL_LIBRARIES - full path to the libraries
+# GSL_LIBRARY_DIRS, the directory where the PLplot library is found.
+
+# CMAKE_GSL_CXX_FLAGS = Unix compiler flags for GSL, essentially "`gsl-config --cxxflags`"
+# GSL_LINK_DIRECTORIES = link directories, useful for rpath on Unix
+# GSL_EXE_LINKER_FLAGS = rpath on Unix
+
+set( GSL_FOUND OFF )
+set( GSL_CBLAS_FOUND OFF )
+
+# Windows, but not for Cygwin and MSys where gsl-config is available
+if( WIN32 AND NOT CYGWIN AND NOT MSYS )
+ # look for headers
+ find_path( GSL_INCLUDE_DIR
+ NAMES gsl/gsl_cdf.h gsl/gsl_randist.h
+ )
+ if( GSL_INCLUDE_DIR )
+ # look for gsl library
+ find_library( GSL_LIBRARY
+ NAMES gsl
+ )
+ if( GSL_LIBRARY )
+ set( GSL_INCLUDE_DIRS ${GSL_INCLUDE_DIR} )
+ get_filename_component( GSL_LIBRARY_DIRS ${GSL_LIBRARY} PATH )
+ set( GSL_FOUND ON )
+ endif( GSL_LIBRARY )
+
+ # look for gsl cblas library
+ find_library( GSL_CBLAS_LIBRARY
+ NAMES gslcblas
+ )
+ if( GSL_CBLAS_LIBRARY )
+ set( GSL_CBLAS_FOUND ON )
+ endif( GSL_CBLAS_LIBRARY )
+
+ set( GSL_LIBRARIES ${GSL_LIBRARY} ${GSL_CBLAS_LIBRARY} )
+ endif( GSL_INCLUDE_DIR )
+
+ mark_as_advanced(
+ GSL_INCLUDE_DIR
+ GSL_LIBRARY
+ GSL_CBLAS_LIBRARY
+ )
+else( WIN32 AND NOT CYGWIN AND NOT MSYS )
+ if( UNIX OR MSYS )
+ find_program( GSL_CONFIG_EXECUTABLE gsl-config
+ /usr/bin/
+ /usr/local/bin
+ )
+
+ if( GSL_CONFIG_EXECUTABLE )
+ set( GSL_FOUND ON )
+
+ # run the gsl-config program to get cxxflags
+ execute_process(
+ COMMAND sh "${GSL_CONFIG_EXECUTABLE}" --cflags
+ OUTPUT_VARIABLE GSL_CFLAGS
+ RESULT_VARIABLE RET
+ ERROR_QUIET
+ )
+ if( RET EQUAL 0 )
+ string( STRIP "${GSL_CFLAGS}" GSL_CFLAGS )
+ separate_arguments( GSL_CFLAGS )
+
+ # parse definitions from cflags; drop -D* from CFLAGS
+ string( REGEX MATCHALL "-D[^;]+"
+ GSL_DEFINITIONS "${GSL_CFLAGS}" )
+ string( REGEX REPLACE "-D[^;]+;" ""
+ GSL_CFLAGS "${GSL_CFLAGS}" )
+
+ # parse include dirs from cflags; drop -I prefix
+ string( REGEX MATCHALL "-I[^;]+"
+ GSL_INCLUDE_DIRS "${GSL_CFLAGS}" )
+ string( REPLACE "-I" ""
+ GSL_INCLUDE_DIRS "${GSL_INCLUDE_DIRS}")
+ string( REGEX REPLACE "-I[^;]+;" ""
+ GSL_CFLAGS "${GSL_CFLAGS}")
+
+ message("GSL_DEFINITIONS=${GSL_DEFINITIONS}")
+ message("GSL_INCLUDE_DIRS=${GSL_INCLUDE_DIRS}")
+ message("GSL_CFLAGS=${GSL_CFLAGS}")
+ else( RET EQUAL 0 )
+ set( GSL_FOUND FALSE )
+ endif( RET EQUAL 0 )
+
+ # run the gsl-config program to get the libs
+ execute_process(
+ COMMAND sh "${GSL_CONFIG_EXECUTABLE}" --libs
+ OUTPUT_VARIABLE GSL_LIBRARIES
+ RESULT_VARIABLE RET
+ ERROR_QUIET
+ )
+ if( RET EQUAL 0 )
+ string(STRIP "${GSL_LIBRARIES}" GSL_LIBRARIES )
+ separate_arguments( GSL_LIBRARIES )
+
+ # extract linkdirs (-L) for rpath (i.e., LINK_DIRECTORIES)
+ string( REGEX MATCHALL "-L[^;]+"
+ GSL_LIBRARY_DIRS "${GSL_LIBRARIES}" )
+ string( REPLACE "-L" ""
+ GSL_LIBRARY_DIRS "${GSL_LIBRARY_DIRS}" )
+ else( RET EQUAL 0 )
+ set( GSL_FOUND FALSE )
+ endif( RET EQUAL 0 )
+
+ MARK_AS_ADVANCED(
+ GSL_CFLAGS
+ )
+ message( STATUS "Using GSL from ${GSL_PREFIX}" )
+ else( GSL_CONFIG_EXECUTABLE )
+ message( STATUS "FindGSL: gsl-config not found.")
+ endif( GSL_CONFIG_EXECUTABLE )
+ endif( UNIX OR MSYS )
+endif( WIN32 AND NOT CYGWIN AND NOT MSYS )
+
+if( GSL_FOUND )
+ if( NOT GSL_FIND_QUIETLY )
+ message( STATUS "FindGSL: Found both GSL headers and library" )
+ endif( NOT GSL_FIND_QUIETLY )
+else( GSL_FOUND )
+ if( GSL_FIND_REQUIRED )
+ message( FATAL_ERROR "FindGSL: Could not find GSL headers or library" )
+ endif( GSL_FIND_REQUIRED )
+endif( GSL_FOUND )
--- /dev/null
+#ifndef SIM_HPP
+#define SIM_HPP
+
+#define BESSELCACHE 10000
+#define WCACHE 1000
+
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_integration.h>
+#include <math.h>
+#include <sys/time.h>
+#include <boost/math/special_functions/bessel.hpp>
+
+#include "latlib/neigh.h"
+
+class sim : public o815::sim {
+public:
+ sim(o815 *_O815);
+ unsigned int lsize4;
+ neigh *nb;
+ int *l, *lp, *p;
+ double kappa[2], lambda[2], Mu[2], beta;
+
+private:
+ void _makeSweep();
+ void _newParas();
+ gsl_rng* rangsl;
+ int lineSweep();
+ int lineUpdate(int x, int vardir);
+ int cubeUpdate(int x, int orient);
+ void getCube(int cube[3], int orient);
+ int plaqUpdate(int x, int mu, int nu, int i);
+ int iPlaq(int i, int j);
+ int lpUpdate(int x, int mu, int i);
+ double I(int p);
+ int wArg(int x, int i);
+ double WF(int nom, int denom, int iflav);
+ static double wkern(double t, void *params);
+ double wGsl(int n, int iflav);
+ struct bcache {
+ double val;
+ bool set;
+ };
+ bcache besselCache[BESSELCACHE];
+ struct kernpara {
+ int n;
+ double kappa;
+ double lambda;
+ };
+ struct wcache {
+ double val;
+ double error;
+ bool set;
+ };
+ gsl_integration_cquad_workspace *w;
+ double rhoLine(int xstart, int vardir, int modifier);
+ double rhoPlaq_Fac(int x, int mu, int i, int modifier);
+ double rhoPlaq(int x, int mu, int nu, int i, int modifier);
+ double rhoLp(int x, int mu, int i, int modifier);
+ bool neverComputedB, neverComputedW[2];
+ double rhoCube(int x, int *cube, int modifier);
+ double wLastLambda[2], wLastKappa[2];
+ double lastBesselBeta;
+ int ipow(int base, int expo);
+ wcache wCache[2][WCACHE*WCACHE];
+ double WFl0(int nom, int denom, int iflav);
+};
+
+sim::sim(o815 *_O815) : o815::sim( _O815,
+ sizeof(int)*
+ _O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[1]*(2*4+2*4+6) ) {
+
+ struct timeval tv;
+
+ lsize4 = _O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[0]*_O815->comargs.lsize[1];
+
+ nb = new neigh(4, _O815->comargs.lsize[0], _O815->comargs.lsize[0], _O815->comargs.lsize[0], _O815->comargs.lsize[1]);
+
+ l = (int*)confMem;
+ lp = (int*)(confMem + sizeof(int)*2*lsize4*4);
+ p = (int*)(confMem + sizeof(int)*2*lsize4*4*2);
+
+ neverComputedB = true;
+ neverComputedW[0] = true;
+ neverComputedW[1] = true;
+
+ w = gsl_integration_cquad_workspace_alloc(10000);
+
+ gettimeofday(&tv, NULL);
+ rangsl = gsl_rng_alloc(gsl_rng_ranlxs0);
+ gsl_rng_set(rangsl, 1000000 * tv.tv_sec + tv.tv_usec);
+}
+
+void sim::_makeSweep()
+{
+ lineSweep();
+
+ for( int x=0; x<lsize4; x++) {
+ for(int orient=0; orient<4; orient++)
+ cubeUpdate(x, orient);
+
+ for(int mu=0; mu<3; mu++)
+ for(int nu=mu+1; nu<4; nu++)
+ for(int i=0; i<2; i++)
+ plaqUpdate(x, mu, nu, i);
+
+ for(int mu=0; mu<4; mu++)
+ for(int i=0; i<2; i++)
+ lpUpdate(x, mu, i);
+ }
+}
+
+void sim::_newParas() {
+ kappa[0] = (*O815->paraQ)["kappaone"];
+ kappa[1] = (*O815->paraQ)["kappatwo"];
+ lambda[0] = (*O815->paraQ)["lambdaone"];
+ lambda[1] = (*O815->paraQ)["lambdatwo"];
+ beta = (*O815->paraQ)["beta"];
+ Mu[0] = (*O815->paraQ)["muone"];
+ Mu[1] = (*O815->paraQ)["mutwo"];
+
+ memset(confMem, 0, sizeof(int)*lsize4*(2*4+2*4+6));
+}
+
+int sim::lineSweep()
+{
+ for (int vardir = 0; vardir < 4; vardir++)
+ for (int x0 = 0; x0 < O815->comargs.lsize[0]; x0++) {
+ for (int x1 = 0; x1 < O815->comargs.lsize[0]; x1++) {
+ for (int x2 = 0; x2 < O815->comargs.lsize[0]; x2++) {
+ for (int x3 = 0; x3 < O815->comargs.lsize[1]; x3++) {
+ lineUpdate(x0
+ + x1 * O815->comargs.lsize[0]
+ + x2 * O815->comargs.lsize[0] * O815->comargs.lsize[0]
+ + x3 * O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0], vardir);
+ if(vardir==3) break;
+ }
+ if(vardir==2) break;
+ }
+ if(vardir==1) break;
+ }
+ if(vardir==0) break;
+ }
+}
+
+int sim::lineUpdate(int x, int vardir)
+{
+ int modifier = 1 - 2 * gsl_rng_uniform_int(rangsl, 2);
+
+ if ( gsl_rng_uniform(rangsl) < rhoLine(x, vardir, modifier) )
+ {
+ const int stepsize = ipow(O815->comargs.lsize[0], vardir);
+
+ for (int xvar = 0; xvar < O815->comargs.lsize[vardir/3]; xvar++)
+ for (int i = 0; i < 2; i++)
+ l[ i*lsize4*4 + (x+xvar*stepsize)*4 + vardir ] += modifier;
+
+ return 1;
+ }
+
+ return 0;
+}
+
+int sim::cubeUpdate(int x, int orient)
+{
+ int modifier = 1 - 2 * gsl_rng_uniform_int(rangsl, 2);
+ int cube[3];
+ getCube(cube, orient);
+
+ if ( gsl_rng_uniform(rangsl) < rhoCube(x, cube, modifier) )
+ {
+ p[ x*6 + iPlaq( cube[0],cube[1] ) ] += modifier;
+ p[ x*6 + iPlaq( cube[0],cube[2] ) ] -= modifier;
+ p[ x*6 + iPlaq( cube[1],cube[2] ) ] += modifier;
+ p[ (*nb)[x*8+cube[0]]*6 + iPlaq( cube[1],cube[2] ) ] -= modifier;
+ p[ (*nb)[x*8+cube[1]]*6 + iPlaq( cube[0],cube[2] ) ] += modifier;
+ p[ (*nb)[x*8+cube[2]]*6 + iPlaq( cube[0],cube[1] ) ] -= modifier;
+
+ return 1;
+ }
+
+ return 0;
+}
+
+void sim::getCube(int cube[3], int orient)
+{
+ int skipped=0;
+ for (int i=0; i<4; i++) {
+ if ( i == orient ) {
+ skipped++;
+ continue;
+ }
+ cube[i-skipped] = i;
+ }
+}
+
+int sim::plaqUpdate(int x, int mu, int nu, int i)
+{
+ int modifier = 1 - 2 * gsl_rng_uniform_int(rangsl, 2);
+ int lmodsign = 1 - 2*i;
+
+ if ( gsl_rng_uniform(rangsl) < rhoPlaq(x, mu, nu, i, modifier) )
+ {
+ p[ x*6 + iPlaq(mu,nu) ] = p[ x*6 + iPlaq(mu,nu) ] + modifier;
+ l[ i*lsize4*4 + x*4 + mu ] = l[ i*lsize4*4 + x*4 + mu ] - lmodsign*modifier;
+ l[ i*lsize4*4 + x*4 + nu ] = l[ i*lsize4*4 + x*4 + nu ] + lmodsign*modifier;
+ l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] = l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] + lmodsign*modifier;
+ l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] = l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] - lmodsign*modifier;
+ return 1;
+ }
+
+ return 0;
+}
+
+int sim::iPlaq(int i, int j)
+{
+ if ( i == 0 ) return j-1;
+ else if ( i == 1 ) return j+1;
+ else return j+2;
+}
+
+int sim::lpUpdate(int x, int mu, int i)
+{
+ int modifier = 1 - 2 * gsl_rng_uniform_int(rangsl, 2);
+
+ if ( (lp[ i*lsize4*4 + x*4 + mu ] + modifier) < 0 )
+ return 0;
+
+ if ( gsl_rng_uniform(rangsl) < rhoLp(x, mu, i, modifier) ) {
+ lp[ i*lsize4*4 + x*4 + mu ] += modifier;
+ return 1;
+ }
+
+ return 0;
+}
+
+double sim::I(int p)
+{
+ /* http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/bessel/mbessel.html */
+
+ p = abs(p);
+
+ if ( beta != lastBesselBeta || neverComputedB )
+ {
+ memset(besselCache, 0, sizeof(bcache)*BESSELCACHE);
+ lastBesselBeta = beta;
+ neverComputedB = false;
+ }
+
+ if ( p < BESSELCACHE )
+ {
+ if ( besselCache[p].set ) return besselCache[p].val;
+ else
+ {
+ besselCache[p].val = boost::math::cyl_bessel_i(p, beta);
+ besselCache[p].set = true;
+ return besselCache[p].val;
+ }
+ }
+ else
+ {
+ return boost::math::cyl_bessel_i(p, beta);
+ }
+}
+
+int sim::wArg(int x, int i)
+{
+ int warg = 0;
+
+ for (int nu=0; nu<4; nu++)
+ warg += abs( l[ i*lsize4*4 + x*4 + nu ] )
+ + 2 * ( lp[ i*lsize4*4 + x*4 + nu ] + lp[ i*lsize4*4 + (*nb)[x*8+nu+4]*4 + nu ] )
+ + abs( l[ i*lsize4*4 + (*nb)[x*8+nu+4]*4 + nu ] );
+
+ return warg;
+}
+
+double sim::WF(int nom, int denom, int iflav)
+{
+ if ( nom == denom ) return 1;
+
+ if ( wLastLambda[iflav] != lambda[iflav] || wLastKappa[iflav] != kappa[iflav] || neverComputedW[iflav] )
+ {
+ memset(wCache[iflav], 0, sizeof(wcache)*WCACHE*WCACHE);
+ wLastLambda[iflav] = lambda[iflav];
+ wLastKappa[iflav] = kappa[iflav];
+ neverComputedW[iflav] = false;
+ }
+
+ if ( nom/2 < WCACHE && denom/2 < WCACHE )
+ {
+ if ( wCache[iflav][ (nom/2)*WCACHE + denom/2 ].set )
+ return wCache[iflav][ (nom/2)*WCACHE + denom/2 ].val;
+ else
+ {
+ /* integral can be solved analytically for lambda==0 */
+ if( lambda == 0 )
+ wCache[iflav][ (nom/2)*WCACHE + denom/2 ].val = WFl0(nom, denom, iflav);
+ else
+ wCache[iflav][ (nom/2)*WCACHE + denom/2 ].val = wGsl(nom, iflav) / wGsl(denom, iflav);
+
+ wCache[iflav][ (nom/2)*WCACHE + denom/2 ].set = true;
+
+ return wCache[iflav][ (nom/2)*WCACHE + denom/2 ].val;
+ }
+ }
+ else
+ {
+ if( lambda == 0 )
+ return WFl0(nom, denom, iflav);
+ else
+ return wGsl(nom, iflav) / wGsl(denom, iflav);
+ }
+}
+
+double sim::wkern(double t, void *params)
+{
+ kernpara *para = (kernpara*)params;
+ double wkern = pow((1-t)/t, para->n + 1) * exp( -para->kappa*pow((1-t)/t, 2) -para->lambda*pow((1-t)/t, 4) ) / pow(t,2);
+ return wkern;
+}
+
+double sim::wGsl(int n, int iflav)
+{
+ double result;
+ gsl_function gslF;
+ kernpara para;
+ int status;
+
+ para.n = n;
+ para.kappa = kappa[iflav];
+ para.lambda = lambda[iflav];
+
+ gslF.function = &wkern;
+ gslF.params = ¶
+
+ status = gsl_integration_cquad(&gslF, 0, 1, 0, 1e-7, w, &result, NULL, NULL);
+
+ return result;
+}
+
+double sim::rhoLine(int xstart, int vardir, int modifier)
+{
+ double rho=1;
+ int tmpwarg;
+
+ const int stepsize = ipow(O815->comargs.lsize[0],vardir);
+ for(int xvar=0; xvar<O815->comargs.lsize[vardir/3]; xvar++)
+ for(int i=0; i<2; i++)
+ {
+ const int x = xstart + xvar*stepsize;
+
+ if( abs( l[ i*lsize4*4 + x*4 + vardir ] + modifier ) > abs( l[ i*lsize4*4 + x*4 + vardir ] ) )
+ rho /= abs( l[ i*lsize4*4 + x*4 + vardir ] + modifier ) + lp[ i*lsize4*4 + x*4 + vardir ];
+ else
+ rho *= abs( l[ i*lsize4*4 + x*4 + vardir ] ) + lp[ i*lsize4*4 + x*4 + vardir ];
+
+ tmpwarg = wArg(x,i);
+ rho *= WF( tmpwarg
+ - abs( l[ i*lsize4*4 + x*4 + vardir ] )
+ + abs( l[ i*lsize4*4 + x*4 + vardir ] + modifier )
+ - abs( l[ i*lsize4*4 + (*nb)[x*8+vardir+4]*4 + vardir ] )
+ + abs( l[ i*lsize4*4 + (*nb)[x*8+vardir+4]*4 + vardir ] + modifier ),
+ tmpwarg, i );
+
+ if(vardir==3)
+ rho *= exp(modifier*Mu[i]);
+ }
+
+ return rho;
+}
+
+double sim::rhoPlaq(int x, int mu, int nu, int i, int modifier)
+{
+ double rho = 1;
+ int lmodsign = 1 - 2*i;
+ int tmpwarg;
+
+ rho *= I( p[ x*6 + iPlaq(mu,nu) ] + modifier ) / I( p[ x*6 + iPlaq(mu,nu) ] );
+
+ tmpwarg = wArg(x,i);
+ rho *= WF( tmpwarg
+ - abs( l[ i*lsize4*4 + x*4 + mu] ) - abs( l[ i*lsize4*4 + x*4 + nu ] )
+ + abs( l[ i*lsize4*4 + x*4 + mu] - lmodsign*modifier ) + abs( l[ i*lsize4*4 + x*4 + nu ] + lmodsign*modifier ),
+ tmpwarg, i );
+
+ tmpwarg = wArg((*nb)[x*8+mu],i);
+ rho *= WF( tmpwarg
+ - abs( l[ i*lsize4*4 + x*4 + mu] ) - abs( l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] )
+ + abs( l[ i*lsize4*4 + x*4 + mu] - lmodsign*modifier ) + abs( l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] - lmodsign*modifier ),
+ tmpwarg, i );
+
+ tmpwarg = wArg((*nb)[(*nb)[x*8+mu]*8+nu],i);
+ rho *= WF( tmpwarg
+ - abs( l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] ) - abs( l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] )
+ + abs( l[ i*lsize4*4 + (*nb)[x*8+mu]*4 + nu ] - lmodsign*modifier ) + abs( l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] + lmodsign*modifier ),
+ tmpwarg, i );
+
+ tmpwarg = wArg((*nb)[x*8+nu],i);
+ rho *= WF( tmpwarg
+ - abs( l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] ) - abs( l[ i*lsize4*4 + x*4 + nu ] )
+ + abs( l[ i*lsize4*4 + (*nb)[x*8+nu]*4 + mu ] + lmodsign*modifier ) + abs( l[ i*lsize4*4 + x*4 + nu ] + lmodsign*modifier ),
+ tmpwarg, i );
+
+ rho *= rhoPlaq_Fac(x, mu, i, -lmodsign*modifier);
+ rho *= rhoPlaq_Fac(x, nu, i, +lmodsign*modifier);
+ rho *= rhoPlaq_Fac((*nb)[x*8+mu], nu, i, -lmodsign*modifier);
+ rho *= rhoPlaq_Fac((*nb)[x*8+nu], mu, i, +lmodsign*modifier);
+
+ return rho;
+}
+
+double sim::rhoPlaq_Fac(int x, int mu, int i, int modifier)
+{
+ if( abs( l[ i*lsize4*4 + x*4 + mu ] + modifier ) > abs( l[ i*lsize4*4 + x*4 + mu] ) )
+ return 1.0 / ( abs( l[ i*lsize4*4 + x*4 + mu ] + modifier ) + lp[ i*lsize4*4 + x*4 + mu ] );
+ else
+ return ( abs( l[ i*lsize4*4 + x*4 + mu ] ) + lp[ i*lsize4*4 + x*4 + mu ] );
+}
+
+double sim::rhoLp(int x, int mu, int i, int modifier)
+{
+ double rho=1;
+ int tmpwarg;
+
+ if(modifier > 0)
+ rho /= ( abs(l[ i*lsize4*4 + x*4 + mu ]) + lp[ i*lsize4*4 + x*4 + mu ] + modifier ) * ( lp[ i*lsize4*4 + x*4 + mu ] + modifier );
+ else
+ rho *= ( abs(l[ i*lsize4*4 + x*4 + mu ]) + lp[ i*lsize4*4 + x*4 + mu ] ) * lp[ i*lsize4*4 + x*4 + mu ];
+
+ tmpwarg = wArg(x,i);
+ rho *= WF( tmpwarg + 2*modifier, tmpwarg, i );
+ tmpwarg = wArg((*nb)[x*8+mu],i);
+ rho *= WF( tmpwarg + 2*modifier, tmpwarg, i );
+
+ return rho;
+}
+
+double sim::rhoCube(int x, int *cube, int modifier)
+{
+ double rho=1;
+ rho *= I( p[ x*6 + iPlaq( cube[0],cube[1] ) ] + modifier ) / I( p[ x*6 + iPlaq( cube[0],cube[1] ) ] );
+ rho *= I( p[ x*6 + iPlaq( cube[0],cube[2] ) ] - modifier ) / I( p[ x*6 + iPlaq( cube[0],cube[2] ) ] );
+ rho *= I( p[ x*6 + iPlaq( cube[1],cube[2] ) ] + modifier ) / I( p[ x*6 + iPlaq( cube[1],cube[2] ) ] );
+ rho *= I( p[ (*nb)[x*8+cube[0]]*6 + iPlaq( cube[1],cube[2] ) ] - modifier ) / I( p[ (*nb)[x*8+cube[0]]*6 + iPlaq( cube[1],cube[2] ) ] );
+ rho *= I( p[ (*nb)[x*8+cube[1]]*6 + iPlaq( cube[0],cube[2] ) ] + modifier ) / I( p[ (*nb)[x*8+cube[1]]*6 + iPlaq( cube[0],cube[2] ) ] );
+ rho *= I( p[ (*nb)[x*8+cube[2]]*6 + iPlaq( cube[0],cube[1] ) ] - modifier ) / I( p[ (*nb)[x*8+cube[2]]*6 + iPlaq( cube[0],cube[1] ) ] );
+ return rho;
+}
+
+int sim::ipow(int base, int expo)
+{
+ int ipow = 1;
+ while( expo>0 )
+ {
+ ipow *= base;
+ expo--;
+ }
+ return ipow;
+}
+
+double sim::WFl0(int nom, int denom, int iflav)
+{
+ double wfl0 = 1;
+
+ if (nom > denom) {
+ for (int i=nom/2; i>denom/2; i--)
+ wfl0 *= i;
+ return wfl0 * pow(kappa[iflav], 0.5*(denom-nom));
+ }
+ else {
+ for (int i=denom/2; i>nom/2; i--)
+ wfl0 *= i;
+ return (1.0 / wfl0) * pow(kappa[iflav], 0.5*(denom-nom));
+ }
+}
+
+#endif
--- /dev/null
+#include <complex>
+#include <cstring>
+
+#include "latlib/o815/o815.h"
+
+#include "sim.hpp"
+
+o815 *O815;
+sim *Sim;
+
+//const complex<double> _i_ = complex<double>(0.0,1.0);
+
+//#include "obs_phi2.hpp"
+//#include "obs_plaq.hpp"
+//#include "obs_corr.hpp"
+
+o815::comoption specOps[] = {
+ { "kappaone", required_argument, NULL, 'r', "set inverse mass kappa_1", "min:max:inc" },
+ { "kappatwo", required_argument, NULL, 's', "set inverse mass kappa_2", "min:max:inc" },
+ { "beta", required_argument, NULL, 'b', "set inverse gauge coupling beta", "min:max:inc" },
+ { "lambdaone", required_argument, NULL, 'p', "set quartic coupling lambda_1", "min:max:inc" },
+ { "lambdatwo", required_argument, NULL, 'q', "set quartic coupling lambda_2", "min:max:inc" },
+ { "kappa", required_argument, NULL, 'k', "set kappa_1 = kappa_2 = kappa", "min:max:inc" },
+ { "lambda", required_argument, NULL, 'l', "set lambda_1 = lambda_2 = lambda", "min:max:inc" },
+ { "muone", required_argument, NULL, 'n', "set chemical potential mu_1", "min:max:inc" },
+ { "mutwo", required_argument, NULL, 'u', "set chemical potential mu_2", "min:max:inc" },
+ { "mu", required_argument, NULL, 'm', "set mu_1 = mu_2 = mu", "min:max:inc" },
+ { "", 0, NULL, 0, "", "" }
+};
+
+void parseSpecOps()
+{
+ for (int isopt = 0; isopt < O815->parsedSpecOps.size(); isopt++)
+ switch(O815->parsedSpecOps[isopt].first) {
+ case 'r':
+ O815->paraQ->addRange("kappaone", O815->parsedSpecOps[isopt].second);
+ break;
+ case 's':
+ O815->paraQ->addRange("kappatwo", O815->parsedSpecOps[isopt].second);
+ break;
+ case 'b':
+ O815->paraQ->addRange("beta", O815->parsedSpecOps[isopt].second);
+ break;
+ case 'p':
+ O815->paraQ->addRange("lambdaone", O815->parsedSpecOps[isopt].second);
+ break;
+ case 'q':
+ O815->paraQ->addRange("lambdatwo", O815->parsedSpecOps[isopt].second);
+ break;
+ case 'l':
+ O815->paraQ->addRange("lambdaone", O815->parsedSpecOps[isopt].second);
+ O815->paraQ->linkParas("lambdatwo", "lambdaone");
+ break;
+ case 'k':
+ O815->paraQ->addRange("kappaone", O815->parsedSpecOps[isopt].second);
+ O815->paraQ->linkParas("kappatwo", "kappaone");
+ break;
+ case 'm':
+ O815->paraQ->addRange("muone", O815->parsedSpecOps[isopt].second);
+ O815->paraQ->linkParas("mutwo", "muone");
+ break;
+ case 'n':
+ O815->paraQ->addRange("muone", O815->parsedSpecOps[isopt].second);
+ break;
+ case 'u':
+ O815->paraQ->addRange("mutwo", O815->parsedSpecOps[isopt].second);
+ break;
+ }
+}
+
+void helpHeader()
+{
+ cout << "Usage: ./u1casc-flux [OPTIONS] [obs1] ... [obsN]" << endl << endl;
+}
+
+void parseLonelyArgs()
+{
+ /*
+ for (vector<char*>::iterator lonit = O815->lonelyArgs.begin(); lonit != O815->lonelyArgs.end(); ++lonit) {
+ if ( strcmp(*lonit, "phi2") == 0 ) {
+ *O815->out->log << "MASTER: registered observable: phi2" << endl << flush;
+ O815->observables.push_back(new obs_phi2(O815));
+ }
+ else if ( strcmp(*lonit, "plaq") == 0 ) {
+ *O815->out->log << "MASTER: registered observable: plaq" << endl << flush;
+ O815->observables.push_back(new obs_plaq(O815));
+ }
+ else if ( strcmp(*lonit, "corr") == 0 ) {
+ *O815->out->log << "MASTER: registered observable: corr" << endl << flush;
+ O815->observables.push_back(new obs_corr(O815));
+ }
+ }
+ */
+}
+
+int main (int argc, char *argv[])
+{
+ O815 = new o815(argc, argv, "u1casc-flux", specOps, &helpHeader);
+
+ O815->addPara("beta", 1);
+ O815->addPara("lambdaone", 1);
+ O815->addPara("lambdatwo", 1);
+ O815->addPara("kappaone", 1);
+ O815->addPara("kappatwo", 1);
+ O815->addPara("muone", 0);
+ O815->addPara("mutwo", 0);
+
+ parseSpecOps();
+ O815->postParaInit();
+ O815->Sim = new sim(O815);
+ parseLonelyArgs();
+ O815->mainLoop();
+
+ delete O815;
+ return 0;
+}
projects / phys / u1casc.git / commitdiff