From: Alexander Schmidt Date: Wed, 26 Feb 2014 20:15:32 +0000 (+0100) Subject: Migrated u1ext-fancy simulation algorithm. X-Git-Url: http://git.treefish.org/~alex/phys/u1casc.git/commitdiff_plain/7a5913c428ac75aa2de2330fa926e8deb3457112?ds=sidebyside Migrated u1ext-fancy simulation algorithm. --- diff --git a/u1casc-flux/.gitignore b/u1casc-flux/.gitignore new file mode 100644 index 0000000..77f2090 --- /dev/null +++ b/u1casc-flux/.gitignore @@ -0,0 +1,8 @@ +CMakeFiles +*~ +Makefile +cmake_install.cmake +CMakeCache.txt +u1casc-flux +data +out diff --git a/u1casc-flux/CMakeLists.txt b/u1casc-flux/CMakeLists.txt new file mode 100644 index 0000000..9b85594 --- /dev/null +++ b/u1casc-flux/CMakeLists.txt @@ -0,0 +1,20 @@ +cmake_minimum_required(VERSION 2.8) + +project(u1casc-flux) + +# #PROFILING +# set(CMAKE_CXX_FLAGS -pg) +# set(CMAKE_EXE_LINKER_FLAGS -pg) +# set(CMAKE_SHARED_LINKER_FLAGS -pg) + +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "./") + +FIND_PACKAGE(GSL REQUIRED) +include_directories(${GSL_INCLUDE_DIRS}) + +SET(CMAKE_BUILD_TYPE Release) + +add_subdirectory(latlib) + +add_executable(u1casc-flux u1casc-flux.cpp) +target_link_libraries(u1casc-flux o815 ${GSL_LIBRARIES} lat_neigh) diff --git a/u1casc-flux/FindGSL.cmake b/u1casc-flux/FindGSL.cmake new file mode 100644 index 0000000..486be73 --- /dev/null +++ b/u1casc-flux/FindGSL.cmake @@ -0,0 +1,135 @@ +# 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 ) diff --git a/u1casc-flux/latlib b/u1casc-flux/latlib new file mode 120000 index 0000000..7076779 --- /dev/null +++ b/u1casc-flux/latlib @@ -0,0 +1 @@ +../latlib \ No newline at end of file diff --git a/u1casc-flux/sim.hpp b/u1casc-flux/sim.hpp new file mode 100644 index 0000000..9eca3d7 --- /dev/null +++ b/u1casc-flux/sim.hpp @@ -0,0 +1,478 @@ +#ifndef SIM_HPP +#define SIM_HPP + +#define BESSELCACHE 10000 +#define WCACHE 1000 + +#include +#include +#include +#include +#include +#include + +#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; xparaQ)["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; xvarcomargs.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 diff --git a/u1casc-flux/u1casc-flux.cpp b/u1casc-flux/u1casc-flux.cpp new file mode 100644 index 0000000..0061116 --- /dev/null +++ b/u1casc-flux/u1casc-flux.cpp @@ -0,0 +1,116 @@ +#include +#include + +#include "latlib/o815/o815.h" + +#include "sim.hpp" + +o815 *O815; +sim *Sim; + +//const complex _i_ = complex(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::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; +}